Unreplicated Two-Way ANOVA Model

Adhesive Example

> adhesive <− data.frame(strength=scan(),

+     pressure=rep(seq(from=120,to=150,by=10),3),

+     temperature=rep(seq(250,270,10),c(4,4,4)))

1: 9.60 9.69 8.43 9.98 11.28 10.10 11.01 10.44

9: 9.00 9.57 9.03 9.80

13:

Read 12 items

The session below 'simulate' hand calculations. It's used only for demonstration.

> adhesive.xtabs <− xtabs(strength~pressure+temperature,data=adhesive)

> adhesive.xtabs

        temperature
pressure   250   260   270
     120  9.60 11.28  9.00
     130  9.69 10.10  9.57
     140  8.43 11.01  9.03
     150  9.98 10.44  9.80

> margin.table(adhesive.xtabs, 1) # Ti.

pressure
  120   130   140   150 
29.88 29.36 28.47 30.22 

> margin.table(adhesive.xtabs, 2) # T.j

temperature
  250   260   270 
37.70 42.83 37.40 

> sum(adhesive.xtabs) # T..

[1] 117.93

> sum(adhesive.xtabs^2) # USS

[1] 1166.349

> sum(margin.table(adhesive.xtabs, 1)^2) # SofTi.^2

[1] 3478.613

> sum(margin.table(adhesive.xtabs,2)^2) # SofT.j^2

[1] 4654.459

> SST <− sum(adhesive.xtabs^2) - sum(adhesive.xtabs)^2/12

> SST

[1] 7.392225

> SSA <− sum(margin.table(adhesive.xtabs, 1)^2)/3 -

+      sum(adhesive.xtabs)^2/12

> SSA

[1] 0.5806917

> SSB <− sum(margin.table(adhesive.xtabs,2)^2)/4 -

+      sum(adhesive.xtabs)^2/12

> SSB

[1] 4.65765

> SSAB <− SST - SSA - SSB

> SSAB

[1] 2.153883

> print( Abar <− rowMeans(adhesive.xtabs) ) # ybar-i.

      120       130       140       150 
 9.960000  9.786667  9.490000 10.073333 

> print( Bbar <− colMeans(adhesive.xtabs) ) # ybar-.j

    250     260     270 
 9.4250 10.7075  9.3500

> print( ylim <− range(pretty(range(Abar,Bbar,ybar))) )

[1]  9.2 10.8

> windows(width=6,height=3,pointsize=10)

> par(mfrow=c(1,2), mar=c(4,4,.5,.5)+.1, oma=c(2,0,3,0))

> plot(seq(120,150,10),Abar,ylim=ylim,

+      xlab=bquote(plain("pressure (")*italic("lb/in")^2*plain(")")),

+      ylab="mean strength",xaxt="n",type="l",lwd=2)

> axis(1, at=seq(120,150,10))

> abline(h=ybar, lty=3, col="darkorange4")

> box(which="figure",col="grey75")

> plot(seq(250,270,10),Bbar,ylim=ylim,

+      xlab=bquote(plain("temperature (")*degree*italic("F")*plain(")")),

+      ylab="mean strength",xaxt="n",type="l",lwd=2)

> axis(1, at=seq(250,270,10))

> abline(h=ybar, lty=3, col="darkorange4")

> box(which="figure",col="grey75")

> mtext(" Adhesive Example", side=3,line=0.5,adj=0, outer=T)

> mtext(bquote(bar(y)[..]==.(ybar)), side=1,line=0.5,outer=T)

> title(main="Main Effect Plots",outer=T)

> box(which="outer", col="grey50")

The main effect plots were produced by above codes.

(adhesive01.png here)

> Alpha <− Abar - ( ybar <− mean(adhesive.xtabs) )

> Alpha # A effects

        120         130         140         150 
 0.13250000 -0.04083333 -0.33750000  0.24583333 

> Beta <− Bbar - ybar

> Beta # B effects

    250     260     270 
-0.4025  0.8800 -0.4775

> num <− sum( (Alpha %o% Beta) * adhesive.xtabs )

> num # numerator of gamma.hat

[1] -0.3321467

> denom <− sum(Alpha^2) * sum(Beta^2)

> denom # denominator of gamma.hat

[1] 0.2253882

> num / denom # gamma.hat

[1] -1.473665

> print( SSgamma <− num^2 / denom ) # SSγ

[1] 0.489473

> print( SSrem <− SSAB - SSgamma ) # SSrem

