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

[Erin Latham]

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

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

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

[Erin Latham]

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"

__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.