diff --git a/R/idiomclassification_mk_pf3.R b/R/idiomclassification_mk_pf3.R
new file mode 100644
index 0000000..7042ad6
--- /dev/null
+++ b/R/idiomclassification_mk_pf3.R
@@ -0,0 +1,361 @@
+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")
+
+
+
+