Loading [MathJax]/jax/output/HTML-CSS/fonts/TeX/fontdata.js

NMSA407 Linear Regression: Tutorial

Two-way Analysis of Variance

Data Howells

Introduction

Full (unbalanced) dataset

data(HowellsAll, package = "mffSM")
head(HowellsAll)

##   gender popul oca gol fgender fpopul fgen.pop fpop.gen
## 1      1     1 123 176       M   BERG   M:BERG   BERG:M
## 2      1     1 115 173       M   BERG   M:BERG   BERG:M
## 3      1     1 117 176       M   BERG   M:BERG   BERG:M
## 4      1     1 113 185       M   BERG   M:BERG   BERG:M
## 5      1     1 121 170       M   BERG   M:BERG   BERG:M
## 6      1     1 118 180       M   BERG   M:BERG   BERG:M

dim(HowellsAll)

## [1] 289   8

summary(HowellsAll)

##      gender           popul            oca             gol        fgender    fpopul   
##  Min.   :0.0000   Min.   :0.000   Min.   :102.0   Min.   :155.0   F:156   AUSTR : 71  
##  1st Qu.:0.0000   1st Qu.:1.000   1st Qu.:113.0   1st Qu.:172.0   M:133   BERG  :109  
##  Median :0.0000   Median :1.000   Median :116.0   Median :178.0           BURIAT:109  
##  Mean   :0.4602   Mean   :1.131   Mean   :115.8   Mean   :178.1                       
##  3rd Qu.:1.0000   3rd Qu.:2.000   3rd Qu.:119.0   3rd Qu.:184.0                       
##  Max.   :1.0000   Max.   :2.000   Max.   :131.0   Max.   :198.0                       
##      fgen.pop      fpop.gen 
##  F:AUSTR :49   AUSTR:F :49  
##  F:BERG  :53   AUSTR:M :22  
##  F:BURIAT:54   BERG:F  :53  
##  M:AUSTR :22   BERG:M  :56  
##  M:BERG  :56   BURIAT:F:54  
##  M:BURIAT:55   BURIAT:M:55

Balanced subset (was artificially created)

data(Howells, package = "mffSM")
head(Howells)

##   gender popul oca gol fgender fpopul fgen.pop fpop.gen
## 1      0     0 110 191       F  AUSTR  F:AUSTR  AUSTR:F
## 2      0     0 117 183       F  AUSTR  F:AUSTR  AUSTR:F
## 3      0     0 119 181       F  AUSTR  F:AUSTR  AUSTR:F
## 4      0     0 108 172       F  AUSTR  F:AUSTR  AUSTR:F
## 5      0     0 113 190       F  AUSTR  F:AUSTR  AUSTR:F
## 6      0     0 110 180       F  AUSTR  F:AUSTR  AUSTR:F

dim(Howells)

## [1] 240   8

summary(Howells)

##      gender        popul        oca             gol        fgender    fpopul       fgen.pop 
##  Min.   :0.0   Min.   :0   Min.   :102.0   Min.   :155.0   F:120   AUSTR :80   F:AUSTR :40  
##  1st Qu.:0.0   1st Qu.:0   1st Qu.:112.0   1st Qu.:173.0   M:120   BERG  :80   F:BERG  :40  
##  Median :0.5   Median :1   Median :116.0   Median :179.0           BURIAT:80   F:BURIAT:40  
##  Mean   :0.5   Mean   :1   Mean   :115.7   Mean   :179.3                       M:AUSTR :40  
##  3rd Qu.:1.0   3rd Qu.:2   3rd Qu.:119.0   3rd Qu.:186.0                       M:BERG  :40  
##  Max.   :1.0   Max.   :2   Max.   :131.0   Max.   :201.0                       M:BURIAT:40  
##      fpop.gen 
##  AUSTR:F :40  
##  AUSTR:M :40  
##  BERG:F  :40  
##  BERG:M  :40  
##  BURIAT:F:40  
##  BURIAT:M:40

Some colors (used in plots below)

genCOL <- rainbow_hcl(2, start = 10, end = 280)
genCOL2 <- c("red3", "darkblue")
names(genCOL) <- names(genCOL2) <- levels(HowellsAll[, "fgender"])

#popCOL <- c(rainbow_hcl(3, start = 5, end = 310), rainbow_hcl(3, start = 25, end = 330))
popCOL <- rainbow_hcl(3, start = 50, end = 250)
popCOL2 <- c("red", "darkgreen", "blue")
names(popCOL) <- names(popCOL2) <- levels(HowellsAll[, "fpopul"])

Frequency tables

Full dataset

with(HowellsAll, table(fpopul, fgender))

##         fgender
## fpopul    F  M
##   AUSTR  49 22
##   BERG   53 56
##   BURIAT 54 55

with(HowellsAll, table(fpop.gen))

## fpop.gen
##  AUSTR:F  AUSTR:M   BERG:F   BERG:M BURIAT:F BURIAT:M 
##       49       22       53       56       54       55

with(HowellsAll, table(fgen.pop))

## fgen.pop
##  F:AUSTR   F:BERG F:BURIAT  M:AUSTR   M:BERG M:BURIAT 
##       49       53       54       22       56       55

margin.table(with(HowellsAll, table(fpopul, fgender)), margin = 1)

## fpopul
##  AUSTR   BERG BURIAT 
##     71    109    109

margin.table(with(HowellsAll, table(fpopul, fgender)), margin = 2)

## fgender
##   F   M 
## 156 133

nrow(HowellsAll)

## [1] 289

Balanced subset

with(Howells, table(fpopul, fgender))

##         fgender
## fpopul    F  M
##   AUSTR  40 40
##   BERG   40 40
##   BURIAT 40 40

with(Howells, table(fpop.gen))

## fpop.gen
##  AUSTR:F  AUSTR:M   BERG:F   BERG:M BURIAT:F BURIAT:M 
##       40       40       40       40       40       40

with(Howells, table(fgen.pop))

## fgen.pop
##  F:AUSTR   F:BERG F:BURIAT  M:AUSTR   M:BERG M:BURIAT 
##       40       40       40       40       40       40

margin.table(with(Howells, table(fpopul, fgender)), margin = 1)

## fpopul
##  AUSTR   BERG BURIAT 
##     80     80     80

margin.table(with(Howells, table(fpopul, fgender)), margin = 2)

## fgender
##   F   M 
## 120 120

nrow(Howells)

## [1] 240

Group sample means

Full dataset

OCA

with(HowellsAll, tapply(oca, list(fpopul, fgender), mean))

##               F        M
## AUSTR  114.6531 113.9545
## BERG   116.9623 117.1607
## BURIAT 117.0370 113.8364

