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?

[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.

[gowerc]

AI review flagged 2 issues for me:

• (major) the LSMs were not invariant to contrasts as it looks like the code was just extracting the beta coeficient which only works for corner-point constraint. Added a bunch more processing code to handle the different contrasts by constructing a small reference matrix and then using that to pick out / sum the relevant beta coeficients.

• (minor) there was an issue where you could get a parameter name clash if you had overlapping visit / treatment labels. I didn't bother fixing this as you have to be going out of your way to cause it so just put an assertion that looks for duplicate labels and then throws an error.

9f6c0b8 @danielinteractive / @luwidmer if you have a chance to see what you think.

Other wise I think the code looks good to me.

[danielinteractive]

Super, thanks @gowerc - I agree that the now general encoding support for the factors is a good idea.

[gowerc]

On reflection think I've confused myself as to what exactly we want or expect the contrast to represent. Apologies if this was in the above conversation, I tried reading but I was struggling to follow :(

So in the current implementation on main we naively force the use corner-point-constraint (contr.treatment) and then just report the treatment contrast effect as the corresponding regression Beta coeficient. In hindsight I'm not even sure if this is correct as this is implicitly equivalent to setting all covariate values + interactions to 0 e.g. ignoring interactions and taking the reference level for categorical covariates. For example is it meaningful to have an age of 0 or assume everyone is Male ? (Edit - I appreciate if theres no treatment interaction terms then the implementation here is fine, but I don't think its right if there is a treatment interaction term).

The current implementation in this PR (as of my most recent commit) slightly addresses this by making the choice of covariates explicit (internally at least) and allowing for non-corner-point-constraint codings. It does this by building a reference dataset for each treatment level where continuous covariates are set to 0 and categorical covariates are set to the reference level using the contrast coding scheme e.g.

corner-point

> mod_treatment <- lm(out ~ rbmiGroup * sex + I(age^2), data = dat)
> design_treatment <- ancova_contrast_design(mod_treatment, dat, "out")
> design_treatment
  (Intercept) rbmiGroupL2 sexMale I(age^2) rbmiGroupL2:sexMale
1           1           0       0        0                   0
2           1           1       0        0                   0

sum-to-zero

> mod_sum <- lm(out ~ rbmiGroup * sex + I(age^2), data = dat)
> design_sum <- ancova_contrast_design(mod_sum, dat, "out")
> design_sum
  (Intercept) rbmiGroup1 sex1 I(age^2) rbmiGroup1:sex1
1           1          1    1        0               1
2           1         -1    1        0              -1

This is then used by the user defined contrasts list to decide which rows / fitted values be calculated and then taking the difference between them.

I'm just a bit lost now as to wether it still makes sense to be setting covariates to 0 / the reference level for categorical values? Would it not be more reasonable to just do the difference of LSMs which then has weighting schemes to determine what values to set the covariates at which is representative ?

@tobiasmuetze / @danielinteractive / @luwidmer / @wwojciech

[wwojciech]

Hi @gowerc

I think these are very valid considerations. A bit aside from the technical discussion, I was wondering whether the current state of the PR would allow us to split this into a separate issue, so that the custom contrast could be considered/implemented in a dedicated PR, while this PR could be completed with the multi-arm support?

If this would not introduce too much mess or require significant additional work, I would be very much in favour of this option. The package I support with @danielinteractive currently requires the multi-arm support first, while, as far as I am aware, the custom contrast is less critical for us (if at all).

@luwidmer , @tobiasmuetze - for your consideration.

[gowerc]

hmm, I mean the concern I raised is tangental to the custom contrast, I think the custom contrast code implemented here by @luwidmer is working as intended. My issue is more around wether it is creating a meaningful metric. Though I think this is an issue with the main version of rbmi as well. Given the code is correct as per its intended purpose I am happy to merge this into main and then create a separate issue for the interactions conversation but I do think the interactions bit needs resolution before I would feel happy pushing to CRAN.

[danielinteractive]

@gowerc yeah makes sense. Maybe @wwojciech you can investigate this concern/question with examples to check if it is a problem?

[gowerc]

I have moved the interaction conversation then to #595 , I am happy to merge but there appears to be a lingering unresolved conversation thread which is blocking the merge (I'd rather not admin force merge unless I really have to)

[danielinteractive]

@gowerc I think it has been resolved now