# Using glmer to estimate treatment interactions

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

• I guess that in the second case you have a more sophisticated random effects structure than you expect, specifically some correlations among the random effects that explain variance. Those effect help to achieve more accurate estimations of the fixed effects, hence the differences. Can you post the random effects output of all the models (i.e., the upper part of summary) so that we can confirm that? – Henrik Apr 7 '14 at 20:37

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))?
• @user43153 Hmm, that is weird. From the output you posted it looks like there should be only 2 random effects parameters for each of the initial models, while, as I expected, model containes the additional correlation parameters which drive the difference. You might try using || instead of | in your call to produce model. Without a reproducible example it now becomes difficult to see what is going on. – Henrik Apr 8 '14 at 7:01