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)+ family = binomial, data = bacteria)
data(bacteria, package = "MASS")
fit1 <- GLMM(y ~ trt, random = ~ 1 | ID,
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
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
