[R-sig-ME] how to get satterthwaite/Kenward-Roger degrees of freedom for heterogeneous-variance one-way ANOVA from an object of nlme::gls

[Maarten Jung]

The emmeans::joint_tests function can give you type III style tests.

library(nlme)
library(emmeans)

d <- ToothGrowth
d$dose <- as.factor(d$dose)

res.gls <- gls(len ~ supp * dose, data = d, weights = varIdent(form = ~ 
1 | supp * dose))

anova(res.gls, type = "marginal")

joint_tests(res.gls, mode = "satterthwaite")
joint_tests(res.gls, mode = "df.error")
joint_tests(res.gls, mode = "asymptotic")

Best,
Maarten

On 16.10.2024 21:42, Qiu, Weiliang /US via R-sig-mixed-models wrote:
> Many thanks for the info!
>
> From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> On Behalf Of Ben Bolker
> Sent: Wednesday, October 16, 2024 12:06 PM
> To: r-sig-mixed-models using r-project.org
> Subject: Re: [R-sig-ME] how to get satterthwaite/Kenward-Roger degrees of freedom for heterogeneous-variance one-way ANOVA from an object of nlme::gls
>
>
>
> emmeans() has some capability for computing Satterthwaite df for gls
> model; it will work for your example, although it wouldn't extend (or I
> don't see how to use it) to testing factors with more than two levels.
>
> https://cran.r-project.org/web/packages/emmeans/vignettes/models.html#K<https://cran.r-project.org/web/packages/emmeans/vignettes/models.html#K>
> help("emmc-functions")
>
> library(emmeans)
> contrast(emmeans(res.gls, ~ group), "trt.vs.ctrl")
>
> contrast estimate SE df t.ratio p.value
> group2 - group1 1.58 0.849 17.8 1.861 0.0794
>
> Degrees-of-freedom method: satterthwaite
>
> I don't see a way to get emmeans to estimate df matching the value
> that nlme::anova() gets (i.e. df = 18)
>
> anova(res.gls, type = "marginal")
> ## df = Inf
> contrast(emmeans(res.gls, ~ group, mode = "asymptotic"), "trt.vs.ctrl")
> ## df = 17
> contrast(emmeans(res.gls, ~ group, mode = "df.error"), "trt.vs.ctrl")
>
> Not sure how this will extend to more complex models (i.e. type-3,
> where you want to test main effects in the presence of interactions).
>
> There is a pull request adding some of this functionality to the
> pbkrtest package (but only does model comparisons, so again may be hard
> to get type-3 Anova out of it), but it looks like it was never merged/is
> now slightly out of sync.
>
> If you can get the Satterthwaite (or Kenward-Roger) df computed, you
> could hack the output of some other method to recompute p-values with
> the new df ...
>
>
>
>
> On 10/16/24 10:42, Qiu, Weiliang /US via R-sig-mixed-models wrote:
>> Greetings. I am using nlme::gls() to perform heterogeneous-variance one-way ANOVA and would like to get type III anova table. It looks like nlme:::anova.gls() does not provide Satterthwaite or Kenward-Roger degree of freedom. Please see below an example.
>>
>>
>> library(nlme)
>>
>> res.gls <- nlme::gls(extra ~ group, data = sleep, weights = varIdent(form = ~1|group))
>>
>> nlme:::anova.gls(res.gls, type = "marginal")
>>
>>
>> Could you suggest how to get type III anova table with Satterthwaite/Kenward-Roger degree of freedom from an object of nlme::gls()? Thanks!
>>
>> Best regards,
>>
>> Weiliang
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models using r-project.org<mailto:R-sig-mixed-models using r-project.org> mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models<https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
> --
> Dr. Benjamin Bolker
> Professor, Mathematics & Statistics and Biology, McMaster University
> Director, School of Computational Science and Engineering
> * E-mail is sent at my convenience; I don't expect replies outside of
> working hours.
>
> _______________________________________________
> R-sig-mixed-models using r-project.org<mailto:R-sig-mixed-models using r-project.org> mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models<https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models