I have data containing 40,000 observations, coming from around 30,000 subjects. Most subjects contribute one observation, some contribute two.
The dependent variable is binary (1/0), and let's say, for the sake of the discussion, that I only have one independent variable, also binary. I want to model the effect of X on Y, and to get an odds ratio.
I started by calculating the odds ratio by a simple 2X2 table, and got a crude OR of 3.18. Since there is dependency between some observations, I ran the 2X2 analysis for the first observation of the subjects and then for the second, separately, and got OR1 = 3.09 and OR2 = 4.36. The Mantel Haenszel is OR_m-h = 3.11. However, when I run a mixed effects logistic regression, I get higher OR's. Using the Laplace approximation, I get an OR of ~4.5. If I use the quadrature method of integration, I get an OR of over 6.
This makes no sense to me. Can someone please explain to me how can this be possible? My random effect (subjectID) has thousands of levels...
My R code of the model is:
glmer(Y ~ X + (1 | SubjectID), data = dataset, family = binomial, control = glmerControl(optimizer = "bobyqa"), nAGQ = 1)