Unit 4: General Factorial Design

1 Packages and Functions

library(agricolae)
# the "nonadditivity" is provided by this package.
anova.R2 <- function (obj.lm)
{
    obj.anova <- anova (obj.lm)
    sq <- obj.anova[,"Sum Sq"]
    R2 <- sq/sum(sq)
    obj.anova[, "R2"] <- R2
    obj.anova
}

2 Two-factor Factorial Design

An engineer is designing a battery for use in a device that will be subjected to some extreme variations in temperature. The only design parameter that he can select at this point is the plate material for the battery, and he has three possible choices. When the devices is manufactured and is shipped to the filed, the engineer has no control over the temperature extremes that the device will encounter, and he knows from experience that temperature will probably affect the effective battery life. However, temperature can be controlled in the product development laboratory for the purpose of a test. The engineer decides to test all three plate materials at three temperature levels – 15, 90, and 1250F because these temperature levels are consistent with the product end-use environment.

The engineer wants to answer the following questions:

• What effects do material type and temperature have on the life of the battery?
• Is there a choice of material that would give uniformly long life regardless of temperature?

battery <- read.table("data/battery.txt", header=T)

## change level order in temperature for better visualization
battery$temperature <- factor (battery$temperature, levels = c("15F", "70F", "125F"))
attach(battery)
battery

material	temperature	response
M1	15F	130
M1	15F	155
M1	15F	74
M1	15F	180
M1	70F	34
M1	70F	40
M1	70F	80
M1	70F	75
M1	125F	20
M1	125F	70
M1	125F	82
M1	125F	58
M2	15F	150
M2	15F	188
M2	15F	159
M2	15F	126
M2	70F	136
M2	70F	122
M2	70F	106
M2	70F	115
M2	125F	25
M2	125F	70
M2	125F	58
M2	125F	45
M3	15F	138
M3	15F	110
M3	15F	168
M3	15F	160
M3	70F	174
M3	70F	120
M3	70F	150
M3	70F	139
M3	125F	96
M3	125F	104
M3	125F	82
M3	125F	60

## look at factor combination
table (temperature, material)

##            material
## temperature M1 M2 M3
##        15F   4  4  4
##        70F   4  4  4
##        125F  4  4  4

##Visualize the data

par(mfrow=c(1,2))
plot(response ~ material, data=battery)
plot(response ~ temperature, data=battery)

par(mfrow=c(1,1))
interaction.plot(temperature,material,response,col = 1:3, ylim = range (response) )
points (temperature, response, col = material, pch = as.numeric(material))

Tabular summary of the dataset:

tapply(response, INDEX = list (material,temperature), mean) -> battery.means
cbind(battery.means, rowMeans(battery.means)) -> battery.means1
rbind(battery.means1, colMeans(battery.means1)) -> battery.means2
battery.means2

##         15F      70F     125F          
## M1 134.7500  57.2500 57.50000  83.16667
## M2 155.7500 119.7500 49.50000 108.33333
## M3 144.0000 145.7500 85.50000 125.08333
##    144.8333 107.5833 64.16667 105.52778

2.1 Fit Additive Models without Interaction Term

# fit the model without interaction (wrong analysis)
g_additive <- lm(response ~ material+temperature, battery)
interaction.plot(temperature,material,  g_additive$fitted.values, col = 1:3, 
                 ylim = range (response), main = "Without Interaction", type = "b")
points (temperature, response,  col = material,   pch = as.numeric(material))

2.2 Fit Interaction Models using “Effect Parametrization”

In this parametrization, the coefficients for "material1, material2, temperature1, temperature2 are main effects (ie, \(\tau_i,\beta_j\)) in the textbook, and the interaction coefficients are \((\tau\beta)_{ij}\) in the textbook.

g <- lm(response ~ material*temperature, battery, 
        contrasts = list(material = contr.sum, temperature=contr.sum))
interaction.plot(temperature,material,  g$fitted.values, col = 1:3, 
                 ylim = range (response),main = "With Interaction",
                 type = "b")
points (temperature, response, col = material, pch = as.numeric(material))

# look at coefficients
summary(g)

