

#Results of 'An ensemble of ordered logistic regression and random forest for child garment size matching'

source('enporffunctions.R') #functions for compting new measures

####################################################
#reading data and description for boys for t-shirt
####################################################


d=read.table("tshirt.txt", sep = ",")


medidas <- c("X10","X0030_mk","X1510","X1520","X1530","X3020_mk",
             "X3030_LR","X4020","X4030","X4040_mk","X4080_LR","X4510","X5020","X5030_mk","X5040_mk",
             "X6010_mk","X6011_mk","X6012_mk","X6510_mk","X7520","X7525","X7540","X8010_LR",
             "X8030_LR","X8040_LR","X8520_LR","X8550_LR")
             #,"X0100","X0110","X9035_LR","X9510_LR","X9520_LR","X9550_LR")

datos=d[,c("User_code","Size_fit","Expert_evaluation", medidas)]



#boys

db=datos[d$SEXO=="V",]
#girls
#db=datos[d$SEXO=="M",]
#db=db[-22,] #nan in expert evaluation 


attach(db)


library(party)
 
ninos=unique(db[,1])



set.seed(2015)

da=db[,4:30] #matrix with the predictors only
#db[,2] is Size_fit, the factor with the classes

####################################################
#mtry with less error rate with oob
####################################################
a=c()
for(i in 5:15){
my_cforest_control <- cforest_control(teststat = "max",
    testtype = "Teststatistic", mincriterion = 0, ntree = 500, mtry = i,
    replace = TRUE) #mimic RF Breiman
    set.seed(2015)
marboli=cforest(factor(Size_fit,ordered=TRUE)+Expert_evaluation~.,data=da,controls=my_cforest_control)

mq=matrix(unlist(predict(marboli,OOB=T)),ncol=4,byrow=TRUE)
cl=max.col(mq[,1:3])-2
a[i-4]= sum(cl==db[,2]) 
}

# a [1] 74 77 76 77 78 76 75 76 76 77 79 boys #i=15 max 
#girls: > a  [1] 88 90 89 88 87 87 90 90 90 90 88 # i=6 max

 
 #eoo=79/109=0.7247706
 #girls=>   90/122= 0.7377049
 
 
####################
#RF with $cforest_unbiased$ as Strobl 2007
###################


set.seed(2015)

marbol=cforest(factor(Size_fit,ordered=TRUE)+Expert_evaluation~.,data=db,controls=cforest_unbiased()) 


qq=predict(marbol,OOB=T)
pm=matrix(unlist(qq),ncol=4,byrow=TRUE)
cl=max.col(pm[,1:3])-2

sum(cl==db[,2]) #73/109= 0.6697248

####################################################
#computing proximities and displaying them 
####################################################

i=15
 my_cforest_control <- cforest_control(teststat = "max",
    testtype = "Teststatistic", mincriterion = 0, ntree = 500, mtry = i,
    replace = TRUE)
    set.seed(2015)
 
marbol=cforest(factor(Size_fit,ordered=TRUE)+Expert_evaluation~.,data=da,controls=my_cforest_control) 

 xxt=proximity(marbol) #these proximities are obtained with training data, not oob
 
 #with OBB
 
  ntetot=marbol@get_where(newdata = da)
  whinbag <- marbol@where #zeros out-of-bag
 
 whout=ntetot
 
 for (i in 1:my_cforest_control@ntree){
 
  aux1=  whout[[i]]
  aux2= whinbag[[i]]
  
  aux1[aux2!=0]=0
   
  whout[[i]]=as.integer(aux1)
  
  }
 
 
rn <- rownames(da)
  prox <- .Call("R_proximity", whout, package = "party")
    prox <- matrix(unlist(prox), ncol = length(prox))
    rownames(prox) <- rn
    colnames(prox) <- rn
   xx= prox #proximities estimated with obb
 
 
 
 #cmds
ramd <- cmdscale(1 - xx, eig=TRUE)


rp=ramd$points

plot(rp,type='n',xlab="Dim 1", ylab="Dim 2")
points(rp[Size_fit==1,],col=4,pch=22) #blue
points(rp[Size_fit==-1,],col=2,pch=21) #red
points(rp[Size_fit==0,],col=3,pch=24) #green

ramd$eig #eigenvalues

cumsum(ramd$eig^2/sum(ramd$eig^2)) #0.4583618 0.6690651 

