Multiple comparisons demo

[msmcfarlin]

Hello,

I recently wrote this comment on #84, though with it being a closed issue I'm thinking maybe it would not get noticed. Sorry, this is my first time using GitHub!

I am interested in the same type of multiple comparisons between groups noted in #84 with betta() and have a question on the demo you detailed. I apologize but I'm getting lost, specifically the step where you calculate the test statistics...

_teststats <- MASS::ginv(expm::sqrtm(A %% bt$cov %% t(A))) %% A %% betahat

What should I be inputting for "betahat"? Also, would you be able to post the adjusted p values for this demo?

Thank you for breakaway, and the other diversity analysis programs from your group! These are really amazing tools.

[ailurophilia]

Hi – thanks for your question!

Breakaway now includes a new function, betta_lincom, specifically for this purpose, so I recommend using that!

Here's an example with some toy data (taken from a response to DivNet issue 37):

set.seed(1)
estimates <- runif(n = 12, min = 1000, max = 1200)
ses <- runif(n = 12, min = 1, max = 5)

# the covariates
types <- rep(c("heart", "lung", "kidney", "liver"), each = 3)
x_default <- model.matrix(~types) # heart is baseline; since alphabetically first
fitted_betta <- betta(estimates, ses, X = x_default) # tests w heart as default

# construct estimates, 95% CIs, hypothesis test for beta_lung - beta_liver
betta_lincom(fitted_betta,
             c(0,0,-1,1))

# construct estimates, 95% CIs, hypothesis test for beta_lung - beta_kidney

betta_lincom(fitted_betta,
             c(0,-1,0,1))

There is more documentation in (the updated version of) Breakaway, including in the hypothesis testing vignette.

As for p-value adjustment, there are a variety of possibilities depending on what exactly your goal is. A variety of options are available in p.adjust in base R.

I hope this helps!

Best, David

[msmcfarlin]

Hi @ailurophilia,

I think this should address my question, though has this function been added to the version of Breakaway that can be downloaded from GitHub? I reinstalled Breakaway and tried the updated vignette but got a could not find function error for betta_lincom. I also didn't see it in the list of functions

[ailurophilia]

Hi @msmcfarlin,

The updated version of breakaway is on the master branch of the github repository. I don't run into the issue you describe if I run devtools::install_github("adw96/breakaway") and then build the vignette. Is this how you are installing breakaway?

Thanks!

David

[ailurophilia]

You can also just directly use devtools::install_github("adw96/breakaway", build_vignettes = TRUE), and then access vignettes with browseVignettes("breakaway"). Let me know if this doesn't work – thanks!

[ailurophilia]

Last thought for the moment: if you don't have devtools installed, the remotes package also includes the install_github function (and I believe this is a more streamlined package than devtools, which is a mildly huge package). So, also should work to call remotes::install_github("adw96/breakaway", build_vignettes = TRUE)!

[msmcfarlin]

Hi David,

I've tried install_github("adw96/breakaway", build_vignettes = TRUE) with both devtools and remotes but still get the error that the function betta_lincom could not be found. I do have the updated vignette with reference to the betta_lincom function though.

[msmcfarlin]

For the moment I've just copied the script for betta_lincom from the github page to try using it in R. Regarding multiple comparisons and p-value adjustments. Would you suggest running betta_lincom for all comparisons, then extract the p-values from each betta_lincom object, and adjust with p.adjust?

Something like this based on the example above? p.adjust(c(p.heart_lung, p.heart_kidney, p.heart_liver, p.lung_kidney, p.lung_liver, p.kidney_liver), method = "BH")

[ailurophilia]

Hmm – I'm glad you found a workaround, though it still doesn't seem ideal that the package isn't installing with vignettes, particularly since betta_lincom is in the namespace and I'm unable to reproduce the error you're encountering.

I'm curious to know if you run the following if you get an error (and also what version of breakaway R tells you you're working with):

remove.packages("breakaway")
remotes::install_github("adw96/breakaway", build_vignettes = TRUE)
packageVersion("breakaway")

And yes – I'm suggesting that if you want to adjust p-values, p.adjust is a reasonable way to do that. As a side note, different adjustments have different interpretations, so it may be helpful to consult the references in the function documentation for p.adjust.

[msmcfarlin]

Hi @ailurophilia,

I got breakaway to install correctly with the betta_lincom function. I had a few packages that needed to be updated and it seems that one of those prevented betta_lincom from loading successfully, though the rest of the package seemed to load. Everything is working now, thanks for your assistance!

[ailurophilia]

@msmcfarlin – awesome! So glad it's working for you now!

[msmcfarlin]

Hi @ailurophilia,

If I am interested in pairwise differences between groups, should I first run betta without an intercept? For example, changing the model matrix in the hypothesis testing vignette to model.matrix(~Day-1). Then run betta, and compare group means with betta_lincom(bt_day_fixed_id_random, linear_com=c(1,-1), signif_cutoff=0.05)?

Thanks for all your assistance!

[ailurophilia]

Hi @msmcfarlin,

I'm not totally sure what model you're fitting, but if your predictor is a categorical variable, it doesn't (or at least shouldn't) matter if you fit your model with or without an intercept. If you fit a model with an intercept, the intercept is the (estimated) mean diversity in an arbitrary reference group, and the other coefficients are estimated differences between mean diversity in each of the other groups and the reference group. So, differences of these coefficients are differences of group means since

(B_1 - B_0) - (B_2 - B_0) = B_1 - B_2

if B_0 is the (estimated) mean in the reference group and B_1 and B_2 are estimated means in two other groups.

Additionally, I'm not sure how many groups you are trying to compare, but note that the linear_com argument of betta_lincom must be the same length as the number of estimated coefficients (including intercept) in your betta model. (This argument tells the function what linear combination of estimated coefficients to compute statistics for -- e.g., c(-1,1,0) asks betta_lincom to produce point estimates, CIs, etc. for the linear combination -1B_0 + 1B_1 + 0*B_2 of estimated coefficients B_0, B_1, B_2 if these are the only coefficients betta returns.)

I hope this helps!