## 
## Call:
## lm(formula = response ~ material * temperature, data = battery, 
##     contrasts = list(material = contr.sum, temperature = contr.sum))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.750 -14.625   1.375  17.937  45.250 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             105.528      4.331  24.367  < 2e-16 ***
## material1               -22.361      6.125  -3.651  0.00111 ** 
## material2                 2.806      6.125   0.458  0.65057    
## temperature1             39.306      6.125   6.418  7.1e-07 ***
## temperature2              2.056      6.125   0.336  0.73975    
## material1:temperature1   12.278      8.662   1.417  0.16778    
## material2:temperature1    8.111      8.662   0.936  0.35735    
## material1:temperature2  -27.972      8.662  -3.229  0.00325 ** 
## material2:temperature2    9.361      8.662   1.081  0.28936    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.98 on 27 degrees of freedom
## Multiple R-squared:  0.7652, Adjusted R-squared:  0.6956 
## F-statistic:    11 on 8 and 27 DF,  p-value: 9.426e-07

# look at model matrix
data.frame(battery, model.matrix(g))

material	temperature	response	X.Intercept.	material1	material2	temperature1	temperature2	material1.temperature1	material2.temperature1	material1.temperature2	material2.temperature2
M1	15F	130	1	1	0	1	0	1	0	0	0
M1	15F	155	1	1	0	1	0	1	0	0	0
M1	15F	74	1	1	0	1	0	1	0	0	0
M1	15F	180	1	1	0	1	0	1	0	0	0
M1	70F	34	1	1	0	0	1	0	0	1	0
M1	70F	40	1	1	0	0	1	0	0	1	0
M1	70F	80	1	1	0	0	1	0	0	1	0
M1	70F	75	1	1	0	0	1	0	0	1	0
M1	125F	20	1	1	0	-1	-1	-1	0	-1	0
M1	125F	70	1	1	0	-1	-1	-1	0	-1	0
M1	125F	82	1	1	0	-1	-1	-1	0	-1	0
M1	125F	58	1	1	0	-1	-1	-1	0	-1	0
M2	15F	150	1	0	1	1	0	0	1	0	0
M2	15F	188	1	0	1	1	0	0	1	0	0
M2	15F	159	1	0	1	1	0	0	1	0	0
M2	15F	126	1	0	1	1	0	0	1	0	0
M2	70F	136	1	0	1	0	1	0	0	0	1
M2	70F	122	1	0	1	0	1	0	0	0	1
M2	70F	106	1	0	1	0	1	0	0	0	1
M2	70F	115	1	0	1	0	1	0	0	0	1
M2	125F	25	1	0	1	-1	-1	0	-1	0	-1
M2	125F	70	1	0	1	-1	-1	0	-1	0	-1
M2	125F	58	1	0	1	-1	-1	0	-1	0	-1
M2	125F	45	1	0	1	-1	-1	0	-1	0	-1
M3	15F	138	1	-1	-1	1	0	-1	-1	0	0
M3	15F	110	1	-1	-1	1	0	-1	-1	0	0
M3	15F	168	1	-1	-1	1	0	-1	-1	0	0
M3	15F	160	1	-1	-1	1	0	-1	-1	0	0
M3	70F	174	1	-1	-1	0	1	0	0	-1	-1
M3	70F	120	1	-1	-1	0	1	0	0	-1	-1
M3	70F	150	1	-1	-1	0	1	0	0	-1	-1
M3	70F	139	1	-1	-1	0	1	0	0	-1	-1
M3	125F	96	1	-1	-1	-1	-1	1	1	1	1
M3	125F	104	1	-1	-1	-1	-1	1	1	1	1
M3	125F	82	1	-1	-1	-1	-1	1	1	1	1
M3	125F	60	1	-1	-1	-1	-1	1	1	1	1

2.3 Fit Interaction Models using “baseline/treatment parametrization”

In this parametrization, the coefficients for “materialM2, materialM3, temperature70F, temperature125F” are NOT difference of main effects to the level 1. The coefficient interpretation is confusing.

g2 <- lm(response ~ material*temperature, battery)
interaction.plot(temperature,material,  g2$fitted.values, col = 1:3, 
                 ylim = range (response),main = "With Interaction",
                 type = "b")
points (temperature, response, col = material, pch = as.numeric(material))

summary(g2)