round(with(HowellsAll, tapply(oca, list(fpopul, fgender), mean)), 1)

##            F     M
## AUSTR  114.7 114.0
## BERG   117.0 117.2
## BURIAT 117.0 113.8

round(with(HowellsAll, tapply(oca, fpopul, mean)), 1)

##  AUSTR   BERG BURIAT 
##  114.4  117.1  115.4

round(with(HowellsAll, tapply(oca, fgender, mean)), 1)

##     F     M 
## 116.3 115.3

round(with(HowellsAll, mean(oca)), 1)

## [1] 115.8

GOL

with(HowellsAll, tapply(gol, list(fpopul, fgender), mean))

##               F        M
## AUSTR  181.1020 190.7727
## BERG   170.5283 180.3214
## BURIAT 171.8333 181.6364

round(with(HowellsAll, tapply(gol, list(fpopul, fgender), mean)), 1)

##            F     M
## AUSTR  181.1 190.8
## BERG   170.5 180.3
## BURIAT 171.8 181.6

round(with(HowellsAll, tapply(gol, fpopul, mean)), 1)

##  AUSTR   BERG BURIAT 
##  184.1  175.6  176.8

round(with(HowellsAll, tapply(gol, fgender, mean)), 1)

##     F     M 
## 174.3 182.6

round(with(HowellsAll, mean(gol)), 1)

## [1] 178.1

Balanced subset

OCA

with(Howells, tapply(oca, list(fpopul, fgender), mean))

##             F       M
## AUSTR  114.80 115.025
## BERG   116.85 116.675
## BURIAT 117.20 113.450

round(with(Howells, tapply(oca, list(fpopul, fgender), mean)), 1)

##            F     M
## AUSTR  114.8 115.0
## BERG   116.8 116.7
## BURIAT 117.2 113.5

round(with(Howells, tapply(oca, fpopul, mean)), 1)

##  AUSTR   BERG BURIAT 
##  114.9  116.8  115.3

round(with(Howells, tapply(oca, fgender, mean)), 1)

##     F     M 
## 116.3 115.0

round(with(Howells, mean(oca)), 1)

## [1] 115.7

GOL

with(Howells, tapply(gol, list(fpopul, fgender), mean))

##              F       M
## AUSTR  181.375 190.375
## BERG   170.450 180.300
## BURIAT 172.175 181.175

round(with(Howells, tapply(gol, list(fpopul, fgender), mean)), 1)

##            F     M
## AUSTR  181.4 190.4
## BERG   170.4 180.3
## BURIAT 172.2 181.2

round(with(Howells, tapply(gol, fpopul, mean)), 1)

##  AUSTR   BERG BURIAT 
##  185.9  175.4  176.7

round(with(Howells, tapply(gol, fgender, mean)), 1)

##     F     M 
## 174.7 183.9

round(with(Howells, mean(gol)), 1)

## [1] 179.3

Boxplots (for full dataset)

OCA

plot(oca ~ fpop.gen, data = HowellsAll, col = genCOL)

plot of chunk fig-ANOVA-01-01

plot(oca ~ fgen.pop, data = HowellsAll, col = popCOL)

plot of chunk fig-ANOVA-01-02

GOL

plot(gol ~ fpop.gen, data = HowellsAll, col = genCOL)

plot of chunk fig-ANOVA-01-03

plot(gol ~ fgen.pop, data = HowellsAll, col = popCOL)

plot of chunk fig-ANOVA-01-04

Spaghetti plots (for full dataset)

with(HowellsAll, interaction.plot(fgender, fpopul, oca, col = popCOL2, lwd = 2))

plot of chunk fig-ANOVA-01-05

with(HowellsAll, interaction.plot(fpopul,  fgender,oca, col = genCOL2, lwd = 2))

plot of chunk fig-ANOVA-01-06

with(HowellsAll, interaction.plot(fgender, fpopul,  gol, col = popCOL2, lwd = 2))

plot of chunk fig-ANOVA-01-07

with(HowellsAll, interaction.plot(fpopul,  fgender, gol, col = genCOL2, lwd = 2))

plot of chunk fig-ANOVA-01-08

F-tests and ANOVA tables with a full (unbalanced) dataset

Models for OCA

oca0All      <- lm(oca ~ 1, data = HowellsAll)
ocaXAll      <- lm(oca ~ fpopul, data = HowellsAll)
ocaZAll      <- lm(oca ~ fgender, data = HowellsAll)
ocaXplusZAll <- lm(oca ~ fpopul + fgender, data = HowellsAll)
ocaXZAll     <- lm(oca ~ fpopul * fgender, data = HowellsAll)
ocaXZAll2    <- lm(oca ~ fgender * fpopul, data = HowellsAll)

Type I ANOVA tables for models `ocaXZAll` and `ocaXZAll2`

anova(ocaXZAll)

## Analysis of Variance Table
## 
## Response: oca
##                 Df Sum Sq Mean Sq F value   Pr(>F)   
## fpopul           2  321.8 160.879  6.3579 0.001991 **
## fgender          1  122.6 122.597  4.8450 0.028534 * 
## fpopul:fgender   2  165.0  82.509  3.2607 0.039807 * 
## Residuals      283 7161.0  25.304                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(ocaXZAll2)

## Analysis of Variance Table
## 
## Response: oca
##                 Df Sum Sq Mean Sq F value    Pr(>F)    
## fgender          1   72.8  72.827  2.8781 0.0908910 .  
## fpopul           2  371.5 185.764  7.3413 0.0007792 ***
## fgender:fpopul   2  165.0  82.509  3.2607 0.0398071 *  
## Residuals      283 7161.0  25.304                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type II ANOVA table for model `ocaXZAll`

Function Anova comes from package car.

library("car")
Anova(ocaXZAll, type = "II")

## Anova Table (Type II tests)
## 
## Response: oca
##                Sum Sq  Df F value    Pr(>F)    
## fpopul          371.5   2  7.3413 0.0007792 ***
## fgender         122.6   1  4.8450 0.0285335 *  
## fpopul:fgender  165.0   2  3.2607 0.0398071 *  
## Residuals      7161.0 283                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type III ANOVA table for model `ocaXZAll`

Anova(ocaXZAll, type = "III")

## Anova Table (Type III tests)
## 
## Response: oca
##                Sum Sq  Df    F value  Pr(>F)    
## (Intercept)    644121   1 25455.4563 < 2e-16 ***
## fpopul            185   2     3.6609 0.02693 *  
## fgender             7   1     0.2928 0.58888    
## fpopul:fgender    165   2     3.2607 0.03981 *  
## Residuals        7161 283                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Rows fpopul, fgender are not really useful.

Models for GOL

