2
$\begingroup$

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!

$\endgroup$
8
  • $\begingroup$ Animal_ID:Animal_Species:Treatment_Group - is this a complete 3 way cross of these three variables? If that's the case, then if an animal ID corresponds to a Bighorn, then that ID crossed with Species = Barbary can't be estimated, because it doesn't exist. How about starting with a simpler random effect like (1|AnimalID)? $\endgroup$
    – Alex J
    Jul 6, 2023 at 0:47
  • $\begingroup$ Hi Alex J, it is not a complete 3-way cross, because each animal (i.e. one level in Animal_ID) is present within only one combination of Animal_Species x Treatment_Group. So Animal_ID is nested within both of these two other predictors. $\endgroup$ Jul 6, 2023 at 1:40
  • $\begingroup$ I'm not super familiar with lme4 syntax, but isn't "\" the nesting operator? cran.r-project.org/web/packages/lme4/vignettes/lmer.pdf Table2 $\endgroup$
    – Alex J
    Jul 6, 2023 at 2:04
  • $\begingroup$ I still don't think the formulation is correct. If Animal ID is nested within Animal Species, that should be (1|Species/ID). If Animal ID is nested within Treatment Group, there should also be (1|Treatment/ID). I am not sure why you would want to do the second though, given Treatment is the thing you are interested in $\endgroup$
    – Alex J
    Jul 6, 2023 at 2:06
  • $\begingroup$ Yes, the "\" operator is one way to do nesting, but another way is to use the ":" operator. Since Animal_ID is nested relative to both Animal_Species and Treatment_Group, I thought I should put all this information together in one term. However, maybe you are right and I should use two separate random effect terms for Animal_ID. $\endgroup$ Jul 6, 2023 at 2:16

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ Hi EdM, thanks for your explanation, its true I am excluding a 2-way interaction between Species x Treatment, and this is because these two predictors are only partially crossed, so I believe it is not possible to include their interaction due to over-specification. With the doubly-nested Animal_ID, maybe there are simply too many coefficients in the model to include Animal_ID as an independent predictor, so I will try your suggestion to start with a simpler model and work upwards. $\endgroup$ Jul 6, 2023 at 17:58
  • $\begingroup$ @ChrisScience you seem to have all combinations of Species x Treatment, so there should be no problem with including their interaction in the model. Again, think back to the simple one-observation-per-individual situation: there's no problem with an interaction in that situation. According to your "3-way contingency table" you have 2 Species, 3 Treatments, and 3 individuals within each of the 6 combinations. That's as crossed as you can be. That wouldn't be "crossed" if Species and Treatment were themselves random effects, but the best model here is for them to represent fixed effects. $\endgroup$
    – EdM
    Jul 6, 2023 at 19:09
  • $\begingroup$ You're right, I had written the 3-way table incorrectly, I've edited my original post with the fix. It should now show missing data in a single combination in Species x Treatment $\endgroup$ Jul 6, 2023 at 21:51
  • $\begingroup$ @ChrisScience with that missing data combination you can't include the interaction and you even have to be careful in how you code and interpret those predictors. If you use T1 as the reference level for treatment and Barbary as the reference for species, then the intercept would be the modeled estimate at the combination of T1 and Barbary that doesn't exist! That would be a problem whether or not you had random effects. That's another advantage o the bottom-up approach I recommend to thinking through the model, as it finds that type of problem very early on. $\endgroup$
    – EdM
    Jul 7, 2023 at 3:02
  • $\begingroup$ @ChrisScience I've edited the answer to provide a potential work-around: code the 5 observed species/treatment combinations as a single 5-level categorical predictor. $\endgroup$
    – EdM
    Jul 7, 2023 at 3:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.