Support for multiarm ancova

[gowerc]

Closes #145

Hey @danielinteractive / @tobiasmuetze,

This contains my initial code proposal for implementing multiarm ancova support (note that it is currently missing documentation and tests).

Some key points:

• To preserve backwards compatibility I've left the original ancova function unchanged. To support multiple groups a new naming scheme is required for the parameters which would break any existing code which is looking for the current names like trt / lsm_rf.

• That is, if users want to use the new multi-group ancova they will have to explicitly set fun=ancova_m_groups in the call to analyse()

• For the new naming convention I've just used the level numbers e.g. lsm_L0 = least squared means for the reference, lsm_L1 = least squared means for the first offset / coefficient. I've also used a more explicit trt_L1_L0 to mean the coefficient of the first offset from the reference.

• As I think we agreed I'm currently only extracting the treatment effects against the reference e.g. L1 - L0 & L2 - L0. I'm not calculating / extracting all pairwise combinations.

• For the naming convention I didn't want to use the user-provided group/level values as that can lead to a ton of special edge cases with whitespaces / special characters that are hard to predict and make value-extraction very error prone. Instead I opted for just renaming the group variable to rbmiGroup and re-leveling it to be L0, L1, etc. Only complication here was that the formula provided in theory can have interaction terms so had to create a function to recursively traverse the AST updating any references of group to rbmiGroup. I will need to put in extensive unit testing to make sure theres no edge cases here with things like I() functions.

If there any no main objections here I'll start trying to add docs and tests.

[danielinteractive]

Hi @gowerc , cool, thanks for working on this!

Sorry if I had a look too late - my only initial high level comment would be that I don't really "like" that we would split out ANCOVA with multiple arms from the version with two arms. I would hope that we could extend the current naming scheme from ancova also to multiple arms. This would make maintenance as well as usage easier long term.

As a possible inspiration, please see https://github.com/johnsonandjohnson/junco/blob/main/R/ancova_rbmi.R where I implemented an ancova function with multiple arms (not focused on backwards compatibility for the naming scheme).

I think in the "worst" case one could just use a condition in the ancova function and have a different naming schedule for 2 arms vs. multiple arms.

What do you think?

[gowerc]

Hi @gowerc , cool, thanks for working on this!

Sorry if I had a look too late - my only initial high level comment would be that I don't really "like" that we would split out ANCOVA with multiple arms from the version with two arms. I would hope that we could extend the current naming scheme from ancova also to multiple arms. This would make maintenance as well as usage easier long term.

As a possible inspiration, please see https://github.com/johnsonandjohnson/junco/blob/main/R/ancova_rbmi.R where I implemented an ancova function with multiple arms (not focused on backwards compatibility for the naming scheme).

I think in the "worst" case one could just use a condition in the ancova function and have a different naming schedule for 2 arms vs. multiple arms.

What do you think?

Interestingly despite what their documentation says I don't think they actually support more than 2 arms, in particular they have near the top of the function checkmate::assert_factor(data[[group]], n.levels = 2L) which would error if there were more than 2 arms.

I am very conflicted here, I agree with everything you said and I am also not happy with splitting it out but that being said I don't see an easy way of merging them without breaking backwards compatibility as the naming scheme has to be different in order for it to make sense we currently have:

lsm_alt_visit_1
lsm_ref_visit_1
trt_visit_1

In hindsight these are already bad names that confuse users. If we were to extend this whilst maintaining backwards compatibility we'd end up with something like:

lsm_alt_visit_1
lsm_alt2_visit_1
lsm_ref_visit_1
trt_visit_1
trt_alt2_visit_1

Which I would argue is even worse.

Alternatively we could just have some code that dispatches to the different ANCOVA function based on the number of levels in group variable but this makes the documentation / explanation clunky e.g.

"If you have 2 groups then it will be formatted like this, O but if you have more than 2 it would be formatted like this even though half the values are the same"

Which is why I settled on just a separate function.

One option though to minimise maintenance is that we could deprecate anocva with a note that it will be removed in the future and replaced with ancova_m_groups then in 6 months to a years time we could fully remove it, we could also create an alias of ancova = ancova_m_groups after removing so we don't have such an awkward name.

