QK Box 8.1

Ercit et al. (2014) focused on whether egg number influenced jump distance of female crickets, taking into account differences in body size (pronotum length) and leg strength (leg mass). They initially fitted a model that included all three predictors plus interaction terms. However, the sample size (n=29) is relatively small for a model to include 7 parameters plus the intercept so we will focus just on the main effects of the three predictors.

Black-horned_Tree_Cricket (Oecanthus nigricornis). Andrew C, CC BY 2.0, via Wikimedia Commons

The paper is here

Ercit, K., Martinez-Novoa, A. & Gwynne, D. T. (2014). Egg load decreases mobility and increases predation risk in female black-horned tree crickets (Oecanthus nigricornis). PLoS One, 9, e110298.

Preliminaries

First, load the required packages (car, lm.beta)

Import ercit data file (ercit.csv)

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

Run diagnostics

Do scatterplotMatrix for all variables

scatterplotMatrix(~distance+eggs+pronotum+legmass, data=ercit, regLine=FALSE, diagonal=list(method='boxplot'))

All variables look symmetrically distributed, no obvious outliers, relationships between distance and each predictor show no signs of non-linearity.

Now check for collinearity among predictors - correlations and scatterplots and VIFs

cor(ercit[,c('eggs','pronotum','legmass')])

          eggs pronotum legmass
eggs     1.000    0.185   0.153
pronotum 0.185    1.000   0.575
legmass  0.153    0.575   1.000

vif(lm(distance~ eggs+pronotum+legmass, data=ercit))

    eggs pronotum  legmass 
    1.04     1.52     1.50 

Some collinearity between pronotum length and egg mass but r < 0.7 and no VIFs >5 so proceed (will do re-analysis later)

Fit model with 3 predictors

ercit.lm <- lm(distance~ eggs+pronotum+legmass, data=ercit)

This command fits this model:

\[ \operatorname{distance} = \alpha + \beta_{1}(\operatorname{eggs}) + \beta_{2}(\operatorname{pronotum}) + \beta_{3}(\operatorname{legmass}) + \epsilon \]

Examine regression diagnostics (residual plot and Cooks D)

plot(ercit.lm)

augment(ercit.lm)

Residual plot OK, no unusually large Cook D values

Display results of model fitting

We’re comfortable with the diagnostics, so we can now assess the parameters:

\[ \operatorname{\widehat{distance}} = -19.06 - 0.11(\operatorname{eggs}) + 18.12(\operatorname{pronotum}) + 0.38(\operatorname{legmass}) \]

options(digits = 3)
tidy(ercit.lm, conf.int = TRUE)

glance(ercit.lm)

Get standarised regression coefficients

Use lm.beta function from lm.beta package

lm.beta.ercit <- lm.beta(ercit.lm)
tidy(lm.beta.ercit, conf.int = TRUE)

glance(lm.beta.ercit)

Show the partial regression (added-variable) plots

avPlots(ercit.lm, ask=F)

Illustrate model comparison approach

There is also an issue with collinearity between pronotum and leg mass. While the collinearity measure does not exceed our threshold, we might also be interested in models with only one of this pair. We’ll drop leg mass and then omit pronotum.

# Omit legg mass from model
ercit1.lm <- lm(distance~ eggs+pronotum, data=ercit)
tidy(ercit1.lm, conf.int = TRUE)

lm.beta.ercit1 <- lm.beta(ercit1.lm)
tidy(lm.beta.ercit1, conf.int = TRUE)

# Omit pronotum from model
ercit2.lm <- lm(distance~ eggs+legmass, data=ercit)
tidy(ercit2.lm, conf.int = TRUE)

lm.beta.ercit2 <- lm.beta(ercit2.lm)
tidy(lm.beta.ercit2, conf.int = TRUE)

Illustration of partial regression coefficient for egg number (Box 8.5)

We need two steps:

1. Variation in egg number that is independent of pronotum and legmass. To get this, we fit a linear model where egg number is the response and pronotum and legmass the predictors. We want the residuals from this regression.

2. Fit a regression model linking jump distance to the residuals from step 1. The slope of this regression is the partial regression coefficient

ercit_egg<-lm(eggs~pronotum+legmass,data=ercit)
model1<-lm(ercit$distance~ercit_egg$residuals)
tidy(model1)

Illustration of model comparison approach (Box 8.6)

First test overall H₀ that all regression coefficients equal zero, then test H₀s about individual regression coefficients (egg number for example)

#Overall test
ercit.lm <- lm(distance~ eggs+pronotum+legmass, data=ercit)
ercit1.lm <- lm(distance~1, data= ercit)
anova(ercit.lm, ercit1.lm)

Analysis of Variance Table

Model 1: distance ~ eggs + pronotum + legmass
Model 2: distance ~ 1
  Res.Df RSS Df Sum of Sq    F  Pr(>F)    
1     25 324                              
2     28 736 -3      -413 10.6 0.00011 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Now drop egg number from model
ercit2.lm <- lm(distance~pronotum+legmass, data= ercit)
#Compare the two models
anova(ercit.lm, ercit2.lm)

Analysis of Variance Table

Model 1: distance ~ eggs + pronotum + legmass
Model 2: distance ~ pronotum + legmass
  Res.Df RSS Df Sum of Sq    F  Pr(>F)    
1     25 324                              
2     26 666 -1      -342 26.4 2.6e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1