#====================================================================================================================
#  Lõputöö "Logistilise regressiooni ja otsustuspuumeetodite kasutamine otsemüügi efektiivsuse suurendamiseks" 
#  praktilise osa tegemisel kasutatud R-kood
#       Hele-Liis Peedosk, mai 2017
#====================================================================================================================
#====================================================================================================================
#  1. Andmestiku sisselugemine
#====================================================================================================================
seeme=199231
algandmed= read.csv("fail.csv", header = TRUE, sep=";")
df.alg<-data.frame(algandmed)

set.seed(seeme+1)
sample <- sample.int(nrow(df.alg), floor(.70*nrow(df.alg)), replace = F)
df.train <- df.alg[sample, ]
df.test <- df.alg[-sample, ]

#====================================================================================================================
# 2.1 Logistilise regressioonimudeli eeltöö ja abifunktsioonid
#====================================================================================================================
# install.packages("car")
# install.packages("pROC")
# install.packages("brglm")
# install.packages("MASS")
# install.packages("gtools")
library(car)
library(pROC)
library(brglm)
library(MASS)
library(gtools)

# 1. Logistiliste mudelite hindamiseks kasutatavad andmed.
df.train.log <- df.alg[sample, ]
df.test.log <- df.alg[-sample, ]


# 2. Funktsioon, mis väljastab mudeli põhjal ROC-kurvi, leitud minimaalse tundlikkuse, mis on üle 95%, vastava klassifitseerimisväärtuse, 
#     vastava spetsiifilisuse ja eksimismaatriksi, lisaks kõik tundlikkused, spetsiifilisused.
rocc.log<-function(model, test, var, minsens=0.95){
  prob=predict(model,type=c("response"), newdata=test)
  g <- roc(var, predictor=prob)
  a<-cbind(g$thresholds, g$sensitivities, g$specificities)
  if (a[2,2]<minsens){
    sens=a[2,2]
    spec=a[2,3]
    threshold=a[2,1]
  }    else{
    ind=max(which(a[,2]>=minsens))
    sens=a[ind,2]
    spec=a[ind,3]
    threshold=a[ind,1]
  }
  plot(g, main=g$auc, print.thres = threshold)
  var[which(var==0)]=2
  klassid<-cbind(var, prob)
  klassid[which(klassid[,2]>threshold),2]=1
  klassid[which(klassid[,2]<=threshold),2]=2
  conf.mat<-table(klassid[,2], klassid[,1])
  var_list <- list("auc" = g$auc, "sens"=sens, "spec"=spec, "threshold"=threshold, "thresholds"= g$thresholds, "sensitivities"=g$sensitivities, "specificities"=g$specificities, "conf.mat"=conf.mat)
  return(var_list)
}


# 3. Funktsioon, mis eemaldab sammuviisi tunnuseid, millel on kõige kõrgem suurus VIF.
#     Allpool toodud funktsioon on modifitseeritud versioon järgneva blogi autori loodud funktisoonist.
#     Allikas: https://beckmw.wordpress.com/2013/02/05/collinearity-and-stepwise-vif-selection/
vif_func<-function(in_frame,val="leping", thresh=10,trace=T,...){
  if(class(in_frame) != 'data.frame') in_frame<-data.frame(in_frame)
  #get initial vif value for all comparisons of variables
  vif_init<-NULL
  var_names <- names(in_frame)
  regressors <- var_names[-which(var_names == val)]
  form <- paste(regressors, collapse = '+')
  form_in <- formula(paste(val, '~', form))
  viff<-vif(glm(form_in, family=binomial(link='logit'),data= in_frame))
  if (class(viff)=="matrix") vif_init<-rbind(vif_init, c(val, viff[,dim(viff)[2]]))
  else vif_init<-rbind(vif_init, c(val, viff))
  vif_max<-max(as.numeric(vif_init), na.rm = TRUE)
  if(vif_max < thresh){
    if(trace==T){ #print output of each iteration
      prmatrix(vif_init,collab=c('var','vif'),rowlab=rep('',nrow(vif_init)),quote=F)
      cat('\n')
      cat(paste('All variables have VIF < ', thresh,', max VIF ',round(vif_max,2), sep=''),'\n\n')
    }
    return(var_names)
  }    else{
    in_dat<-in_frame
    #backwards selection of explanatory variables, stops when all VIF values are below 'thresh'
    while(vif_max >= thresh){
      vif_vals<-NULL
      var_names <- names(in_dat)
      regressors <- var_names[-which(var_names == val)]
      form <- paste(regressors, collapse = '+')
      form_in <- formula(paste(val, '~', form))
      vif_add<-vif(glm(form_in, family=binomial(link='logit'),data= in_frame))
      if (class(vif_add)=="matrix") vif_vals<-c(val,  vif_add[,dim(vif_add)[2]])
      else vif_vals<-c(val,vif_add)
      vif_max<-max(as.numeric(vif_vals), na.rm = TRUE)
      if(vif_max<thresh) break
      if(trace==T){ #print output of each iteration
        prmatrix(vif_vals,collab=c('var','vif'),quote=F)
        cat('\n')
        cat('removed: ',names(vif_vals)[which(vif_vals==vif_max)],vif_max,'\n\n')
        flush.console()
      }
      in_dat<-in_dat[,!names(in_dat) %in% names(vif_vals)[which(vif_vals==vif_max)]]
    }
    var_names <- names(in_dat)
    regressors <- var_names[-which(var_names == val)]
    form <- paste(regressors, collapse = '+')
    print(formula(paste(val, '~', form)))
    return(formula(paste(val, '~', form)))
  }
}



#====================================================================================================================
# 2.2 Logistilised regressioonimudelid algsete väärtustega
#====================================================================================================================
# 1. Kasutades automaatset multikollineaarsete tunnuste eemaldamist, leitakse tunnused, mille põhjal edaspidi modelleerida. 
keep.dat<-vif_func(in_frame=df.train.log,val="leping", thresh=3,trace=TRUE)
lm1<-glm(keep.dat, family=binomial(link='logit'),data=df.train.log)

# 2. Tehakse sammuviisi valik multikollineaarsuse probleemita andmete tunnuste seast.
lm2.step<-step(lm1, direction = 'both')
lm2.step.sum<-summary(lm2.step)

