etlp_predict_using_glm.rmd

---
title: "Interpretable Analysis of School Policy Decisions, linear models"
author: "Charles Saluski"
# date: "1/4/2022"
output: pdf_document
---

```{r}
library(glmnet)
library(mlr3)
library(mlr3learners)
library(data.table)
library(mlr3extralearners)
library(mlr3tuning)

csv.data.loc <- "./Data Sources CSV"
ic.joined.dt.loc <- paste(csv.data.loc, "/ic.cwis.nces.computed.combined.csv", sep = "")
cwis.joined.dt.loc <- paste(csv.data.loc, "/cwis.nces.computed.combined.csv", sep = "")
cl.joined.dt.loc <- paste(csv.data.loc, "/ic.cwis.nces.cl.computed.combined.csv", sep = "")

ic.joined.dt <- as.data.table(read.csv(ic.joined.dt.loc))
ic.joined.dt <- ic.joined.dt[complete.cases(ic.joined.dt)]
cwis.joined.dt <- as.data.table(read.csv(cwis.joined.dt.loc))
cl.joined.dt <- as.data.table(read.csv(cl.joined.dt.loc))
```

```{r}
exclude.cols <- c("X", "State.District.ID", "CWIS_session", "IC_NCES.District.Name..to.check.", "IC_School.District", "IC_Teacher_leader_More_than_6", "IC_Total_more_than_10")

ic.predict.dt <- ic.joined.dt[, !..exclude.cols]
ic.predict.no.cfa.dt <- ic.joined.dt[, !c("CWIS_CFA_avg", ..exclude.cols)]
ic.predict.cfa.dt <- ic.joined.dt[, !c("CWIS_ETLP_avg", ..exclude.cols)]

cwis.predict.dt <- cwis.joined.dt[, !..exclude.cols]
cwis.predict.no.cfa.dt <- cwis.joined.dt[, !c("CWIS_CFA_avg", ..exclude.cols)]
cwis.predict.cfa.dt <- cwis.joined.dt[, !c("CWIS_ETLP_avg", ..exclude.cols)]

cl.predict.dt <- cl.joined.dt[, !..exclude.cols]
cl.predict.no.cfa.dt <- cl.joined.dt[, !c("CWIS_CFA_avg", ..exclude.cols)]
cl.predict.cfa.dt <- cl.joined.dt[, !c("CWIS_ETLP_avg", ..exclude.cols)]
```

```{r}
set.seed(123)
num.folds <- 10

task.full.regr <- TaskRegr$new(id = "ic.etlp", backend = ic.predict.dt, target = "CWIS_ETLP_avg")

task.no.cfa.regr <- TaskRegr$new(id = "ic.etlp.no.cfa", backend = ic.predict.no.cfa.dt, target = "CWIS_ETLP_avg")

task.cfa.regr <- TaskRegr$new(id = "ic.cfa", backend = ic.predict.cfa.dt, target = "CWIS_CFA_avg")


task.cwis.full.regr <- TaskRegr$new(id = "cwis.etlp", backend = cwis.predict.dt, target = "CWIS_ETLP_avg")

task.cwis.no.cfa.regr <- TaskRegr$new(id = "cwis.etlp.no.cfa", backend = cwis.predict.no.cfa.dt, target = "CWIS_ETLP_avg")

task.cwis.cfa.regr <- TaskRegr$new(id = "cwis.cfa", backend = cwis.predict.cfa.dt, target = "CWIS_CFA_avg")


task.cl.regr <- TaskRegr$new(id = "cl", backend = cl.predict.dt, target = "CWIS_ETLP_avg")

task.cl.no.cfa.regr <- TaskRegr$new(id = "cl.no.cfa", backend = cl.predict.no.cfa.dt, target = "CWIS_ETLP_avg")

task.cl.cfa.regr <- TaskRegr$new(id = "cl.cfa", backend = cl.predict.cfa.dt, target = "CWIS_CFA_avg")

task.name.vec <- c("ic.etlp", "ic.etlp.no.cfa", "ic.cfa", "cwis.etlp", "cwis.etlp.no.cfa", "cwis.cfa", "cl", "cl.no.cfa", "cl.cfa")

task.list <- list(task.full.regr, task.no.cfa.regr, task.cfa.regr, task.cwis.full.regr, task.cwis.no.cfa.regr, 
task.cwis.cfa.regr, 
task.cl.regr, task.cl.no.cfa.regr, task.cl.cfa.regr)

learner.list <- list()
learner.list[["regr.featureless"]] <- LearnerRegrFeatureless$new()
learner.list[["regr.ctree"]] <- LearnerRegrCTree$new()
# cv_glmnet returns 2 models, one with s1 and one with minimum
learner.list[["regr.cv_glmnet"]] <- LearnerRegrCVGlmnet$new()
learner.list[["regr.cforest"]] <- LearnerRegrCForest$new()

learner.list[["regr.xgboost.tuned"]] <- AutoTuner$new(
  learner = lrn("regr.xgboost"),
  resampling = rsmp("cv", folds = 3),
  measure = msr("regr.mse"),
  search_space = ps(
    eta = p_dbl(lower = 0, upper = 1),
    nrounds = p_int(lower = 1, upper = 16)
  ),
  terminator = trm("none"),
  tuner = tnr("grid_search", resolution = 5),
  store_tuning_instance = TRUE
)

learner.name.vec <- names(learner.list)

resampling <- rsmp("cv", folds = num.folds)

benchmark.obj <- benchmark_grid(
  # tasks, learners, and resamplings
  # we'll only give a learner vector, same tasks and resamplings
  task = task.list,
  learners = learner.list,
  resamplings = list(resampling)
)

benchmark.res <- benchmark(benchmark.obj, store_models = TRUE)

measure <- msr("regr.mse")
```

