Two-Way ANOVA Model, #1

Factory Example

> factory <− data.frame(output=scan(),

+   shift=rep(rep(paste("Shift",1:3),c(2,2,2)),2),

+   gender=factor(rep(x<−c("Male","Female"),c(6,6)),

+                 levels=x))

1: 8 6 2 4 4 6 10 12 6 8 6 8

13:

Read 12 items

> factory

   output   shift gender
1       8 Shift 1   Male
2       6 Shift 1   Male
3       2 Shift 2   Male
4       4 Shift 2   Male
5       4 Shift 3   Male
6       6 Shift 3   Male
7      10 Shift 1 Female
8      12 Shift 1 Female
9       6 Shift 2 Female
10      8 Shift 2 Female
11      6 Shift 3 Female
12      8 Shift 3 Female

> levels(factory$gender)

[1] "Male"   "Female"

> x

[1] "Male"   "Female"

# The rhs (right hand side) of the formula in

# next command means all columns in data but the lhs.

print(T.. <− xtabs(output~.,data=factory)) # cell totals

         gender
shift     Male Female
  Shift 1   14     22
  Shift 2    6     14
  Shift 3   10     14

print(nij <− xtabs(rep(1,12)~shift+gender,data=factory)) # cell sizes

         gender
shift     Male Female
  Shift 1    2      2
  Shift 2    2      2
  Shift 3    2      2

T.. / nij # cell means

         gender
shift     Male Female
  Shift 1    7     11
  Shift 2    3      7
  Shift 3    5      7

> x <− paste(factory$gender, unclass(factory$shift), sep=".")

> x

 [1] "Male.1"   "Male.1"   "Male.2"   "Male.2"   "Male.3"   "Male.3"  
 [7] "Female.1" "Female.1" "Female.2" "Female.2" "Female.3" "Female.3"

> unique(x)

[1] "Male.1"   "Male.2"   "Male.3"   "Female.1" "Female.2" "Female.3"

> factory$comb <− factor(x, levels=unique(x))

> factory

   output   shift gender     comb
1       8 Shift 1   Male   Male.1
2       6 Shift 1   Male   Male.1
3       2 Shift 2   Male   Male.2
4       4 Shift 2   Male   Male.2
5       4 Shift 3   Male   Male.3
6       6 Shift 3   Male   Male.3
7      10 Shift 1 Female Female.1
8      12 Shift 1 Female Female.1
9       6 Shift 2 Female Female.2
10      8 Shift 2 Female Female.2
11      6 Shift 3 Female Female.3
12      8 Shift 3 Female Female.3

> levels(factory$comb)

[1] "Male.1"   "Male.2"   "Male.3"   "Female.1" "Female.2" "Female.3"

> with(factory,tapply(output,comb,mean)) # cell means

  Male.1   Male.2   Male.3 Female.1 Female.2 Female.3 
       7        3        5       11        7        7

> summary(aov(output~comb,data=factory)) # Hcells

            Df Sum Sq Mean Sq F value  Pr(>F)  
comb         5 70.667  14.133  7.0667 0.01689 *
Residuals    6 12.000   2.000                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> summary(factory.aov<−aov(output~shift*gender, data=factory)) # ANOVA

             Df Sum Sq Mean Sq F value   Pr(>F)   
shift         2 34.667  17.333  8.6667 0.017003 * 
gender        1 33.333  33.333 16.6667 0.006484 **
shift:gender  2  2.667   1.333  0.6667 0.547708   
Residuals     6 12.000   2.000                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The interaction is insignificant which can be observed in the interaction plots below.

> windows(width=6, # 6in wide

+         height=3, # 3in tall

+         pointsize=10) # 10-point font

> par(mfrow=c(1,2), # 1 by 2 multiple figures, row-wise

+     mar=c(4,4,.5,3)+.1, # four margins (b,l,t,r) for each figure

+     oma=c(2,0,3,0), # outer margins (b,l,t,r)

+     bty="l") # box type (surrounding each plot region) is l (L-shape)

> # Note that each margin is in number of text lines.

> # Colors for distinct profiles will be set next.

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

> cols # colors used for distinct profiles

[1] "darkmagenta" "darkgreen"   "darkblue" 

> # To get a list of available colors, use command

> # colors()

> # For more info, do help(colors), help(palette), etc.

> # R function interaction.plot is used

> # to create interaction plot. Specify profile-factor,

> # trace-factor, and response in this order.