(Despite my tone I am not confident on the best path forward here as all options appear bad to me so please do challenge if you still disagree)

[danielinteractive]

Interestingly despite what their documentation says I don't think they actually support more than 2 arms, in particular they have near the top of the function checkmate::assert_factor(data[[group]], n.levels = 2L) which would error if there were more than 2 arms.

Hmm... ok weird. But I am pretty sure we wanted this to work with multiple arms, that is also what the other code suggests... I will need to check in September when I am back.

I am very conflicted here, I agree with everything you said and I am also not happy with splitting it out but that being said I don't see an easy way of merging them without breaking backwards compatibility as the naming scheme has to be different in order for it to make sense we currently have:

lsm_alt_visit_1
lsm_ref_visit_1
trt_visit_1

Yeah I understand the need for backwards compatibility. But I wonder if an alternative solution here could be to have this old naming scheme used for 2 arms, and the new naming scheme used for more than 2 arms?

In hindsight these are already bad names that confuse users. If we were to extend this whilst maintaining backwards compatibility we'd end up with something like:

lsm_alt_visit_1
lsm_alt2_visit_1
lsm_ref_visit_1
trt_visit_1
trt_alt2_visit_1

Which I would argue is even worse.

Personally I think that would also be acceptable. The user will not see these names too much anyway if they just use the rbmi functions.

Alternatively we could just have some code that dispatches to the different ANCOVA function based on the number of levels in group variable but this makes the documentation / explanation clunky e.g.

Yeah I would not have two functions necessarily but just two naming schemes. Or just go with the naming scheme just described which I think is fine.

One option though to minimise maintenance is that we could deprecate anocva with a note that it will be removed in the future and replaced with ancova_m_groups then in 6 months to a years time we could fully remove it, we could also create an alias of ancova = ancova_m_groups after removing so we don't have such an awkward name.

(Despite my tone I am not confident on the best path forward here as all options appear bad to me so please do challenge if you still disagree)

Personally I would just go with the naming scheme you suggested, I think that is reasonable. Just for the sake of a better naming scheme I would not go the route of two functions or function deprecation. Just my 2 cents 😄

[tobiasmuetze]

Hi @gowerc / @danielinteractive,
I finally got around to looking at this.

I'd prefer to have a single ancova() function.

As for the naming convention, I agree with Daniel that the following

lsm_alt_visit_1
lsm_alt2_visit_1
lsm_ref_visit_1
trt_visit_1
trt_alt2_visit_1

is acceptable.

I am wondering if long term, it would be better to allow the user to match the group levels to ref, alt, etc.

[wwojciech]

Hi @gowerc @tobiasmuetze @danielinteractive,

This functionality would be really helpful for the project I'm currently working on.

Is there any chance this PR might get merged?

Thanks!
Wojtek

[danielinteractive]

@wwojciech I think we were not yet fully aligned on the right API here for the function, and @gowerc is still on leave at the moment... @luwidmer what is your call on this one from maintainer perspective?

[wwojciech]

Hi @danielinteractive @luwidmer,

To make downstream parsing easier and more robust, have you considered exposing the components currently concatenated into the parameter column as separate columns in the output of analyse() -> pool()?

For example:

  estimate_type  group_var group_level_1 group_level_2 visit    est
  <chr>          <chr>     <chr>         <chr>         <chr>    <dbl> 
1 lsm            TRT01P    Placebo       NA            Week16   55.8
2 lsm            TRT01P    A             NA            Week16   64.2
3 lsm            TRT01P    B             NA            Week16   66.7   ......
4 contrast       TRT01P    A             Placebo       Week16   8.4
5 contrast       TRT01P    B             Placebo       Week16   10.9 

for var$group having 3 levels: "Placebo", "A", "B" (in this order).

[wwojciech]

Hi @gowerc , @luwidmer - can you please update on this PR?

[luwidmer]

@tobiasmuetze, from your experience implementing this ad-hoc, what are your thoughts on what this should look like in rbmi?

[tobiasmuetze]

