[R] help with glmmPQL
John Fox
jfox at mcmaster.ca
Sat Nov 27 03:10:19 CET 2004
Dear Doug,
I assume that in the absence of a suitable anova() method, there's nothing
wrong with testing the difference between two GLMM objects nested in their
random effects by manually comparing the log-likelihoods as printed or
returned by logLik(). Is that correct?
Regards,
John
--------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
--------------------------------
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Douglas Bates
> Sent: Friday, November 26, 2004 3:13 PM
> To: Spencer Graves
> Cc: R-help; bates at cs.wisc.edu; Jose.Pinheiro at pharma.novartis.com
> Subject: Re: [R] help with glmmPQL
>
> 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.
>
> I don't think that would be appropriate. AFAIK the logLik
> generic applied to a glmmPQL object does not return the
> likelihood of the fitted model or an approximation to that likelihood.
>
> That's why there is a difference between the values of the
> likelihood for the same model fit by glmmPQL and fit using
> the PQL method in GLMM from the lme4 package.
>
> I will add anova methods for GLMM objects whenever I manage
> to dig myself out of the current rewrite of the
> representation of linear mixed-effects models.
>
> > 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.
> > I looked briefly at the code for "anova.lme", and it
> looks to me
> > like "glmmPQL" could be added to the list of allowed options in the
> > following test roughly half way through the code:
> > if (!all(match(termsClass, c("gls", "gnls", "lm", "lmList",
> > "lme", "nlme", "nlsList", "nls"), 0))) {
> > stop(paste("Objects must inherit from classes \"gls\",
> > \"gnls\"",
> > "\"lm\",\"lmList\", \"lme\",\"nlme\",\"nlsList\", or
> > \"nls\""))
> > }
> >
> > hope this helps. spencer graves
> >
> > A.J. Rossini wrote:
> >
> >> For better or worse, it's holidays in the states. Very
> amusing for
> >> me being in a non-Thanksgiving celebrating country.
> >>
> >> In addition, it's not a problem. The complaint is valid.
> Probably
> >> no one has coded up the right solution yet for comparison.
> >>
> >> I can't recall if one would want those statistics for a binomial
> >> random effects model, but I do recall some issues with model
> >> comparison in that setting, though they are a bit dated
> (say, 2 years
> >> or so).
> >>
> >> On Fri, 26 Nov 2004 09:31:40 +0700, Andrew Criswell
> >> <r-stats at arcriswell.com> wrote:
> >>
> >>
> >>> Hello:
> >>>
> >>> Will someone PLEASE help me with this problem. This is the third
> >>> time I've posted it.
> >>>
> >>> When I appply anova() to two equations estimated using glmmPQL, I
> >>> get a complaint,
> >>>
> >>>
> >>>
> >>>> anova(fm1, fm2)
> >>>>
> >>>
> >>> Error in anova.lme(fm1, fm2) : Objects must inherit from classes
> >>> "gls", "gnls" "lm","lmList", "lme","nlme","nlsList", or "nls"
> >>>
> >>>
> >>> The two equations I estimated are these:
> >>>
> >>>
> >>>
> >>>> fm1 <- glmmPQL(choice ~ day + stereotypy,
> >>>>
> >>>
> >>> + random = ~ 1 | bear, data = learning, family =
> >>> binomial)
> >>>
> >>>
> >>>
> >>>> fm2 <- glmmPQL(choice ~ day + envir + stereotypy,
> >>>>
> >>>
> >>> + random = ~ 1 | bear, data = learning, family =
> >>> binomial)
> >>>
> >>> Individually, I get results from anova():
> >>>
> >>>
> >>>
> >>>> anova(fm1)
> >>>>
> >>>
> >>> numDF denDF F-value p-value
> >>> (Intercept) 1 2032 7.95709 0.0048
> >>> day 1 2032 213.98391 <.0001
> >>> stereotypy 1 2032 0.42810 0.5130
> >>>
> >>>
> >>>
> >>>> anova(fm2)
> >>>>
> >>>
> >>> numDF denDF F-value p-value
> >>> (Intercept) 1 2031 5.70343 0.0170
> >>> day 1 2031 213.21673 <.0001
> >>> envir 1 2031 12.50388 0.0004
> >>> stereotypy 1 2031 0.27256 0.6017
> >>>
> >>>
> >>> I did look through the archives but didn't finding
> anything relevant
> >>> to my problem.
> >>>
> >>> Hope someone can help.
> >>>
> >>> ANDREW
> >>> ____________________________
> >>> _
> >>> platform i586-mandrake-linux-gnu
> >>> arch i586
> >>> os linux-gnu
> >>> system i586, linux-gnu
> >>> status
> >>> major 2
> >>> minor 0.0
> >>> year 2004
> >>> month 10
> >>> day 04
> >>> language R
> >>>
> >>> --
> >>> Andrew R. Criswell, Ph.D.
> >>> Graduate School, Bangkok University
> >>>
> >>> ______________________________________________
> >>> R-help at stat.math.ethz.ch mailing list
> >>> https://stat.ethz.ch/mailman/listinfo/r-help
> >>> PLEASE do read the posting guide!
> >>> http://www.R-project.org/posting-guide.html
> >>>
> >>>
> >>
> >>
> >>
> >>
> >>
> >
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
More information about the R-help
mailing list