Nesting in mixed effect models

Question

I am interested in how a response varies between time periods, and regions. My question revolves around nesting of random effects, specifically putting site (unique within each region) as a random effect.

Simulated data set for mixed effects modelling

factor1 <- c(rep("a", 100), rep("b", 100))#time period
factor2 <- rep(rep(1:5, 20), 2) #5 unique regions
factor3 <- c(1:100, 1:100)#100 unique sites

response <- c(rnorm(100, 10, 3), rnorm(100, 15, 3))

df <- data.frame(time.period = factor1, region = as.factor(factor2), site = as.factor(factor3), response = response)
df %>% str() #structure of the data

#quick visual
ggplot(df, aes(time.period, response)) + geom_boxplot() + facet_grid(~region)

lm1 <- lm(response ~ time.period*region, data=df)#fixed effects model
anova(lm1)
summary(lm1)

I received comments from a reviewer which stated that the "among site variation should be controlled for as random effects". Specifically, they suggested nesting the random effect of site within the fixed effect of region. They also suggested that including a site effect would potentially remove a region effect.

I'm having trouble conceptually with this, and have the following questions

1) Would you even need to nest site within region? Each site has a unique identifier, so wouldn't nesting be unnecessary?

2) How would including a random site effect "remove a region effect"? Doesn't that imply making region a random effect?

*Note that I actually have multiple observations within each time period (2-4).

Below is the code I've written for implementing the suggested approach:

library(lme4)

lme1 <- lmer(response ~ time.period*region*(1|site), data=df)
anova(lme1)
summary(lme1)

Answer 1

I am interested in how a response varies between time periods, and regions

As they are central to your research question, the main effects of time.period and region, along with their interaction, should be fixed effects in the model formulation.

Since each site is sampled several times, in periods 2-4, the reviewer is correct to point out that among-site variation needs to be accounted for, and one way to do that is by using a mixed effects model with random intercepts for site. So, the proposed model:

lme1 <- lmer(response ~ time.period*region*(1|site), data=df)

does make sense. However there are further considerations. When we have repeated measures over time, not only do we expect there to be correlation within sites, but we also expect there to be autocorrelation. This can be handled implicitly in lmer by fitting random slopes for time. Such a model would be specified as:

lmer(response ~ time.period*region*(1+time.period|site), data=df)

Some other packages in R, such as glmmTMB and mmrm can handle autocorrelation explicitly via the estimation of an AR(1) parameter).

As for the other questions:

1. Would you even need to nest site within region? Each site has a unique identifier, so wouldn't nesting be unnecessary?

Your intuition is correct, and it would not make sense to have site nested within region

2. How would including a random site effect "remove a region effect"? Doesn't that imply making region a random effect?

Including a random site effect doesn't remove the region effect. It separates the variation due to region from the variation due to specific sites within each region. This allows for a more accurate estimate of the true region effect.