APTS Design of Experiments 2022

Helicopter experiment

library(apts.doe)
library(FrF2) # for the half-normal plot and Plackett-Burman designs
library(effects) # for factorial effects plots
library(fractional) # to simplify display of Alias matrices

The experiment

Video: School children

Website: Paper helicopter experiment

An example helicopter

We will have four factors with the following levels:

w <- c(0.03, 0.09) # rotor width [0.03, 0.09]m
r <- c(0.07, 0.12) # rotor length [0.07, 0.12]m
t <- c(0.07, 0.12) # tail length [0.07, 0.12]m
d <- c(0.06, 0.12) # paper density [0.06, 0.12] kgm^-2

Full factorial design

All combinations of the two levels of the four factors. Flighttime is simulated using the function flighttime from the apts.doe package (uses a loosely physics-based model - see Shen et al., 2014, JQT)

X <- expand.grid(w = w, r = r, t = t, d = d) # form the factorial design by taking all 16 combinations
X <- X[sample(1:16), ] # randomise the run order
y <- flighttime(w = X[, 1], r = X[, 2], t = X[, 3], d = X[, 4], theta1 = 1.9) # simulate the data
y

 [1]  1.69025135  3.66005203  0.18266155  0.74274840  0.01720983  9.65204331  1.04978019  2.30328202  6.11983986
[10]  6.77142980  1.93361926  0.72682565 25.72002448  0.97101753 36.72006858  6.69900887

Modelling

Set up the four columns of X as factors, and recode the levels as -1, 1

```r
X <- lapply(X, factor)
for(i in 1:4) contrasts(X[[i]]) <- contr.twolevel() # function from apts.doe to recode levels
heli <- data.frame(X, y)

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Fit the linear model, and get the coefficients


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuY29lZihoZWxpLmxtKVxuXG5gYGAifQ== -->

```r
coef(heli.lm)

(Intercept)          w1          r1          t1          d1       w1:r1       w1:t1       w1:d1       r1:t1 
  6.5599914  -4.2674372   5.3532943  -1.0485400  -3.8050235  -3.5351687   0.8945174   2.4684047  -0.7893593 
      r1:d1       t1:d1    w1:r1:t1    w1:r1:d1    w1:t1:d1    r1:t1:d1 w1:r1:t1:d1 
 -3.0964263   0.5989819   0.7377477   2.1343005  -0.5471642   0.4579221  -0.4257895 

Produce a half-normal plot to decide which effects are important (which estimated effects look like outliers)

par(pty = "s", mar = c(8, 4, 1, 2))
DanielPlot(heli.lm, main = "", datax = F, half = T, autolab = F) # from FrF2

Produce effects plots for important main effects and interactions

helip.lm <- lm(y ~ (.) ^ 3, data = heli)
plot(Effect("w", helip.lm), main = "", rug = F, aspect = 1)

plot(Effect("r", helip.lm), main = "", rug = F, aspect = 1)

plot(Effect("d", helip.lm), main = "", rug = F, aspect = 1)

## Effect is from effects package

plot(Effect(c("w", "r"), helip.lm), main = "", rug = F, x.var = "w")

plot(Effect(c("w", "d"), helip.lm), main = "", rug = F, x.var = "w")

plot(Effect(c("r", "d"), helip.lm), main = "", rug = F, x.var = "r")

plot(Effect(c("w", "r", "d"), helip.lm), main = "", rug = F, x.var = "w")

Half replicate fractional factorial

What happens if we only look at a half-fraction?

coef(heliFF1.lm) 

(Intercept)          w1          r1          t1          d1       w1:r1       w1:t1       w1:d1       r1:t1 
   6.134202   -3.809515    4.806130    1.085761   -3.067276   -2.936187   -2.201909    1.679045          NA 
      r1:d1       t1:d1    w1:r1:t1    w1:r1:d1    w1:t1:d1    r1:t1:d1 w1:r1:t1:d1 
         NA          NA          NA          NA          NA          NA          NA 

Compare the estimated coefficients with the above for the full factorial - e.g. betahat1(frac) = betahat1(full) + betahat234(full)