## 
## Call:
## lm(formula = response ~ material * temperature, data = battery)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -60.750 -14.625   1.375  17.938  45.250 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  134.75      12.99  10.371 6.46e-11 ***
## materialM2                    21.00      18.37   1.143 0.263107    
## materialM3                     9.25      18.37   0.503 0.618747    
## temperature70F               -77.50      18.37  -4.218 0.000248 ***
## temperature125F              -77.25      18.37  -4.204 0.000257 ***
## materialM2:temperature70F     41.50      25.98   1.597 0.121886    
## materialM3:temperature70F     79.25      25.98   3.050 0.005083 ** 
## materialM2:temperature125F   -29.00      25.98  -1.116 0.274242    
## materialM3:temperature125F    18.75      25.98   0.722 0.476759    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 25.98 on 27 degrees of freedom
## Multiple R-squared:  0.7652, Adjusted R-squared:  0.6956 
## F-statistic:    11 on 8 and 27 DF,  p-value: 9.426e-07

data.frame(model.matrix(g2))

X.Intercept.	materialM2	materialM3	temperature70F	temperature125F	materialM2.temperature70F	materialM3.temperature70F	materialM2.temperature125F	materialM3.temperature125F
1	0	0	0	0	0	0	0	0
1	0	0	0	0	0	0	0	0
1	0	0	0	0	0	0	0	0
1	0	0	0	0	0	0	0	0
1	0	0	1	0	0	0	0	0
1	0	0	1	0	0	0	0	0
1	0	0	1	0	0	0	0	0
1	0	0	1	0	0	0	0	0
1	0	0	0	1	0	0	0	0
1	0	0	0	1	0	0	0	0
1	0	0	0	1	0	0	0	0
1	0	0	0	1	0	0	0	0
1	1	0	0	0	0	0	0	0
1	1	0	0	0	0	0	0	0
1	1	0	0	0	0	0	0	0
1	1	0	0	0	0	0	0	0
1	1	0	1	0	1	0	0	0
1	1	0	1	0	1	0	0	0
1	1	0	1	0	1	0	0	0
1	1	0	1	0	1	0	0	0
1	1	0	0	1	0	0	1	0
1	1	0	0	1	0	0	1	0
1	1	0	0	1	0	0	1	0
1	1	0	0	1	0	0	1	0
1	0	1	0	0	0	0	0	0
1	0	1	0	0	0	0	0	0
1	0	1	0	0	0	0	0	0
1	0	1	0	0	0	0	0	0
1	0	1	1	0	0	1	0	0
1	0	1	1	0	0	1	0	0
1	0	1	1	0	0	1	0	0
1	0	1	1	0	0	1	0	0
1	0	1	0	1	0	0	0	1
1	0	1	0	1	0	0	0	1
1	0	1	0	1	0	0	0	1
1	0	1	0	1	0	0	0	1

2.4 ANOVA

# ANOVA table for additive models
anova.R2(g_additive)

	Df	Sum Sq	Mean Sq	F value	Pr(>F)	R2
material	2	10683.72	5341.8611	5.947226	0.0065146	0.1375935
temperature	2	39118.72	19559.3611	21.775920	0.0000012	0.5038023
Residuals	31	27844.53	898.2106	NA	NA	0.3586042

# ANOVA table for interaction models
anova.R2(g)

	Df	Sum Sq	Mean Sq	F value	Pr(>F)	R2
material	2	10683.722	5341.861	7.911372	0.0019761	0.1375935
temperature	2	39118.722	19559.361	28.967692	0.0000002	0.5038023
material:temperature	4	9613.778	2403.444	3.559535	0.0186112	0.1238139
Residuals	27	18230.750	675.213	NA	NA	0.2347902

# ANOVA for g2 for interaction models

anova.R2(g2)

	Df	Sum Sq	Mean Sq	F value	Pr(>F)	R2
material	2	10683.722	5341.861	7.911372	0.0019761	0.1375935
temperature	2	39118.722	19559.361	28.967692	0.0000002	0.5038023
material:temperature	4	9613.778	2403.444	3.559535	0.0186112	0.1238139
Residuals	27	18230.750	675.213	NA	NA	0.2347902

# The interaction SS is the difference of RSS from g_additive to g
anova (g_additive, g)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
31	27844.53	NA	NA	NA	NA
27	18230.75	4	9613.778	3.559535	0.0186112

2.5 Effect Comparison with Tukey Method

