[R] GLMMPQL and negbinomial: trouble with the X-axis in PREDICT

Erin Latham els.latham at gmail.com
Mon Oct 26 20:26:47 CET 2009

I'm having some difficulty with graphing outputs of a GLM model I've
been working. I have count data for both my predictor (only 1) and
response variables, and I have pseudoreplication which I've modeled as
a random effect. The odTest() from pscl:: indicated that the negative
binomial distribution fit better than Poisson, and I then proceeded by
estimating theta from glm.nb.
My residual plot of m1 looks great, and now I'd like an output showing
the relationship between BERRIES and the respose, but I can't seem to
get "BERRIES" on the x-axis. It keeps spitting out "INDEX" which is
just all my numbered observations from 1 to 121, and therefore not to
I did try creating a new vector to plot the data, but I got an error
message (see below).
Does anyone have suggestions on where to look to solve this problem?

Erin Latham,
M.GIS Candidate, Geography
University of Calgary, AB, Canada

> m0 <- glm.nb(MANUAL ~ BERRIES + BUSHID , data=a)
> m1 <- glmmPQL(MANUAL ~ BERRIES, random = ~ 1|BUSHID, data=a, family=negative.binomial(theta))
> summary(m1)
Linear mixed-effects model fit by maximum likelihood
 Data: a
  AIC BIC logLik
   NA  NA     NA

Random effects:
 Formula: ~1 | BUSHID
        (Intercept)  Residual
StdDev:    1.081722 0.6577578

Variance function:
 Structure: fixed weights
 Formula: ~invwt
Fixed effects: MANUAL ~ BERRIES
               Value  Std.Error DF  t-value p-value
(Intercept) 7.602133 0.23105575 91 32.90172       0
BERRIES     0.001369 0.00023395 28  5.85324       0
BERRIES -0.485

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-2.14692244 -0.52715093 -0.02833408  0.48002141  2.70931877

Number of Observations: 121
Number of Groups: 30

> MyData<-data.frame(BERRIES = seq(from =20, to=3600, by=100))
> e<-predict(m1,MyData,type="response")
Error in predict.lme(object, newdata, level = level, na.action = na.action) :
  Cannot evaluate groups for desired levels on "newdata"

More information about the R-help mailing list