#h-plot
di=1-xx

windows()
#hplot, with 2 dimensions
n=dim(da)[1]
p=prcomp(di)
v=p$rotation
pr=p$sdev
pv=cbind(pr[1]*v[1:n,1],pr[2]*v[1:n,2])
#pv=cbind(pr[1]*v[1:n,1],pr[2]*v[1:n,2],pr[3]*v[1:n,3])


cumsum(c(pr[1]^4,pr[2]^4))/sum(pr^4) #first and second 0.7075972 0.8840140

plot(pv,type='n',xlab="Dim 1", ylab="Dim 2")
points(pv[Size_fit==1,],col=4,pch=22) #blue
points(pv[Size_fit==-1,],col=2,pch=21) #red
points(pv[Size_fit==0,],col=3,pch=24) #green
points(pv[c(4,54,83),],col=1,pch=3) #black cross

#plot(pv[,2:3],type='n',xlab="Dim 2", ylab="Dim 3")
#points(pv[Size_fit==1,2:3],col=4,pch=22) #blue
#points(pv[Size_fit==-1,2:3],col=2,pch=21) #red
#points(pv[Size_fit==0,2:3],col=3,pch=24) #green
#points(pv[c(4,54,83),2:3],col=1,pch=3) #black cross

#plot(pv[,c(1,3)],type='n',xlab="Dim 1", ylab="Dim 3")
#points(pv[Size_fit==1,c(1,3)],col=4,pch=22) #blue
#points(pv[Size_fit==-1,c(1,3)],col=2,pch=21) #red
#points(pv[Size_fit==0,c(1,3)],col=3,pch=24) #green
#points(pv[c(4,54,83),c(1,3)],col=1,pch=3) #black cross

#spl=scatterplot3d(pv[,1],pv[,2],pv[,3],type='n',xlab="Dim 1", ylab="Dim 2",zlab="Dim 3",angle=a)
#spl$points3d(pv[Size_fit==1,1:3],col=4,pch=22) #blue
#spl$points3d(pv[Size_fit==-1,],col=2,pch=21) #red
#spl$points3d(pv[Size_fit==0,],col=3,pch=24) #green
#spl$points3d(pv[c(4,54,83),],col=1,pch=3) #black cross
 
######################################################################
#computing outliers, class probabilities, Q_i^k, TM_i^k and prototypes
######################################################################
 
 library(randomForest)

o=outlier(xx, cls=db[,2]) #outlier measure
 
# > max(o) 6.496779

#class probabilities
 mq=matrix(unlist(predict(marbol,OOB=T)),ncol=4,byrow=TRUE)
 
#Qik
k=21 #the number of members in class -1 for boys is 22
#min(table(db[,2]))-1 #it's the number by defect in classCenter
set.seed(2015)
situa=qik(xx,db[,2],k) #Q_i^k computing

protu=c() #probability given by RF for true class
 
 for(i in 1: dim(da)[1]){
 protu[i]= mq[i,db[i,2]+2]
 }
summary((abs(situa-protu))) 

situall=tmik(xx,db[,2],k) #TM_i^k computing

#representants
which(situall==1)
# 51 54 83


#cbind(pv[db[,2]==-1,],situa[db[,2]==-1],situall[db[,2]==-1])
#mq[db[,2]==-1,1]


gr=db[,2]
for (k in seq(21,1,-1)){
set.seed(2015)
situa=qik(xx,gr,k) #Q_i^k computing
 situall=tmik(xx,gr,k)
 wll=which(situall==1)
 wh=situa[wll]==1
 if(any(wh)){cat("individuals", wll[wh], "classes", gr[wll[wh]]  ,"and k", k, "\n")}
 }


#individuals 4 classes 1 and k 19 
#individuals 4 classes 1 and k 18 
#individuals 4 classes 1 and k 17 
#individuals 39 classes 1 and k 16 
#individuals 39 classes 1 and k 15 
#individuals 39 classes 1 and k 14 
#individuals 39 classes 1 and k 13 
#individuals 39 classes 1 and k 12 
#individuals 15 classes 1 and k 11 
#individuals 15 classes 1 and k 10 
#individuals 15 classes 1 and k 9 
#individuals 15 classes 1 and k 8 
#individuals 15 54 83 classes 1 0 -1 and k 7 
#individuals 15 54 83 classes 1 0 -1 and k 6 
#individuals 15 54 83 classes 1 0 -1 and k 5 
#individuals 3 15 83 classes 0 1 -1 and k 4 
#individuals 13 15 99 classes 0 1 -1 and k 3 
#individuals 13 15 80 classes 0 1 -1 and k 2 
#individuals 3 13 15 80 classes 0 0 1 -1 and k 1


