CODEWITHSUNDEEP
rstatisticslme4nlme
mfbeuq
mfbeuq

Reputation: 23

GLS / GLM nested design with autocorrelation over time

Still fairly new to GLM and a bit confused about how to establish my model.

About my project:

I sampled the microbiome (and measured a diversity index value = Shannon) from the root system of a sample of 9 trees (=tree1_cat).

In each tree I sampled fine and thick roots (=rootpart) and each tree was sampled four times (=days) over the course of one season. Thus I have a nested design but have to keep time in mind for autocorrelation. Also not all values are present, thus I have a few missing values). So far I have tried and tested the following:

Model <- gls(Shannon ~ tree1_cat/rootpart + tree1_cat + days,  
     na.action = na.omit, data = psL.meta, 
     correlation = corAR1(form =~ 1|days), 
     weights = varIdent(form= ~ 1|days))

Furthermore I've tried to get more insight and used anova(Model) to get the p-values of those factors. Am I allowed to use those p-values? Also I've used emmeans(Model, specs = pairwise ~ rootpart)for pairwise comparisons but since rootpart was entered as nested factor it only gives me the paired interactions.

It all works, but I am not sure, whether this is the right model! Any help would be highly appreciated!

Upvotes: 0

Views: 777

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226547

It would be helpful to know your scientific question, but let's suppose you're interested in differences in Shannon diversity between fine and thick roots and in time trends. A model you could use would be:

library(lmerTest)
lmer(Shannon ~ rootpart*days + (rootpart*days|tree1_cat), data = ...)

The fixed-effect component rootpart*days can be expanded into 1 + rootpart + days + rootpart:days (where 1 signifies the intercept)

  • intercept: SD in fine roots on day 0 (hopefully that's the beginning of the season)
  • rootpart: difference between fine and thick roots on day 0
  • days: change per day in SD in fine roots (slope)
  • rootpart:days difference in slope between thick roots and fine roots

The random-effect component (rootpart*days|tree1_cat) measures how all four of these effects vary across trees, and their correlations (e.g. do trees with a larger-than-average difference between fine and thick roots on day 0 also have a larger-than-average change over time in fine root SD?)

This 'maximal' random effects model is almost certainly too complex for your data; a rough rule of thumb says you should have 10-20 data points per parameter estimated, the fixed-effect model takes 4 parameters. A full model with 4 random effects requires the estimate of a 4×4 covariance matrix, which has (4*5)/2 = 10 parameters all by itself. I might just try (1+days|tree1_cat) (random slopes) or (rootpart|tree_cat) (among-tree difference in fine vs. thick differences), with a bias towards allowing for the variation in the effect that is your primary interest (e.g. if your primary question is about fine vs. thick then go with (rootpart|tree_cat).

I probably wouldn't worry about autocorrelation at all, nor heteroscedasticity by day (your varIdent(~1|days) term) unless those patterns are very strongly evident in the data.

If you want to allow for autocorrelation you'll need to fit the model with nlme::lme or glmmTMB (lmer still doesn't have machinery for autocorrelation models); something like

library(nlme)
lme(Shannon ~ rootpart*days,
    random = ~days|tree1_cat,
    data = ...,
    correlation = corCAR1(form = ~days|tree1_cat)
)

You need to use corCAR1 (continuous-time autoregressive order-1) rather than the more common corAR1 for unevenly sampled data. Be aware that lme is more finicky/worse at dealing with singular models, so you may discover you need to simplify your model before you can actually get this model to run.

Upvotes: 1

Related Questions