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
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
mydata2. If anyone can help me understand this phenomenon, I'll be more than grateful.