# Sammuviisi regressiooni tulemus, kust on lisaks eemaldatud ebaolulised tunnused olulisusnivool 0,1.
lm2<-glm(leping ~ L4 + L6 + T1 + 
           L8 + kõne + I1 + I3 + K8 + KV7 + 
           K12 + M3 + M11 + 
           M18, family=binomial(link='logit'),data=df.train.log)
lm2.sum<-summary(lm2)
lm2.roc<-rocc.log(lm2, df.test.log, df.test.log$leping)
c(lm2.roc$threshold, lm2.roc$sens, lm2.roc$spec)


# 3. Tehakse sammuviisi valik, kasutades multikollineaarsuse probleemita andmete tunnuseid ja nende koosmõjusid tunnusega 'kõne'.
lm3.step<-step(lm1, scope = leping ~ L1 + L2 + L3 + L4 + 
                 L5 + L6 + L7 + T1 + kõne + L8 + T2 + T3 + L9 + 
                 L10 + I1 + I2 + I3 + K1 +  K4 + K9 + K10 + K11 + 
                 KV3 + KV7 + KV8 + KV10 + K12 + K13 + 
                 K14 + M3 + M4 + M7 + M8 + M11 + 
                 M14 + M16 + M17 +  M18 + MV1 + MV3 + 
                 MV7 + MV8 + MV10 +  MV12 + MV14 + 
                 MV17 + MV18 + MV20+
                 L1:kõne + L2:kõne + L3:kõne + L4:kõne + 
                 L5:kõne + L6:kõne + L7:kõne + T1:kõne + 
                 L8:kõne + T2:kõne + T3:kõne + L9:kõne + 
                 L10:kõne + I1:kõne + I2:kõne + I3:kõne + K1:kõne + 
                 K4:kõne + K9:kõne + K10:kõne + K11:kõne + 
                 KV3:kõne + KV7:kõne + KV8:kõne + 
                 KV10:kõne + K12:kõne + K13:kõne + 
                 K14:kõne + M3:kõne + M4:kõne + 
                 M7:kõne + M8:kõne + M11:kõne + 
                 M14:kõne + M16:kõne + M17:kõne + 
                 M18:kõne + MV1:kõne + MV3:kõne + 
                 MV7:kõne + MV8:kõne + MV4:kõne + 
                 MV12:kõne + MV14:kõne + 
                 MV17:kõne + MV18:kõne + MV20
               , direction = 'both')
lm3.step.sum<-summary(lm3.step)

# Sammuviisi valiku põhjal saadud mudel, kust on eemaldatud ka ebaolulised tunnused.
lm3<-glm(leping ~ L4 + L6 + T1 + L8 + T3 + I1 + I3 + K9 + 
           KV7 + K12 + M3 + M8 + M11 +  
           kõne + M8:kõne + T3:kõne + L6:kõne
         , family=binomial(link='logit'),data=df.train.log)
lm3.sum<-summary(lm3)
lm3.roc<-rocc.log(lm3, df.test.log, df.test.log$leping)
c(lm3.roc$threshold, lm3.roc$sens, lm3.roc$spec)


# 4 Logistiline regressioonimudel, kasutades ekspertvaliku põhjal valitud tunnuseid.      
lm5<-glm(leping ~ kõne + K12 + KV1_prots_kaart_sum + L6 + K8 + M6 
         + L8 +KV7 +KV6 +L10 + T3 +L4
         , family=binomial(link='logit'),data=df.train.log)
lm5.sum<-summary(lm5)
lm5.roc<-rocc.log(lm5, df.test.log, df.test.log$leping)
c(lm5.roc$threshold, lm5.roc$sens, lm5.roc$spec)


#====================================================================================================================
# 3.1 Logistilised regressioonimudelid grupeerimisega - eeltöö
#====================================================================================================================
# 1. Tehakse grupeeritud tunnustega andmestik.
df.train.gr<-df.train.log
df.test.gr<-df.test.log


# 2.1 Grupeerimise funktsioon mittenegatiivsete tunnuste jaoks, mis tagastab ka vastavad gruppide piirid:
#     1) jagab gruppidesse, kus eraldi on objektid, millel väärtus 0, ülejäänud jagatud tertsiilidesse;
#     2) kui valitud gruppidesse jääb vähem kui 100 objekti, teeb tunnuse binaarseks.   
group.pos<-function (x){
  q<-unique(quantile(x[x!=0], probs=seq(0,1,1/3)))  
  if (quantile(x[x!=0], probs=seq(0,1,1/3))[1]==quantile(x[x!=0], probs=seq(0,1,1/3))[2]) breax<-c(-Inf,0, q[1], q[2], Inf)
  else breax<-c(-Inf,0, q[2], q[3], Inf)
  ret<-cut(as.matrix(x), breaks = breax, labels<-FALSE,  right = TRUE)-1
  if (min(table(ret))<100) {
    ret[x==0]=0
    ret[x!=0]=1
    cuts<-rbind(c(0, 1), c(0,Inf))
  } else {cuts<-rbind(names(table(ret)), c(0, q[2], q[3], Inf))}
  ret_list=list(vals=ret, cuts=cuts )
  return(ret_list)
}

# 2.2 Grupeerimise funktsioon, mis tagastab ka vastavad gruppide piirid:
#     1) jagab viieks grupiks: 0-ga võrdsed, väiksemad kui positiivsete mediaan, suuremad kui positiivsete mediaan ning sama negatiivsete korral;
#     2) kui valitud gruppidesse jääb vähem kui 100 objekti, teeb tunnuse kolmetasemeliseks: negatiivsed, 0-id ja positiivsed.
group.neg<-function (x){
  q.neg<-unique(quantile(x[x<0], probs=seq(0,1,1/2)))
  q.pos<-unique(quantile(x[x>0], probs=seq(0,1,1/2)))
  breax<-c(-Inf, q.neg[2], -0.000001,0, q.pos[2], Inf)
  ret<-cut(as.matrix(x), breaks = breax, labels<-FALSE,  right = TRUE)-3
  if (min(table(ret))<100) {
    ret[x<0]=-1
    ret[x==0]=0
    ret[x>0]=1
    cuts<-rbind(c(-1, 0, 1), c(-0.000001, 0,Inf))
  } else{ cuts<-rbind(names(table(ret)), c(q.neg[2], -0.000001,0, q.pos[2], Inf))}
  ret_list=list(vals=ret, cuts=cuts)
  return(ret_list)
}

