simulation.Rnw

\documentclass{article}
\usepackage[left=20mm, right=15mm, top=15mm, bottom=20mm]{geometry}

\begin{document}
%\VignetteEngine{knitr::knitr}

<<setup>>=
rm(list = ls())

library("knitr")
opts_chunk$set(cache = FALSE, echo = TRUE, message = FALSE)

library("partykit")
library("sandwich")
library("plyr")
library("mvtnorm")
library("ggplot2")
theme_set(theme_bw(base_size = 18))


source("basis/personalised_models.R")
source("basis/variable_importance.R")
source("basis/dependence_plots.R")
source("basis/pruning.R")

set.seed(123)
@

\section{Simulation of data}
<<datasimulation>>=
sim_data <- function(n = 500, p = 10, beta = 3, sd = 1){
  
  lev <- c("C", "A")
  a <- rep(factor(lev, labels = lev, levels = lev), length = n)
  
  ## z variables are correlated
  sigma <- diag(p) 
  sigma[sigma == 0] <- 0.2
  ztemp <- rmvnorm(n, sigma = sigma)
  z <- (pnorm(ztemp) * 2 * pi) - pi  
  
  #   z <- matrix(runif(p * n, min = -pi, max = pi), nrow = n)
  colnames(z) <- paste0("z", 1:ncol(z))
  z1 <- z[,1]
  
  y <- 1.9 + 0.2 * (a %in% "A") + beta * cos(z1) * (a %in% "A") + rnorm(n, 0, sd)
  
  data.frame(y = y, a = a, z)
}

@


\section{The function needed for fitting the cforest}
<<lm_fittingfunction>>=
##' fitting function for cforest with linear model
lmf <- function(data, weights, parm = c(1,2)) {
  
  tb <- table(data[["a"]][weights > 0])
  ## only one treatment arm left (or only one observation in one treatment arm); 
  ## we don't want to split further...
  if (any(tb < 2)) return(matrix(0, nrow = nrow(data), ncol = length(parm)))
  
  mod <- lm(y ~ a, weights = weights, data = data, subset = weights > 0)
  ef <- as.matrix(estfun(mod)[, parm])
  ret <- matrix(0, nrow = nrow(data), ncol = ncol(ef))
  ret[weights > 0,] <- ef
  ret
}
@

 
 \section{Computing individualised models}
 
 \Sexpr{knit_child("personalised_models.Rnw")}
 
<<sind_mods>>=
 beta <- 3
 train <- sim_data(p = 10, beta = beta, n = 600)
 test <- sim_data(p = 10, beta = beta, n = 600)
 Z <- names(train)[!(names(train) %in% c("y", "a"))]
 cfm <- as.formula(paste("y + a ~", paste(Z, collapse = "+")))
 
 lmforest <-  cforest(cfm, data = train, ytrafo = lmf, ntree = 100, cores = 1, 
                      perturb = list(replace = FALSE))
 pmods <- person_mods(object = lmforest, basemod = "lm", newdata = test, 
                      OOB = FALSE, parallel = TRUE)
@
 
 \section{Variable importance}
 
 \Sexpr{knit_child("variable_importance.Rnw")}
 
<<svarimp, out.width = "0.7\\textwidth", out.height = "0.7\\textwidth">>=
 ###### variable importance ######
 lmforest_pruned <- prune_forest(lmforest, endpoint = "other", 
                                 treatmentvar = "a", minobs = 2)
 VI <- varimp(forest = lmforest_pruned, loglik = comp_loglik.lm, 
              basemod = "lm", OOB = TRUE, parallel = TRUE)
 VI
 
 text.z <- sapply(Z, function(x) bquote(z[.(gsub("z", "", x))]))
 VI$variable <- factor(VI$variable, levels = Z)
 
 save(VI, file = "sim_varimp.rda")
 
 ggplot(VI, aes(y = VI, x = variable)) + geom_bar(stat = "identity", width = .1) + 
   coord_flip() + scale_x_discrete(labels = text.z) +
   theme(panel.grid.major.y = element_blank())
 
@
 
 
 \section{Comparison to regular forest}
