I am examining a temporal dataset of contaminant concentrations from individual specimens over a 100 year time span, comparing mercury (THg) concentrations over time (in years). For some of the years, I have single samples, whereas for others I have over 10 samples in a given year. How best could I control for the differences in sampling effort among years?
This is one of the benefits of using a multilevel or mixed model. By specifying a random intercept for specimen, the model gives greater weight to specimens with more measurements and less weight to those with fewer. In the case of a specimen with a single measurement, the model would likely predict that such a specimen would have a value on the outcome that is essentially equivalent to the mean outcome across all specimens over time.
In lme4, an example of this model with a single predictor (time) would be the following:
m <- lmer(y ~ 1 + x + (1|specimen), data=df)
You can use
ranef(m) to see exactly what the model predicts for each of your specimens. For those specimens with very few observations or just 1 observation, the prediction would be 0, which is understood relative to the
(Intercept) from the summary of your model. A positive prediction means that for that specimen, the model predicts their outcome value to be higher than the intercept and vice-versa for a negative prediction.
You can further plot these using
dotplot(ranef(m,condVar=TRUE)) to see the distribution of predictions.
Based on the comments to the post, I'm adding a bit more information. The model I suggested above,
m, treats the effect of time (
x) on y as linear and as "fixed," meaning the trend is the same across all specimens. This is an assumption and it can easily be relaxed, allowing for unique linear trends for each specimen (i.e., a random slope on
x) and tested as follows:
m1 <- lmer(y ~ 1 + x + (x|specimen), data=df) anova(m, m1) #likelihood ratio test comparing fixed time effect to random time effect
If the $\chi^2$ test is significant, then the random effect model is preferred; the increased complexity of
m1 is a better fit to the data. It may be worth exploring different functional forms of x given the long measurement period here. That can be done with, e.g.:
m2 <- lmer(y ~ 1 + poly(x, 2) + (x|specimen), data=df) #fixed only m3 <- lmer(y ~ 1 + poly(x, 2) + (poly(x, 2)|specimen), data=df) #fixed + random
These models can again be likelihood ratio tested against simpler models. Note that if all you are doing is changing a fixed effect predictor and want to run a likelihood ratio test, you should be using full maximum likelihood (
REML=FALSE in your
lmer model). The
anova command will do this automatically.
Because of the unique nature of your data, it may be worth exploring generalized additive mixed models to allow for the most flexible modeling of the effect of time on the outcome. This has an introduction to GAMs and then some specific information on the mixed version of GAMs.
Here is a description from The Analysis Factor article on Approaches to Repeated Measures Data: Repeated Measures ANOVA, Marginal, and Mixed Models by Karen Grace-MAartin. To quote some comments:
In a marginal model (AKA, the population averaged model), the model equation is written just like any linear model. There is a single response and a single residual. The difference between the marginal model and a linear model is that the residuals are not assumed to be independent with constant variance.
In a marginal model, we can directly estimate the correlations among each individual’s residuals. (We do assume the residuals across different individuals are independent of each other). We can specify that they are equally correlated, as in the RM ANOVA, but we’re not limited to that assumption. Each correlation can be unique, or measurements closer in time can have higher correlations than those farther away. There are a number of common patterns that the residuals tend to take.
And, with respect to the option of employing a Linear Mixed Model:
It too controls for non-independence among the repeated observations for each individual, but it does so in a conceptually different way. Rather than just estimate the correlation among an individual’s repeated observations, it actually adds one or more random effects for Individuals to the model.
The model equation therefore includes extra parameters to include any random effects. They take the form of additional residual terms, each of which has its own variance to be estimated.
This literally means the model is controlling for the effects of individual. The simplest mixed model, the random intercept model, controls for the fact that some individuals always have higher values than others. By controlling for this variation, we’ve taken it out of the original residual.
Individual growth curve models are a specific type of mixed model that uniquely models each individual’s value of the outcome over time. They are particularly useful when the research question is about how covariates affect not only the value of the dependent variable, but its change over time.
The biggest advantage of mixed models is their incredible flexibility. They can handle clustered individuals as well as repeated measures (even in the same model). They can handle crossed random effects, where there are repeated measures not only on an individual, but also on each stimulus.
Time can easily be considered continuous or categorical, and covariates can be measured just once per individual or repeatedly at each observation. Unbalanced data are no problem, and even if some outcomes are missing for some individuals, they won’t be dropped from the model.
The biggest disadvantage of mixed models, at least for someone new to them, is their incredible flexibility. It’s easy to mis-specify a mixed model, and this is a place where a little knowledge is definitely dangerous.
In the case that your data contains individuals with say measurements of heavy metals in their blood over time, but not necessarily the same individuals (given that is a 100 years of data), I suspect that with the simplest mixed model, namely the random intercept model, one may have a model specification issue. For example, in measuring blood levels of say Lead which actually relates to water and air exposure to all individuals over time in a particular locale, the process for controlling the effects of the individual (who happen to change over time), there could be a misstatement in the individual effect (the intercept). The latter may actually be part of the time trend (from say, mining activity and associated pollution over years with a corresponding impact on the general population). In such a case, working with the marginal model may be better.
Bottom line, more clarity on what is being modeled (like pollution related to a particular locale) and is the selected model appropriate and likely could provide the best fit for a time trend (or forecast, which may be of interest).
On the last comment, to quote from a prior thread: "However, the problem is that so far in my experience population-level predictions based on mixed models are significantly worse than predictions based on standard regression models with fixed effects only." Supplied answers seem to suggest, it depends (hence, my call for clarity).