TukeyHSD(aov(response ~ material*temperature, battery))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = response ~ material * temperature, data = battery)
## 
## $material
##           diff       lwr      upr     p adj
## M2-M1 25.16667 -1.135677 51.46901 0.0627571
## M3-M1 41.91667 15.614323 68.21901 0.0014162
## M3-M2 16.75000 -9.552344 43.05234 0.2717815
## 
## $temperature
##               diff        lwr       upr     p adj
## 70F-15F  -37.25000  -63.55234 -10.94766 0.0043788
## 125F-15F -80.66667 -106.96901 -54.36432 0.0000001
## 125F-70F -43.41667  -69.71901 -17.11432 0.0009787
## 
## $`material:temperature`
##                    diff         lwr        upr     p adj
## M2:15F-M1:15F     21.00  -40.823184  82.823184 0.9616404
## M3:15F-M1:15F      9.25  -52.573184  71.073184 0.9998527
## M1:70F-M1:15F    -77.50 -139.323184 -15.676816 0.0065212
## M2:70F-M1:15F    -15.00  -76.823184  46.823184 0.9953182
## M3:70F-M1:15F     11.00  -50.823184  72.823184 0.9994703
## M1:125F-M1:15F   -77.25 -139.073184 -15.426816 0.0067471
## M2:125F-M1:15F   -85.25 -147.073184 -23.426816 0.0022351
## M3:125F-M1:15F   -49.25 -111.073184  12.573184 0.2016535
## M3:15F-M2:15F    -11.75  -73.573184  50.073184 0.9991463
## M1:70F-M2:15F    -98.50 -160.323184 -36.676816 0.0003449
## M2:70F-M2:15F    -36.00  -97.823184  25.823184 0.5819453
## M3:70F-M2:15F    -10.00  -71.823184  51.823184 0.9997369
## M1:125F-M2:15F   -98.25 -160.073184 -36.426816 0.0003574
## M2:125F-M2:15F  -106.25 -168.073184 -44.426816 0.0001152
## M3:125F-M2:15F   -70.25 -132.073184  -8.426816 0.0172076
## M1:70F-M3:15F    -86.75 -148.573184 -24.926816 0.0018119
## M2:70F-M3:15F    -24.25  -86.073184  37.573184 0.9165175
## M3:70F-M3:15F      1.75  -60.073184  63.573184 1.0000000
## M1:125F-M3:15F   -86.50 -148.323184 -24.676816 0.0018765
## M2:125F-M3:15F   -94.50 -156.323184 -32.676816 0.0006078
## M3:125F-M3:15F   -58.50 -120.323184   3.323184 0.0742711
## M2:70F-M1:70F     62.50    0.676816 124.323184 0.0460388
## M3:70F-M1:70F     88.50   26.676816 150.323184 0.0014173
## M1:125F-M1:70F     0.25  -61.573184  62.073184 1.0000000
## M2:125F-M1:70F    -7.75  -69.573184  54.073184 0.9999614
## M3:125F-M1:70F    28.25  -33.573184  90.073184 0.8281938
## M3:70F-M2:70F     26.00  -35.823184  87.823184 0.8822881
## M1:125F-M2:70F   -62.25 -124.073184  -0.426816 0.0474675
## M2:125F-M2:70F   -70.25 -132.073184  -8.426816 0.0172076
## M3:125F-M2:70F   -34.25  -96.073184  27.573184 0.6420441
## M1:125F-M3:70F   -88.25 -150.073184 -26.426816 0.0014679
## M2:125F-M3:70F   -96.25 -158.073184 -34.426816 0.0004744
## M3:125F-M3:70F   -60.25 -122.073184   1.573184 0.0604247
## M2:125F-M1:125F   -8.00  -69.823184  53.823184 0.9999508
## M3:125F-M1:125F   28.00  -33.823184  89.823184 0.8347331
## M3:125F-M2:125F   36.00  -25.823184  97.823184 0.5819453

par (mfrow=c(1,3), mar=c(4,4,4,4))
plot(TukeyHSD(aov(response ~ material*temperature, battery)), las=2)

2.5.1 Computing Tukey’s CIs

\[se(\hat \mu_i-\hat \mu_j)=\hat\sigma\times \sqrt{1/n_i+1/n_j}\] where \(\hat\sigma=\sqrt{MSE}\) from the interaction model, where \(n_i\) and \(n_j\) are the sample sizes in \(\hat \mu_i\) and \(\hat \mu_j\). This SE will be multiplied by the t quantile or Tukey quantiles to obtain the LSD CI or family-wise CI. # Comparing Main Effects of Materials

