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

U-boats in WW-II

This is the time when we celebrate the end of the second world war in The Netherlands, so I thought to do somwthig with data from that era. One of the things I enjoyed were books on the sea warfare, such as 'The Cruel Sea'by Nicholas Monsarrat. In that book it was told that in 1943 the tide turned, and hence I started to look for data which could support that narrative.

Data

Data which I found are from www.uboat.net, a site with quite a lot of info on the topic. I extracted two data sets from the website, ships hit by U-boats and U-boats lost. In addition on wikipedia I found a table with ships lost, broken down by category. Wikipedia also has a good entry for the battle of the Atlantic.

Plots

The number of allied boats lost increases quite a lot in 1942. What happened is the USA entering the war, but having no protection against German U-boats. The second half of 1942 has technological advances for Germany. In 1943 to quote Wikipedia The Battle of the Atlantic was won by the Allies in two months. There was no single reason for this; what had changed was a sudden convergence of technologies, combined with an increase in Allied resources'. The figures show number of ships hit (data from uboat.net) and in Gross Registered Tons (GRT). The pattern is reasonable similar.

The counter plot is U-boats lost. Here wikipedia and uboat.net have slightly different numbers. What can be seen is an increase in losses in 1942. Then there is 1943 a big peak with over 40 U-boats lost.
To quote wikipedia: In all, 43 U-boats were destroyed in May, 34 in the Atlantic. This was 25% of German U-boat arm (U-Bootwaffe) (UBW)'s total operational strength. This was the turning point in the battle for the Atlantic.

Code
Sys.setlocale(category='LC_TIME',locale='C')
library(ggplot2)

