# 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()

> # 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)

> 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).

> 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
```

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
```

> 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)

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

> dev.off()

```null device
1
```

> 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
```

> 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
```