gr=db[,2]
for (k in seq(21,1,-1)){
set.seed(2015)
situa=qik(xx,gr,k) #Q_i^k computing
  whs=situa==1
 if(any(whs)){cat("individuals,", which(whs)," classes", gr[which(whs)], "and k" ,k, "\n")}
 }

#individuals, 59 106  classes 1 1 and k 20 
#individuals, 4 59 106  classes 1 1 1 and k 19 
#individuals, 4 31 51 59 106  classes 1 1 1 1 1 and k 18 
#individuals, 4 31 51 59 106  classes 1 1 1 1 1 and k 17 
#individuals, 4 31 39 51 59 70 106  classes 1 1 1 1 1 1 1 and k 16 
#individuals, 4 31 39 51 57 59 70 98 106  classes 1 1 1 1 1 1 1 1 1 and k 15 
#individuals, 4 31 39 51 57 59 70 98 106  classes 1 1 1 1 1 1 1 1 1 and k 14 
#individuals, 4 6 31 39 51 57 59 70 98 106  classes 1 1 1 1 1 1 1 1 1 1 and k 13 
#individuals, 4 6 27 31 39 47 51 57 59 70 76 98 106  classes 1 1 1 1 1 1 1 1 1 1 1 1 1 and k 12 
#individuals, 4 6 15 27 31 39 47 51 57 59 70 76 79 98 106  classes 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 and k 11 
#individuals, 4 6 15 27 31 39 47 51 55 57 59 70 76 79 98 106  classes 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 and k 10 
#individuals, 4 6 15 27 31 39 41 47 51 55 57 59 70 76 79 98 106  classes 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 and k 9 
#individuals, 1 4 6 15 27 31 39 41 47 51 55 57 59 70 76 79 98 106  classes 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 and k 8 
#individuals, 1 4 6 15 27 31 39 41 47 51 54 55 57 59 70 76 79 83 98 106 109  classes 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 -1 1 1 1 and k 7 
#individuals, 1 4 6 15 27 31 39 41 42 47 51 54 55 57 59 70 76 79 83 98 106 109  classes 0 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 -1 1 1 1 and k 6 
#individuals, 1 4 6 15 27 31 34 39 41 42 45 47 51 54 55 57 59 70 76 79 83 84 98 106 109  classes 0 1 1 1 1 1 0 1 1 0 1 1 1 0 1 1 1 1 1 1 -1 0 1 1 1 and k 5 
#individuals, 1 3 4 6 7 15 27 31 34 39 41 42 44 45 47 51 54 55 57 59 70 76 79 83 84 90 98 106 109  classes 0 0 1 1 0 1 1 1 0 1 1 0 0 1 1 1 0 1 1 1 1 1 1 -1 0 1 1 1 1 and k 4 
#etc

###################################################################
#Computing importance measures
###################################################################

#removing nans for comparing with those by randomForest

#removing nans, boys

kk=db[-24,1:30]
kk=kk[-24,]

kk=kk[-62,]

db=kk
###end removing nans


####################################################
#mtry with less error rate with oob
####################################################
da=db[,4:30]
attach(db)
a=c()
for(i in 5:15){
my_cforest_control <- cforest_control(teststat = "max",
    testtype = "Teststatistic", mincriterion = 0, ntree = 500, mtry = i,
    replace = TRUE)
    set.seed(2015)
marboli=cforest(factor(Size_fit,ordered=TRUE)+Expert_evaluation~.,data=da,controls=my_cforest_control)

mq=matrix(unlist(predict(marboli,OOB=T)),ncol=4,byrow=TRUE)
cl=max.col(mq[,1:3])-2
a[i-4]= sum(cl==db[,2]) 
}
#> a  [1] 75 73 76 74 72 76 75 76 75 75 76 #i=7 biggest accuracy and lowest number, parsimony
#> 76/106= 0.7169811

