Fig. 5.5 Boxplot of standardised residuals obtained by a linear regression model applied on the log-transformed sibling negotiation data. The y-axis shows the values of the residuals and the horizontal axis the nests. Note that some nests have residuals that are above or below the zero line, indicating the need for a random effect a fixed term or as a random term, but we already discussed that this has to be as a random term.

5.10.2 Step 2 of the Protocol: Fit the Model with GLS

In this step we fit the model using the gls function. It allows us to compare the linear regression model with the mixed effects model that we will calculate using the lme function in a moment.

> library(nlme)

> Form <- formula(LogNeg ~ SexParent * FoodTreatment +

SexParent * ArrivalTime)

To reduce the code, we have used the formula expression. The numerical output in the object M.gls is identical to that of the lm function.

5.10.3 Step 3 of the Protocol: Choose a Variance Structure

In Chapter 4, this step consisted of finding the optimal variance structure in terms of heterogeneity. We can still do that here, but adding the random component nest is our first priority. Note that the random intercept is also part of the 'choose a variance structure' process. This means that the following random intercept mixed effects model is fitted.

LogNegj = a + jj1 x SexParentij + j2 x Foodtreatmentij

+ jj3 x ArrivalTimeij + j4 x SexParentij x FoodTreatmentj + j5 x SexParentj x ArrivalTimeij + ai + £„

LogNegj is the log-10 transformed sibling negotiation for observation j at nest i. SexParentij and FoodTreatmentij are nominal variables with two levels, and ArrivalTimeij is a continuous variable. The second line contains interactions. The term ai is a random intercept and is assumed to be normally distributed with mean 0 and variance d2. The residual eij is assumed to be normally distributed with mean 0 and variance a2. Both random terms are assumed to be independent of each other.

5.10.4 Step 4: Fit the Model

The linear mixed effects model is applied in R with the following code.

> M1.lme <- lme(Form, random = ~ 1 | Nest, method = "REML", data = Owls)

5.10.5 Step 5 of the Protocol: Compare New Model with Old Model

We use the anova command to compare the models M.gls and Mi.lme. Note that the models were estimated with REML, which allows us to apply the likelihood ratio test to see whether we need the random intercept.

Model df AIC BIC logLik Test L.Ratio p-value

M.gls 1 7 64.37422 95.07058 -25.18711

M1.lme 2 8 37.71547 72.79702 -10.85773 1 vs 2 28.65875 <.0001

The likelihood ratio test indicates that the model with the random intercept is considerably better. You would quote this statistic in a paper as L = 28.65 (df = 1, p < 0.001). Recall from Section 5.8 that we are testing on the boundary here. If we did the correction for testing on the boundary, the p-value would get even smaller. Because the random intercept is highly significant, testing on the boundary is not a problem here.

The AIC of the model with the random intercept is also considerably smaller, confirming the results of the likelihood ratio test. As well as the random intercept, it is also an option to use a random intercept and random slope model. In this case, you assume that the strength of the relationship between sibling negotiation and arrival time changes randomly between the nests. We leave this as an exercise to the reader.

0.2 0.3 0.4 Fitted values Food treatment

0.2 0.3 0.4 Fitted values Food treatment

Was this article helpful?

## Post a comment