inst/india_summary.R

In mboost: Model-Based Boosting

library("mboost")
load("india.Rdata")

inds <- paste("w", 1:44, sep="")
quantiles <- c(0.05, 0.1, 0.5)
pred1 <- pred2 <- pred3 <- matrix(0, ncol=4, nrow=length(inds))
colnames(pred1) <- colnames(pred2) <- colnames(pred3) <- c("add","vcm","stump","tree")
add1.selected <- add2.selected <- add3.selected <- list()
vcm1.selected <- vcm2.selected <- vcm3.selected <- list()
add1.stabsel <- add2.stabsel <- add3.stabsel <- matrix(NA, ncol=length(inds), nrow=210000)
add1.its <- add2.its <- add3.its <- rep(0, length(inds))
vcm1.its <- vcm2.its <- vcm3.its <- rep(0, length(inds))
vcm1.stabsel <- vcm2.stabsel <- vcm3.stabsel <- matrix(NA, ncol=length(inds), nrow=210000)

loss <- function(y, f, tau)
  {
  sum(tau*(y - f)*((y - f) >= 0) - (1-tau)*(y - f)*((y - f) < 0))
  }

# read the predictions
  
for(i in 1:length(inds))
  {
  print(i)
  indiapred <- india[india[,inds[i]]==3,]
  indiapred <- cbind(indiapred, ff = 0)
  names(indiapred)[71] <- "FALSE"

  load(paste("additive/model",i,"_0.05.Rdata",sep=""))
  pred <- predict(model, newdata=indiapred)
  pred1[i,1] <- loss(indiapred$stunting, pred, tau=0.05)
  add1.selected[[i]] <- factor(selected(model), levels=1:28, labels=names(model$baselearner))
  add1.stabsel[1:mstop(model),i] <- selected(model) 
  add1.its[i] <- mstop(model)

  load(paste("additive/model",i,"_0.1.Rdata",sep=""))
  pred <- predict(model, newdata=indiapred)
  pred2[i,1] <- loss(indiapred$stunting, pred, tau=0.1)
  add2.selected[[i]] <- factor(selected(model), levels=1:28, labels=names(model$baselearner))
  add2.stabsel[1:mstop(model),i] <- selected(model) 
  add2.its[i] <- length(selected(model))

  load(paste("additive/model",i,"_0.5.Rdata",sep=""))
  pred <- predict(model, newdata=indiapred)
  pred3[i,1] <- loss(indiapred$stunting, pred, tau=0.5)
  add3.selected[[i]] <- factor(selected(model), levels=1:28, labels=names(model$baselearner))
  add3.stabsel[1:mstop(model),i] <- selected(model) 
  add3.its[i] <- length(selected(model))


  load(paste("VCM/model",i,"_0.05.Rdata",sep=""))
  pred <- predict(model, newdata=indiapred)
  pred1[i,2] <- loss(indiapred$stunting, pred, tau=0.05)
  vcm1.selected[[i]] <- factor(selected(model), levels=1:54, labels=names(model$baselearner))
  vcm1.stabsel[1:mstop(model),i] <- selected(model) 
  vcm1.its[i] <- length(selected(model))

  load(paste("VCM/model",i,"_0.1.Rdata",sep=""))
  pred <- predict(model, newdata=indiapred)
  pred2[i,2] <- loss(indiapred$stunting, pred, tau=0.1)
  vcm2.selected[[i]] <- factor(selected(model), levels=1:54, labels=names(model$baselearner))
  vcm2.stabsel[1:mstop(model),i] <- selected(model) 
  vcm2.its[i] <- length(selected(model))

  load(paste("VCM/model",i,"_0.5.Rdata",sep=""))
  pred <- predict(model, newdata=indiapred)
  pred3[i,2] <- loss(indiapred$stunting, pred, tau=0.5)
  vcm3.selected[[i]] <- factor(selected(model), levels=1:54, labels=names(model$baselearner))
  vcm3.stabsel[1:mstop(model),i] <- selected(model) 
  vcm3.its[i] <- length(selected(model))


  load(paste("stumps/model",i,"_0.05.Rdata",sep=""))
  pred1[i,3] <- loss(indiapred$stunting, pred, tau=0.05)

  load(paste("stumps/model",i,"_0.1.Rdata",sep=""))
  pred2[i,3] <- loss(indiapred$stunting, pred, tau=0.1)

  load(paste("stumps/model",i,"_0.5.Rdata",sep=""))
  pred3[i,3] <- loss(indiapred$stunting, pred, tau=0.5)


  load(paste("trees/model",i,"_0.05.Rdata",sep=""))
  pred1[i,4] <- loss(indiapred$stunting, pred, tau=0.05)

  load(paste("trees/model",i,"_0.1.Rdata",sep=""))
  pred2[i,4] <- loss(indiapred$stunting, pred, tau=0.1)

  load(paste("trees/model",i,"_0.5.Rdata",sep=""))
  pred3[i,4] <- loss(indiapred$stunting, pred, tau=0.5)
  }
  
  

# plotting of effects in the additive model

preddata <- preddatafm <- list()
plotdata.add1 <- plotdata.add2 <- plotdata.add3 <- list()
plotdata.vcm1 <- plotdata.vcm2 <- plotdata.vcm3 <- list()

varnames <- c("cage", "breastfeeding", "mbmi", "mage", "medu", "edupartner",
            "csex", "ctwin", "cbirthorder", "munemployed", "mreligion",
            "mresidence", "deadchildren", "wealth", "electricity", "radio",
            "television", "refrigerator", "bicycle", "motorcycle", "car")