<<cforest>>=
 
 set.seed(22)
 fm_cf <- as.formula("y ~ . + .:a")
 condforest <- cforest(fm_cf, data = train, ntree = 100, cores = 50, 
                       perturb = list(replace = FALSE))
 
 t0 <- t1 <- test
 t1$a[t1$a == "C"] <- "A"
 t0$a[t1$a == "A"] <- "C"
 
 trt_eff_condforest <- predict(condforest, newdata = t1, type = "response") -
   predict(condforest, newdata = t0, type = "response")
 
 
 trt_eff_pmods <- sapply(pmods, coef)["aA", ]
 
 
 trt_eff <- data.frame(trt_eff = c(trt_eff_pmods, trt_eff_condforest), 
                       type = rep(c("pmods", "cforest"), each = nrow(train)),
                       z1 = rep(test$z1, times = 2))
 
 save(trt_eff, file = "sim_compareforests.rda")
 
 
 p <- ggplot(data = trt_eff) +
   geom_line(aes(x = z1, y = (0.2 + beta * cos(z1)))) + 
   geom_point(aes(x = z1, y = trt_eff, color = type), alpha = 0.6) +
   ylab(expression(beta(z[1]))) +  xlab(bquote(z[.(gsub("z", "", Z[1]))])) 
 p
@
 
 
<<pruned_forest>>=
 pmods_pruned <- person_mods(object = lmforest_pruned, basemod = "lm", newdata = test, 
                             OOB = FALSE, parallel = TRUE)
 
 trt_eff_pruned <- data.frame(trt_eff = sapply(pmods_pruned, coef)["aA", ], 
                              type = rep(c("pmods_pruned"), each = nrow(train)),
                              z1 = test$z1)
 p + geom_point(data = trt_eff_pruned, aes(x = z1, y = trt_eff, color = type), alpha = 0.6)
@
 

\section{Function which computes log-likelihoods for simulated data}

<<sim_complogliks>>=
##' simulation for cforest check
##' comparison with two linear models
##' (1) lm(y ~ a, data = train)
##' (2) lm(y ~ a * I(z1 > 0), data = train)
##' @param i unused
##' @param data list of 2 list(train, test). If NULL training data of size n and test data of size m are obtained.
sim_wi <- function(i, data = NULL, n = 600, m = 600, p = 10, beta = 3,
                   person.mods = FALSE, comp.cforest = FALSE){
  
  if(is.null(data)) {
    data <- sim_data(n = (n + m), p = p, beta = beta)
    train <- data[1:n, ]
    test <- data[(n+1):(n+m), ]
  } else {
    train <- data$train
    test <- data$test
  }
  
  Z <- names(train)[!(names(train) %in% c("y", "a"))]
  
  
  ### linear model 1 ############################
  lm1 <- lm(y ~ a, data = train)
  loglik1 <- sum(comp_loglik.lm(lm1, test, response = "y"))
  
  
  ### linear model 1 ############################
  lm2 <- lm(y ~ a * cos(z1), data = train)
  loglik2 <- sum(comp_loglik.lm(lm2, test, response = "y"))
  
  ### linear model random forest ################
  cfm <- as.formula(paste("y + a ~", paste(Z, collapse = "+")))
  lmforest <-  cforest(cfm, data = train, ytrafo = lmf, ntree = 100, cores = 50, 
                       perturb = list(replace = FALSE))
  
  mods <- person_mods(object = lmforest, basemod = "lm", newdata = test)
  logLiksf <- sapply(1:nrow(test), comp_loglik, mods = mods, dat = test, 
                     basemod = "lm", loglik = NULL)
  loglikf <- sum(logLiksf)
  
  
  if (person.mods) {
    pm <- cbind(t(sapply(mods, coef)), test)
    
    if(comp.cforest) {
      condforest1 <- cforest(y ~ ., data = train, ntree = 100, cores = 50, 
                            perturb = list(replace = FALSE))
      fm_cf <- as.formula("y ~ . + .:a")
      condforest2 <- cforest(fm_cf, data = train, ntree = 100, cores = 50, 
                            perturb = list(replace = FALSE))
      t0 <- t1 <- test
      t1$a[t1$a == "C"] <- "A"
      t0$a[t1$a == "A"] <- "C"
      
      trt_eff_condforest1 <- predict(condforest1, newdata = t1, type = "response") -
        predict(condforest1, newdata = t0, type = "response")
      trt_eff_condforest2 <- predict(condforest2, newdata = t1, type = "response") -
        predict(condforest2, newdata = t0, type = "response")
      
      trt.eff_cforest1 <- cbind(trt_eff = trt_eff_condforest1, test)
      trt.eff_cforest2 <- cbind(trt_eff = trt_eff_condforest2, test)

      return(list(logliks = c(loglik1, loglik2, loglikf),
                  person.mods = pm,
                  trt.eff_cforest1 = trt.eff_cforest1,
                  trt.eff_cforest2 = trt.eff_cforest2))
    } else{
      return(list(logliks = c(loglik1, loglik2, loglikf),
                  person.mods = pm))
    }
  } else {
    return(c(loglik1, loglik2, loglikf))
  }
  
  
  
}
@

