CODEWITHSUNDEEP
rmodelmixedglmm
Jessica Girona
Jessica Girona

Reputation: 21

Random effects nested

I'm pretty new to the world of mixed models, and I'm looking for validation from people who are more knowledgeable than I am.

I have a randomized complete block design with different 8 treatments repeated 4 times. I've been monitoring aphid populations on plants in the field (count) at different dates. I want to know if one of the treatments was more effective than another (in other words, in which treatment did I get fewer aphids?). The trial was carried out 7 times (2 sites in 2022, 3 sites in 2023 and 2 sites in 2024).

To sum up, I have the following variables :

Aphids : Count 
Plant : 10 Plant per plot (new at each day of assessment)
Plot : 32 per site 
Block : 4 levels per site
Treatment : 8 levels per site 
Site : 7 levels (field where is located the trial)
Date : the day of the assessment
Year : 3 levels

I guess I have to put in random effect:

These random effects are nested, so I tried to respond at my question with the following model :

modele.abond <- glmmTMB(Aphid ~ Treatment + (1|Year) + (1|Year/Site) + (1|Year/Site/Block) + 
                          (1|Year:Site:Block), 
                        data = DF, 
                        family = nbinom2)

Given the hierarchical structure, does this model seem relevant to you?

Many thanks in advance !

Upvotes: 1

Views: 67

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226742

I would suggest something like

Aphid ~ Treatment*Year +  cs(1+Treatment|Year:(Site/Block))
  • it's conceptually sensible to treat year as a random effect, but with only three years you're unlikely to be able to fit the model effectively (or, you could equivalently say that it's hard to generalize to other years if you only have a sample of 3 years), so I'm suggesting treating it as a fixed effect (and allowing for the effect of treatment to vary across years)
  • as @Axeman says, a syntax like a/b/c expands to a + a:b + a:b:c, so your specification is redundant. Here I've written Year:(Site/Block) which will expand to Year:Site + Year:Site:Block (since the top-level effect of Year is already in the fixed effect specification)
  • by specifying cs(1 + Treatment|...) we're allowing the effect of treatment to vary across blocks, with each treatment having its own variance but every pair of treatments having the same correlation. If you want to restrict this to a homogeneous compound symmetric matrix (every treatment has the same between-site and between-block variance), you can use the map argument, but it's a little bit tricky:
ntreat <- 8
mapvec <- list(theta = factor(rep(1:4, c(ntreat, 1, ntreat, 1))))

Since the parameters for each random effect term (Y:S and Y:S:B) are specified with the ntreat log-SDs followed by one correlation parameter, this says that the ntreat log-SDs for each term should all be estimated as the same value.

Upvotes: 7

Related Questions