cloglog and lmer
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
> GPG key:
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