R-scripts/bootstrap.R

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