R-scripts/bootstrap.R

   1 ## PsN R script for plotting results of bootstrap\r
   2 ## Justin Wilkins\r
   3 ## October 2005\r
   4 \r
   5 ## set basic options\r
   6 ## autogenerated code begins\r
   7 ## we can set other things, like colours, line types, layout, etc here if we wish\r
   8 \r
   9 ## This script assumes "ofv" is the last column before we get into the parameter estimates\r
  10 \r
  11 min.failed    <- FALSE     # do we want to omit minimization failed runs?\r
  12 cov.failed    <- FALSE      # do we want to omit covariance failed runs?\r
  13 cov.warnings  <- TRUE      # do we want to omit covariance failed runs?\r
  14 boundary      <- FALSE      # do we want to omit boundary runs?\r
  15 showmean      <- FALSE      # show line for mean\r
  16 showquart     <- FALSE      # show line for quartiles\r
  17 \r
  18 ## autogenerated code ends\r
  19 \r
  20 ## read file\r
  21 bootstrap.data <- read.csv("raw_results1.csv", header=T)\r
  22 \r
  23 ## replace underscores\r
  24 for (i in 1:length(names(bootstrap.data))) {\r
  25   names(bootstrap.data)[i] <- gsub("\_", "\.", names(bootstrap.data)[i])\r
  26 }\r
  27 \r
  28 ## find ofv column index\r
  29 index <- 0\r
  30 seen  <- FALSE\r
  31 \r
  32 for (i in names(bootstrap.data)) {\r
  33   if (!seen) {\r
  34     index <- index + 1\r
  35   }\r
  36   if (i == "ofv") {\r
  37     seen <- TRUE\r
  38   }\r
  39 }\r
  40 \r
  41 ## get number of parameters\r
  42 n       <- length(colnames(bootstrap.data)) - index\r
  43 nparams <- length(colnames(bootstrap.data))\r
  44 \r
  45 cols <- c()\r
  46 \r
  47 ## separate out original model fit\r
  48 p1 <- subset(bootstrap.data, bootstrap.data$model != 0)\r
  49 o1 <- subset(bootstrap.data, bootstrap.data$model == 0)\r
  50 \r
  51 #names(p1)[2] <- "minimization.successful"\r
  52 #names(p1)[3] <- "covariance.step.successful"\r
  53 #names(p1)[4] <- "covariance.step.warnings"\r
  54 #names(p1)[5] <- "estimate.near.boundary"\r
  55 \r
  56 #cat(nrow(p1))\r
  57 if (min.failed) {\r
  58   p1 <- subset(p1, minimization.successful == 1)\r
  59 } \r
  60 if (cov.failed) {\r
  61   p1 <- subset(p1, covariance.step.successful == 1)\r
  62 }\r
  63 if (cov.warnings) {\r
  64   p1 <- subset(p1, covariance.step.warnings == 0)\r
  65 }\r
  66 if (boundary) {\r
  67   p1 <- subset(p1, estimate.near.boundary == 0)\r
  68 }\r
  69 \r
  70 ## stats and plots for each- single\r
  71 \r
  72 for (i in index:nparams) {\r
  73   if (mode(p1[[i]]) == "numeric") {\r
  74     sp <- summary(p1[[i]])\r
  75     # IQR <- diff(summary(p1[[i]])[c(5,2)])\r
  76     dp <- density(p1[[i]], na.rm=T)\r
  77     parmlabel <- names(p1)[i]\r
  78   \r
  79     postscript(file=paste("bootstrap.", parmlabel, ".ps", sep=""), paper="a4",\r
  80       title=paste("Bootstrap results - ", parmlabel, sep=""),\r
  81       horizontal=T) \r
  82   \r
  83     qu <- quantile(p1[[i]], c(0.025, 0.975), na.rm=T)\r
  84     \r
  85     hist(p1[[i]], \r
  86           main = paste("Bootstrap results - ", parmlabel, sep=""),\r
  87           xlab = parmlabel,\r
  88           # ylim = c(0, max(dp$y)),\r
  89           # ylim = c(0, 1),\r
  90           xlim = c(min(dp$x), max(dp$x)),\r
  91           breaks = 20,\r
  92           # xlim = c(min(p1[[i]]) - min(p1[[i]])*0.15,max(p1[[i]]) + max(p1[[i]])*0.15),\r
  93           probability = T,\r
  94           sub=paste(paste(paste("n = ", nrow(p1), sep=""), "; Median = ", sp[4], sep=""), "; Orig = ", o1[[i]], sep="") )\r
  95           \r
  96     #h <- hist(p1[[i]], prob=T, plot=F)\r
  97         \r
  98     lines(dp, lwd=2, lty=3, col="red")\r
  99   \r
 100     if (showquart) {\r
 101       abline(v=sp[2], lwd= 1, lty=3, col="red") ## 1st quartile\r
 102       abline(v=sp[5], lwd= 1, lty=3, col="red") ## 3rd quartile\r
 103     }\r
 104     if (showmean) {\r
 105       abline(v=sp[3], lty=2, lwd=1, col="red") ## mean\r
 106     }\r
 107     abline(v=sp[4], lty=1, lwd=2, col="red") ## median\r
 108     abline(v=o1[[i]], lty=2, lwd=1, col="red") ## original\r
 109     abline(v=qu[1], lty=4, lwd=1, col="red") ## 2.5% CL\r
 110     abline(v=qu[2], lty=4, lwd=1, col="red") ## 97.5% CL\r
 111   \r
 112     text(qu[1], max(dp$y), labels=signif(qu[1], digits = 3), cex = .8, adj = c(0,0), pos='2')\r
 113     text(qu[2], max(dp$y), labels=signif(qu[2], digits = 3), cex = .8, adj = c(0,0), pos='4')\r
 114     #text(sp[4], max(h$density), labels=paste("Med: ", signif(sp[4], digits = 3), sep=""), adj = c(-1,0), cex = .8, pos='2')\r
 115   }\r
 116 }\r
 117 \r
 118 ## stats and plots for each - 6 per sheet\r
 119 \r
 120 total  <- 0\r
 121 bspage <- 0\r
 122 \r
 123 for (i in index:nparams) {\r
 124   if (mode(p1[[i]]) == "numeric") {\r
 125     sp <- summary(p1[[i]])\r
 126     # IQR <- diff(summary(p1[[i]])[c(5,2)])\r
 127     dp <- density(p1[[i]], na.rm=T)\r
 128     parmlabel <- names(p1)[i]\r
 129   \r
 130     if (total == 0) {\r
 131       bspage <- bspage + 1\r
 132       postscript(file=paste("bootstrap.page", bspage, ".ps", sep=""), paper="a4",\r
 133         title="Bootstrap results",\r
 134         horizontal=T) \r
 135       par(mfrow = c(3,3))\r
 136     }\r
 137     total <- total + 1\r
 138   \r
 139     qu <- quantile(p1[[i]], c(0.025, 0.975), na.rm=T)\r
 140     \r
 141     hist(p1[[i]], \r
 142           main = paste("Bootstrap results - ", parmlabel, sep=""),\r
 143           xlab = parmlabel,\r
 144           # ylim = c(0, max(dp$y)),\r
 145           # ylim = c(0, 1),\r
 146           xlim = c(min(dp$x), max(dp$x)),\r
 147           breaks = 20,\r
 148           # xlim = c(min(p1[[i]]) - min(p1[[i]])*0.15,max(p1[[i]]) + max(p1[[i]])*0.15),\r
 149           probability = T,\r
 150           sub=paste(paste(paste("n = ", nrow(p1), sep=""), "; Median = ", sp[4], sep=""), "; Orig = ", o1[[i]], sep="") )\r
 151     \r
 152     #h <- hist(p1[[i]], prob=T, plot=F)\r
 153         \r
 154     lines(dp, lwd=2, lty=3, col="red")\r
 155   \r
 156     if (showquart) {\r
 157       abline(v=sp[2], lwd= 1, lty=3, col="red") ## 1st quartile\r
 158       abline(v=sp[5], lwd= 1, lty=3, col="red") ## 3rd quartile\r
 159     }\r
 160     if (showmean) {\r
 161       abline(v=sp[3], lty=2, lwd=1, col="red") ## mean\r
 162     }\r
 163     abline(v=sp[4], lty=1, lwd=2, col="red") ## median\r
 164     abline(v=o1[[i]], lty=2, lwd=1, col="red") ## original\r
 165     abline(v=qu[1], lty=4, lwd=1, col="red") ## 2.5% CL\r
 166     abline(v=qu[2], lty=4, lwd=1, col="red") ## 97.5% CL\r
 167   \r
 168     text(qu[1], max(dp$y), labels=signif(qu[1], digits = 3), cex = .8, adj = c(0,0), pos='2')\r
 169     text(qu[2], max(dp$y), labels=signif(qu[2], digits = 3), cex = .8, adj = c(0,0), pos='4')\r
 170     #text(sp[4], max(h$density), labels=paste("Med: ", signif(sp[4], digits = 3), sep=""), adj = c(-1,0), cex = .8, pos='2')\r
 171   \r
 172     if (total == 9) {\r
 173       total <- 0\r
 174       dev.off()\r
 175     }\r
 176   }\r
 177 }