R - unexpected probabilities obtained with lsmeans using a glmer (lme4) mixed effect model

This is a follow-up question to this one.

I am working on survey data which contains multiple-response types of questions. More specifically, respondents were asked to "check all that applies" on 2 sections of the questionnaire. My client needs to compare the response proportions for each of the options listed. At first I was looking into pairwise tests for proportions, but as was pointed out, the dependency (each subject has a tendency to select few, or many items) should be taken into account, such as in a mixed effects model.

I managed to get results that make a lot of sense. After restructuring the data, I started by just using a logistic regression, combined with lsmeans to make sure that the probabilities did reflect the raw proportions for each item. And it works as expected.

library(lme4)
library(lsmeans)
mod.fixed <- glm(valBin ~ var, family=binomial, data=mydata)
lsmeans(mod.fixed, "var", type="response")

# var       prob         SE df asymp.LCL asymp.UCL
# Q11a 0.5222222 0.03723097 NA 0.4492960 0.5942133
# Q11b 0.5977654 0.03665039 NA 0.5243325 0.6670602
# Q11c 0.6500000 0.03555121 NA 0.5775582 0.7161250
# Q11d 0.5810056 0.03687803 NA 0.5074972 0.6510856
# Q11e 0.5786517 0.03701001 NA 0.5049313 0.6490249
#
# Confidence level used: 0.95
# Intervals are back-transformed from the logit scale


These correspond exactly to the raw proportions of positive answers to each item in the survey.

Going a step further, I used glmer allowing each subject to have its own intercept.

mod.mixed <- glmer(valBin ~ var + (1|ID), family=binomial, data=mydata)
lsmeans(mod.mixed, "var", type="response")
# var       prob         SE df asymp.LCL asymp.UCL
# Q11a 0.5650139 0.07870159 NA 0.4094875 0.7087161
# Q11b 0.7163331 0.06658291 NA 0.5705559 0.8275796
# Q11c 0.7983241 0.05408306 NA 0.6720504 0.8843446
# Q11d 0.6831449 0.07048872 NA 0.5324562 0.8032161
# Q11e 0.6802639 0.07094074 NA 0.5289000 0.8012684
#
# Confidence level used: 0.95
# Intervals are back-transformed from the logit scale


Up to there, the difference in probabilities between the fixed and mixed models make sense (the prob's from the mixed model are not completely off the map). However, I have a second set of items that follow the exact same logic, so I followed the same steps:

mod2.fixed <- glm(valBin ~ var, family=binomial, data=mydata2)
lsmeans(mod.fixed, "var", type="response")
# var       prob         SE df asymp.LCL asymp.UCL
# Q12a 0.8222222 0.02849686 NA 0.7593968 0.8714205
# Q12b 0.8722222 0.02488306 NA 0.8150493 0.9135951
# Q12c 0.8722222 0.02488306 NA 0.8150493 0.9135951
# Q12d 0.8388889 0.02740178 NA 0.7777517 0.8856801
# Q12e 0.7333333 0.03296086 NA 0.6640175 0.7928108
#
# Confidence level used: 0.95
# Intervals are back-transformed from the logit scale


The prob values still correspond precisely to the raw observed proportions.

Now the mystery comes in in the next step, where I use glmer and lsmeans with that second set of data.

mod2.mixed <- glmer(valBin ~ var + (1|ID), family=binomial, data=mydata2)
lsmeans(mod.mixed, "var", type="response")
# var       prob           SE df asymp.LCL asymp.UCL
# Q12a 0.9995412 0.0004351133 NA 0.9970612 0.9999285
# Q12b 0.9998568 0.0001464262 NA 0.9989384 0.9999807
# Q12c 0.9998568 0.0001464273 NA 0.9989384 0.9999807
# Q12d 0.9996788 0.0003118857 NA 0.9978485 0.9999521
# Q12e 0.9977529 0.0017893754 NA 0.9893514 0.9995290
#
# Confidence level used: 0.95
# Intervals are back-transformed from the logit scale


I can't explain why all the probabilities are near 1, whereas in the data, they range from .73 to .87. And also, why are the SE's so low? I'm linking sample, anonymous data here. (Just ignore the smok column)

I've checked for anomalies in the data and I don't see anything dramatically different from mydata1 in mydata2. If anyone can help me understand this phenomenon, I'll be more than grateful.

I verify the same results. Notable is the comparison of estimates (i'm calling the models fitted to the second dataset mod2.fixed and mod2.mixed):

> mod.mixed
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial  ( logit )
Formula: valBin ~ var + (1 | ID)
Data: mydata
AIC       BIC    logLik  deviance  df.resid
941.4224  970.2100 -464.7112  929.4224       890
Random effects:
Groups Name        Std.Dev.
ID     (Intercept) 2.853
Number of obs: 896, groups:  ID, 180
Fixed Effects:
(Intercept)      varQ11b      varQ11c      varQ11d      varQ11e
0.2615       0.6648       1.1143       0.5067       0.4934

> mod2.mixed
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: binomial  ( logit )
Formula: valBin ~ var + (1 | ID)
Data: mydata2
AIC       BIC    logLik  deviance  df.resid
556.1389  584.9532 -272.0694  544.1389       894
Random effects:
Groups Name        Std.Dev.
ID     (Intercept) 8.02
Number of obs: 900, groups:  ID, 180
Fixed Effects:
(Intercept)      varQ12b      varQ12c      varQ12d      varQ12e
7.6865       1.1645       1.1645       0.3566      -1.5906


The model for the second dataset has a much greater estimated SD for subjects -- and an absolutely huge intercept of 7.69 compared with 0.26. I don't think the small SEs are a concern, as they are a byproduct of the distorted intercept. To see this, I changed just the intercept and redid the summary:

> lsm2= lsmeans(mod2.mixed, "var")
> lsm2@bhat
(Intercept)     varQ12b     varQ12c     varQ12d     varQ12e
7.6864525   1.1644807   1.1644800   0.3566375  -1.5905781
> lsm2@bhat[1] = .26
> summary(lsm2, type="resp")  ### These are fake results for testing purposes
var       prob        SE df  asymp.LCL asymp.UCL
Q12a 0.5646363 0.2332415 NA 0.16802950 0.8927991
Q12b 0.8060399 0.1598540 NA 0.35903794 0.9685836
Q12c 0.8060398 0.1598552 NA 0.35903453 0.9685840
Q12d 0.6494534 0.2211250 NA 0.21635511 0.9255534
Q12e 0.2090638 0.1319715 NA 0.05240947 0.5581554

Confidence level used: 0.95
Intervals are back-transformed from the logit scale


You can see that now the SEs are now much larger than those of mod.mixed -- reflecting the very much larger SD of subjects.

I am not sure yet why this is happening, but I do note that the mod2.fixed results look kind of suspicious too, in that two of the treatments have the same mean and all of the other means are round-looking numbers, in a way -- like something you'd get from a very small dataset. It makes me think that there are really only a few distinct subjects, and many of them are copies of one another.

• This is very helpful and you're right, the proportions seem suspicious, I'll look into that. – Dominic Comtois Oct 27 '16 at 3:35
• I think I know what happens... Out of 180 subjects, 119 (66%) have 1's everywhere. That must be what makes the model go crazy... – Dominic Comtois Oct 27 '16 at 19:21