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]

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.