11 Correlation structures
NoteIn this chapter
In this chapter, you will learn
- how to distinguish a trend from a true residual autocorrelation, and how to add a correlation structure to
glsorglmmTMB(Section 11.1); - how to detect and model temporal, spatial and phylogenetic correlation (Section 11.2, Section 11.3, Section 11.4);
- how multivariate GLMs describe the correlation between several responses (Section 11.6).
Workflow steps (Chapter 12): step 3 (correlation structures from the design), step 5 (residual checks: autocorrelation).
Prerequisites: random effects (Chapter 4), modeling the variance (Section 10.1.2), GLMs (Chapter 9).
11.1 General idea
Except for the random effects (Chapter 4), we have so far assumed that residual errors are independent. However, this is not always the case - we may find that residuals show autocorrelation.
Note
Correlation means that one variable correlates with another. Autocorrelation means that data points of one variable that are close to each other are more similar (positive autocorrelation, the usual case) or more different (negative autocorrelation) than expected for independent data. This implies that autocorrelation is only defined if there is a “distance relationship” between observations.
Autocorrelation can always occur if we have a distance relationship between observations. In random effects, distance is expressed by group, and all members of a group are equally correlated. Common examples of continuous distances, where the correlation decays with distance, include:
- Spatial distance.
- Temporal distance.
- Phylogenetic distance.
Here a visualization from Roberts et al. (2017) (reproduced as OA, copyright: the authors).
Roberts, D.R. et al. (2017) “Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure,” Ecography, 40(8), pp. 913–929.
Tip 11.1: The “autoregressive” model zoo: SAR, CAR, geostatistical, and phylogenetic models
The word “autoregressive” is used for several related, but mathematically distinct, constructions. It’s worth knowing the differences, even though the practical recipe (de-trend first, then add a correlation structure) is the same for all of them.
(Simultaneous) autoregressive models (SAR): the classic time-series AR(1) model, \(x_{t+1} = a \cdot x_t + \epsilon\), defines each value recursively from the previous one plus noise. This is what
corAR1/ar1()implement, and the same idea generalizes to space as spatial SAR models (\(Y = BY + \epsilon\)).Conditional autoregressive models (CAR) (Besag, 1974): instead of a recursive equation, you specify the distribution of each point conditional on its neighbors, \(Y_i \mid Y_{-i} \sim N(\sum_{j \sim i} c_{ij} Y_j,\ \tau_i^2)\). This requires a discrete neighborhood structure (e.g. an adjacency matrix for areal/lattice data, as used in packages like
CARBayesor INLA’sbesagmodel) and is common in spatial epidemiology. CAR and SAR are both Markovian and closely related, but they are built differently and generally do not induce the same joint distribution.Geostatistical / covariance-function models: the spatial example used later in this chapter (Section 11.3;
corExp,exp()in glmmTMB) doesn’t use a discrete neighborhood at all. Instead, the correlation between any two points is a direct function of their (continuous) distance, e.g. exponential decay. This “kriging” / Gaussian-process style is a third, separate family - despite being spatial, it is technically neither SAR nor CAR.Phylogenetic models (Brownian motion, Ornstein-Uhlenbeck): perhaps surprisingly, these are close relatives of the SAR family, not CAR. Brownian motion is a continuous-time random walk (the \(a \to 1\) limit of AR(1)), and the Ornstein-Uhlenbeck process used for
corMartins/corPagelis exactly the continuous-time analog of AR(1): its exact discrete-time solution is \(x_{t+1} = \mu + e^{-\theta \Delta t}(x_t - \mu) + \epsilon_t\), an AR(1) recursion with \(a = e^{-\theta \Delta t}\).
You don’t need to memorize this taxonomy for the course, but if you go looking specifically for “CAR models” in the literature or in packages like INLA, be aware that this refers to the Besag-style conditional construction - not to the temporal, phylogenetic, or corExp-style spatial models we actually use in this chapter.
Besag, J. (1974) “Spatial interaction and the statistical analysis of lattice systems,” Journal of the Royal Statistical Society: Series B (Methodological), 36(2), pp. 192–225.
11.1.1 Models to deal with autocorrelation
If we find autocorrelation in the residuals of our model, there can be several reasons, which we can address by different structures.
Important
In the context of regression models, we are never interested in the autocorrelation of the response / predictors per se, but only in the residuals. Thus, it doesn’t make sense to assume that you need a spatial model only because you have a spatially autocorrelated signal.
Autocorrelation can occur because we have a correlated misfit, i.e. there is a trend along the distance axis (time, space, or phylogeny). If this is the case, de-trending the model (with a linear regression term or a spline) will remove the residual autocorrelation. We should always de-trend first before we consider moving to a model with a residual correlation structure.
Only after accounting for the trend, we should test if there is a residual spatial / temporal / phylogenetic autocorrelation. If that is the case, we would usually use a so-called autoregressive structure (see Tip 11.1 for the different flavors of these models). In these models, we make parametric assumptions for how the correlation between data points falls off with distance. When we speak about spatial / temporal / phylogenetic regressions, we usually refer to these or similar models.
11.1.2 R implementation
De-trending needs no new tools, only regression terms or splines on time or space, so the rest of this chapter concentrates on how to specify “real” correlation structures.
To account for “real” autocorrelation of residuals, similarly to the variance modeling (Section 10.1.2), we can add correlation structures
- for normal responses in
nlme::gls, see https://stat.ethz.ch/R-manual/R-devel/library/nlme/html/corClasses.html - for GLMs using
glmmTMB, see https://cran.r-project.org/web/packages/glmmTMB/vignettes/covstruct.html.
The following sections (Section 11.2, Section 11.3, Section 11.4) provide examples and further comments on how to do this.
Note
Especially for spatial models, both nlme and glmmTMB are relatively slow. There are a large number of specialized packages that deal in particular with the problem of spatial models, including MASS::glmmPQL, BRMS, INLA, spaMM, and many more. Bayesian autoregressive models, and state-space models that separate process and observation error, are covered in Hartig (2024) (autoregressive models, state-space models). To keep things simple and concentrate on the principles, however, we will stick with the packages you already know.
Hartig, F. (2024) “Introduction to Bayesian statistics with R.” Available at: https://florianhartig.github.io/LearningBayes/.
11.2 Temporal correlation structures
To introduce temporal autoregressive models, let’s simulate some data first. The most simple (and common) temporal structure is the AR1 model, aka autoregressive model with lag 1. The AR1 assumes that the next data point (or residual) originates from the last data point, multiplied by a factor a (with |a| < 1), plus a normally distributed residual, in the form
\[x_{t+1} = a \cdot x_t + \epsilon\]
Let’s simulate some data according to this model
# simulate temporally autocorrelated data
AR1sim<-function(n, a){
x = rep(NA, n)
x[1] = 0
for(i in 2:n){
x[i] = a * x[i-1] + rnorm(1)
}
return(x)
}
set.seed(123)
obs = AR1sim(1000, 0.9)
plot(obs)
As we can see, we have a temporal correlation here. As we have not modeled / specified any further predictors, the correlation in the signal will transform directly in the correlations of the residuals if we fit a model:
fit = lm(obs~1)
summary(fit)
Call:
lm(formula = obs ~ 1)
Residuals:
Min 1Q Median 3Q Max
-5.8485 -1.4315 0.0559 1.3729 6.6317
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.17019 0.06745 2.523 0.0118 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.133 on 999 degrees of freedomNote that the estimate of the intercept is significant, although we started the simulation at zero: because neighboring points are similar, the data contain less independent information than n = 1000 suggests, so the standard error of the lm is too small and its p-value too optimistic. Let’s look at the residuals, which have the same autocorrelation as the data.
plot(residuals(fit))
We can quantify the autocorrelation by the acf function, which quantifies correlation between observations as a function of their lag (temporal distance). Note that although we modeled only a lag of 1, we will get correlations with many lags, because the correlation effect “trickles down”.
acf(residuals(fit))
To check what the actual underlying model is, it may be useful to plot the partial correlation coefficient, i.e. the correlation at lag k that remains after accounting for lags 1 to k-1 (the last coefficient of an AR(k) model).
pacf(residuals(fit))
Here, we see that we actually only have a correlation with lag 1. You can also check for temporal correlation with the DHARMa package
library(DHARMa)
testTemporalAutocorrelation(fit, time = 1:1000)
Durbin-Watson test
data: simulationOutput$scaledResiduals ~ 1
DW = 0.25973, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is not 0Remember that a misfit trend in time or space also creates residual autocorrelation, and the line between a large-scale trend and autocorrelation is fluid, so in real data always de-trend first (Section 11.1.1).
Let’s see how we can fit the AR1 model to data. First, with nlme
library(nlme)
fitGLS = gls(obs~1, correlation = corAR1(0.5, form = ~ 1)) # 0.5 is only a start value
summary(fitGLS)Generalized least squares fit by REML
Model: obs ~ 1
Data: NULL
AIC BIC logLik
2827.461 2842.181 -1410.73
Correlation Structure: AR(1)
Formula: ~1
Parameter estimate(s):
Phi
0.8865404
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.1620833 0.2741234 0.5912785 0.5545
Standardized residuals:
Min Q1 Med Q3 Max
-2.72610072 -0.66439352 0.02987808 0.64460056 3.09922868
Residual standard error: 2.142402
Degrees of freedom: 1000 total; 999 residualSecond, with glmmTMB
library(glmmTMB)
time <- factor(1:1000) # time variable
group = factor(rep(1,1000)) # group (for multiple time series)
fitGLMMTMB = glmmTMB(obs ~ ar1(time + 0 | group))
summary(fitGLMMTMB) Family: gaussian ( identity )
Formula: obs ~ ar1(time + 0 | group)
AIC BIC logLik -2*log(L) df.resid
2828.4 2848.1 -1410.2 2820.4 996
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
group time1 4.49382 2.1199 0.89 (ar1)
Residual 0.02047 0.1431
Number of obs: 1000, groups: group, 1
Dispersion estimate for gaussian family (sigma^2): 0.0205
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.1619 0.2740 0.591 0.555If you check the results, you can see that
- Both models correctly estimate the AR1 parameter
- The p-value for the intercept is in both models n.s., as expected
NoteResidual checks conditional on the fitted random effects
glmmTMB fits the AR1 structure as a random effect, with one value per time point. When DHARMa simulates new data to calculate residuals, the question is whether these random effects are re-simulated (unconditional simulations) or kept at their fitted values (conditional simulations, see also Tip 4.3). Since DHARMa 0.5.0, simulateResiduals() simulates conditionally by default. For glmmTMB models, DHARMa does this via glmmTMB::set_simcodes(fit$obj, val = "fix"); for a few covariance structures, this is not yet possible, and DHARMa then falls back to unconditional simulations with a message. With simulateREs = "unconditional", all random effects are re-simulated. To simulate conditionally directly with glmmTMB:
glmmTMB::set_simcodes(fitGLMMTMB$obj, val = "fix") # note: this modifies the model object
simulate(fitGLMMTMB, nsim = 10)Let’s compare the two options for our AR1 model:
resU = simulateResiduals(fitGLMMTMB, simulateREs = "unconditional")
testTemporalAutocorrelation(resU, time = 1:1000)
Durbin-Watson test
data: simulationOutput$scaledResiduals ~ 1
DW = 0.24923, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is not 0resC = simulateResiduals(fitGLMMTMB) # conditional (default)
testTemporalAutocorrelation(resC, time = 1:1000)
Durbin-Watson test
data: simulationOutput$scaledResiduals ~ 1
DW = 2.8049, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is not 0The unconditional residuals still contain the full autocorrelation, because the AR1 process is re-simulated, so they cannot tell us whether the AR1 model is adequate. The conditional residuals no longer show the positive autocorrelation, but a negative one (Durbin-Watson statistic > 2): with one random effect per time point, the fitted random effects also absorb part of the independent noise, so the conditional residuals are not independent either. Therefore, use DHARMa’s autocorrelation test mainly on the model without a correlation structure, to decide whether you need one. To check the fitted AR1 model, a cleaner option for normal responses are the normalized residuals of gls, from which the estimated correlation has been removed:
acf(residuals(fitGLS, type = "normalized"))
NoteWhy de-trending is important
As discussed in Section 11.1.1, a misfit trend also creates residual autocorrelation. Let’s simulate data with a nonlinear time trend, but with independent residuals, i.e. without any autocorrelation:
set.seed(123)
time = 1:1000/100
y = time + 2*sin(time/0.4) + rnorm(1000)
data = data.frame(y = y, time = time)
plot(y ~ time, data = data)
If we only remove a linear trend, the residuals are strongly autocorrelated, because the sine wave that we did not model is still in them:
fit1 = lm(y ~ time, data = data)
res1 = simulateResiduals(fit1, plot = TRUE)
testTemporalAutocorrelation(res1, time = data$time)
Durbin-Watson test
data: simulationOutput$scaledResiduals ~ 1
DW = 0.62908, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is not 0The correct reaction is not to add an AR1 structure, but to de-trend properly, e.g. with a spline (Section 2.5.3):
library(mgcv)
fit2 = gam(y ~ s(time), data = data)
res2 = simulateResiduals(fit2, plot = TRUE)
testTemporalAutocorrelation(res2, time = data$time)
Durbin-Watson test
data: simulationOutput$scaledResiduals ~ 1
DW = 1.9629, p-value = 0.5574
alternative hypothesis: true autocorrelation is not 0After de-trending, the residuals look fine and there is no autocorrelation left, so we don’t need a correlation structure.
What would have happened if we had added an AR1 structure to the linear model instead?
fit3 = gls(y ~ time, correlation = corAR1(form = ~ 1), data = data, method = "ML")
summary(fit3)Generalized least squares fit by maximum likelihood
Model: y ~ time
Data: data
AIC BIC logLik
3361.755 3381.386 -1676.877
Correlation Structure: AR(1)
Formula: ~1
Parameter estimate(s):
Phi
0.6570029
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.4899957 0.23786121 2.060007 0.0397
time 0.9050320 0.04112851 22.004980 0.0000
Correlation:
(Intr)
time -0.865
Standardized residuals:
Min Q1 Med Q3 Max
-2.654859670 -0.770146615 -0.004566779 0.758276379 2.686580562
Residual standard error: 1.716295
Degrees of freedom: 1000 total; 998 residualAIC(fit1, fit2, fit3) df AIC
fit1 3.00000 3925.468
fit2 10.99608 2867.488
fit3 4.00000 3361.755The AR1 model estimates a considerable autocorrelation (Phi), although there is none in the data: the correlation structure absorbs the misfit of the trend. In this case, the AR1 model is inappropriate. It interprets a systematic trend as random correlated noise, it fits much worse than the GAM (AIC), and its normalized residuals are still not independent (check acf(residuals(fit3, type = "normalized"))). An AR1 structure is meant for the autocorrelation that remains after the trend has been removed, not as a replacement for de-trending.
11.2.1 Exercise — Hurricanes revisited?
CautionExercise
Look at the hurricane study that we also analyze in the case studies (Section C.1), which is, after all, temporal data. This data set is located in DHARMa.
library(DHARMa)
hurricanes$sMasFem = scale(hurricanes$MasFem)
hurricanes$sMinpressure = scale(hurricanes$Minpressure_Updated_2014)
hurricanes$sNDAM = scale(hurricanes$NDAM)
fit = glmmTMB(alldeaths ~ sMasFem * (sMinpressure + sNDAM),
data = hurricanes, family = nbinom2)
# Residual checks with DHARMa.
res = simulateResiduals(fit)
# Checking for temporal autocorrelation
res2 = recalculateResiduals(res, group = hurricanes$Year)
testTemporalAutocorrelation(res2, time = unique(hurricanes$Year))
Durbin-Watson test
data: simulationOutput$scaledResiduals ~ 1
DW = 2.5518, p-value = 0.04758
alternative hypothesis: true autocorrelation is not 0We see that there is a tiny bit of autocorrelation. Note that this autocorrelation is negative, so if we have a higher residual in one year, it is lower in the next year. Normally, we would expect temporal autocorrelation to be positive, so this is a bit weird. It could be an effect of people being more careful if there were many fatalities in the previous year, but it could also be a statistical fluke.
In any case, as an exercise, we want to address the issue - add an AR1 term to the model!
11.3 Spatial correlation structures
Spatial models work very similarly to the temporal models. This time, we start directly with an example, using a data set with the thickness of coal seams, that we try to predict with a spatial (soil) predictor.
library(EcoData)
plot(thick ~ soil, data = thickness)
Let’s fit a simple LM to this
fit = lm(thick ~ soil, data = thickness)
summary(fit)
Call:
lm(formula = thick ~ soil, data = thickness)
Residuals:
Min 1Q Median 3Q Max
-6.0414 -1.1975 0.0876 1.4836 4.9584
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 31.9420 3.1570 10.118 1.54e-15 ***
soil 2.2552 0.8656 2.605 0.0111 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.278 on 73 degrees of freedom
Multiple R-squared: 0.08508, Adjusted R-squared: 0.07254
F-statistic: 6.788 on 1 and 73 DF, p-value: 0.01111DHARMa checks:
res = simulateResiduals(fit)
testSpatialAutocorrelation(res, x = thickness$north, y = thickness$east)
DHARMa Moran's I test for distance-based autocorrelation
data: res
observed = 0.210870, expected = -0.013514, sd = 0.021940, p-value <
2.2e-16
alternative hypothesis: Distance-based autocorrelationFor spatial data, we often look at spatial variograms, which are similar to an acf but in spatial directions
library(gstat)
tann.dir.vgm = variogram(residuals(fit) ~ 1,
loc =~ east + north, data = thickness,
alpha = c(0, 45, 90, 135))
plot(tann.dir.vgm)
How to read a variogram: the x-axis is the distance between pairs of points, the y-axis the semivariance (half the mean squared difference between the residuals of pairs at that distance); for independent residuals the curve is flat, autocorrelation shows as a low semivariance at short distances that increases with distance and then levels off. A curve that keeps rising without leveling off, or that differs between the four directions, indicates a large-scale trend. Both the DHARMa plots and the variograms are more indicative of a spatial trend. Let’s remove this with a 2d-spline, called a tensor spline:
library(mgcv)
fit1 = gam(thick ~ soil + te(east, north) , data = thickness)
summary(fit1)
Family: gaussian
Link function: identity
Formula:
thick ~ soil + te(east, north)
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 39.68933 0.26498 149.780 <2e-16 ***
soil 0.12363 0.07275 1.699 0.0952 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df F p-value
te(east,north) 21.09 22.77 721.3 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.996 Deviance explained = 99.7%
GCV = 0.033201 Scale est. = 0.022981 n = 75plot(fit1, pages = 0, lwd = 2)
We can check the model again, and the problem is gone (the test is no longer significant)
res = simulateResiduals(fit1)
testSpatialAutocorrelation(res, x = thickness$north, y = thickness$east)
DHARMa Moran's I test for distance-based autocorrelation
data: res
observed = -0.024242, expected = -0.013514, sd = 0.021860, p-value =
0.6236
alternative hypothesis: Distance-based autocorrelationIf we had still seen a signal, we should have fit an autoregressive model. Here it’s not necessary, but just to show you the syntax - first nlme:
fit2 = gls(thick ~ soil , correlation = corExp(form =~ east + north) , data = thickness)
summary(fit2)Generalized least squares fit by REML
Model: thick ~ soil
Data: thickness
AIC BIC logLik
164.3474 173.5092 -78.17368
Correlation Structure: Exponential spatial correlation
Formula: ~east + north
Parameter estimate(s):
range
719.4121
Coefficients:
Value Std.Error t-value p-value
(Intercept) 42.81488 5.314542 8.056176 0.0000
soil 0.02662 0.199737 0.133289 0.8943
Correlation:
(Intr)
soil -0.12
Standardized residuals:
Min Q1 Med Q3 Max
-1.5811122 -0.7276873 -0.5028102 -0.2092991 0.3217326
Residual standard error: 5.573087
Degrees of freedom: 75 total; 73 residualThe range parameter is the distance scale of the exponential decay (correlation = exp(-distance / range)). Note the huge range here (about 700, larger than the study area): without a trend term, the correlation structure absorbs the large-scale trend, and the soil effect changes accordingly. This is only to show the syntax; in a real analysis, de-trend first.
Second, for glmmTMB. Here, we again have to prepare the data first
thickness$pos <- numFactor(thickness$east,
thickness$north)
thickness$group <- factor(rep(1, nrow(thickness)))
fit3 = glmmTMB(thick ~ soil + exp(pos + 0 | group) , data = thickness)The output of summary is a bit chunky, which is why I suppress it here
summary(fit3)If you wonder why there is such a large correlation matrix displayed: both the AR1 and the exp(pos + 0 | group) structure impose a particular correlation structure on the random effects. Per default, glmmTMB shows correlations of random effects if they are estimated. In the case of the AR1 structure, the programmers apparently suppressed this, and just showed the estimate of the AR1 parameter. Here, however, they didn’t implement this feature, so you see the entire correlation structure, which is, admittedly, less helpful and should be changed.
11.3.1 Exercise — Does agriculture influence plant species richness?
TipSolution
?EcoData::plantcounts
plants_sf <- plantcounts
str(plants_sf)'data.frame': 285 obs. of 6 variables:
$ tk : int 65341 65342 65343 65344 65351 65352 65353 65354 65361 65362 ...
$ area : num 33.6 33.6 33.6 33.6 33.6 ...
$ richness: int 767 770 741 756 550 434 433 448 527 505 ...
$ agrarea : num 0.488 0.431 0.484 0.598 0.422 ...
$ lon : num 11.4 11.5 11.4 11.5 11.5 ...
$ lat : num 49.5 49.5 49.4 49.4 49.5 ...plants_sf$sAgrarea <- scale(plants_sf$agrarea)
plants_sf$longitude <- plants_sf$lon
plants_sf$latitude <- plants_sf$lat
library(sf)
plants_sf <- sf::st_as_sf(plants_sf, coords = c('longitude', 'latitude'), crs
= st_crs("+proj=longlat +ellps=bessel
+towgs84=606,23,413,0,0,0,0 +no_defs"))
library(mapview)
mapview(plants_sf["richness"], map.types = "OpenTopoMap")fit <- glmmTMB::glmmTMB(richness ~ sAgrarea + offset(log(area)),
family = nbinom1, data = plants_sf)
summary(fit) Family: nbinom1 ( log )
Formula: richness ~ sAgrarea + offset(log(area))
Data: plants_sf
AIC BIC logLik -2*log(L) df.resid
3348.8 3359.8 -1671.4 3342.8 282
Dispersion parameter for nbinom1 family (): 14.3
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.66825 0.01047 254.79 < 2e-16 ***
sAgrarea -0.03316 0.01021 -3.25 0.00117 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(DHARMa)
res <- simulateResiduals(fit)
plot(res)
testSpatialAutocorrelation(res, x = plants_sf$lon, y = plants_sf$lat)
DHARMa Moran's I test for distance-based autocorrelation
data: res
observed = 0.0958792, expected = -0.0035211, sd = 0.0047788, p-value <
2.2e-16
alternative hypothesis: Distance-based autocorrelationfit2<-mgcv::gam(richness ~ sAgrarea + te(lon, lat) + offset(log(area)),
family = nb, data = plants_sf)
summary(fit2)
Family: Negative Binomial(61.871)
Link function: log
Formula:
richness ~ sAgrarea + te(lon, lat) + offset(log(area))
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.66098 0.00800 332.634 <2e-16 ***
sAgrarea -0.02419 0.01013 -2.389 0.0169 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Approximate significance of smooth terms:
edf Ref.df Chi.sq p-value
te(lon,lat) 18.39 21.48 198 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.401 Deviance explained = 48.4%
-REML = 1624.3 Scale est. = 1 n = 285plot(fit2)
library(mgcViz)
b <- getViz(fit2)
print(plot(b, allTerms = F), pages = 1) # Calls print.plotGam()
#plotRGL(sm(b, 1), residuals = TRUE)
res <- simulateResiduals(fit2)
plot(res)
testSpatialAutocorrelation(res, x = plants_sf$lon, y = plants_sf$lat)
DHARMa Moran's I test for distance-based autocorrelation
data: res
observed = -0.0013879, expected = -0.0035211, sd = 0.0047792, p-value =
0.6553
alternative hypothesis: Distance-based autocorrelation11.4 Phylogenetic structures (PGLS)
Closely related species share evolutionary history, so their trait values are not independent. Under Brownian motion, the covariance of two species is proportional to the length of the branch they share (vcv(tree)); a PGLS is simply a gls with this covariance structure, where the distance in the correlation structure is the phylogenetic distance. We fit with method = "ML" so that models with different fixed effects or correlation structures can be compared by AIC / LRT (Section 4.2.4).
This section is mostly taken from https://lukejharmon.github.io/ilhabela/instruction/2015/07/03/PGLS/. The two datasets associated with this example are in the EcoData package.
Perform analysis:
library(EcoData)
library(ape)
library(geiger)
library(nlme)
library(phytools)
library(DHARMa)
library(effects)To plot the phylogenetic tree, use
plot(anolisTree)Regress species traits
# Check whether names are matching in both files.
name.check(anolisTree, anolisData)$tree_not_data
[1] "ahli" "alayoni" "alfaroi" "aliniger"
[5] "allisoni" "allogus" "altitudinalis" "alumina"
[9] "alutaceus" "angusticeps" "argenteolus" "argillaceus"
[13] "armouri" "bahorucoensis" "baleatus" "baracoae"
[17] "barahonae" "barbatus" "barbouri" "bartschi"
[21] "bremeri" "breslini" "brevirostris" "caudalis"
[25] "centralis" "chamaeleonides" "chlorocyanus" "christophei"
[29] "clivicola" "coelestinus" "confusus" "cooki"
[33] "cristatellus" "cupeyalensis" "cuvieri" "cyanopleurus"
[37] "cybotes" "darlingtoni" "distichus" "dolichocephalus"
[41] "equestris" "etheridgei" "eugenegrahami" "evermanni"
[45] "fowleri" "garmani" "grahami" "guafe"
[49] "guamuhaya" "guazuma" "gundlachi" "haetianus"
[53] "hendersoni" "homolechis" "imias" "inexpectatus"
[57] "insolitus" "isolepis" "jubar" "krugi"
[61] "lineatopus" "longitibialis" "loysiana" "lucius"
[65] "luteogularis" "macilentus" "marcanoi" "marron"
[69] "mestrei" "monticola" "noblei" "occultus"
[73] "olssoni" "opalinus" "ophiolepis" "oporinus"
[77] "paternus" "placidus" "poncensis" "porcatus"
[81] "porcus" "pulchellus" "pumilis" "quadriocellifer"
[85] "reconditus" "ricordii" "rubribarbus" "sagrei"
[89] "semilineatus" "sheplani" "shrevei" "singularis"
[93] "smallwoodi" "strahmi" "stratulus" "valencienni"
[97] "vanidicus" "vermiculatus" "websteri" "whitemani"
$data_not_tree
[1] "1" "10" "100" "11" "12" "13" "14" "15" "16" "17" "18" "19"
[13] "2" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29" "3"
[25] "30" "31" "32" "33" "34" "35" "36" "37" "38" "39" "4" "40"
[37] "41" "42" "43" "44" "45" "46" "47" "48" "49" "5" "50" "51"
[49] "52" "53" "54" "55" "56" "57" "58" "59" "6" "60" "61" "62"
[61] "63" "64" "65" "66" "67" "68" "69" "7" "70" "71" "72" "73"
[73] "74" "75" "76" "77" "78" "79" "8" "80" "81" "82" "83" "84"
[85] "85" "86" "87" "88" "89" "9" "90" "91" "92" "93" "94" "95"
[97] "96" "97" "98" "99" # Plot traits.
plot(anolisData[, c("awesomeness", "hostility")])
fit = lm(hostility ~ awesomeness, data = anolisData)
summary(fit)
Call:
lm(formula = hostility ~ awesomeness, data = anolisData)
Residuals:
Min 1Q Median 3Q Max
-0.7035 -0.3065 -0.0416 0.2440 0.7884
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.10843 0.03953 2.743 0.00724 **
awesomeness -0.88116 0.03658 -24.091 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3807 on 98 degrees of freedom
Multiple R-squared: 0.8555, Adjusted R-squared: 0.8541
F-statistic: 580.4 on 1 and 98 DF, p-value: < 2.2e-16plot(allEffects(fit, partial.residuals = TRUE))
Check for phylogenetic signal in residuals.
# Calculate weight matrix for phylogenetic distance.
w = 1/cophenetic(anolisTree)
diag(w) = 0
Moran.I(residuals(fit), w)$observed
[1] 0.05067625
$expected
[1] -0.01010101
$sd
[1] 0.00970256
$p.value
[1] 3.751199e-10Conclusion: signal in the residuals, a normal lm will not work.
You can also check with DHARMa; this also works for GLMMs.
res = simulateResiduals(fit)
testSpatialAutocorrelation(res, distMat = cophenetic(anolisTree))
DHARMa Moran's I test for distance-based autocorrelation
data: res
observed = 0.0509093, expected = -0.0101010, sd = 0.0097304, p-value =
3.609e-10
alternative hypothesis: Distance-based autocorrelationAn old-school method to deal with the problem are the so-called Phylogenetically Independent Contrasts (PICs) (Felsenstein, 1985). The idea here is to transform your data in a way that an lm is still appropriate. For completeness, I show the method here.
Felsenstein, J. (1985) “Phylogenies and the comparative method,” The American Naturalist, 125(1), pp. 1–15.
# Extract columns.
host = anolisData[, "hostility"]
awe = anolisData[, "awesomeness"]
# Give them names.
names(host) = names(awe) = anolisData$species
# Calculate PICs.
hPic = pic(host, anolisTree)
aPic = pic(awe, anolisTree)
# Make a model.
picModel = lm(hPic ~ aPic - 1)
summary(picModel) # Yes, significant.
Call:
lm(formula = hPic ~ aPic - 1)
Residuals:
Min 1Q Median 3Q Max
-2.1051 -0.4188 0.0103 0.3137 4.9991
Coefficients:
Estimate Std. Error t value Pr(>|t|)
aPic -0.97758 0.04516 -21.65 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8967 on 98 degrees of freedom
Multiple R-squared: 0.827, Adjusted R-squared: 0.8253
F-statistic: 468.6 on 1 and 98 DF, p-value: < 2.2e-16# plot results.
plot(hPic ~ aPic)
abline(a = 0, b = coef(picModel))
Now, new school, with a PGLS
pglsModel = gls(hostility ~ awesomeness,
correlation = corBrownian(phy = anolisTree, form =~ species),
data = anolisData, method = "ML")
summary(pglsModel)Generalized least squares fit by maximum likelihood
Model: hostility ~ awesomeness
Data: anolisData
AIC BIC logLik
190.9775 198.793 -92.48875
Correlation Structure: corBrownian
Formula: ~species
Parameter estimate(s):
numeric(0)
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.1505613 0.26262709 0.573289 0.5678
awesomeness -0.9775847 0.04515861 -21.647804 0.0000
Correlation:
(Intr)
awesomeness -0.042
Standardized residuals:
Min Q1 Med Q3 Max
-0.76019997 -0.39056977 -0.04941607 0.19596725 1.07373699
Residual standard error: 0.8877372
Degrees of freedom: 100 total; 98 residualcoef(pglsModel)(Intercept) awesomeness
0.1505613 -0.9775847 plot(allEffects(pglsModel, partial.residuals = TRUE),
partial.residuals = list(smooth = FALSE))
In the effects plot of a gls, the points are the raw data, which are correlated through the phylogeny. We therefore switch off the smoother, which would ignore this correlation and could deviate from the PGLS line even if the model is fine.
OK, same result (the PIC slope equals the PGLS slope, here -0.978), but PGLS is WAY more flexible than PICs. For example, we can include a discrete predictor:
pglsModel2 = gls(hostility ~ ecomorph,
correlation = corBrownian(phy = anolisTree, form =~ species),
data = anolisData, method = "ML")
summary(pglsModel2)Generalized least squares fit by maximum likelihood
Model: hostility ~ ecomorph
Data: anolisData
AIC BIC logLik
374.9473 395.7887 -179.4737
Correlation Structure: corBrownian
Formula: ~species
Parameter estimate(s):
numeric(0)
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.4843515 0.9852947 0.4915803 0.6242
ecomorphGB -0.6315992 1.0376354 -0.6086909 0.5442
ecomorphT -1.0585278 1.3848069 -0.7643865 0.4466
ecomorphTC -0.8558138 0.8919905 -0.9594427 0.3398
ecomorphTG -0.4085610 1.0433839 -0.3915731 0.6963
ecomorphTW -0.4039460 1.0409616 -0.3880508 0.6989
ecomorphU -0.7021719 0.8534461 -0.8227490 0.4128
Correlation:
(Intr) ecmrGB ecmrpT ecmrTC ecmrTG ecmrTW
ecomorphGB -0.636
ecomorphT -0.463 0.451
ecomorphTC -0.551 0.553 0.432
ecomorphTG -0.641 0.723 0.451 0.540
ecomorphTW -0.602 0.572 0.419 0.551 0.549
ecomorphU -0.648 0.615 0.483 0.553 0.606 0.611
Standardized residuals:
Min Q1 Med Q3 Max
-1.40279955 -0.29332546 -0.02361606 0.24737549 1.40998056
Residual standard error: 2.11863
Degrees of freedom: 100 total; 93 residualanova(pglsModel2)Denom. DF: 93
numDF F-value p-value
(Intercept) 1 0.01986379 0.8882
ecomorph 6 0.23482069 0.9641# We can even include multiple predictors:
pglsModel3 = gls(hostility ~ ecomorph * awesomeness,
correlation = corBrownian(phy = anolisTree, form =~ species),
data = anolisData, method = "ML")
summary(pglsModel3)Generalized least squares fit by maximum likelihood
Model: hostility ~ ecomorph * awesomeness
Data: anolisData
AIC BIC logLik
184.7095 223.7871 -77.35476
Correlation Structure: corBrownian
Formula: ~species
Parameter estimate(s):
numeric(0)
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.4944955 0.3719683 1.329402 0.1872
ecomorphGB -0.4531317 0.3941838 -1.149544 0.2535
ecomorphT 0.0742925 0.5481020 0.135545 0.8925
ecomorphTC -0.4676392 0.3497486 -1.337072 0.1847
ecomorphTG -0.3580159 0.3938436 -0.909031 0.3659
ecomorphTW -0.4587856 0.3932411 -1.166678 0.2466
ecomorphU -0.3378731 0.3315315 -1.019128 0.3110
awesomeness -1.1761120 0.0716546 -16.413639 0.0000
ecomorphGB:awesomeness 0.3437046 0.1686570 2.037891 0.0446
ecomorphT:awesomeness 0.2067933 0.1646899 1.255653 0.2126
ecomorphTC:awesomeness 0.1694683 0.1610109 1.052527 0.2955
ecomorphTG:awesomeness 0.5616740 0.1137243 4.938911 0.0000
ecomorphTW:awesomeness 0.0637364 0.2112562 0.301702 0.7636
ecomorphU:awesomeness 0.0801693 0.1282734 0.624988 0.5336
Correlation:
(Intr) ecmrGB ecmrpT ecmrTC ecmrTG ecmrTW ecmrpU awsmns
ecomorphGB -0.637
ecomorphT -0.453 0.438
ecomorphTC -0.548 0.543 0.410
ecomorphTG -0.644 0.721 0.442 0.533
ecomorphTW -0.607 0.571 0.402 0.547 0.551
ecomorphU -0.644 0.606 0.459 0.548 0.600 0.603
awesomeness -0.088 0.082 0.059 0.081 0.081 0.081 0.105
ecomorphGB:awesomeness 0.054 -0.146 -0.028 -0.021 -0.040 -0.042 -0.058 -0.425
ecomorphT:awesomeness 0.059 -0.056 -0.303 -0.033 -0.062 -0.033 -0.029 -0.435
ecomorphTC:awesomeness 0.092 -0.085 -0.065 -0.286 -0.075 -0.090 -0.104 -0.450
ecomorphTG:awesomeness 0.061 -0.079 -0.037 -0.052 -0.102 -0.054 -0.067 -0.630
ecomorphTW:awesomeness 0.061 -0.023 -0.026 -0.021 -0.035 -0.100 -0.038 -0.339
ecomorphU:awesomeness 0.090 -0.075 -0.080 -0.092 -0.078 -0.068 -0.254 -0.551
ecmGB: ecmrT: ecmTC: ecmTG: ecmTW:
ecomorphGB
ecomorphT
ecomorphTC
ecomorphTG
ecomorphTW
ecomorphU
awesomeness
ecomorphGB:awesomeness
ecomorphT:awesomeness 0.184
ecomorphTC:awesomeness 0.187 0.189
ecomorphTG:awesomeness 0.269 0.273 0.285
ecomorphTW:awesomeness 0.165 0.148 0.146 0.212
ecomorphU:awesomeness 0.221 0.229 0.255 0.347 0.184
Standardized residuals:
Min Q1 Med Q3 Max
-1.22551732 -0.37996923 -0.03423687 0.30758384 1.16666459
Residual standard error: 0.7630593
Degrees of freedom: 100 total; 86 residualanova(pglsModel3)Denom. DF: 86
numDF F-value p-value
(Intercept) 1 0.1416 0.7076
ecomorph 6 1.6740 0.1371
awesomeness 1 549.8314 <.0001
ecomorph:awesomeness 6 4.5226 0.0005The gls function also allows other models of trait evolution than Brownian motion, e.g. an Ornstein-Uhlenbeck process. When trying this, however, I noted that the model does not converge due to a scaling problem.
tempTree = anolisTree
pglsModelLambda = gls(hostility ~ awesomeness,
correlation = corPagel(0.9, phy = tempTree,
fixed = FALSE,
form =~ species),
data = anolisData, method = "ML")In this case, if we insist on an OU process, we would have to search for a more stable numeric implementation, either using a Bayesian approach or the caper package.
11.4.1 Exercise — Does bird song change with altitude?
CautionExercise
The following exercise uses data from the 2017 Latin American Macroevolution Workshop in Córdoba (http://www.phytools.org/Cordoba2017/), which examined the influence of environmental factors on the evolution of song in a group of Asian bird species called “barbets.” The following code cleans the raw data available from the EcoData package:
library(ape)
library(EcoData)
dat = barbetData
tree = barbetTree
dat$species = row.names(dat)
plot(tree)
# dropping species in the phylogeny for which we don't have data
obj<-geiger::name.check(tree,dat)
reducedTree<-drop.tip(tree, obj$tree_not_data)
geiger::name.check(reducedTree,dat)[1] "OK"# Moran.I (see below) matches residuals and distance matrix by position,
# so we sort the data in the order of the tips of the tree
dat = dat[reducedTree$tip.label, ]Task: Check if there is a relationship between altitude at which a species is found and the length of the note in its song, which uses the variables Lnote~Lnalt
TipSolution
dat$sLnalt = scale(dat$Lnalt)
plot(Lnote~Lnalt, data = dat)
fit <- lm(Lnote~ sLnalt, data = dat)
summary(fit)
Call:
lm(formula = Lnote ~ sLnalt, data = dat)
Residuals:
Min 1Q Median 3Q Max
-0.11597 -0.06798 -0.03097 0.01019 0.34676
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.11652 0.01994 5.844 1.91e-06 ***
sLnalt 0.02870 0.02025 1.418 0.166
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1145 on 31 degrees of freedom
Multiple R-squared: 0.06087, Adjusted R-squared: 0.03058
F-statistic: 2.009 on 1 and 31 DF, p-value: 0.1663library(effects)
plot(allEffects(fit,partial.residuals = T))
Bit of a misfit, to get a good fit, after playing around, I logged the response and added a quadratic and a cubic effect - you can probably also find other solutions.
fit <- lm(log(Lnote) ~ sLnalt + I(sLnalt^2) + I(sLnalt^3), data = dat)
summary(fit)
Call:
lm(formula = log(Lnote) ~ sLnalt + I(sLnalt^2) + I(sLnalt^3),
data = dat)
Residuals:
Min 1Q Median 3Q Max
-1.54583 -0.34238 -0.05857 0.41449 1.81049
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.04435 0.19067 -15.967 6.64e-16 ***
sLnalt -0.02132 0.22516 -0.095 0.92520
I(sLnalt^2) 0.76725 0.21638 3.546 0.00135 **
I(sLnalt^3) 0.27479 0.10296 2.669 0.01233 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7305 on 29 degrees of freedom
Multiple R-squared: 0.3376, Adjusted R-squared: 0.2691
F-statistic: 4.928 on 3 and 29 DF, p-value: 0.006911plot(allEffects(fit,partial.residuals = T))
Now, with a more complex polynomial for Lnalt, how do we see if there is an overall effect of Lnalt? Easiest option is to do a LRT:
fit0 = lm(log(Lnote)~ 1, data = dat)
anova(fit0, fit)Analysis of Variance Table
Model 1: log(Lnote) ~ 1
Model 2: log(Lnote) ~ sLnalt + I(sLnalt^2) + I(sLnalt^3)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 32 23.366
2 29 15.477 3 7.8893 4.9276 0.006911 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Check residuals for phylogenetic correlation. Note that Moran.I() matches the residuals to the distance matrix by position. This is why we sorted the data in the order of the tips of the tree above; otherwise, the test silently compares the wrong species.
w = 1/cophenetic(reducedTree)
diag(w) = 0
Moran.I(residuals(fit), w)$observed
[1] 0.03858967
$expected
[1] -0.03125
$sd
[1] 0.03427974
$p.value
[1] 0.04161546The residuals are phylogenetically correlated (p = 0.04), so we should not trust the p-value of the lm. Let’s fit the same model with a PGLS:
fit <- gls(log(Lnote) ~ sLnalt + I(sLnalt^2) + I(sLnalt^3),
correlation = corBrownian(phy = reducedTree,
form =~ species), data = dat,
method = "ML")
summary(fit)Generalized least squares fit by maximum likelihood
Model: log(Lnote) ~ sLnalt + I(sLnalt^2) + I(sLnalt^3)
Data: dat
AIC BIC logLik
70.00218 77.48472 -30.00109
Correlation Structure: corBrownian
Formula: ~species
Parameter estimate(s):
numeric(0)
Coefficients:
Value Std.Error t-value p-value
(Intercept) -2.3159881 0.6851223 -3.380401 0.0021
sLnalt 0.2257658 0.2502888 0.902021 0.3745
I(sLnalt^2) 0.2145252 0.1960182 1.094415 0.2828
I(sLnalt^3) 0.0626423 0.0997766 0.627826 0.5350
Correlation:
(Intr) sLnalt I(L^2)
sLnalt 0.224
I(sLnalt^2) -0.286 -0.565
I(sLnalt^3) -0.230 -0.791 0.861
Standardized residuals:
Min Q1 Med Q3 Max
-1.557262740 -0.758489071 -0.349734902 -0.001682088 1.067069864
Residual standard error: 1.101971
Degrees of freedom: 33 total; 29 residualNone of the three terms is significant on its own, but because we have a linear, a quadratic and a cubic term, the best way to test for an overall effect of Lnalt is again an LRT against the intercept-only model:
fit0 <- gls(log(Lnote) ~ 1,
correlation = corBrownian(phy = reducedTree,
form =~ species), data = dat,
method = "ML")
anova(fit0, fit) Model df AIC BIC logLik Test L.Ratio p-value
fit0 1 2 69.91762 72.91064 -32.95881
fit 2 5 70.00218 77.48472 -30.00109 1 vs 2 5.915445 0.1158After accounting for the phylogeny, the evidence for an effect of altitude on note length is much weaker (p = 0.12, compared to p = 0.007 in the lm). The AIC also prefers the PGLS (70.0, compared to 78.7 for the lm, AIC() of both models). Closely related species have similar songs and live at similar altitudes, so the lm treated them as independent replicates and overstated the evidence.
Addition: what would happen if we do the same with a misspecified model (only a linear effect, response not logged)? Have a look at the p-values of the fitted models. Can you explain what’s going on here?
fit <- lm(Lnote~ sLnalt, data = dat)
summary(fit)
Call:
lm(formula = Lnote ~ sLnalt, data = dat)
Residuals:
Min 1Q Median 3Q Max
-0.11597 -0.06798 -0.03097 0.01019 0.34676
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.11652 0.01994 5.844 1.91e-06 ***
sLnalt 0.02870 0.02025 1.418 0.166
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1145 on 31 degrees of freedom
Multiple R-squared: 0.06087, Adjusted R-squared: 0.03058
F-statistic: 2.009 on 1 and 31 DF, p-value: 0.1663plot(allEffects(fit, partial.residuals = T))
Moran.I(residuals(fit), w)$observed
[1] 0.09520617
$expected
[1] -0.03125
$sd
[1] 0.0334855
$p.value
[1] 0.000159082fit <- gls(Lnote~ sLnalt,
correlation = corBrownian(phy = reducedTree,
form =~ species), data = dat,
method = "ML")
summary(fit)Generalized least squares fit by maximum likelihood
Model: Lnote ~ sLnalt
Data: dat
AIC BIC logLik
-77.67763 -73.1881 41.83881
Correlation Structure: corBrownian
Formula: ~species
Parameter estimate(s):
numeric(0)
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.18020988 0.07194706 2.504757 0.0177
sLnalt 0.03941982 0.01556802 2.532103 0.0166
Correlation:
(Intr)
sLnalt 0.15
Standardized residuals:
Min Q1 Med Q3 Max
-1.5260202 -1.0620664 -0.7067039 -0.4679653 2.1556803
Residual standard error: 0.1249469
Degrees of freedom: 33 total; 31 residualHere, it is the other way round: the PGLS effect estimate is significant (p = 0.017), while the normal lm is not (p = 0.17). The residuals of the lm are strongly phylogenetically correlated (Moran’s I, p < 0.001), and the AIC clearly prefers the PGLS (-77.7, compared to -45.4 for the lm), so here, too, the lm is the model that violates its assumptions. The PGLS changes the slope because it effectively compares differences between related species (as the PICs above) instead of treating all species as independent.
What we learn from these examples is that phylogenetic correlation can make the p-values of an lm too small or too large, so it should not be ignored. Instead of only testing the residuals, you can also estimate the strength of the phylogenetic signal together with the model, e.g. with Pagel’s \(\lambda\) (corPagel(), \(\lambda = 0\) means no phylogenetic signal, \(\lambda = 1\) corresponds to Brownian motion):
fitPagel <- gls(Lnote~ sLnalt,
correlation = corPagel(1, phy = reducedTree,
form =~ species), data = dat,
method = "ML")
coef(fitPagel$modelStruct$corStruct, unconstrained = FALSE) lambda
1.014796 Here, \(\lambda\) is estimated to be about 1, i.e. a strong phylogenetic signal, consistent with the Brownian motion model that we used in the PGLS.
11.5 Take-home messages
Tip
- Only the autocorrelation of the residuals matters, not that of the response or the predictors. Ignoring it makes p-values unreliable, so check the residuals of a model without a correlation structure (
acf, variogram, DHARMa) to see whether you need one. - De-trend first (e.g. with a spline), and add a correlation structure only if autocorrelation remains. A correlation structure fitted on top of a misspecified trend absorbs the misfit and wrongly suggests a random correlation.
- To check a fitted correlation structure, use the normalized residuals of
gls(acf(residuals(fit, type = "normalized"))); DHARMa’s autocorrelation test is mainly useful for the model without correlation structure. - Phylogenetic correlation can make the p-values of an
lmtoo small or too large. Check the residuals for phylogenetic signal and use a PGLS if there is some. - Do not interpret a model that has not converged; take warnings seriously and compare several starting values. Residual correlations between responses are not evidence of causal interactions.
11.6 Multivariate GLMs
In recent years, multivariate GLMs, in particular the multivariate probit model, and latent variable versions thereof, have become popular for the analysis of community data. The keyword here is “joint species distribution models” (jSDMs).
Briefly, what we want a jSDM to do is to fit a vector of responses (could be abundance or presence / absence data) as a function of environmental predictors and a covariance between the responses. The model is thus
\[y \sim f(x) + \Sigma\] where y is a response vector, f(x) is a matrix with effects, and \(\Sigma\) is a covariance matrix that we want to estimate.
You can fit these models directly in lme4 or glmmTMB via implementing an RE per species
lme4::glmer(abundance ~ 0 + species + env:species + (0 + species | site), family = poisson, data = dat)
glmmTMB(abundance ~ 0 + species + env:species + (0 + species | site), family = poisson, data = dat)however, these so-called full-rank models have a lot of degrees of freedom and are slow to compute if the number of species gets large. It has therefore become common to fit rank-reduced versions of these models using a latent-variable reparameterization of the model above. A latent-variable version of this model can be fit in glmmTMB via this syntax
glmmTMB(abundance ~ 0 + species + env:species + rr(0 + species | site, d = 2),
family = poisson, data = dat)The parameter d controls the number of latent variables. With d latent variables, all correlations between the responses are generated by d hidden gradients across the sampling units (here: sites), and each response has one loading per gradient; with small d, correlations close to ±1 therefore occur by construction.
Of course, there are many more specialized packages for fitting latent-variable jSDMs in R right now, including hmsc, gllvm or sjSDM, but I find it nice to place these models within the general topic of correlations, and to realize that we can fit them with the same methods and packages as all other correlation structures. The rr() structure in glmmTMB is described by McGillycuddy et al. (2025).
McGillycuddy, M. et al. (2025) “Parsimoniously fitting large multivariate random effects in glmmTMB,” Journal of Statistical Software, 112(1).
11.6.1 Case study: cover of understory species
If a species has more cover than expected in a plot, do other species have more cover too, or less? Separate models for each species cannot answer this (Section 9.2.4.1). We use the varespec data from the plant cover case study in the GLM chapter (Section 9.7.2). To keep the model estimable with only 24 plots, we keep the 21 species that occur in at least 14 plots; the rarer half of the species is dropped, so all correlations refer to the common species. The cover is absolute (plot totals range from 64 to 126%), so a Dirichlet model does not apply. With 24 plots, this is a small data set: take the analysis as a demonstration, not a template.
library(vegan)
data(varespec)
data(varechem)
summary(rowSums(varespec)) Min. 1st Qu. Median Mean 3rd Qu. Max.
64.11 89.67 100.53 100.74 111.39 125.61 keep <- names(which(colSums(varespec > 0) >= 14))
wide <- data.frame(site = factor(rownames(varespec)), varespec[, keep] / 100,
Al = scale(varechem$Al)[, 1], N = scale(varechem$N)[, 1]) # scaled; plain vectors for reshape()
long <- reshape(wide, direction = "long", varying = keep, v.names = "cover",
timevar = "species", times = keep, idvar = "site")
long$species <- factor(long$species)As a baseline, we fit an ordered beta model in which each species has its own intercept, its own response to aluminum and nitrogen, and its own precision (dispformula = ~ 0 + species, see Section 10.1.2; with a common precision for all species, the residual checks fail clearly). The plot effect (1 | site) accounts for plots with more or less vegetation overall. We use the ordered beta rather than the hurdle beta because the hurdle would double the number of parameters per species, too many for 24 plots.
ctl <- glmmTMBControl(optCtrl = list(iter.max = 2000, eval.max = 2000))
m_base <- glmmTMB(cover ~ 0 + species + species:Al + species:N + (1 | site),
dispformula = ~ 0 + species, family = ordbeta(),
data = long, control = ctl)
plot(simulateResiduals(m_base))qu = 0.75, log(sigma) = -2.4794 : outer Newton did not converge fully.
Now we add two latent variables. The larger number of iterations (ctl) is needed for convergence:
m_rr <- update(m_base, . ~ . + rr(0 + species | site, d = 2))
m_rr$sdr$pdHess[1] TRUEm_rr$fit$message[1] "relative convergence (4)"Do the species co-vary beyond what the predictors and the plot effect (1 | site) explain? We compare the two models with AIC and BIC (Note 8.1):
ll <- sapply(list(m_base, m_rr), logLik)
k <- sapply(list(m_base, m_rr), function(m) attr(logLik(m), "df"))
data.frame(model = c("baseline", "baseline + rr(d = 2)"), df = k,
AIC = -2 * ll + 2 * k,
BIC_obs = -2 * ll + log(nrow(long)) * k,
BIC_plots = -2 * ll + log(nlevels(long$site)) * k) model df AIC BIC_obs BIC_plots
1 baseline 87 -2745.017 -2377.653 -2642.526
2 baseline + rr(d = 2) 128 -2843.759 -2303.270 -2692.968AIC clearly favors the joint model. BIC depends on what we count as \(n\): R uses the 504 observations (21 species × 24 plots), but the latent variables are estimated from only 24 plots. With \(n\) = 24, BIC agrees with AIC. We conclude that the species co-vary beyond the common plot effect. We use two latent variables because they can be displayed; more would improve AIC further, but 24 plots cannot support many more parameters.
TipChecking convergence
Latent-variable models often have several local optima, and optimizers can stop early. Refit the model from random starting values of the variance parameters (theta) and compare the log-likelihoods of the fits that converged:
set.seed(4)
nth <- length(m_rr$obj$env$parList()$theta)
sapply(1:6, function(i) {
m <- update(m_rr, start = list(theta = rnorm(nth, 0, 0.5)))
if (m$sdr$pdHess && m$fit$convergence == 0) as.numeric(logLik(m)) else NA
})[1] NA NA 1549.88 NA NA 1549.88as.numeric(logLik(m_rr))[1] 1549.88Two of six starts reach the same optimum as our model; the others stop without convergence (NA). Rules: count a fit only if pdHess is TRUE and the message reports convergence; if converged fits reach different optima, take the best and be cautious; if none converge, reduce d or the number of species. Never interpret a model that has not converged.
We check the residuals by simulating new latent variables (simulateREs = "unconditional", see the explanation of conditional and unconditional simulations in Section 11.2), which checks the overall fit of the model:
plot(simulateResiduals(m_rr, simulateREs = "unconditional"))
Finally, the estimated correlations. In VarCorr(), site is the plot effect and site.1 the rr() term. For readability, we show only species with a latent standard deviation above 0.5:
S <- VarCorr(m_rr)$cond$site.1
dimnames(S) <- list(keep, keep)
round(sort(sqrt(diag(S)), decreasing = TRUE), 2)Callvulg Stersp Cetreric Cladunci Empenigr Cladrang Cladcris Vaccviti
2.79 1.72 1.32 1.02 0.88 0.65 0.57 0.56
Dicrfusc Cladcocc Cladstel Cladcorn Pinusylv Cladarbu Claddefo Cladfimb
0.48 0.45 0.43 0.39 0.39 0.32 0.32 0.30
Cladgrac Pleuschr Polyjuni Pohlnuta Ptilcili
0.18 0.18 0.14 0.13 0.10 sel <- sqrt(diag(S)) > 0.5
heatmap(cov2cor(S[sel, sel]), symm = TRUE, scale = "none")
These are correlations of the latent plot effects on the logit scale, after removing the plot’s overall cover level; they are not correlations of the observed cover. Read them as groups of species that tend to have more (or less) cover in the same plots, not as single values. The fixed effects are now conditional on the latent variables, so some of them change compared with the baseline model (e.g. for Calluna).
ImportantResidual correlation is not a causal link
A residual correlation between two responses means that they are higher (or lower) in the same units than the predictors explain. This can be caused by unmeasured conditions of the units, by spatial structure (Section 11.3), or by a direct effect of one response on the other; the model cannot tell these apart. For species, this means that co-occurrence is not evidence of interactions (Blanchet, Cazelles and Gravel, 2020).
Blanchet, F.G., Cazelles, K. and Gravel, D. (2020) “Co-occurrence is not evidence of ecological interactions,” Ecology Letters, 23(7), pp. 1050–1063.
NoteCommunity data in practice
In ecology, data on many species are often analyzed with ordination (NMDS or PCoA with PERMANOVA, RDA, CCA), which describes the species composition as a whole but gives no model for single species. Model-based alternatives are mvabund::manyglm, gllvm (which also offers ordered and hurdle beta for cover data, as compared by Korhonen et al. (2024)), Hmsc (Bayesian) and sjSDM. If the hypothesis is about diversity, species richness may be the simpler response.
Korhonen, P. et al. (2024) “A comparison of joint species distribution models for percent cover data,” Methods in Ecology and Evolution, 15(12), pp. 2359–2372.