refval <- list()
for(i in 1:length(varnames))
  {
  print(i)
  if(is.numeric(india[,varnames[i]]))
    refval[[i]] <- mean(india[,varnames[i]])
  if(is.factor(india[,varnames[i]]))
    refval[[i]] <- unique(india[,varnames[i]])[which.max(table(india[,varnames[i]]))]
  }

for(i in 1:length(varnames))
  {
  print(i)
  vals <- unique(india[,varnames[i]])
  preddata[[i]] <- india[1:length(vals),]
  for(j in 1:length(varnames))
    {
    if(j==i)
      preddata[[i]][,varnames[j]] <- vals
    else
      preddata[[i]][,varnames[j]] <- rep(refval[[j]], length(vals))
    }
  preddata[[i]] <- cbind(preddata[[i]], ff = 0)
  names(preddata[[i]])[71] <- "FALSE"
    
  preddatafm[[i]] <- preddata[[i]]
  preddatafm[[i]][,"csex"] <- rep(unique(india[,"csex"])[which.min(table(india[,"csex"]))], length(vals))


  plotdata.add1[[i]] <- data.frame(vals, matrix(0, nrow=length(vals), ncol=length(inds)))
  names(plotdata.add1[[i]]) <- c(varnames[i],paste("f", 1:length(inds), sep=""))

  plotdata.add2[[i]] <- data.frame(vals, matrix(0, nrow=length(vals), ncol=length(inds)))
  names(plotdata.add1[[i]]) <- c(varnames[i],paste("f", 1:length(inds), sep=""))

  plotdata.add3[[i]] <- data.frame(vals, matrix(0, nrow=length(vals), ncol=length(inds)))
  names(plotdata.add1[[i]]) <- c(varnames[i],paste("f", 1:length(inds), sep=""))


  plotdata.vcm1[[i]] <- data.frame(vals, matrix(0, nrow=length(vals), ncol=length(inds)))
  names(plotdata.vcm1[[i]]) <- c(varnames[i],paste("fm", 1:length(inds), sep=""))
  plotdata.vcm1[[length(varnames)+i]] <- data.frame(vals, matrix(0, nrow=length(vals), ncol=length(inds)))
  names(plotdata.vcm1[[length(varnames)+i]]) <- c(varnames[i],paste("ff", 1:length(inds), sep=""))

  plotdata.vcm2[[i]] <- data.frame(vals, matrix(0, nrow=length(vals), ncol=length(inds)))
  names(plotdata.vcm2[[i]]) <- c(varnames[i],paste("fm", 1:length(inds), sep=""))
  plotdata.vcm2[[length(varnames)+i]] <- data.frame(vals, matrix(0, nrow=length(vals), ncol=length(inds)))
  names(plotdata.vcm2[[length(varnames)+i]]) <- c(varnames[i],paste("ff", 1:length(inds), sep=""))

  plotdata.vcm3[[i]] <- data.frame(vals, matrix(0, nrow=length(vals), ncol=length(inds)))
  names(plotdata.vcm3[[i]]) <- c(varnames[i],paste("fm", 1:length(inds), sep=""))
  plotdata.vcm3[[length(varnames)+i]] <- data.frame(vals, matrix(0, nrow=length(vals), ncol=length(inds)))
  names(plotdata.vcm3[[length(varnames)+i]]) <- c(varnames[i],paste("ff", 1:length(inds), sep=""))
  }

for(i in 1:length(inds))
  {
  print(i)
  load(paste("additive/model",i,"_0.05.Rdata",sep=""))
  for(j in 1:length(varnames))
    {
    plotdata.add1[[j]][,(i+1)] <- predict(model, newdata=preddata[[j]])
    }

  load(paste("additive/model",i,"_0.1.Rdata",sep=""))
  for(j in 1:length(varnames))
    {
    plotdata.add2[[j]][,(i+1)] <- predict(model, newdata=preddata[[j]])
    }

  load(paste("additive/model",i,"_0.5.Rdata",sep=""))
  for(j in 1:length(varnames))
    {
    plotdata.add3[[j]][,(i+1)] <- predict(model, newdata=preddata[[j]])
    }


  load(paste("VCM/model",i,"_0.05.Rdata",sep=""))
  for(j in 1:length(varnames))
    {
    plotdata.vcm1[[j]][,(i+1)] <- predict(model, newdata=preddata[[j]])
    plotdata.vcm1[[length(varnames)+j]][,(i+1)] <- predict(model, newdata=preddatafm[[j]])
    }

  load(paste("VCM/model",i,"_0.1.Rdata",sep=""))
  for(j in 1:length(varnames))
    {
    plotdata.vcm2[[j]][,(i+1)] <- predict(model, newdata=preddata[[j]])
    plotdata.vcm2[[length(varnames)+j]][,(i+1)] <- predict(model, newdata=preddatafm[[j]])
    }

  load(paste("VCM/model",i,"_0.5.Rdata",sep=""))
  for(j in 1:length(varnames))
    {
    plotdata.vcm3[[j]][,(i+1)] <- predict(model, newdata=preddata[[j]])
    plotdata.vcm3[[length(varnames)+j]][,(i+1)] <- predict(model, newdata=preddatafm[[j]])
    }
  }
  
save("pred1", "pred2", "pred3", "plotdata.add1", "plotdata.add2",
    "plotdata.add3", "plotdata.vcm1", "plotdata.vcm2", "plotdata.vcm3",
    "add1.selected", "add2.selected", "add3.selected", 
    "vcm1.selected", "vcm2.selected", "vcm3.selected", 
    "add1.its", "add2.its", "add3.its", "vcm1.its", "vcm2.its", "vcm3.its", 
    "add1.stabsel", "add2.stabsel", "add3.stabsel", 
    "vcm1.stabsel", "vcm2.stabsel", "vcm3.stabsel", 
    "varnames", file="plotresults.Rdata")