# www.uboat.net
data1 <- read.delim(text=
'year Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1939 2 5 1 1
1940 2 5 3 4 1 1 2 2 1 1 2 0
1941 0 0 5 2 1 4 0 4 2 2 5 10
1942 3 2 6 3 4 3 12 9 10 16 13 5
1943 7 18 16 16 40 17 38 24 11 26 19 7
1944 15 20 25 21 23 25 21 31 21 10 8 14
1945 11 22 27 52 22 ')
data1 <- reshape(data1,
    direction='long',
    varying=list(names(data1)[-1]),
    idvar='year',
    timevar='month',
    v.names='Uboat',
    times=names(data1)[-1])
data1$date <- as.Date(paste(data1$year,data1$month,'01',sep='-'),format='%Y-%b-%d')


data2 <- read.delim(text=
'year Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1939 52 37 30 47
1940 58 54 26 11 17 66 43 66 64 72 36 50
1941 15 47 50 48 66 65 21 32 60 51 17 28
1942 66 82 99 89 146 146 109 131 116 120 142 76
1943 49 88 131 57 49 27 62 30 25 31 16 17
1944 20 28 22 12 17 22 25 37 15 10 11 25
1945 24 24 21 23 6
')
data2 <- reshape(data2,
    direction='long',
    varying=list(names(data2)[-1]),
    idvar='year',
    timevar='month',
    v.names='Convoy',
    times=names(data2)[-1])
data2$date <- as.Date(paste(data2$year,data2$month,'01',sep='-'),format='%Y-%b-%d')

both <- merge(data1,data2)
both <- both[order(both$date),]
both <- both[complete.cases(both),]
qplot(y=Convoy,
    x=date,
    data=both,
    geom='line',
    xlab='',
    ylab='Convoy ships hit')
r2 <- read.delim(text=
'Monthyear "Sunk by U-Boat" "Sunk by aircraft" "Sunk by warship or raider" "Sunk by mines" "German submarines lost"
Sep 1939 153879 0 5051 29537 2
Oct 1939 134807 0 32058 29490 5
Nov 1939 51589 0 1722 120958 1
Dec 1939 80881 2949 22506 82712 1
Jan 1940 111263 23693 0 77116 1
Feb 1940 169566 853 1761 54740 6
Mar 1940 62781 8694 0 35051 1
Apr 1940 32467 13409 5358 19799 5
May 1940 55580 158348 6893 47716 1
Jun 1940 284113 105193 61857 86087 0
Jul 1940 195825 70193 80796 33598 2
Aug 1940 267618 53283 63350 11433 3
Sep 1940 295335 56328 96288 8269 1
Oct 1940 352407 8752 32134 32548 1
Nov 1940 146613 66438 123671 46672 2
Dec 1940 212590 14890 55728 54331 0
Jan 1941 126782 78597 80796 17107 0
Feb 1941 196783 89305 89096 16507 0
Mar 1941 243020 113314 138906 23585 5
Apr 1941 249375 323454 91579 24888 2
May 1941 325492 146302 15002 23194 1
Jun 1941 310143 61414 17759 15326 4
Jul 1941 94209 9275 5792 8583 0
Aug 1941 80310 23862 24897 1400 4
Sep 1941 202820 40812 22910 14948 2
Oct 1941 156554 35222 3305 19737 2
Nov 1941 62196 23015 17715 1714 5
Dec 1941 124070 72850 6661 63853 10
Jan 1942 327357 57086 3275 10079 3
Feb 1942 476451 133746 0 7242 2
Mar 1942 537980 55706 25614 16862 6
Apr 1942 431664 82924 131188 15002 3
May 1942 607247 59041 19363 18795 4
Jun 1942 700235 54769 48474 19936 3
Jul 1942 476065 74313 54358 8905 11
Aug 1942 544410 60532 50516 0 9
Sep 1942 485413 57526 24388 0 10
Oct 1942 619417 5686 7576 5157 16
Nov 1942 729160 53868 19178 992 13
Dec 1942 330816 4853 12312 1618 4
Jan 1943 203128 25503 7040 18475 6
Feb 1943 359328 75 4858 34153 19
Mar 1943 627377 65128 0 884 15
Apr 1943 327943 3034 1742 11961 15
May 1943 264853 20942 0 1568 41
Jun 1943 97753 6083 17655 4334 17
Jul 1943 242145 106005 7176 72 37
Aug 1943 86579 14133 0 19 25
Sep 1943 118841 22905 9977 4396 10
Oct 1943 97407 22680 0 19774 26
Nov 1943 66585 62452 8538 6666 19
Dec 1943 86967 75471 0 6086 8
Jan 1944 92278 24237 6420 7176 16
Feb 1944 92923 21616 2085 0 20
Mar 1944 142944 0 7840 7176 25
Apr 1944 62149 19755 0 0 21
May 1944 24424 2873 0 0 24
Jun 1944 57875 9008 1812 24654 25
Jul 1944 63351 0 7219 8114 24
Aug 1944 98729 0 7176 7194 35
Sep 1944 43368 0 0 1437 22
Oct 1944 7176 0 0 4492 13
Nov 1944 29592 7247 1141 0 8
Dec 1944 58518 35920 0 35612 15
Jan 1945 56988 7176 2365 16368 14
Feb 1945 65233 7177 3899 18076 21
Mar 1945 65077 0 3968 36064 33
Apr 1945 72439 22822 0 8733 53
May 1945 11439 7176 0 0 35
')
r2$date <- as.Date(paste(r2$Monthyear,'01',sep='-'),format='%b %Y-%d')
qplot(y=Sunk.by.U.Boat,
    x=date,
    data=r2,
    geom='line',
    xlab='',
    ylab='Sunk by U Boat (GRT)')

qplot(y=Uboat,
    x=date,
    data=both,
    geom='line',
    xlab='',
    ylab='U boats sunk')+
    geom_line(aes(y=German.submarines.lost,x=date),data=r2)

Should I use premium Diesel? Result: No

A while ago I had a post: 'Should I use premium Diesel? Setup. Since that time the data has been acquired. This post describes the results.

Data

Data is registered by me in 2014 and 2015. 2014 has standard Diesel, while 2015 has premium. Both are from the same months of the year, so should have approximately the same weather and traffic.

Results

The figure shows the main result: While the 2015 has some better results, it is not by much. It is not in the ballpark of 5% which was originally defined.

Results from distribution fitting and t.test confirm the visual observation.

fitdistr

Standard
      mean          sd    
  3.59517971   0.19314598 
 (0.04828649) (0.03414371)
Premium 
      mean          sd    
  3.54015573   0.23905047 
 (0.05976262) (0.04225855)

t.test

Welch Two Sample t-test

data:  usage by kind
t = -0.6934, df = 28.732, p-value = 0.4936
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.2173822  0.1073343
sample estimates:
 mean in group Premium mean in group Standard 
              3.540156               3.595180 

Posterior density for premium Diesel

Previously I defined a prior for the analysis of the data. This was a mixed distribution from three parts. Like Standard fuel, 5% better and something else in the neighborhood of 3.6. Given the simplicity of the model and the non-standard prior I used MCMCpack's MCMCmetrop1R function. This allowed me to write it all as R code.
The result is a density for mean usage quite close to 3.6.

Iterations = 501:20500
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 20000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

       Mean      SD  Naive SE Time-series SE
[1,] 3.5409 0.06413 0.0004535       0.001722
[2,] 0.2699 0.05353 0.0003785       0.001448

2. Quantiles for each variable:

       2.5%   25%    50%    75%  97.5%
var1 3.4147 3.494 3.5465 3.5884 3.6523
var2 0.1873 0.232 0.2629 0.2998 0.3953

Conclusion

While premium Diesel had a few occasions with remarkable fuel usage, it did not result in a relevant decrease in fuel usage overall. In addition, while the advised prices for standard and premium are 5% apart, the achievable price difference is 10%. The real price fighters do not do premium, hence the pressure to discount is less. This 10% price increase is more than I initially expected. Hence the data do not show benefit for using premium Diesel.

Data and code

library(MASS)

library(ggplot2)

library(MCMCpack)

r1 <- read.table(text='l km

28.25 710.6

22.93 690.4

28.51 760.5

23.22 697.9

31.52 871.2

24.68 689.6

30.85 826.9

23.04 699

29.96 845.3

30.16 894.7

25.71 696

23.6 669.8

28.57 739

27.23 727.4

18.31 499.9

24.28 689.5',header=TRUE)

r1$usage=100*r1$l/r1$km

r1$Fuel='Standard'

r2 <- read.table(text='l km

23.97 696

25.33 693.5

24.19 698.5

23.27 689.6

25.78 707.8

15.84 510.1

32.23 860.2

23.29 674.6

25.33 678.5

26.44 686

36.43 913.2

26.46 850.1

24.38 678.8

34.85 1001.3

28.34 858.5

25.85 698.8'

,header=TRUE)

r2$usage=100*r2$l/r2$km

r2$Fuel='Premium'

rr <- rbind(r1,r2)

rr$Fuel <- factor(rr$Fuel)

ggplot(rr,aes(Fuel,usage))+geom_violin()+scale_y_continuous('usage (l/100 km)')

fitdistr(r1$usage,'normal')

fitdistr(r2$usage,'normal')

t.test(usage ~ kind, data=rr)

priormeandens <- function(usage) {

  dens <- (dnorm(usage,3.6,.05)+

      dnorm(usage,3.6/1.05,.05)+

      dnorm(usage,(3.6+3.6/1.05)/2,.15))/3

   log(dens)

 }

 mydens <- function(theta,x) {

   musage=theta[1]

   sdusage=theta[2]

   datdens <- sum(dnorm(x,musage,sdusage,log=TRUE))

   priormean <- priormeandens(musage)

   priorsd <- dunif(sdusage,0,1,log=TRUE)

   datdens+priormean+priorsd

 } 

mcmc1 <- MCMCmetrop1R(mydens,

    theta.init=c(3.8,.2),

    seed=as.integer((Sys.time())),

    x=r2$usage)

plot(mcmc1)

summary(mcmc1)