[R-sig-ME] Logistic Regression for Matched Case-Control Data

[Viechtbauer Wolfgang (STAT)]

Thanks for looking into this. Yes, for nAGQ values above ~10, it's indeed pretty stable, regardless of the version, and it does give the right value. I was just worried about the behavior for a smaller number of quadrature points. But maybe this is nothing unusual when using glmer() for this application. However, the SE is always off, regardless of the nAGQ value.

I also tried fitting the same model with the glmmML() function from the package with the same name:

library(glmmML)
glmmML(event ~ group, cluster=id, family=binomial, method="ghq", n.points=21)

This gives b1=.8109 with SE=.4249 as expected. So, at least with respect to the SE, glmer() is doing something different here.

Best,
Wolfgang

> -----Original Message-----
> From: r-sig-mixed-models-bounces at r-project.org [mailto:r-sig-mixed-models-
> bounces at r-project.org] On Behalf Of Ben Bolker
> Sent: Tuesday, June 18, 2013 16:49
> To: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] Logistic Regression for Matched Case-Control Data
> 
> Viechtbauer Wolfgang (STAT <wolfgang.viechtbauer at ...> writes:
> 
> >
> > Dear All,
> 
> > I am trying to wrap my head around using logistic random-effects
> > regression models for the analysis of matched data and the results I
> > obtain with lmer() when using this approach. So, let's say we have
> > dichotomous outcomes for matched subjects in 2 groups (could also be
> > repeated measurements on a single subject -- the idea is the
> > same). We can write the subject-specific model as:
> 
> > P(Y_ij = 1) = alpha_i + beta x_ij,
> 
> > where Y_ij is the observed outcome (either 1 or 0) for subject j
> > (either 1 or 2) for pair i (j = 1, ..., n), x_ij=0 for group 1,
> > x_ij=1 for group 2, alpha_i is the intercept for pair i, and beta is
> > the log(OR). Estimating this model with fixed intercepts is
> > problematic, so the usual approach is to use a conditional logistic
> > model (conditioning on Y_i1 + Y_i2).
> 
>  [snip snip snip to make Gmane happy]
> 
> > library(lme4)
> > lmer(event ~ group + (1 | id), family=binomial, nAGQ=21)
> 
> > According to Agresti (2002; see also Neuhaus et al., 1994), this
> > should yield the same estimate of beta for matched pairs with a
> > non-negative sample log odds ratio (as in this example). With the
> > right number of quadrature points (e.g., nAGQ=21), this does indeed
> > yield the same value. But this appears to depend highly on what
> > value of nAGQ is chosen. Also, the SE is different. I tried fitting
> > the same model in some other software to compare results. For
> > example, in Stata using xtmelogit, I get results that match up
> > completely with those from the other approaches, including the SE of
> > hat(beta).
> 
> > So, I am wondering if the high dependence of lmer() on the number of
> > quadrature points in this particular application is expected. And is
> > there an explanation for why the SE of hat(beta) is different in
> > lmer()? I realize that clusters of size 2 is probably not something
> > that lmer() was meant for, but the theory apparently says that the
> > results should match up, so I am wondering if there is an intuitive
> > explanation for what is going on here.
> 
>   With the development version of lme4 (which might ?? be more stable
> ?) I get values of the fixed effect of 0.8109 (plus or minus about
> 3e-5) as long as nAGQ >= 9 ...  for lme4.0 (which should be equivalent
> to the CRAN version) it seems to take _slightly_ higher nAGQ to get
> precise answers, but it still doesn't seem that sensitive to me.
> (Might vary across platforms?)  However, I can't speak to the standard
> error of hat(beta) ...
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models