I am generally fine with the idea that ref refers to the first factor level of vars$group, alt1 to the second factor level, etc. However, I think the output of pool() should be more explicit about the underlying study arms and the corresponding comparisons. In that respect, I agree with @wwojciech’s proposal, but I would suggest extending the pool() output while retaining the current parameter variable for backwards compatibility.

I would also suggest allowing users to specify the contrasts of interest, since there may be cases where two alt levels need to be compared directly. The default could remain all comparisons of each alt level against the ref level.

[wwojciech]

Thanks, @tobiasmuetze.

I have a positive feeling about keeping parameter column for backward compatibility. In that case, the parameter column in the pool() output could look like this:

parameter
<chr>
lsm_L1_Week16
lsm_L2_Week16
lsm_L3_Week16
contr_L1_L3_Week16
contr_L2_L3_Week16
contr_L2_L1_Week16

as a result of the following settings:

1. var$group has 3 levels: "A", "B", and "Placebo" (in this order).
2. The user requested the following contrasts:

• "A" - "Placebo",
• "B" - "Placebo",
• "B" - "A".

One additional suggestion: if possible, I would propose numbering the levels from L1 rather than L0. This would be more consistent with R conventions, where (for example, in lm()) the reference level is simply the first factor level rather than "level 0." I think starting with L1 would therefore be more intuitive for R users.

[tobiasmuetze]

One additional suggestion: if possible, I would propose numbering the levels from L1 rather than L0. This would be more consistent with R conventions, where (for example, in lm()) the reference level is simply the first factor level rather than "level 0." I think starting with L1 would therefore be more intuitive for R users.

I would agree in principle if we were to start building this from scratch. However, giving that {rbmi} is used in a regulated environment for decision making, I'd prefer to prioritize backward compatibility in this case and the usage of ref, alt, alt2, etc., and keeping the parameter structure the same.

Any suggestions how to specify the contrast would be appreciated.

[wwojciech]

One additional suggestion: if possible, I would propose numbering the levels from L1 rather than L0. This would be more consistent with R conventions, where (for example, in lm()) the reference level is simply the first factor level rather than "level 0." I think starting with L1 would therefore be more intuitive for R users.

I would agree in principle if we were to start building this from scratch. However, giving that {rbmi} is used in a regulated environment for decision making, I'd prefer to prioritize backward compatibility in this case and the usage of ref, alt, alt2, etc., and keeping the parameter structure the same.

Any suggestions how to specify the contrast would be appreciated.

Sure, thank you @tobiasmuetze.

Then, there must either be a clear definition of the relationship between ref, alt, et2, etc., and the levels of group, or the output table of pool() should follow the structure I proposed above (of course, with the parameter column retained for backward compatibility). Otherwise, group-level (label) identification may become error-sensitive.

Constrast specification

One approach that comes to mind naturally would be to specify contrasts through an additional group_contrasts argument in set_vars(). This could be a list of pairs of levels of group, where the first element in each pair is the minuend. For example:

group_contrasts = list(
  c("A", "Placebo"),
  c("B", "Placebo")
)

This would represent the contrasts A - Placebo and B - Placebo, respectively.

Alternatively, this could be represented as a matrix or a named list, if there is an opportunity for {rbmi} to make use of the names. For example, names could be propagated to outputs such as the table produced by pool().

I proposed the group_ prefix so that the argument is conceptually connected to the group argument of set_vars().

[tobiasmuetze]

Then, there must either be a clear definition of the relationship between ref, alt, et2, etc., and the levels of group, or the output table of pool() should follow the structure I proposed above (of course, with the parameter column retained for backward compatibility). Otherwise, group-level (label) identification may become error-sensitive.

Agree!

Constrast specification

At first sight, I am fine with the group_contrasts proposal, but need to think the implications fully through.

[wwojciech]

@tobiasmuetze - Can we please remove this from the assertion:

length(unique(data[[group]])) == length(levels(data[[group]]))

and then, of course, modify the code below that computes trt to:

