[R] lmer and mixed effects logistic regression

[Spencer Graves]

'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