# 2.3 Grupeerimise funktsioon testandmestiku tunnuse 'x' jaoks, mis võtab lõikekohtadeks sisendina saadud treeningandmestiku gruppide piirid.
group.test<-function(x, cuts.m){
  n=dim(cuts.m)[2]
  vals<-rep(NA, length(x))
  vals[x<=as.numeric(cuts.m[2,1])]=cuts.m[1,1]
  for (j in 2:n) vals[x<=as.numeric(cuts.m[2,j]) & x>as.numeric(cuts.m[2,(j-1)])]=cuts.m[1,j]
  return(as.numeric(vals))
}


# 3. Vektor selliste tunnustega, mida tuleb grupeerida.
pos<-c("K1", "K2", "K3", "K5", "K6", "K7", "K10", "KV1", 
       "KV2", "KV3", "KV5", "KV6", "KV10", "K12", "K13", 
       "M1", "M2", "M5", "M6", "M11", "M12", "M13", 
       "M15", "M16", "MV1", "MV2", "MV5","MV11", "MV12", 
       "MV13", "MV15", "MV16", "MV6")

neg<-c("K4", "K11", "M3", "M7", "M14", 
       "M17", "MV3", "MV7", "MV14", "MV17")


# 4.1 Treeningandmestiku grupeerimine.
cuts <- vector(mode="list", length=length(pos)+length(neg))
names(cuts)<-c(pos, neg)
for (i in 1:length(colnames(df.train.gr))){
  if (is.element(colnames(df.train.gr)[i],pos)==TRUE){
    g<-group.pos(df.train.gr[i])
    df.train.gr[i]<-g$vals
    cuts[[colnames(df.train.gr)[i]]]<-g$cuts
  }
  if (is.element(colnames(df.train.gr)[i], neg)==TRUE){
    g<-group.neg(df.train.gr[i])
    df.train.gr[i]<-g$vals
    cuts[[colnames(df.train.gr)[i]]]<-g$cuts
  }
}


# 4.2 testandmestiku grupeerimine
for (i in 1:length(colnames(df.test.gr))){
  if (colnames(df.test.gr)[i] %in% names(cuts)){
    df.test.gr[i]<-group.test(df.test.gr[i],cuts[[which(names(cuts)==colnames(df.test.gr)[i])]])
  }
}


# 5.1 Vektor, mis sisaldab mittelineaarsete gruppidega tunnuseid.
mittelin<-c("KV3", "KV5", "KV6", "KV10", "M5", "M6", "M16")


# 5.2 Mittelineaarsete tunnuste faktoriks muutmine treening- ja testandmetikus.
for (i in 1:length(mittelin)){
  b<-mittelin[i]
  df.train.gr[[b]]<-as.factor(df.train.gr[[b]])
} 

for (i in 1:length(mittelin)){
  b<-mittelin[i]
  df.test.gr[[b]]<-as.factor(df.test.gr[[b]])
} 


# 6. Grupeerimisel tekkinud aliaste ehk kattuvate tunnuste eemaldamine nii test- kui treeningandmestikus.
lmg0<-glm(leping ~ ., family=binomial(link='logit'),data=df.train.gr)
alias(lmg0)
eemaldada<-c("MV5", "MV6", "MV15", "MV16"
             , "MV1", "M5", "KV10", "KV6")
for (i in 1:length(eemaldada)){
  b<-eemaldada[i]
  df.train.gr<-df.train.gr[,-which(colnames(df.train.gr)==b)]
} 
for (i in 1:length(eemaldada)){
  b<-eemaldada[i]
  df.test.gr<-df.test.gr[,-which(colnames(df.test.gr)==b)]
} 


#====================================================================================================================
# 3.2 Logistilised regressioonimudelid grupeerimisega - modelleerimine
#====================================================================================================================
# 1. Sammuviisi multikollineaarsete eemaldamine ning tulemuse põhjal mudeli tegemine.
keep.dat.gr<-vif_func(in_frame=df.train.gr,var="leping", thresh=3,trace=TRUE)
lmg1<-glm(keep.dat.gr, family=binomial(link='logit'),data=df.train.gr)


# 2. Sammuviisi AIC valik
lmg2.step<-step(lmg1, direction = 'both')
lmg2.step.sum<-summary(lmg2.step)
lmg2<-glm(leping ~ L4 + L6 + 
            T1 + L8 + kõne + I1 + 
            I3 + K10 + KV4
          , family=binomial(link='logit'),data=df.train.gr)
lmg2.sum<-summary(lmg2)
lmg2.roc<-rocc.log(lmg2, df.test.gr, df.test.gr$leping)
c(lmg2.roc$threshold, lmg2.roc$sens, lmg2.roc$spec)


# 4. Sammuviisi AIC, kasutades multikollineaarsuse probleemist vabasid tunnuseid ning arvestades ka nende koosmõjusid tunnusega 'kõne'.
lmg3.step<-step(lmg1, scope = leping ~ L1 + L2 + L3 + L4 + 
                  L5 + L6 + L7 + T1 + L8 + T2 + T3 + L9 + 
                  L10 + kõne + I1 + I2 + I3 + K1 +  K4 + K6 + K8 + K10 + 
                  K11 + KV3 + KV4 +  KV5 + KV7 + KV8 + 
                  K12 + K13 + K14 +  M3 + M4 + 
                  M6 + M7 + M8 + M11 + M12 + M14 + 
                  M16 + M17 + M18 + MV2 + MV3 + 
                  MV6 + MV7 + MV8 + MV10 + 
                  MV11 + MV12 + MV14 + MV17 + MV18 + MV19 + MV20 + 
                  L1:kõne + L2:kõne + L3:kõne + L4:kõne + 
                  L5:kõne + L6:kõne + L7:kõne + T1:kõne + 
                  L8:kõne + T2:kõne + T3:kõne + L9:kõne + 
                  L10:kõne + kõne:kõne + I1:kõne + I2:kõne + I3:kõne + K1:kõne + 
                  K4:kõne + K6:kõne + K8:kõne + K10:kõne + 
                  K11:kõne + KV3:kõne + KV4:kõne + 
                  KV5:kõne + KV7:kõne + KV8:kõne + 
                  K12:kõne + K13:kõne + K14:kõne + 
                  M3:kõne + M4:kõne + M6:kõne + M7:kõne + M8:kõne + 
                  M11:kõne + M12:kõne + M14:kõne + 
                  M16:kõne + M17:kõne + M18:kõne + 
                  MV2:kõne + MV3:kõne + MV4:kõne + MV7:kõne + 
                  MV8:kõne + MV10:kõne + MV11:kõne + MV12:kõne + MV14:kõne + 
                  MV17:kõne + MV18:kõne + MV19:kõne + 
                  MV20:kõne
                , direction = 'both')
