[R-sig-ME] cloglog and lmer
Luca Borger
lborger at uoguelph.ca
Fri Mar 19 20:46:32 CET 2010
Hello,
in case this is useful, I tried it out (exact same code, only added
"verbose=TRUE") and got the same error message as Sasha (using the latest
released lme4 version on CRAN, I believe):
####################################################################
> library(lme4)
> (g1 <- glmer(y~x+(1|f),data=dat,
+ family=binomial(link="cloglog"),verbose=TRUE))
0: 267.35545: 0.163299 0.728061 2.49517
Error in mer_finalize(ans) :
mu[i] must be in the range (0,1): mu = 5.70618e-316, i = 115478304
>
> (g2 <- glmer(y~x+(1|f1)+(1|f2),data=dat,
+ family=binomial(link="cloglog"),verbose=TRUE))
0: 693.88107: 0.115470 0.115470 0.364870 1.73946
Error in mer_finalize(ans) :
mu[i] must be in the range (0,1): mu = 9.34227e-313, i = 115155608
>
> sessionInfo()
R version 2.10.1 (2009-12-14)
i386-pc-mingw32
locale:
[1] LC_COLLATE=English_United Kingdom.1252 LC_CTYPE=English_United
Kingdom.1252 LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C LC_TIME=English_United
Kingdom.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lme4_0.999375-32 Matrix_0.999375-37 lattice_0.18-3
loaded via a namespace (and not attached):
[1] grid_2.10.1
>
####################################################################
Cheers,
Luca
---------------------------
Luca Börger, PhD
Postdoctoral Research Fellow
Department of Integrative Biology
University of Guelph
Guelph, Ontario, Canada N1G 2W1
----- Original Message -----
From: "Ben Bolker" <bolker at ufl.edu>
To: "Sasha Goodman" <sashag at stanford.edu>; "R Mixed Models"
<r-sig-mixed-models at r-project.org>
Sent: Friday, March 19, 2010 3:25 PM
Subject: Re: [R-sig-ME] 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
>> <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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list