R help - Oct 2009 - GLMMPQL and negbinomial: trouble with the X-axis in PREDICT

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
scale.
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)
>theta<-m0$theta
> 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
 Correlation:
        (Intr)
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"

Ahh, yes. The level made all the difference. I originally thought the
level only applied to binomial response variables, so I disregarded
it.
And I'll just add that >plot(e) below still produced an indexed x
axis, so I just combined the vectors and then was able to plot the
scaled variable properly.
> MyData<-data.frame(BERRIES = seq(from =20, to=3600, by=100))
> e<-predict(m1,MyData,type="response", level=0)
> q2<-cbind(MyData,e)
> plot(q2$BERRIES,q2$e)
Thanks for the fast responses!

~erin

Erin Latham, MGIS Candidate
Geography Dept, University of Calgary, AB, Canada