lmg3.step.sum<-summary(lmg3.step)
lmg3<-glm(leping ~ L4 + L6 + T1 + L8 + T3 + kõne + 
            I1 + I3 + K4 + K10 +  KV4 + KV5 + 
            KV7 + M11 + MV14 +  kõne:M11 + kõne:KV5 + 
            kõne:K4 + kõne:MV14 + 
            T3:kõne + kõne:I3, family=binomial(link='logit'),data=df.train.gr)
lmg3.sum<-summary(lmg3)
lmg3.roc<-rocc.log(lmg3, df.test.gr, df.test.gr$leping)
c(lmg3.roc$threshold, lmg3.roc$sens, lmg3.roc$spec)


# 5. Grupeeritud tunnuste põhjal ekspertvalikul koostatud mudel.    
lmg5<-glm(leping ~ kõne +K10 + KV4 + T3 + T1
          + L4 + I3 + L6 +L10, family=binomial(link='logit'),data=df.train.gr)
lmg5.sum<-summary(lmg5)
lmg5.roc<-rocc.log(lmg5, df.test.gr, df.test.gr$leping)
c(llmg5.roc$threshold, lmg5.roc$sens, lmg5.roc$spec)


#====================================================================================================================
# 4.1 Logistilised regressioonimudelid kuupjuurteisendusega - eeltöö
#====================================================================================================================

# 1. Tehakse teisendatud tunnustega andmestik.
df.train.tr<-df.train.log
df.test.tr<-df.test.log

# 2.Vektor selliste tunnustega, millest võetakse kuupjuur.
kuubiks<-c("K1", "K2", "K3", "K4", "K5", 
           "K6", "K7", "K10", "K11", "KV1", 
           "KV2", "KV3", "KV5", "KV6", 
           "KV10", "K12", "K13",  
           "M1", "M2", "M3", "M5", 
           "M6", "M7", "M11", "M12", "M13", 
           "M14", "M15", "M16", "MV1", 
           "MV2", "MV3", "MV11", "MV12", 
           "MV13", "MV14")


# 3. Teeb kuupjuurteisenduse valitud tunnustele nii test- kui treeningandmestikul. 
for (i in 1:length(colnames(df.train.tr))){
  if (is.element(colnames(df.train.tr)[i], kuubiks)==TRUE){
    df.train.tr[i]<-sign(df.train[i])*abs((df.train[i]))^(1/3)
  }
}

for (i in 1:length(colnames(df.test.tr))){
  if (is.element(colnames(df.test.tr)[i], kuubiks)==TRUE){
    df.test.tr[i]<-sign(df.test[i])*abs((df.test[i]))^(1/3)
  }
}


#====================================================================================================================
# 4.2 Logistilised regressioonimudelid kuupjuurteisendusega - modelleerimine
#====================================================================================================================
# 1. sammuviisi multikollineaarsete eemaldamine
keep.dat.tr<-vif_func(in_frame=df.train.tr,var="leping", thresh=3,trace=TRUE)
lmt1<-glm(keep.dat.tr, family=binomial(link='logit'),data=df.train.tr)

# 2. Sammuviisi AIC, kasutades multikollineaarsuse probleemist vabu andmeid.
lmt2.step<-step(lmt1, direction = 'both')
lmt2.step.sum<-summary(lmt2.step)

# Sammuviisi valiku põhjal saadud mudel, kust eemaldati ka ebaolulised tunnused.
lmt2<-glm(leping ~ L4 +  L6 + 
            T1 + L8 + kõne + I1 + I3 + 
            K10 + KV7 +
            M18
          , family=binomial(link='logit'),data=df.train.tr)
lmt2.sum<-summary(lmt2)
lmt2.roc<-rocc.log(lmt2, df.test.tr, df.test.tr$leping)
c(lmt2.roc$threshold, lmt2.roc$sens, lmt2.roc$spec)


# 3. Sammuviisi AIC, kus lisaks on valitud tunnuste juurde arvestatud ka koosmõjusid tunnusega 'kõne'.
lmt3.step<-step(lmt1, scope = leping ~ L1 + L2 + L3 + L4 + 
                  L5 + L6 + L7 + T1 + L8 + T2 + T3 + L9 + 
                  L10 + kõne + I1 + I2 + I3 + K4 + K8 + K10 + K11 + KV7 + 
                  KV8 + KV10 + K13 + K14 + M3 + M6+ M7 + M8 + M14 + M17 + 
                  M18 + MV1 + MV3 + MV7 + MV8 + MV10 + 
                  MV12 + MV14 + MV17 + MV18 + MV20 +
                  L1:kõne  + L2:kõne  + L3:kõne  + L4:kõne  + 
                  L5:kõne  + L6:kõne  + L7:kõne  + T1:kõne  + 
                  L8:kõne  + T2:kõne  + T3:kõne  + L9:kõne  + 
                  L10:kõne  + kõne:kõne  + I1:kõne  + I2:kõne  + I3:kõne  + K4:kõne  + 
                  K8:kõne  + K10:kõne  + K11:kõne  + KV7:kõne  + 
                  KV8:kõne  + KV10:kõne  + 
                  K13:kõne  + K14:kõne  + M3:kõne  + 
                  M14:kõne  + M7:kõne  + 
                  M8:kõne  + M14:kõne  + M17:kõne  + 
                  M18:kõne  + MV1:kõne  + MV3:kõne  + 
                  MV7:kõne  + MV8:kõne  + MV10  + 
                  MV12:kõne  + MV14:kõne  + 
                  MV17:kõne  + MV18:kõne  + MV20:kõne 
                , direction = 'both')
lmt3.step.sum<-summary(lmt3.step)