gol0All      <- lm(gol ~ 1, data = HowellsAll)
golXAll      <- lm(gol ~ fpopul, data = HowellsAll)
golZAll      <- lm(gol ~ fgender, data = HowellsAll)
golXplusZAll <- lm(gol ~ fpopul + fgender, data = HowellsAll)
golXZAll     <- lm(gol ~ fpopul * fgender, data = HowellsAll)
golXZAll2    <- lm(gol ~ fgender * fpopul, data = HowellsAll)

Type I ANOVA tables for models `golXZAll` and `golXZAll2`

anova(golXZAll)

## Analysis of Variance Table
## 
## Response: gol
##                 Df  Sum Sq Mean Sq  F value Pr(>F)    
## fpopul           2  3448.1  1724.1  43.3542 <2e-16 ***
## fgender          1  6649.7  6649.7 167.2172 <2e-16 ***
## fpopul:fgender   2     0.2     0.1   0.0024 0.9976    
## Residuals      283 11254.0    39.8                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(golXZAll2)

## Analysis of Variance Table
## 
## Response: gol
##                 Df  Sum Sq Mean Sq  F value Pr(>F)    
## fgender          1  4937.1  4937.1 124.1509 <2e-16 ***
## fpopul           2  5160.7  2580.4  64.8873 <2e-16 ***
## fgender:fpopul   2     0.2     0.1   0.0024 0.9976    
## Residuals      283 11254.0    39.8                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type II ANOVA table for model `golXZAll`

Anova(golXZAll, type = "II")

## Anova Table (Type II tests)
## 
## Response: gol
##                 Sum Sq  Df  F value Pr(>F)    
## fpopul          5160.7   2  64.8873 <2e-16 ***
## fgender         6649.7   1 167.2172 <2e-16 ***
## fpopul:fgender     0.2   2   0.0024 0.9976    
## Residuals      11254.0 283                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type III ANOVA table for model `golXZAll`

Anova(golXZAll, type = "III")

## Anova Table (Type III tests)
## 
## Response: gol
##                 Sum Sq  Df    F value    Pr(>F)    
## (Intercept)    1607100   1 40413.1028 < 2.2e-16 ***
## fpopul            3350   2    42.1161 < 2.2e-16 ***
## fgender           1420   1    35.7071 6.879e-09 ***
## fpopul:fgender       0   2     0.0024    0.9976    
## Residuals        11254 283                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Rows fpopul, fgender are not really useful.

F-tests and ANOVA tables with a balanced dataset

Models for OCA

oca0      <- lm(oca ~ 1, data = Howells)
ocaX      <- lm(oca ~ fpopul, data = Howells)
ocaZ      <- lm(oca ~ fgender, data = Howells)
ocaXplusZ <- lm(oca ~ fpopul + fgender, data = Howells)
ocaXZ     <- lm(oca ~ fpopul * fgender, data = Howells)
ocaXZ2    <- lm(oca ~ fgender * fpopul, data = Howells)

Type I ANOVA tables for models `ocaXZ` and `ocaXZ2`

anova(ocaXZ)

## Analysis of Variance Table
## 
## Response: oca
##                 Df Sum Sq Mean Sq F value  Pr(>F)  
## fpopul           2  150.9  75.454  3.0497 0.04926 *
## fgender          1   91.3  91.267  3.6888 0.05599 .
## fpopul:fgender   2  191.6  95.804  3.8722 0.02216 *
## Residuals      234 5789.6  24.742                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(ocaXZ2)

## Analysis of Variance Table
## 
## Response: oca
##                 Df Sum Sq Mean Sq F value  Pr(>F)  
## fgender          1   91.3  91.267  3.6888 0.05599 .
## fpopul           2  150.9  75.454  3.0497 0.04926 *
## fgender:fpopul   2  191.6  95.804  3.8722 0.02216 *
## Residuals      234 5789.6  24.742                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

With balanced dataset, type I ANOVA table in a model with two categorical covariates that includes their interaction does not depend on the ordering of the rows with the main effects. This can be generally proved.

Type II ANOVA table for model `ocaXZ`

Anova(ocaXZ, type = "II")

## Anova Table (Type II tests)
## 
## Response: oca
##                Sum Sq  Df F value  Pr(>F)  
## fpopul          150.9   2  3.0497 0.04926 *
## fgender          91.3   1  3.6888 0.05599 .
## fpopul:fgender  191.6   2  3.8722 0.02216 *
## Residuals      5789.6 234                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type III ANOVA table for model `ocaXZ`

Anova(ocaXZ, type = "III")

## Anova Table (Type III tests)
## 
## Response: oca
##                Sum Sq  Df    F value  Pr(>F)    
## (Intercept)    527162   1 21306.6325 < 2e-16 ***
## fpopul            134   2     2.7174 0.06813 .  
## fgender             1   1     0.0409 0.83986    
## fpopul:fgender    192   2     3.8722 0.02216 *  
## Residuals        5790 234                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Rows fpopul, fgender are not really useful.

Models for GOL

gol0      <- lm(gol ~ 1, data = Howells)
golX      <- lm(gol ~ fpopul, data = Howells)
golZ      <- lm(gol ~ fgender, data = Howells)
golXplusZ <- lm(gol ~ fpopul + fgender, data = Howells)
golXZ     <- lm(gol ~ fpopul * fgender, data = Howells)
golXZ2    <- lm(gol ~ fgender * fpopul, data = Howells)

Type I ANOVA tables for models `golXZ` and `golXZ`

anova(golXZ)

## Analysis of Variance Table
## 
## Response: gol
##                 Df Sum Sq Mean Sq  F value Pr(>F)    
## fpopul           2 5242.1  2621.1  65.1743 <2e-16 ***
## fgender          1 5170.8  5170.8 128.5753 <2e-16 ***
## fpopul:fgender   2    9.6     4.8   0.1198 0.8872    
## Residuals      234 9410.6    40.2                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(golXZ2)

## Analysis of Variance Table
## 
## Response: gol
##                 Df Sum Sq Mean Sq  F value Pr(>F)    
## fgender          1 5170.8  5170.8 128.5753 <2e-16 ***
## fpopul           2 5242.1  2621.1  65.1743 <2e-16 ***
## fgender:fpopul   2    9.6     4.8   0.1198 0.8872    
## Residuals      234 9410.6    40.2                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type II ANOVA table for model `golXZ`

Anova(golXZ, type = "II")

## Anova Table (Type II tests)
## 
## Response: gol
##                Sum Sq  Df  F value Pr(>F)    
## fpopul         5242.1   2  65.1743 <2e-16 ***
## fgender        5170.8   1 128.5753 <2e-16 ***
## fpopul:fgender    9.6   2   0.1198 0.8872    
## Residuals      9410.6 234                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Type III ANOVA table for model `golXZ`