trt <- lapply(
  levels(data2[["rbmiGroup"]])[-1],
  function(x) {
    term <- sprintf("rbmiGroup%s", x)
    if (term %in% names(coef(mod))) {
      list(
        est = coef(mod)[[term]],
        se = sqrt(vcov(mod)[term, term]),
        df = df.residual(mod)
      )
    } else {
      list(est = NA, se = NA, df = NA)
    }
  }
)

This is because there are some cases where we have no records for a certain level of data[[vars$group]], and this should not cause the function to stop.

[tobiasmuetze]

Is my understanding correct that if this assertion is removed, the group variable in the input data set could have a certain level but no subjects with this level? In which clinical trial scenario would this be relevant?

[wwojciech]

Is my understanding correct that if this assertion is removed, the group variable in the input data set could have a certain level but no subjects with this level? In which clinical trial scenario would this be relevant?

Yes, it is. There might be no data for a given level of data[[vars$group]].

This is relevant in scenarios where we fit an ANCOVA model using a subset of the full dataset defined by a subgroup variable (such as geographic region). In that case, for some levels of the subgroup variable, certain levels of data[[vars$group]] may be absent from the subset.

For example, in a given geographic region, we may not have data for all treatment groups. See the example table below for a dataset with three treatment groups (i.e., data[[vars$group]] has three levels), where the table presents results for only two treatment groups of interest.

[tobiasmuetze]

I’d suggest keeping the assertion, because this scenario seems to point to an inadequate model/data setup rather than something the function should accommodate.

For subgroup analyses based on multiply imputed data, the treatment-by-subgroup interaction should be included in the imputation model. If a treatment-by-subgroup combination is completely absent in the data, then fitting such an imputation model should already be problematic.

In that sense, the assertion length(unique(data[[group]])) == length(levels(data[[group]])) at least ensures that all declared treatment levels are actually represented in the dataset being analyzed. If a level is absent after subsetting, that should probably be treated as a signal that the analysis model is not adequate for the intended subgroup analysis.

The statistical concern is that if the imputation model does not include the treatment-by-subgroup interaction, subsequent subgroup-specific treatment effect estimates, or analyses of treatment-by-subgroup interaction, may be biased toward the null when interactions are present; see Schafer. More generally, Austin et al. state that “it is important to include in the imputation model all the variables that will be included in the analysis model.”

[wwojciech]

I’d suggest keeping the assertion, because this scenario seems to point to an inadequate model/data setup rather than something the function should accommodate.

For subgroup analyses based on multiply imputed data, the treatment-by-subgroup interaction should be included in the imputation model. If a treatment-by-subgroup combination is completely absent in the data, then fitting such an imputation model should already be problematic.

In that sense, the assertion length(unique(data[[group]])) == length(levels(data[[group]])) at least ensures that all declared treatment levels are actually represented in the dataset being analyzed. If a level is absent after subsetting, that should probably be treated as a signal that the analysis model is not adequate for the intended subgroup analysis.

The statistical concern is that if the imputation model does not include the treatment-by-subgroup interaction, subsequent subgroup-specific treatment effect estimates, or analyses of treatment-by-subgroup interaction, may be biased toward the null when interactions are present; see Schafer. More generally, Austin et al. state that “it is important to include in the imputation model all the variables that will be included in the analysis model.”

Thank you @tobiasmuetze

My understanding is as follows.

Assume that the imputation model was correctly specified (following the recommendations of Schafer and Austin), i.e., the treatment-by-subgroup interaction is included in the imputation model. Even in this case, it may still happen that, for certain levels of the subgroup variable, not all treatment groups are represented in the observed data. As a result, after imputation there will still be no observations for that particular treatment-by-subgroup combination, since multiple imputation imputes missing values for existing subjects rather than creating new observations.

However, I do not think this necessarily prevents a valid analysis within that subgroup. If the estimand of interest does not involve the missing treatment level, then contrasts among the remaining observed treatment groups remain estimable. In that situation, I would expect the ANCOVA model to fit successfully and to return valid estimates for the estimable contrasts, while returning NA for treatment effects involving the absent treatment level.

For this reason, I think the assertion may be too restrictive.

@luwidmer, @danielinteractive - what do you think?

[luwidmer]