# Koosmõjudega sammuviisi AIC tulemus, millest eemaldatud mitteolulised tunnused.
lmt3<-glm(leping ~ L4 + L6 + T1 + 
            L8 + T3 + kõne + I1 + I3 + 
            K10 + KV7 + K13 + M7 + 
            M8 + M18 + MV18 + 
            kõne:M8 + T3:kõne + kõne:K4 + 
            L6:kõne, family=binomial(link='logit'),data=df.train.tr)
lmt3.sum<-summary(lmt3)
lmt3.roc<-rocc.log(lmt3, df.test.tr, df.test.tr$leping)
c(lmt3.roc$threshold, lmt3.roc$sens, lmt3.roc$spec)


# 4. Ekspertvaliku põhjal valitud teisendatud tunnustega log. regressioonimudel. 
lmt5<-glm(leping ~ kõne + KV3 + KV4 + KV7 + I1
          + K10 + L4 +L6 + L8 + M19
          , family=binomial(link='logit'),data=df.train.tr)
(lmt5.sum<-summary(lmt5))
lmt5.roc<-rocc.log(lmt5, df.test.tr, df.test.tr$leping)
c(lmt5.roc$threshold, lmt5.roc$sens, lmt5.roc$spec)


#====================================================================================================================
# 6. Puude eeltöö
#====================================================================================================================
# 1. Kasutatavad paketid.
# install.packages("rpart")
# install.packages("rpartScore")
# install.packages("randomForest")
# install.packages("caret")
# install.packages("CORElearn")
# install.packages("e1071")
# install.packages("plyr")
library(rpart)
library(rpartScore)
library(randomForest)
library(caret)
require(e1071)
library(pROC)
library(plyr)
library(DMwR)


# 1. Puumeetodite jaoks tehakse test- ja treeningandmestiud.
df.train.puu<- df.alg[sample, ]
df.test.puu<- df.alg[-sample, ]

# 2. Kuna puudel võib binaarne 0/1 tekkitada mõningaid probleeme, muudetakse uuritav tunnus sobivaks.
df.train.puu$leping<-as.factor(df.train.puu$leping)
a<-factor(df.train.puu$leping, labels=c("ei", "jah"))
df.train.puu$leping <- factor(a, levels = c("jah", "ei"))

df.test.puu$leping<-as.factor(df.test.puu$leping)
a<-factor(df.test.puu$leping, labels=c("ei", "jah"))
df.test.puu$leping <- factor(a, levels = c("jah", "ei"))  

# 3. Faktortunnused muudetakse faktoriteks.
for (i in 2:81){
  if (any(c(2:8, 10, 12, 15, 17, 81)==i)){
    print(i)
    df.train.puu[, i]<-as.factor(df.train.puu[, i])
  }
}

for (i in 2:81){
  if (any(c(2:8, 10, 12, 15, 17, 81)==i)){
    print(i)
    df.test.puu[, i]<-as.factor(df.test.puu[, i])
  }
}


# 4.1 Funktsioon alavaliku tegemiseks. Alusena kasutati paketi 'caret' koodi, lisati juurde positiivse klassi osakaalu 
#     muutmine parameetriga 'rate'. Võetakse tasakaalustatud otsutusmetsa algoritmi põhiselt bootstrap-valim positiivsetest ning 
#     lih. juh. tagasipanekuga valik negatiivsetest.

downsample.tr.bs=function (x, y, rate=0.5, list=FALSE) {
  xc <- class(x)
  if (!is.data.frame(x)) 
    x <- as.data.frame(x)
  if (!is.factor(y)) {
    y<-as.factor(y)  }
  minClass <- min(table(y))
  maxClass<- round(((1-rate)*minClass)/rate)
  x$.outcome <- y
  mincl<-if ((table(y))[1]<(table(y))[2]) names(table(y))[1] else names(table(y))[2] 
  maxcl<-if ((table(y))[1]>(table(y))[2]) names(table(y))[1] else names(table(y))[2] 
  ind<-c(sample(which(x$.outcome==mincl), min(table(y)), replace=TRUE), sample(which(x$.outcome==maxcl), min(maxClass, max(table(y))), replace=TRUE))
  x<-x[ind,]
  y <- x$.outcome
  x <- x[, !(colnames(x) %in% c("y", ".outcome")), drop = FALSE]
  if (list) {
    if (xc[1] == "matrix") 
      x <- as.matrix(x)
    out <- list(x = x, y = y)
  }
  else {
    out <- cbind(y,x)
    colnames(out)[1] <- "leping"
    rownames(out)<-seq(1:length(out[,1]))
  }
  out}


# 4.2 Funktsioon alavaliku tegemiseks. Analoogne eelnevale, kuid siin võetakse valimisse kõik positiivsed objektid ning 
#     lih. juh. tagasipanekuta valik negatiivsetest.
downsample.tr=function (x, y, rate=0.5, list=FALSE) {
  xc <- class(x)
  if (!is.data.frame(x)) 
    x <- as.data.frame(x)
  if (!is.factor(y)) {
    y<-as.factor(y)  }
  minClass <- min(table(y))
  #print(minClass)
  maxClass<- round(((1-rate)*minClass)/rate)
  #print(maxClass)
  x$.outcome <- y
  mincl<-if ((table(y))[1]<(table(y))[2]) names(table(y))[1] else names(table(y))[2] 
  maxcl<-if ((table(y))[1]>(table(y))[2]) names(table(y))[1] else names(table(y))[2] 
  ind<-c(which(x$.outcome==mincl), sample(which(x$.outcome==maxcl), min(maxClass, max(table(y))), replace=FALSE))
  x<-x[ind,]
  y <- x$.outcome
  x <- x[, !(colnames(x) %in% c("y", ".outcome")), drop = FALSE]
  if (list) {
    if (xc[1] == "matrix") 
      x <- as.matrix(x)
    out <- list(x = x, y = y)
  }  else {
    out <- cbind(y,x)
    colnames(out)[1] <- "leping"
    rownames(out)<-seq(1:length(out[,1]))
  }
  out}  