i=7
my_cforest_control <- cforest_control(teststat = "max",
    testtype = "Teststatistic", mincriterion = 0, ntree = 500, mtry = i,
    replace = TRUE)
    set.seed(2015)

marbol=cforest(factor(Size_fit,ordered=TRUE)+Expert_evaluation~.,data=da,controls=my_cforest_control) 

#######################
#IPM with party, two responses
#######################

ntree=500
pupf=ipmparty(marbol,da,ntree) #IPM case-wise computing with OBB with party

#global IPM
pua=apply(pupf,2,mean) 

#rank global IPM  
 rank(pua)
 #28- rank(pua)
 #X10 X0030_mk    X1510    X1520    X1530 X3020_mk X3030_LR    X4020    X4030 X4040_mk X4080_LR    X4510    X5020 X5030_mk X5040_mk X6010_mk X6011_mk 
#  1        2       16       10       24       22       17       26       23       11       13        9       18       27       25        8        7 
#X6012_mk X6510_mk    X7520    X7525    X7540 X8010_LR X8030_LR X8040_LR X8520_LR X8550_LR 
#      19       20        4        3       14        5        6       15       21       12 

 
  
#IPM by classes

puac=matrix(0,nrow=3,ncol=dim(da)[2])
puac[1,]=apply(pupf[db[,2]==-1,],2,mean) 
puac[2,]=apply(pupf[db[,2]==0,],2,mean) 
puac[3,]=apply(pupf[db[,2]==1,],2,mean) 

colnames(puac)=colnames(da)
apply(puac,1,rank)



da1=as.data.frame(matrix(rep(0,27),nrow=1,ncol=27))
pupfn=ipmpartynew(marbol,da1,ntree)#IPM case-wise computing for new cases for party library
da2=rbind(da1,da1)
pupfn=ipmpartynew(marbol,da2,ntree)#IPM case-wise computing for new cases for party library
#colnames(pupfn)=colnames(da)

#######################################
#Comparing IPM with importance measures of RF Breiman. Only one response
########################################

#determining mtry for one response with RF Breiman, less OBB error rate
library(randomForest)
er=c()
for(i in 5:15){
set.seed(2015)
marbolr=randomForest(factor(Size_fit,ordered=TRUE)~.,data=da,ntree=500,mtry=i,importance=T,localImp=T,keep.inbag=T) 
er[i-4]=marbolr$err.rate[500,1]
}
#er: 0.2830189 0.3113208 0.2924528 0.2735849 0.3018868 0.2830189 0.2830189 0.2924528 0.3113208 0.2924528 0.3113208 #mtry=8

set.seed(2015)
marbolr=randomForest(factor(Size_fit,ordered=TRUE)~.,data=da,ntree=500,mtry=8,importance=T,localImp=T,keep.inbag=T) 

marbolr$importance
marbolr$localImportance

pupf=ipmrf(marbolr,da,ntree) #IPM case-wise computing with OBB with RandomForest

#global IPM
pua=apply(pupf,2,mean) 

#rank global IPM  
 rank(pua)
 
  
#IPM by classes

puac=matrix(0,nrow=3,ncol=dim(da)[2])
puac[1,]=apply(pupf[db[,2]==-1,],2,mean) 
puac[2,]=apply(pupf[db[,2]==0,],2,mean) 
puac[3,]=apply(pupf[db[,2]==1,],2,mean) 

colnames(puac)=colnames(da)
apply(puac,1,rank)



da1=as.data.frame(matrix(rep(0,27),nrow=1,ncol=27))
pupfn=ipmrfnew(marbolr,da1,ntree)#IPM case-wise computing for new cases for randomForest library
da2=rbind(da1,da1)
pupfn=ipmrfnew(marbolr,da2,ntree)#IPM case-wise computing for new cases for randomForest library


library(irr)
pim=rank(pua)
im=importance(marbolr)
ii=apply(im,2,rank)
 mmm=cbind(ii[,4],pim)
 kendall(mmm) #concordance global IPM and MeanDecreaseAccuracy  0.956 
 
 mmm=cbind(ii[,5],pim) 
 kendall(mmm) #concordance IPM and MeanDecreaseGini  0.995 
 
 mmm=cbind(ii[,5],ii[,4])
 kendall(mmm) #concordance MeanDecreaseAccuracy and MeanDecreaseGini  0.954 

piac=apply(puac,1,rank)
iiic=ii[,1:3]

