Specifying a complicated Mixed Effects Model

Question

I am trying to translate the design information of an experiment into a mixed effects model, but I don't have any experience with mixed effects models. I was hoping someone could help me with this problem.

---

I have an experiment with 5 predictors:

• Animal_ID :: [Categorical, Unordered, 18 levels] --> (1, 2, ..., 18)
• Animal_Species :: [Categorical, Unordered, 2 levels] --> (Bighorn, Barbary)
• Treatment_Group :: [Categorical, Unordered, 3 levels] --> (T1, T2, T3)
• Sample_Type :: [Categorical, Unordered, 2 levels] --> (Blood, Tissue)
• Collection_Time :: [Numerical, Continuous] --> (-Inf -> +Inf)

---

The Animal_ID predictor is the identifier for a single animal, and these animals are randomly chosen from larger populations, therefore Animal_ID is a Random Effect Predictor. All the other predictors are Fixed Effect Predictors. Thus:

• Animal_ID --> Random
• Animal_Species --> Fixed
• Treatment_Group --> Fixed
• Sample_Type --> Fixed
• Collection_Time --> Fixed

---

Looking at contingency tables of the data set, I made the following Crossed/Nested observations for each pair of categorical predictors:

• Animal_ID is Nested within Animal_Species
• Animal_ID is Nested within Treatment_Group
• Animal_ID is Partially Crossed with Sample_Type (i.e. crossed but with some missing data)
• Animal_Species is Partially Crossed with Treatment_Group
• Animal_Species is Fully Crossed with Sample_Type (i.e. crossed with no missing data)
• Treatment_Group is Partially Crossed with Sample_Type

---

Furthermore, looking at a 3-way contingency table, Animal_ID is Nested within the Crossed predictors Animal_Species and Treatment_Group, which can be depicted in a table:

Animal_Species : Treatment_Group : Animal_IDs

Bighorn :: T1 :: 1, 2, 3, 4

Bighorn :: T2 :: 5, 6, 7, 8

Bighorn :: T3 :: 9, 10, 11, 12

Barbary :: T1 :: missing data

Barbary :: T2 :: 13, 14, 15

Barbary :: T3 :: 16, 17, 18

---

Before attempting to translate the design information into a mixed model, I make several assumptions:

• Since Animal_ID is Partially Crossed with Sample_Type, their 2-way interaction cannot be added to the model
• Since Animal_Species is Partially Crossed with Treatment_Group, their 2-way interaction cannot be added to the model
• Since Animal_Species is Fully Crossed with Sample_Type, their 2-way interaction can be included in the model
• Since Treatment_Group is Partially Crossed with Sample_Type, their 2-way interaction cannot be added to the model

---

Finally, writing a mixed model formula (in LME4 notation) based on the given assumptions, and the fact that Animal_ID is nested, I am inclined to write down something like this:

Response ~ (1|Animal_ID:Animal_Species:Treatment_Group) + Animal_Species + Treatment_Group + Sample_Type + Collection_Time + Animal_Species:Sample_Type

But this throws the error "boundary (singular) fit".

---

In conclusion, despite my observations of the design information, the mixed model I think should work does not, and I don't know what I'm missing. If anyone can help point me in the right direction, I would appreciate any suggestions. Thanks!

Answer 1

When you get an error like that, you probably are trying to estimate more coefficients than the data will allow. You might be over-thinking this into a model that's more complicated than necessary (or even possible).

In particular, the : separator you use is for different levels of grouping of random effects (see the GLMM FAQ page) so your invocation of it to include your fixed effects Animal_Species and Treatment_Group might not be doing what you think. If you want the associations of those fixed-effect predictors with outcome to vary among individuals, then you use random slopes among individuals for the fixed-effect predictors.

Instead of trying to figure out the structure of the entire data set, build up from a simple situation in which you only have one type of sample and one observation time. Then the ID values wouldn't matter, with only one observation per animal.* That could be a simple linear model of a continuous outcome based on two categorical predictors and their interaction (also known as a two-way analysis of variance):

Response ~ Animal_Species + Treatment_Group + Animal_Species:Treatment_Group

Your model didn't include that interaction term. Was that a deliberate choice?

The edited question explains what's going on, bringing up a potentially bigger problem: there are NO observations that include both T1 and Barbary, so that a simple interaction setup isn't possible and you need to exert some extra effort. In that case, you could code the 5 separate observed treatment/species combinations as a single categorical predictor, and use post-modeling analysis to compare treatment/species combinations of interest that actually exist.

Now, say that you have multiple samples over time and you want to include Collection_Time in the model. This is where you have to start making choices with respect to random effects, representing differences among individuals the share the same values of Animal_Species and Treatment_Group.

If you think that baseline outcome values differ among individuals but the associations of other predictors with outcome are the same beyond that, then you need a random intercept: a term of (1|ID) added to the model. For example:

Response ~ Collection_Time + Animal_Species + Treatment_Group + Animal_Species:Treatment_Group + (1|ID)

If you think that the association of Treatment_Group with outcome might differ among individuals, then you could include a random slope for Treatment_Group. As shown on this site's lmer cheat sheet, you could allow those random individual-specific slopes to be uncorrelated with the random intercepts by adding a term (0 + Treatment_Group|ID) to the above model. If you want to impose a correlation between the intercepts and slopes instead (often a reasonable choice), replace the (1|ID) above with (1 + Treatment_Group|ID).

Proceed similarly to consider what you think is going on with respect to the fixed-effect predictors Animal_Species, Collection_Time, Sample_Type, and interactions among any or all fixed-effect predictors. In principle, each of those individual effects and interactions could be included as a random slope among ID values.

In practice, you typically run out of data before you can fit that many predictors, interactions, and random effects. You have to apply your knowledge of the subject matter to make choices about what to include and what to exclude.

---

*Although the animals are still "randomly chosen from larger populations" in this situation, with only 1 observation per individual there are no intra-individual correlations to take into account so no random-effect term is involved.