| blob | 0102d081bf6abcecc530aa730f36efbceac8dab4 |
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 }