cbind(iiic,piac) #per class importance MeanDecreaseAccuracy variable

mmm=cbind(ii[,1],piac[,1])
 kendall(mmm) #concordance class -1   0.955 


mmm=cbind(ii[,2],piac[,2])
 kendall(mmm) #concordance class 0    0.586 p-value = 0.248


mmm=cbind(ii[,3],piac[,3])
 kendall(mmm) #concordance class 1    0.905 

#28-cbind(ii,piac,pim)
#mmm=cbind(ii[,1],ii[,2])
#kendall(mmm) # 0.545 
#mmm=cbind(ii[,2],ii[,3])
#kendall(mmm)
#mmm=cbind(piac[,2],piac[,3])
# kendall(mmm) 0.982 
# mmm=cbind(piac[,1],piac[,2])
# kendall(mmm) 0.981
# mmm=cbind(piac[,1],piac[,3])
# kendall(mmm) 0.951





#################################################


###########################
#Classification boys Leave One Out individual
###########################

###############
#POLR
###############
#regression

library(MASS)

addterm_start2 <- function(model1, datos, ind)
{
  
  c <- model1$coefficients
  
  faltan <- datos[,c(-1,-which(names(datos)%in%names(c)))]
  #grepl(names(datos)[2], names(c))
  #TRUE%in%grepl(names(datos)[2], names(c))
  
  estan <- datos[,which(names(datos)%in%names(c))]
  
  min <- extractAIC(model1)[2]
  aic <- rep(NA, dim(faltan)[2])
  for(i in 1:dim(faltan)[2])
  {
    # print(names(faltan)[i])
    
    model2 <- polr(reformulate(termlabels = c(names(c), names(faltan)[i]), response = "Size_fit")
      ,data = datos, subset = ind,start = c(model1$coefficients, 0, model1$zeta))
    aic[i] <- extractAIC(model2)[2]
    # print(aic[i])
  }
  
  idmin <- which.min(aic)
  if(aic[idmin] < min)
  {
    return(names(faltan)[idmin])
  }
  else
  {
    return("")
  }
  
}

step_start <- function(model1, dat, ind)
{
  nuevoterm <- addterm_start2(model1, dat, ind)
  # print(nuevoterm)
  while(nuevoterm != "")
  {
    c <- model1$coefficients
    
    model2 <- polr(reformulate(termlabels = c(names(c), nuevoterm), response = "Size_fit")
                   ,data = dat, subset = ind,start = c(model1$coefficients, 0, model1$zeta))
    # print(extractAIC(model2))
    
    model1 <- model2
    nuevoterm <- addterm_start2(model1, dat, ind)
    # print(nuevoterm)
  }
  
  return(model1)
  
}


#size fit
ninos=unique(db[,1])

set.seed(2015)
preaire=c()
da=db[,c(2,4:30)]
da$Size_fit <- as.factor(da$Size_fit) 

for (i in 1:length(ninos)){

TR=db[,1] %in% ninos[-i]


nullmodel <- polr(as.factor(Size_fit) ~ 1, data = da, subset = TR)
stepmodel <- step_start(nullmodel, da,TR)

TE=db[,1] %in% ninos[i]


preaire=rbind(preaire,predict(stepmodel, da[TE,],type='probs'))

}


clre=max.col(preaire[,1:3])-2


sum(clre==db[,2]) #76 correct from 106


###########
#RF selecting mtry each time LOU
############
preai=c()


da=db[,-1]
attach(da)
mtv=c()

who=1:length(ninos)

for (i in who){

TR=db[,1] %in% ninos[-i]

aooo=c()

for (j in 5:15){
my_cforest_control <- cforest_control(teststat = "max",
    testtype = "Teststatistic", mincriterion = 0, ntree = 500, mtry = j,
    replace = TRUE)
set.seed(2015)
marboli=cforest(factor(Size_fit,ordered=TRUE)+Expert_evaluation~.,data=da,subset=TR,controls=my_cforest_control) 



qqi=predict(marboli,OOB=T)
pmi=matrix(unlist(qqi),ncol=4,byrow=TRUE)
cli=max.col(pmi[,1:3])-2

aooo[j-4]=sum(cli==db[TR,2])



}
mtv[i]=which.max(aooo)+4


my_cforest_control <- cforest_control(teststat = "max",
    testtype = "Teststatistic", mincriterion = 0, ntree = 500, mtry = which.max(aooo)+4,
    replace = TRUE)
set.seed(2015)
marboli=cforest(factor(Size_fit,ordered=TRUE)+Expert_evaluation~.,data=da,subset=TR,controls=my_cforest_control) 


TE=db[,1] %in% ninos[i]

qq=predict(marboli,db[TE,],OOB=FALSE)


preai=rbind(preai,matrix(unlist(qq),ncol=4,byrow=TRUE))


}


