[This article was first published on Wiekvoet, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
The paper helicopter is one of the devices to explain about design of experiments. The aim is to create the longest flying paper helicopter by means of experimental design.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Paper helicopters are a nice example, because they are cheap to make, easy to test landing time and sufficient variables to make it non obvious.
Rather than make and measure my own helicopters, I decided to use data from the internet. In this post I use data from williamghunter.net and http://www.rose-hulman.edu. There is more data on the internet, but these two are fairly similar. Both use a fractional factorial design of 16 runs and they have the same variables. However, a quick check showed that these were different results and, very important, the aliasing structure was different.
Data
Data were taken from the above given locations. Rather than using the coded units, the data was converted to sizes in cm. Time to land was converted to seconds.Since these were separate experiments, it has to be assumed that they used different paper, different heights to drop helicopters from. It even seems, that different ways were found to attach a paperclip to the helicopters.
Simple analysis
To confirm the data an analysis on coded units was performed. These results were same as given by the websites, results not shown here. My own analysis starts with real world units and is by regression. Disadvantage or real world units is that one cannot compare the size of the effects, however, given the designs used, the t-statistic can be used for this purpose.The first data set shows WingLength and BodyLength to have the largest effects.
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.92798 0.54903 3.512 0.009839 **
PaperTyperegular1 -0.12500 0.13726 -0.911 0.392730
WingLength 0.17435 0.03088 5.646 0.000777 ***
BodyLength -0.08999 0.03088 -2.914 0.022524 *
BodyWidth 0.01312 0.07205 0.182 0.860634
PaperClipYes 0.05000 0.13726 0.364 0.726403
FoldYes -0.10000 0.13726 -0.729 0.489918
TapedBodyYes -0.15000 0.13726 -1.093 0.310638
TapedWingYes 0.17500 0.13726 1.275 0.242999
The second data set shows WingLength, PaperClip and PaperType to have the largest effects.
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.73200 0.21737 3.368 0.01196 *
PaperTyperegular2 0.28200 0.06211 4.541 0.00267 **
WingLength 0.16654 0.01223 13.622 2.7e-06 ***
BodyLength -0.02126 0.01630 -1.304 0.23340
BodyWidth -0.03307 0.04890 -0.676 0.52058
PaperClipYes -0.35700 0.06211 -5.748 0.00070 ***
FoldYes 0.04500 0.06211 0.725 0.49222
TapedBodyYes -0.14700 0.06211 -2.367 0.04983 *
TapedWingYes 0.06600 0.06211 1.063 0.32320
Combined analysis
The combination analysis is programmed in Jags. To capture a different falling distance, a Mul parameter is used, which defines a multiplicative effect between the two experiments. In addition, both sets have their own measurement error. There are four types of paper, two from each data set, and three levels of paperclip, no paperclip assumed same for both experiments. In addition to the parameters given earlier, residuals are estimated, in order to have some idea about the quality of fit.
The model then, looks like this
jmodel <- function() {
for (i in 1:n) {
premul[i] <- (test[i]==1)+Mul*(test[i]==2)
mu[i] <- premul[i] * (
WL*WingLength[i]+
BL*BodyLength[i] +
PT[PaperType[i]] +
BW*BodyWidth[i] +
PC[PaperClip[i]] +
FO*Fold[i]+
TB*TapedBody[i]+
TW*TapedWing[i]
)
Time[i] ~ dnorm(mu[i],tau[test[i]])
residual[i] <- Time[i]-mu[i]
}
for (i in 1:2) {
tau[i] <- pow(StDev[i],-2)
StDev[i] ~dunif(0,3)
}
for (i in 1:4) {
PT[i] ~ dnorm(PTM,tauPT)
}
tauPT <- pow(sdPT,-2)
sdPT ~dunif(0,3)
PTM ~dnorm(0,0.01)
WL ~dnorm(0,0.01)
BL ~dnorm(0,1000)
BW ~dnorm(0,1000)
PC[1] <- 0
PC[2]~dnorm(0,0.01)
PC[3]~dnorm(0,0.01)
FO ~dnorm(0,1000)
TB ~dnorm(0,0.01)
TW ~dnorm(0,0.01)
Mul ~ dnorm(1,1) %_% I(0,2)
}
for (i in 1:n) {
premul[i] <- (test[i]==1)+Mul*(test[i]==2)
mu[i] <- premul[i] * (
WL*WingLength[i]+
BL*BodyLength[i] +
PT[PaperType[i]] +
BW*BodyWidth[i] +
PC[PaperClip[i]] +
FO*Fold[i]+
TB*TapedBody[i]+
TW*TapedWing[i]
)
Time[i] ~ dnorm(mu[i],tau[test[i]])
residual[i] <- Time[i]-mu[i]
}
for (i in 1:2) {
tau[i] <- pow(StDev[i],-2)
StDev[i] ~dunif(0,3)
}
for (i in 1:4) {
PT[i] ~ dnorm(PTM,tauPT)
}
tauPT <- pow(sdPT,-2)
sdPT ~dunif(0,3)
PTM ~dnorm(0,0.01)
WL ~dnorm(0,0.01)
BL ~dnorm(0,1000)
BW ~dnorm(0,1000)
PC[1] <- 0
PC[2]~dnorm(0,0.01)
PC[3]~dnorm(0,0.01)
FO ~dnorm(0,1000)
TB ~dnorm(0,0.01)
TW ~dnorm(0,0.01)
Mul ~ dnorm(1,1) %_% I(0,2)
}
Inference for Bugs model at “C:/Users/Kees/AppData/Local/Temp/Rtmp4o0rhh/model16f468e854ce.txt”, fit using jags,
4 chains, each with 3000 iterations (first 1500 discarded)
n.sims = 6000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
BL -0.029 0.014 -0.056 -0.038 -0.028 -0.019 -0.001 1.001 4400
BW -0.005 0.025 -0.052 -0.023 -0.006 0.011 0.044 1.002 1900
FO 0.005 0.028 -0.050 -0.014 0.005 0.023 0.058 1.001 6000
Mul 1.166 0.149 0.819 1.087 1.176 1.254 1.433 1.028 130
PC[1] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 1
PC[2] 0.066 0.141 -0.208 -0.021 0.061 0.147 0.360 1.002 2300
PC[3] -0.362 0.070 -0.501 -0.404 -0.362 -0.319 -0.225 1.001 6000
PT[1] 1.111 0.397 0.516 0.864 1.059 1.286 2.074 1.021 150
PT[2] 1.019 0.379 0.437 0.783 0.974 1.186 1.925 1.019 160
PT[3] 0.728 0.170 0.397 0.615 0.728 0.840 1.068 1.002 2900
PT[4] 0.991 0.168 0.655 0.885 0.993 1.103 1.309 1.002 1600
StDev[1] 0.133 0.039 0.082 0.108 0.127 0.150 0.225 1.005 540
StDev[2] 0.304 0.075 0.192 0.251 0.292 0.343 0.488 1.003 1300
TB -0.144 0.059 -0.264 -0.181 -0.144 -0.108 -0.025 1.001 4100
TW 0.084 0.059 -0.033 0.045 0.084 0.122 0.203 1.001 4400
WL 0.164 0.013 0.138 0.156 0.164 0.172 0.188 1.004 810
residual[1] 0.174 0.146 -0.111 0.079 0.173 0.268 0.464 1.002 1700
residual[2] 0.466 0.158 0.162 0.361 0.463 0.567 0.780 1.004 730
residual[3] 0.150 0.170 -0.173 0.041 0.147 0.253 0.499 1.003 1100
residual[4] -0.416 0.162 -0.733 -0.523 -0.418 -0.308 -0.099 1.001 3800
residual[5] -0.087 0.168 -0.419 -0.198 -0.084 0.026 0.238 1.005 560
residual[6] -0.085 0.156 -0.397 -0.184 -0.084 0.016 0.221 1.003 1200
residual[7] -0.056 0.159 -0.371 -0.156 -0.055 0.047 0.251 1.003 910
residual[8] -0.203 0.157 -0.527 -0.304 -0.198 -0.100 0.095 1.001 6000
residual[9] 0.150 0.150 -0.139 0.052 0.148 0.247 0.451 1.001 6000
residual[10] 0.103 0.156 -0.200 0.003 0.101 0.206 0.415 1.004 720
residual[11] 0.133 0.160 -0.176 0.027 0.131 0.237 0.454 1.002 2100
residual[12] 0.335 0.177 -0.006 0.218 0.332 0.451 0.689 1.004 830
residual[13] -0.436 0.156 -0.747 -0.536 -0.436 -0.337 -0.128 1.002 2100
residual[14] 0.098 0.162 -0.227 -0.007 0.099 0.205 0.410 1.004 670
residual[15] -0.018 0.160 -0.340 -0.118 -0.015 0.084 0.292 1.003 920
residual[16] -0.127 0.155 -0.441 -0.224 -0.125 -0.027 0.173 1.001 3600
residual[17] 0.037 0.088 -0.135 -0.018 0.037 0.093 0.215 1.002 1600
residual[18] -0.088 0.090 -0.274 -0.141 -0.086 -0.031 0.081 1.002 2500
residual[19] -0.074 0.088 -0.248 -0.129 -0.072 -0.018 0.100 1.002 1900
residual[20] -0.079 0.088 -0.259 -0.133 -0.076 -0.023 0.091 1.001 3800
residual[21] -0.037 0.087 -0.201 -0.093 -0.039 0.016 0.141 1.002 3000
residual[22] 0.051 0.087 -0.128 -0.001 0.053 0.107 0.221 1.001 4800
residual[23] -0.008 0.084 -0.177 -0.061 -0.009 0.046 0.159 1.001 5500
residual[24] 0.129 0.086 -0.047 0.076 0.130 0.185 0.294 1.002 1900
residual[25] 0.196 0.087 0.030 0.141 0.196 0.249 0.370 1.003 1400
residual[26] -0.027 0.084 -0.195 -0.081 -0.026 0.029 0.138 1.001 6000
residual[27] 0.070 0.088 -0.101 0.016 0.070 0.124 0.247 1.001 3700
residual[28] -0.166 0.089 -0.355 -0.221 -0.163 -0.108 0.004 1.001 3700
residual[29] -0.052 0.087 -0.223 -0.107 -0.053 0.002 0.124 1.001 4300
residual[30] 0.039 0.089 -0.139 -0.016 0.038 0.095 0.218 1.002 2500
residual[31] -0.079 0.089 -0.245 -0.135 -0.080 -0.026 0.103 1.002 2300
residual[32] 0.048 0.085 -0.122 -0.006 0.049 0.102 0.214 1.002 2300
deviance -15.555 7.026 -26.350 -20.655 -16.487 -11.540 0.877 1.004 750
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = var(deviance)/2)
pD = 24.6 and DIC = 9.0
DIC is an estimate of expected predictive error (lower deviance is better).
Striking in the results is big residuals, for instance for observations 2, 4 and 13. The residuals for observations 4 and 13 are also big when a similar classical model is used, hence this is a clear indication of some kind of interaction.
4 chains, each with 3000 iterations (first 1500 discarded)
n.sims = 6000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
BL -0.029 0.014 -0.056 -0.038 -0.028 -0.019 -0.001 1.001 4400
BW -0.005 0.025 -0.052 -0.023 -0.006 0.011 0.044 1.002 1900
FO 0.005 0.028 -0.050 -0.014 0.005 0.023 0.058 1.001 6000
Mul 1.166 0.149 0.819 1.087 1.176 1.254 1.433 1.028 130
PC[1] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 1
PC[2] 0.066 0.141 -0.208 -0.021 0.061 0.147 0.360 1.002 2300
PC[3] -0.362 0.070 -0.501 -0.404 -0.362 -0.319 -0.225 1.001 6000
PT[1] 1.111 0.397 0.516 0.864 1.059 1.286 2.074 1.021 150
PT[2] 1.019 0.379 0.437 0.783 0.974 1.186 1.925 1.019 160
PT[3] 0.728 0.170 0.397 0.615 0.728 0.840 1.068 1.002 2900
PT[4] 0.991 0.168 0.655 0.885 0.993 1.103 1.309 1.002 1600
StDev[1] 0.133 0.039 0.082 0.108 0.127 0.150 0.225 1.005 540
StDev[2] 0.304 0.075 0.192 0.251 0.292 0.343 0.488 1.003 1300
TB -0.144 0.059 -0.264 -0.181 -0.144 -0.108 -0.025 1.001 4100
TW 0.084 0.059 -0.033 0.045 0.084 0.122 0.203 1.001 4400
WL 0.164 0.013 0.138 0.156 0.164 0.172 0.188 1.004 810
residual[1] 0.174 0.146 -0.111 0.079 0.173 0.268 0.464 1.002 1700
residual[2] 0.466 0.158 0.162 0.361 0.463 0.567 0.780 1.004 730
residual[3] 0.150 0.170 -0.173 0.041 0.147 0.253 0.499 1.003 1100
residual[4] -0.416 0.162 -0.733 -0.523 -0.418 -0.308 -0.099 1.001 3800
residual[5] -0.087 0.168 -0.419 -0.198 -0.084 0.026 0.238 1.005 560
residual[6] -0.085 0.156 -0.397 -0.184 -0.084 0.016 0.221 1.003 1200
residual[7] -0.056 0.159 -0.371 -0.156 -0.055 0.047 0.251 1.003 910
residual[8] -0.203 0.157 -0.527 -0.304 -0.198 -0.100 0.095 1.001 6000
residual[9] 0.150 0.150 -0.139 0.052 0.148 0.247 0.451 1.001 6000
residual[10] 0.103 0.156 -0.200 0.003 0.101 0.206 0.415 1.004 720
residual[11] 0.133 0.160 -0.176 0.027 0.131 0.237 0.454 1.002 2100
residual[12] 0.335 0.177 -0.006 0.218 0.332 0.451 0.689 1.004 830
residual[13] -0.436 0.156 -0.747 -0.536 -0.436 -0.337 -0.128 1.002 2100
residual[14] 0.098 0.162 -0.227 -0.007 0.099 0.205 0.410 1.004 670
residual[15] -0.018 0.160 -0.340 -0.118 -0.015 0.084 0.292 1.003 920
residual[16] -0.127 0.155 -0.441 -0.224 -0.125 -0.027 0.173 1.001 3600
residual[17] 0.037 0.088 -0.135 -0.018 0.037 0.093 0.215 1.002 1600
residual[18] -0.088 0.090 -0.274 -0.141 -0.086 -0.031 0.081 1.002 2500
residual[19] -0.074 0.088 -0.248 -0.129 -0.072 -0.018 0.100 1.002 1900
residual[20] -0.079 0.088 -0.259 -0.133 -0.076 -0.023 0.091 1.001 3800
residual[21] -0.037 0.087 -0.201 -0.093 -0.039 0.016 0.141 1.002 3000
residual[22] 0.051 0.087 -0.128 -0.001 0.053 0.107 0.221 1.001 4800
residual[23] -0.008 0.084 -0.177 -0.061 -0.009 0.046 0.159 1.001 5500
residual[24] 0.129 0.086 -0.047 0.076 0.130 0.185 0.294 1.002 1900
residual[25] 0.196 0.087 0.030 0.141 0.196 0.249 0.370 1.003 1400
residual[26] -0.027 0.084 -0.195 -0.081 -0.026 0.029 0.138 1.001 6000
residual[27] 0.070 0.088 -0.101 0.016 0.070 0.124 0.247 1.001 3700
residual[28] -0.166 0.089 -0.355 -0.221 -0.163 -0.108 0.004 1.001 3700
residual[29] -0.052 0.087 -0.223 -0.107 -0.053 0.002 0.124 1.001 4300
residual[30] 0.039 0.089 -0.139 -0.016 0.038 0.095 0.218 1.002 2500
residual[31] -0.079 0.089 -0.245 -0.135 -0.080 -0.026 0.103 1.002 2300
residual[32] 0.048 0.085 -0.122 -0.006 0.049 0.102 0.214 1.002 2300
deviance -15.555 7.026 -26.350 -20.655 -16.487 -11.540 0.877 1.004 750
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = var(deviance)/2)
pD = 24.6 and DIC = 9.0
DIC is an estimate of expected predictive error (lower deviance is better).
Striking in the results is big residuals, for instance for observations 2, 4 and 13. The residuals for observations 4 and 13 are also big when a similar classical model is used, hence this is a clear indication of some kind of interaction.
Model with interactions
Adding the most obvious interactions, such as WingLength*TapedBody did not really provide a suitable answer. Indeed, large residuals at observations 4 and 13, which are at opposite sides in the fractional factorial design, can not be resolved with one interaction.
Hence I proceeded with adding all two way interactions. Since this was expected to result in a model without clear estimates, all interactions had a strong prior; mean was 0 and precision (tau) was 1000. This model was subsequently reduced by giving the interactions which clearly differed from 0 a lesser precision while interactions which where clearly zero were removed. During this process the parameter Fold was removed from the parameter set. Finally, quadratic effects were added. There is one additional parameter, other, it has no function in the model, but tells what the properties of the prior for the interactions are. Parameters with a standard deviation less than other have information added from the data.
jmodel <- function() {
for (i in 1:n) {
premul[i] <- (test[i]==1)+Mul*(test[i]==2)
mu[i] <- premul[i] * (
WL*WingLength[i]+
BL*BodyLength[i] +
PT[PaperType[i]] +
BW*BodyWidth[i] +
PC[PaperClip[i]] +
TB*TapedBody[i]+
TW*TapedWing[i]+
WLBW*WingLength[i]*BodyWidth[i]+
WLPC[1]*WingLength[i]*(PaperClip[i]==2)+
WLPC[2]*WingLength[i]*(PaperClip[i]==3)+
BLPT[1]*BodyLength[i]*(PaperType[i]==2)+
BLPT[2]*BodyLength[i]*(PaperType[i]==3)+
BLPC[1]*BodyLength[i]*(PaperClip[i]==2)+
BLPC[2]*BodyLength[i]*(PaperClip[i]==3)+
BWPC[1]*BodyWidth[i]*(PaperClip[i]==2)+
BWPC[2]*BodyWidth[i]*(PaperClip[i]==3) +
WLWL*WingLength[i]*WingLength[i]+
BLBL*BodyLength[i]*BodyLength[i]+
BWBW*BodyWidth[i]*BodyWidth[i]
)
Time[i] ~ dnorm(mu[i],tau[test[i]])
residual[i] <- Time[i]-mu[i]
}
for (i in 1:2) {
tau[i] <- pow(StDev[i],-2)
StDev[i] ~dunif(0,3)
WLPC[i] ~dnorm(0,1)
BLPT[i] ~dnorm(0,1)
BLPC[i] ~dnorm(0,1)
BWPC[i] ~dnorm(0,1)
}
for (i in 1:3) {
PT[i] ~ dnorm(PTM,tauPT)
}
tauPT <- pow(sdPT,-2)
sdPT ~dunif(0,3)
PTM ~dnorm(0,0.01)
WL ~dnorm(0,0.01)
BL ~dnorm(0,0.01)
BW ~dnorm(0,0.01)
PC[1] <- 0
PC[2]~dnorm(0,0.01)
PC[3]~dnorm(0,0.01)
TB ~dnorm(0,0.01)
TW ~dnorm(0,0.01)
WLBW~dnorm(0,1)
WLTW~dnorm(0,1)
WLWL~dnorm(0,1)
BLBL~dnorm(0,1)
BWBW~dnorm(0,1)
other~dnorm(0,1)
Mul ~ dnorm(1,1) %_% I(0,2)
}
Inference for Bugs model at “C:/Users/Kees/AppData/Local/Temp/Rtmp4o0rhh/model16f472b05364.txt”, fit using jags,
5 chains, each with 4000 iterations (first 2000 discarded), n.thin = 2
n.sims = 5000 iterations saved
mu.vect sd.vect 2.5% 25% 50% 75% 97.5% Rhat n.eff
BL 0.021 0.197 -0.367 -0.080 0.027 0.121 0.396 1.021 590
BLBL -0.001 0.015 -0.027 -0.009 -0.003 0.006 0.031 1.015 1200
BLPC[1] -0.099 0.105 -0.295 -0.125 -0.086 -0.053 0.021 1.100 560
BLPC[2] -0.110 0.111 -0.334 -0.134 -0.094 -0.060 0.018 1.130 250
BLPT[1] -0.038 0.190 -0.503 -0.124 0.001 0.069 0.286 1.005 600
BLPT[2] 0.058 0.038 -0.031 0.045 0.063 0.078 0.113 1.063 400
BW -0.430 0.558 -1.587 -0.657 -0.389 -0.143 0.463 1.045 960
BWBW 0.009 0.094 -0.160 -0.031 0.009 0.052 0.176 1.053 1300
BWPC[1] -0.224 0.173 -0.615 -0.295 -0.209 -0.133 0.064 1.011 5000
BWPC[2] -0.093 0.101 -0.285 -0.137 -0.091 -0.044 0.085 1.040 5000
Mul 1.053 0.145 0.680 0.997 1.069 1.139 1.281 1.098 290
PC[1] 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000 1
PC[2] 1.459 2.367 -3.571 0.333 1.565 2.617 6.138 1.019 420
PC[3] 0.401 0.732 -0.619 0.032 0.309 0.629 1.954 1.074 320
PT[1] 1.353 1.437 -1.364 0.556 1.318 2.088 4.128 1.032 480
PT[2] 1.906 1.767 -1.087 0.828 1.726 2.814 5.879 1.013 1300
PT[3] 0.731 1.419 -1.864 -0.058 0.682 1.444 3.535 1.032 520
StDev[1] 0.108 0.082 0.045 0.067 0.088 0.120 0.302 1.023 450
StDev[2] 0.267 0.156 0.122 0.177 0.229 0.301 0.706 1.021 390
TB -0.146 0.051 -0.247 -0.172 -0.145 -0.119 -0.048 1.011 5000
TW 0.086 0.054 -0.007 0.055 0.082 0.112 0.204 1.010 1700
WL 0.209 0.380 -0.496 0.007 0.188 0.394 1.035 1.014 670
WLBW 0.051 0.062 -0.013 0.026 0.043 0.062 0.167 1.159 220
WLPC[1] 0.057 0.210 -0.304 -0.063 0.024 0.152 0.556 1.004 1600
WLPC[2] 0.020 0.027 -0.031 0.010 0.021 0.033 0.066 1.044 2400
WLWL -0.014 0.026 -0.072 -0.026 -0.011 0.001 0.032 1.014 5000
other 0.002 1.007 -1.973 -0.680 0.000 0.684 1.952 1.002 2200
residual[1] 0.227 0.272 -0.178 0.066 0.190 0.334 0.935 1.041 390
residual[2] 0.035 0.231 -0.447 -0.084 0.037 0.160 0.503 1.007 2500
residual[3] 0.026 0.269 -0.404 -0.118 -0.002 0.131 0.587 1.039 430
residual[4] -0.123 0.279 -0.542 -0.276 -0.157 -0.018 0.530 1.053 370
residual[5] -0.046 0.241 -0.535 -0.168 -0.043 0.083 0.422 1.008 5000
residual[6] -0.094 0.241 -0.568 -0.221 -0.095 0.035 0.390 1.005 2600
residual[7] 0.284 0.268 -0.139 0.140 0.263 0.392 0.861 1.046 430
residual[8] 0.018 0.240 -0.460 -0.107 0.022 0.144 0.494 1.006 5000
residual[9] 0.121 0.299 -0.310 -0.042 0.079 0.223 0.827 1.054 300
residual[10] 0.038 0.237 -0.428 -0.086 0.034 0.155 0.518 1.006 3100
residual[11] -0.077 0.251 -0.562 -0.204 -0.073 0.046 0.401 1.020 5000
residual[12] 0.153 0.262 -0.286 0.013 0.133 0.267 0.711 1.035 610
residual[13] -0.024 0.244 -0.466 -0.160 -0.035 0.095 0.493 1.008 5000
residual[14] -0.019 0.244 -0.537 -0.140 -0.013 0.111 0.456 1.006 5000
residual[15] -0.159 0.250 -0.663 -0.281 -0.156 -0.038 0.302 1.026 860
residual[16] 0.034 0.273 -0.531 -0.076 0.056 0.178 0.486 1.037 410
residual[17] 0.001 0.115 -0.185 -0.057 -0.008 0.047 0.232 1.047 890
residual[18] 0.016 0.105 -0.187 -0.038 0.017 0.067 0.211 1.014 3300
residual[19] -0.068 0.108 -0.262 -0.118 -0.068 -0.017 0.127 1.036 5000
residual[20] 0.067 0.114 -0.138 0.017 0.067 0.115 0.270 1.046 4500
residual[21] 0.003 0.117 -0.223 -0.046 0.007 0.057 0.203 1.044 3200
residual[22] -0.004 0.113 -0.202 -0.059 -0.007 0.044 0.211 1.035 2000
residual[23] -0.039 0.134 -0.313 -0.081 -0.023 0.027 0.145 1.097 300
residual[24] 0.009 0.114 -0.197 -0.042 0.009 0.061 0.223 1.039 5000
residual[25] 0.045 0.110 -0.170 -0.005 0.046 0.095 0.248 1.028 5000
residual[26] -0.044 0.108 -0.252 -0.096 -0.043 0.007 0.165 1.024 4000
residual[27] 0.046 0.112 -0.164 -0.005 0.046 0.100 0.264 1.022 3600
residual[28] -0.062 0.115 -0.296 -0.104 -0.053 -0.004 0.112 1.047 1400
residual[29] -0.025 0.143 -0.321 -0.064 -0.006 0.042 0.153 1.110 230
residual[30] -0.016 0.118 -0.228 -0.066 -0.015 0.037 0.196 1.042 1400
residual[31] -0.025 0.115 -0.239 -0.072 -0.021 0.028 0.174 1.047 1300
residual[32] 0.020 0.111 -0.176 -0.033 0.017 0.066 0.233 1.041 2600
deviance -32.864 19.923 -62.354 -46.843 -35.763 -22.807 16.481 1.014 420
For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
DIC info (using the rule, pD = var(deviance)/2)
pD = 196.7 and DIC = 163.8
DIC is an estimate of expected predictive error (lower deviance is better).
Model discussion
This model does not have the big residuals. In addition it seems that some parameters, e.g. WLWL and WLBW have small mean values and small standard deviations. To me this suggests that they are indeed estimated and found to be close to 0. After all, if the data contained no information, their standard deviation would be similar to the prior, which is much larger, as seen from the other parameter.The quadratic effects were added to allow detection of a maximum. There is not much presence of these effects, except perhaps in WingLength (parameter WLWL).
For descriptive purposes, I will leave these parameters in. However, for predictive purposes, it may be better to remove them or shrink them closer to zero.
Given the complex way in which the parameters are chosen, it is very well possible that a different model would be better. In hindsight, I might have used the BMA function to do a more thorough selection. Thus the model needs to be validated some more. Since I found two additional data sets on line, these might be used for this purpose.
Code
h1 <- read.table(sep=’t’,header=TRUE,text=’
PaperType WingLength BodyLength BodyWidth PaperClip Fold TapedBody TapedWing Time
regular1 7.62 7.62 3.175 No No No No 2.5
bond 7.62 7.62 3.175 Yes No Yes Yes 2.9
regular1 12.065 7.62 3.175 Yes Yes No Yes 3.5
bond 12.065 7.62 3.175 No Yes Yes No 2.7
regular1 7.62 12.065 3.175 Yes Yes Yes No 2
bond 7.62 12.065 3.175 No Yes No Yes 2.3
regular1 12.065 12.065 3.175 No No Yes Yes 2.9
bond 12.065 12.065 3.175 Yes No No No 3
regular1 7.62 7.62 5.08 No Yes Yes Yes 2.4
bond 7.62 7.62 5.08 Yes Yes No No 2.6
regular1 12.065 7.62 5.08 Yes No Yes No 3.2
bond 12.065 7.62 5.08 No No No Yes 3.7
regular1 7.62 12.065 5.08 Yes No No Yes 1.9
bond 7.62 12.065 5.08 No No Yes No 2.2
regular1 12.065 12.065 5.08 No Yes No No 3
bond 12.065 12.065 5.08 Yes Yes Yes Yes 3
‘)
h2 <- read.table(sep=’t’,header=TRUE,text=’
PaperType BodyWidth BodyLength WingLength PaperClip Fold TapedBody TapedWing Time
regular2 2.54 3.81 5.08 No No No No 1.74
construction 2.54 3.81 5.08 No Yes Yes Yes 1.296
regular2 3.81 3.81 5.08 Yes No Yes Yes 1.2
construction 3.81 3.81 5.08 Yes Yes No No 0.996
regular2 2.54 7.62 5.08 Yes Yes Yes No 1.056
construction 2.54 7.62 5.08 Yes No No Yes 1.104
regular2 3.81 7.62 5.08 No Yes No Yes 1.668
construction 3.81 7.62 5.08 No No Yes No 1.308
regular2 2.54 3.81 10.16 Yes Yes No Yes 2.46
construction 2.54 3.81 10.16 Yes No Yes No 1.74
regular2 3.81 3.81 10.16 No Yes Yes No 2.46
construction 3.81 3.81 10.16 No No No Yes 2.184
regular2 2.54 7.62 10.16 No No Yes Yes 2.316
construction 2.54 7.62 10.16 No Yes No No 2.208
regular2 3.81 7.62 10.16 Yes No No No 1.98
construction 3.81 7.62 10.16 Yes Yes Yes Yes 1.788
‘)
l1 <- lm(Time ~ PaperType + WingLength + BodyLength + BodyWidth +
PaperClip + Fold + TapedBody + TapedWing, data=h1)
summary(l1)
residuals(l1)
l2 <- lm(Time ~ PaperType + WingLength + BodyLength + BodyWidth +
PaperClip + Fold + TapedBody + TapedWing, data=h2)
summary(l2)
h1$test <- ‘WH’
# WingLength, BodyLength
h2$test <- ‘RH’
#WhingLegnth, PaperClip, PaperType
helis <- rbind(h1,h2)
helis$test <- factor(helis$test)
helis$PaperClip2 <- factor(ifelse(helis$PaperClip==’No’,’No’,as.character(helis$test)),
levels=c(‘No’,’WH’,’RH’))
library(R2jags)
datain <- list(
PaperType=c(1:4)[helis$PaperType],
WingLength=helis$WingLength,
BodyLength=helis$BodyLength,
BodyWidth=helis$BodyWidth,
PaperClip=c(1,2,3)[helis$PaperClip2],
Fold=c(0,1)[helis$Fold],
TapedBody=c(0,1)[helis$TapedBody],
TapedWing=c(0,1)[helis$TapedWing],
test=c(1,2)[helis$test],
Time=helis$Time,
n=nrow(helis))
parameters <- c(‘Mul’,’WL’,’BL’,’PT’,’BW’,’PC’,’FO’,’TB’,’TW’,’StDev’,’residual’)
jmodel <- function() {
for (i in 1:n) {
premul[i] <- (test[i]==1)+Mul*(test[i]==2)
mu[i] <- premul[i] * (
WL*WingLength[i]+
BL*BodyLength[i] +
PT[PaperType[i]] +
BW*BodyWidth[i] +
PC[PaperClip[i]] +
FO*Fold[i]+
TB*TapedBody[i]+
TW*TapedWing[i]
)
Time[i] ~ dnorm(mu[i],tau[test[i]])
residual[i] <- Time[i]-mu[i]
}
for (i in 1:2) {
tau[i] <- pow(StDev[i],-2)
StDev[i] ~dunif(0,3)
}
for (i in 1:4) {
PT[i] ~ dnorm(PTM,tauPT)
}
tauPT <- pow(sdPT,-2)
sdPT ~dunif(0,3)
PTM ~dnorm(0,0.01)
WL ~dnorm(0,0.01)
BL ~dnorm(0,1000)
BW ~dnorm(0,1000)
PC[1] <- 0
PC[2]~dnorm(0,0.01)
PC[3]~dnorm(0,0.01)
FO ~dnorm(0,1000)
TB ~dnorm(0,0.01)
TW ~dnorm(0,0.01)
Mul ~ dnorm(1,1) %_% I(0,2)
}
jj <- jags(model.file=jmodel,
data=datain,
parameters=parameters,
progress.bar=’gui’,
n.chain=4,
n.iter=3000,
inits=function() list(Mul=1.3,WL=0.15,BL=-.08,PT=rep(1,4),
BW=0,PC=c(NA,0,0),FO=0,TB=0,TW=0))
jj
#################################
datain <- list(
PaperType=c(2,1,3,1)[helis$PaperType],
WingLength=helis$WingLength,
BodyLength=helis$BodyLength,
BodyWidth=helis$BodyWidth,
PaperClip=c(1,2,3)[helis$PaperClip2],
TapedBody=c(0,1)[helis$TapedBody],
TapedWing=c(0,1)[helis$TapedWing],
test=c(1,2)[helis$test],
Time=helis$Time,
n=nrow(helis))
parameters <- c(‘Mul’,’WL’,’BL’,’PT’,’BW’,’PC’,’TB’,’TW’,’StDev’,’residual’,
‘WLBW’,’WLPC’, ‘WLWL’,
‘BLPT’ ,’BLPC’, ‘BLBL’,
‘BWPC’, ‘BWBW’, ‘other’)
jmodel <- function() {
for (i in 1:n) {
premul[i] <- (test[i]==1)+Mul*(test[i]==2)
mu[i] <- premul[i] * (
WL*WingLength[i]+
BL*BodyLength[i] +
PT[PaperType[i]] +
BW*BodyWidth[i] +
PC[PaperClip[i]] +
TB*TapedBody[i]+
TW*TapedWing[i]+
WLBW*WingLength[i]*BodyWidth[i]+
WLPC[1]*WingLength[i]*(PaperClip[i]==2)+
WLPC[2]*WingLength[i]*(PaperClip[i]==3)+
BLPT[1]*BodyLength[i]*(PaperType[i]==2)+
BLPT[2]*BodyLength[i]*(PaperType[i]==3)+
BLPC[1]*BodyLength[i]*(PaperClip[i]==2)+
BLPC[2]*BodyLength[i]*(PaperClip[i]==3)+
BWPC[1]*BodyWidth[i]*(PaperClip[i]==2)+
BWPC[2]*BodyWidth[i]*(PaperClip[i]==3) +
WLWL*WingLength[i]*WingLength[i]+
BLBL*BodyLength[i]*BodyLength[i]+
BWBW*BodyWidth[i]*BodyWidth[i]
)
Time[i] ~ dnorm(mu[i],tau[test[i]])
residual[i] <- Time[i]-mu[i]
}
for (i in 1:2) {
tau[i] <- pow(StDev[i],-2)
StDev[i] ~dunif(0,3)
WLPC[i] ~dnorm(0,1)
BLPT[i] ~dnorm(0,1)
BLPC[i] ~dnorm(0,1)
BWPC[i] ~dnorm(0,1)
}
for (i in 1:3) {
PT[i] ~ dnorm(PTM,tauPT)
}
tauPT <- pow(sdPT,-2)
sdPT ~dunif(0,3)
PTM ~dnorm(0,0.01)
WL ~dnorm(0,0.01)
BL ~dnorm(0,0.01)
BW ~dnorm(0,0.01)
PC[1] <- 0
PC[2]~dnorm(0,0.01)
PC[3]~dnorm(0,0.01)
TB ~dnorm(0,0.01)
TW ~dnorm(0,0.01)
WLBW~dnorm(0,1)
WLTW~dnorm(0,1)
WLWL~dnorm(0,1)
BLBL~dnorm(0,1)
BWBW~dnorm(0,1)
other~dnorm(0,1)
Mul ~ dnorm(1,1) %_% I(0,2)
}
jj <- jags(model.file=jmodel,
data=datain,
parameters=parameters,
progress.bar=’gui’,
n.chain=5,
n.iter=4000,
inits=function() list(Mul=1.3,WL=0.15,BL=-.08,PT=rep(1,3),
PC=c(NA,0,0),TB=0,TW=0))
jj
To leave a comment for the author, please follow the link and comment on their blog: Wiekvoet.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.