Anova(golXZ, type = "III")

## Anova Table (Type III tests)
## 
## Response: gol
##                 Sum Sq  Df    F value    Pr(>F)    
## (Intercept)    1315876   1 32720.0068 < 2.2e-16 ***
## fpopul            2760   2    34.3097 8.577e-14 ***
## fgender           1620   1    40.2822 1.128e-09 ***
## fpopul:fgender      10   2     0.1198    0.8872    
## Residuals         9411 234                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Rows fpopul, fgender are not really useful.

Levels of the two factors

(lpop <- levels(HowellsAll[, "fpopul"]))

## [1] "AUSTR"  "BERG"   "BURIAT"

(lgen <- levels(HowellsAll[, "fgender"]))

## [1] "F" "M"

outer(lpop, lgen, paste, sep = ":")

##      [,1]       [,2]      
## [1,] "AUSTR:F"  "AUSTR:M" 
## [2,] "BERG:F"   "BERG:M"  
## [3,] "BURIAT:F" "BURIAT:M"

(lpopgen <- as.character(outer(lpop, lgen, paste, sep = ":")))

## [1] "AUSTR:F"  "BERG:F"   "BURIAT:F" "AUSTR:M"  "BERG:M"   "BURIAT:M"

Reference group pseudocontrast parameterization `contr.treatment`

Pseudocontrast matrices

(CX.trt <- contr.treatment(3))

##   2 3
## 1 0 0
## 2 1 0
## 3 0 1

(CZ.trt <- contr.treatment(2))

##   2
## 1 0
## 2 1

(CXZ.trt <- kronecker(CZ.trt, CX.trt))

##      [,1] [,2]
## [1,]    0    0
## [2,]    0    0
## [3,]    0    0
## [4,]    0    0
## [5,]    1    0
## [6,]    0    1

rownames(CX.trt) <- lpop
rownames(CZ.trt) <- lgen
rownames(CXZ.trt) <- lpopgen

print(CX.trt)

##        2 3
## AUSTR  0 0
## BERG   1 0
## BURIAT 0 1

print(CZ.trt)

##   2
## F 0
## M 1

print(CXZ.trt)

##          [,1] [,2]
## AUSTR:F     0    0
## BERG:F      0    0
## BURIAT:F    0    0
## AUSTR:M     0    0
## BERG:M      1    0
## BURIAT:M    0    1

Matrix that (together with the intercept term) parameterizes the two-way classified group means

(D.trt <- cbind(kronecker(rep(1, 2), CX.trt), kronecker(CZ.trt, rep(1, 3)), CXZ.trt))

##          [,1] [,2] [,3] [,4] [,5]
## AUSTR:F     0    0    0    0    0
## BERG:F      1    0    0    0    0
## BURIAT:F    0    1    0    0    0
## AUSTR:M     0    0    1    0    0
## BERG:M      1    0    1    1    0
## BURIAT:M    0    1    1    0    1

OCA (full dataset)

Interaction model

ocaXZAll.trt <- lm(oca ~ fpopul * fgender, data = HowellsAll)
summary(ocaXZAll.trt)

## 
## Call:
## lm(formula = oca ~ fpopul * fgender, data = HowellsAll)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.1607  -3.1607   0.0455   3.1636  13.8393 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           114.6531     0.7186 159.548   <2e-16 ***
## fpopulBERG              2.3092     0.9969   2.316   0.0213 *  
## fpopulBURIAT            2.3840     0.9925   2.402   0.0169 *  
## fgenderM               -0.6985     1.2910  -0.541   0.5889    
## fpopulBERG:fgenderM     0.8970     1.6112   0.557   0.5782    
## fpopulBURIAT:fgenderM  -2.5022     1.6110  -1.553   0.1215    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.03 on 283 degrees of freedom
## Multiple R-squared:  0.07842,    Adjusted R-squared:  0.06214 
## F-statistic: 4.816 on 5 and 283 DF,  p-value: 0.0003046

drop1(ocaXZAll.trt, test = "F")

## Single term deletions
## 
## Model:
## oca ~ fpopul * fgender
##                Df Sum of Sq  RSS    AIC F value  Pr(>F)  
## <none>                      7161 939.68                  
## fpopul:fgender  2    165.02 7326 942.27  3.2607 0.03981 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

GOL (full dataset)

Interaction model

golXZAll.trt <- lm(gol ~ fpopul * fgender, data = HowellsAll)
summary(golXZAll.trt)

## 
## Call:
## lm(formula = gol ~ fpopul * fgender, data = HowellsAll)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5283  -4.3214  -0.3214   4.4717  17.6786 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           181.1020     0.9009 201.030  < 2e-16 ***
## fpopulBERG            -10.5737     1.2498  -8.461 1.45e-15 ***
## fpopulBURIAT           -9.2687     1.2442  -7.450 1.14e-12 ***
## fgenderM                9.6707     1.6184   5.976 6.88e-09 ***
## fpopulBERG:fgenderM     0.1224     2.0198   0.061    0.952    
## fpopulBURIAT:fgenderM   0.1323     2.0196   0.066    0.948    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.306 on 283 degrees of freedom
## Multiple R-squared:  0.4729, Adjusted R-squared:  0.4636 
## F-statistic: 50.79 on 5 and 283 DF,  p-value: < 2.2e-16

drop1(golXZAll.trt, test = "F")

## Single term deletions
## 
## Model:
## gol ~ fpopul * fgender
##                Df Sum of Sq   RSS    AIC F value Pr(>F)
## <none>                      11254 1070.3               
## fpopul:fgender  2   0.19404 11254 1066.3  0.0024 0.9976

Additive model

golXplusZAll.trt <- lm(gol ~ fpopul + fgender, data = HowellsAll)
summary(golXplusZAll.trt)

## 
## Call:
## lm(formula = gol ~ fpopul + fgender, data = HowellsAll)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5400  -4.3103  -0.3103   4.4600  17.6897 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  181.0712     0.7814 231.724   <2e-16 ***
## fpopulBERG   -10.5311     0.9706 -10.850   <2e-16 ***
## fpopulBURIAT  -9.2213     0.9695  -9.511   <2e-16 ***
## fgenderM       9.7703     0.7529  12.977   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.284 on 285 degrees of freedom
## Multiple R-squared:  0.4729, Adjusted R-squared:  0.4674 
## F-statistic: 85.24 on 3 and 285 DF,  p-value: < 2.2e-16

drop1(golXplusZAll.trt, test = "F")

