Solved – Nested random effects and interaction terms in lme4

lme4-nlmenested datarrandom-effects-model

I have a data set containing various vegetation and geomorphic variables sampled in 3 distances on both sides of 43 drainage ditches (Location). Roughly half of these ditches are occupied by a beaver, the other half is empty. Now I want to run a model with the binomial response variable Status ("beaver == 1" / "beaver == 0")

I'm struggling with the order and layout of the nested and interaction effects using glmer. So far I've got

fit <- glmer(Status ~ BankslopeScaled + Connectivity + 
                      Canal_width + Distance:Food_crops + 
                      Distance:Edible_trees + 
                (1 | Distance/Side/Location), 
              data, family=binomial(link="logit")

but I'm not sure if I still have pseudoreplication in my data or whether I correctly applied the formula in order to estimate the influence of the predictors in every distance on both sides in each Location.

Like, if food_crops in the 3rd distance on the left side is lower than edible_trees in the 2nd distance on the right side, then …

I feel like there's something wrong with my random effects-term.

My output looks like this:

summary(fit)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial  ( logit )
Formula: Status ~ BankslopeScaled + Connectivity + Canal_width + Distance:Food_crops +  
Distance:Edible_trees + (1 | Distance/Side/Location)
Data: Satz

     AIC      BIC   logLik deviance df.resid 
   314.6    360.8   -144.3    288.6      245 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.18541 -0.71205  0.07243  0.82483  1.75303 

Random effects:
 Groups                   Name        Variance  Std.Dev. 
 Location:(Side:Distance) (Intercept) 2.834e-02 1.683e-01
 Side:Distance            (Intercept) 2.074e-10 1.440e-05
 Distance                 (Intercept) 2.085e-10 1.444e-05
 Number of obs: 258, groups:  Location:(Side:Distance), 258; Side:Distance, 6; Distance, 3

 Fixed effects:
                        Estimate Std. Error z value Pr(>|z|)    
 (Intercept)            -2.86517    0.79747  -3.593 0.000327 ***
 BankslopeScaled         1.76475    0.62541   2.822 0.004776 ** 
 Connectivity            0.10394    0.02729   3.809 0.000140 ***
 Canal_width             0.19138    0.11089   1.726 0.084364 .  
 Distance1:Food_crops    0.03667    0.09366   0.391 0.695441    
 Distance2:Food_crops    0.10852    0.08996   1.206 0.227694    
 Distance3:Food_crops    0.06303    0.08502   0.741 0.458510    
 Distance1:Edible_trees  0.02273    0.01327   1.712 0.086818 .  
 Distance2:Edible_trees -0.01750    0.02992  -0.585 0.558738    
 Distance3:Edible_trees  0.09769    0.07986   1.223 0.221201    
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 [correlation of fixed effects snipped]

A point into the right direction is much appreciated!

Best Answer

A few things to keep in mind:

You probably shouldn't treat Distance as a random effect in any case (estimates of random effects with only 3 levels are very unreliable)
You already implicitly have the main effects of Distance in your model through the Distance:Food_crops and Distance:Edible_trees interactions.
... for both these reasons, you should leave Distance out of your random-effects specification: (1 | Distance:Side) + (1 | Distance:Side:Location) should work to include the two lower nested levels.
you will still probably find that the (1|Distance:Side) random effect has variance of effectively zero ... but you might as well leave it in.

Related Solutions

R – Using Nested Random Effects in LME4

Let me try to answer Q3 first:

Is it justifiable to break up the models into two? I don't think so. Instead of finding a way around the message fixed-effect model matrix is rank deficient so dropping 1 column / coefficient, we should try to figure out why that happened. Usually, what this means is that you have not enough information to estimate the specified model.

I created an example that matches your experimental design above:

set.seed(23)
## YOUR DESIGN AS DEPITCTED IN THE IMAGE
surv <- data.frame(
  Gradient = rep(c("In", "Out"), each = 30),
  Site = rep(paste("Site", c(1:4), sep = ""), each = 15),
  Transect = rep(paste("T", c(1:12), sep = ""), each = 5),
  Plate = factor(sample(c(1:60))),
  # GET SOME RANDOM NUMBERS FROM POISSON DISTRIBUTION
  Response = rpois(60, 6)
)

Now using your model a for example, we can reproduce the message fixed-effect model matrix is rank deficient so dropping 1 column / coefficient (actually that's also the case for your models b, c, d, and e):

library(lme4)
model_a <- glmer(Response ~ Gradient + Site + (1|Site),
                 data = surv, family = "poisson", nAGQ = 7)

When you look at the fixed effects output (summary(model_a), also check: model.matrix(model_a)), you can see that Site4 was dropped. This has most likely to do with the fact that your fixed effects (factors and levels) cannot be represented as unique combinations of each other. Or in other words, you have levels of Site which are not represented in levels of Gradient. This is nothing new since you already specified that Site is nested within Gradient. Now, this also means that you cannot test whether levels of Site are significantly different. Instead, you would want to quantify and account for the (random-)effect of Site (and Transect) on the fixed-effect Gradient.

Here's how I would specify the model:

model_b <- glmer(Response ~ Gradient + (1|Site) + (1|Transect), 
                 data=surv, family = "poisson")
summary(model_b)

[showing only relevant output]

Random effects:
 Groups   Name        Variance Std.Dev.
 Transect (Intercept) 0        0       
 Site     (Intercept) 0        0       
Number of obs: 60, groups:  Transect, 12; Site, 4

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.81916    0.07352  24.743   <2e-16 ***
GradientOut -0.04987    0.10530  -0.474    0.636    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

When you look at the part of the output where is says Number of obs: 60, groups: Transect, 12; Site, 4, you can see that the grouping matches your design. Also since your nesting is explicit, you don't need to specify nesting in your random-effects statement (using the \ sign) since it's implied by your data. The variance estimates of the random effects will probably change when you enter your real data.

As for Q1, I think that Plate itself isn't a level since it has no replicates and is also not repeatedly measured over time. The interaction of Transect:Plate would be the lowest level in your case, which is represented by (1|Transect) due to the explicit nesting structure in your data.

As for Q2, since Site seems to be of interest to you (included as fixed effect in your models), it could be treated as such, if you reduce its levels from four to two and replicate it in the Gradient levels In and Out. However, I don't know whether this is possible because I don't know enough about your Site factor. If that's reasonable, you could do this:

surv2 <- data.frame(
  Gradient = rep(c("In","Out"), each = 30),
  Site = rep(paste("Site", c(1:2), 2, sep = ""), each = 15),
  Transect = rep(paste("T", c(1:12), sep = ""), each = 5),
  Plate = factor(sample(c(1:60))),
  # GET SOME RANDOM NUMBERS FROM POISSON DISTRIBUTION
  Response = rpois(60, 6)
)

model_c <- glmer(Response ~ Gradient + Site  + (1|Transect), 
                 data = surv2, family = "poisson")
summary(model_c)

and potentially also allowing for random slopes between levels of Site:

model_d <- glmer(Response ~ Gradient + Site  + (Site|Transect), 
                 data = surv2, family = "poisson")
summary(model_d)

Note: This won't give you the fixed-effect model matrix is rank deficient so dropping 1 column / coefficient message any more.

I don't think that you need to choose between models by comparing the AICs in this case here. I think the most useful model given the data is model_b and perhaps model_c or model_d if this is feasible. But maybe someone else can comment if that's incorrect.

Solved – Specification of nested random effects in lme4

Yes, : is the interaction operator. For example it is specifically mentioned as such in the lme4: Mixed-effects modeling with R book draft by D. Bates (see Sect. 2.2.1.1 on Nested Factors)
It effectively says that instead of having an intercept varying among schools and among pupils within schools ((1|schools)+(1|school:pupil)), you have intercept variations only among pupils within schools. It is a rather restrictive model and I think is hard to justify because you would expect at least some minimal effect to be due to school as a factor.
The interpretation would be that you are having an intercept varying among schools, among pupils and among pupils within schools. I strongly suspect this would be an overdetermined model. It could potentially work if you have multiple individual students who attended multiple schools in your sample but even then it would be probably computationally intractable. Conceptually, you would have to say that you are able to differentiate between the variation due to school, pupil and school:pupil; which I think is a bit of a tall order for not exceptionally large and well-behaved samples.

Best Answer

Related Solutions

R – Using Nested Random Effects in LME4

Solved – Specification of nested random effects in lme4

Related Question