cl=max.col(preai[,1:3])-2

conjunto=db[,1] %in% ninos

sum(cl==db[conjunto,2])  #77 correct from 106

#mtv: 12  7 13  8 13 14 14 12 10 13  7  8 12  5 10  8 12 13 11  5 12 15 10 12  8 10  8 10  8  7 15  8  7 15 15 15 10  7 12 11 13 14 11 11  7 15  9  8  8 14
table(mtv)
# 5  7  8  9 10 11 12 13 14 15 
# 2  6  9  1  6  4  7  5  4  6 

# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#  5.00    8.00   11.00   10.64   13.00   15.00 

summary(abs(preai[,4]-db[conjunto,3])) #expert evaluation
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#0.06366 0.85260 1.31900 1.63600 1.98800 6.75600 

sum((abs(preai[,4]-db[conjunto,3]))>3) #13
sum((abs(preai[,4]-db[conjunto,3]))>4) # 4

#evl1=db[db[,2]==-1,3]
# evb1=db[db[,2]==1,3]
# ev0=db[db[,2]==0,3]
# table(evl1)
#evl1
#1 2 4 5 6 7 9 
#2 1 7 3 2 4 2 
# table(evb1)
#evb1
# 1  2  3  4  5  6  7  8  9 
# 5  4  4  4  5 12  5  3  1 
# table(ev0)
#ev0
# 6  7  8  9 10 
# 1  5 11 22  3 
 


############
#combining: ensemble
##########
#weighting

w=seq(0,1,0.01)


prc=c()
clw=c()
aw=c()

for (i in 1:length(w)){
prc=w[i]*preai[,1:3]+(1-w[i])*preaire
clw=max.col(prc[,1:3])-2
aw[i]=sum(clw==db[,2])
}

#aw: 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 76 77 77 77 77 77 77 77 77
# [51] 77 77 77 77 77 77 77 77 77 77 77 78 78 77 78 76 77 76 77 77 78 78 77 77 77 78 78 78 79 80 80 79 79 79 79 78 78 77 77 76 76 76 76 76 77 76 78 78 78 77
#[101] 77

eb=1-aw/106



are=clre==db[,2] #aciertos de regresion
 aa=cl==db[,2] #aciertos de forest 
mm=table(are,aa)
#       aa
#are     FALSE TRUE
#  FALSE    14   16
#  TRUE     15   61

yule= (mm[1,1]*mm[2,2]-mm[1,2]*mm[2,1])/(mm[1,1]*mm[2,2]+mm[1,2]*mm[2,1]) #0.5612431

#i=85 #obtained in comparative.R
#prc=w[i]*preai[,1:3]+(1-w[i])*preaire
#cbind(db[,1:3],prc,preai[,4], d[as.numeric(rownames(db)),42])

###########################
#Classification girls Leave One Out individual
###########################

###############
#POLR
###############
#regression



db=datos[d$SEXO=="M",]

#removing nan

kk=db[-22,1:30] #nan de expert evaluation
kk=kk[-22,]

kk=kk[-64,]

kk=kk[-64,]


db=kk
###end removing nans

#size fit
ninos=unique(db[,1])

set.seed(2015)
preaire=c()
da=db[,c(2,4:30)]
da$Size_fit <- as.factor(da$Size_fit) 

for (i in 1:length(ninos)){

TR=db[,1] %in% ninos[-i]


nullmodel <- polr(as.factor(Size_fit) ~ 1, data = da, subset = TR)
stepmodel <- step_start(nullmodel, da,TR)

TE=db[,1] %in% ninos[i]


preaire=rbind(preaire,predict(stepmodel, da[TE,],type='probs'))

}


clre=max.col(preaire[,1:3])-2


sum(clre==db[,2]) # 88 correct from 119 =0,73949579831932773109243697478992

###########
#RF selecting mtry each time LOU
############
preai=c()


