# Mixed model that tests for within-subject differences between two conditions while including covariates

I would like to use a mixed model to test if the values in a response variable differ between two conditions that each subject underwent. On top, I would like to control for the influence of other continuous and categorical variables (without having any hypothesis on the specific influence of these covariates on the value differences between the conditions, aka. just include them as covariates). The "minimum" model would be to check if values for condition A are greater then for condition B and to include the variable "site" as a covariate of no interest to control for batch-effects (subjects are coming from different sites, and the difference between condition A and B might differ between sites).

Questions:

1. What would be the corresponding lmer formula for this? I wrote value ~ condition + (1|subject) + (1 | site), but I am not sure if this is right?

2. I have the directed hypothesis, that for condition A values will be higher than for condition B. How do include this within my model? Do I simply divide the p-value in half after the model has been fitted?

I assume that the data snippet is representative of the full dataset: every participant provides one response under each condition and that site, age and sex are participant-level covariates.

To test the hypothesis that responses are higher under A than under B, let's first estimate the condition effect: the difference between A and B adjusting for age and sex and taking into account the "subjects nested with sites" structure. So the fixed component of the linear mixed model is y ~ age + sex + condition and the random component is 1 | site/subject. This specification works "out of the box" so to say, whether there is one observation under each condition per subject or multiple — not necessarily balanced — observations per condition & subject. We just plug in the data into lme4::lmer (or nlme::lme or glmmTMB::glmmTMB or MASS::glmmPQL. For example:

m1 <- lme4::lmer(
y ~ age + sex + condition + (1 | site/subject)
)


And if we declare the condition as factor with ordered levels c("B", "A"), with summary(m1) we get an estimate for the effect of condition A compared to that of condition B, i.e. an estimate for the difference A - B. We also get a t-statistic.

The goal of this analysis is to calculate a p-value when the alternative hypothesis is A > B. The challenge is not that the alternative is one-sided (it's easy to calculate the p-value by hand) but to determine the degrees of freedom of the t-statistic. For a linear mixed model this is not straightforward at all: getting degrees of freedom from lmer, How to determine Degrees of Freedom in Linear (Mixed Effect) Regression, Degrees freedom reported by the lmer model don't seem plausible, How to obtain the p-value (check significance) of an effect in a lme4 mixed model?. And unfortunately no agreement on the best practice emerges from those discussions either.

However, if the dataset is large enough (and the model fits the data well), the t-statistic is approximately normal. Then the choice for degrees of freedom formula doesn't make a practical difference. In a very imprecise way, by "large enough" I mean the number of subjects n is at least 50 more than the number of sites k. (Note #1: I increased this from 30 to 50; there isn't a hard-set rule anyway for when the $$t$$ is "sufficiently" close to the normal; see also Why does the t-distribution become more normal as sample size increases?. Note #2: We also need degrees of freedom for the covariates, so let's say n ≥ 50 + k + 2 as the age and sex covariates each have one degree of freedom.)

Finally compute the p-value for the greater-than test as:

# p-value = Pr{T < statistic}
# DF can vary depending on the method used to fit the model
pt(statistic, df, lower.tail = FALSE)

# or

# p-value = Pr{Z < statistic}  (normal approximation)
pnorm(statistic, lower.tail = FALSE)


The model y ~ age + sex + condition + (1 | site/subject) works well indeed, and since it doesn't need the data to be balanced, is probably the way to go about the analysis. Nevertheless, I also show how to estimate the difference between the two conditions with a model for the within-subject differences. This analysis has more steps and originally I had missed one of those steps — now fixed — which can be instructive. (Or at least it was instructive for me to make the mistake and then figure out how to fix it.)

First let's re-arrange the data so that each row corresponds to a subject:

#>   subject site    age   sex     A     B
#>         1 NY       23     1  17.1  20.6
#>         2 LA       17     0  13.1  14.3
#>         3 CHI      32     1  15.5  11.2
#>         4 CHI      26     1  21.2  20.8


This format highlights that the experimental design is paired as each participant provides two responses, one under A and another under B. Since the interest is whether A results in higher responses and the experimental condition varies within subject, the "basic estimators" of the difference between responses under A and B, $$\Delta = \mu_A - \mu_B$$, are the within-participant differences $$\delta_i = Y_{i,A} - Y_{i,B}$$. The hypothesis $$\mu_A > \mu_B$$ becomes $$\Delta > 0$$.

Let's add these differences to the data table:

#>   subject site    age   sex     A     B  delta
#>         1 NY       23     1  17.1  20.6 -3.5
#>         2 LA       17     0  13.1  14.3 -1.20
#>         3 CHI      32     1  15.5  11.2  4.3
#>         4 CHI      26     1  21.2  20.8  0.400


Now we can fit a mixed model for the difference between the two conditions, with covariates age and sex and a random site effect.

lmer(
delta ~ age + sex + (1 | site)
)


In this formulation, condition doesn't appear in the model formula and the intercept term corresponds to the difference $$\Delta = \mu_A - \mu_B$$ when the covariates are "set to 0". However, age = 0 doesn't make sense and in any case we want an estimate of the condition effect that's representative of the entire population under study (young adults).

So we center the continuous covariate age to its mean and use the sum contrast for the categorical covariate sex. In R formula notation:

m2 <- lmer(
delta ~ I(age - mean(age)) + sex + (1 | site),
contrasts = list(sex = "contr.sum")
)


Note: Centering the covariates is the step I had omitted in my original answer.

• Perfect, thanks so much! Is there a specific reason why you computed delta, instead of including the condition as a categorical predictor? Or is it just a matter of taste and both approaches would lead to the same result? Commented Mar 11 at 15:20
• I find this formulation more elegant: We want to estimate the average difference between A and B, and the within-participant differences are the "natural" estimators. The mixed-effects model with (1 | subject) component should give the same result. (Something to try for fun?) Make sure the subject IDs are unique or use (1 | site/subject) otherwise. Also, note the similarity with a paired t-test. Commented Mar 11 at 17:52
• I decided to simulate some data so that I can compare the models M1: lme(delta ~ age + sex, random = ~ 1 | site) and M2: y ~ condition + age + sex, random = ~ 1 | site / subject. The estimates (intercept in M1, beta_condition in M2) are not the same but both are valid estimators of the true difference. However, the std. error is smaller in M2. The smaller std. error the better; I don't understand why yet. Commented Mar 11 at 21:20
• Okay, I know why: the presence of covariates changes the meaning of the intercept in M1 in a subtle way. This model specification still works, just have to be more careful about the covariates (age and sex). Will revise my answer; M1 and M2 give the same answer as expected (if specified correctly). Thanks for asking the followup question! Commented Mar 11 at 21:57