[1] 1.664410

> print( MSrem <− SSrem / 5 ) # MSrem

[1] 0.3328821

> print( Fgamma <− SSgamma / MSrem ) # Fγ

[1] 1.470409

> pf(Fgamma, 1, 5, lower=F) # Pvalγ

[1] 0.2794440

> cols <− paste("dark",c("magenta","green","blue","orange4"),sep="")

> windows(width=6,height=3,pointsize=10)

> par(mfrow=c(1,2),mar=c(4,4,.5,3)+.1,oma=c(2,0,3,0),bty="l")

> with(data=adhesive, {

+   interaction.plot(pressure,temperature,strength,type="b",

+     col=cols[1:3], pch=1:3)

+   box(which="figure",col="grey80")

+   interaction.plot(temperature,pressure,strength,type="b",

+     col=cols[1:4], pch=1:4)

+   box(which="figure",col="grey80")

+ } )

> mtext("Adhesive Example",side=1, outer=T, adj=0, line=0.5)

> title(main="Interaction Plots", outer=T)

> box(which="outer",col="grey50")

The above codes generate the interaction plots below:

(adhesive02.png here)

> adhesive.aov <− aov(strength~factor(pressure)+factor(temperature),

+     data=adhesive)

> summary(adhesive.aov)

                    Df Sum Sq Mean Sq F value  Pr(>F)  
factor(pressure)     3 0.5807  0.1936  0.5392 0.67270  
factor(temperature)  2 4.6576  2.3288  6.4873 0.03162 *
Residuals            6 2.1539  0.3590                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> require(multcomp)

Loading required package: multcomp
Loading required package: mvtnorm
[1] TRUE

> adhesive.tukey <− glht(adhesive.aov,

+     linfct=mcp("factor(temperature)"="Tukey"))

> summary(adhesive.tukey)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = strength ~ factor(pressure) + factor(temperature), 
    data = adhesive)

Linear Hypotheses:
               Estimate Std. Error t value p value  
260 - 250 == 0   1.2825     0.4237   3.027  0.0526 .
270 - 250 == 0  -0.0750     0.4237  -0.177  0.9829  
270 - 260 == 0  -1.3575     0.4237  -3.204  0.0423 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
(Adjusted p values reported)

> confint(adhesive.tukey,level=.9)

         Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = strength ~ factor(pressure) + factor(temperature), 
    data = adhesive)

Estimated Quantile = 2.5165

Linear Hypotheses:
               Estimate lwr     upr    
260 - 250 == 0  1.2825   0.2164  2.3486
270 - 250 == 0 -0.0750  -1.1411  0.9911
270 - 260 == 0 -1.3575  -2.4236 -0.2914

90% family-wise confidence level

> confint(adhesive.tukey)

         Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = strength ~ factor(pressure) + factor(temperature), 
    data = adhesive)

Estimated Quantile = 3.068

Linear Hypotheses:
               Estimate lwr      upr     
260 - 250 == 0  1.28250 -0.01731  2.58231
270 - 250 == 0 -0.07500 -1.37481  1.22481
270 - 260 == 0 -1.35750 -2.65731 -0.05769

95% family-wise confidence level

> windows(width=5,height=3,pointsize=10)

> par(mar=c(5,5,4,1)+.1)

> plot(adhesive.tukey, sub="Adhesive Data",

+     ylim=c(.5,3.5))

> mtext("Tukey's Method for Factor temperature",side=3,line=0.5)

The codes above generate a plot of Tukey's 95% confidence intervals.

(adhesive03.png here)

> summary(lm(strength~factor(pressure)+poly(temperature,2),

+         data=adhesive))

Call:
lm(formula = strength ~ factor(pressure) + poly(temperature, 
    2), data = adhesive)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6575 -0.4902  0.1233  0.3067  0.6400 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)             9.9600     0.3459  28.793 1.16e-07 ***
factor(pressure)130    -0.1733     0.4892  -0.354   0.7352    
factor(pressure)140    -0.4700     0.4892  -0.961   0.3738    
factor(pressure)150     0.1133     0.4892   0.232   0.8245    
poly(temperature, 2)1  -0.1061     0.5991  -0.177   0.8653    
poly(temperature, 2)2  -2.1556     0.5991  -3.598   0.0114 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 0.5991 on 6 degrees of freedom
Multiple R-Squared: 0.7086,     Adjusted R-squared: 0.4658 
F-statistic: 2.918 on 5 and 6 DF,  p-value: 0.1124