da=db[,-1]
attach(da)
mtv=c()

who=1:length(ninos)

for (i in who){

TR=db[,1] %in% ninos[-i]

aooo=c()

for (j in 5:15){
my_cforest_control <- cforest_control(teststat = "max",
    testtype = "Teststatistic", mincriterion = 0, ntree = 500, mtry = j,
    replace = TRUE)
set.seed(2015)
marboli=cforest(factor(Size_fit,ordered=TRUE)+Expert_evaluation~.,data=da,subset=TR,controls=my_cforest_control) 


qqi=predict(marboli,OOB=T)
pmi=matrix(unlist(qqi),ncol=4,byrow=TRUE)
cli=max.col(pmi[,1:3])-2

aooo[j-4]=sum(cli==db[TR,2])

}
mtv[i]=which.max(aooo)+4


my_cforest_control <- cforest_control(teststat = "max",
    testtype = "Teststatistic", mincriterion = 0, ntree = 500, mtry = which.max(aooo)+4,
    replace = TRUE)
set.seed(2015)
marboli=cforest(factor(Size_fit,ordered=TRUE)+Expert_evaluation~.,data=da,subset=TR,controls=my_cforest_control) 


TE=db[,1] %in% ninos[i]

qq=predict(marboli,db[TE,],OOB=FALSE)

preai=rbind(preai,matrix(unlist(qq),ncol=4,byrow=TRUE))
}

cl=max.col(preai[,1:3])-2

conjunto=db[,1] %in% ninos

sum(cl==db[conjunto,2])  #86 correct from 119

#mtv: 12 15  5 13 11 13 11 14  8  5 14  8  5  9  5 14 15 12 14 12 11  5  7  7  5 10 11 11  6 10 13  7  7 11 11 14 15
#[38]  9  5 11 13 15  6  8  9 14  5 14  6  9  6  7 12 15  6 13 13 12
table(mtv)
# 5  6  7  8  9 10 11 12 13 14 15 
# 8  5  5  3  4  2  8  5  6  7  5
# 

# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 5.00    7.00   11.00   10.07   13.00   15.00 


summary(abs(preai[,4]-db[conjunto,3])) #expert evaluation
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#0.02992 0.54180 1.28000 1.43300 1.95900 5.35500

sum((abs(preai[,4]-db[conjunto,3]))>3) #9

 sum((abs(preai[,4]-db[conjunto,3]))>4) # 6


############
#combining: ensemble
##########
#weighting

w=seq(0,1,0.01)


prc=c()
clw=c()
aw=c()

for (i in 1:length(w)){
prc=w[i]*preai[,1:3]+(1-w[i])*preaire
clw=max.col(prc[,1:3])-2
aw[i]=sum(clw==db[,2])
}

#aw:  [1] 88 88 88 89 89 89 89 89 89 89 89 89 89 89 89 89 90 90 91 92 92 92 92 92 92 92 92 92 92 92 92 92 92 93 93 93 93 93 92 92 92 92 92 92 92 92 91 91 90 90
 #[51] 90 90 90 90 89 89 89 89 89 89 89 89 89 89 89 89 90 90 89 89 89 90 89 89 89 89 89 89 89 89 89 89 89 89 89 88 88 88 88 88 88 88 88 88 88 87 87 87 86 86
#[101] 86
# 93/119= 0.7815126



eg=1-aw/119



are=clre==db[,2] #aciertos de regresion
 aa=cl==db[,2] #aciertos de forest 
mm=table(are,aa)
#       aa
#are     FALSE TRUE
#  FALSE   23    8
#  TRUE    10   78


yule= (mm[1,1]*mm[2,2]-mm[1,2]*mm[2,1])/(mm[1,1]*mm[2,2]+mm[1,2]*mm[2,1]) # 0.9146211
plot(w,aw,type='l')


plot(w,eb,xlab="Weight",ylab="Error rate", type='n',ylim=c(0.215,0.285),cex.axis=1.5,cex.lab=1.5)

lines(w,eb)#solid black boys
lines(w,eg,lty=2,col=4)#blue dashed girl


# d[as.numeric(rownames(db)),c(1,42)]
#i=53 #obtained in comparative.R
#prc=w[i]*preai[,1:3]+(1-w[i])*preaire
#cbind(db[,1:3],prc,preai[,4], d[as.numeric(rownames(db)),42])



