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

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

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




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

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

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




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

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