R/idiomclassification_mk_pf3.R - private/idioms - Gitiles

 library(caret)
 library(tidyverse)
 library(DMwR)
 library(randomForest)
 library(FSelector)
 library(MLmetrics)
 # library(randomForestExplainer)
 # may need to: options(expressions = 5e5) to avoid stackoverflow for installing package

 set.seed(42)

 # Test

 ngramfile2<-"gold03_anno_ml_synfeat_nstopw" # 2nd dataset
 ngramfile1 <-"goldstandard01_anno_ml_synfeat_nstop1" # 1st dataset

 setwd(dirname(rstudioapi::getSourceEditorContext()$path))
 stopwords <- readLines(con = "../data/stopwords.txt",encoding="UTF-8")
 oringramme1 <- read.csv(paste("../data/",ngramfile1,".csv",sep=""), header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)
 oringramme2 <- read.csv(paste("../data/",ngramfile2,".csv",sep=""), header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)


 syfeaturenames <- read.csv("../data/syfeatures.tsv", header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)

 featurenames <- read.csv("../data/features.tsv", header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)

 deleteStopwords =  function(wl, stopwords = NULL) {
   wl[!(wl %in% stopwords)]
 }

 oringramme1 <- oringramme1 %>%
   mutate(CO_IDIOM = ifelse(is.na(CO_IDIOM),0,CO_IDIOM)) %>%  # treat NAs as 0
   filter(CO_IDIOM < 2)   # just two classes: 0 no idiom, 1 idiom

 oringramme2 <- oringramme2 %>%
   mutate(CO_IDIOM = ifelse(is.na(CO_IDIOM),0,CO_IDIOM)) %>%  # treat NAs as 0
   filter(CO_IDIOM < 2)   # just two classes: 0 no idiom, 1 idiom

 # Reduce number of classes, treat null values, add NSTOPW, change names for SY features and rename all features
 # new featurenames.

 ngramme1 <- oringramme1 %>%
   add_column(NSTOPW = sapply(oringramme1$tokens,function(x) length(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords)))) %>%
   add_column(nstokens = sapply(oringramme1$tokens, function(x) paste(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords),collapse=" "))) %>%
   # select(-matches("CO_TOKEN.*"), -tokens) %>%
   add_column(CO_GRAM = sapply(oringramme1$tokens, function(x) length(unlist(strsplit(x," ")))))%>%
   select(-matches("CO_TOKEN.*")) %>% # keep tokens for interpretability
   mutate(across(matches(".rank.*"), ~ replace_na(.x, 1000))) %>%
   mutate(across(c("dice", "lfmd", "llr", "ld", "pmi"), ~ replace_na(.x, min(.x) - 1))) %>%
   rename_at(syfeaturenames$innames, ~ syfeaturenames[syfeaturenames$innames==.x,]$synames ) %>%
   rename_at(featurenames$oldnames, ~ featurenames[featurenames$oldnames==.x,]$newnames ) %>%
   mutate(across(everything(), ~ replace_na(.x, 0))) %>%
   mutate(CO_IDIOM = as.factor(if_else(CO_IDIOM == 1, "idiom", "no_idiom"))) # just two classes: 0 no idiom, 1 idiom

 ngramme2 <- oringramme2 %>%
   add_column(NSTOPW = sapply(oringramme2$tokens,function(x) length(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords)))) %>%
   add_column(nstokens = sapply(oringramme2$tokens, function(x) paste(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords),collapse=" "))) %>%
   # select(-matches("CO_TOKEN.*"), -tokens) %>%
   select(-matches("CO_TOKEN.*")) %>% # keep tokens for interpretability
   mutate(across(matches(".rank.*"), ~ replace_na(.x, 1000))) %>%
   mutate(across(c("dice", "lfmd", "llr", "ld", "pmi"), ~ replace_na(.x, min(.x) - 1))) %>%
   rename_at(syfeaturenames$innames, ~ syfeaturenames[syfeaturenames$innames==.x,]$synames ) %>%
   rename_at(featurenames$oldnames, ~ featurenames[featurenames$oldnames==.x,]$newnames ) %>%
   mutate(across(everything(), ~ replace_na(.x, 0))) %>%
   mutate(CO_IDIOM = as.factor(if_else(CO_IDIOM == 1, "idiom", "no_idiom"))) # just two classes: 0 no idiom, 1 idiom

 # combine

 ngramme1<-subset(ngramme1,select=-c(CO_SONGS))
 ngramme <- rbind(ngramme1,ngramme2[colnames(ngramme1)])
 # ngramme<-ngramme1

 ngramme<-ngramme1%>% distinct(nstokens,.keep_all=T)
 ngramme1<-ngramme1%>% distinct(nstokens,.keep_all=T)

 # Optional
 write.table(ngramme,file=paste("../data/","combined_noduplicates.tsv",sep=""), sep = "\t", quote=F)

 write.table(ngramme1,file=paste("../data/","dataset1_noduplicates.tsv",sep=""), sep = "\t", quote=F)


 # featuresets

 o_vars <- c("O_C2_N", "O_C2_SGT", "O_DEREKO", "O_GRAM", "O_NSTOPW")
 o_vars_1 <- c("O_DEREKO", "O_GRAM", "O_NSTOPW")
 o_vars_2 <- c("O_C2_N", "O_C2_SGT")
 co_vars <- c("CO_VEC","CO_VEC_LEX")
 sy_c1_vars <- c("SY_C1_C_L","SY_C1_C_R","SY_C1_DICE","SY_C1_LD","SY_C1_LDAF","SY_C1_LL","SY_C1_MI","SY_C1_MI_L","SY_C1_MI_R","SY_C1_MI2","SY_C1_MI3","SY_C1_NMI")
 sy_c1_vars_1 <- c("SY_C1_LD","SY_C1_LDAF","SY_C1_LL","SY_C1_MI","SY_C1_MI3")
 # sy_c1_vars_1 <- c("SY_C1_LD","SY_C1_LL","SY_C1_MI","SY_C1_MI3")
 sy_c2_vars <- c("SY_C2_EXP","SY_C2_G","SY_C2_K","SY_C2_LD","SY_C2_LL","SY_C2_LMI","SY_C2_MI","SY_C2_MI3","SY_C2_T","SY_C2_Z")
 sy_c2_vars_1 <- c("SY_C2_G","SY_C2_LD","SY_C2_LL","SY_C2_MI","SY_C2_MI3")
 sy_w_vars <- c("SY_W_AVG","SY_W_CON","SY_W_MAX","SY_W_NSUM","SY_W_NSUM_AF")
 # sy_w_vars_1 <- c("SY_W_AVG","SY_W_CON","SY_W_MAX","SY_W_NSUM")
 sy_r_vars <- c("SY_C1_R","SY_W_R1","SY_W_R2","SY_R_D")

 all_vars <- c(o_vars,co_vars,sy_c1_vars,sy_c2_vars,sy_w_vars,sy_r_vars)
 all_vars_1 <- c(o_vars, co_vars, sy_c1_vars_1, sy_c2_vars_1, sy_w_vars, sy_r_vars)


 # formulae for training and testing rf

 all_fml <- as.formula(paste("CO_IDIOM ~ ", paste(all_vars, collapse= "+")))
 all_fml_1 <- as.formula(paste("CO_IDIOM ~ ", paste(all_vars_1, collapse= "+")))

 sy_c1_fml_1 <- as.formula(paste("CO_IDIOM ~ ", paste(sy_c1_vars_1, collapse= "+")))
 sy_c2_fml_1 <- as.formula(paste("CO_IDIOM ~ ", paste(sy_c2_vars_1, collapse= "+")))
 sy_w_fml <- as.formula(paste("CO_IDIOM ~ ", paste(sy_c2_vars_1, collapse= "+")))
 sy_r_fml <- as.formula(paste("CO_IDIOM ~ ", paste(sy_r_vars, collapse= "+")))
 co_fml <- as.formula(paste("CO_IDIOM ~ ", paste(co_vars, collapse= "+")))
 o_fml_1 <- as.formula(paste("CO_IDIOM ~ ", paste(o_vars_1, collapse= "+")))
 o_fml <- as.formula(paste("CO_IDIOM ~ ", paste(o_vars, collapse= "+")))
 o_fml_2 <- as.formula(paste("CO_IDIOM ~ ", paste(o_vars_2, collapse= "+")))

 # Simple train/test split

 set.seed(111)
 trainRows <- sample(nrow(ngramme), nrow(ngramme)*0.8, replace = FALSE)
 train <- ngramme[trainRows,]
 test <- ngramme[setdiff(1:nrow(ngramme),trainRows),]

 cat("Random Forest\n")

 rf_classifier = randomForest(all_fml, train, importance=TRUE)

 # only SY features
 # rf_classifier = randomForest(fmlasy, train, importance=TRUE)

 prediction_for_table <- predict(rf_classifier, test %>% select(-CO_IDIOM))

 res <- confusionMatrix(prediction_for_table, test$CO_IDIOM, positive= "idiom",mode="everything")
 print(res)
 collected_results <- bind_cols("rf" = res$byClass)

 # Sensitivity is recall of class 1
 # Pos Pred Value is precision
 varImpPlot(rf_classifier)

 cat("Random Forest with cutoff\n")
 prediction_for_table <- predict(rf_classifier,test %>% select(-CO_IDIOM), cutoff = c(0.3, 0.7))
 res <- confusionMatrix(prediction_for_table,test$CO_IDIOM, positive = "idiom",mode="everything")
 collected_results <- bind_cols(collected_results, "rf with cutoff" = res$byClass)
 print(res)

 cat("With SMOTE resampled training data\n")
 smoted.data <- SMOTE(all_fml, subset(train, select = c("CO_IDIOM", all_vars)), perc.over = 1200, perc.under = 100)
 rf_classifier = randomForest(all_fml, smoted.data, importance=TRUE)
 prediction_for_table <- predict(rf_classifier,test %>% select(-CO_IDIOM))
 res <- confusionMatrix(prediction_for_table,test$CO_IDIOM, positive = "idiom",mode="everything")
 collected_results <- bind_cols(collected_results, "rf with SMOTE" = res$byClass)
 print(res)

 cat("With SMOTE and cutoff\n")
 prediction_for_table <- predict(rf_classifier,test %>% select(-CO_IDIOM), cutoff = c(0.3, 0.7))
 res <- confusionMatrix(prediction_for_table,test$CO_IDIOM, positive = "idiom",mode="everything")
 collected_results <- bind_cols(collected_results, "rf with SMOTE and cutoff" = res$byClass)
 print(res)

 collected_results <- collected_results %>%
   round(3) %>%
   add_column(measure = names(res$byClass)) %>%
   column_to_rownames("measure")

 View(collected_results)

 # Analysing tradeoff between Fscore, Recall, Precision for various cutoffs
 # full range from precision almost 100% to recall almost 100%
 rf_classifier = randomForest(all_fml, train, importance=TRUE)
 cvalues<-tibble()
 for (c in c(seq(from=0.4,to=0.99,by=0.025),0.999)) {
   prediction_for_table <- predict(rf_classifier, test %>% select(-CO_IDIOM), cutoff = c(1-c, c))
   conf<-confusionMatrix(prediction_for_table, test$CO_IDIOM, positive = "idiom")
   cvalues <-bind_rows(cvalues, c(cutoff=c, conf$byClass))
 }
 cvalues %>%
   select(c("cutoff", "Recall", "Precision", "F1", "Specificity", "Balanced Accuracy")) %>%
   pivot_longer(!cutoff, names_to=c("measure")) %>%
   ggplot(aes(cutoff, value, colour=measure)) + geom_line(size=1) +
   scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
   scale_y_continuous(breaks = scales::pretty_breaks(n = 10))

 # cross validation

 featuresets<-data.frame(read.csv("../data/featuresets.tsv", header = FALSE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE))

 collected_results<-tibble()
 for (i in c(1:nrow(featuresets))) {
   set.seed(325)
   train_control <- trainControl(method="repeatedcv", number=5, repeats=5, verboseIter=T,classProbs = TRUE,
                                 savePredictions = "final")
   vars<-unlist(strsplit(featuresets[i,2],","))
   fml <- as.formula(paste("CO_IDIOM ~ ", paste(vars, collapse= "+")))
   rf <- train(fml,
               data=ngramme1,
               method='rf',
               tuneGrid = data.frame(.mtry = floor(sqrt(length(vars)))),
               cutoff=c(0.3,0.7),
               trControl=train_control)
   res<-confusionMatrix(rf$pred$pred, rf$pred$obs,mode="everything")
   collected_results <-bind_rows(collected_results, c(features=featuresets[i,1], res$byClass))
 }

 collected_results1<-collected_results[c(2:16),c(1,6,7,8,12)]
 collected_results1<-collected_results1%>%
   mutate(across(everything(), ~ replace_na(.x, 0))) %>%
   mutate(across(c(2:5), ~ round(as.numeric(.x),digits=3)))

 # Analysing tradeoff between Fscore, Recall, Precision for various cutoffs
 # full range from precision almost 100% to recall almost 100%


 vars<-unlist(strsplit(featuresets[2,2],","))
 fml <- as.formula(paste("CO_IDIOM ~ ", paste(vars, collapse= "+")))
 cvalues<-tibble()
 for (c in c(seq(from=0.4,to=0.99,by=0.025),0.999)) {
   set.seed(325)
   train_control <- trainControl(method="repeatedcv", number=5, repeats=5, verboseIter=T,classProbs = TRUE,
                                 savePredictions = "final")
   rf <- train(fml,
               data=ngramme1,
               method='rf',
               tuneGrid = data.frame(.mtry = floor(sqrt(length(vars)))),
               cutoff=c(1-c,c),
               trControl=train_control)
   res<-confusionMatrix(rf$pred$pred, rf$pred$obs,mode="everything")
   cvalues <-bind_rows(cvalues, c(cutoff=c, res$byClass))
 }
 cvalues %>%
   select(c("cutoff", "Recall", "Precision", "F1", "Specificity", "Balanced Accuracy")) %>%
   pivot_longer(!cutoff, names_to=c("measure")) %>%
   ggplot(aes(cutoff, value, colour=measure)) + geom_line(size=1) +
   scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
   scale_y_continuous(breaks = scales::pretty_breaks(n = 10))

 cvaluesall<-cvalues


 # Using estimates by random forest on entire dataset

 library(randomForest)
 rf_classifier_full = randomForest(all_fml_1, data=ngramme1, importance=TRUE, cutoff=c(0.3,0.7))
 rf_classifier_full
 # class.error is 1 - recall
 varImpPlot(rf_classifier_full)

 # Feature ranking

 # rf features as table

 # correlated features seem to split their rankings

 rfranks<-importance(rf_classifier_full)[,3:4]

 # ttest

 idioms<-ngramme1 %>% filter(CO_IDIOM == "idiom")
 nonidioms<-ngramme1 %>% filter(CO_IDIOM != "idiom")

 m1<-mean(unlist(nonidioms["SY_C1_LD"]),na.rm=T)

 nonidioms1<-nonidioms%>%filter(SY_C1_LD > m1)

 mean(unlist(idioms["SY_C1_LD"]),na.rm=T)

 idioms1<-ngramme1 %>% filter(CO_IDIOM == "idiom")
 idioms2<-ngramme2 %>% filter(CO_IDIOM == "idiom")

 ttestPvalues<-sapply(all_vars_1,
                      function(sel) t.test(idioms[sel],nonidioms[sel])$p.value)

 t.test(idioms["SY_C1_LL"],nonidioms["SY_C1_LL"])

 ttestSignificance<-sapply(all_vars_1,
                      function(sel) {
                        p<-t.test(idioms[sel],nonidioms[sel])$p.value
                        if (p < 0.001) {
                          return("***")
                        }
                        if (p < 0.01) {
                          return("**")
                        }
                        if (p < 0.05) {
                          return ("*")
                        }
                        return(" ")
                      })


 ttestSignificance1<-sapply(all_vars_1,
                           function(sel) {
                             p<-t.test(idioms[sel],nonidioms1[sel])$p.value
                             if (p < 0.001) {
                               return("***")
                             }
                             if (p < 0.01) {
                               return("**")
                             }
                             if (p < 0.05) {
                               return ("*")
                             }
                             return(" ")
                           })


 # information gain
 # multiply by 1000 to avoid undersized bins
 # features are ranked individually no matter their correlation
 igain<-information.gain(all_fml_1, data=ngramme1%>%mutate_at(all_vars_1, ~ . * 1000),unit="log2")

 # difference between means (positive or negative?)

 diffMeans<-sapply(all_vars_1,function(sel) mean(unlist(idioms[sel]),na.rm=T)-mean(unlist(nonidioms[sel]),na.rm=T))

 diffMeansSign<-sapply(all_vars_1,function(sel) ifelse(mean(unlist(idioms[sel]),na.rm=T)-mean(unlist(nonidioms[sel]),na.rm=T)>0,"+","-"))

 diffMeansSign1<-sapply(all_vars_1,function(sel) ifelse(mean(unlist(idioms[sel]),na.rm=T)-mean(unlist(nonidioms1[sel]),na.rm=T)>0,"+","-"))


 featurenames[,c("newnames","explanation")]

 featureRanks<-cbind(rownames(rfranks),rfranks,igain,ttestPvalues,ttestSignificance,diffMeans,diffMeansSign,ttestSignificance1,diffMeansSign1)
 colnames(featureRanks)[1]<-"newnames"

 featureRanks<-merge(featureRanks,featurenames[,c("newnames","explanation")],by="newnames")


 #randomForestExplainer::explain_forest(rf_classifier )

 # averate estimates and feature ranks over 10 runs

 errrate<-0
 conf<-matrix(0,2,3)
 featureRanks<-matrix(0,4,length(vars))
 for (i in 1:10) {
   rfc =randomForest(fmla, data=ngramme, importance=TRUE)
   #rfc =randomForest(fmla, data=ngramme, importance=TRUE, cutoff=c(0.2, 0.8))
   errrate<-errrate+rfc$err.rate[100,1]
   conf<-conf+rfc$confusion
   featureRanks<-featureRanks+
     cbind(importance(rfc)[,3:4],
           information.gain(fmla, data=ngramme%>%mutate_at(vars, ~ . * 1000),unit="log2"),
           sapply(vars,
                  function(sel) t.test(idioms[sel],nonidioms[sel])$p.value))
   print(errrate/i)
   conf1<-round(
     rbind(
       cbind(conf[,1:2]/i,(1-conf[,3]/i)*100),
       c(100*diag(conf[,1:2])/colSums(conf[,1:2]),NA),
       c(rowSums(conf[,1:2]/i),NA)),digits=2)
   colnames(conf1)<-c("1","0","rec")
   rownames(conf1)<-c("1","0","prec","sum")
   print(conf1)
 }
 featureRanks<-featureRanks/10
 colnames(featureRanks)<-c("MeanDecreaseAccuracy","MeanDecreaseGini","InformationGain","Ttest")
	library(caret)
	library(tidyverse)
	library(DMwR)
	library(randomForest)
	library(FSelector)
	library(MLmetrics)
	# library(randomForestExplainer)
	# may need to: options(expressions = 5e5) to avoid stackoverflow for installing package

	set.seed(42)

	# Test

	ngramfile2<-"gold03_anno_ml_synfeat_nstopw" # 2nd dataset
	ngramfile1 <-"goldstandard01_anno_ml_synfeat_nstop1" # 1st dataset

	setwd(dirname(rstudioapi::getSourceEditorContext()$path))
	stopwords <- readLines(con = "../data/stopwords.txt",encoding="UTF-8")
	oringramme1 <- read.csv(paste("../data/",ngramfile1,".csv",sep=""), header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)
	oringramme2 <- read.csv(paste("../data/",ngramfile2,".csv",sep=""), header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)


	syfeaturenames <- read.csv("../data/syfeatures.tsv", header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)

	featurenames <- read.csv("../data/features.tsv", header = TRUE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE)

	deleteStopwords = function(wl, stopwords = NULL) {
	wl[!(wl %in% stopwords)]
	}

	oringramme1 <- oringramme1 %>%
	mutate(CO_IDIOM = ifelse(is.na(CO_IDIOM),0,CO_IDIOM)) %>% # treat NAs as 0
	filter(CO_IDIOM < 2) # just two classes: 0 no idiom, 1 idiom

	oringramme2 <- oringramme2 %>%
	mutate(CO_IDIOM = ifelse(is.na(CO_IDIOM),0,CO_IDIOM)) %>% # treat NAs as 0
	filter(CO_IDIOM < 2) # just two classes: 0 no idiom, 1 idiom

	# Reduce number of classes, treat null values, add NSTOPW, change names for SY features and rename all features
	# new featurenames.

	ngramme1 <- oringramme1 %>%
	add_column(NSTOPW = sapply(oringramme1$tokens,function(x) length(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords)))) %>%
	add_column(nstokens = sapply(oringramme1$tokens, function(x) paste(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords),collapse=" "))) %>%
	# select(-matches("CO_TOKEN.*"), -tokens) %>%
	add_column(CO_GRAM = sapply(oringramme1$tokens, function(x) length(unlist(strsplit(x," ")))))%>%
	select(-matches("CO_TOKEN.*")) %>% # keep tokens for interpretability
	mutate(across(matches(".rank.*"), ~ replace_na(.x, 1000))) %>%
	mutate(across(c("dice", "lfmd", "llr", "ld", "pmi"), ~ replace_na(.x, min(.x) - 1))) %>%
	rename_at(syfeaturenames$innames, ~ syfeaturenames[syfeaturenames$innames==.x,]$synames ) %>%
	rename_at(featurenames$oldnames, ~ featurenames[featurenames$oldnames==.x,]$newnames ) %>%
	mutate(across(everything(), ~ replace_na(.x, 0))) %>%
	mutate(CO_IDIOM = as.factor(if_else(CO_IDIOM == 1, "idiom", "no_idiom"))) # just two classes: 0 no idiom, 1 idiom

	ngramme2 <- oringramme2 %>%
	add_column(NSTOPW = sapply(oringramme2$tokens,function(x) length(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords)))) %>%
	add_column(nstokens = sapply(oringramme2$tokens, function(x) paste(deleteStopwords(tolower(unlist(strsplit(x," "))),stopwords),collapse=" "))) %>%
	# select(-matches("CO_TOKEN.*"), -tokens) %>%
	select(-matches("CO_TOKEN.*")) %>% # keep tokens for interpretability
	mutate(across(matches(".rank.*"), ~ replace_na(.x, 1000))) %>%
	mutate(across(c("dice", "lfmd", "llr", "ld", "pmi"), ~ replace_na(.x, min(.x) - 1))) %>%
	rename_at(syfeaturenames$innames, ~ syfeaturenames[syfeaturenames$innames==.x,]$synames ) %>%
	rename_at(featurenames$oldnames, ~ featurenames[featurenames$oldnames==.x,]$newnames ) %>%
	mutate(across(everything(), ~ replace_na(.x, 0))) %>%
	mutate(CO_IDIOM = as.factor(if_else(CO_IDIOM == 1, "idiom", "no_idiom"))) # just two classes: 0 no idiom, 1 idiom

	# combine

	ngramme1<-subset(ngramme1,select=-c(CO_SONGS))
	ngramme <- rbind(ngramme1,ngramme2[colnames(ngramme1)])
	# ngramme<-ngramme1

	ngramme<-ngramme1%>% distinct(nstokens,.keep_all=T)
	ngramme1<-ngramme1%>% distinct(nstokens,.keep_all=T)

	# Optional
	write.table(ngramme,file=paste("../data/","combined_noduplicates.tsv",sep=""), sep = "\t", quote=F)

	write.table(ngramme1,file=paste("../data/","dataset1_noduplicates.tsv",sep=""), sep = "\t", quote=F)



	# featuresets

	o_vars <- c("O_C2_N", "O_C2_SGT", "O_DEREKO", "O_GRAM", "O_NSTOPW")
	o_vars_1 <- c("O_DEREKO", "O_GRAM", "O_NSTOPW")
	o_vars_2 <- c("O_C2_N", "O_C2_SGT")
	co_vars <- c("CO_VEC","CO_VEC_LEX")
	sy_c1_vars <- c("SY_C1_C_L","SY_C1_C_R","SY_C1_DICE","SY_C1_LD","SY_C1_LDAF","SY_C1_LL","SY_C1_MI","SY_C1_MI_L","SY_C1_MI_R","SY_C1_MI2","SY_C1_MI3","SY_C1_NMI")
	sy_c1_vars_1 <- c("SY_C1_LD","SY_C1_LDAF","SY_C1_LL","SY_C1_MI","SY_C1_MI3")
	# sy_c1_vars_1 <- c("SY_C1_LD","SY_C1_LL","SY_C1_MI","SY_C1_MI3")
	sy_c2_vars <- c("SY_C2_EXP","SY_C2_G","SY_C2_K","SY_C2_LD","SY_C2_LL","SY_C2_LMI","SY_C2_MI","SY_C2_MI3","SY_C2_T","SY_C2_Z")
	sy_c2_vars_1 <- c("SY_C2_G","SY_C2_LD","SY_C2_LL","SY_C2_MI","SY_C2_MI3")
	sy_w_vars <- c("SY_W_AVG","SY_W_CON","SY_W_MAX","SY_W_NSUM","SY_W_NSUM_AF")
	# sy_w_vars_1 <- c("SY_W_AVG","SY_W_CON","SY_W_MAX","SY_W_NSUM")
	sy_r_vars <- c("SY_C1_R","SY_W_R1","SY_W_R2","SY_R_D")

	all_vars <- c(o_vars,co_vars,sy_c1_vars,sy_c2_vars,sy_w_vars,sy_r_vars)
	all_vars_1 <- c(o_vars, co_vars, sy_c1_vars_1, sy_c2_vars_1, sy_w_vars, sy_r_vars)


	# formulae for training and testing rf

	all_fml <- as.formula(paste("CO_IDIOM ~ ", paste(all_vars, collapse= "+")))
	all_fml_1 <- as.formula(paste("CO_IDIOM ~ ", paste(all_vars_1, collapse= "+")))

	sy_c1_fml_1 <- as.formula(paste("CO_IDIOM ~ ", paste(sy_c1_vars_1, collapse= "+")))
	sy_c2_fml_1 <- as.formula(paste("CO_IDIOM ~ ", paste(sy_c2_vars_1, collapse= "+")))
	sy_w_fml <- as.formula(paste("CO_IDIOM ~ ", paste(sy_c2_vars_1, collapse= "+")))
	sy_r_fml <- as.formula(paste("CO_IDIOM ~ ", paste(sy_r_vars, collapse= "+")))
	co_fml <- as.formula(paste("CO_IDIOM ~ ", paste(co_vars, collapse= "+")))
	o_fml_1 <- as.formula(paste("CO_IDIOM ~ ", paste(o_vars_1, collapse= "+")))
	o_fml <- as.formula(paste("CO_IDIOM ~ ", paste(o_vars, collapse= "+")))
	o_fml_2 <- as.formula(paste("CO_IDIOM ~ ", paste(o_vars_2, collapse= "+")))

	# Simple train/test split

	set.seed(111)
	trainRows <- sample(nrow(ngramme), nrow(ngramme)*0.8, replace = FALSE)
	train <- ngramme[trainRows,]
	test <- ngramme[setdiff(1:nrow(ngramme),trainRows),]

	cat("Random Forest\n")

	rf_classifier = randomForest(all_fml, train, importance=TRUE)

	# only SY features
	# rf_classifier = randomForest(fmlasy, train, importance=TRUE)

	prediction_for_table <- predict(rf_classifier, test %>% select(-CO_IDIOM))

	res <- confusionMatrix(prediction_for_table, test$CO_IDIOM, positive= "idiom",mode="everything")
	print(res)
	collected_results <- bind_cols("rf" = res$byClass)

	# Sensitivity is recall of class 1
	# Pos Pred Value is precision
	varImpPlot(rf_classifier)

	cat("Random Forest with cutoff\n")
	prediction_for_table <- predict(rf_classifier,test %>% select(-CO_IDIOM), cutoff = c(0.3, 0.7))
	res <- confusionMatrix(prediction_for_table,test$CO_IDIOM, positive = "idiom",mode="everything")
	collected_results <- bind_cols(collected_results, "rf with cutoff" = res$byClass)
	print(res)

	cat("With SMOTE resampled training data\n")
	smoted.data <- SMOTE(all_fml, subset(train, select = c("CO_IDIOM", all_vars)), perc.over = 1200, perc.under = 100)
	rf_classifier = randomForest(all_fml, smoted.data, importance=TRUE)
	prediction_for_table <- predict(rf_classifier,test %>% select(-CO_IDIOM))
	res <- confusionMatrix(prediction_for_table,test$CO_IDIOM, positive = "idiom",mode="everything")
	collected_results <- bind_cols(collected_results, "rf with SMOTE" = res$byClass)
	print(res)

	cat("With SMOTE and cutoff\n")
	prediction_for_table <- predict(rf_classifier,test %>% select(-CO_IDIOM), cutoff = c(0.3, 0.7))
	res <- confusionMatrix(prediction_for_table,test$CO_IDIOM, positive = "idiom",mode="everything")
	collected_results <- bind_cols(collected_results, "rf with SMOTE and cutoff" = res$byClass)
	print(res)

	collected_results <- collected_results %>%
	round(3) %>%
	add_column(measure = names(res$byClass)) %>%
	column_to_rownames("measure")

	View(collected_results)

	# Analysing tradeoff between Fscore, Recall, Precision for various cutoffs
	# full range from precision almost 100% to recall almost 100%
	rf_classifier = randomForest(all_fml, train, importance=TRUE)
	cvalues<-tibble()
	for (c in c(seq(from=0.4,to=0.99,by=0.025),0.999)) {
	prediction_for_table <- predict(rf_classifier, test %>% select(-CO_IDIOM), cutoff = c(1-c, c))
	conf<-confusionMatrix(prediction_for_table, test$CO_IDIOM, positive = "idiom")
	cvalues <-bind_rows(cvalues, c(cutoff=c, conf$byClass))
	}
	cvalues %>%
	select(c("cutoff", "Recall", "Precision", "F1", "Specificity", "Balanced Accuracy")) %>%
	pivot_longer(!cutoff, names_to=c("measure")) %>%
	ggplot(aes(cutoff, value, colour=measure)) + geom_line(size=1) +
	scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
	scale_y_continuous(breaks = scales::pretty_breaks(n = 10))

	# cross validation

	featuresets<-data.frame(read.csv("../data/featuresets.tsv", header = FALSE, sep = "\t", dec=".", quote="", encoding="UTF-8",stringsAsFactors=FALSE))

	collected_results<-tibble()
	for (i in c(1:nrow(featuresets))) {
	set.seed(325)
	train_control <- trainControl(method="repeatedcv", number=5, repeats=5, verboseIter=T,classProbs = TRUE,
	savePredictions = "final")
	vars<-unlist(strsplit(featuresets[i,2],","))
	fml <- as.formula(paste("CO_IDIOM ~ ", paste(vars, collapse= "+")))
	rf <- train(fml,
	data=ngramme1,
	method='rf',
	tuneGrid = data.frame(.mtry = floor(sqrt(length(vars)))),
	cutoff=c(0.3,0.7),
	trControl=train_control)
	res<-confusionMatrix(rf$pred$pred, rf$pred$obs,mode="everything")
	collected_results <-bind_rows(collected_results, c(features=featuresets[i,1], res$byClass))
	}

	collected_results1<-collected_results[c(2:16),c(1,6,7,8,12)]
	collected_results1<-collected_results1%>%
	mutate(across(everything(), ~ replace_na(.x, 0))) %>%
	mutate(across(c(2:5), ~ round(as.numeric(.x),digits=3)))

	# Analysing tradeoff between Fscore, Recall, Precision for various cutoffs
	# full range from precision almost 100% to recall almost 100%


	vars<-unlist(strsplit(featuresets[2,2],","))
	fml <- as.formula(paste("CO_IDIOM ~ ", paste(vars, collapse= "+")))
	cvalues<-tibble()
	for (c in c(seq(from=0.4,to=0.99,by=0.025),0.999)) {
	set.seed(325)
	train_control <- trainControl(method="repeatedcv", number=5, repeats=5, verboseIter=T,classProbs = TRUE,
	savePredictions = "final")
	rf <- train(fml,
	data=ngramme1,
	method='rf',
	tuneGrid = data.frame(.mtry = floor(sqrt(length(vars)))),
	cutoff=c(1-c,c),
	trControl=train_control)
	res<-confusionMatrix(rf$pred$pred, rf$pred$obs,mode="everything")
	cvalues <-bind_rows(cvalues, c(cutoff=c, res$byClass))
	}
	cvalues %>%
	select(c("cutoff", "Recall", "Precision", "F1", "Specificity", "Balanced Accuracy")) %>%
	pivot_longer(!cutoff, names_to=c("measure")) %>%
	ggplot(aes(cutoff, value, colour=measure)) + geom_line(size=1) +
	scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
	scale_y_continuous(breaks = scales::pretty_breaks(n = 10))

	cvaluesall<-cvalues


	# Using estimates by random forest on entire dataset

	library(randomForest)
	rf_classifier_full = randomForest(all_fml_1, data=ngramme1, importance=TRUE, cutoff=c(0.3,0.7))
	rf_classifier_full
	# class.error is 1 - recall
	varImpPlot(rf_classifier_full)

	# Feature ranking

	# rf features as table

	# correlated features seem to split their rankings

	rfranks<-importance(rf_classifier_full)[,3:4]

	# ttest

	idioms<-ngramme1 %>% filter(CO_IDIOM == "idiom")
	nonidioms<-ngramme1 %>% filter(CO_IDIOM != "idiom")

	m1<-mean(unlist(nonidioms["SY_C1_LD"]),na.rm=T)

	nonidioms1<-nonidioms%>%filter(SY_C1_LD > m1)

	mean(unlist(idioms["SY_C1_LD"]),na.rm=T)

	idioms1<-ngramme1 %>% filter(CO_IDIOM == "idiom")
	idioms2<-ngramme2 %>% filter(CO_IDIOM == "idiom")

	ttestPvalues<-sapply(all_vars_1,
	function(sel) t.test(idioms[sel],nonidioms[sel])$p.value)

	t.test(idioms["SY_C1_LL"],nonidioms["SY_C1_LL"])

	ttestSignificance<-sapply(all_vars_1,
	function(sel) {
	p<-t.test(idioms[sel],nonidioms[sel])$p.value
	if (p < 0.001) {
	return("***")
	}
	if (p < 0.01) {
	return("**")
	}
	if (p < 0.05) {
	return ("*")
	}
	return(" ")
	})


	ttestSignificance1<-sapply(all_vars_1,
	function(sel) {
	p<-t.test(idioms[sel],nonidioms1[sel])$p.value
	if (p < 0.001) {
	return("***")
	}
	if (p < 0.01) {
	return("**")
	}
	if (p < 0.05) {
	return ("*")
	}
	return(" ")
	})


	# information gain
	# multiply by 1000 to avoid undersized bins
	# features are ranked individually no matter their correlation
	igain<-information.gain(all_fml_1, data=ngramme1%>%mutate_at(all_vars_1, ~ . * 1000),unit="log2")

	# difference between means (positive or negative?)

	diffMeans<-sapply(all_vars_1,function(sel) mean(unlist(idioms[sel]),na.rm=T)-mean(unlist(nonidioms[sel]),na.rm=T))

	diffMeansSign<-sapply(all_vars_1,function(sel) ifelse(mean(unlist(idioms[sel]),na.rm=T)-mean(unlist(nonidioms[sel]),na.rm=T)>0,"+","-"))

	diffMeansSign1<-sapply(all_vars_1,function(sel) ifelse(mean(unlist(idioms[sel]),na.rm=T)-mean(unlist(nonidioms1[sel]),na.rm=T)>0,"+","-"))


	featurenames[,c("newnames","explanation")]

	featureRanks<-cbind(rownames(rfranks),rfranks,igain,ttestPvalues,ttestSignificance,diffMeans,diffMeansSign,ttestSignificance1,diffMeansSign1)
	colnames(featureRanks)[1]<-"newnames"

	featureRanks<-merge(featureRanks,featurenames[,c("newnames","explanation")],by="newnames")



	#randomForestExplainer::explain_forest(rf_classifier )

	# averate estimates and feature ranks over 10 runs

	errrate<-0
	conf<-matrix(0,2,3)
	featureRanks<-matrix(0,4,length(vars))
	for (i in 1:10) {
	rfc =randomForest(fmla, data=ngramme, importance=TRUE)
	#rfc =randomForest(fmla, data=ngramme, importance=TRUE, cutoff=c(0.2, 0.8))
	errrate<-errrate+rfc$err.rate[100,1]
	conf<-conf+rfc$confusion
	featureRanks<-featureRanks+
	cbind(importance(rfc)[,3:4],
	information.gain(fmla, data=ngramme%>%mutate_at(vars, ~ . * 1000),unit="log2"),
	sapply(vars,
	function(sel) t.test(idioms[sel],nonidioms[sel])$p.value))
	print(errrate/i)
	conf1<-round(
	rbind(
	cbind(conf[,1:2]/i,(1-conf[,3]/i)*100),
	c(100*diag(conf[,1:2])/colSums(conf[,1:2]),NA),
	c(rowSums(conf[,1:2]/i),NA)),digits=2)
	colnames(conf1)<-c("1","0","rec")
	rownames(conf1)<-c("1","0","prec","sum")
	print(conf1)
	}
	featureRanks<-featureRanks/10
	colnames(featureRanks)<-c("MeanDecreaseAccuracy","MeanDecreaseGini","InformationGain","Ttest")