A <- Amatrix(model.matrix(heliFF1.lm)[, 1:8], model.matrix(heliFF1.lm)[, 9:16])
# Amatrix is a simple function from apts.doe for calculating the alias matrix
fractional(A)

            r1:t1 r1:d1 t1:d1 w1:r1:t1 w1:r1:d1 w1:t1:d1 r1:t1:d1 w1:r1:t1:d1
(Intercept) .     .     .     .        .        .        .        1          
w1          .     .     .     .        .        .        1        .          
r1          .     .     .     .        .        1        .        .          
t1          .     .     .     .        1        .        .        .          
d1          .     .     .     1        .        .        .        .          
w1:r1       .     .     1     .        .        .        .        .          
w1:t1       .     1     .     .        .        .        .        .          
w1:d1       1     .     .     .        .        .        .        .          

coef(heli.lm)[1:8] + A %*% coef(heli.lm)[-(1:8)]

                 [,1]
(Intercept)  6.134202
w1          -3.809515
r1           4.806130
t1           1.085761
d1          -3.067276
w1:r1       -2.936187
w1:t1       -2.201909
w1:d1        1.679045

coef(heliFF1.lm)

(Intercept)          w1          r1          t1          d1       w1:r1       w1:t1       w1:d1       r1:t1 
   6.134202   -3.809515    4.806130    1.085761   -3.067276   -2.936187   -2.201909    1.679045          NA 
      r1:d1       t1:d1    w1:r1:t1    w1:r1:d1    w1:t1:d1    r1:t1:d1 w1:r1:t1:d1 
         NA          NA          NA          NA          NA          NA          NA 

coef(heli.lm)

(Intercept)          w1          r1          t1          d1       w1:r1       w1:t1       w1:d1       r1:t1 
  6.5599914  -4.2674372   5.3532943  -1.0485400  -3.8050235  -3.5351687   0.8945174   2.4684047  -0.7893593 
      r1:d1       t1:d1    w1:r1:t1    w1:r1:d1    w1:t1:d1    r1:t1:d1 w1:r1:t1:d1 
 -3.0964263   0.5989819   0.7377477   2.1343005  -0.5471642   0.4579221  -0.4257895 

Pick out the other half-fraction with x1x2x3x4 = -1; equivalent to assigning x4 = -x1x2x3

