# Make some random error
error <- rnorm(100)
# Make a time-series starting at 0
ts <- rep(0, 100)
# Make it an AR(1) number
for(i in seq(2, length(ts)))
ts[i] <- ts[i-1] + error[i]
# Estimate the AR(1) coefficients
arima(ts, c(1,0,0))
#>
#> Call:
#> arima(x = ts, order = c(1, 0, 0))
#>
#> Coefficients:
#> ar1 intercept
#> 0.9931 9.4306
#> s.e. 0.0082 7.0824
#>
#> sigma^2 estimated as 0.9251: log likelihood = -140.15, aic = 286.29
# Estimate the MA(1) coefficients
arima(ts, c(0,0,1))
#>
#> Call:
#> arima(x = ts, order = c(0, 0, 1))
#>
#> Coefficients:
#> ma1 intercept
#> 0.8907 9.9899
#> s.e. 0.0330 0.6799
#>
#> sigma^2 estimated as 13.05: log likelihood = -271.13, aic = 548.26
# ... the AIC and log-likelihood is better for AR(1)!
arima(ts, c(2,0,3))
#> Warning in arima(ts, c(2, 0, 3)): possible convergence problem: optim gave code
#> = 1
#>
#> Call:
#> arima(x = ts, order = c(2, 0, 3))
#>
#> Coefficients:
#> Warning in sqrt(diag(x$var.coef)): NaNs produced
#> ar1 ar2 ma1 ma2 ma3 intercept
#> 0.167 0.821 0.8547 0.0686 0.0119 9.5044
#> s.e. NaN NaN NaN 0.1363 0.1042 7.7260
#>
#> sigma^2 estimated as 0.9192: log likelihood = -139.91, aic = 293.83
# ... how to fit an ARM(2,3) model (if you want)6 Space, time, and deep-time
6.1 Time-series analysis (temporal autocorrelation)
Overview
This is the first session in the ‘autocorrelation’ series, in which we’re going to be controlling for the kinds of problems that can appear when you collect non-independent—autocorrelated (correlated with itself)—data. In this first session we’re going to cover time series analysis, which could broadly be defined as any set of data that is collected over time. Things like species’ population estimates, stock prices, and temperature are all the kinds of data that are classically called time series data. Of course, it can be a little difficult (indeed, impossible!) to think of data that isn’t collected at some point in time, so perhaps time series data is better thought of as data where the temporal component is important to the questions you care about. The confusing thing about time series data is that most statisticians are interested in removing the influence of time on your dataset, whereas most biologists are interested in quantifying change over time. Those two objectives are not always as easy to reconcile as you might first think! To begin, we’re going to cover the classic models of time series data—ARIMA models—and then move on to how you can remove the influence of time (make your data ‘stationary’), measure the magnitude of an effect (regress the data), and finally do everything at once (using hierarchical models).
6.1.1 AutoRegressive Moving Average (ARMA) models
AutoRegressive Moving Average (ARMA) models are the fundamental building blocks of all time series analyses. The best way to approach them is to build them up in stages, tackling their two separate components, and then finally putting them all together in a single model.
6.1.2 AutoRegressive models
Imagine you have some time series data on a response variable, which we’ll call \(y\), and how it varies over time (\(t\)). We could model \(y\) at each point in time as a function of its value in the previous timestep like this:
\[ y_t = \mu + a y_{t-1} \]
Where \(y_t\) is our value of \(y\) at time \(t\), \(\mu\) is some overall mean1 that describes \(y\), and \(a\) is some measure of how important \(y_t\) is (it’s a parameter/coefficient in our model). This would be called a first-order autoregressive model: first order because it references only one time-step (the previous time-step), and auto-regressive because we’re essentially regressing our response variable against itself2. This is sometimes written as AR(1) in statistics journals. We could generalize such a model to deal with any amount of auto-regression like this:
\[ y_t = \mu + \Sigma{a_iy_{t-i}} \tag{6.1}\]
Where now we’re summing up over \(i\) orders of autoregression, and so this is the general equation for an AR(i) model. Notice that fitting this model to data would be straightforward: it’s a standard linear model, only now instead of comparing \(y\) to some explanatory variables, we’re now also introducing \(y\) as an explanatory variable of itself. Other than this weirdness, the fitting procedure is exactly the same as before. Because we fit this using least squares3, we would call such a model an example of a Generalized Least Squares model, but there’s no more complexity to it than that. This sort of model would be good for something like stock prices or population size estimates: the previous time-step(s) have an impact on the next time-step, and this model captures that well.
6.1.3 Moving Average models
Rarely, if ever, do our models fit our data perfectly, but so far in this course we haven’t really fit models were we have to directly think about the error associated with our measurements4. That all changes today, when we consider moving average models, the first-order definition of which can be written as:
\[ y_t = \mu + \epsilon_t + a \cdot \epsilon_{t-1} \]
Where the only new terms are \(\epsilon_t\) and \(\epsilon_{t-1}\), which are random errors at times \(t\) and \(t-1\) respectively. What makes this moving average model different from that of an auto-regressive model is its focus on averaging the random error that is introduced into the data. This is importantly different from the autoregressive model above, in that now it’s the underlying random noise that is driving a sort-of moving window of impact on the variable. Mathematicians care a great deal about the difference between moving average and autoregressive models because AR models involve the deterministic part of the response variable, whereas MA models involve the stochastic part. We can generalize this equation, in much the same was as we did before, beyond a first order MA model to consider additional levels like this:
\[ y_t = \mu + \Sigma{a_i\epsilon_{t-i}} \tag{6.2}\]
Where now we’re summing up over \(i\) orders of moving average, and so this is the general equation for an MA(i) model. This sort of model would be good for capturing something like blood concentrations of hormones, where recent events matter for predicting the future but, ultimately, those effects wear off through time. In an AR model, because everything depends on the previous value, the effect of early events propagate through the system; in an MA model, the effect is more limited.
6.1.4 Linking it all together
An ARMA model is the combination of Equation 6.1 and Equation 6.2 into one coherent package. Formally, we can write this as:
\[ y_t = \mu + \sum\limits_{i=1}^p{a_iy_{t-i}} + \sum\limits_{j=1}^q{b_j\epsilon_{t-j}} \]
Where \(p\) is the order of the AR part, \(q\) the order of the MA part, and we have to-rename our coefficients so that we have separate multipliers for the AR (\(a\), indexed by \(i\)) and the MA (\(b\), indexed by \(j\)). I’m sorry to have to write out the Sigma notation in full, because I think it looks confusing, but it’s necessary to keep track of the fact that our ARMA process has both \(p\) and \(q\) orders for each of its two components. Indeed, we could write this as \(ARMA(p,q)\), where \(ARMA(1,0)\) would be the same as \(AR(1)\) and \(ARMA(0,1)\) the same as \(MA(1)\).
To make things a little clearer, let’s simulate an ARMA(1,0) (i.e., an AR(1)) process and then see whether an AR or MA model fits it best.
Here the arima function (more on why the function is called arima and not arma in a moment) fits both kinds of model using generalized least squares. You will rarely, if ever, find yourself fitting ARMA models in this way, but hopefully this shows that there’s no magic to it. If you want to figure out what kind of model best fits your data, simply compare the AICs of different structures, using the code above.
6.1.5 Differencing and Integration
A key assumption of ARMA models is that their underlying data are stationary. Stationarity means that it doesn’t matter where in the time series you are: the underlying process is exactly the same. The problem is, of course, that most time series are not like that whatsoever. Without going into a huge amount of detail, there are weaker forms of stationarity and we can fit modified ARMA models—called Autoregressive integrated moving average models—in the case that we find non-stationary data.
The most common form of non-stationarity is trended or seasonal data. For example, it might be the case that, over time, your response variable is increasing. In that case you can use differencing to detrend the data. First-order differencing is literally subtracting \(y_{t-1}\) from \(y_t\), and second-order differencing is when you do the differencing a second time after the first differencing. All of this subtraction means you end up losing data (one data point per difference), and while it sounds complicated hopefully the R code below makes it clear that it’s not as bad as it sounds.
# difference <- ignore the first - ignore the last
first.diff.ts <- ts [-1] - ts [-length(ts)]
second.diff.ts <- first.diff.ts[-1] - first.diff.ts[-length(first.diff.ts)]
# There is a built-in function for this
identical(diff(ts), first.diff.ts)
#> [1] TRUEIt’s difficult to write out differencing in a general equation that doesn’t look scary, so I hope this operational definition is enough for you. There is also something called seasonal differencing, where you calculate the differences among comparable seasons: January 2017 subtracted from January 2018, for example.
Differencing is so useful that it’s the foundation of ARIMA: the opposite of differencing is integrating, and its order is given the term \(d\), making ARIMA models of the form ARMIMA(p,d,q). Thus the model below is ARIMA(1,1,0)—a first-order autoregressive model with a first-order degree of differencing.
# Add trend
trend.ts <- ts + seq(0,1,length.out=100)
# Fit with differencing
arima(trend.ts, c(1,1,0))
#>
#> Call:
#> arima(x = trend.ts, order = c(1, 1, 0))
#>
#> Coefficients:
#> ar1
#> 0.0172
#> s.e. 0.1007
#>
#> sigma^2 estimated as 0.9307: log likelihood = -136.92, aic = 277.84
# ...which gives a better fit than AR(1)
arima(trend.ts, c(1,0,0))
#>
#> Call:
#> arima(x = trend.ts, order = c(1, 0, 0))
#>
#> Coefficients:
#> ar1 intercept
#> 0.9939 9.8962
#> s.e. 0.0075 7.6238
#>
#> sigma^2 estimated as 0.9293: log likelihood = -140.43, aic = 286.85So this gives us our first test of whether or not a time series dataset has a trend: whether an ARIMA model with a non-zero differencing term is statistically supported according to AIC. Remember that we only cared about models within 4 \(\delta\)AIC units of the top model, so a difference of about 4 units is pretty definitive evidence for something.
6.1.6 Estimating the significance of a time series trend
These ARIMA models are all very well and good, but in practice we need to fit models that account not just for change through time, but also for other explanatory variables. So let’s now talk about how to figure out whether the change in a variable in a dataset is statistically significant or not.
6.1.7 The simplest case—time is the only explanatory variable
While it is possible to test whether something is changing through time by making use of an ARIMA model’s differencing term, I suggest that, in the simplest case, you should regress your response variable against time. If the regression is statistically significant, congratulations! Your effect is significant. End of story.
6.1.8 A tricky case—several continuous explanatory variables
Imagine you are modeling someone’s mood as a function of both the time and the local temperature (some like it hot, some like it cold). What do you do now? The problem is that you must judge the significance of the relationship of both of these variables, while also accounting for the temporal autocorrelation of both the explanatory variables (weather tends to be auto-correlated) and the response variable itself (which is the signal we are trying to measure).
The classical solution is, in actual fact, quite difficult. We can explicitly model the expect auto-correlation among the experimental units by specifying something called a covariance matrix based on an ARIMA model of our choosing. A covariance matrix is a way of specifying how similar (or independent) we expect difference things to be from each other, and in this case it can be specified based on the expectations of our ARIMA model. It’s a way of telling R how much each observation in the time series we should expect to be linked to every other observation, on the basis of the ARIMA equations I’ve already shown you. Luckily, we can fit these quite simply in R, even if the underlying math might seem opaque:
library(nlme)
temperature <- cumsum(rnorm(100))
time <- seq_len(100)
mood <- temperature + 1.5 * time + rnorm(100)
gls <- gls(mood ~ time + temperature, correlation=corAR1(form=~time))You might, at first glance, stare at these results and despair, because they imply that time has no effect on mood even though we know it does. The reason is that we have specified an AR(1) form of covariance, and so the effect of time has been soaked up into that. Luckily, we can read from the output that Phi1 is roughly equal to \(0.8\)—that means that, with each year, there is a correlation of \(0.8\) (out of a maximum of 1!) between one year and the next. Thus even though there is no “significant” effect of year, we have a number we can report that shows there is a positive correlation between mood and year. Indeed, we could even perform a test of the significance of that extra year term, proving our point if we needed to in front of a reviewer (remember that our AR(1) term has soaked up a degree of freedom).
null <- gls(mood ~ time + temperature)
anova(gls, null)
#> Model df AIC BIC logLik Test L.Ratio p-value
#> gls 1 5 307.0889 319.9625 -148.5445
#> null 2 4 305.5756 315.8744 -148.7878 1 vs 2 0.4866202 0.4854Note that, if you wished, you could use this sort of model testing to figure out the “significant” ARIMA structure for your data (check out the corARMA function). I would honestly advise you to go easy on such tests, use diagnostic plots, or, perhaps the best, use a hierarchical model.
6.1.9 The hardest case that turns out to be easy—hierarchical models to the rescue
What would we do if we were monitoring more than one person over the course of our study? In other words, what would we do if we had more than one experimental grouping?
The obvious solution to this would be to fit a hierarchical model with an ARIMA model framework sitting inside it. If you would like to do that, you have the following options:
- Fit a mixed effects model (or Bayesian hierarchical model using
rstan) with year “nested” within each observational group. Something like(year|person)would do the trick. - Use
nlme, notlme4, for fitting mixed effects models. In this case, everything works as you would expect. The way of specifying a mixed effects model innlmeis slightly different, and is something like the below, but otherwise it’s straightforward.
library(nlme)
model <- nlme(response ~ explanatory, random = person ~ 1)- Go to the
lme4GitHub repository, find theflexLambdabranch (https://github.com/lme4/lme4/tree/flexLambda), and use that. While this might sound terrifying, you can install directly from GitHub, and the code is quite stable. This is also a great source for a half-finished Phylogenetic Generalized Linear Mixed Model implementation (that may not be finished, but is much faster than mine!). The problem is the person maintaining this has now left academia, I think, and so it’s unlikely to be finished. - Use
rstanorJAGSto build your own ARIMA model. Again, this sounds terrifying, but isn’t that bad once you start. There is also a wrapper function inrstanarmcalledstan_jm, which is remarkably user friendly for the amount of stuff that it can do… - Use the
mgcvlibrary to fit a Generalized Additive Mixed Effects model where you add a smoothed term for time. I don’t cover such models in this class, but the option is there and the package is very nice indeed.
Alternatively, of course, you could difference your response/explanatory variable and then fit a standard regression. Sorry to leave you on something of a cliff-hanger; I’m afraid this is what the cutting edge of statistics feels like!
6.1.10 A note on plotting
You can use something called a partial auto-correlation function to get a clue as to what kind of ARIMA model you might need to fit to your data. The two plots below will give you a feeling for what’s going on.
# Auto-correlation function (nearly useless)
acf(ts, main="Auto-correlation")
# Partial auto-correlation function (useful)
pacf(ts, main="Partial auto-correlation")
The auto-correlation function plot shows the average correlation of all the parts of the response variable at a given time offset. So, for our data with a strong autoregressive component, you can see a strong initial correlation that tails off. The partial autocorrelation function, however, accounts for the correlation across the lags: it is run across the differences of the appropriate variable. The reason this is useful is for autoregressive models, however many lags are significantly correlated is the amount of differencing you need to remove the autocorrelation. So, for example, the single peak at lag 0 means the data is AR(1).
6.1.11 Exercises
As ever, below is some code to load today’s dataset into R, and do a few modifications to it. You must attempt two exercises from the three below. So, for example, you could attempt both data exercises, or one data exercise and one programming exercise.
# Load in the temperature data like this
data <- read.delim("global-temp.tsv", as.is=TRUE)
# Plot the data out
with(data, plot(mean.temp ~ year, type="l"))
# Create an extra explanatory variable
data$recent <- data$year > 1960
# Load in the kick data like this
data <- read.csv("ramona-kicks.csv", as.is=TRUE)
# Plot the data out
with(data, plot(duration, type="l"))
home.stretch <- c(rep(FALSE,40), rep(TRUE,40))- One of today’s datasets is of global temperature averages for the last hundred-or-so years. The data come from NASA (https://climate.nasa.gov/vital-signs/global-temperature/).
- From the partial autocorrelation function, does it appear as though there is temporal autocorrelation in these data?
- Use
arimato test whether these data are best explained by an ARIMA(0,0,1), ARIMA(0,1,0), ARIMA(1,0,0), or ARIMA(1,1,1) model. - Use GLS regression to test whether there is a difference between the mean global temperature now or in the recent (see code above).
- Does your result mean that it’s warmer, colder, or the same temperature as it has always been right now? Why?
- One of today’s datasets is a time-series of my daughter’s kicks while in the womb5. Each day, my wife recorded how long it took her to feel ten kicks in her womb. We’re going to analyze the data to see if there’s a pattern within it.
- From the partial autocorrelation function, does it appear as though there is temporal autocorrelation in these data?
- Use
arimato test whether these data are best explained by an ARIMA(0,0,1), ARIMA(0,1,0), ARIMA(1,0,0), or ARIMA(1,1,1) model. - Use GLS regression to test whether there is a significant difference between kick-counts in the last 30 days (the variable
home.stretch). - Does the presence or absence of a difference in the length of time it takes for ten kicks indicate anything that might be of clinical use?
- Today’s code challenge is to make it a little clearer to you what differencing means. Fill in the blank(s) to write a function that takes, as input, a time series and the order of differences to be applied.
my.diff <- function(ts, n.diff){
for(i in 1:n.diff){
ts <- ____(____)
}
return(ts)
}6.2 Spatial analysis and autocorrelation
Overview
This is the second session in the ‘autocorrelation’ series, in which we’re going to cover the kinds of problems that can arise when we collect data at different points in space. As with time series analysis, it can be difficult to think of how we could collect data outside of space, so perhaps it’s best to think of these techniques as important for data where your explanatory variables change across space. We will only be covering two-dimensional space in this lecture series, and it is very rare for biologists (or really any scientists) to need to control for spatial autocorrelation in a three-dimensional setting. Typically, if you have variation in height (elevation) in your data, this is fitted as an explanatory variable. Perhaps the only exception to this in a biological setting is marine data, and while I have only rarely seen studies of pelagic data that need to account for all three dimensions, in such cases the extension of these techniques (if not these precise packages) to such data is trivial. To begin with, we’re going to extend the classic concept of correlation to account for spatial autocorrelation, which will give us a way to detect whether spatial autocorrelation exists in our datasets. We’ll then cover two common methods for working with spatial data (autoregressive models and Generalized Least Squares (GLS) models), and then finish up with some hands-on experience working with spatial data.
6.2.1 Detecting spatial autocorrelation
While we have spent a great deal of time in this course talking about the square of the correlation coefficient (\(r^2\); see Equation 2.9), but we have not considered its square-root in much detail. You have encountered the correlation coefficient in previous statistical classes, but let’s review it here as we can detect spatial autocorrelation using an extension of it. The correlation coefficient, \(r\), measures whether or not two vectors of data are correlated. There are many different definitions of it but, importantly, the Pearson’s correlation coefficient we are focusing on today doesn’t care which of them is the response or explanatory variable—they’re all the same to it. It is defined as:
\[ r = \frac{\text{Outliers in x} \cdot \text{Outliers in y}}{\text{Variance in x and y}} = \frac{\sum((x - \bar{x})(y-\bar{y}))}{\sqrt{\sum(x-\bar{x})^2}\sqrt{\sum(y-\bar{y})^2}} \tag{6.3}\]
Where \(x\) and \(y\) are the two variables we are testing for correlation, and \(\bar{x}\) and \(\bar{y}\) are the means of these two variables. The idea of a correlation test is to see whether values that are greater than the mean in one variable are paired with values that are above or below the mean in the other variable (that’s our observation, and our expectation is they’re not), in the context of how much background variation there is in the two variables6. If there were lots of unusually large values in \(x\) that were paired with unusually low values in \(y\), we would tend to be multiplying lots of large positive values with lots of negative values, and would sum these up to get a very large negative value. This would be divided by the overall variation in the two variables, and we would get a negative number that was close to \(-1\). The correlation coefficient can only range between \(-1\) and \(+1\): values closer to \(-1\) mean there is a negative correlation between the values, values closer to \(+1\) mean a positive value, and values closer to \(0\) tend to mean there’s little correlation between the variables.
So how can we use this to measure the correlation of values in the same variable through space? Using something called Moran’s I, which is remarkably similar to a Pearson’s correlation coefficient, and is defined as:
\[ I = \text{How many things} \cdot \frac{\text{Outliers here} \cdot \text{Outliers there}}{\text{Variation}} = \frac{n}{\sum w_{i,j}} \cdot \frac{\sum w_{i,j}(x_i-\bar{x})(x_j-\bar{x})}{\sum(x_i-\bar{x})^2} \]
Where \(x\) is our only variable, \(\bar{x}\) is still the mean of \(x\), \(n\) is our number of data points, \(w\) is a matrix (table) of spatial weights among the points and has dimensions \(i\) (rows) and \(j\) (columns) each corresponding to a particular observation of \(x\) such that it has \(n\) columns and \(n\) rows. The spatial weights, \(w\), are directly proportional to the distances among points: closer points are weighted higher (you can decide how much higher; see GLS below). Thus Moran’s I is just the same, really, as a Pearson’s correlation coefficient, only now instead of looking to see whether outliers in \(x\) and \(y\) are correlated, we’re looking to see whether outliers in \(x\) are in roughly the same area in space. In order to bring that concept of space into this, we need to weight our points according to how close they are, as otherwise we would simply be asking how correlated all points within a variable are with each other, which isn’t useful for our purposes7. There are two (equally fine) way to view the standardization term \(\frac{n}{\sum w_{i,j}}\) in all of this: you can view it as dividing through by \(\frac{\sum w_{i,j}}{n}\), in which case it’s part of the variance component (take the overall variation in the weights matrix given how many points we have), or you can view it as you viewed the degrees of freedom in \(F\)-ratios, in which case we’re (again) controlling for how many observations we have and also where they are in space. Either is fine.
Moran’s I values can be compared with an expected distribution of Moran’s I values if there were no aggregation of similar values across space in exactly the same way as any other test statistics. It’s quite common to generate a null distribution of I values through some kind of simulation procedure, however, because there’s a dependence on the locations of your data points in space that can be tricky to solve analytically. Either way, the end result is the same: you get a distribution of expected Moran’s I values, and you compare your observed value with that to find out what fraction of the distribution is greater than your observation8. Note that you can also have spatial disaggregation: values of Moran’s I that are less than \(0\) and reflect points that are close to one-another in space being different from each other.
6.2.2 Accounting for spatial autocorrelation
Imagine you’ve calculated Moran’s I and it appears to be significantly different from zero: you’ve got spatial autocorrelation of some kind. Or, perhaps, you have a strong a priori reason to suppose you’ve got spatial autocorrelation and you don’t want to test for it because you know it’s there9. What now? There is a huge body of literature devoted to modeling spatial data, and I’m going to present to you two reasonably straightforward methods for dealing with it. If you are interested, I cannot recommend Dormann et al. (2007) in Ecography highly enough. It’s showing its age (it tests its models on a Pentium 4!), but the general concepts involved really haven’t changed that much with the exception of recent focuses on Wavelet and spatial point-pattern analysis10.
6.2.3 Simultaneous AutoRegressive models
Simultaneous AutoRegressive models (SARs) are generalizations of the autoregressive models we met last time in a temporal context. The “simultaneous” component here accounts for the fact that we can regress the data against nearby values of the response variable, the explanatory variable, or both during model-fitting. You pick between those options depending on the kind of autocorrelation problem you think, or test to find, that you have. If you think the underlying variable that is driving the autocorrelation is in the explanatory variable (e.g., temperature is driving mood but temperature is spatially autocorrelated) then you would regress your data against nearby temperature values. If you think it’s in the response variable, (e.g., happy people make others around them happy), then you might regress against nearby response variables. If you’re not sure… Go for both.
6.2.4 Generalized Least Squares
We met Generalized Least Squares (GLS) models last time as a way to flexibly add any kind of ARIMA structure into a model. The good news is they can work in exactly the same way for spatial data. Where this can be useful is if you found the decision about what kind of SAR error structure to fit uncomfortable in the last section: with GLS, you have a bit more flexibility, and can also have autocorrelation in your error component. This is analogous to the decision we made in the time series session about whether we wanted to account for autocorrelation of the response variable (autoregressive models) or the error underlying the response variable (moving average models), although the specifics are a little different.
To work with GLS models, you need to specify the functional form of the autocorrelation you think is driving your data. For example, do you think there’s an exponential (decaying) effect of space, where nearby points matter a lot and far away points matter very little? Then you can specify an exponential model of the form \(e^{-\frac{a}{d}}\), where \(e\) is the exponential function, \(d\) is the distance among points, and \(a\) is an estimated parameter that determines how fast the influence of nearby points decays. This is exactly analogous to the moving average models from last time, only now instead of an order (first, second, third, etc.) determining how many timesteps we “listen” to from earlier in the time series, now we have a continuous \(d\)istance and an amount we listen (the exponential function with its slope determined by \(a\)).
6.2.5 Hands on with spatial data in R
I’m not going to sugar-coat it: the hardest part of spatial analysis in R is the actual R. For lots of very good computational reasons11 spatial analysis in R can be a real pain. So please view what follows as a cookbook: when you come to having to do your own spatial analysis in R, take these instructions and modify them to suit your own purposes.
First off, as always, I’m going to simulate some data. This requires the use of a rather fancy package to generate spatially-autocorrelated errors. I’m happy to give a very cursory introduction to how the package works if you’re interested, but please just treat it as a magic black box unless you’re willing to (1) take a serious class in spatial analysis12 or (2) read the vignette of the package carefully. I’m constrained in this class by how much I can teach in half a semester: using the analysis tools I’m teaching you in a real analysis is fine, but using these simulation tools without understanding them more deeply could lead to Very Bad Things.
library(fields)
#> Loading required package: spam
#> Spam version 2.11-4 (2026-05-28) is loaded.
#> Type 'help( Spam)' or 'demo( spam)' for a short introduction
#> and overview of this package.
#> Help for individual functions is also obtained by adding the
#> suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
#>
#> Attaching package: 'spam'
#> The following objects are masked from 'package:base':
#>
#> backsolve, forwardsolve
#> Loading required package: viridisLite
#> Loading required package: RColorBrewer
#>
#> Try help(fields) to get started.
f.obj <- list(x= seq_len(10), y=seq_len(10))
f.obj <- circulantEmbeddingSetup(f.obj)
temperature <- circulantEmbedding(f.obj)
error <- circulantEmbedding(f.obj)
data <- data.frame(
error=as.numeric(error), temperature=as.numeric(temperature),
y=as.numeric(as.numeric(row(error))), x=as.numeric(col(error))
)
data$mood <- with(data, 3 + temperature*.5 + error)We now have a data.frame, called data, that contains some information about people’s average moods through space as a function of temperature. In our simulated data people are strange, and seem to prefer a warmer environment. OK, so how do we go about telling R that we have spatial data, and then how do we use Moran’s I to test whether there’s spatial autocorrelation in this dataset?
# Load useful packages
library(sp)
library(spdep)
#> Loading required package: spData
#> To access larger datasets in this package, install the spDataLarge
#> package with: `install.packages('spDataLarge',
#> repos='https://nowosad.github.io/drat/', type='source')`
#> Loading required package: sf
#> Linking to GEOS 3.14.1, GDAL 3.12.2, PROJ 9.7.1; sf_use_s2() is TRUE
# Copy our data
sp.data <- data
# Tell R our new data structure is spatial
coordinates(sp.data) <- ~x+y
# Build a spatial weight matrix
neighbors <- knn2nb(knearneigh(sp.data, k=8))
weights <- nb2listw(neighbors)
moran.test(data$mood, weights)
#>
#> Moran I test under randomisation
#>
#> data: data$mood
#> weights: weights
#>
#> Moran I statistic standard deviate = 4.7771, p-value = 8.893e-07
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic Expectation Variance
#> 0.214692701 -0.010101010 0.002214327
# For good measure, check our residuals for autocorrelation
naive.model <- lm(mood ~ temperature, data=data)
moran.test(residuals(naive.model), weights)
#>
#> Moran I test under randomisation
#>
#> data: residuals(naive.model)
#> weights: weights
#>
#> Moran I statistic standard deviate = 3.44, p-value = 0.0002909
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic Expectation Variance
#> 0.151813170 -0.010101010 0.002215417
# ...and this has the consequence that our coefficients are wrong
summary(naive.model)
#>
#> Call:
#> lm(formula = mood ~ temperature, data = data)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.64286 -0.76471 -0.01161 0.56784 2.29128
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 3.10888 0.10096 30.793 < 2e-16 ***
#> temperature 0.39558 0.09467 4.179 6.38e-05 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 0.9985 on 98 degrees of freedom
#> Multiple R-squared: 0.1512, Adjusted R-squared: 0.1426
#> F-statistic: 17.46 on 1 and 98 DF, p-value: 6.378e-05
# ... the slope should be 0.5 and the intercept 3A few things about what R is doing “behind the scenes” here. By giving (one) of our datasets coordinates, we’re implicitly telling R to create a spatial data object for us. This is, for reasons I alluded to above, a quite tricky thing to do, and if you have a very large dataset then make sure you have a very large computer. Second, we build our weights matrix using (in this case) the eight nearest neighbors of each point. R can do a lot with that, and it might surprise you that this is enough to detect autocorrelation: you can use more neighbors if you want, but remember that the importance of space tails off with distance and setting this number too high can take up a lot of memory very quickly.
Most importantly, from the perspective of the statistics, we can see that both our response variable and our residuals in our model show spatial autocorrelation. This is a problem: it means our data, and our inferences from that data, are pseudoreplicated, and so we need to account for this in our model. So let’s fit an SAR model to account for this:
# Load yet more packages
library(spatialreg)
#> Loading required package: Matrix
#>
#> Attaching package: 'Matrix'
#> The following object is masked from 'package:spam':
#>
#> det
#>
#> Attaching package: 'spatialreg'
#> The following objects are masked from 'package:spdep':
#>
#> get.ClusterOption, get.coresOption, get.mcOption,
#> get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
#> set.coresOption, set.mcOption, set.VerboseOption,
#> set.ZeroPolicyOption
# Fit the model
sar.model <- errorsarlm(mood ~ temperature, data=data, weights)
# Check the output (it's correct!)
summary(sar.model)
#>
#> Call:errorsarlm(formula = mood ~ temperature, data = data, listw = weights)
#>
#> Residuals:
#> Min 1Q Median 3Q Max
#> -2.56858 -0.68272 -0.04057 0.45409 2.18510
#>
#> Type: error
#> Coefficients: (asymptotic standard errors)
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 3.09957 0.17704 17.5081 < 2.2e-16
#> temperature 0.35219 0.10147 3.4708 0.0005188
#>
#> Lambda: 0.46998, LR test value: 7.8148, p-value: 0.0051821
#> Asymptotic standard error: 0.14372
#> z-value: 3.2702, p-value: 0.0010748
#> Wald statistic: 10.694, p-value: 0.0010748
#>
#> Log likelihood: -136.8246 for error model
#> ML residual variance (sigma squared): 0.87555, (sigma: 0.93571)
#> Number of observations: 100
#> Number of parameters estimated: 4
#> AIC: 281.65, (AIC for lm: 287.46)
# Check for residual autocorrelation (there's none!)
moran.test(residuals(sar.model), weights)
#>
#> Moran I test under randomisation
#>
#> data: residuals(sar.model)
#> weights: weights
#>
#> Moran I statistic standard deviate = 0.48946, p-value = 0.3123
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic Expectation Variance
#> 0.012930766 -0.010101010 0.002214195Voilà! We’ve managed to get the correct parameter estimates out of our model, and removed the residual autocorrelation from our model, all with a SAR model. This particular kind of model only accounts for autocorrelation in our errors, but this was still sufficient to get accurate parameter estimates despite me simulating my explanatory variable and my errors to have spatial autocorrelation. So the take-home from this is not to worry too much if you don’t know what kind of SAR model to fit: one with autocorrelation in the errors is normally sufficient. Let’s move onto a GLS model to see how we fare there…
# Load the package
library(nlme)
# Fit the model (note we're using sp.data now)
gls <- gls(mood ~ temperature, data=sp.data, corr=corExp())
# Check the output - it's correct!
# - p-values are rounded so "0" means "very small"
summary(gls)
#> Generalized least squares fit by REML
#> Model: mood ~ temperature
#> Data: sp.data
#> AIC BIC logLik
#> 293.7738 304.1137 -142.8869
#>
#> Correlation Structure: Exponential spatial correlation
#> Formula: ~1
#> Parameter estimate(s):
#> range
#> 0.4722517
#>
#> Coefficients:
#> Value Std.Error t-value p-value
#> (Intercept) 3.1116643 0.11380528 27.342003 0e+00
#> temperature 0.3867136 0.09876006 3.915688 2e-04
#>
#> Correlation:
#> (Intr)
#> temperature -0.135
#>
#> Standardized residuals:
#> Min Q1 Med Q3 Max
#> -2.65118802 -0.76370733 -0.01176372 0.56295312 2.27901877
#>
#> Residual standard error: 1.000554
#> Degrees of freedom: 100 total; 98 residual
# Check for residual autocorrelation (there's none!)
moran.test(residuals(gls), weights)
#>
#> Moran I test under randomisation
#>
#> data: residuals(gls)
#> weights: weights
#>
#> Moran I statistic standard deviate = 3.4607, p-value = 0.0002694
#> alternative hypothesis: greater
#> sample estimates:
#> Moran I statistic Expectation Variance
#> 0.152793167 -0.010101010 0.002215512The nice thing about the GLS is we have a lot of flexibility about the form of the autocorrelation we think is going on—there are many options beyond corExp for exponential decay. The downside of this flexibility is it can be mind-numbingly slow to fit a model. Even in this reasonably small test-case, it took a very long time. There are other packages that offer spatial GLS regression, and many of those are much faster, but there is an advantage to having nlme::gls in your back-pocket: it can fit mixed effects models. Such models can be fit in Bayesian land, and there are plenty of tutorials for stan and BUGS out there13, but perhaps we can talk about those some other time.
6.2.6 Coda: space isn’t always time, and a statistical fix isn’t always the right answer
I can’t leave you without mentioning something that crops up a lot in my field: the modeling of space as if it were evolutionary time. A long time ago, in a galaxy far, far away, someone thought it would be a very good idea to account for correlation through time on a phylogenetic tree as if it were correlation through space [see Diniz‐Filho et al. (1998) in Evolution]. It was a very good idea, and under some circumstances it can be useful14, but in general it is absolutely not a good idea. This was definitively shown quite some time ago [see Freckleton et al. (2011) in The American Naturalist], but every now and again someone chances upon the idea of using some sort of spatial method in a new context. If you are going to use spatial methods for something that isn’t space, do make sure they’re mathematically simple to describe15, or else you will end up over-engineering something.
Equally, and I might add that this applies to everything I have taught you in this course, do remember that a statistical fix is no substitute for biological insight. There are very few, if any, cases where having all of the important variables that drive a system in a model will go wrong in a straight-forward multiple regression of the type you learned in lecture 4. Everything else, from hierarchical models to spatial autocorrelation, is only needed if you don’t have all the right information to hand. A good example of this is the Mantel test in spatial statistics, which I have deliberately not taught you because it can so easily lead to problems in spatial analyses [see Legendre et al. (2015) in Methods in Ecology & Evolution—the answer to the title of their paper is “no”]. Anything that “cancels out” or “controls for” something inevitably means that there is some process that is being ignored. Modeling spatial autocorrelation with a specific autoregressive or GLS function is not “canceling out” the autocorrelation, it’s measuring its impact. An even better thing to do, of course, would be to find out what is driving that autocorrelation and measure it.
6.2.7 Exercises
As ever, below is some code to load today’s dataset into R, and do a few modifications to it. This dataset is a little different, in that it’s loaded into your R session by the code you’re going to execute below. This code is very useful, not just because the dataset is useful, but also because it shows you how to extract information about particular points on the Earth’s surface (in this case random ones). I can assure you that there will come a day when you will be very grateful for remembering this exercise…
library(terra)
#> terra 1.9.34
#>
#> Attaching package: 'terra'
#> The following object is masked from 'package:fields':
#>
#> describe
library(geodata)
#>
#> Attaching package: 'geodata'
#> The following object is masked from 'package:fields':
#>
#> world
r <- worldclim_global(var = "bio", res = 10, path = tempdir())
#> Cached as: /tmp/RtmpnoARbu/climate/wc2.1_10m//wc2.1_10m_bio.zip
points <- expand.grid(lat = seq(0, 50), long = seq(-50, -150))
v.points <- vect(points, geom = c("long", "lat"), crs = crs(r))
data <- extract(r, v.points, ID = FALSE)
names(data) <- c("temp.mean", "diurnal.range", "isothermality",
"temp.season", "max.temp", "min.temp", "temp.range", "temp.wettest",
"temp.driest", "temp.mean.warmest", "temp.mean.coldest", "precip",
"precip.wettest", "precip.driest", "precip.season", "precip.wettest",
"precip.driest", "precip.warmest", "precip.coldest")
data$lat <- points$lat
data$long <- points$long
data <- na.omit(data)- Today’s dataset is of global weather/climate conditions in roughly the present day. The code above shows you how to download that, and then extract information about these environmental variables at a set of points on the Earth’s surface. You can adapt this code to study some particular points you’re interested in if you wish!
- Regress two of the environmental variables against one-another in a standard regression. What does it show?
- Is there significant spatial autocorrelation in either variable or the model residuals?
- Use either GLS regression or SAR regression to account for potential biases in the data. What does this show?
- Today’s code challenge is all about making it easier for you to grab climate information for any point on the planet’s surface. Do pay attention to the code I have given you above…
grab.cliamte <- function(lat, long){
raster <- getData(____, var="bio", res=10)
points <- data.frame(lat=___, long=___)
coordinates(____) <- ___ ___ ___ ___
return(____(____, points))
}6.3 Deep time (evolutionary) considerations
Overview
This is the second of four sessions in your options series, and is itself the last of two related sessions that form a sort of “what to do with community data” series. If there is a single take-home from the literature that has developed around functional trait and phylogenetic metrics, and their examination in the context of null models, I think it is that concepts matter more than mathematics. In the nearly ten years since the first true eco-phylogenetic metrics were published, a number of criticisms have been raised about various aspects of the metrics. I strongly believe that debate would have been better advanced by distinguishing clearly between the properties of metrics and their performance in measuring a purportedly important property. I would argue that almost every metric has performed perfectly with respect to what it was intended to do, but as each was designed to have different properties they clearly differ. There is no ‘best’ flavor of ice cream, simply flavors that tend to be preferable to many people, and in much the same way there is no ‘right’ eco-phylogenetic metric. There is a parallel, in my opinion, with the literature surrounding conservation prioritization and phylogenetics16
6.3.1 What is functional ecology? Why care about phylogeny?
I think functional ecology is such an important field that I don’t think the field, as it is so-named, really needs a name. Functional ecologists study the ‘functional traits’ of species in order to understand how they function within their environment. They focus on traits that matter (“let the concept of trait be functional!” is their battle-cry), and I am always a bit confused when I meet ecologists who claim not to care about using functional traits. They are a way to generalize across species and ecosystems, and to make sense of the chaos of the millions of species and countless individuals that make up the biosphere. We would be lost without them.
Thus community ecologists, who try to understand the processes by which assemblages of organisms come to be, naturally want to make use of functional traits. One way is to use those traits as predictors in models (see the Fourth Corner analyses we discussed last time), but another is through metrics of the trait structures of assemblages. Ecologists tend to be interested in one of two fundamentally different aspects of ecological structure: the scale and variance of the traits of species in an assemblage. Scale is a measure of the central tendency (often mean) of traits within an assemblage. Typically, community ecologists talk about how environmental conditions ‘filter’ down a ‘pool’ of species that could potentially co-occur within a given assemblage17. Such processes lead to species resembling one-another in terms of trait(s)—perhaps at the top of the hill there are lots of trees, and so the mean height of hill-top assemblages is higher than at the bottom of the hill. Variance is something different: it’s how much species tend to resemble that central tendency, how variable all the species within an assemblage are. Some people prefer the term dispersion for such patterns, particularly in the phylogenetics literature. The problem that has dogged the field, in my opinion, is that so rarely do we distinguish between these two aspects of structure. Environmental filtering, as I have described it, would tend to move the scale of an assemblage in any direction, but would always reduce the variance. This isn’t the case for all processes.
The problem is that it can be time-consuming to collect functional traits data, and we don’t always know what functional traits to collect until after we have finished a study. Through great pain and effort, we have a reasonable set of functional traits that are meaningful for plants, mammals, and birds, but even within those taxa the set of ‘minimal’ traits is far from universally agreed-upon18 and outside those groups… Good luck. Thus the idea came about that maybe species’ evolutionary history—their phylogeny—could be used to act as a proxy for (missing) functional trait information. If species’ evolution was shaped by functional traits (…or shaped it…), then ideally the record of that evolution should be a good substitute for their functional traits. Thus if closely related species resemble one-another, then species co-occur and are closely related must be environmentally filtered on the basis of shared traits, while distantly related species co-occurring must be the result of competition on the basis of those same shared traits. The idea was a good one, and became a victim of its own success: calculating the phylogenetic structure of ecological assemblages became a very popular thing to do, but thinking about why you might want to do it was not so popular. The most important critiques are (1) why bother using phylogeny if you already have the functional trait data19, (2) phylogeny isn’t always a good proxy for functional traits20, and (3) you cannot unambiguously map a (univariate) pattern onto a process. This last point was articulated wonderfully by Mayfield & Levine21 in a seminal paper whose central point seems to have been, in my opinion, entirely missed by the ecological community: environmental filtering and competitive exclusion do not unambiguously onto a single metric of phylogenetic or functional trait structure. This is something we will be addressing only in passing, but I cannot recommend enough the importance of reading both this seminal paper and the work by Peter(s) Chesson and Adler that underlies this point…
6.3.2 Functional trait metrics
We will discuss three fundamental functional trait metrics, but I must emphasize that there are many more: the community-weighted mean, convex hull volume (often called functional richness), and functional dispersion. There are almost more names for metrics than there are metrics themselves, so I would encourage you to use the citations I provide when describing them in manuscripts.
The community-weighted mean (CWM) of sites can either incorporate species’ abundances or not. It is, bluntly, the mean value of a given trait of the species or individuals within an assemblage22. If the weighting is by individual, then it is an abundance-weighted measure, but if it is not, then it ignores the (relative) abundances of species within a given assemblage. Calculating it within R is reasonably trivial; I will use this opportunity to simulate some species, their phylogeny, their functional traits, and then their assemblages. Please don’t worry about tracking all of this code; feel free to treat it as a black-box.
# Simulate phylogeny
library(geiger)
#> Loading required package: ape
#> Registered S3 method overwritten by 'ape':
#> method from
#> plot.mst spdep
#>
#> Attaching package: 'ape'
#> The following objects are masked from 'package:terra':
#>
#> rotate, trans, zoom
#> Loading required package: phytools
#> Loading required package: maps
#>
#> Attaching package: 'phytools'
#> The following object is masked from 'package:terra':
#>
#> rescale
tree <- sim.bdtree(n=25)
# Simulate traits
traits <- sim.char(tree, 1, "BM", nsim=1)[,,1]
# Simulate environmental filtering along a gradient
env <- seq(min(traits),max(traits),length.out=20)
comm <- abs(outer(env, traits, `-`))
comm[comm > 1] <- 1; comm <- 1 - comm
abund <- comm * 5
for(i in seq_len(length(comm)))
comm[i] <- rbinom(1, 1, comm[i])
for(i in seq_len(length(abund)))
abund[i] <- rpois(1, abund[i])
rownames(comm) <- rownames(abund) <- letters[seq_along(env)]
# Sometimes sites with no species are simulated; let's remove those for simplicity
comm <- comm[rowSums(comm)>0, ]
# Calculate CWM (p/a and abundance-weighted) 'by hand'
quick.cwm <- function(mat, trt){
for(i in seq_len(nrow(mat)))
mat[i,] <- mat[i,] * trt
mat[mat==0] <- NA
return(apply(mat, 1, mean, na.rm=TRUE))
}
quick.cwm(comm, traits)
#> a b c d e f
#> -1.58963245 -1.38947443 -1.43825851 -1.35695414 -0.97309209 -0.79906296
#> g h i j k l
#> -0.65548029 -0.62899310 -0.43920736 -0.34315628 -0.14068926 0.02873733
#> m n o p q r
#> 0.52745425 0.52745425 1.64732597 2.12943822 1.80803005 1.74858377
#> s t
#> 1.64732597 1.93901100
quick.cwm(abund, traits)
#> a b c d e f
#> -5.87207923 -5.12468548 -5.27743916 -3.85793918 -3.50209938 -2.68203830
#> g h i j k l
#> -2.86772767 -2.43746936 -1.86751363 -0.56973273 -0.68159761 -0.05054772
#> m n o p q r
#> 1.42429179 1.27403540 2.31684570 5.19966085 6.91071205 4.94197791
#> s t
#> 3.74301159 5.74549833
# Using an R package if you must...
library(FD)
#> Loading required package: ade4
#>
#> Attaching package: 'ade4'
#> The following object is masked from 'package:spdep':
#>
#> mstree
#> Loading required package: geometry
#> Loading required package: vegan
#> Loading required package: permute
#>
#> Attaching package: 'vegan'
#> The following object is masked from 'package:phytools':
#>
#> scores
dbFD(traits, comm, w.abun=FALSE)$CWM # Could use 'comm' if wanted
#> FEVe: Could not be calculated for communities with <3 functionally singular species.
#> FDis: Equals 0 in communities with only one functionally singular species.
#> FRic: Only one continuous trait or dimension in 'x'. FRic was measured as the range, NOT as the convex hull volume.
#> FDiv: Cannot not be computed when 'x' contains one single continuous trait or dimension.
#> Trait
#> a -1.58963245
#> b -1.38947443
#> c -1.43825851
#> d -1.35695414
#> e -0.97309209
#> f -0.79906296
#> g -0.65548029
#> h -0.62899310
#> i -0.43920736
#> j -0.34315628
#> k -0.14068926
#> l 0.02873733
#> m 0.52745425
#> n 0.52745425
#> o 1.64732597
#> p 2.12943822
#> q 1.80803005
#> r 1.74858377
#> s 1.64732597
#> t 1.93901100
dbFD(traits, abund)$CWM
#> FEVe: Could not be calculated for communities with <3 functionally singular species.
#> FRic: Only one continuous trait or dimension in 'x'. FRic was measured as the range, NOT as the convex hull volume.
#> FDiv: Cannot not be computed when 'x' contains one single continuous trait or dimension.
#> Trait
#> a -1.55437391
#> b -1.50726044
#> c -1.40013692
#> d -1.21453641
#> e -1.03340638
#> f -0.86370725
#> g -0.75219086
#> h -0.57352220
#> i -0.47773605
#> j -0.20990153
#> k -0.25965623
#> l -0.02808207
#> m 0.31650929
#> n 1.27403540
#> o 0.99293387
#> p 1.48561738
#> q 1.72767801
#> r 1.64732597
#> s 1.87150580
#> t 1.91516611
# ...as a coding bonus, think about why 'comm' and 'abund' have different
# ...presence/absence CWMsCWMs are good for measuring the average kind of thing you find in a community. Perhaps the tops of hills have tall trees, etc. That doesn’t tell you, however, the variation within a community. This variation is important because it can either (potentially) reveal something about competition23, or it can reveal something about the relative functioning of a community. A community filled with very dissimilar species might be that way because excluding competition has driven similar species to drive each other locally extinct; it could also be the case that such a community provides more ecosystem functions because it contains more variable species. Perhaps it might even be more resilient to change because of the variation it contains: it’s difficult to tell on the basis of these single pieces of data, but that doesn’t stop people trying!
Convex hull volume measures the volume of the space taken up by all the species within a community. It is a powerful metric of the variation in a community, but it is totally ignorant of relative density within the space, and is easily biased by a single outlier species. It’s commonly used in the ecosystem function literature, but these caveats should be kept in mind. A solution to these problems is functional dispersion, which is the mean distance of each species in trait space from the mean (centroid) of all the species in that space. Both these metrics can be calculated with single dimensions of trait data, as I show below, but frankly they’re much more meaningful when you have more than one dimension of trait data. Figure 6.1 shows how they are both calculated from the papers that first defined them (in ecology); notice how only functional dispersion can account for species’ abundances, and it does so by making use of a weighted mean.
# Function to calculate a convex hull volume
conv.hull <- function(comm, trt){
vols <- numeric(nrow(comm))
for(i in seq_along(vols)){
curr.spp <- which(comm[i,] > 0)
if(length(curr.spp) > 2){
vols[i] <- convhulln(trt[curr.spp,], "FA")$vol
} else vols[i] <- NA
}
return(vols)
}
# Calculating convex hull volume with an extra piece of (fake) trait data
conv.hull(comm, cbind(traits,rnorm(25)))
#> [1] 0.6366659 1.4511003 1.2857660 1.3242029 1.7396994 1.9734327 2.2220920
#> [8] 1.7464464 2.2164698 1.7000760 1.1864008 0.7965363 NA NA
#> [15] NA NA 0.7533160 NA NA NA
# ... volume does sort of exist for a single trait (it's the range)
# ... but that's boring so we'll ignore it
# Calculate dispersion
dbFD(traits, comm)$FDis
#> FEVe: Could not be calculated for communities with <3 functionally singular species.
#> FDis: Equals 0 in communities with only one functionally singular species.
#> FRic: Only one continuous trait or dimension in 'x'. FRic was measured as the range, NOT as the convex hull volume.
#> FDiv: Cannot not be computed when 'x' contains one single continuous trait or dimension.
#> a b c d e f g
#> 0.26537124 0.25381972 0.27149465 0.19877344 0.30696828 0.30732234 0.25427547
#> h i j k l m n
#> 0.24900786 0.23938293 0.23266182 0.33006708 0.31200425 0.00000000 0.00000000
#> o p q r s t
#> 0.09502243 0.00000000 0.20107743 0.00000000 0.09502243 0.17870089
# ... notice is warns us that dispersion makes no sense when there's only
# ... one species because that species is the mean (centroid) of the community
6.3.3 Eco-phylogenetic structure
As I have mentioned, it is possible to use species’ evolutionary history as a proxy for functional trait distances. This is an area I have spent some time studying, and so I follow my own review24 in describing these metrics. I would strongly encourage you, if you are interested, to read Caroline Tucker’s excellent review of this topic25 if you are interested in it. That review, frankly, is much more popular than mine and much better, although obviously they both cover much the same material.
Figure 6.2 shows what I consider to be the four kinds of eco-phylogenetic metrics on the basis of the data used to calculate them: shape, evenness, dispersion, and dissimilarity. Shape metrics contain only information about the species present in an assemblage, while evenness metrics incorporate abundance information. One describes the shape of the phylogeny of the species in an assemblage, the other how evenly the individuals within a species are distributed across them. There are many different kinds of these metrics: perhaps the two most popular are Mean Phylogenetic Distance (MPD; also sometimes called the Mean Pairwise Distance) and the Mean Nearest Taxon Distance (MNTD). MPD is the mean of the distances between all the species on a phylogeny, and MNTD is the mean of the distances between each species and its single closest relative. Another popular metric is Faith’s PD, which is the sum of the branch lengths on a phylogeny. This often correlates so strongly with species richness that practicing eco-phylogeneticists rarely use it, but it is a vital metric to know about26. Calculating these metrics is trivial in my package pez.
# Create a 'comparative community' object to match all your data together
# - this can also hold trait and environmental information, but let's ignore that
library(pez)
c.data <- comparative.comm(tree, comm)
abun.c.data <- comparative.comm(tree, abund)
# Calculate specific metrics
.mpd(c.data)
.pd(abun.c.data) # See my footnote for what the extra column is
# Calculate lots of metrics
pez.shape(c.data)
pez.evenness(abun.c.data)Dispersion metrics are the only metrics that consider context. They ask whether the species within an assemblage are more, or less, closely related to one-another than some other set of species—the source pool of species that could be in a region. This kind of question has dominated community ecology for decades—it is one thing to say that an assemblage has low functional dispersion, and quite another to say whether that dispersion is lower than you would expect by chance. The definition of the source pool is one of the most critical components of the entire enterprise. The most common dispersion metrics are \(SES_{MPD}\) and \(SES_{MNTD}\), so-named because they are the Standard Effect Sizes of MPD and MNTD respectively. They involve comparing each metric with a series of randomizations, where you randomly permute either your real data, or some other set of simulated assemblages, lots of times and calculate the MPD and MNTD values under those randomizations. An SES is a test statistic: the observation minus the mean under permutation (the expectation), divided by the standard deviation of the permutations (the variation). Thus, like all test statistics, they’re sensitive to the variation in the simulations as well as the magnitude of the departure from what we would expect by chance. I go over this in some detail in my other statistical espresso class, but remember—\(t = \frac{o-e}{v}\)—and thus it is possible to get a large test statistic value by being very certain (low \(v\)) or having a big effect size (\(o-e\)). I have argued that most eco-phylogenetic studies, by ignoring this, may be flawed27.
Finally, dissimilarity compares the phylogenetic structure of each assemblage with each other assemblage, to see how similar they are. The most common metric of this is UniFrac, which sums the branch lengths two assemblages have in common and divides it by the total branch length of all the species within two assemblages. By repeating this process across all the assemblages in an ecosystem, it is possible to estimate site-level distances. Indeed, this can even be used as the input for a hierarchical cluster analysis!
# Dispersion
.ses.mpd(c.data)
pez.dispersion(abun.c.data)
# Hierarchical cluster analysis of UniFrac
plot(hclust(.unifrac(c.data)))
6.3.4 Exercises
- Below is code to load an example dataset shipped with
pez, which describes some invertebrate data. If you’re interested, read the help file entry (?laja).
library(pez)
data(laja)
# You may need to remove species that are missing from every site
river.sites <- river.sites[,colSums(river.sites)>0]
c.data <- comparative.comm(invert.tree, river.sites, invert.traits)
#> Warning in comparative.comm(invert.tree, river.sites, invert.traits): Mismatch
#> between phylogeny and other data, dropping 2 tips
#> Warning in comparative.comm(invert.tree, river.sites, invert.traits): Mismatch
#> between traits and other data, dropping 2 columns- Calculate some functional diversity metrics for this data. Pick whichever traits you want, but make sure you calculate at least one example of each diversity metric.
- Describe, briefly, what you think these diversity metrics tell you.
- Calculate one phylogenetic shape, one phylogenetic evenness, and one phylogenetic dispersion metric for each of these sites.
- Describe, briefly, what you think these phylogenetic diversity metrics tell you.
or equation!↩︎
‘auto’ is Ancient Greek ‘himself/herself/itself’↩︎
There’s nothing in principle to stop us using this formulation for a response variable that’s binary or count data, I should add↩︎
In our GLMs we did think about our error distributions, but we didn’t really have to work directly with terms for error that much.↩︎
Recording how long it takes your child to kick ten times is one way to track their health while developing, and is an absolutely fascinating way to spend an evening.↩︎
Surprise! It’s another test statistic!↩︎
Indeed, it would just be the variance of \(x\)↩︎
If you have ever taken a course in bootstrapping, you might be a bit concerned that I’ve glossed over quite a lot of details about how to compare an observed value to a bootstrapped distribution. Yet you might also note that, from the definitions of frequentist distributions I gave you earlier in the class, frequentists treat probability as a set of observations in the long-run and so there’s no difference really between bootstrapping and analytical expectations. If you haven’t taken a class in bootstrapping, move on! Nothing to see here…↩︎
Good for you! You have been paying attention!↩︎
Which really, truly are beyond the scope of this course, but I can point you in the direction of good introductions if you are interested.↩︎
Most of which involve the fact that you have to work with spatial distance matrices, because matrices tend to take up lots of memory and are slow to manipulate↩︎
…and I do mean a serious class in spatial analysis↩︎
JAGS is not so spatial-friendly, sorry!↩︎
See, for example, a paper I thought I would hate being involved with but that actually changed my mind on a lot of these eigenvector approaches; Morales‐Castilla et al. (2017) in Global Ecology and Biogeography.↩︎
A single equation, with a handful of parameters, for the correlation structure.↩︎
See Isaac & Pearse (2018) in Phylogenetic Diversity (pp. 27–39). Perhaps this is just the way (applied) phylogeneticists tend to be about things…↩︎
They almost invariably then say that there is no such thing as community ‘filtering’ processes, but given every single community ecologist I have ever met uses the term at some point I’m going to carry on the tradition here.↩︎
If you enjoy watching other people’s arguments, see Mikula et al. (2018) Evolution 72(12): 2832-2835.↩︎
The ‘phylogenetic middleman’ problem that many claim to have thought of first (and may well have done) but only Swenson (2013; Ecography 36(3): 264–276) had the sense to write up and name…↩︎
True (Mazel et al. (2018) Nature communications 9(1): 2888) but it is more often than not.↩︎
2010; Ecology letters 13(9): 1085–1093.↩︎
I prefer to use the term ‘assemblage’ to describe a taxonomically, spatially, and temporally defined unit of ecological measurement, rather than ‘community’. This is because the term ‘community’ means many things to many people (books have literally been written about it), and thus ‘assemblage’ is much less ambiguous and controversial a term. Yet the metric is called CWM, and so that is the term I will use.↩︎
See my overview for why this is risky↩︎
Pearse et al. (2014) in Modern phylogenetic comparative methods and their application in evolutionary biology (pp. 451–464).↩︎
Tucker et al. (2017)Biological Reviews 92(2): 698–715.↩︎
I added a ‘new’ metric, which was Faith’s PD corrected for overall species richness, in the first version of my package
pez. I am convinced I read about it somewhere, but I have no idea where, and so if you find where I first read it I will buy you a Gatorade.↩︎Pearse et al. (2013) Ecology 94 (12): 2861–2872)↩︎