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, April 2007
   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 excl.id <- c()              # exclude samples that have this individual
  24
  25 ## autogenerated code ends
  26
  27 ## read files
  28 bootstrap.data <- read.csv("raw_results1.csv", header=T)
  29 incl.ids       <- read.csv("included_individuals1.csv", header=F)
  30
  31 ## replace underscores
  32 for (i in 1:length(names(bootstrap.data))) {
  33   names(bootstrap.data)[i] <- gsub("\_", "\.", names(bootstrap.data)[i])
  34 }
  35
  36 ## find ofv column index
  37 index <- 0
  38 seen  <- FALSE
  39
  40 for (i in names(bootstrap.data)) {
  41   if (!seen) {
  42     index <- index + 1
  43   }
  44   if (i == "ofv") {
  45     seen <- TRUE
  46   }
  47 }
  48
  49 ## get number of parameters
  50 n       <- length(colnames(bootstrap.data)) - index
  51 nparams <- length(colnames(bootstrap.data))
  52
  53 ## separate out original model fit
  54 p1 <- subset(bootstrap.data, bootstrap.data$model != 0)
  55 o1 <- subset(bootstrap.data, bootstrap.data$model == 0)
  56
  57 incl.flag <- rep(0,length(rownames(p1)))
  58 for( i in excl.id ) {
  59   incl.flag <- incl.flag + rowSums( incl.ids == i )
  60 }
  61
  62 p1 <- p1[(incl.flag==0),]
  63
  64 #names(p1)[2] <- "minimization.successful"
  65 #names(p1)[3] <- "covariance.step.successful"
  66 #names(p1)[4] <- "covariance.step.warnings"
  67 #names(p1)[5] <- "estimate.near.boundary"
  68
  69 #cat(nrow(p1))
  70 if (min.failed) {
  71   p1 <- subset(p1, minimization.successful == 1)
  72 }
  73 if (cov.failed) {
  74   p1 <- subset(p1, covariance.step.successful == 1)
  75 }
  76 if (cov.warnings) {
  77   p1 <- subset(p1, covariance.step.warnings == 0)
  78 }
  79 if (boundary) {
  80   p1 <- subset(p1, estimate.near.boundary == 0)
  81 }
  82
  83 ## stats and plots for each- single
  84
  85 for (i in index:nparams) {
  86   if (mode(p1[[i]]) == "numeric" &&
  87       sum(p1[[i]],na.rm=T)) {
  88     sp <- summary(p1[[i]])
  89     # IQR <- diff(summary(p1[[i]])[c(5,2)])
  90     dp <- density(p1[[i]], na.rm=T)
  91     parmlabel <- names(p1)[i]
  92
  93     pdf(file=paste("bootstrap.", parmlabel, ".pdf", sep=""), paper="special",
  94       title=paste("Bootstrap results - ", parmlabel, sep=""),width=10,height=7 )
  95
  96     qu <- quantile(p1[[i]], c(0.025, 0.975), na.rm=T)
  97
  98     legend=paste("n = ", nrow(p1), sep="")
  99     if (showmean) {
 100       legend=paste(legend, "; Mean = ", sp[4], sep="")
 101     }
 102     if (showmedian) {
 103       legend=paste(legend, "; Median = ", sp[3], sep="")
 104     }
 105     if (showoriginal) {
 106       legend=paste(legend, "; Orig = ", o1[[i]], sep="")
 107     }
 108
 109     hist(p1[[i]],
 110           main = paste("Bootstrap results - ", parmlabel, sep=""),
 111           xlab = parmlabel,
 112           # ylim = c(0, max(dp$y)),
 113           # ylim = c(0, 1),
 114           xlim = c(min(dp$x), max(dp$x)),
 115           breaks = 20,
 116           # xlim = c(min(p1[[i]]) - min(p1[[i]])*0.15,max(p1[[i]]) + max(p1[[i]])*0.15),
 117           probability = T,
 118           sub=legend )
 119 #          sub=paste(paste(paste("n = ", nrow(p1), sep=""), "; Median = ", sp[4], sep=""), "; Orig = ", o1[[i]], sep="") )
 120
 121     #h <- hist(p1[[i]], prob=T, plot=F)
 122
 123     lines(dp, lwd=2, lty=3, col="red")
 124
 125     if (showquart) {
 126       abline(v=sp[2], lwd= 1, lty=3, col="red") ## 1st quartile
 127       abline(v=sp[5], lwd= 1, lty=3, col="red") ## 3rd quartile
 128     }
 129     if (showmean) {
 130       abline(v=sp[4], lty=2, lwd=1, col="red") ## mean
 131     }
 132     if (showmedian) {
 133       abline(v=sp[3], lty=1, lwd=2, col="red") ## median
 134     }
 135     if (showoriginal) {
 136       abline(v=o1[[i]], lty=2, lwd=1, col="red") ## original
 137     }
 138     if (show95CI) {
 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       text(qu[1], max(dp$y), labels=signif(qu[1], digits = 3), cex = .8, adj = c(0,0), pos='2')
 142       text(qu[2], max(dp$y), labels=signif(qu[2], digits = 3), cex = .8, adj = c(0,0), pos='4')
 143     }
 144 #    abline(v=sp[4], lty=1, lwd=2, col="red") ## median
 145 #    abline(v=o1[[i]], lty=2, lwd=1, col="red") ## original
 146 #    abline(v=qu[1], lty=4, lwd=1, col="red") ## 2.5% CL
 147 #    abline(v=qu[2], lty=4, lwd=1, col="red") ## 97.5% CL
 148
 149 #    text(qu[1], max(dp$y), labels=signif(qu[1], digits = 3), cex = .8, adj = c(0,0), pos='2')
 150 #    text(qu[2], max(dp$y), labels=signif(qu[2], digits = 3), cex = .8, adj = c(0,0), pos='4')
 151     #text(sp[4], max(h$density), labels=paste("Med: ", signif(sp[4], digits = 3), sep=""), adj = c(-1,0), cex = .8, pos='2')
 152   }
 153 }
 154
 155 ## stats and plots for each - 6 per sheet
 156
 157 total  <- 0
 158 bspage <- 0
 159
 160 for (i in index:nparams) {
 161   if (mode(p1[[i]]) == "numeric" &&
 162       sum(p1[[i]],na.rm=T)) {
 163     sp <- summary(p1[[i]])
 164     # IQR <- diff(summary(p1[[i]])[c(5,2)])
 165     dp <- density(p1[[i]], na.rm=T)
 166     parmlabel <- names(p1)[i]
 167
 168     if (total == 0) {
 169       bspage <- bspage + 1
 170       pdf(file=paste("bootstrap.page", bspage, ".pdf", sep=""), paper="special",
 171         title="Bootstrap results",width=10,height=7)
 172       par(mfrow = c(3,3))
 173     }
 174     total <- total + 1
 175
 176     qu <- quantile(p1[[i]], c(0.025, 0.975), na.rm=T)
 177
 178     legend=paste("n = ", nrow(p1), sep="")
 179     if (showmean) {
 180       legend=paste(legend, "; Mean = ", sp[3], sep="")
 181     }
 182     if (showmedian) {
 183       legend=paste(legend, "; Median = ", sp[4], sep="")
 184     }
 185     if (showoriginal) {
 186       legend=paste(legend, "; Orig = ", o1[[i]], sep="")
 187     }
 188
 189     hist(p1[[i]],
 190           main = paste("Bootstrap results - ", parmlabel, sep=""),
 191           xlab = parmlabel,
 192           # ylim = c(0, max(dp$y)),
 193           # ylim = c(0, 1),
 194           xlim = c(min(dp$x), max(dp$x)),
 195           breaks = 20,
 196           # xlim = c(min(p1[[i]]) - min(p1[[i]])*0.15,max(p1[[i]]) + max(p1[[i]])*0.15),
 197           probability = T,
 198           sub=legend )
 199 #=paste(paste(paste("n = ", nrow(p1), sep=""), "; Median = ", sp[4], sep=""), "; Orig = ", o1[[i]], sep="") )
 200
 201     #h <- hist(p1[[i]], prob=T, plot=F)
 202
 203     lines(dp, lwd=2, lty=3, col="red")
 204
 205     if (showquart) {
 206       abline(v=sp[2], lwd= 1, lty=3, col="red") ## 1st quartile
 207       abline(v=sp[5], lwd= 1, lty=3, col="red") ## 3rd quartile
 208     }
 209     if (showmean) {
 210       abline(v=sp[3], lty=2, lwd=1, col="red") ## mean
 211     }
 212     if (showmedian) {
 213       abline(v=sp[4], lty=1, lwd=2, col="red") ## median
 214     }
 215     if (showoriginal) {
 216       abline(v=o1[[i]], lty=2, lwd=1, col="red") ## original
 217     }
 218     if (show95CI) {
 219       abline(v=qu[1], lty=4, lwd=1, col="red") ## 2.5% CL
 220       abline(v=qu[2], lty=4, lwd=1, col="red") ## 97.5% CL
 221       text(qu[1], max(dp$y), labels=signif(qu[1], digits = 3), cex = .8, adj = c(0,0), pos='2')
 222       text(qu[2], max(dp$y), labels=signif(qu[2], digits = 3), cex = .8, adj = c(0,0), pos='4')
 223     }
 224
 225 #    abline(v=sp[4], lty=1, lwd=2, col="red") ## median
 226 #    abline(v=o1[[i]], lty=2, lwd=1, col="red") ## original
 227 #    abline(v=qu[1], lty=4, lwd=1, col="red") ## 2.5% CL
 228 #    abline(v=qu[2], lty=4, lwd=1, col="red") ## 97.5% CL
 229
 230 #    text(qu[1], max(dp$y), labels=signif(qu[1], digits = 3), cex = .8, adj = c(0,0), pos='2')
 231 #    text(qu[2], max(dp$y), labels=signif(qu[2], digits = 3), cex = .8, adj = c(0,0), pos='4')
 232     #text(sp[4], max(h$density), labels=paste("Med: ", signif(sp[4], digits = 3), sep=""), adj = c(-1,0), cex = .8, pos='2')
 233
 234     if (total == 9) {
 235       total <- 0
 236       dev.off()
 237     }
 238   }
 239 }