> with(data=factory,

+      {interaction.plot(shift,gender,output,type='b',

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

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

+       interaction.plot(gender,shift,output,type='b',

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

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

+      }

+     )

> # Note that type='b' above means both (points and lines).

> title(sub=" Factory Example", # place substitle

+       line=0.5, # start at line 0.5 outward

+       outer=TRUE, # at outer margin (bottom)

+       adj=0) # left justified

> title(main="Interaction Plot", # place main title

+       outer=TRUE) # at outer margin (top)

Make the graphics windows current by clicking inside it. Then on R console, do FileSave asPng (or Pdf, Bmp, etc). Make sure you include the extension in the file name (factory01.png in my example).

(factory01.png here)

> dev.off() # Turn off (windows) device

null device 
          1 

The session below generates main effect plots.

> print( ybar.A <− with(factory,tapply(output,shift,mean)) )

Shift 1 Shift 2 Shift 3 
      9       5       6    ⇐ means for shifts

> print( ybar.B <− with(factory,tapply(output,gender,mean)) )

    Male   Female 
5.000000 8.333333    ⇐ means for gender

> print( ybar <− mean(factory$output) )

[1] 6.666667    ⇐ grand mean

> ylim <− range(pretty(range(ybar.A,ybar.B,ybar)))

> 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(1:3, ybar.A,

+      xlab="shift", ylab="mean output",

+      ylim=ylim, # use the pre-calculated y ranges

+      xaxt="n", # no x-axis line

+      type='l', # plot type is line

+      lwd=2) # line width is 2

> axis(1, # side 1 axis (viz. x-axis)

+      at=1:3, # tick marks's positions

+      labels=levels(factory$shift)) # tick labels

> abline(h=ybar, # horizontal reference line

+        lty="dotted", # line type

+        col='darkorange4') # of this color

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

> plot(1:2, ybar.B, ylim=ylim,

+      xlim=c(.7,2.3), # use this ranges for x

+      xlab="gender", ylab="mean output",

+      xaxt='n', type='l', lwd=2)

> axis(1, at=1:2, lab=levels(factory$gender))

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

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

> title(outer=T,

+       main="Main Effect Plots\nfactory example")

> title(line=0.5,outer=TRUE,

+       sub=bquote(bar(y)[...]==.(round(ybar,2)))

+      )

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

> dev.off()

null device 
          1 

(factory02.png here)

As expected, the main effect plots do reflect the test results (that both main effects are significant). The follow-up study can be focused on the main effects.

> #

> x<−glht(factory.aov,linfct=mcp(shift="Tukey"))

> summary(x)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = output ~ shift * gender, data = factory)

Linear Hypotheses:
                       Estimate Std. Error t value p value  
Shift 2 - Shift 1 == 0       -4          1      -4  0.0171 *
Shift 3 - Shift 1 == 0       -3          1      -3  0.0544 .
Shift 3 - Shift 2 == 0        1          1       1  0.6035  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
(Adjusted p values reported)

> confint(x)

Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = output ~ shift * gender, data = factory)

Estimated Quantile = 3.0707  ≅ q.05;3,6/√2

Linear Hypotheses:
                       Estimate lwr      upr     
Shift 2 - Shift 1 == 0 -4.00000 -7.07068 -0.92932
Shift 3 - Shift 1 == 0 -3.00000 -6.07068  0.07068
Shift 3 - Shift 2 == 0  1.00000 -2.07068  4.07068

95% family-wise confidence level

> qtukey(.05,3,6,lower=F) # q.05;3,6

[1] 4.339195

> qtukey(.05,3,6,lower=F)/sqrt(2) # Tukey's critical value

[1] 3.068274

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

> par(mar=c(5,7,4,1)+.1) # leave 7 lines for left margin

> plot(x)

> title(main="factor shift",line=1)

> title(sub="Factory Example", line=3, adj=1)

> box(which="figure", col="cyan", lwd=2)

> dev.off()

null device 
          1

(factory03.png here)

> x <− glht(factory.aov,linfct=mcp(gender="Tukey"))

> summary(x)

         Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = output ~ shift * gender, data = factory)

Linear Hypotheses:
                   Estimate Std. Error t value p value   
Female - Male == 0   3.3333     0.8165   4.082 0.00648 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
(Adjusted p values reported)

> confint(x)

Simultaneous Confidence Intervals for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = output ~ shift * gender, data = factory)

Estimated Quantile = 2.4469

Linear Hypotheses:
                   Estimate lwr    upr   
Female - Male == 0 3.3333   1.3354 5.3312

95% family-wise confidence level

Since there is only one interval, the above Estimated Quantile is actually t.025,6.

> qt(.025,6,lower=F) # t.025,6

[1] 2.446912