Sunday, May 17, 2015

The paper helicopter experiment

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.
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
It seems then, that the two experiments show somewhat different effects. WingLength is certainly important. BodyLength maybe. Regarding paper, both have regular paper, but one has bond paper and the other construction paper. It is not difficult to imagine these are quite different.

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)
}

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. 

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

No comments:

Post a Comment