# 5.  Funktsioon, mis väljastab puumudeli põhjal ROC-kurvi, leitud minimaalse tundlikkuse, mis on üle 95%, vastava klassifitseerimisväärtuse, 
#     vastava spetsiifilisuse ja eksimismaatriksi, lisaks kõik tundlikkused, spetsiifilisused.
rocc.tree<-function(model, test, var, minsens=0.95){
  prob <- predict(model, newdata = test, type = "prob")
  g <- roc(response = var,
           predictor = prob[, "jah"],
           levels = rev(levels(var)))
  a<-cbind(g$thresholds, g$sensitivities, g$specificities)
  if (a[2,2]<minsens){
    sens=a[2,2]
    spec=a[2,3]
    threshold=a[2,1]
  }  else{
    ind=max(which(a[,2]>=minsens))
    sens=a[ind,2]
    spec=a[ind,3]
    threshold=a[ind,1]
  }
  plot(g, print.thres=threshold)
  klassid<-cbind(var, prob[,1])
  klassid[which(klassid[,2]>threshold),2]=1
  klassid[which(klassid[,2]<=threshold),2]=2
  conf.mat<-table(klassid[,2], klassid[,1])
  var_list <- list("auc" = g$auc, "sens"=sens, "spec"=spec, "threshold"=threshold, "thresholds"= g$thresholds, "sensitivities"=g$sensitivities, "specificities"=g$specificities, "conf.mat"=conf.mat)
  return(var_list)
}


#====================================================================================================================
# 7.1 Üksik klassifitseerimispuu
#====================================================================================================================
attach(df.train.puu)

# 1. Klassifitseerimispuu treenimine.   
ctrl.ct <- trainControl(method = "none", seeds=seeme)
ct <- train(leping ~ ., data = df.train.puu,
            method = "rpart1SE",
            trControl =  ctrl.ct)
# Protseduuri tulemusena et tekkinud otsustuspuule ühtegi lehte. Ainult juurtipp mis prognoosib alati negatiivsesse klassi.

#====================================================================================================================
# 7.2 Taasvalikuga klassifitseerimispuud
#====================================================================================================================

# 1. Tehakse esmalt erinevate osakaaludega alavalikud ning SMOTE taasvalik treeningandmestiku pealt.
set.seed(seeme+1)
df.train.puu.down2 <- downsample.tr(x = df.train.puu[, -1],
                                    y = df.train.puu$leping, rate=0.2, list=FALSE)
set.seed(seeme+2)
df.train.puu.down3 <- downsample.tr(x = df.train.puu[, -1],
                                    y = df.train.puu$leping, rate=0.3, list=FALSE)
set.seed(seeme+3)
df.train.puu.down4 <- downsample.tr(x = df.train.puu[, -1],
                                    y = df.train.puu$leping, rate=0.4, list=FALSE)
set.seed(seeme+4)
df.train.puu.down5 <- downsample.tr(x = df.train.puu[, -1],
                                    y = df.train.puu$leping, list=FALSE)
set.seed(seeme)
df.train.smote<- SMOTE(leping ~ ., data = df.train.puu)   

# 2. Klassifitseerimispuu, kasutades valimit positiivse klassi osakaaluga 0,2.
detach(df.train.puu)
attach(df.train.puu.down2)
ctrl.ct.down<- trainControl(method = "none", seeds=seeme, classProbs = TRUE,
                            summaryFunction = twoClassSummary.cust)
ct.down2 <- train(leping ~ ., data = df.train.puu.down2,
                  method = "rpart1SE",
                  trControl =  ctrl.ct.down)
# Tulemusena tekkis ainult juurtipust koosnev puu.


# 3. Klassifitseerimispuu, kasutades valimit positiivse klassi osakaaluga 0,3.
detach(df.train.puu.down2)
attach(df.train.puu.down3)
ct.down3 <- train(leping ~ ., data = df.train.puu.down3,
                  method = "rpart1SE",
                  trControl =  ctrl.ct.down)
# Tulemusena tekkis ainult juurtipust koosnev puu.


# 4. Klassifitseerimispuu, kasutades valimit positiivse klassi osakaaluga 0,4.  
detach(df.train.puu.down3)
attach(df.train.puu.down4)
ct.down4 <- train(leping ~ ., data = df.train.puu.down4,
                  method = "rpart1SE",
                  trControl =  ctrl.ct.down)
# Puu graafiline kuju
plot(ct.down4$finalModel)
text(ct.down4$finalModel)

ctd4.roc<-rocc.tree(ct.down4, df.test.puu, df.test.puu$leping)
c(ctd4.roc$threshold, ctd4.roc$sens, ctd4.roc$spec)


# 5. Klassifitseerimispuu, kasutades valimit positiivse klassi osakaaluga 0,5.  
detach(df.train.puu.down4)
attach(df.train.puu.down5)
ct.down5 <- train(leping ~ ., data = df.train.puu.down5,
                  method = "rpart1SE",
                  trControl =  ctrl.ct.down)

# Puu graafiline kuju
plot(ct.down5$finalModel)
text(ct.down5$finalModel)

ctd5.roc<-rocc.tree(ct.down5, df.test.puu, df.test.puu$leping)
c(ctd5.roc$threshold, ctd5.roc$sens, ctd5.roc$spec)
detach(df.train.puu.down5)


# 6. Klassifitseerimispuu, kasutades taasvaliku meetodit SMOTE.     
attach(df.train.smote)
ct.smote<- train(leping ~ ., data = df.train.smote,
                 method = "rpart1SE",
                 trControl =  ctrl.ct.down)

ctsm.roc<-rocc.tree(ct.smote, df.test.puu, df.test.puu$leping)
c(ctsm.roc$threshold, ctsm.roc$sens, ctsm.roc$spec)
detach(df.train.smote)


#====================================================================================================================
# 8.1 Otsustusmets
#====================================================================================================================
# 1. Tehakse raamistik tulemuste kogumise jaoks.
tulem<-data.frame()
tulem<-rbind(tulem, rep(NA, 8))
colnames(tulem)<- c("Tüüp", "ntree", "mtry", "m", "nodesize", "cutoff", "sens", "spec") 
tulem<-tulem[-1,]

ntrees<-c(500, 1000)                           # puude arv
Nss<-c(1,5)                                    # minimaalne lehtede arv
ms<-c(9, 40, 70)                               # kirjeldavate tunnuste arv, mida kaasata

attach(df.train.puu)

# 2. Treeni otsustusmetsad erinevate parameetrite kombinatsioonide põhjal.    
ctrl.rf <- trainControl(method = "oob", classProbs = TRUE, seeds=seeds)

