Describing data structure and specifying a linear mixed model in nlme with nested and crossed effects

Question

I am trying to specify a linear mixed model to analyse data with the following structure and have several questions about correctly describing the structure of the data and how to specify the model. I often have similar questions when trying to specify mixed models and this example gives me the opportunity to query some of my (mis)understandings.

A response variable - lets call it metabolic rate (MR) - of animals was measured under 6 different conditions (A through F). Animals were split into two groups (G1 and G2) with each animal tested under 3 of the conditions only. The order of testing (conditions) was dependent on the group the animal was in. The whole experiment was run twice (R1 and R2). Unfortunately, the order of testing within each group was different between Run 1 and Run 2 such that animals were tested in the following order within each group:

• Run 1, Group 1: Conditions A, then B then C,
• Run 1, Group 2: Conditions D, E, F,
• Run 2, Group 1: Conditions F, B, A,
• Run 2, Group 2: Conditions C, E, D.

Each Run and Group combination used a different set of 10 animals, each with a unique ID (1-40 across the four groups).

I have made an example dataset with this structure in R which may help:

set.seed(1)

df <- data.frame(Run = rep(c("R1", "R2"), each = 60),
                 Group = rep(c("G1", "G2"), each = 30, times = 2),
                  Condition = rep(c(rep(c("A", "B", "C"), each = 10),
                                  rep(c("D", "E", "F"), each = 10),
                                  rep(c("F", "B", "A"), each = 10),
                                  rep(c("C", "E", "D"), each = 10))),
                 ID = c(rep(1:10, 3), rep(11:20, 3),
                        rep(21:30, 3), rep(31:40, 3)),
                 MR = rnorm(120, 1, 0.1))

In describing the data structure are the following correct:

Animal ID is nested within Group?

Is Group nested within Run (as Group 1 and 2 are different between the two Runs)?

Also, is Condition is fully crossed with Run but partially crossed within Group?

I am mainly interested in whether the metabolic rate for the different conditions differed between Run 1 and Run 2 i.e. a significant Run * Condition interaction or main effect of Run. I would ideally like to specify the model using nlme in R if possible.

This is the thought process I have been working through so far:

My first attempt to model the data was:

library(nlme)

m1 <- lme(MR ~ Condition * Run, random = ~1|ID, data = df)
summary(m1)
anova(m1)

however, I think that this won’t capture the fact that MR may also depend on the order in which the animals were tested (i.e. Group). Is this true?

As Group 1 in Run 1 has a different sequence of conditions from Group 1 in Run 2 (with the same being true for group 2) I specified a new column of df, giving each group a unique identifier.

df$NewGroup <- interaction(df$Run, df$Group)

And as I think ID is nested in NewGroup I then tried the following model:

m2 <- lme(MR ~ Condition * Run, random = ~1|NewGroup/ID, data = df)
summary(m2)
anova(m2)

Is this a correct model specification and does it capture the effect of Group?

Two concerns I have are:

1. Is it a problem that the conditions tested within each group are different?
2. NewGroup has only four levels – is this a problem for a random effect?

Given ID is also nested within Run, why does that not also need to be specified (or is this implicit from the fact that Run is included as a main effect?).

The original plan was to collect all the data in just two runs, with each animal being tested under all six conditions, just with a different treatment order (order of conditions) between the runs. There would be only one group of 20 animals per run. In this case, although I think animal is still nested within run, would the model m1 (above) be appropriate?

I think a summary of where I am unsure how to proceed may be that although there is a mixture of crossing and nesting of effects (both random/fixed), nesting only needs to be specified for a random effect when it involves a second random effect - is this true here, and is this always the case?

Edit

Following suggestions in the comments I also tried some models with NewGroup as a main effect, but the model won't run with both Run and NewGroup, regardless of whether I include the interaction between these two factors or not.

e.g.

m3 <- lme(MR ~ Condition * Run + NewGroup, random = ~1|ID, data = df)

Gives the following error:

Error in MEEM(object, conLin, control$niterEM) : Singularity in backsolve at level 0, block 1

I am wondering if this is because of the incomplete nature of the design (i.e. within Run 1, Group 1 and Group 2 only have half of the conditions each, and the same within Run 2).

Answer 1

There are quite a lot of questions and issues here. I will do my best to work through them one by one.

Understanding the data structure.

First, I would like to point out, and you may be well aware of it, but it is always good to remember that nesting or crossing of random factors is a property of the experimental design, and therefore the data. It is not a property of the model.

• ID: A unique identifier for each animal which is is tested under three conditions and belongs to one specific group in one specific run.
• Group: There are two groups (G1 and G2) per run.
• Run: There are two runs (R1 and R2). The order of conditions tested within each group is different in each run.
• Condition: There are six different conditions (A through F), but each animal is tested under only three conditions, depending on its group and run.

Nested vs Crossed random effects.

(For some detailed background on this very issue see my answer here: Crossed vs nested random effects: how do they differ and how are they specified correctly in lme4?)

• Since each animal belongs to only one group, ID is nested within Group.
• Since each Group (G1 and G2) is specific to each Run, Group is nested in Run.
• since each Condition appears in both Runs, these factors are crossed (and appear to be fully crossed).
• since each Condition appears in both runs but in different groups and sequences, Condition is crossed with Run
• since each Group only experiences three of the six Conditions, Condition is (partially) crossed with Group

However, most of that is not important.

Now to address your specific questions

Is it a problem that the conditions tested within each group are different?

No, this is not a problem. Each animal is exposed to different conditions based on their Group and Run, which implies that the response variable (MR) could vary due to these differences. By including Condition as a fixed effect and properly structuring the random effects, you can account for the variability introduced by these different conditions.

NewGroup has only four levels – is this a problem for a random effect?

Yes and I will talk about that below

Error in MEEM(object, conLin, control$niterEM) : Singularity in backsolve at level 0, block 1

I am wondering if this is because of the incomplete nature of the design (i.e. within Run 1, Group 1 and Group 2 only have half of the conditions each, and the same within Run 2).

The error arises due to collinearity or redundancy between the Run and NewGroup variables, as NewGroup already contains information about Run and Group. Including both Run and NewGroup in the fixed effects leads to multicollinearity, causing singularity issues.

To avoid this, you can't include Run and NewGroup as fixed effects in the model. Nor can you include just Group along with the interaction between Run and Condition (which you need to answer your research question)

The issue can be seen from this cross-tabulation between Condition and Group:

> table(df$Condition, df$Group)
   
    G1 G2
  A 20  0
  B 20  0
  C 10 10
  D  0 20
  E  0 20
  F 10 10

Given your research question:

I am mainly interested in whether the metabolic rate for the different conditions differed between Run 1 and Run 2 i.e. a significant Run * Condition interaction or main effect of Run. I would ideally like to specify the model using nlme in R if possible.

Given the details above this implies the following models:

m1 < lme(MR ~ Run*Condition, random = ~1|NewGroup/ID, data = df)

or

m2 <- lme(MR ~ Run*Condition, random = ~1|ID, data = df)

The problem with m1 is that NewGroup has only 4 levels.

The problem with m2 is that it does not take into account the variability in the outcome due to the order in which the animals were tested.

So, neither model is ideal, but I am out of ideas for any alternative. Perhaps others can chime in with ideas ?