\section{Log-likelihoods for simulated data for the three methods}

%\Sexpr{knit_child("log_likelihoods.Rnw")}

<<slogliks, out.width = "0.6\\textwidth", out.height = "0.6\\textwidth">>=
set.seed(223)
logliks.pm <- lapply(1:100, sim_wi, person.mods = TRUE, comp.cforest = TRUE)
logliks <- t(sapply(logliks.pm, function(x) x$logliks))
llnames <- c("incorrect model", "correct model", "model based forest")
colnames(logliks) <- llnames
logliks

logliks_gg <- data.frame(id = rep(1:nrow(logliks), times = 3),
                         nam = rep(llnames, each = nrow(logliks)), 
                         loglik = unlist(as.list(logliks)))

save(logliks_gg, file = "sim_logliks.rda")

ggplot(logliks_gg, aes(x = nam, y = loglik)) + 
  geom_boxplot() +
  geom_line(aes(group = id), alpha = 0.2) +
  xlab("") + ylab("log-likelihood") +
  theme(panel.grid.major.x = element_blank())
@


\section{Comparison forests in simulation}
<<scompare_error>>=
msefun <- function(error) {
  mean(error^2)
}

get_error <- function(x, df = FALSE, beta = 3) {
  
  trt_eff <- x$person.mods[, c("aA", "z1")]
  names(trt_eff)[names(trt_eff) == "aA"] <- "trt_eff"
  trt_eff$type <- "pmods"
  trt_eff <- rbind(trt_eff,
                   cbind(x$trt.eff_cforest1[, c("trt_eff", "z1")], type = "cforest1"),
                   cbind(x$trt.eff_cforest2[, c("trt_eff", "z1")], type = "cforest2"))
  
  trt_eff$true <- (0.2 + beta * cos(trt_eff$z1))
  trt_eff$err <- trt_eff$trt_eff - trt_eff$true
  
  if(df) {
    return(trt_eff)
  } else {
    c(mse_cforest1 = msefun(trt_eff$err[trt_eff$type == "cforest1"]),
      mse_cforest2 = msefun(trt_eff$err[trt_eff$type == "cforest2"]),
      mse_pmods = msefun(trt_eff$err[trt_eff$type == "pmods"]))
  }
}


mse <- sapply(logliks.pm, get_error)
save(mse, file = "sim_comparemse.rda")
@


\section{Partial dependence plots}
With out of sample data.
\Sexpr{knit_child("partial_dependence_plots.Rnw")}

<<spdplots, out.width = "0.5\\textwidth", out.height = "0.4\\textwidth", fig.show='hold', message=FALSE>>=
###### partial dependence plots ######
ggplot(data = train, aes(x = z1, y = (0.2 + beta * cos(z1)) * (a %in% "A"), group = a)) +
  geom_line() + ylab(expression(beta(z[1]))) +  xlab(bquote(z[.(gsub("z", "", Z[1]))])) 



# personalized models
library(plyr)
pd <- ldply(1:length(logliks.pm), 
            function(i) {
              logliks.pm[[i]]$person.mods$iteration <- i
              return(logliks.pm[[i]]$person.mods)
            })

pd$iteration <- as.factor(pd$iteration)

save(train, file = "sim_traindat.rda")
save(Z, pd, file = "sim_dependenceplotdat.rda")

dependenceplot.sim <- function(variable, treatment) {
  ggplot(data = pd, aes_string(x = variable, y = treatment)) + 
    geom_point(alpha = 0.1) + 
    xlab(bquote(z[.(gsub("z", "", variable))])) + 
    ylab(bquote(bar(beta)(z[.(gsub("z", "", variable))])))+ 
    xlim(-pi, pi) + ylim(-pi, pi)              
}
p <- lapply(Z, dependenceplot.sim, treatment = "aA")

p[[1]] + geom_line(aes(x = z1, y = (0.2 + beta * cos(z1)))) 
p[[2]]

@



\end{document}