CODEWITHSUNDEEP
rlme4
hgecol
hgecol

Reputation: 21

Using glmer.nb(), the error message:(maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate is returned

When using glmer.nb, we just get error message

> glm1 <- glmer.nb(Jul ~ scale(I7)+ Maylg+(1|Year), data=bph.df)

Error: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate In addition: Warning message: In theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached

Who can help me? Thanks very much!

My data listed below.

   Year    Jul A7    Maylg L7b
331 1978   1948  6 1.322219   4
343 1979   8140 32 2.678518   2
355 1980 106896 26 2.267172   2
367 1981  36227 25 4.028205   2
379 1982  19085 18 2.752816   2
391 1983  26010 32 2.086360   3
403 1984   1959  1 2.506505   4
415 1985   8025 18 2.656098   0
427 1986   9780 20 1.939519   0
439 1987  48235 29 4.093912   1
451 1988   7473 30 2.974972   2
463 1989   2850 25 2.107210   2
475 1990  10555 18 2.557507   3
487 1991  70217 30 4.843563   0
499 1992   2350 31 1.886491   2
511 1993   3363 32 2.956649   4
523 1994   5140 37 1.934498   4
535 1995  14210 36 2.492760   1
547 1996   3644 27 1.886491   1
559 1997   9828 29 1.653213   1
571 1998   3119 41 2.535294   4
583 1999   5382 10 2.472756   3
595 2000    690  5 1.886491   2
607 2001    871 13       NA   2
619 2002  12394 27 0.845098   5
631 2003   4473 36 1.342423   2

Upvotes: 2

Views: 4758

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226577

You're going to have a lot of problems with this data set, among other things, because you have an observation-level random effect (you only have one data point per Year) and are trying to fit a negative binomial model. That essentially means you're trying to fit the overdispersion in two different ways at the same time.

If you fit the Poisson model, you can see that the results are strongly underdispersed (for a Poisson model, the residual deviance should be approximately equal to the residual degrees of freedom).

library("lme4")
glm0 <- glmer(Jul ~ scale(A7)+ Maylg+(1|Year), data=bph.df,
              family="poisson")
print(glm0)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: poisson  ( log )
Formula: Jul ~ scale(A7) + Maylg + (1 | Year)
   Data: bph.df
      AIC       BIC    logLik  deviance  df.resid 
 526.4904  531.3659 -259.2452  518.4904        21 
Random effects:
 Groups Name        Std.Dev.
 Year   (Intercept) 0.9555  
Number of obs: 25, groups:  Year, 25
Fixed Effects:
(Intercept)    scale(A7)        Maylg  
     7.3471       0.3363       0.6732  
deviance(glm0)/df.residual(glm0)
## [1] 0.0003479596

Or alternatively:

library("aods3")
gof(glm0)
## D  = 0.0073, df = 21, P(>D) = 1
## X2 = 0.0073, df = 21, P(>X2) = 1

glmmADMB does manage to fit it, but I don't know how far I would trust the results (the dispersion parameter is very large, indicating that the model has basically converged to a Poisson distribution anyway).

bph.df <- na.omit(transform(bph.df,Year=factor(Year)))
glmmadmb(Jul ~ scale(A7)+ Maylg+(1|Year), data=bph.df,
         family="nbinom")
GLMM's in R powered by AD Model Builder:

  Family: nbinom 
  alpha = 403.43 
  link = log 

Fixed effects:
  Log-likelihood: -259.25 
  AIC: 528.5 
  Formula: Jul ~ scale(A7) + Maylg + (1 | Year) 
(Intercept)   scale(A7)       Maylg 
  7.3628472   0.3348105   0.6731953 

Random effects:
Structure: Diagonal matrix
Group=Year
            Variance StdDev
(Intercept)   0.9105 0.9542

Number of observations: total=25, Year=25

The results are essentially identical to the Poisson model from lme4 above.

Upvotes: 2

Related Questions