[R-sig-ME] bootstrap confidence interval and p-value for treatment vs control comparisons

Thu Nov 14 22:28:46 CET 2024

Dear Ben, all

Something went wrong when I replied to Ben's answer, so I am hereby reposting our two last mails that weren't sent to the entire mailing list along with (below) my new reply.
My apologies for this messy approach, I don't know how else to remedy this. Sorry!

My newest reply:
Dear Ben

Alright, thank you very much for sharing your opinion!
I agree that 88 observations is an unfortunate low number...
We have exactly 4 land use types.

Thank you again!

All the best
Lieke

Ben wrote:
> I would probably choose option 1.  With only 88 observations, I
> would take *any* conclusions with a big grain of salt (how many land use
> types do you have?); according to Frank Harrell's rules of thumb, you
> can't even do a reasonable job estimating parameters if you have more
> than about 4 land use types (which would be few enough that I wouldn't
> bother with multiple comparisons corrections anyway).
>
>    Or I might just show everything unadjusted and take the whole thing
> with a slightly bigger grain of salt ...
>
>   While I wouldn't disagree with someone who figured out how to do both
> MC-adjusted p-values and CIs (using some defensible approach), it also
> wouldn't bother me very much if unadjusted results were presented.
>
>   How many land use types *do* you have?
>
>   cheers
>    Ben Bolker

On 11/14/24 12:54, Lieke Moereels wrote:
> Dear Ben
>
> Thank you very much for your quick reply and your help!
> I really appreciate it a lot.
>
> I see, I didn't think of this issue with the Dunnett adjustment's use
> of the residual degrees of freedom. Thank you for pointing that out!
>
> I thought I'd best adjust the confidence intervals for the multiple
> comparisons as well, and not just the p-values (in this case
> afterwards), but perhaps you wouldn't think so?
> It seems like there is still some discussion on these multiplicity
> adjusted CIs, but it does seem a bit odd to report adjusted p-values but
> unadjusted CIs?
>
> This idea was among others discussed here: Benjamini, Y., and D.
> Yekutieli. “False Discovery Rate-Adjusted Multiple Confidence Intervals
> for Selected Parameters <https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.tandfonline.com%2Fdoi%2F&data=05%7C02%7CLieke.Moereels%40UGent.be%7Cb29f47a11366421b25a708dd04e2c8eb%7Cd7811cdeecef496c8f91a1786241b99c%7C1%7C0%7C638672094274001449%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=0Lly7dHi0RV%2FLZzbPA7gV8y4R3yA%2B2pxVHWrv%2F2m%2BUI%3D&reserved=0
> pdf/10.1198/016214504000001907>.” /Journal of the American Statistical
> Association/ 100, no. 469 (2005): 71–81.
> However, if using their approach to control the False Coverage Rate, I
> would in practice not perform any adjustment I believe, since I am in
> fact reporting the CIs of all 4 parameter estimates (intercept and 3
> beta-parameters), not just selecting some of them.
>
> Unfortunately, there seem to not be many options to adjust constructed
> CIs for multiplicity, besides using Bonferroni or Sidak adjustments that
> will probably be considered too conservative or the Dunnett adjustment
> that wouldn't be an option here because of what you pointed out. (I
> would of course use the same adjustment method for the p-values and the
> confidence intervals.)
>
> I'm sorry if this is bothersome, but if I could just ask this, I would
> be really happy to hear what (if hopefully any) you yourself  would
> consider the most appropriate approach out of the following options (in
> my case of a sample size of at least 88 observations and the model Y ~
> landuse + (1|location)):
>
>  1.
>     report Holm-adjusted p-values obtained from parametric bootstrap CIs
>     together with these unadjusted bootstrap CIs
>  2.
>     report Sidak- or Bonferroni-adjusted p-values obtained from
>     parametric bootstrap CIs together with the corresponding Sidak- or
>     Bonferroni-adjusted bootstrap CIs
>  3.
>     report Dunnett-like-adjusted p-values obtained from Wald tests
>     together with the corresponding multiplicity corrected Wald tests
>     (emmeans-approach; Wald-intervals not much more narrow than
>     bootstrap intervals in most cases)
>
> Thank you again a lot, apologies for all the questions and have a great
> autumn day
> Lieke

________________________________
From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> on behalf of Ben Bolker <bbolker using gmail.com>
Sent: Wednesday, November 13, 2024 4:21 PM
To: r-sig-mixed-models using r-project.org <r-sig-mixed-models using r-project.org>
Subject: Re: [R-sig-ME] bootstrap confidence interval and p-value for treatment vs control comparisons

   I think the "trt.vs.ctrl" part will be relatively easy, because this
