dat <- data.frame(
run = 1:16,
A = c(-1, 1, 1,-1, 1,-1,-1, 1, 1,-1,-1, 1,-1, 1, 1,-1),
B = c(-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1,-1, 1),
C = c(-1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1,-1,-1, 1, 1),
D = c(-1,-1,-1,-1, 1, 1, 1, 1,-1,-1,-1,-1, 1, 1, 1, 1),
E = c(-1,-1,-1,-1,-1,-1,-1,-1, 1, 1, 1, 1, 1, 1, 1, 1),
y1 = c(1330,1935,1770,1275,1880,1385,1220,2155,
1715,1385,1000,1990,1275,1660,1880,1275),
y2 = c(1330,1935,1770,1275,1935,1440,1165,2100,
1715,1550,1165,1990,1660,1605,1935,1220),
y3 = c(1165,1880,1770,1275,1880,1495,1440,2100,
1660,1550,1495,1990,1550,1660,1935,1275)
)
#Run means and within-run variances
dat$ybar <- rowMeans(dat[, c("y1", "y2", "y3")])
dat$s2 <- apply(dat[, c("y1", "y2", "y3")], 1, var)
# for zero variances
minpos <- min(dat$s2[dat$s2 > 0])
dat$s2star <- ifelse(dat$s2 == 0, 0.5 * minpos, dat$s2)
# Dispersion response
dat$z <- log(dat$s2star)
# Two-factor interactions
fac <- c("A", "B", "C", "D", "E")
for(i in 1:(length(fac)-1)){
for(j in (i+1):length(fac)){
dat[[paste0(fac[i], fac[j])]] <- dat[[fac[i]]] * dat[[fac[j]]]
}
}
# replicate columns to long format for location analysis
long <- reshape(dat,
varying = c("y1", "y2", "y3"),
v.names = "y",
timevar = "rep",
times = 1:3,
direction = "long")
# Full location model
loc.full <- lm(y ~ A+B+C+D+E+
AB+AC+AD+AE+BC+BD+BE+CD+CE+DE,
data = long)
summary(loc.full)
##
## Call:
## lm(formula = y ~ A + B + C + D + E + AB + AC + AD + AE + BC +
## BD + BE + CD + CE + DE, data = long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -220.00 -36.67 0.00 18.33 275.00
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.605e+03 1.408e+01 113.990 < 2e-16 ***
## A 2.635e+02 1.408e+01 18.717 < 2e-16 ***
## B 3.667e+01 1.408e+01 2.604 0.013853 *
## C -2.292e+00 1.408e+01 -0.163 0.871733
## D 2.521e+01 1.408e+01 1.790 0.082866 .
## E -1.604e+01 1.408e+01 -1.139 0.263033
## AB 1.146e+01 1.408e+01 0.814 0.421781
## AC 8.250e+01 1.408e+01 5.859 1.63e-06 ***
## AD 6.995e-14 1.408e+01 0.000 1.000000
## AE -4.125e+01 1.408e+01 -2.930 0.006211 **
## BC 2.063e+01 1.408e+01 1.465 0.152727
## BD -5.271e+01 1.408e+01 -3.743 0.000716 ***
## BE -2.979e+01 1.408e+01 -2.116 0.042232 *
## CD 1.375e+01 1.408e+01 0.977 0.336119
## CE 9.167e+00 1.408e+01 0.651 0.519675
## DE -3.667e+01 1.408e+01 -2.604 0.013853 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 97.55 on 32 degrees of freedom
## Multiple R-squared: 0.9313, Adjusted R-squared: 0.8991
## F-statistic: 28.93 on 15 and 32 DF, p-value: 2.168e-14
# Selected location model
loc.mod <- lm(y ~ A + B + AC + AE + BD + BE + DE, data = long)
summary(loc.mod)
##
## Call:
## lm(formula = y ~ A + B + AC + AE + BD + BE + DE, data = long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -277.292 -41.250 -6.875 57.292 217.708
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1605.00 14.21 112.984 < 2e-16 ***
## A 263.54 14.21 18.552 < 2e-16 ***
## B 36.67 14.21 2.581 0.013620 *
## AC 82.50 14.21 5.808 8.81e-07 ***
## AE -41.25 14.21 -2.904 0.005975 **
## BD -52.71 14.21 -3.710 0.000629 ***
## BE -29.79 14.21 -2.097 0.042344 *
## DE -36.67 14.21 -2.581 0.013620 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 98.42 on 40 degrees of freedom
## Multiple R-squared: 0.9126, Adjusted R-squared: 0.8973
## F-statistic: 59.69 on 7 and 40 DF, p-value: < 2.2e-16
# Dispersion model
disp.mod <- lm(z ~ A + B + AB, data = dat)
summary(disp.mod)
##
## Call:
## lm(formula = z ~ A + B + AB, data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3438 -0.5199 0.1733 0.2303 1.5466
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.8076 0.1998 39.086 5.07e-14 ***
## A -1.0648 0.1998 -5.331 0.000179 ***
## B -0.6529 0.1998 -3.268 0.006724 **
## AB 0.6529 0.1998 3.268 0.006724 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.799 on 12 degrees of freedom
## Multiple R-squared: 0.8058, Adjusted R-squared: 0.7572
## F-statistic: 16.59 on 3 and 12 DF, p-value: 0.0001438
#Predicted dispersion for A and B combinations
ABgrid <- expand.grid(A = c(-1, 1),
B = c(-1, 1))
ABgrid$AB <- ABgrid$A * ABgrid$B
ABgrid$zhat <- predict(disp.mod, newdata = ABgrid)
ABgrid
## A B AB zhat
## 1 -1 -1 1 10.178165
## 2 1 -1 -1 6.742767
## 3 -1 1 -1 7.566727
## 4 1 1 1 6.742767
#Generating all 32 factor combinations
grid <- expand.grid(A = c(-1, 1),
B = c(-1, 1),
C = c(-1, 1),
D = c(-1, 1),
E = c(-1, 1))
#interaction columns needed in selected location model
grid$AC <- grid$A * grid$C
grid$AE <- grid$A * grid$E
grid$BD <- grid$B * grid$D
grid$BE <- grid$B * grid$E
grid$DE <- grid$D * grid$E
# Predicted tensile strength
grid$yhat <- predict(loc.mod, newdata = grid)
# Best setting
grid[which.max(grid$yhat), ]
## A B C D E AC AE BD BE DE yhat
## 8 1 1 1 -1 -1 1 -1 -1 -1 1 2074.792