```{r}
result.dt <- benchmark.res$score(measure)

score.result.list <- list()
# we can't do a for loop over the learner.name.vec because the cv_glmnet needs
# to be run twice, once for s1 and once for minimum
for (method in learner.name.vec) {
  for (task.name in task.name.vec) {
    curr.dt <- result.dt[learner_id == method & task_id == task.name]
    method.learner.list <- curr.dt$learner
    for (i in 1:num.folds) {
      if (method == "regr.cv_glmnet") {
        curr.model <- method.learner.list[[i]]$model
        lambda.min.index <- curr.model$index[1]
        mse.min <- curr.model$cvm[lambda.min.index]
        lambda.1se.index <- curr.model$index[2]
        mse.1se <- curr.model$cvm[lambda.1se.index]
        score.result.list[[paste(method, task.name, "1se", i, sep = ".")]] <- data.table(
          method = paste(method, "1se", sep = "."),
          fold = i,
          mse.loss = mse.1se,
          task.name
        )
        score.result.list[[paste(method, task.name, "min", i, sep = ".")]] <- data.table(
          method = paste(method, "min", sep = "."),
          fold = i,
          mse.loss = mse.min,
          task.name
        )
      } else {
        score.result.list[[paste(method, task.name, i, sep = ".")]] <- data.table(
          method = paste(method, sep = "."),
          fold = i,
          mse.loss = curr.dt[i, "regr.mse"][[1]],
          task.name
        )
      }
    }
  }
}
err.dt <- do.call(rbind, score.result.list)
```


LASSO models prove to be much more accurate than the featureless model, disproving the null hypothesis. 
```{r}
library(ggplot2)
method.levels <- err.dt[, .(mean = mean(mse.loss)), by = method][order(-mean), method]
err.dt[, Method := factor(method, method.levels)]
err.plot <- ggplot() +
  geom_point(data = err.dt, aes(x = mse.loss, y = Method)) +
  ggtitle(paste("MSE loss by method and task")) +
  facet_grid(task.name ~ .)

png(filename = "./img_out/glmnet/regr.loss.mse.all.png", width = 6, height = 12, unit = "in", res = 200)
print(err.plot)
dev.off()


err.plot <- ggplot() +
  geom_point(data = err.dt[Method %in% c("regr.cv_glmnet.min", "regr.cv_glmnet.1se", "regr.featureless")], aes(x = mse.loss, y = Method)) +
  ggtitle(paste("MSE loss by method and task")) +
  facet_grid(task.name ~ .)

png(filename = "./img_out/glmnet/regr.loss.mse.both.png", width = 6, height = 8, unit = "in", res = 200)
print(err.plot)
dev.off()


err.plot <- ggplot() +
  geom_point(data = err.dt[Method %in% c("regr.cv_glmnet.min", "regr.cv_glmnet.1se", "regr.featureless") & task.name == "ic.etlp"], aes(x = mse.loss, y = Method)) +
  ggtitle(paste("MSE loss by method and task")) +
  facet_grid(task.name ~ ., labeller = label_both)

png(filename = "./img_out/glmnet/regr.loss.mse.etlp.png", width = 6, height = 4, unit = "in", res = 200)
print(err.plot)
dev.off()
```

Now we examine the factors that are found to be imporant in the models.
```{r}
# we want a dt with each model's coefficients
# then count and display which coefficients are important
cv.glm.dt <- result.dt[learner_id == "regr.cv_glmnet"]
glm.method.v <- c("lambda.min", "lambda.1se")
glm.coef.list <- list()
for (task.name in task.name.vec) {
  for (method in glm.method.v) {
    curr.dt <- cv.glm.dt[task_id == task.name]
    for (fold in 1:num.folds) {
      curr.coef.mat <- as.matrix(
        coef(curr.dt[iteration == fold]$learner[[1]]$model, s = method)[-1, ]
      )
      glm.coef.list[[paste(method, task.name, fold)]] <- data.table(
        method,
        var = rownames(curr.coef.mat),
        coef = as.numeric(curr.coef.mat),
        task_id = task.name
      )
    }
  }
}

# this dt has columns of coefs of each var and a column with the method
glm.coef.dt <- do.call(rbind, glm.coef.list)

# dt with var method coef
# make count
glm.coef.dt[, count := sum(coef != 0), by = .(method, task_id, var)]

```

```{r}
coef.file.dest <- "./Data Sources CSV/regr.glm.coef.csv"

write.csv(glm.coef.dt, file = coef.file.dest, row.names=FALSE)
```


```{r}
for (method.select in glm.method.v) {
  for (task.name in task.name.vec) {
    var.coef.plot <- ggplot() +
      geom_point(data = glm.coef.dt[method.select == method & task_id == task.name & count > 0], aes(x = coef, y = var)) +
      facet_grid(count ~ ., scales = "free", space = "free") +
      ggtitle(paste("Coefficients of model ", method.select, " in task ", task.name)) 
    # scale_y_continuous(breaks=1:num.folds)
    filename <- paste("img_out/glmnet/glm_results/",method.select, task.name, ".png", sep = "")
    print(filename)
    png(filename = filename, width = 12, height = 8, unit = "in", res = 200)
    print(var.coef.plot)
    dev.off()
  }
}
```