## Single term deletions
## 
## Model:
## gol ~ fpopul + fgender
##         Df Sum of Sq   RSS    AIC F value    Pr(>F)    
## <none>               11254 1066.3                      
## fpopul   2    5160.7 16415 1171.4  65.345 < 2.2e-16 ***
## fgender  1    6649.7 17904 1198.5 168.396 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimates of all (partial) effects of `fpopul` given `fgender` in the additive model

Function LSest comes from package mffSM.

library("mffSM")
L.effX.trt <- cbind(0, rbind(CX.trt[2,] - CX.trt[1,], CX.trt[3,] - CX.trt[1,], CX.trt[3,] - CX.trt[2,]), 0)
colnames(L.effX.trt) <- names(coef(golXplusZAll.trt))
rownames(L.effX.trt) <- paste(lpop[c(2, 3, 3)], "-", lpop[c(1, 1, 2)], sep = "")
print(L.effX.trt)

##              (Intercept) fpopulBERG fpopulBURIAT fgenderM
## BERG-AUSTR             0          1            0        0
## BURIAT-AUSTR           0          0            1        0
## BURIAT-BERG            0         -1            1        0

LSest(golXplusZAll.trt, L = L.effX.trt)

##                Estimate Std. Error    t value P value       Lower     Upper
## BERG-AUSTR   -10.531148  0.9705782 -10.850385 < 2e-16 -12.4415591 -8.620737
## BURIAT-AUSTR  -9.221329  0.9695097  -9.511332 < 2e-16 -11.1296364 -7.313021
## BURIAT-BERG    1.309819  0.8512377   1.538723 0.12498  -0.3656911  2.985330

Estimates of all (partial) effects of `fgender` given `fpopul` in the additive model

L.effZ.trt <- cbind(0, 0, 0, CZ.trt[2,] - CZ.trt[1,])
colnames(L.effZ.trt) <- names(coef(golXplusZAll.trt))
rownames(L.effZ.trt) <- paste(lgen[2], "-", lgen[1], sep = "")
print(L.effZ.trt)

##     (Intercept) fpopulBERG fpopulBURIAT fgenderM
## M-F           0          0            0        1

LSest(golXplusZAll.trt, L = L.effZ.trt)

##     Estimate Std. Error  t value    P value    Lower    Upper
## M-F 9.770313  0.7529092 12.97675 < 2.22e-16 8.288345 11.25228

Sum contrast parameterization `contr.sum`

Contrast matrices

(CX.sum <- contr.sum(3))

##   [,1] [,2]
## 1    1    0
## 2    0    1
## 3   -1   -1

(CZ.sum <- contr.sum(2))

##   [,1]
## 1    1
## 2   -1

(CXZ.sum <- kronecker(CZ.sum, CX.sum))

##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
## [3,]   -1   -1
## [4,]   -1    0
## [5,]    0   -1
## [6,]    1    1

rownames(CX.sum) <- lpop
rownames(CZ.sum) <- lgen
rownames(CXZ.sum) <- lpopgen

print(CX.sum)

##        [,1] [,2]
## AUSTR     1    0
## BERG      0    1
## BURIAT   -1   -1

print(CZ.sum)

##   [,1]
## F    1
## M   -1

print(CXZ.sum)

##          [,1] [,2]
## AUSTR:F     1    0
## BERG:F      0    1
## BURIAT:F   -1   -1
## AUSTR:M    -1    0
## BERG:M      0   -1
## BURIAT:M    1    1

Matrix that (together with the intercept term) parameterizes the two-way classified group means

(D.sum <- cbind(kronecker(rep(1, 2), CX.sum), kronecker(CZ.sum, rep(1, 3)), CXZ.sum))

##          [,1] [,2] [,3] [,4] [,5]
## AUSTR:F     1    0    1    1    0
## BERG:F      0    1    1    0    1
## BURIAT:F   -1   -1    1   -1   -1
## AUSTR:M     1    0   -1   -1    0
## BERG:M      0    1   -1    0   -1
## BURIAT:M   -1   -1   -1    1    1

OCA (full dataset)

Interaction model

ocaXZAll.sum <- lm(oca ~ fpopul * fgender, data = HowellsAll,
                   contrasts = list(fpopul = "contr.sum", fgender = "contr.sum"))
summary(ocaXZAll.sum)

## 
## Call:
## lm(formula = oca ~ fpopul * fgender, data = HowellsAll, contrasts = list(fpopul = "contr.sum", 
##     fgender = "contr.sum"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.1607  -3.1607   0.0455   3.1636  13.8393 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      115.6007     0.3129 369.455  < 2e-16 ***
## fpopul1           -1.2969     0.4866  -2.665 0.008138 ** 
## fpopul2            1.4608     0.4187   3.489 0.000563 ***
## fgender1           0.6168     0.3129   1.971 0.049671 *  
## fpopul1:fgender1  -0.2675     0.4866  -0.550 0.582896    
## fpopul2:fgender1  -0.7160     0.4187  -1.710 0.088376 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.03 on 283 degrees of freedom
## Multiple R-squared:  0.07842,    Adjusted R-squared:  0.06214 
## F-statistic: 4.816 on 5 and 283 DF,  p-value: 0.0003046

drop1(ocaXZAll.sum, test = "F")

## Single term deletions
## 
## Model:
## oca ~ fpopul * fgender
##                Df Sum of Sq  RSS    AIC F value  Pr(>F)  
## <none>                      7161 939.68                  
## fpopul:fgender  2    165.02 7326 942.27  3.2607 0.03981 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimates of all $\alpha^X$ parameters

L.alphaX <- cbind(0, CX.sum, 0, 0, 0)
colnames(L.alphaX) <- names(coef(ocaXZAll.sum))
rownames(L.alphaX) <- lpop
print(L.alphaX)

##        (Intercept) fpopul1 fpopul2 fgender1 fpopul1:fgender1 fpopul2:fgender1
## AUSTR            0       1       0        0                0                0
## BERG             0       0       1        0                0                0
## BURIAT           0      -1      -1        0                0                0

LSest(ocaXZAll.sum, L = L.alphaX)

##         Estimate Std. Error    t value    P value      Lower      Upper
## AUSTR  -1.296861  0.4866057 -2.6651167 0.00813815 -2.2546868 -0.3390351
## BERG    1.460825  0.4187409  3.4886132 0.00056284  0.6365830  2.2850669
## BURIAT -0.163964  0.4186786 -0.3916225 0.69563187 -0.9880833  0.6601554

alphaX <- LSest(ocaXZAll.sum, L = L.alphaX)[["Estimate"]]
names(alphaX) <- lpop
print(alphaX)

##     AUSTR      BERG    BURIAT 
## -1.296861  1.460825 -0.163964

Estimates of all $\alpha^Z$ parameters

