[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
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"
More information about the R-help
mailing list