```r
heli.FF2 <- heli[X.full[, 16] == -1,]
heliFF2.lm <- lm(y ~ (.) ^ 4, data = heli.FF2)
coef(heliFF2.lm) 

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiKEludGVyY2VwdCkgICAgICAgICAgdzEgICAgICAgICAgcjEgICAgICAgICAgdDEgICAgICAgICAgZDEgICAgICAgdzE6cjEgICAgICAgdzE6dDEgICAgICAgdzE6ZDEgICAgICAgcjE6dDEgXG4gICA2Ljk3NzQxNSAgIC00Ljg0Mzc5MyAgICA1Ljg3NDkwMSAgIC0zLjEyNzU3NCAgIC00LjQ4MDc2NyAgIC00LjIyOTA3NSAgICAzLjk5ODUyMCAgICAzLjMxMTgzNiAgICAgICAgICBOQSBcbiAgICAgIHIxOmQxICAgICAgIHQxOmQxICAgIHcxOnIxOnQxICAgIHcxOnIxOmQxICAgIHcxOnQxOmQxICAgIHIxOnQxOmQxIHcxOnIxOnQxOmQxIFxuICAgICAgICAgTkEgICAgICAgICAgTkEgICAgICAgICAgTkEgICAgICAgICAgTkEgICAgICAgICAgTkEgICAgICAgICAgTkEgICAgICAgICAgTkEgXG4ifQ== -->

(Intercept) w1 r1 t1 d1 w1:r1 w1:t1 w1:d1 r1:t1 6.977415 -4.843793 5.874901 -3.127574 -4.480767 -4.229075 3.998520 3.311836 NA r1:d1 t1:d1 w1:r1:t1 w1:r1:d1 w1:t1:d1 r1:t1:d1 w1:r1:t1:d1 NA NA NA NA NA NA NA

<!-- rnb-output-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


Compare the estimated coefficients with the above for the full factorial 
- e.g. betahat1(frac) = betahat1(full) - betahat234(full) 


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuQSA8LSBBbWF0cml4KG1vZGVsLm1hdHJpeChoZWxpRkYyLmxtKVssIDE6OF0sIG1vZGVsLm1hdHJpeChoZWxpRkYyLmxtKVssIDk6MTZdKVxuZnJhY3Rpb25hbChBKVxuYGBgXG5gYGAifQ== -->

```r
```r
A <- Amatrix(model.matrix(heliFF2.lm)[, 1:8], model.matrix(heliFF2.lm)[, 9:16])
fractional(A)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiICAgICAgICAgICAgcjE6dDEgcjE6ZDEgdDE6ZDEgdzE6cjE6dDEgdzE6cjE6ZDEgdzE6dDE6ZDEgcjE6dDE6ZDEgdzE6cjE6dDE6ZDFcbihJbnRlcmNlcHQpICAuICAgICAuICAgICAuICAgICAuICAgICAgICAuICAgICAgICAuICAgICAgICAuICAgICAgIC0xICAgICAgICAgXG53MSAgICAgICAgICAgLiAgICAgLiAgICAgLiAgICAgLiAgICAgICAgLiAgICAgICAgLiAgICAgICAtMSAgICAgICAgLiAgICAgICAgIFxucjEgICAgICAgICAgIC4gICAgIC4gICAgIC4gICAgIC4gICAgICAgIC4gICAgICAgLTEgICAgICAgIC4gICAgICAgIC4gICAgICAgICBcbnQxICAgICAgICAgICAuICAgICAuICAgICAuICAgICAuICAgICAgIC0xICAgICAgICAuICAgICAgICAuICAgICAgICAuICAgICAgICAgXG5kMSAgICAgICAgICAgLiAgICAgLiAgICAgLiAgICAtMSAgICAgICAgLiAgICAgICAgLiAgICAgICAgLiAgICAgICAgLiAgICAgICAgIFxudzE6cjEgICAgICAgIC4gICAgIC4gICAgLTEgICAgIC4gICAgICAgIC4gICAgICAgIC4gICAgICAgIC4gICAgICAgIC4gICAgICAgICBcbncxOnQxICAgICAgICAuICAgIC0xICAgICAuICAgICAuICAgICAgICAuICAgICAgICAuICAgICAgICAuICAgICAgICAuICAgICAgICAgXG53MTpkMSAgICAgICAtMSAgICAgLiAgICAgLiAgICAgLiAgICAgICAgLiAgICAgICAgLiAgICAgICAgLiAgICAgICAgLiAgICAgICAgIFxuIn0= -->

        r1:t1 r1:d1 t1:d1 w1:r1:t1 w1:r1:d1 w1:t1:d1 r1:t1:d1 w1:r1:t1:d1

(Intercept) . . . . . . . -1
w1 . . . . . . -1 .
r1 . . . . . -1 . .
t1 . . . . -1 . . .
d1 . . . -1 . . . .
w1:r1 . . -1 . . . . .
w1:t1 . -1 . . . . . .
w1:d1 -1 . . . . . . .

<!-- rnb-output-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuY29lZihoZWxpLmxtKVsxOjhdICsgQSAlKiUgY29lZihoZWxpLmxtKVstKDE6OCldXG5gYGBcbmBgYCJ9 -->

```r
```r
coef(heli.lm)[1:8] + A %*% coef(heli.lm)[-(1:8)]

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiICAgICAgICAgICAgICAgICBbLDFdXG4oSW50ZXJjZXB0KSAgNi45Nzc0MTVcbncxICAgICAgICAgIC00Ljg0Mzc5M1xucjEgICAgICAgICAgIDUuODc0OTAxXG50MSAgICAgICAgICAtMy4xMjc1NzRcbmQxICAgICAgICAgIC00LjQ4MDc2N1xudzE6cjEgICAgICAgLTQuMjI5MDc1XG53MTp0MSAgICAgICAgMy45OTg1MjBcbncxOmQxICAgICAgICAzLjMxMTgzNlxuIn0= -->

             [,1]

(Intercept) 6.977415 w1 -4.843793 r1 5.874901 t1 -3.127574 d1 -4.480767 w1:r1 -4.229075 w1:t1 3.998520 w1:d1 3.311836

<!-- rnb-output-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->

## Plackett-Burman design

Now pick out the runs corresponding to a 12 run Plackett-Burman design

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxucGIxMiA8LSBwYigxMiwgNCkgIyBmcm9tIEZyRjJcbnBiMTIgPC0gc2FwcGx5KHBiMTIsIGZ1bmN0aW9uKHgpIGFzLm51bWVyaWMoYXMuY2hhcmFjdGVyKHgpKSlcbmhlbGkuUEIgPC0gaGVsaVttYXRjaChkYXRhLmZyYW1lKHQocGIxMikpLCBkYXRhLmZyYW1lKHQoWC5mdWxsWywyOjVdKSkpLCBdXG5oZWxpLlBCXG5gYGAifQ== -->

```r
pb12 <- pb(12, 4) # from FrF2
pb12 <- sapply(pb12, function(x) as.numeric(as.character(x)))
heli.PB <- heli[match(data.frame(t(pb12)), data.frame(t(X.full[,2:5]))), ]
heli.PB

