11 Correlation structures
11.1 General Idea
Except for the random effects, we have so far assumed that residual errors are independent. However, that must not always be 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 have similar values. 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. Apart from random effects, where distance is expressed by group, common examples of continuous distance relationships include:
- Random effects (distance = group)
- Spatial distance.
- Temporal distance.
- Phylogenetic distance.
Here a visualization from Roberts et al., 2016 (reproduced as OA, copyright: the authors).
TipThe “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 + (1-a) \cdot \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 (
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.
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 spatially correlated misfit, i.e. there is a trend in the given space (e.g. time, space, 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 because 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 the box above 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
To de-trend, you can just use standard regression terms or splines on time or space. For the rest of this chapter, we will concentrate on how to specify “real” correlation structures. However, in the case studies, you should always de-trend first.
To account for “real” autocorrelation of residuals, similar as for the variance modelling, 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 pages 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. To keep things simple and concentrate on the principles, however, we will stick with the packages you already know.
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 a weighted mean of the last data point and a residual normal distribution in the form
\[x_{t+1} = a \cdot x_t + (1-a) \cdot \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] + (1-a) * 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
-0.58485 -0.14315 0.00559 0.13729 0.66317
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.017019 0.006745 2.523 0.0118 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2133 on 999 degrees of freedomNote that the estimate of the intercept is significant, although we started the simulation at zero. 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, which is estimated by fitting autoregressive models of successively higher orders and checking their residuals.
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)This is DHARMa 0.5.0. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
Note that the default setting in simulateResiduals was changed to conditional simulations since version 0.5.0. This is likely to change residual calculations for all hierarchical (in particular random effect) models. If you want to switch back to the old package version defaults, please use the argument simulateREs = "user-specified" in simulateResiduals(). For more details, see ?simulateREsiduals.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 0
Note
Remember: in general, for spatial / temporal data, there are two processes that can created residual autocorreation:
- There is a spatial misfit trend in time / space, which creates a correlation in space / time.
- There truly is a spatial correlation, after accounting for the trend.
Unfortunately, the distinction between a larger trend and a correlation is quite fluid. Nevertheless, one should always first check for and remove the trend, typically by including time/space as a predictor, potentially in a flexible way (GAMs come in handy). After this is done, we can fit a model with a temporally/spatially correlated error.
Let’s see how we can fit the AR1 model to data. First, with nlme
library(nlme)
Attaching package: 'nlme'The following object is masked from 'package:DHARMa':
getDatafitGLS = gls(obs~1, corr = corAR1(0.771, form = ~ 1))
summary(fitGLS)Generalized least squares fit by REML
Model: obs ~ 1
Data: NULL
AIC BIC logLik
-1773.104 -1758.384 889.552
Correlation Structure: AR(1)
Formula: ~1
Parameter estimate(s):
Phi
0.8865404
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.01620833 0.02741235 0.5912784 0.5545
Standardized residuals:
Min Q1 Med Q3 Max
-2.72610052 -0.66439347 0.02987808 0.64460052 3.09922846
Residual standard error: 0.2142402
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))Warning in glmmTMB(obs ~ ar1(time + 0 | group)): use of the 'data' argument is
recommendedsummary(fitGLMMTMB) Family: gaussian ( identity )
Formula: obs ~ ar1(time + 0 | group)
AIC BIC logLik -2*log(L) df.resid
-1776.7 -1757.1 892.4 -1784.7 996
Random effects:
Conditional model:
Groups Name Variance Std.Dev. Corr
group time1 0.0449381 0.21199 0.89 (ar1)
Residual 0.0002047 0.01431
Number of obs: 1000, groups: group, 1
Dispersion estimate for gaussian family (sigma^2): 0.000205
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.01619 0.02740 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
NoteTrend and autocorrelation with glmmTMB
As I mentioned earlier, first detrend and then add correlation structure if there is autocorrelation. After both steps we should no longer see any pattern in the conditional residuals. Unfortunately, checking the conditional residuals is a bit complicated because glmmTMB does not support conditional simulations, while lme4 does, but it does not support correlation structures. However, there is a workaround.
Let’s start with a small simulation with a time trend and autocorrelation:
time = 1:1000/100
y = time +2*(sin(time/0.4)) + rnorm(1000)
data =
data.frame(y = y, time = time, timeF = as.factor(1:1000), group = as.factor(1))
plot(y ~time)
- Detrend
fit1 = glmmTMB(y~time, data = data)
res = simulateResiduals(fit1, plot = TRUE)
testTemporalAutocorrelation(res, time = data$time)
Durbin-Watson test
data: simulationOutput$scaledResiduals ~ 1
DW = 0.66181, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is not 0Test for temporal autocorrelation is siginificant -> add autoregressive structure
- Add autoregressive structure
fit2 = glmmTMB(y~time + ar1(0+timeF|group), data = data)
res = simulateResiduals(fit2, plot = TRUE)
The residual plot did not change because glmmTMB:::simulate.glmmTMB does not generate conditional predictions. But we can generate them ourselves:
- Create conditional predictions and simulations
We can create a custom DHARMa object with our own simulations:
pred = predict(fit2, re.form = NULL)
simulations = sapply(1:250, function(i) rnorm(1000, pred, summary(fit2)$sigma))
res = createDHARMa(simulations, data$y, pred)
plot(res)
Voila, the residuals look good now!
11.2.1 Exercise - hurricanes revisited?
CautionExercise
Look at the hurricane study that we used before, which is, after all, temporal data. This data set is located in DHARMa.
library(DHARMa)
fit = glmmTMB(alldeaths ~ scale(MasFem) *
(scale(Minpressure_Updated_2014) + scale(NDAM)),
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 mere 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 similar 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)
Both the DHARMa plots and the variograms are more indicative of a spatial trend. Let’s remove this with a 2d-spine, 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
res = simulateResiduals(fit1)Registered S3 method overwritten by 'mgcViz':
method from
+.gg ggplot2testSpatialAutocorrelation(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 autocorrelationAlmost the same, but simpler:
fit = lm(thick ~ soil + north + I(north^2), data = thickness)If we would have 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 residualSecond, 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$agrarea_scaled <- 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 ~ agrarea_scaled + offset(log(area)),
family = nbinom1, data = plants_sf)
summary(fit) Family: nbinom1 ( log )
Formula: richness ~ agrarea_scaled + 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 ***
agrarea_scaled -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 ~ agrarea_scaled + te(lon, lat),
offset(log(area)), family = nb, data = plants_sf)
summary(fit2)
Family: Negative Binomial(67.736)
Link function: log
Formula:
richness ~ agrarea_scaled + te(lon, lat)
Parametric coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 6.183373 0.004096 1509.72 < 2e-16 ***
agrarea_scaled -0.024366 0.005355 -4.55 5.37e-06 ***
---
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) 22.53 23.76 850.3 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R-sq.(adj) = 0.413 Deviance explained = 50.2%
-REML = 5622.6 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.0030357, expected = -0.0035211, sd = 0.0047800, p-value =
0.9191
alternative hypothesis: Distance-based autocorrelation11.4 Phylogenetic Structures (PGLS)
This 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)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")])
plot(hostility ~ awesomeness, data = anolisData)
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-16abline(fit)
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, using this works also 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, J. (1985) “Phylogenies and the comparative method”. American Naturalist, 125, 1–15.). The idea here is to transform your data in a way that an lm is still appropriate. For completeness, I show the method here.
# Extract columns.
host = anolisData[, "hostility"]
awe = anolisData[, "awesomeness"]
# Give them names.
names(host) = names(awe) = rownames(anolisData)
# Calculate PICs.
hPic = pic(host, anolisTree)Warning in pic(host, anolisTree): the names of argument 'x' and the tip labels
of the tree did not match: the former were ignored in the analysis.aPic = pic(awe, anolisTree)Warning in pic(awe, anolisTree): the names of argument 'x' and the tip labels
of the tree did not match: the former were ignored in the analysis.# Make a model.
picModel = lm(hPic ~ aPic - 1)
summary(picModel) # Yes, significant.
Call:
lm(formula = hPic ~ aPic - 1)
Residuals:
Min 1Q Median 3Q Max
-1.30230 -0.23485 0.06003 0.34772 0.92222
Coefficients:
Estimate Std. Error t value Pr(>|t|)
aPic -0.91964 0.03887 -23.66 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4263 on 98 degrees of freedom
Multiple R-squared: 0.851, Adjusted R-squared: 0.8495
F-statistic: 559.9 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(hostility ~ awesomeness, data = anolisData)
abline(pglsModel, col = "red")
OK, same result, 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?
11.4.1.1 Exercise
Caution
The following exercise uses data from a study by Corboda et al., 2017, which examined the influence of environmental factors on the evolution of song in an 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"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
plot(Lnote~Lnalt, data = dat)
fit <- lm(Lnote~ scale(Lnalt), data = dat)
summary(fit)
Call:
lm(formula = Lnote ~ scale(Lnalt), 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 ***
scale(Lnalt) 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)Loading required package: carDatalattice theme set by effectsTheme()
See ?effectsTheme for details.plot(allEffects(fit,partial.residuals = T))Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictor scale(Lnalt) is a one-column matrix that was converted to a vector
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) ~ scale(Lnalt) + I(scale(Lnalt)^2) + I(scale(Lnalt)^3), data = dat)
summary(fit)
Call:
lm(formula = log(Lnote) ~ scale(Lnalt) + I(scale(Lnalt)^2) +
I(scale(Lnalt)^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 ***
scale(Lnalt) -0.02132 0.22516 -0.095 0.92520
I(scale(Lnalt)^2) 0.76725 0.21638 3.546 0.00135 **
I(scale(Lnalt)^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))Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors scale(Lnalt), I(scale(Lnalt)^2), I(scale(Lnalt)^3) are one-column
matrices that were converted to vectors
Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors scale(Lnalt), I(scale(Lnalt)^2), I(scale(Lnalt)^3) are one-column
matrices that were converted to vectors
Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictors scale(Lnalt), I(scale(Lnalt)^2), I(scale(Lnalt)^3) are one-column
matrices that were converted to vectors
Now, with a more complex polynomial for Lnalt, how to 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) ~ scale(Lnalt) + I(scale(Lnalt)^2) + I(scale(Lnalt)^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
w = 1/cophenetic(reducedTree)
diag(w) = 0
Moran.I(residuals(fit), w)$observed
[1] -0.06427062
$expected
[1] -0.03125
$sd
[1] 0.03427974
$p.value
[1] 0.3354125Nothing! So we could leave the model as it is. Just for completeness, fit the same comparison with a PGLS, effect remains significant, but p-value a bit larger.
fit <- gls(Lnote~ scale(Lnalt) * I(scale(Lnalt)^2),
correlation = corBrownian(phy = reducedTree,
form =~ species), data = dat,
method = "ML")
summary(fit)Generalized least squares fit by maximum likelihood
Model: Lnote ~ scale(Lnalt) * I(scale(Lnalt)^2)
Data: dat
AIC BIC logLik
-76.95977 -69.47723 43.47988
Correlation Structure: corBrownian
Formula: ~species
Parameter estimate(s):
numeric(0)
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.14433555 0.07391392 1.9527520 0.0606
scale(Lnalt) 0.02481932 0.02700222 0.9191586 0.3656
I(scale(Lnalt)^2) 0.03434214 0.02114728 1.6239510 0.1152
scale(Lnalt):I(scale(Lnalt)^2) 0.01160067 0.01076432 1.0776963 0.2901
Correlation:
(Intr) scl(L) I((L)^
scale(Lnalt) 0.224
I(scale(Lnalt)^2) -0.286 -0.565
scale(Lnalt):I(scale(Lnalt)^2) -0.230 -0.791 0.861
Standardized residuals:
Min Q1 Med Q3 Max
-1.5858745 -0.8709001 -0.7491330 -0.3869937 2.0469866
Residual standard error: 0.1188854
Degrees of freedom: 33 total; 29 residualBtw, did you notice that neither Lnalt has a significant effect - so is there no dependence on Lnalt? Because we have the linear and quadaratic effect, the best way to answer this is to run an ANOVA:
fit0 <- gls(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 -73.47388 -70.48087 38.73694
fit 2 5 -76.95977 -69.47723 43.47988 1 vs 2 9.485883 0.0235which tells us that Lnalt has a significant effect, even though we can’t reject that either the linear or the quadrart term is zero.
Addition: what would happen if we do the same with a misspecified model? Have a look at the p-values of the fitted models. Can you explain what’s going on here?
fit <- lm(Lnote~ scale(Lnalt), data = dat)
summary(fit)
Call:
lm(formula = Lnote ~ scale(Lnalt), 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 ***
scale(Lnalt) 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))Warning in Analyze.model(focal.predictors, mod, xlevels, default.levels, : the
predictor scale(Lnalt) is a one-column matrix that was converted to a vector
w = 1/cophenetic(reducedTree)
diag(w) = 0
Moran.I(residuals(fit), w)$observed
[1] -0.05831654
$expected
[1] -0.03125
$sd
[1] 0.0334855
$p.value
[1] 0.4189142fit <- gls(Lnote~ scale(Lnalt),
correlation = corBrownian(phy = reducedTree,
form =~ species), data = dat,
method = "ML")
summary(fit)Generalized least squares fit by maximum likelihood
Model: Lnote ~ scale(Lnalt)
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
scale(Lnalt) 0.03941982 0.01556802 2.532103 0.0166
Correlation:
(Intr)
scale(Lnalt) 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 residualThe observation is that the PGLS effect estimate is significant while normal lm is not. The reason is probably that the PGLS is re-weighting residuals, and it seems that in this case, the re-weighting is changing the slope. What we learn by this example is that a PGLS can increase significance, and in this case I would argue wrongly so, as we have no indication that there is a phylogenetic signal. I would therefore NOT recommend to blindly fit PGLS, but rather test first if a PGLS is needed, and only then apply.
11.5 Multivariate GLMs
In the 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(abundance ~ 0 + species + env:species + (0+species|site)
glmmTMB(abundance ~ 0 + species + env:species + (0+species|site)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 the 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(Species + 0|id, d = 2))The parameter d controlls the number of latent variables.
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 set this models in the general topic of correlations, and realize that we can fit them with the same methods and packages as for all other correlation structures.