Tartu et al. (2016) collected data on blood mercury concentrations in male and female Arctic black-legged kittiwakes (Rissa tridactyla) (Factor Sex, with two levels) that were collected during the incubation period and while chicks were being reared (Factor Breeding Stage, 2 levels). During incubation, they collected 48 females and 44 males, with 17 and 21, respectively, during chick rearing. The authors analyzed these data with a generalized linear model based on a normal distribution and identity link. In contrast, keeping with the theme for this chapter, we fitted a two-factor OLS linear model including the fixed main effects of site and season and their interaction. With these unequal sample sizes, we have the options of Type I, II or III SS.

Anderson, Brian, Public domain, via Wikimedia Commons

The paper is here

Tartu, S., Bustamante, P., Angelier, F., Lendvai, A. Z., Moe, B., Blevin, P., Bech, C., Gabrielsen, G. W., Bustnes, J. O. & Chastel, O. (2016). Mercury exposure, stress and prolactin secretion in an Arctic seabird: an experimental study. Functional Ecology, 30, 596-604.

Preliminaries

First, load the required packages (car)

Import tartu_hg data file (tartu_hg.csv)

tartu_hg <- read.csv("../data/tartu_hg.csv")
head(tartu_hg,10)

Examine boxplots by cell

boxplot(hg~sex*stage,data=tartu_hg)

No particular issues here, so happy to proceed

Fit OLS model (Type I SS) to untransformed data

tartu.aov <- aov(hg~sex*stage, data=tartu_hg)
summary(tartu.aov)
             Df Sum Sq Mean Sq F value   Pr(>F)    
sex           1  6.360   6.360  44.680 6.74e-10 ***
stage         1  6.978   6.978  49.022 1.34e-10 ***
sex:stage     1  0.394   0.394   2.766   0.0988 .  
Residuals   126 17.936   0.142                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
plot(tartu.aov)

Note pattern in residuals but untransformed data will be analysed

Refit model so stage goes first (still Type I SS)

tartu1.aov <- aov(hg~stage*sex, data=tartu_hg)
summary(tartu1.aov)
             Df Sum Sq Mean Sq F value   Pr(>F)    
stage         1  6.076   6.076  42.685 1.44e-09 ***
sex           1  7.262   7.262  51.017 6.48e-11 ***
stage:sex     1  0.394   0.394   2.766   0.0988 .  
Residuals   126 17.936   0.142                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
plot(tartu1.aov)

Now get Type II and Type III SS using car package

tartu.lm1 <- lm(hg~stage*sex, data=tartu_hg, contrasts=list(stage=contr.sum, sex=contr.sum))
Anova(tartu.lm1, type='III')
Anova Table (Type III tests)

Response: hg
             Sum Sq  Df   F value    Pr(>F)    
(Intercept) 167.105   1 1173.8984 < 2.2e-16 ***
stage         6.748   1   47.4012 2.438e-10 ***
sex           4.768   1   33.4957 5.345e-08 ***
stage:sex     0.394   1    2.7661   0.09877 .  
Residuals    17.936 126                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Anova(tartu.lm1, type='II')
Anova Table (Type II tests)

Response: hg
           Sum Sq  Df F value    Pr(>F)    