is what you get from the default treatment contrasts in R anyway (that
is, if `landuse` is a factor, then the beta (fixed-effect) parameters
corresponding to `landuse` will be equal to the contrasts between each
level of landuse and the reference level.

   So

bb <- bootMer(fitted_model,
    FUN = function(x) {
       ff <- fixef(f)
       ff[startsWith(names(ff), "landuse")]
    },
    nsim = 1000  ## or whatever
}

and confint(bb) will give you the confidence intervals.
bb$t will give you the ensemble of bootstrap values

You could get the two-tailed p-values something like this:

boot_pval <- function(x) { 2*min(mean(x<0), mean(x>0)) }
apply(bb$t, 2, boot_pval)

   However, I think a Dunnett adjustment will be harder.  Based on this
<https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstats.stackexchange.com%2Fquestions%2F631129%2Fmulticomp-package-and-emmeans-package-produce-different-adjust-pvalues-for-dunne&data=05%7C02%7CLieke.Moereels%40UGent.be%7Ca41701cf99c54430abf308dd03f6e5f0%7Cd7811cdeecef496c8f91a1786241b99c%7C1%7C0%7C638671081155038608%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=ydOeNey4swzgaO0feSYadgXZD4S6oMRtVAIC1jHXlNw%3D&reserved=0<https://stats.stackexchange.com/questions/631129/multicomp-package-and-emmeans-package-produce-different-adjust-pvalues-for-dunne>>,
the Dunnett adjustment relies the t distribution (or multivariate t). In
particular, if you look at
<https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fen.wikipedia.org%2Fwiki%2FDunnett%2527s_test&data=05%7C02%7CLieke.Moereels%40UGent.be%7Ca41701cf99c54430abf308dd03f6e5f0%7Cd7811cdeecef496c8f91a1786241b99c%7C1%7C0%7C638671081155061334%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=dlNvVEJ5rakYyhjZ0IhN9ZoCZwfeal52eF%2FtWv8%2BsRU%3D&reserved=0<https://en.wikipedia.org/wiki/Dunnett%27s_test>>, you can see that step
5 includes knowing the residual degrees of freedom, which are hard to
define for (G)LMMs ...

   You could choose one of the other methods implemented in base R's
p.adjust() ...

   If you do want to apply Dunnett's correction on bootstrap replicates,
you might need to ask on Cross Validated
(https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstats.stackexchange.com%2F&data=05%7C02%7CLieke.Moereels%40UGent.be%7Ca41701cf99c54430abf308dd03f6e5f0%7Cd7811cdeecef496c8f91a1786241b99c%7C1%7C0%7C638671081155076909%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=Id6%2BDKtfLna1iUqjZkf5roWZ8uFTc%2FQKO9u%2F4sd62wo%3D&reserved=0<https://stats.stackexchange.com/>), although I don't know if you'll get a
satisfactory answer there either ...

   good luck,
    Ben Bolker

On 11/13/24 07:58, Lieke Moereels via R-sig-mixed-models wrote:
> Dear all
> I hope you�re doing well. I posted a similar question earlier but didn�t get a reply and because I�m under a bit of pressure to advance on this matter, I hope it�s fine I ask this modified, shorter question on the same topic.
> If I have a simple GLMM specified as: diversity ~ land use type + (1|location) where land use type has 4 levels (let�s say 1 control and 3 treatments) and location at least 8.
> How could I obtain parametric bootstrap confidence intervals and associated p-values for the difference in mean diversity between each of the treatments and the control, that are corrected for these multiple comparisons by a Dunnett correction?
> (So I'm basically looking for a parametric bootstrapping alternative to emmeans' "confint(contrast(emmeans(fitted_model, ~ landuse), "trt.vs.ctrl", ref = 1))")
>
> Any help would be greatly appreciated!
>
> Have a very nice day!
> Lieke
>
>        [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=05%7C02%7CLieke.Moereels%40UGent.be%7Ca41701cf99c54430abf308dd03f6e5f0%7Cd7811cdeecef496c8f91a1786241b99c%7C1%7C0%7C638671081155092051%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=9a5Zr2kHSV8wFr6ZHUSat52HOKm%2FyoRvQP48d3BIwMk%3D&reserved=0<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 mailing list
https://eur03.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=05%7C02%7CLieke.Moereels%40UGent.be%7Ca41701cf99c54430abf308dd03f6e5f0%7Cd7811cdeecef496c8f91a1786241b99c%7C1%7C0%7C638671081155107178%7CUnknown%7CTWFpbGZsb3d8eyJFbXB0eU1hcGkiOnRydWUsIlYiOiIwLjAuMDAwMCIsIlAiOiJXaW4zMiIsIkFOIjoiTWFpbCIsIldUIjoyfQ%3D%3D%7C0%7C%7C%7C&sdata=ys9GuIEtxEM6RYJVrqsfaXw7y%2FE%2FfUq7iqGWBa%2BVQg4%3D&reserved=0<https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models>

	[[alternative HTML version deleted]]