Paper Helicopter experiment, part II

Last week I created a JAGS model combining data from two paper helicopter datasets. This week, I will use the model to find the longest flying one.

Predicting

The JAGS/RJAGS system has no predict() function that I know of. What I therefore did is adapt the model so during estimation of the parameters the predictions were made. Using this adapted model, two prediction steps were made.
In step one predictions from the whole design space were combined. To keep the number of predictions at least somewhat limited, only a few levels were used for the continuous variables. This step was used to select the best region within  the whole space. Step two focuses on the best region and provides more detailed predictions.

Step 1 

After predicting from the whole experimental space, the mean and lower 5% limits of predicted times were plotted.

It was decided to focus on the region top right. At least 2.7 for the lower 5% limit, and at least 3.7 for the mean time. The associated settings are summarized below.
        PaperType    WingLength       BodyWidth       BodyLength     TapedBody
 bond        :72   Min.   : 7.408   Min.   :2.540   Min.   : 3.810   No :114  
 regular1    :72   1st Qu.:12.065   1st Qu.:3.387   1st Qu.: 6.562   Yes: 54  
 construction:24   Median :12.065   Median :4.233   Median : 6.562            
                   Mean   :11.871   Mean   :4.163   Mean   : 6.955            
                   3rd Qu.:12.065   3rd Qu.:5.080   3rd Qu.: 9.313            
                   Max.   :12.065   Max.   :5.080   Max.   :12.065            
 TapedWing PaperClip PaperClip2 Fold     test       Time        
 No : 68   No :84    No: 20     No:168   WH:168   Mode:logical  
 Yes:100   Yes:84    WH:148                       NA's:168      
                     RH:  0                                     
                                                                
                                                                
                                                                
      Mean            u95             l05       
 Min.   :3.203   Min.   :3.574   Min.   :2.701  
 1st Qu.:3.327   1st Qu.:3.808   1st Qu.:2.776  
 Median :3.465   Median :3.975   Median :2.847  
 Mean   :3.486   Mean   :4.108   Mean   :2.873  
 3rd Qu.:3.636   3rd Qu.:4.388   3rd Qu.:2.952  
 Max.   :3.877   Max.   :5.044   Max.   :3.165  

Phase 2

The second prediction only varied PaperType, BodyWidth, BodyLength and TapedWing. All others were set at their most occurring setting. As can be seen, there is a bit of a trade-off. It is possible to select the longest time, but that incurs some chance of a much lower time, because of model uncertainty. On the other hand, for a slightly lesser mean time, we can have the certainty.

It is my choice to avoid the more uncertain region. Hence I will base my choice on the lower limit. Here we can see that there is a tradeoff. The bond paper needs a slightly longer BodyLength, while Regular paper can have a short BodyLength. BodyWidth should be maximized, but that is not a sensitive parameter.

For completeness, the mean prediction. This shows hardly any interaction. Hence the need for higher BodyLength in bond type paper is due to lack of experiments in this region. A few confirming final experiments seem to be in order. Within those, we could also include a low BodyWidth, since the models are unclear if this should be maximized or minimized.

Code used

Code for actual data are in previous post. This code starts after reading in those data.
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)
library(ggplot2)

helispred <- expand.grid(
    PaperType=c('bond','regular1','construction'),
    WingLength=seq(min(helis$WingLength),max(helis$WingLength),length.out=4),
    BodyWidth=seq(min(helis$BodyWidth),max(helis$BodyWidth),length.out=4),
    BodyLength=seq(min(helis$BodyLength),max(helis$BodyLength),length.out=4),
    TapedBody=c('No','Yes'),
    TapedWing=c('No','Yes'),
    PaperClip=c('No','Yes'),
    PaperClip2=c('No','WH','RH'),
    Fold='No',
    test='WH',
    Time=NA)

helisboth <- rbind(helis,helispred)


#################################
datain <- list(
    PaperType=c(2,1,3,1)[helisboth$PaperType],
    WingLength=helisboth$WingLength,
    BodyLength=helisboth$BodyLength,
    BodyWidth=helisboth$BodyWidth,
    PaperClip=c(1,2,3)[helisboth$PaperClip2],
    TapedBody=c(0,1)[helisboth$TapedBody],
    TapedWing=c(0,1)[helisboth$TapedWing],
    test=c(1,2)[helisboth$test],
    Time=helisboth$Time,
    n=nrow(helis),
    m=nrow(helispred))

parameters <- c('Mul','WL','BL','PT','BW','PC','TB','TW','StDev',
    'WLBW','WLPC',            'WLWL',
    'BLPT'       ,'BLPC',     'BLBL',
    'BWPC',                   'BWBW',  'other','pred')

jmodel <- function() {
  for (i in 1:(n+m)) {     
    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]       
          )
  }
  for (i in 1:n) {
    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)
  for (i in 1:m) {
    pred[i] <- mu[i+n]
  }
}

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

predmat <- jj$BUGSoutput$sims.matrix[,grep('pred',dimnames(jj$BUGSoutput$sims.matrix)[[2]],value=TRUE)]
helispred$Mean <- colMeans(predmat)
helispred$u95 <- apply(predmat,2,function(x) quantile(x,.95))
helispred$l05 <- apply(predmat,2,function(x) quantile(x,.05))
png('select1.png')
qplot(y=Mean,x=l05,data=helispred)
dev.off()
select <- helispred[helispred$Mean>3.2 & helispred$l05>2.7,]
summary(select)

########

helispred <- expand.grid(
    PaperType=c('bond','regular1'),
    WingLength=12.065,
    BodyWidth=seq(2.5,5,length.out=11),
    BodyLength=seq(3.8,12,length.out=11),
    TapedBody=c('No'),
    TapedWing=c('No','Yes'),
    PaperClip='No',
    PaperClip2=c('WH'),
    Fold='No',
    test='WH',
    Time=NA)

helisboth <- rbind(helis,helispred)
datain <- list(
    PaperType=c(2,1,3,1)[helisboth$PaperType],
    WingLength=helisboth$WingLength,
    BodyLength=helisboth$BodyLength,
    BodyWidth=helisboth$BodyWidth,
    PaperClip=c(1,2,3)[helisboth$PaperClip2],
    TapedBody=c(0,1)[helisboth$TapedBody],
    TapedWing=c(0,1)[helisboth$TapedWing],
    test=c(1,2)[helisboth$test],
    Time=helisboth$Time,
    n=nrow(helis),
    m=nrow(helispred))

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

predmat <- jj$BUGSoutput$sims.matrix[,grep('pred',dimnames(jj$BUGSoutput$sims.matrix)[[2]],value=TRUE)]
helispred$Mean <- colMeans(predmat)
helispred$u95 <- apply(predmat,2,function(x) quantile(x,.95))
helispred$l05 <- apply(predmat,2,function(x) quantile(x,.05))

#
png('select2.png')
qplot(y=Mean,x=l05,data=helispred)
dev.off()


png('l05.png')
v <- ggplot(helispred, aes(BodyLength, BodyWidth, z = l05))
v + stat_contour(aes(colour= ..level.. )) +
    scale_colour_gradient(name='Time' )+
    facet_grid(PaperType ~ TapedWing )+
    ggtitle('Lower 95% predicion') 
dev.off()

png('mean.png')

v <- ggplot(helispred, aes(BodyLength, BodyWidth, z = Mean))
v + stat_contour(aes(colour= ..level.. )) +
    scale_colour_gradient(name='Time' )+
    facet_grid(PaperType ~ TapedWing ) +
    ggtitle('Mean prediction')
dev.off()