Fit the model, and get the coefficients

```r
heliPB.lm <- lm(y ~ (.) ^ 4, data = heli.PB)
coef(heliPB.lm)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiKEludGVyY2VwdCkgICAgICAgICAgdzEgICAgICAgICAgcjEgICAgICAgICAgdDEgICAgICAgICAgZDEgICAgICAgdzE6cjEgICAgICAgdzE6dDEgICAgICAgdzE6ZDEgICAgICAgcjE6dDEgXG4gIDUuNDgyMTE1NCAgLTMuMzQ4NDkzNCAgIDQuMzc5NjAxMSAgLTEuNjMyMjc0NSAgLTIuOTg1NDY3NiAgLTMuNzA2MjUxNiAgLTAuNTIyNjg1NiAgIDIuNDQzMDQ2OSAgIDAuNjI2NTExMCBcbiAgICAgIHIxOmQxICAgICAgIHQxOmQxICAgIHcxOnIxOnQxICAgIHcxOnIxOmQxICAgIHcxOnQxOmQxICAgIHIxOnQxOmQxIHcxOnIxOnQxOmQxIFxuIC0zLjAyNTkwNTYgIC0wLjk3MjQ3NTggICAgICAgICAgTkEgICAgICAgICAgTkEgICAgICAgICAgTkEgICAgICAgICAgTkEgICAgICAgICAgTkEgXG4ifQ== -->

(Intercept) w1 r1 t1 d1 w1:r1 w1:t1 w1:d1 r1:t1 5.4821154 -3.3484934 4.3796011 -1.6322745 -2.9854676 -3.7062516 -0.5226856 2.4430469 0.6265110 r1:d1 t1:d1 w1:r1:t1 w1:r1:d1 w1:t1:d1 r1:t1:d1 w1:r1:t1:d1 -3.0259056 -0.9724758 NA NA NA NA NA

<!-- rnb-output-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


And then the alias matrix; notice all the partial aliasing between the main effect of a factor and all the two-factor interactions not involving that factor


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuWC5QQiA8LSBtb2RlbC5tYXRyaXgoaGVsaVBCLmxtKVxuQSA8LSBBbWF0cml4KFguUEJbLCAxOjVdLCBYLlBCWywgNjoxNl0pIFxuZnJhY3Rpb25hbChBKVxuYGBgXG5gYGAifQ== -->

```r
```r
X.PB <- model.matrix(heliPB.lm)
A <- Amatrix(X.PB[, 1:5], X.PB[, 6:16]) 
fractional(A)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiICAgICAgICAgICAgdzE6cjEgdzE6dDEgdzE6ZDEgcjE6dDEgcjE6ZDEgdDE6ZDEgdzE6cjE6dDEgdzE6cjE6ZDEgdzE6dDE6ZDEgcjE6dDE6ZDEgdzE6cjE6dDE6ZDFcbihJbnRlcmNlcHQpICAgIC4gICAgIC4gICAgIC4gICAgIC4gICAgIC4gICAgIC4gIC0xLzMgICAgIC0xLzMgICAgICAxLzMgICAgIC0xLzMgICAgIC0xLzMgICAgICAgXG53MSAgICAgICAgICAgICAuICAgICAuICAgICAuICAtMS8zICAtMS8zICAgMS8zICAgICAuICAgICAgICAuICAgICAgICAuICAgICAtMS8zICAgICAtMS8zICAgICAgIFxucjEgICAgICAgICAgICAgLiAgLTEvMyAgLTEvMyAgICAgLiAgICAgLiAgLTEvMyAgICAgLiAgICAgICAgLiAgICAgLTEvMyAgICAgICAgLiAgICAgIDEvMyAgICAgICBcbnQxICAgICAgICAgIC0xLzMgICAgIC4gICAxLzMgICAgIC4gIC0xLzMgICAgIC4gICAgIC4gICAgIC0xLzMgICAgICAgIC4gICAgICAgIC4gICAgIC0xLzMgICAgICAgXG5kMSAgICAgICAgICAtMS8zICAgMS8zICAgICAuICAtMS8zICAgICAuICAgICAuICAtMS8zICAgICAgICAuICAgICAgICAuICAgICAgICAuICAgICAtMS8zICAgICAgIFxuIn0= -->

        w1:r1 w1:t1 w1:d1 r1:t1 r1:d1 t1:d1 w1:r1:t1 w1:r1:d1 w1:t1:d1 r1:t1:d1 w1:r1:t1:d1

(Intercept) . . . . . . -1/3 -1/3 1/3 -1/3 -1/3
w1 . . . -1/3 -1/3 1/3 . . . -1/3 -1/3
r1 . -1/3 -1/3 . . -1/3 . . -1/3 . 1/3
t1 -1/3 . 1/3 . -1/3 . . -1/3 . . -1/3
d1 -1/3 1/3 . -1/3 . . -1/3 . . . -1/3

<!-- rnb-output-end -->

<!-- rnb-chunk-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuaGVsaVBCMi5sbSA8LSBsbSh5IH4gKC4pLCBkYXRhID0gaGVsaS5QQilcbmNvZWYoaGVsaVBCMi5sbSlcbmBgYFxuYGBgIn0= -->

```r
```r
heliPB2.lm <- lm(y ~ (.), data = heli.PB)
coef(heliPB2.lm)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiKEludGVyY2VwdCkgICAgICAgICAgdzEgICAgICAgICAgcjEgICAgICAgICAgdDEgICAgICAgICAgZDEgXG4gICA1LjQ4MjExNSAgIC0yLjg3Mjg1NCAgICA0LjA2MzYzOSAgICAxLjQyNjEyNyAgIC0yLjEzMzExNiBcbiJ9 -->

(Intercept) w1 r1 t1 d1 5.482115 -2.872854 4.063639 1.426127 -2.133116

<!-- rnb-output-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuY29lZihoZWxpLmxtKVsxOjVdICsgQSAlKiUgY29lZihoZWxpLmxtKVstKDE6NSldXG5gYGBcbmBgYCJ9 -->

```r
```r
coef(heli.lm)[1:5] + A %*% coef(heli.lm)[-(1:5)]

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiICAgICAgICAgICAgICAgICBbLDFdXG4oSW50ZXJjZXB0KSAgNS40ODIxMTVcbncxICAgICAgICAgIC0yLjg3Mjg1NFxucjEgICAgICAgICAgIDQuMDYzNjM5XG50MSAgICAgICAgICAgMS40MjYxMjdcbmQxICAgICAgICAgIC0yLjEzMzExNlxuIn0= -->

             [,1]

(Intercept) 5.482115 w1 -2.872854 r1 4.063639 t1 1.426127 d1 -2.133116 ```