Tabular summary of the dataset:

tapply(response, INDEX = list ( material, temperature), mean) -> battery.means
cbind(battery.means, rowMeans(battery.means)) -> battery.means1
rbind(battery.means1, colMeans(battery.means1)) -> battery.means2
battery.means2

##         15F      70F     125F          
## M1 134.7500  57.2500 57.50000  83.16667
## M2 155.7500 119.7500 49.50000 108.33333
## M3 144.0000 145.7500 85.50000 125.08333
##    144.8333 107.5833 64.16667 105.52778

table(material,temperature)

##         temperature
## material 15F 70F 125F
##       M1   4   4    4
##       M2   4   4    4
##       M3   4   4    4

Computing SE and ME:

se <- sigma (g) * sqrt(1/12+1/12); se

## [1] 10.60827

# number of means to select is 3
ME <- se * qtukey(0.95, 3, g$df.residual) / sqrt(2)

# Comparing means between temperatures

# CI for main effects between temperature 15F to temperature 70F
battery.means2[4,1] -battery.means2[4,2]  + c (-ME, +ME)

## [1] 10.94766 63.55234

# Comparing means between materials

# CI for main effects between material 1 and material 2
battery.means2[1,4] -battery.means2[2,4]  + c (-ME, +ME)

## [1] -51.469011   1.135677

# for comparing means for all combinations of "material" and "temperature"
# 4 is number of replicates
se2 <- sigma (g) * sqrt(1/4+1/4); se2

## [1] 18.37407

ME2 <- se2 * qtukey (0.95, 9, g$df.residual) / sqrt(2)


battery.means[2,3]-battery.means[1,1] + c(-ME2, ME2)

## [1] -147.07318  -23.42682

2.6 Model Checking

par(mfrow=c(1,2))
qqnorm(g$residuals)
qqline(g$residuals)
shapiro.test(g$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  g$residuals
## W = 0.97606, p-value = 0.6117

plot(g$fit,g$res,xlab="fitted vlaue", ylab="Residuals", 
     col = temperature,  pch = as.numeric(material),
     main = "Residuals vs Fitted in Interaction Model"); 
abline (h = 0)

detach ("battery")

3 Tukey’s Nonadditive Test for Factorial Experiments with a Single Replicate

tukeydata <- read.table("data/tukeydata.txt", header=T)
tukeydata

temperature	pressure	y
t100	p25	5
t100	p30	4
t100	p35	6
t100	p40	3
t100	p45	5
t125	p25	3
t125	p30	1
t125	p35	4
t125	p40	2
t125	p45	3
t150	p25	1
t150	p30	1
t150	p35	3
t150	p40	1
t150	p45	2

attach(tukeydata)

## The following object is masked from package:datasets:
## 
##     pressure

par (mfrow=c(1,1))
interaction.plot(temperature, pressure,y, type = "b")

table (temperature, pressure)

##            pressure
## temperature p25 p30 p35 p40 p45
##        t100   1   1   1   1   1
##        t125   1   1   1   1   1
##        t150   1   1   1   1   1

names(tukeydata)

## [1] "temperature" "pressure"    "y"

g<-lm(y~temperature+pressure)
anova(g)

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
temperature	2	23.33333	11.66667	46.66667	0.0000388
pressure	4	11.60000	2.90000	11.60000	0.0020634
Residuals	8	2.00000	0.25000	NA	NA

library(agricolae)
#We need to supply the degree freedom and values of MSE from the additive model, which are 8 and 0.25 respectively. 
nonadditivity(y, temperature, pressure, 8, 0.25)

## 
## Tukey's test of nonadditivity
## y 
## 
## P : 2.666667
## Q : 72.17778 
## 
## Analysis of Variance Table
## 
## Response: residual
##               Df  Sum Sq  Mean Sq F value Pr(>F)
## Nonadditivity  1 0.09852 0.098522  0.3627  0.566
## Residuals      7 1.90148 0.271640

## $P
## Nonadditivity 
##      2.666667 
## 
## $Q
## Nonadditivity 
##      72.17778 
## 
## $ANOVA
## Analysis of Variance Table
## 
## Response: residual
##               Df  Sum Sq  Mean Sq F value Pr(>F)
## Nonadditivity  1 0.09852 0.098522  0.3627  0.566
## Residuals      7 1.90148 0.271640

detach ("tukeydata")