[R] help with glmmPQL
Spencer Graves
spencer.graves at pdf.com
Fri Nov 26 20:48:25 CET 2004
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
More information about the R-help
mailing list