[R-sig-ME] cloglog and lmer
Ben Bolker
bolker at ufl.edu
Fri Mar 19 20:25:44 CET 2010
What are the results of sessionInfo() ?
I am using a slightly hacked version of lme4, but I don't *think* that
any of my hacks should affect the functioning of cloglog.
You ran exactly the toy code example that I provided, right?
It's possible that this problem is right on the edge of numerical
tractability and thus that there are some tiny platform-dependent
differences that lead to a probability estimate that is very small on
one platform and exactly 0 (due to underflow) on another platform,
triggering the error.
In fact the error message you got
... mu = 9.35952e-313, i = 1068030290 ...
suggests this may be the case -- although I don't see how mu =
9.3592e-313 is triggering a (mui<=0 || mui>1) flag ...
It also suggests there may be another bug in the reporting because I
doubt there are really 1068030290 cases in the data.
However, there were some old bugs of this type that do not appear in
the current (at least in the development) version of lme4, so I strongly
recommend that you at least check to see if you have the latest released
version, and perhaps switch to the development version
(install.packages("lme4",repos="http://r-forge.r-project.org")).
It might also be nice if (but this is definitely a can of worms with
reasonable arguments in either direction) lme4 could be tolerant of mu
values on the borders of the allowed region or (???) slightly beyond it,
so that optimizations could proceed to a final answer that might not be
at those boundaries.
cheers
Ben Bolker
>
> sessionInfo()
R version 2.10.1 (2009-12-14)
i486-pc-linux-gnu
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_0.999375-32-2 Matrix_0.999375-38 lattice_0.18-3
loaded via a namespace (and not attached):
[1] grid_2.10.1
Sasha Goodman wrote:
> I just tried your toy code with crossed random effects, and got the
> following error when using glmer
>
> "Error in mer_finalize(ans) :
> mu[i] must be in the range (0,1): mu = 9.35952e-313, i = 1068030290"
>
> No model was ever fit, so the following happened:
>
> fixef(g2)
> Error: object 'g2' not found
>
> This is exactly the problem I was having before.
>
>
>
> On Fri, Mar 19, 2010 at 11:46 AM, Sasha Goodman <sashag at stanford.edu
> <mailto:sashag at stanford.edu>> wrote:
>
> Ben, thanks for getting back to me.
> The problem was with crossed random effects and cloglog.
>
>
> There was a strange error:
> "mu[i] must be in the range (0,1)"
>
> Here is the original post:
>
>
>
> ..........
> We have been using lme4 with a logit link and crossed random effects.
> It works very nicely. However, because our outcome is very rare, we
> are trying the cloglog link with lmer. The model simply does not run,
>
>
> however. The errors is "mu[i] must be in the range (0,1)". A google
> search reveals no one else has ever posted this particular problem.
>
> I'm sending the following details in case it helps the developers.
>
>
> Advice is also welcomed.
>
> Descriptives for the two variables:
> Y: [0,1], ∈ R, E=0.007238727, SE=0.08477285, the binary response
> X : [0.4318617,0.998886], ∈ R, E=0.9886799, SE=0.03572924, continuous
> in the range [0.4318617,0.998886]
>
>
>
> ## Crossed with cloglog fails
> h = lmer(Y ~ X + (1 | i) + (1 | j) + (1|t), D2 ,family =
> binomial(link="cloglog"),control = list(msVerbose = 1))
> 0: 5328.8713: 0.100868 0.0552843 0.00896829 -35.7191 31.0045
>
>
> Error in mer_finalize(ans) :
> mu[i] must be in the range (0,1): mu = 0, i = 253517704
>
> ## Crossed with logit works
> h = lmer(Y ~ X + (1 | i) + (1 | j) + (1|t), D2 ,family =
> binomial(link="logit"),control = list(msVerbose = 1))
>
>
> summary(h)
> Generalized linear mixed model fit by the Laplace approximation
> Formula: Y ~ X + (1 | i) + (1 | j) + (1 | t)
> Data: D2
> AIC BIC logLik deviance
> 1797 1843 -893.5 1787
> Random effects:
>
>
> Groups Name Variance Std.Dev.
> j (Intercept) 17.9346 4.2349
> i (Intercept) 126.9971 11.2693
> t (Intercept) 2.6644 1.6323
> Number of obs: 66310, groups: j, 253; i, 76; t, 2
>
> Fixed effects:
>
>
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) -81.38 19.29 -4.219 2.46e-05 ***
> X 60.91 19.26 3.162 0.00157 **
> ---
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
>
>
>
>
> ## Simple GLM with the same data runs perfectly
> h = glm(Y ~ X , D2 ,family = binomial(link="cloglog"))
>
> glm(formula = Y ~ X, family = binomial(link = "cloglog"),
> data = D2)
>
>
>
> Deviance Residuals:
> Min 1Q Median 3Q Max
> -0.1316 -0.1271 -0.1270 -0.1212 3.9149
>
> Coefficients:
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) -35.719 6.833 -5.228 1.72e-07 ***
>
>
> X 31.004 6.869 4.514 6.37e-06 ***
> ---
> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>
> (Dispersion parameter for binomial family taken to be 1)
>
> Null deviance: 5687.7 on 66309 degrees of freedom
>
>
> Residual deviance: 5643.9 on 66308 degrees of freedom
> AIC: 5647.9
>
> Number of Fisher Scoring iterations: 10
>
>
>
> On Fri, Mar 19, 2010 at 11:03 AM, Ben Bolker <bolker at ufl.edu
> <mailto:bolker at ufl.edu>> wrote:
>
> Hmmm. Can you remind me what didn't work?
> Examples below give at least vaguely plausible answers ...
>
> cheers
> Ben Bolker
>
> set.seed(1001)
>
> N <- 20
> n <- 100
> x <- runif(N*n)
> f <- factor(rep(LETTERS[1:N],each=n))
> dat <- data.frame(x,f)
>
> beta <- c(1,3)
> alpha <- rnorm(N,sd=0.5)
>
> mm <- cbind(model.matrix(~x,data=dat),
> model.matrix(~f-1,data=dat))
>
> eta <- mm %*% c(beta,alpha)
> y <- rbinom(N*n,prob=1-exp(-exp(eta)),size=1)
> dat <- data.frame(dat,y)
>
> library(lme4)
> (g1 <- glmer(y~x+(1|f),data=dat,
> family=binomial(link="cloglog")))
>
> fixef(g1) ## matches (1,3) reasonably well
> VarCorr(g1) ## matches sd=0.5 reasonably well
>
> ## now try crossed random effects
> set.seed(1001) ## reset
> N1 <- 10
> N2 <- 10
> n <- 20
> ntot <- n*N1*N2
> dat <- expand.grid(f1=LETTERS[1:N1],f2=letters[1:N2],rep=1:n)
> dat <- data.frame(dat,x=runif(ntot))
> alpha1 <- rnorm(N1,sd=0.5)
> alpha2 <- rnorm(N2,sd=1)
>
> mm <- cbind(model.matrix(~x,data=dat),
> model.matrix(~f1-1,data=dat),
> model.matrix(~f2-1,data=dat))
>
> eta <- mm %*% c(beta,alpha1,alpha2)
> y <- rbinom(ntot,prob=1-exp(-exp(eta)),size=1)
>
> (g2 <- glmer(y~x+(1|f1)+(1|f2),data=dat,
> family=binomial(link="cloglog")))
>
> ## warning about false convergence
> fixef(g2) ## still OK.
> VarCorr(g2) ## reasonable but ??
>
>
>
> Sasha Goodman wrote:
> > Bump. So what was your solution to the cloglog issue? Has it
> been fixed
> > in lme4?
> >
> > Other people keep asking me about this issue.
> >
> > On Tue, Nov 3, 2009 at 4:48 PM, Sasha Goodman
> <sashag at stanford.edu <mailto:sashag at stanford.edu>
> > <mailto:sashag at stanford.edu <mailto:sashag at stanford.edu>>> wrote:
> >
> > So, I want to use the cloglog link. Is your solution
> working now?
> > Does it worked with crossed random effects?
> >
> >
> > On Thu, Oct 1, 2009 at 1:46 PM, Sasha Goodman
> <sashag at stanford.edu <mailto:sashag at stanford.edu>
> > <mailto:sashag at stanford.edu <mailto:sashag at stanford.edu>>>
> wrote:
> >
> > Very interested! Never came up with a solution. Even
> retried
> > with a recent version of lme4 recently.
> >
> >
> > On Thu, Oct 1, 2009 at 2:29 PM, Ben Bolker
> <bolker at ufl.edu <mailto:bolker at ufl.edu>
> > <mailto:bolker at ufl.edu <mailto:bolker at ufl.edu>>> wrote:
> >
> >
> > I came across your post to the r-sig-mixed list
> from May.
> > Did you ever resolve your problem? I may have a
> solution
> > (of sorts)
> > if you're still interested.
> >
> > Ben Bolker
> >
> > --
> > Ben Bolker
> > Associate professor, Biology Dep't, Univ. of Florida
> > bolker at ufl.edu <mailto:bolker at ufl.edu>
> <mailto:bolker at ufl.edu <mailto:bolker at ufl.edu>> /
> > www.zoology.ufl.edu/bolker
> <http://www.zoology.ufl.edu/bolker>
> <http://www.zoology.ufl.edu/bolker>
> > GPG key:
> www.zoology.ufl.edu/bolker/benbolker-publickey.asc
> <http://www.zoology.ufl.edu/bolker/benbolker-publickey.asc>
> >
> <http://www.zoology.ufl.edu/bolker/benbolker-publickey.asc>
> >
> >
> >
> >
> > --
> > Sasha Goodman
> > Doctoral Candidate, Organizational Behavior
> >
> >
> >
> >
> > --
> > Sasha Goodman
> > Doctoral Candidate, Organizational Behavior
> >
> >
> >
> >
> > --
> > Sasha Goodman
> > Doctoral Candidate, Organizational Behavior
>
>
> --
> Ben Bolker
> Associate professor, Biology Dep't, Univ. of Florida
> bolker at ufl.edu <mailto:bolker at ufl.edu> /
> people.biology.ufl.edu/bolker <http://people.biology.ufl.edu/bolker>
> GPG key: people.biology.ufl.edu/bolker/benbolker-publickey.asc
> <http://people.biology.ufl.edu/bolker/benbolker-publickey.asc>
>
>
>
>
> --
> Sasha Goodman
> Doctoral Candidate, Organizational Behavior
>
>
>
>
> --
> Sasha Goodman
> Doctoral Candidate, Organizational Behavior
--
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bolker at ufl.edu / people.biology.ufl.edu/bolker
GPG key: people.biology.ufl.edu/bolker/benbolker-publickey.asc
More information about the R-sig-mixed-models
mailing list