stage      6.9783   1 49.0222 1.341e-10 ***
sex        7.2624   1 51.0173 6.477e-11 ***
stage:sex  0.3938   1  2.7661   0.09877 .  
Residuals 17.9362 126                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
LS0tCnRpdGxlOiAiUUsgQm94IDcuMyIKCm91dHB1dDogCiAgaHRtbF9ub3RlYm9vazoKICAgIHRoZW1lOiBmbGF0bHkKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFKQpgYGAKClRhcnR1IGV0IGFsLiAoMjAxNikgY29sbGVjdGVkIGRhdGEgb24gYmxvb2QgbWVyY3VyeSBjb25jZW50cmF0aW9ucyBpbiBtYWxlIGFuZCBmZW1hbGUgQXJjdGljIGJsYWNrLWxlZ2dlZCBraXR0aXdha2VzICgqUmlzc2EgdHJpZGFjdHlsYSopIChGYWN0b3IgU2V4LCB3aXRoIHR3byBsZXZlbHMpIHRoYXQgd2VyZSBjb2xsZWN0ZWQgZHVyaW5nIHRoZSBpbmN1YmF0aW9uIHBlcmlvZCBhbmQgd2hpbGUgY2hpY2tzIHdlcmUgYmVpbmcgcmVhcmVkIChGYWN0b3IgQnJlZWRpbmcgU3RhZ2UsIDIgbGV2ZWxzKS4gRHVyaW5nIGluY3ViYXRpb24sIHRoZXkgY29sbGVjdGVkIDQ4IGZlbWFsZXMgYW5kIDQ0IG1hbGVzLCB3aXRoIDE3IGFuZCAyMSwgcmVzcGVjdGl2ZWx5LCBkdXJpbmcgY2hpY2sgcmVhcmluZy4gVGhlIGF1dGhvcnMgYW5hbHl6ZWQgdGhlc2UgZGF0YSB3aXRoIGEgZ2VuZXJhbGl6ZWQgbGluZWFyIG1vZGVsIGJhc2VkIG9uIGEgbm9ybWFsIGRpc3RyaWJ1dGlvbiBhbmQgaWRlbnRpdHkgbGluay4gSW4gY29udHJhc3QsIGtlZXBpbmcgd2l0aCB0aGUgdGhlbWUgZm9yIHRoaXMgY2hhcHRlciwgd2UgZml0dGVkIGEgdHdvLWZhY3RvciBPTFMgbGluZWFyIG1vZGVsIGluY2x1ZGluZyB0aGUgZml4ZWQgbWFpbiBlZmZlY3RzIG9mIHNpdGUgYW5kIHNlYXNvbiBhbmQgdGhlaXIgaW50ZXJhY3Rpb24uIFdpdGggdGhlc2UgdW5lcXVhbCBzYW1wbGUgc2l6ZXMsIHdlIGhhdmUgdGhlIG9wdGlvbnMgb2YgVHlwZSBJLCBJSSBvciBJSUkgU1MuCgpbIVtdKC4uL21lZGlhLzUxMnB4LUJsYWNrLWxlZ2dlZF9LaXR0aXdha2VfYW5kX0NoaWNrLmpwZyldKGh0dHBzOi8vY29tbW9ucy53aWtpbWVkaWEub3JnL3dpa2kvRmlsZTpCbGFjay1sZWdnZWRfS2l0dGl3YWtlX2FuZF9DaGljay5qcGcpCgpBbmRlcnNvbiwgQnJpYW4sIFB1YmxpYyBkb21haW4sIHZpYSBXaWtpbWVkaWEgQ29tbW9ucwoKVGhlIHBhcGVyIGlzIFtoZXJlXShodHRwczovL2RvaS5vcmcvMjAxMC4xMTExLzEzNjUtMjQzNS4xMjUzNCkKClRhcnR1LCBTLiwgQnVzdGFtYW50ZSwgUC4sIEFuZ2VsaWVyLCBGLiwgTGVuZHZhaSwgQS4gWi4sIE1vZSwgQi4sIEJsZXZpbiwgUC4sIEJlY2gsIEMuLCBHYWJyaWVsc2VuLCBHLiBXLiwgQnVzdG5lcywgSi4gTy4gJiBDaGFzdGVsLCBPLiAoMjAxNikuIE1lcmN1cnkgZXhwb3N1cmUsIHN0cmVzcyBhbmQgcHJvbGFjdGluIHNlY3JldGlvbiBpbiBhbiBBcmN0aWMgc2VhYmlyZDogYW4gZXhwZXJpbWVudGFsIHN0dWR5LiAqRnVuY3Rpb25hbCBFY29sb2d5KiwgMzAsIDU5Ni02MDQuCgojIyMgUHJlbGltaW5hcmllcwoKRmlyc3QsIGxvYWQgdGhlIHJlcXVpcmVkIHBhY2thZ2VzIChjYXIpCgpgYGB7ciBpbmNsdWRlPUZBTFNFLCByZXN1bHRzPSdoaWRlJ30Kc291cmNlKCIuLi9SL2xpYnJhcmllcy5SIikgICAjVGhpcyBpcyB0aGUgY29tbW9uIGxpYnJhcnkKYGBgCgpJbXBvcnQgdGFydHVfaGcgZGF0YSBmaWxlICh0YXJ0dV9oZy5jc3YpCgpgYGB7cn0KdGFydHVfaGcgPC0gcmVhZC5jc3YoIi4uL2RhdGEvdGFydHVfaGcuY3N2IikKaGVhZCh0YXJ0dV9oZywxMCkKYGBgCgojIyMgRXhhbWluZSBib3hwbG90cyBieSBjZWxsCgpgYGB7ciB9CmJveHBsb3QoaGd+c2V4KnN0YWdlLGRhdGE9dGFydHVfaGcpCmBgYAoKTm8gcGFydGljdWxhciBpc3N1ZXMgaGVyZSwgc28gaGFwcHkgdG8gcHJvY2VlZAoKIyMgRml0IE9MUyBtb2RlbCAoVHlwZSBJIFNTKSB0byB1bnRyYW5zZm9ybWVkIGRhdGEKCmBgYHtyIH0KdGFydHUuYW92IDwtIGFvdihoZ35zZXgqc3RhZ2UsIGRhdGE9dGFydHVfaGcpCnN1bW1hcnkodGFydHUuYW92KQpwbG90KHRhcnR1LmFvdikKYGBgCgpOb3RlIHBhdHRlcm4gaW4gcmVzaWR1YWxzIGJ1dCB1bnRyYW5zZm9ybWVkIGRhdGEgd2lsbCBiZSBhbmFseXNlZAoKIyMjIFJlZml0IG1vZGVsIHNvIHN0YWdlIGdvZXMgZmlyc3QgKHN0aWxsIFR5cGUgSSBTUykKCmBgYHtyIH0KdGFydHUxLmFvdiA8LSBhb3YoaGd+c3RhZ2Uqc2V4LCBkYXRhPXRhcnR1X2hnKQpzdW1tYXJ5KHRhcnR1MS5hb3YpCnBsb3QodGFydHUxLmFvdikKYGBgCgojIyMgTm93IGdldCBUeXBlIElJIGFuZCBUeXBlIElJSSBTUyB1c2luZyBjYXIgcGFja2FnZQoKYGBge3IgfQp0YXJ0dS5sbTEgPC0gbG0oaGd+c3RhZ2Uqc2V4LCBkYXRhPXRhcnR1X2hnLCBjb250cmFzdHM9bGlzdChzdGFnZT1jb250ci5zdW0sIHNleD1jb250ci5zdW0pKQpBbm92YSh0YXJ0dS5sbTEsIHR5cGU9J0lJSScpCkFub3ZhKHRhcnR1LmxtMSwgdHlwZT0nSUknKQpgYGAK