From: Douglas Bates <bates_at_stat.wisc.edu>

Date: Sun 18 Jun 2006 - 21:58:51 EST

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 22:03:06 2006

Date: Sun 18 Jun 2006 - 21:58:51 EST

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.

On 6/17/06, Spencer Graves <spencer.graves@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@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 22:03:06 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 - 02:11:15 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.
*