What would you think of the following: relax the assertion, but not with a silent NA to address these two concerns:

• Imputation / analysis model alignment (@tobiasmuetze): the treatment-by-subgroup interaction must be in the imputation model, or subgroup effects bias toward the null (Schafer; Austin et al.). Agreed.
• Structural empty cells (@wwojciech): MI imputes outcomes for existing subjects, never creates new ones. So a congenial model can still leave an arm×subgroup cell empty by design (e.g. an arm that never opened in a region).

However, then length(unique(data[[group]])) == length(levels(data[[group]])) is the wrong check to do?
It is neither necessary (estimable contrasts don't need every level) nor sufficient (one subject per level does not guarantee the treatment-by-subgroup interaction being in the imputation model) for the above criteria I think?
If the estimand excludes the absent arm, the remaining contrasts could still be identified?

But I'd keep @tobiasmuetze's signal instead of suppressing it, e.g., replace stop() with a warning() naming the empty level, return estimable contrasts + NA for the rest. That distinguishes "empty by design" from a data-prep bug, which a silent NA can't?

[wwojciech]

Thank you @ludmier. I agree with your points. Having said that, I would be slightly cautious about issuing a warning here. A user could repeatedly get warnings for perfectly valid analyses, making the output noisy.

Is there any way to be more specific, i.e. identify a data-preparation issue such as imputation model misspecification and issue the warning only in such cases?

---

Apart from data-preparation issues, I would also distinguish between the following two cases:

1. An empty treatment level, but the requested contrasts are still estimable.
   -> I would not warn. Return the estimable contrasts and NA for the non-estimable ones.
2. The user explicitly requests a contrast involving the absent treatment.
   -> This is more actionable. A warning would make sense, since the requested estimand cannot be estimated.

That said, I am not trying to push strongly for one solution here. Whether you decide to keep the warning or not is ultimately fine with me, as long as we don't stop() :)

[wwojciech]

Also, @tobiasmuetze , please revisit the existing implementation of ancova_single_m_group(). The current code implicitly assumes that treatment (contr.treatment()) contrasts are used for the grouping factor. If a different contrast coding is active (e.g. contr.sum() or another custom contrast), the coefficient extraction no longer corresponds to treatment-vs-reference comparisons, which could lead to incorrect or misleading results without necessarily producing an error.

[tobiasmuetze]

Also, @tobiasmuetze , please revisit the existing implementation of ancova_single_m_group(). The current code implicitly assumes that treatment (contr.treatment()) contrasts are used for the grouping factor. If a different contrast coding is active (e.g. contr.sum() or another custom contrast), the coefficient extraction no longer corresponds to treatment-vs-reference comparisons, which could lead to incorrect or misleading results without necessarily producing an error.

I agree that the current coefficient extraction only supports the originally discussed setup, where each group level is compared against the first/reference level. However, based on our recent discussion, we now want to support user-specified contrasts, rather than assuming comparisons to the first level.

[luwidmer]

@tobiasmuetze & @wwojciech: any objections to merging into this branch the v2 version that implements user-specified contrasts (#590) for further discussion / iteration (if needed)?

[wwojciech]

@tobiasmuetze & @wwojciech: any objections to merging into this branch the v2 version that implements user-specified contrasts (#590) for further discussion / iteration (if needed)?

Hi @luwidmer - I have no objections. That said, I must admit I didn’t have a chance to check all the details and only had a very brief look at it.

[munoztd0]

@luwidmer and @tobiasmuetze sorry to bother but any plans to merge it into main now that the new API has been merged in ?

[danielinteractive]

I understand that we have an ooo situation so need to wait a couple of days - that said, if there is any support needed here happy to jump in any time from my side.

[danielinteractive]

Thanks again @luwidmer for the nice proposal, I just had a quick look through and it looks really nice. Also the testing seems to be very thorough which is great to see, especially that we get full backwards compatibility. Minor things below.

[luwidmer]

@gowerc could you take a look as to whether you're fine with what I did to your branch 😅 ?

[tobiasmuetze]

I reviewed the functionality from a user perspective. I have no suggestions for changes. Great work and very useful.