L.alphaZ <- cbind(0, 0, 0, CZ.sum, 0, 0)
colnames(L.alphaZ) <- names(coef(ocaXZAll.sum))
rownames(L.alphaZ) <- lgen
print(L.alphaZ)

##   (Intercept) fpopul1 fpopul2 fgender1 fpopul1:fgender1 fpopul2:fgender1
## F           0       0       0        1                0                0
## M           0       0       0       -1                0                0

LSest(ocaXZAll.sum, L = L.alphaZ)

##     Estimate Std. Error   t value  P value         Lower         Upper
## F  0.6167898  0.3128953  1.971234 0.049671  0.0008924028  1.2326872759
## M -0.6167898  0.3128953 -1.971234 0.049671 -1.2326872759 -0.0008924028

alphaZ <- LSest(ocaXZAll.sum, L = L.alphaZ)[["Estimate"]]
names(alphaZ) <- lgen
print(alphaZ)

##          F          M 
##  0.6167898 -0.6167898

Estimates of all $\alpha^{XZ}$ parameters

L.alphaXZ <- cbind(0, 0, 0, 0, CXZ.sum)
colnames(L.alphaXZ) <- names(coef(ocaXZAll.sum))
rownames(L.alphaXZ) <- lpopgen
print(L.alphaXZ)

##          (Intercept) fpopul1 fpopul2 fgender1 fpopul1:fgender1 fpopul2:fgender1
## AUSTR:F            0       0       0        0                1                0
## BERG:F             0       0       0        0                0                1
## BURIAT:F           0       0       0        0               -1               -1
## AUSTR:M            0       0       0        0               -1                0
## BERG:M             0       0       0        0                0               -1
## BURIAT:M           0       0       0        0                1                1

LSest(ocaXZAll.sum, L = L.alphaXZ)

##            Estimate Std. Error    t value  P value      Lower      Upper
## AUSTR:F  -0.2675320  0.4866057 -0.5497921 0.582896 -1.2253578  0.6902939
## BERG:F   -0.7160149  0.4187409 -1.7099236 0.088376 -1.5402569  0.1082270
## BURIAT:F  0.9835469  0.4186786  2.3491692 0.019503  0.1594275  1.8076662
## AUSTR:M   0.2675320  0.4866057  0.5497921 0.582896 -0.6902939  1.2253578
## BERG:M    0.7160149  0.4187409  1.7099236 0.088376 -0.1082270  1.5402569
## BURIAT:M -0.9835469  0.4186786 -2.3491692 0.019503 -1.8076662 -0.1594275

alphaXZ <- LSest(ocaXZAll.sum, L = L.alphaXZ)[["Estimate"]]
names(alphaXZ) <- lpopgen
print(alphaXZ)

##    AUSTR:F     BERG:F   BURIAT:F    AUSTR:M     BERG:M   BURIAT:M 
## -0.2675320 -0.7160149  0.9835469  0.2675320  0.7160149 -0.9835469

Estimates of $\alpha^{XZ}$ parameters stored in a matrix

(alphaXZm <- matrix(alphaXZ, nrow = length(lpop), ncol = length(lgen)))

##            [,1]       [,2]
## [1,] -0.2675320  0.2675320
## [2,] -0.7160149  0.7160149
## [3,]  0.9835469 -0.9835469

rownames(alphaXZm) <- lpop
colnames(alphaXZm) <- lgen
print(alphaXZm)

##                 F          M
## AUSTR  -0.2675320  0.2675320
## BERG   -0.7160149  0.7160149
## BURIAT  0.9835469 -0.9835469

Estimated group means

coef(ocaXZAll.sum)["(Intercept)"] + alphaXZm + cbind(alphaX, alphaX) + rbind(alphaZ, alphaZ, alphaZ)

##               F        M
## AUSTR  114.6531 113.9545
## BERG   116.9623 117.1607
## BURIAT 117.0370 113.8364

(ocaMean <- with(HowellsAll, tapply(oca, list(fpopul, fgender), mean)))

##               F        M
## AUSTR  114.6531 113.9545
## BERG   116.9623 117.1607
## BURIAT 117.0370 113.8364

Estimates of the differences between the means of the means of groups by `fpopul`

Those are equal to differences of $\alpha^X$ parameters.

L.effX.sumXZ <- cbind(0, rbind(CX.sum[2,] - CX.sum[1,], CX.sum[3,] - CX.sum[1,], CX.sum[3,] - CX.sum[2,]), 0, 0, 0)
colnames(L.effX.sumXZ) <- names(coef(ocaXZAll.sum))
rownames(L.effX.sumXZ) <- paste(lpop[c(2, 3, 3)], "-", lpop[c(1, 1, 2)], sep = "")
print(L.effX.sumXZ)

##              (Intercept) fpopul1 fpopul2 fgender1 fpopul1:fgender1 fpopul2:fgender1
## BERG-AUSTR             0      -1       1        0                0                0
## BURIAT-AUSTR           0      -2      -1        0                0                0
## BURIAT-BERG            0      -1      -2        0                0                0

LSest(ocaXZAll.sum, L = L.effX.sumXZ)

##               Estimate Std. Error   t value    P value      Lower      Upper
## BERG-AUSTR    2.757686  0.8055844  3.423211 0.00071032  1.1719880  4.3433837
## BURIAT-AUSTR  1.132897  0.8054873  1.406474 0.16068053 -0.4526097  2.7184037
## BURIAT-BERG  -1.624789  0.6815323 -2.384023 0.01778371 -2.9663047 -0.2832731

Estimated means of the means and their differences

apply(ocaMean, 1, mean)

##    AUSTR     BERG   BURIAT 
## 114.3038 117.0615 115.4367

outer(apply(ocaMean, 1, mean), apply(ocaMean, 1, mean), "-")

##           AUSTR      BERG    BURIAT
## AUSTR  0.000000 -2.757686 -1.132897
## BERG   2.757686  0.000000  1.624789
## BURIAT 1.132897 -1.624789  0.000000

Estimates of the differences between the means of the means of groups by `fgender`

Those are equal to differences of $\alpha^Z$ parameters.

L.effZ.sumXZ <- cbind(0, 0, 0, CZ.sum[2,] - CZ.sum[1,], 0, 0)
colnames(L.effZ.sumXZ) <- names(coef(ocaXZAll.sum))
rownames(L.effZ.sumXZ) <- paste(lgen[2], "-", lgen[1], sep = "")
print(L.effZ.sumXZ)

##     (Intercept) fpopul1 fpopul2 fgender1 fpopul1:fgender1 fpopul2:fgender1
## M-F           0       0       0       -2                0                0

LSest(ocaXZAll.sum, L = L.effZ.sumXZ)

