QK Box 8.16

Yavno and Fox (2013) collected introduced pumpkinseed sunfish from four sites in Canada, two flowing water sites and two still water sites, and returned them to a holding facility. Progeny from mating between individuals within each population were then allocated to artificial channels and subjected to either still or flowing water for 80 days. Various morphological features for each fish were measured at the end of the experiment. Each morphological measurement was analyzed in a two-factor linear model, including the value of the centroid of nine points on the body as a measure of body size. Each morphological variable and the value of the centroid was log transformed before analysis.

Pumpkinseed sunfish, Lepomis gibbosus. © Smithsonian Environmental Research Center, CC BY 4.0

The paper is here and the data are in the supplementary material S2

Yavno, S. & Fox, M. G. (2013). Morphological change and phenotypic plasticity in native and non-native pumpkinseed sunfish in response to sustained water velocities. Journal of Evolutionary Biology, 26, 2383-95.

Preliminaries

First, load the required packages (car, interactions)

Import yavno data file (yavno.csv)

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

Check that the covariate is not affected by factors

boxplot(log10(centroid)~pops*watervel, data=yavno)

cov.aov <- aov(log10(centroid)~pops*watervel, data=yavno)
summary(cov.aov)

               Df Sum Sq Mean Sq F value   Pr(>F)    
pops            3 0.4780 0.15932  34.481  < 2e-16 ***
watervel        1 0.0730 0.07301  15.801 9.23e-05 ***
pops:watervel   3 0.0461 0.01538   3.329   0.0202 *  
Residuals     248 1.1459 0.00462                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Covariate affected by both factors and interaction using the anova result, but we’re dealing with quite a sensitive test, give the residual df. Looking at the boxplots, however, the covariate ranges are broadly overlapping, and one group is a bit higher than the rest. It’s not a clear-cut decision, but we’ll do ancova anyway to match original paper.

First bodydep

boxplot(log10(bodydep)~pops*watervel, data=yavno)

scatterplot(log10(bodydep)~log10(centroid)|pops, data=yavno)

scatterplot(log10(bodydep)~log10(centroid)|watervel, data=yavno)

all looks OK

Fit model using lm to get Type III SS

bodydep.lm1 <- lm(log10(bodydep)~pops*watervel*log10(centroid), contrasts=list(pops=contr.sum,watervel=contr.sum), data=yavno)
plot(bodydep.lm1)

Anova(bodydep.lm1, type='III')

Anova Table (Type III tests)

Response: log10(bodydep)
                               Sum Sq  Df   F value    Pr(>F)    
(Intercept)                   0.05549   1  519.7533 < 2.2e-16 ***
pops                          0.00019   3    0.5874 0.6238187    
watervel                      0.00143   1   13.3754 0.0003133 ***
log10(centroid)               0.34148   1 3198.3933 < 2.2e-16 ***
pops:watervel                 0.00023   3    0.7295 0.5353057    
pops:log10(centroid)          0.00017   3    0.5310 0.6614451    
watervel:log10(centroid)      0.00142   1   13.2907 0.0003270 ***
pops:watervel:log10(centroid) 0.00025   3    0.7924 0.4991730    
Residuals                     0.02562 240                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Combine covariate by factor terms into a single term

bodydep.lm2 <- lm(log10(bodydep)~pops*watervel+log10(centroid), contrasts=list(pops=contr.sum,watervel=contr.sum), data=yavno)
anova(bodydep.lm1, bodydep.lm2)

Analysis of Variance Table

Model 1: log10(bodydep) ~ pops * watervel * log10(centroid)
Model 2: log10(bodydep) ~ pops * watervel + log10(centroid)
  Res.Df      RSS Df Sum of Sq      F  Pr(>F)  
1    240 0.025624                              
2    247 0.027431 -7 -0.001807 2.4178 0.02074 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Anova(bodydep.lm2, type='III')

Anova Table (Type III tests)

Response: log10(bodydep)
                 Sum Sq  Df   F value    Pr(>F)    
(Intercept)     0.06383   1  574.7058 < 2.2e-16 ***
pops            0.00420   3   12.6181 1.058e-07 ***
watervel        0.00008   1    0.7637    0.3830    
log10(centroid) 0.38831   1 3496.4440 < 2.2e-16 ***
pops:watervel   0.00050   3    1.5072    0.2132    
Residuals       0.02743 247                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Plot of slopes

We’ll create log10(centroid) as a new variable, rather than transforming it repeatedly in models.

yavno$lcentroid <- log10(yavno$centroid)
bodydep.lm3 <- lm(log10(bodydep)~pops*watervel*lcentroid, data=yavno)
interact_plot(bodydep.lm3, pred=lcentroid, modx=pops, plot.points=TRUE, data=yavno)

interact_plot(bodydep.lm3, pred=lcentroid, modx=watervel, plot.points=TRUE, data=yavno)

Unequal slopes apparent for water velocity groups, but the sample size is large and the difference in slopes isn’t great, so fit homogeneous slopes ancova model anyway

bodydep.lm4 <- lm(log10(bodydep)~pops+watervel+pops*watervel+log10(centroid), contrasts=list(pops=contr.sum,watervel=contr.sum), data=yavno)
Anova(bodydep.lm4, type='III')

Anova Table (Type III tests)

Response: log10(bodydep)
                 Sum Sq  Df   F value    Pr(>F)    
(Intercept)     0.06383   1  574.7058 < 2.2e-16 ***
pops            0.00420   3   12.6181 1.058e-07 ***
watervel        0.00008   1    0.7637    0.3830    
log10(centroid) 0.38831   1 3496.4440 < 2.2e-16 ***
pops:watervel   0.00050   3    1.5072    0.2132    
Residuals       0.02743 247                        
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Get adjusted means

adjmeans <- effect("pops*watervel", bodydep.lm4, se=TRUE)
summary(adjmeans)

 pops*watervel effect
          watervel
pops           Flow    Still
  Otonabee 1.195001 1.191510
  Rice     1.203540 1.205033
  Susqueda 1.198707 1.201057
  Ter      1.200549 1.204955

 Lower 95 Percent Confidence Limits
          watervel
pops           Flow    Still
  Otonabee 1.191188 1.187786
  Rice     1.199925 1.201094
  Susqueda 1.194843 1.197484
  Ter      1.196344 1.201321

 Upper 95 Percent Confidence Limits
          watervel
pops           Flow    Still
  Otonabee 1.198813 1.195233
  Rice     1.207155 1.208971
  Susqueda 1.202571 1.204631
  Ter      1.204754 1.208590