Solved – Using glmer to estimate treatment interactions

Question

In my data, I have two treatment conditions with repeated measures for each subject. I would like to run a mixed logistic regression separately for each of my two conditions where my binary outcome DV (dependent variable) is regressed on my IV (independent variable) and also have a random slope and intercept fitted for each subject.

So, I run the following:

modelT0 <- glmer(DV ~ IV + (1|subject) + (0 + IV|subject), data = D0, family = binomial)
modelT1 <- glmer(DV ~ IV + (1|subject) + (0 + IV|subject), data = D1, family = binomial)

In the above, D0 and D1 are data sets restricted to treatment conditions 0 and 1, respectively. What I would like to do is compare the estimated fixed effects coefficient on IV across conditions to see if it significantly changes. To do this, I pool D0 and D1 into a single data set, D, and create a treatment indicator that takes value 0 in D0 and 1 in D1. I then run:

model <- glmer(DV ~ IV + treatment + treatment:IV + (1 + treatment|subject:treatment) 
               + (0 + IV + treatment:IV|subject:treatment), data = D, family = binomial)

I should be able to look at the fixed effects coefficient on treatment:IV to get my answer, but the issue is that for whatever combination of random effects I seem to specify, the coefficients from the pooled regression are slightly different from the regressions specified by treatment. So for instance, the fixed effect coefficient on treatment:IV plus the one on IV in model is not equal to the coefficient on IV in model1.

Any idea about what I might be doing wrong or how to answer the question I have? Thanks!

EDIT:

As per Henrik's suggestion, I'm copying the random effects output of the models below:

summary(modelT0):

    Random effects:
    Groups    Name        Variance  Std.Dev. 
    subject   (Intercept) 1.412e-07 0.0003758
    subject.1 IV          1.650e+00 1.2844341

summary(modelT1):

    Random effects:
    Groups    Name        Variance Std.Dev.
    subject   (Intercept) 0.00378  0.06148 
    subject.1 IV          0.26398  0.51379 

summary(model):

    Random effects:
    Groups              Name         Variance  Std.Dev. Corr 
    subject.treatment   (Intercept)  0.0005554 0.02357       
                        treatment    0.0066042 0.08127  -0.88
    subject.treatment.1 IV           1.6500112 1.28453       
                        IV:treatment 1.0278663 1.01384  -0.93

Answer 1

As I said in my comment, I expect the problem to be one of additional random effects parameters in model compared to in modelT0 + modelT1 (specifically corraletions).

Hence I would first check for the number of random effects parameters. Is attr(logLik(model), "df") - length(fixef(model)) equal to 2 * attr(logLik(modelT0), "df") - length(fixef(modelT0))?

If not, the additional random effects parameters explain error variance which helps to more precisely estimate the fixed effects.