##     Estimate Std. Error   t value  P value     Lower        Upper
## M-F -1.23358  0.6257906 -1.971234 0.049671 -2.465375 -0.001784806

Estimated means of the means and their differences

apply(ocaMean, 2, mean)

##        F        M 
## 116.2175 114.9839

outer(apply(ocaMean, 2, mean), apply(ocaMean, 2, mean), "-")

##          F       M
## F  0.00000 1.23358
## M -1.23358 0.00000

GOL (full dataset)

Interaction model

golXZAll.sum <- lm(gol ~ fpopul * fgender, data = HowellsAll,
                   contrasts = list(fpopul = "contr.sum", fgender = "contr.sum"))
summary(golXZAll.sum)

## 
## Call:
## lm(formula = gol ~ fpopul * fgender, data = HowellsAll, contrasts = list(fpopul = "contr.sum", 
##     fgender = "contr.sum"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5283  -4.3214  -0.3214   4.4717  17.6786 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      179.36570    0.39225 457.271  < 2e-16 ***
## fpopul1            6.57168    0.61002  10.773  < 2e-16 ***
## fpopul2           -3.94083    0.52494  -7.507 7.91e-13 ***
## fgender1          -4.87781    0.39225 -12.435  < 2e-16 ***
## fpopul1:fgender1   0.04246    0.61002   0.070    0.945    
## fpopul2:fgender1  -0.01876    0.52494  -0.036    0.972    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.306 on 283 degrees of freedom
## Multiple R-squared:  0.4729, Adjusted R-squared:  0.4636 
## F-statistic: 50.79 on 5 and 283 DF,  p-value: < 2.2e-16

drop1(golXZAll.sum, test = "F")

## Single term deletions
## 
## Model:
## gol ~ fpopul * fgender
##                Df Sum of Sq   RSS    AIC F value Pr(>F)
## <none>                      11254 1070.3               
## fpopul:fgender  2   0.19404 11254 1066.3  0.0024 0.9976

Additive model

golXplusZAll.sum <- lm(gol ~ fpopul + fgender, data = HowellsAll,
                   contrasts = list(fpopul = "contr.sum", fgender = "contr.sum"))                       
summary(golXplusZAll.sum)

## 
## Call:
## lm(formula = gol ~ fpopul + fgender, data = HowellsAll, contrasts = list(fpopul = "contr.sum", 
##     fgender = "contr.sum"))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5400  -4.3103  -0.3103   4.4600  17.6897 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 179.3722     0.3797 472.421  < 2e-16 ***
## fpopul1       6.5842     0.5811  11.330  < 2e-16 ***
## fpopul2      -3.9470     0.5157  -7.654 3.03e-13 ***
## fgender1     -4.8852     0.3765 -12.977  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.284 on 285 degrees of freedom
## Multiple R-squared:  0.4729, Adjusted R-squared:  0.4674 
## F-statistic: 85.24 on 3 and 285 DF,  p-value: < 2.2e-16

drop1(golXplusZAll.sum, test = "F")

## Single term deletions
## 
## Model:
## gol ~ fpopul + fgender
##         Df Sum of Sq   RSS    AIC F value    Pr(>F)    
## <none>               11254 1066.3                      
## fpopul   2    5160.7 16415 1171.4  65.345 < 2.2e-16 ***
## fgender  1    6649.7 17904 1198.5 168.396 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimates of all $\alpha^X$ parameters in the additive model

L.alphaX <- cbind(0, CX.sum, 0)
colnames(L.alphaX) <- names(coef(golXplusZAll.sum))
rownames(L.alphaX) <- lpop
print(L.alphaX)

##        (Intercept) fpopul1 fpopul2 fgender1
## AUSTR            0       1       0        0
## BERG             0       0       1        0
## BURIAT           0      -1      -1        0

LSest(golXplusZAll.sum, L = L.alphaX)

##         Estimate Std. Error   t value    P value     Lower     Upper
## AUSTR   6.584159  0.5811231 11.330059 < 2.22e-16  5.440321  7.727997
## BERG   -3.946989  0.5156772 -7.653992 3.0336e-13 -4.962008 -2.931970
## BURIAT -2.637170  0.5150067 -5.120651 5.6141e-07 -3.650869 -1.623470

alphaX <- LSest(golXplusZAll.sum, L = L.alphaX)[["Estimate"]]
names(alphaX) <- lpop
print(alphaX)

##     AUSTR      BERG    BURIAT 
##  6.584159 -3.946989 -2.637170

Estimates of all $\alpha^Z$ parameters in the additive model

L.alphaZ <- cbind(0, 0, 0, CZ.sum)
colnames(L.alphaZ) <- names(coef(golXplusZAll.sum))
rownames(L.alphaZ) <- lgen
print(L.alphaZ)

##   (Intercept) fpopul1 fpopul2 fgender1
## F           0       0       0        1
## M           0       0       0       -1

LSest(golXplusZAll.sum, L = L.alphaZ)

##    Estimate Std. Error   t value    P value     Lower     Upper
## F -4.885157  0.3764546 -12.97675 < 2.22e-16 -5.626141 -4.144173
## M  4.885157  0.3764546  12.97675 < 2.22e-16  4.144173  5.626141

alphaZ <- LSest(golXplusZAll.sum, L = L.alphaZ)[["Estimate"]]
names(alphaZ) <- lgen
print(alphaZ)

##         F         M 
## -4.885157  4.885157

Estimated group means: model based

(golModelMean <- coef(golXplusZAll.sum)["(Intercept)"] + cbind(alphaX, alphaX) + rbind(alphaZ, alphaZ, alphaZ))

##          alphaX   alphaX
## AUSTR  181.0712 190.8415
## BERG   170.5400 180.3103
## BURIAT 171.8498 181.6202

colnames(golModelMean) <- names(alphaZ)
print(golModelMean)

##               F        M
## AUSTR  181.0712 190.8415
## BERG   170.5400 180.3103
## BURIAT 171.8498 181.6202

Estimated group means: sample means

(golMean <- with(HowellsAll, tapply(gol, list(fpopul, fgender), mean)))

##               F        M
## AUSTR  181.1020 190.7727
## BERG   170.5283 180.3214
## BURIAT 171.8333 181.6364

Estimates of the differences between the means of the means of groups by `fpopul`

Those are partial effects of fpopul given fgender (if the additive model can be assumed). They are equal to differences in $\alpha^X$ 's.

L.peffX.sumXZ <- cbind(0, rbind(CX.sum[2,] - CX.sum[1,], CX.sum[3,] - CX.sum[1,], CX.sum[3,] - CX.sum[2,]), 0)
colnames(L.peffX.sumXZ) <- names(coef(golXplusZAll.sum))
rownames(L.peffX.sumXZ) <- paste(lpop[c(2, 3, 3)], "-", lpop[c(1, 1, 2)], sep = "")
print(L.peffX.sumXZ)

