[FEATURE REQUEST] testing linear combinations

[adw96]

One of the great things about consulting is that you realize what people need... and don't yet have.

I was chatting with someone today who would love to test whether a linear combination of two radEmu parameters equals zero. Specifically, they had covariate information for each month, and want to know what taxa are the most abundant in April compared to March, May compared to April, and so on. I think of this as wanting to test $\beta_{j, march} - \beta_{j, april} = 0$ and so on.

Of course, with some elaborate design matrix manipulation, this can be done (at least, for categorical contrasts) by repeatedly modifying the baseline category and refitting the model. But, it shouldn't be necessary, because the model already has all the information it needs... for a Wald test, at least.

This is a feature request for testing linear combinations of radEmu parameters. Wald tests are top priority, score tests second priority.

Here is an example with data that we could play with, of course there are billions of others.

library(tidyverse)
library(radEmu)
data("wirbel_sample")
data("wirbel_otu")
wirbel_sample
mOTU_names <- colnames(wirbel_otu)
set.seed(1)
restricted_mOTU_names <- mOTU_names %>% sample(80)
restricted_samples <- which(wirbel_otu[, restricted_mOTU_names] %>% apply(1, sum) > 0) %>% sample(100)

all_fit <- emuFit(formula = ~ Group + Study + Age_spline.1 + Age_spline.2,
                  data = wirbel_sample[restricted_samples, ], 
                  Y = wirbel_otu[restricted_samples, restricted_mOTU_names],
                  run_score_tests = FALSE)

# Does StudyFR-CRC - StudyFR-CRC = 0 for Bacteroides fragilis [ref_mOTU_v2_0286]?

rigr has a nice format for doing this:

mod_combo <- c(0, -10, 1)
lincom(mod_rigr, mod_combo, null.hypoth = -1)

and I would love to have something similar. I think it would be reasonable to assume for now that only hypotheses for the same taxa will be tested, so mod_combo can remain a vector (and doesn't have to have taxon indices eg as rows).

@svteichman would you mind looking into this? (after you're back!!!!) For the Wald I believe it's just a matter of a new $h(\cdot)$ in the relevant test statistics:

Lets start here and once that's implemented we can discuss score tests. We may have to ask @ailurophilia for some advice fitting under the null... but maybe parametrization invariance will help us?! Some thinking might save some computation.

@ailurophilia just checking you haven't done this before / have code already? I know how much you love testing a good linear combination too, so just checking ❤

[svteichman]

Sounds good! I agree that robust wald tests should be pretty straightforward, and then we'll have to think harder about the score tests.

Is the hypothesis that we're testing (in your example) $\beta_{j,march} - g_{march}(\beta_{march}) - (\beta_{k,april} - g_{april}(\beta_{april})) = 0$?

[adw96]

such a good point -- and great catch, thanks @svteichman . I think I now just live and breathe the constraint 😹 . Indeed the hypothesis that we would be testing is $$\beta_{j, March} - g_{March}(\beta_{\cdot, March}) - (\beta_{j, April} - g_{April}(\beta_{\cdot, April})) = 0$$

(= as you wrote but with $j = k$)