help with glmmPQL
Hi, Deepayan:
Thanks much. That's very helpful.
ANDREW: How difficult would it be for you to generate a Monte Carlo to
simulate data according to your two models, e.g., "y ~ trt" and "y ~ trt
+ I(week > 2)"? If you did that, you could then fit both models to each
set of simulated data, and compute and store logLR <-
fit2$logLik-fit1$logLik for each one. This will give you a reference
distribution, from which you can estimate both the p-value and the
statistical power of this analysis against your chosen alternatives.
If you do that, I suggest you use "GLMM" in library(lme4), not
glmmPQL. The "logLik" produced by glmmPLQ for model 2 was LESS THAN
that for model 1. If the function were maximizing the likelihood,
fit2$logLik should be greater than fit1$logLik, not less.
Of course, of you simulate both models and compute the
distribution of your favorite test statistic, you can get a p-value that
is as good as your simulation. I've done this kind of thing before, and
it should be relatively easy in R using rnorm and rbinom.
hope this helps.
spencer graves
Deepayan Sarkar wrote:
On Friday 26 November 2004 12:35, Spencer Graves wrote:
HI, DOUG & JOSE:
Is there some reason that "anova.lme" should NOT accept an
object of class "glmmPQL" in the example below? If you don't see one
either, then I suggest you consider modifying the code as described
below.
HI, ANDREW:
I couldn't find your data "learning" in my Windows installation
of R 2.0.0pat, which meant that I had to take the time to find
another example before I could get the error message you described.
I got it from modifying the example in the documentation for
"glmmPQL" as follows:
fit1 <- glmmPQL(y ~ trt, random = ~ 1 | ID,
family = binomial, data = bacteria)
fit2 <- glmmPQL(y ~ trt + I(week > 2), random = ~ 1 | ID,
family = binomial, data = bacteria)
anova(fit1, fit2)
Error in anova.lme(fit1, fit2) : Objects must inherit from
classes "gls", "gnls" "lm","lmList", "lme","nlme","nlsList", or "nls"
I then checked the class of fit1 and fit2:
class(fit1)
[1] "glmmPQL" "lme"
class(fit2)
[1] "glmmPQL" "lme"
As an experiment, I changed the class of fit1 and fit2:
class(fit1) <- "lme" class(fit2) <- "lme" anova(fit1, fit2)
Model df AIC BIC logLik Test L.Ratio p-value
fit1 1 5 1054.623 1071.592 -522.3117
fit2 2 6 1113.622 1133.984 -550.8111 1 vs 2 56.99884 <.0001
Unless someone like Doug or Jose tells us, "Don't do that", I
would use these answers.
These likelihoods are (AFAIK) NOT the likelihoods of the models fitted, they are likelihoods of an lme model that approximates it. Thus, the test may not be appropriate. Having 'anova(fit1, fit2)' silently producing an answer would certainly be misleading. E.g., lme4 produces the following:
library(lme4) data(bacteria, package = "MASS") fit1 <- GLMM(y ~ trt, random = ~ 1 | ID,
+ family = binomial, data = bacteria)
fit2 <- GLMM(y ~ trt + I(week > 2), random = ~ 1 | ID,
+ family = binomial, data = bacteria)
logLik(fit1)
`log Lik.' -104.3780 (df=5)
logLik(fit2)
`log Lik.' -97.68162 (df=6) Deepayan
Spencer Graves, PhD, Senior Development Engineer O: (408)938-4420; mobile: (408)655-4567