for (ntree in ntrees){
  for (Ns in Nss)  {
    for (m in ms){
      rf <- train(x=df.train.puu[,-1], y=df.train.puu$leping, data = df.train.puu,
                  method = "rf",
                  ntree = ntree,
                  tuneGrid=expand.grid(mtry=c(m)),
                  trControl =  ctrl.rf,
                  nodesize=Ns)
      rf.roc<-rocc.tree(rf, df.test.puu, df.test.puu$leping)
      print(c("rf", ntree, 1, m, Ns, rf.roc$threshold, rf.roc$sens, rf.roc$spec))
      de<- c("rf", ntree, 1, m, Ns, rf.roc$threshold, rf.roc$sens, rf.roc$spec)
      names(de)=c("Tüüp", "ntree", "mtry", "m", "nodesize", "cutoff", "sens", "spec") 
      tulem<-rbind(tulem,de)
    }
  }
}


# 3. Mudeli stabiliseerumise graafik.
layout(matrix(c(1,2),nrow=1),
       width=c(4,1)) 
par(mar=c(5,4,4,0),col.lab = "white") 
plot(rf$finalModel, main="", log="y")
par(col.lab = "black")
title(xlab = "Puude arv", ylab = "Veamäär") 
par(mar=c(5,0,4,2)) 
plot(c(0,1),type="n", axes=F, xlab="", ylab="")
legend("top", colnames(rf$finalModel$err.rate),col=1:4,cex=0.8,fill=1:4)
layout(matrix(c(1,1),nrow=1),
       width=c(4,1)) 

#====================================================================================================================
# 8.2 Taasvalikuga otsustusmetsad
#====================================================================================================================
# 1. Parameetrite ja seemnete tegemine. 
M<-40                           # otsib 40 juhusliku arvu seast parimat arvu m, mille puhul Kappa suurim
Tr<-1500                        # puude arv
Ns<-1                           # minimaalne lehtede arv
seeds <- vector(mode = "list", length = 1+1)
set.seed(seeme)
seeds[[1]] <- sample.int(1000, M)
seeds[[2]]<-seeme


# 2. Otsustusmetsa konstrueerimine, kasutades iga puu korral erinevaid SMOTE taasvaliku meetodil saadud valimeid. 
ctrl.rfs <- trainControl(method = "oob", sampling="smote",  classProbs = TRUE, seeds=seeds, search="random")
rfs <- train(leping ~ ., data = df.train.puu,
             method = "rf",
             ntree = Tr,
             tuneLength=M,
             metric = "Kappa",
             trControl =  ctrl.rfs)

plot(rfs$finalModel)
rfs.roc<-rocc.tree(rfs, df.test.puu, df.test.puu$leping)
c("rfs", Tr, M, rfs$bestTune$mtry, Ns, rfs.roc$threshold, rfs.roc$sens, rfs.roc$spec)
(de<- c("rfs", Tr, M, rfs$bestTune$mtry, Ns, rfs.roc$threshold, rfs.roc$sens, rfs.roc$spec))
names(de)=c("Tüüp", "ntree", "mtry", "m", "nodesize", "cutoff", "sens", "spec") 
tulem<-rbind(tulem,de)


# 3. Otsustusmetsa konstrueerimine, kasutades iga puu korral erinevaid 20% alavaliku meetodil saadud valimeid. 
ctrl.rfd2 <- trainControl(method = "oob", sampling=function(x,y) downsample.tr.bs(x,y, rate=0.2, list=TRUE),  classProbs = TRUE, seeds=seeds, search="random")
rfd2 <- train(leping ~ ., data = df.train.puu,
              method = "rf",
              ntree = Tr,
              tuneLength=M,
              metric = "Kappa",
              trControl =  ctrl.rfd2)

plot(rfd2$finalModel)
rfd2.roc<-rocc.tree(rfd2, df.test.puu, df.test.puu$leping)
c("rfd2", Tr, M, rfd2$bestTune$mtry, Ns, rfd2.roc$threshold, rfd2.roc$sens, rfd2.roc$spec)
de<-c("rfd2", Tr, M, rfd2$bestTune$mtry, Ns, rfd2.roc$threshold, rfd2.roc$sens, rfd2.roc$spec)
names(de)=c("Tüüp", "ntree", "mtry", "m", "nodesize", "cutoff", "sens", "spec") 
tulem<-rbind(tulem,de)


# 3. Otsustusmetsa konstrueerimine, kasutades iga puu korral erinevaid 30% alavaliku meetodil saadud valimeid. 
ctrl.rfd3 <-  trainControl(method = "oob", sampling=function(x,y) downsample.tr.bs(x,y, rate=0.3, list=TRUE),  classProbs = TRUE, seeds=seeds, search="random")
rfd3 <- train(leping ~ ., data = df.train.puu,
              method = "rf",
              ntree = Tr,
              tuneLength=M,
              metric = "Kappa",
              trControl =  ctrl.rfd3)

plot(rfd3$finalModel)
rfd3.roc<-rocc.tree(rfd3, df.test.puu, df.test.puu$leping)
c("rfd3", Tr, M, rfd3$bestTune$mtry, Ns, rfd3.roc$threshold, rfd3.roc$sens, rfd3.roc$spec)
de<-c("rfd3", Tr, M, rfd3$bestTune$mtry, Ns, rfd3.roc$threshold, rfd3.roc$sens, rfd3.roc$spec)
names(de)=c("Tüüp", "ntree", "mtry", "m", "nodesize", "cutoff", "sens", "spec") 
tulem<-rbind(tulem,de)



# 4.  Otsustusmetsa konstrueerimine, kasutades iga puu korral erinevaid 40% alavaliku meetodil saadud valimeid. 
ctrl.rfd4 <- trainControl(method = "oob", sampling=function(x,y) downsample.tr.bs(x,y, rate=0.4, list=TRUE),  classProbs = TRUE, seeds=seeds, search="random")
rfd4 <- train(leping ~ ., data = df.train.puu,
              method = "rf",
              ntree = Tr,
              tuneLength=M,
              metric = "Kappa",
              trControl =  ctrl.rfd4)

plot(rfd4$finalModel)
rfd4.roc<-rocc.tree(rfd4, df.test.puu, df.test.puu$leping)
c("rfd4", Tr, M, rfd4$bestTune$mtry, Ns, rfd4.roc$threshold, rfd4.roc$sens, rfd4.roc$spec)
de<-c("rfd4", Tr, M, rfd4$bestTune$mtry, Ns, rfd4.roc$threshold, rfd4.roc$sens, rfd4.roc$spec)
names(de)=c("Tüüp", "ntree", "mtry", "m", "nodesize", "cutoff", "sens", "spec") 
tulem<-rbind(tulem,de)


