Re: [R] lmer and mixed effects logistic regression

From: Spencer Graves <spencer.graves_at_pdf.com>
Date: Sun 18 Jun 2006 - 02:46:41 EST

'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@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@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


R-help@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 Received on Sun Jun 18 02:54:23 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Tue 20 Jun 2006 - 00:11:01 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.