[R] help with glmmPQL

[Spencer Graves]

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