##              (Intercept) fpopul1 fpopul2 fgender1
## BERG-AUSTR             0      -1       1        0
## BURIAT-AUSTR           0      -2      -1        0
## BURIAT-BERG            0      -1      -2        0

LSest(golXplusZAll.sum, L = L.peffX.sumXZ)

##                Estimate Std. Error    t value P value       Lower     Upper
## BERG-AUSTR   -10.531148  0.9705782 -10.850385 < 2e-16 -12.4415591 -8.620737
## BURIAT-AUSTR  -9.221329  0.9695097  -9.511332 < 2e-16 -11.1296364 -7.313021
## BURIAT-BERG    1.309819  0.8512377   1.538723 0.12498  -0.3656911  2.985330

Estimated means of the means (model-based) and their differences

apply(golModelMean, 1, mean)

##    AUSTR     BERG   BURIAT 
## 185.9563 175.4252 176.7350

outer(apply(golModelMean, 1, mean), apply(golModelMean, 1, mean), "-")

##             AUSTR      BERG    BURIAT
## AUSTR    0.000000 10.531148  9.221329
## BERG   -10.531148  0.000000 -1.309819
## BURIAT  -9.221329  1.309819  0.000000

Estimates of the differences between the means of the means of groups by `fgender`

Those are partial effects of fgender given fpopul (if the additive model can be assumed). They are equal to differences in $\alpha^Z$ 's.

L.peffZ.sumXZ <- cbind(0, 0, 0, CZ.sum[2,] - CZ.sum[1,])
colnames(L.peffZ.sumXZ) <- names(coef(golXplusZAll.sum))
rownames(L.peffZ.sumXZ) <- paste(lgen[2], "-", lgen[1], sep = "")
print(L.peffZ.sumXZ)

##     (Intercept) fpopul1 fpopul2 fgender1
## M-F           0       0       0       -2

LSest(golXplusZAll.sum, L = L.peffZ.sumXZ)

##     Estimate Std. Error  t value    P value    Lower    Upper
## M-F 9.770313  0.7529092 12.97675 < 2.22e-16 8.288345 11.25228

Estimated means of the means (model-based) and their differences

apply(golModelMean, 2, mean)

##        F        M 
## 174.4870 184.2573

outer(apply(golModelMean, 2, mean), apply(golModelMean, 2, mean), "-")

##          F         M
## F 0.000000 -9.770313
## M 9.770313  0.000000

NMSA407 Linear Regression: Tutorial

Introduction

Full (unbalanced) dataset

Balanced subset (was artificially created)

Some colors (used in plots below)

Frequency tables

Full dataset

Balanced subset

Group sample means

Full dataset

OCA

GOL

Balanced subset

OCA

GOL

Boxplots (for full dataset)

OCA

GOL

Spaghetti plots (for full dataset)

F-tests and ANOVA tables with a full (unbalanced) dataset

Models for OCA

Type I ANOVA tables for models ocaXZAll and ocaXZAll2

Type II ANOVA table for model ocaXZAll

Type III ANOVA table for model ocaXZAll

Models for GOL

Type I ANOVA tables for models golXZAll and golXZAll2

Type II ANOVA table for model golXZAll

Type III ANOVA table for model golXZAll

F-tests and ANOVA tables with a balanced dataset

Models for OCA

Type I ANOVA tables for models ocaXZ and ocaXZ2

Type II ANOVA table for model ocaXZ

Type III ANOVA table for model ocaXZ

Models for GOL

Type I ANOVA tables for models golXZ and golXZ

Type II ANOVA table for model golXZ

Type III ANOVA table for model golXZ

Levels of the two factors

Reference group pseudocontrast parameterization contr.treatment

Pseudocontrast matrices

Matrix that (together with the intercept term) parameterizes the two-way classified group means

OCA (full dataset)

Interaction model

GOL (full dataset)

Interaction model

Additive model

Estimates of all (partial) effects of fpopul given fgender in the additive model

Estimates of all (partial) effects of fgender given fpopul in the additive model

Sum contrast parameterization contr.sum

Contrast matrices

Matrix that (together with the intercept term) parameterizes the two-way classified group means

OCA (full dataset)

Interaction model

Estimates of all αX\alpha^X parameters

Estimates of all αZ\alpha^Z parameters

Estimates of all αXZ\alpha^{XZ} parameters

Estimates of αXZ\alpha^{XZ} parameters stored in a matrix

Estimated group means

Estimates of the differences between the means of the means of groups by fpopul

Estimated means of the means and their differences

Estimates of the differences between the means of the means of groups by fgender

Estimated means of the means and their differences

GOL (full dataset)

Interaction model

Additive model

Estimates of all αX\alpha^X parameters in the additive model

Estimates of all αZ\alpha^Z parameters in the additive model

Estimated group means: model based

Estimated group means: sample means

Estimates of the differences between the means of the means of groups by fpopul

Estimated means of the means (model-based) and their differences

Estimates of the differences between the means of the means of groups by fgender

Estimated means of the means (model-based) and their differences

Type I ANOVA tables for models `ocaXZAll` and `ocaXZAll2`

Type II ANOVA table for model `ocaXZAll`

Type III ANOVA table for model `ocaXZAll`

Type I ANOVA tables for models `golXZAll` and `golXZAll2`

Type II ANOVA table for model `golXZAll`

Type III ANOVA table for model `golXZAll`

Type I ANOVA tables for models `ocaXZ` and `ocaXZ2`

Type II ANOVA table for model `ocaXZ`

Type III ANOVA table for model `ocaXZ`

Type I ANOVA tables for models `golXZ` and `golXZ`

Type II ANOVA table for model `golXZ`

Type III ANOVA table for model `golXZ`

Reference group pseudocontrast parameterization `contr.treatment`

Estimates of all (partial) effects of `fpopul` given `fgender` in the additive model

Estimates of all (partial) effects of `fgender` given `fpopul` in the additive model

Sum contrast parameterization `contr.sum`

Estimates of all $\alpha^X$ parameters

Estimates of all $\alpha^Z$ parameters

Estimates of all $\alpha^{XZ}$ parameters

Estimates of $\alpha^{XZ}$ parameters stored in a matrix

Estimates of the differences between the means of the means of groups by `fpopul`

Estimates of the differences between the means of the means of groups by `fgender`

Estimates of all $\alpha^X$ parameters in the additive model

Estimates of all $\alpha^Z$ parameters in the additive model

Estimates of the differences between the means of the means of groups by `fpopul`

Estimates of the differences between the means of the means of groups by `fgender`