[R-sig-ME] glmer.nb warning and prediction problem
Adriaan de Jong
Adr|@@n@de@Jong @end|ng |rom @|u@@e
Sat Mar 5 12:29:21 CET 2022
Hi Ben,
Thank you very much for your comments and suggestions.
As for point 2.
Adding quotation marks didn't solve the problem. But from the lapply(model.frame(mod5), unique) printout (listed below) I understand that the Dec_Hour = 7.0 was not represented among the unique values, but 7.00 was. I changed the newd vector (newd<-data.frame(Years=c(1,29),Dag="34",Dec_Hour="7.00"), but this didn't solve the problem either. The Dag variable does include the value 34 so that couldn't cause the problem. So yes, an error message with more information wouldn't hurt :)
From a biological point of view, the decimal-hours ("Dec_Hour") form of time-of-day doesn't make much more sense than using just full hours ("Hour") in the random effects part. When I changed "Dec_Hour" to "Hour" in the original glmer.nb, things went smoothly (at least after the iteration limit warning). The prediction function
(newd4<-data.frame(Years=c(1,29),Dag="34",Hour="7")
pred4<-predict(mod4,newd4)).
also worked well. With and without levels in quotation marks!
As for point 1.
This is great news, because negative binomial modelling has been troublesome for me. Being a non-professional statistician, it makes me wonder, though, when and how to apply the var/mean ratio (often seen as a "law" in family selection).
Never mind, I changed to a poisson glmer (glmer(Total~ Years + (1|Dag/Dec_Hour), family=poisson)) and (obviously) received no iteration limit warning, but again, the prediction function didn't work. Similar to the above, this problem was solved by changing from decimal to Hours in the random effects part.
From this I conclude that the use of Dec_Hour (decimal hour based on a combination of hour and minutes data) is inappropriate in both glmer and glmer.nb when predictions are wanted. Is this information provided anywhere in the documentation of these functions? From my point of view, the interesting thing is that Dec_Hour worked flawless in gam and the predict.gam functions.
This leaves me in a model selection situation. Based on AIC (list below), the glmer.nb with just full hours seems to be the model of choice, especially because it allows for predictions. Does this choice make sense to you?
Cheers,
Adjan
Function and formulaAIC
glmer(Total~ Years + (1|Dag/Hour), family=poisson)4973
glmer(Total~ Years + (1|Dag/Dec_Hour), family=poisson)4920 (prediction failure)
glmer.nb(Total~ Years + (1|Dag/Dec_Hour))4922 (iteration limit warning, prediction failure)
glmer.nb(Total~ Years + (1|Dag/Hour))4914 (iterarion limit warning)
lapply(model.frame(mod5)
$Total
[1] 9 3 7 14 10 8 12 5 1 2 11 6 17 4 19 15 13 16 18 21 20 0 22
$Years
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
$Dag
[1] 33 41 42 47 48 53 59 62 66 67 68 72 77 78 79 31 32 37 38 40 43 46 54 57 61 70 35 50 51 52 56 63 65 75
[35] 36 39 44 49 58 69 73 34 60 80 45 64 55 74 76 71
$Dec_Hour
[1] 11.42 11.17 15.42 5.33 7.92 15.92 16.92 17.17 10.92 20.42 17.42 14.92 8.42 11.67 16.42 12.17 16.17
[18] 17.33 10.17 14.17 17.67 15.67 19.50 13.42 9.50 16.25 10.50 12.50 11.75 6.42 19.08 8.17 14.00 16.00
[35] 20.17 7.67 11.08 17.25 11.00 12.42 10.00 13.00 17.83 6.58 13.17 15.08 7.83 14.42 8.00 9.58 14.50
[52] 9.25 12.25 18.83 6.67 16.67 8.33 12.00 9.67 12.33 7.50 14.67 15.17 20.92 11.25 19.00 9.00 11.92
[69] 12.83 17.50 18.00 18.67 12.58 8.67 20.67 8.50 17.92 13.67 18.25 18.17 11.50 16.33 11.83 13.08 15.50
[86] 19.17 18.50 10.75 9.83 13.83 14.75 8.58 9.08 10.83 11.33 15.25 17.75 12.92 18.58 19.67 13.50 13.92
[103] 10.42 17.00 18.33 5.92 6.33 7.70 15.00 15.83 13.33 6.75 8.83 18.42 19.42 15.05 15.58 5.17 12.67
[120] 13.75 10.33 12.08 6.83 9.17 10.58 20.58 20.75 6.92 16.83 5.42 12.75 8.92 5.67 14.08 18.08 16.58
[137] 14.33 5.08 7.17 6.17 9.75 17.08 10.08 5.25 5.27 10.67 20.33 19.92 17.48 6.50 15.33 14.83 16.50
[154] 20.00 9.92 16.75 5.75 5.83 20.83 7.00 19.33 17.58 18.75 20.25 5.00 19.58 15.75 7.25 16.08 7.75
[171] 7.33 7.42 5.50 20.08 8.25 8.08 6.00 7.58 13.58 18.92 14.58 14.25 19.75 6.08 19.83 7.08 10.25
[188] 9.42 9.33 6.38 13.25 6.25 19.25 8.75
-----Original Message-----
From: Ben Bolker <bbolker using gmail.com>
Sent: den 4 mars 2022 21:02
To: Adriaan de Jong <Adriaan.de.Jong using slu.se>
Cc: r-sig-mixed-models using r-project.org
Subject: Re: [R-sig-ME] glmer.nb warning and prediction problem
(It wouldn't hurt for us to make the error message slightly more informative by listing *which* new levels were found ...)
On Fri, Mar 4, 2022 at 3:02 PM Ben Bolker <bbolker using gmail.com> wrote:
>
> (1) tl;dr don't worry about it. It's very, very likely that you have
> data whose conditional distribution is indistinguishable from a
> Poisson. I know that the *marginal* variance/mean ratio is large, but
> once we put in the covariates and the random effects, what's left is
> presumably equi- or underdispersed (variance ~ mean or variance <
> mean). The mean counts at the beginning of the time series are about 8
> (exp(2)), decreasing over time; your dispersion parameter is 798, >>
> mean - that means the distribution is effectively Poisson. The
> iteration limit warning comes from the initial call to MASS::glm.nb,
> which tries to determine an initial guess for the dispersion parameter
> via an iterative algorithm - it gives up after a while.
> You could try a Poisson model - my guess is that it will give
> practically indistinguishable results (and a likelihood ratio test
> probably won't reject the null hypothesis that the conditional
> distribution is Poisson).
> Or you could ignore the warning.
>
> 2. You might need to put the levels of `Dag` and `Dec_Hour` in
> quotation marks. What is `lapply(model.frame(mod5), unique)` ?
>
>
> On Fri, Mar 4, 2022 at 7:02 AM Adriaan de Jong <Adriaan.de.Jong using slu.se> wrote:
> >
> > Dear list members,
> >
> > What I’m trying to do is to model a linear trend over the years (Years = Year – 1992, 1993 was the starting year of the 29 year data series) of bird numbers (“Total”) with a var/mean ratio of 2.85. The variables “Dag” (= day of the study period starting on April 1 = 1, in this subset ranging from 31 to 80) and “Dec_Hour” (= decimal hour between 4 AM and 9 PM) are nested random effects. The subset for this analysis (Trend1) contains 993 rows (=counting trips).
> >
> > The mod5<-glmer.nb(Total~ Years + (1|Dag/Dec_Hour),data=Trend1) command under lme4 gives me the following warning message:
> >
> > In theta.ml(Y, mu, weights = object using resp$weights, limit = limit, :
> > iteration limit reached
> >
> > I presume I need to provide glmer.nb with a lmerControl string, but have no clue what it should contain.
> >
> > I guess that the Dec_Hour may be causing the iteration limit
> > problem. The model seems to converge, though, and the summary looks
> > like this
> >
> >
> > Generalized linear mixed model fit by maximum likelihood (Laplace
> > Approximation) [ glmerMod]
> > Family: Negative Binomial(798.6507) ( log )
> > Formula: Total ~ Years + (1 | Dag/Dec_Hour)
> > Data: Trend1
> >
> > AIC BIC logLik deviance df.resid
> > 4921.7 4946.2 -2455.8 4911.7 988
> >
> > Scaled residuals:
> > Min 1Q Median 3Q Max
> > -1.9784 -0.6199 -0.0885 0.4511 3.8851
> >
> > Random effects:
> > Groups Name Variance Std.Dev.
> > Dec_Hour:Dag (Intercept) 0.15182 0.3896
> > Dag (Intercept) 0.09292 0.3048
> > Number of obs: 993, groups: Dec_Hour:Dag, 945; Dag, 50
> >
> > Fixed effects:
> > Estimate Std. Error z value Pr(>|z|)
> > (Intercept) 2.09539 0.06025 34.78 <2e-16 ***
> > Years -0.04078 0.00248 -16.44 <2e-16 ***
> > ---
> > Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> >
> > Correlation of Fixed Effects:
> > (Intr)
> > Years -0.594
> >
> > This makes sense to me, but I would like to predict the 1993 and 2021 endpoint "Total" values of the yearly trend (numerically Years = 1 and 29) given the peak values of observed birds found to be Dag=34 and Dec_Hour=7. For this I tried:
> >
> > newd<-data.frame(Years=c(1,29),Dag=34,Dec_Hour=7.0)
> > pred<-predict(mod5,newd)
> >
> > This rendered the message:
> > Error in levelfun(r, n, allow.new.levels = allow.new.levels) :
> > new levels detected in newdata
> >
> > Could anyone give me a hint on (a) how to avoid the iteration problem and (b) how to adjust the predict function? Thanks in advance for your help.
> >
> > Cheers,
> > Adjan
> >
> > Adriaan "Adjan" de Jong
> > Associate professor
> > Swedish University of Agricultural Sciences
> > ---
> > När du skickar e-post till SLU så innebär detta att SLU behandlar
> > dina personuppgifter. För att läsa mer om hur detta går till, klicka
> > här <https://www.slu.se/om-slu/kontakta-slu/personuppgifter/>
> > E-mailing SLU will result in SLU processing your personal data. For
> > more information on how this is done, click here
> > <https://www.slu.se/en/about-slu/contact-slu/personal-data/>
> > _______________________________________________
> > R-sig-mixed-models using r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
---
När du skickar e-post till SLU så innebär detta att SLU behandlar dina personuppgifter. För att läsa mer om hur detta går till, klicka här <https://www.slu.se/om-slu/kontakta-slu/personuppgifter/>
E-mailing SLU will result in SLU processing your personal data. For more information on how this is done, click here <https://www.slu.se/en/about-slu/contact-slu/personal-data/>
More information about the R-sig-mixed-models
mailing list