[R] lmer and mixed effects logistic regression

[Douglas Bates]

If I understand correctly Rick it trying to fit a model with random
effects on a binary response when there are either 1 or 2 observations
per group.  I think that is very optimistic because there is so little
information available per random effect (exactly 1 or 2 bits of
information per random effect). I'm not surprised that there are
difficulties in fitting such a model.

Rick should try to monitor the iterations to see what is happening to
the parameter estimates and perhaps should try the Laplace
approximation to the log-likelihood.  This is done by adding method =
"Laplace" and control = list(msVerbose = TRUE) to the call to lmer.
This causes the PQL iterations to be followed by direct optimization
of the Laplace approximation.  Because the problem is probably in the
PQL iterations Rick may want to turn those off and do direct
optimization of the Laplace approximation only.  In that case the
control argument should be list(usePQL=FALSE, msVerbose = TRUE)


On 6/17/06, Spencer Graves <spencer.graves at pdf.com> wrote:
>
>
> 'lmer' RETURNS AN ERROR WHEN SAS NLMIXED RETURNS AN ANSWER
>
>           Like you, I would expect lmer to return an answer when SAS
> NLMIXED does, and I'm concerned that it returns an error message instead.
>
>           Your example is not self contained, and I've been unable to
> get the result you got.  This could occur for several reasons.  First,
> what do you get for "sessionInfo()"?  I got the following:
>
> sessionInfo()
> Version 2.3.1 (2006-06-01)
> i386-pc-mingw32
>
> attached base packages:
> [1] "methods"   "stats"     "graphics"  "grDevices" "utils"     "datasets"
> [7] "base"
>
> other attached packages:
>        lme4    lattice     Matrix
>   "0.995-2"   "0.13-8" "0.995-10"
>
>           If you are using different versions of lmer, lattice and / or
> Matrix, please update.packages and try again.  If that does not fix the
> problem, could you please try to find a self-contained example that is
> as simple as possible and still generates the same error message?  Then
> send that to this list.
>
>           When I ran your example limited only to the observations you
> listed, I got a different result:
>
> Warning messages:
> 1: fitted probabilities numerically 0 or 1 occurred in: glm.fit(X, Y,
> weights = weights, offset = offset, family = family,
> 2: IRLS iterations for PQL did not converge
>
>           We get this these "Warning messages" because all the values
> of "response" are constant within each level of "id".  With data like
> this, the likelihood will be maximized by sending Std.Dev.(id) to Inf;
> when 'lmer' quit, Std.Dev.(id) was already at 58779, which
> computationally is moderately close to Inf for the current "lmer"
> algorithm.
>
>           How many different patterns of "response" by "id" do you have?
> With all 0's, the random adjustment to "(Intercept)" for that "id",
> ignoring the Bayesian shrinkage, is "-Inf".  With all 1's, this random
> adjustment would be "+Inf".  With two observations for a subject, if
> they are different, this random adjustment would exactly cancel the
> estimate of the fixed-effect for "(Intercept").  Do you have only these
> three patterns, all 0's, all 1's, and half 0's, half 1's?  If yes, that
> might explain why Std.Dev.(id) wants to go to Inf.
>
>           I don't know exactly how many different patterns you need, but
> you should be able to have some levels of "id" with all 0's or all 1's
> as long as those did not dominate, as they do in this cut down example.
>
>           If this is the problem, then SAS NLMIXED could be achieving false
> convergence and either not telling you about it or reporting it so
> quietly you have not reported it here.
>
>
> CREATING A REPRODUCIBLE EXAMPLE
>
>           I'm not eager to experiment with a large complicated example.
> You could help us help you by trying to produce a simpler,
> self-contained example.  For example, do you get the same answer when
> you delete 'age' from the model:
>
>           lmer(response~(1|id),data=gdx,family=binomial)
>
>           If yes, then you could easily just count the number of levels of "id"
> that have all 0's, all 1's, (0, 1) = (1, 0), etc.  Then write one or two
> lines that would generate that special case and submit that to this
> listserve.  Before you do, however, I encourage you to experiment to
> find small numbers of replicates that reproduce the error you get.  Then
> submit that to this list.
>
> STARTING VALUES?
>
>           The help page for "lmer" describes the following argument:
>
>     start: a list of relative precision matrices for the random effects.
>             This has the same form as the slot '"Omega"' in a fitted
>            model.  Only the upper triangle of these symmetric matrices
>            should be stored.
>
>           To understand this, it might help to experiment with one of
> the examples on this help page:
>
> > fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)
> > fm1 at Omega
> $Subject
> 2 x 2 Matrix of class "dpoMatrix"
>              (Intercept)       Days
> (Intercept)   1.0746247 -0.2942832
> Days         -0.2942832 18.7549595
>
>           Just taking a guess, I tried the following:
>
> fm1. <- lmer(Reaction ~ Days+(Days|Subject), sleepstudy,
>               start=list(Subject=diag(1:2)))
>
> fm1.c <- lmer(Reaction ~ Days+(Days|Subject), sleepstudy,
>               start=list(Subject=array(1, .1, .1, 2),
>                       dim=c(2,2)))
>
>            Both of these returned the same answer.  Starting values are
> only required for the random effects, because these models are all
> "plinear" in the sense described in the "nls" documentation.
>
>           For "glmmPQL", I was able to get the "update" function to work.
> Thus, if you can get an answer with one model, you can use that as an
> input to "update".  I would expect "update" to try to use the parameter
> estimates from one model fit as starting values for the modification,
> where it can find a way to do that.
>
>           Hope this helps.
>           Spencer Graves
>
> Rick Bilonick wrote:
> > I'm using FC4 and R 2.3.1 to fit a mixed effects logistic regression.
> > The response is 0/1 and both the response and the age are the same for
> > each pair of observations for each subject (some observations are not
> > paired). For example:
> >
> > id response age
> > 1    0      30
> > 1    0      30
> >
> > 2    1      55
> > 2    1      55
> >
> > 3    0      37
> >
> > 4    1      52
> >
> > 5    0      39
> > 5    0      39
> >
> > etc.
> >
> > I get the following error:
> >
> >> (lmer(response~(1|id)+age,data=gdx,family=binomial))
> > Error in deviance(.Call("mer_coefGets", x, pars, 2, PACKAGE =
> > "Matrix")) :
> >         Leading minor of order 2 in downdated X'X is not positive
> > definite
> >
> > Similar problem if I use quasibinomial.
> >
> > If I use glm, of course it thinks I have roughly twice the number of
> > subjects so the standard errors would be smaller than they should be.
> >
> > I used SAS's NLMIXED and it converged without problems giving me
> > parameter estimates close to what glm gives, but with larger standard
> > errors. glmmPQL(MASS) gives very different parameter estimates.
> >
> > Is it reasonable to fit a mixed effects model in this situation?
> >
> > Is there some way to give starting values for lmer and/or glmmPQL?
> >
> > Rick B.
> >
> > ______________________________________________
> > R-help at stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
>