# 5.  Otsustusmetsa konstrueerimine, kasutades iga puu korral erinevaid 50% alavaliku meetodil saadud valimeid. 
ctrl.rfd5 <- trainControl(method = "oob", sampling=function(x,y) downsample.tr(x,y, rate=0.5, list=TRUE),  classProbs = TRUE, seeds=seeds, search="random")
rfd5 <- train(leping ~ ., data = df.train.puu,
              method = "rf",
              ntree = Tr,
              tuneLength=M,
              metric = "Kappa",
              trControl =  ctrl.rfd5)

plot(rfd5$finalModel)
rfd5.roc<-rocc.tree(rfd5, df.test.puu, df.test.puu$leping)
cbind("rfd5", Tr, M, rfd5$bestTune, Ns, rfd5.roc$threshold, rfd5.roc$sens, rfd5.roc$spec)
names(de)=c("Tüüp", "ntree", "mtry", "m", "nodesize", "cutoff", "sens", "spec") 
tulem<-rbind(tulem,de)


# 6. 40% alavaliku meetodiga otsustusmetsa stabiliseerumise graafik.   
layout(matrix(c(1,2),nrow=1),
       width=c(4,1)) 
par(mar=c(5,4,4,0),col.lab = "white") 
plot(rfdx$finalModel, main="", log="y")
par(col.lab = "black")
title(xlab = "Puude arv", ylab = "Veamäär") 
par(mar=c(5,0,4,2)) 
plot(c(0,1),type="n", axes=F, xlab="", ylab="")
legend("top", colnames(rfdx$finalModel$err.rate),col=1:4,cex=0.8,fill=1:4)


detach(df.train.puu)
#====================================================================================================================
# 9.1 Agregeerimata tunnustega puumudelite eeltöö
#====================================================================================================================

# 1. Algandmestiku lugemine.
set.seed(seeme)
algandmed_teh= read.csv("fail2.csv", header = TRUE, sep=";")
df.alg.v<-data.frame(algandmed_teh)

set.seed(seeme+1)
sample <- sample.int(nrow(df.alg.v), floor(.70*nrow(df.alg.v)), replace = F)
df.train.v <- df.alg.v[sample, ]
df.test.v <- df.alg.v[-sample, ]

# 2. Tunnuste sobivaks tegemine.  
# uuritav tunnus
df.train.v$leping<-as.factor(df.train.v$leping)
a<-factor(df.train.v$leping, labels=c("ei", "jah"))
df.train.v$leping <- factor(a, levels = c("jah", "ei"))

df.test.v$leping<-as.factor(df.test.v$leping)
a<-factor(df.test.v$leping, labels=c("ei", "jah"))
df.test.v$leping <- factor(a, levels = c("jah", "ei"))  

# kvalitatiivsed kirjeldavad tunnused 
for (i in 2:16){
  if (any(c(2:8, 10, 12, 14, 16)==i)){
    print(i)
    df.train.v[, i]<-as.factor(df.train.v[, i])
  }
}

for (i in 2:16){
  if (any(c(2:8, 10, 12, 14, 16)==i)){
    print(i)
    df.test.v[, i]<-as.factor(df.test.v[, i])
  }
}

#====================================================================================================================
# 9.2 Agregeerimata tunnustega puumudelite konstreueerimine
#====================================================================================================================
# 1. Klassifitseerimispuu tasakaalustatud valimiga.
set.seed(seeme)
df.train.v.d5<- downsample.tr(x = df.train.v[, -1],
                              y = df.train.v$leping, list=FALSE) 
attach(df.train.v.d5)
ctd5.v<- train(leping ~ ., data = df.train.v.d5,
               method = "rpart1SE",
               trControl =  ctrl.ct.down)

ctd5.v.roc<-rocc.tree(ctd5.v, df.test.v, df.test.v$leping)
c(ctd5.v.roc$threshold, ctd5.v.roc$sens, ctd5.v.roc$spec)

plot(ctd5.v$finalModel)
text(ctd5.v$finalModel)
detach(df.train.v.d5)


# 2. Klassifitseerimispuu konstrueerimine, kasutades 40% alavaliku meetodil saadud valimit. 
set.seed(seeme)
df.train.v.d4<- downsample.tr(x = df.train.v[, -1],
                              y = df.train.v$leping, rate=0.4, list=FALSE) 
attach(df.train.v.d4)
ctd4.v <- train(leping ~ ., data = df.train.v.d4,
                method = "rpart1SE",
                trControl =  ctrl.ct.down)
ctd4.v.roc<-rocc.tree(ctd4.v, df.train.v.d4, df.train.v.d4$leping)
c(ctd4.v.roc$auc, ctd4.v.roc$threshold, ctd4.v.roc$sens,ctd4.v.roc$spec)

plot(ctd4.v$finalModel)
text(ctd4.v$finalModel)

detach(df.train.v.4)



# 3.  Otsustusmetsa konstrueerimine, kasutades iga puu korral erinevaid 40% alavaliku meetodil saadud valimeid. 
ctrl.rfd4 <- trainControl(method = "oob", sampling=function(x,y) downsample.tr.bs(x,y, rate=0.4, list=TRUE),  classProbs = TRUE, seeds=seeds, search="random")
rfd4.v <- train(leping ~ ., data = df.train.v,
                method = "rf",
                ntree = Tr,
                tuneLength=M,
                metric = "Kappa",
                trControl =  ctrl.rfd4)

plot(rfd4.v$finalModel)
rfd4.roc.v<-rocc.tree(rfd4.v, df.test.v, df.test.v$leping)
c("rfd4.v", Tr, M, rfd4.v$bestTune$mtry, Ns, rfd4.roc.v$threshold, rfd4.roc.v$sens, rfd4.roc.v$spec)
de<-c("rfd4.v", Tr, M, rfd4.v$bestTune$mtry, Ns, rfd4.roc.v$threshold, rfd4.roc.v$sens, rfd4.roc.v$spec)
names(de)=c("Tüüp", "ntree", "mtry", "m", "nodesize", "cutoff", "sens", "spec") 
tulem<-rbind(tulem,de)