Chapter 13 Analysis of covariance, nested data and mixed effects
13.1 Introduction
We’ve seen that regression is often not the best available technique to use for bivariate analysis.The books written by Zuur and associates show some of the big challenges involved in analysing real life ecological data. The idealised assumptions of linear regression are very rarely met in full.
Zuur writes Always apply the simplest statistical technique on your data, but ensure it is applied correctly! And here is a crucial problem. In ecology, the data are seldom modelled adequately by linear regression models. If they are, you are lucky. If you apply a linear regression model on your data, then you are implicitly assuming a whole series of assumptions, and once the results are obtained, you need to verify all of them. What should we do if we violate all the assumptions? The answer is simple: reject the model. But what do we do if we only violate one of the assumptions? And how much can we violate the assumptions before we are in trouble?
This sounds frightening. If model assumptions are rarely met, how can we ever trust them? Do we need the help of an expert statistician such as Alain Zuur to analyse all ecological data? Perhaps all we can use are simple tests. But do we really know that their assumptions hold?
The best advice to a student aiming to analyse data for an MSc dissertation is simply to always make the best possible attempt at diagnosing issues with the chosen analysis. Be aware of all the assumptions used and be critical at all times. But don’t despair if some assumptions are not met. Always point out the problems when discussing the analysis. In many cases the issues will in fact turn out to be unimportant. The key message held within most data sets can often be extracted, even if models have minor violations of assumptions. In other cases there may be more serious problems. Sometimes a more advanced analysis than was orginally planned may be necessary. If you have spotted a problem and understand its implications, it could be easier than you think to build a slightly more complex model. You will not be the first person to have come across the issue. The biggest mistakes usually involve failing to account for lack of independence in the data. This chapter runs through an example of this in practice. The data initially seem simple enough.
An important issue arises when the same relationship between two variables is repeated multiple times in a data set. In the past you may have handled this situation by subsetting the data various times and repeating the analysis for each subset. However it is possible to build more oomprehensive models that look at a population of relationships.
13.2 Whale teeth and isotope ratios
In Zuur et al. (2009) a rather complex analysis involving generalised additive models and correlated residuals is developed and justified as a means of ensuring that all model assumptions are met. The results of the analysis are also presented in a paper that shows how the statistical analysis was used to answer an interesting scientific question. It is well worth reading this paper before analysing the data.
The researchers aimed to piece together information regarding feeding behaviour and migration of Sperm whales, based on analysis of the composition of their teeth. Specifically they looked at carbon and nitrogen isotope ratios. N isotope ratios can be indicative of trophic position in marine communities. Trophic level is expected to be higher for larger individuals of the same species. Sperm whales are approximately 4m long at birth, with males growing to 18 m.The relative trophic position of the different stages in the life of several male sperm whales can be investigated through the use of N stable isotope ratios. It was expected that the trophic position would increase with age, as the animals become larger.
Whale teeth have growth rings. To obtain data the researchers analysed isotope ratios from the bands of individual teeth.
The first task in the analysis could be simply to establish that the expectation of increased N isotope ratio with age is supported by the data. We will concentrate on that before looking at more complex issues.
13.2.1 Moby’s tooth
Let’s first look at a single tooth that was obtained from a famous whale called Moby that was washed up on the shores of the Firth of Forth in 1997 after many rescue attempts. Moby was 15.2 metres in length and weighed 38.5 tons.
Moby’s bones went on display in the Royal Scottish Museum, Chamber’s street soon after death. The original display failed to attract the public due to the appalling smell.
Moby’s skull has since been de-oderised and is now redisplayed. It has even been associated with a Turner prize.
Moby is clearly a star turn. But, how does his life history compare to that of other Sperm Whales? This is the topic we will analyse in this class.
Whales<-read.csv("https://tinyurl.com/aqm-data/whaleteeth.csv")
Moby<-subset(Whales,Whale=="Moby")13.3 Plotting the data
We should first plot the data to look at the pattern. Notice the way a more complex label is set up for the y axis. There are many tricks such as this in R that are rather hard to discover. If you need to use Greek symbols on the axes of your own figures try adapting this code.
ylabel <-expression(paste(delta^{15}, "N"))
xlabel<-"Estimated age"
library(ggplot2)
theme_set(theme_bw())
g0<-ggplot(data=Moby,aes(x=Age,y=X15N))
g1<-g0+geom_point() +xlab(xlabel) +ylab(ylabel)
g1
13.4 Fitting a regression
OK, so this all looks quite simple at this stage. The data points lie more or less along a straight line. So we could fit a regression.
mod<-lm(data=Moby,X15N~Age)
summary(mod)##
## Call:
## lm(formula = X15N ~ Age, data = Moby)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.07102 -0.28706 0.04346 0.33820 1.13724
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.748940 0.163559 71.83 <2e-16 ***
## Age 0.113794 0.006186 18.40 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4859 on 40 degrees of freedom
## Multiple R-squared: 0.8943, Adjusted R-squared: 0.8917
## F-statistic: 338.4 on 1 and 40 DF, p-value: < 2.2e-16g1+stat_smooth(method="lm")+labs(y = ylabel,x=xlabel)
13.4.1 Diagnostics
Diagnostics show that most of the assumptions of a regression appear to be met. However an issue shows up when we look at the first diagnostic plot.
plot(mod,which=1)
There is a clear pattern in the residuals.
13.4.2 Testing for serial correlation
In fact the residuals are serially correlated. This can be checked using a Durbin Watson test.
library(lmtest)
dwtest(Moby$X15N~Moby$Age)##
## Durbin-Watson test
##
## data: Moby$X15N ~ Moby$Age
## DW = 1.1458, p-value = 0.0009539
## alternative hypothesis: true autocorrelation is greater than 0We can confirm this using the acf function.
acf(residuals(mod))
If any of the lines apart from the first extends above or below the blue line the data are significantly serially correlated. In other words the last value can be used to partially predict the next value. If one residual is negative, the next is also likely to be negative and vice versa. The correlation is significant at lag 1 and lag 2.
13.5 Interpreting the results
Does the serial correlation matter? Well in a purely statistical sense it could matter quite a lot. Zuur et al take the effect into account through the use of a rather cunning statistical trick using an ARIMA model. However most MSc students would be unlikely to be aware of this rather advanced method. In this context it is not really worth the effort. The only real difference between a model that includes autocorrelation and one that does not involves the size of the p-value and the width of the confidence bands. Providing that a straight line is a reasonable description of the underlying pattern (and we will come on to that later), the model is not changed.
There is something much more important to consider however. We have spotted the serial correlation (lack of independence) between years that occurs in the case of Moby. However, you may have realised that in reality all the data we have just analysed are totally lacking in independence anyway! They have been taken from a single whale. So, in one sense we have been conducting statistics with a sample size of one.
In effect, the relationship that we have established applies only to Moby. It would also be really nice to have data from some other teeth from Moby to check even this. We certainly cannot generalise it to other whales. We may have found out some interesting information along the way. The researchers were interested in understanding the trophic ecology of mature sperm whales in some detail. If they are right in assuming that high levels of delta N are associated with feeding on larger prey then there may be some indication that Moby spent several years in a row feeding on prey that may have been larger than expected for his age, followed by a switch to rather smaller prey. This is quite a speculative interpretation, but it may tell us something.
13.6 Exercise
Now try implementing an appropriate method to look at alternatives to a straight line to capture the empirical relationship between age and delta N measurements in Moby’s tooth. Comment on your findings.
13.7 Finding a general pattern
We actually have data for fifteen whales. So we could compare the pattern shown by Moby with the population of whale teeth. So let’s plot out the data from all fifteen whales at once.
g0<-ggplot(data=Whales,aes(x=Age,y=X15N,col=Whale))
g1<-g0+geom_point() + xlab(xlabel) +ylab(ylabel)
g1
There is still a clear pattern, but with much more scatter this time. This is to be expected, as each whale has its own life history that contributes variability.
We can look at this by plotting out the data for each individual whale using a facet wrap to split the data.
g0<-ggplot(data=Whales,aes(x=Age,y=X15N))
g1<-g0+geom_point()+geom_line()+ labs(y = ylabel,x=xlabel)
g1+facet_wrap("Whale")
Each whale seems to have its own pattern. Most are approximately linear, like Moby, but some show a very different pattern
For the moment we will ignore this and continue to fit a simple regression model to all the data. Note THIS IS THE WRONG WAY TO WORK WITH THE DATA!
mod1<-lm(X15N~Age,data=Whales)
summary(mod1)##
## Call:
## lm(formula = X15N ~ Age, data = Whales)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5041 -0.7013 -0.1141 0.6209 3.3444
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.228411 0.108598 112.6 <2e-16 ***
## Age 0.095715 0.005631 17.0 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9953 on 305 degrees of freedom
## Multiple R-squared: 0.4865, Adjusted R-squared: 0.4848
## F-statistic: 288.9 on 1 and 305 DF, p-value: < 2.2e-16confint(mod1)## 2.5 % 97.5 %
## (Intercept) 12.01471479 12.4421063
## Age 0.08463424 0.1067958The problem with this model is that we have far too many degrees of freedom. The model has assumed that each data point is independent, which clearly is not true.
13.8 Analysis of covariance
The data that we have available to predict N15 levels consists of two variables. One is Age and the other is the identity of the individual whale from which the data was obtained. This is a categorical variable. If we continue to put to one side the issue of independence of the data points within the regression for each whale we can analyse the data using a technique known as analysis of covariance.
Analysis of covariance takes into account the differences between each group of observations. It is used to answer the following questions.
- Does the level of a categorical factor affect the response variable?
- Does the slope of a regression line differ between levels of the categorical variable?
In our case the first question implies that each whale may have a different mean value for \(\delta^{15}N\). If you read the paper carefully you will see that the researchers were aware that this is a rather trivial question. They looked at more subtle differences related to time of birth and life history. However we will ignore this for the moment too and carry out the analysis in order to illustrate the point.
As mentioned previously, R has a very consistent syntax that makes model fitting very easy. All we need to do to add in an effect for each whale to the model is to write lm(X15N~Age+Whale). The line of R code preceding this sets the contrasts for the model. This determines the output from the summary of the model. In this case we want to contrast the results for other whales with Moby (number 11 in the list).
contrasts(Whales$Whale)<-contr.treatment(levels(Whales$Whale),base=11)
mod2<-lm(X15N~Age+Whale,data=Whales)
anova(mod2)## Analysis of Variance Table
##
## Response: X15N
## Df Sum Sq Mean Sq F value Pr(>F)
## Age 1 286.21 286.215 431.039 < 2.2e-16 ***
## Whale 10 106.27 10.627 16.004 < 2.2e-16 ***
## Residuals 295 195.88 0.664
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod2)##
## Call:
## lm(formula = X15N ~ Age + Whale, data = Whales)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.70699 -0.51434 -0.07552 0.47030 2.64139
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.256788 0.174661 70.175 < 2e-16 ***
## Age 0.092183 0.005159 17.869 < 2e-16 ***
## WhaleI1/98 -0.010461 0.181371 -0.058 0.95404
## WhaleM143/96D -0.229444 0.220575 -1.040 0.29909
## WhaleM143/96E 0.385063 0.235813 1.633 0.10355
## WhaleM2583/94(1) -0.936579 0.216095 -4.334 2.01e-05 ***
## WhaleM2583/94(10) -0.357032 0.209785 -1.702 0.08983 .
## WhaleM2583/94(7) -0.901441 0.223214 -4.038 6.87e-05 ***
## WhaleM2679/93 -0.059622 0.185189 -0.322 0.74772
## WhaleM2683/93 0.267238 0.227314 1.176 0.24069
## WhaleM447/98 1.210339 0.192499 6.288 1.16e-09 ***
## WhaleM546/95 0.685273 0.219544 3.121 0.00198 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8149 on 295 degrees of freedom
## Multiple R-squared: 0.6671, Adjusted R-squared: 0.6547
## F-statistic: 53.73 on 11 and 295 DF, p-value: < 2.2e-16The anova result is very clear. There is a significant additive effect that can be attributed to the individual. So we can see that whales do differ in their baseline N15 levels. In order to interpret the model summary you need to be aware that the intercept represents the value at age zero for Moby, which was set to be the “control” or baseline for the contrasts. The other whale values are compared to the baseline. Thus if a coefficient for an individual whale is not significant it means that the intercept for that whale was not different to that of Moby.
So, the answer to the first question, just as we would expect, is a definite yes. Six of the ten whales that we compared with Moby have significantly different levels of N15 after allowing for age.
How do we answer the second question? This involves adding an interaction term to the model. The interaction implies that each whale could have it’s own intercept AND it’s own slope. If the interaction term is significant then we have a very complex model that cannot be easily generalised.
mod3<-lm(X15N~Age*Whale,data=Whales)
anova(mod3)## Analysis of Variance Table
##
## Response: X15N
## Df Sum Sq Mean Sq F value Pr(>F)
## Age 1 286.215 286.215 734.509 < 2.2e-16 ***
## Whale 10 106.266 10.627 27.271 < 2.2e-16 ***
## Age:Whale 10 84.828 8.483 21.769 < 2.2e-16 ***
## Residuals 285 111.055 0.390
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The interaction term is highly significant. This is not suprising when we consider the lattice plot, which showed a series of very different patterns for each whale.
If we summarise the model now we will obtain even more complex output which shows how the coefficients depend on the identity of each whale.
summary(mod3)##
## Call:
## lm(formula = X15N ~ Age * Whale, data = Whales)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.58547 -0.40951 -0.01552 0.37746 2.41087
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.748940 0.210125 55.914 < 2e-16 ***
## Age 0.113794 0.007947 14.320 < 2e-16 ***
## WhaleI1/98 1.464565 0.303845 4.820 2.34e-06 ***
## WhaleM143/96D 1.905086 0.376915 5.054 7.73e-07 ***
## WhaleM143/96E 0.720710 0.364725 1.976 0.049115 *
## WhaleM2583/94(1) -1.517527 0.336649 -4.508 9.58e-06 ***
## WhaleM2583/94(10) 0.747337 0.328174 2.277 0.023512 *
## WhaleM2583/94(7) -0.559979 0.346500 -1.616 0.107179
## WhaleM2679/93 -0.750516 0.296705 -2.529 0.011962 *
## WhaleM2683/93 0.509286 0.412830 1.234 0.218351
## WhaleM447/98 1.557803 0.305940 5.092 6.45e-07 ***
## WhaleM546/95 3.269282 0.341382 9.577 < 2e-16 ***
## Age:WhaleI1/98 -0.065573 0.011918 -5.502 8.36e-08 ***
## Age:WhaleM143/96D -0.142106 0.022432 -6.335 9.23e-10 ***
## Age:WhaleM143/96E -0.004390 0.027327 -0.161 0.872476
## Age:WhaleM2583/94(1) 0.065493 0.020050 3.267 0.001222 **
## Age:WhaleM2583/94(10) -0.065797 0.018155 -3.624 0.000343 ***
## Age:WhaleM2583/94(7) -0.007142 0.022432 -0.318 0.750421
## Age:WhaleM2679/93 0.041481 0.012471 3.326 0.000996 ***
## Age:WhaleM2683/93 -0.001922 0.025478 -0.075 0.939929
## Age:WhaleM447/98 -0.012176 0.013906 -0.876 0.381993
## Age:WhaleM546/95 -0.194624 0.021171 -9.193 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6242 on 285 degrees of freedom
## Multiple R-squared: 0.8112, Adjusted R-squared: 0.7973
## F-statistic: 58.33 on 21 and 285 DF, p-value: < 2.2e-16Now we can see that six of the whales do not just have a different baseline level of \(\delta^{15}N\) to Moby, the change over time when modelled as a linear response is also different for these individuals.
13.8.1 More about interactions
The model above could also be written in R longhand as:
mod3<-lm(X15N~Age+Whale+Age:Whale,data=Whales)The {*} symbol is simply short hand for this. It does not imply multiplication. The interaction component is Age:Whale. Don’t confuse statistical interactions with ecological interactions such so those that occur in a food web. A statistical interaction implies that we need to know additional information in order to predict something. For example if an experiment showed an interaction between nitrogen fertilizer level and soil type we could not be sure that adding nitrogen always enhances growth in the same way. Perhaps it has no effect at all when applied to heavy clay soils, although it increases growth generally. Finding interactions suggest complexity. This can often be interesting in itself.
13.8.2 Plotting the model
The model with interactions can be plotted using ggplots as a series of individual regression lines for each whale. Notice that the slope differs between each of the panels.
g0<-ggplot(data=Whales,aes(x=Age,y=X15N,group=Whale))
g1<-g0+geom_point()+ labs(y = ylabel,x=xlabel) +geom_smooth(method="lm")
g1+facet_wrap("Whale")