Treatment-group-specific variances

[wlandau]

The interface could be brm_formula(..., variance = "heterogeneous_group"). (Homogeneous group-specific variances would be an odd choice, and I see no reason for it.) Then in brm_marginal_draws(), the calculation of effect size would need to use group-specific variances.

[wlandau]

Glancing at https://opensource.nibr.com/bamdd/src/02h_mmrm.html, I see sigma ~ 1 + AVISIT + TRT01P + AVISIT:TRT01P. This makes me think users are going to want the full range of options for sigma, including maybe subgroup-specific variances. So this is more complicated than I originally thought.

Here is what I propose:

1. Create a new function brm_formula_sigma() to define a sigma formula with the same kind of guardrails and formal arguments as brm_formula().
2. Add a new sigma argument to brm_formula() which accepts a formula from brm_formula_sigma(). Use this instead of the current variance argument of brm_formula().
3. In get_draws_sigma() (part of brm_draws_marginal()), use brms::make_standata(formula, data)$X_sigma derive time-specific variances from whatever distributional parameterization for sigma that the user throws at us. This is similar to how the custom marginal mean estimation works (through brm_transform_marginal()).

[andrew-bean]

So this is more complicated than I originally thought.

One interesting wrinkle from the vignette is: the contrasts(AVISIT) <- MASS::contr.sdif(...) setting affects the parametrization of the (log linear) model for sigma as well as the model for the mean.

In our package ... should the brms.mmrm archetype chosen by the user affect both the mean parametrization and the variance parametrization? I guess yes?

I looks like sigma ~ 1 + AVISIT + TRT01P + AVISIT:TRT01P would explicitly model subject-specific effects.

I don't think there is any subject-specific effect here. Suppose there are 4 treatment groups and 4 visits like in my colleagues' vignette. The formula above allows 16 different variances in these group-by-visit cells. This is more flexible than what is assumed in the current brms.mmrm model, which is 4 variances (one per visit) assumed to be shared across treatment group within visits. But it is not so flexible as to allow different residual variances at the patient level.

[chstock]

I don't think there is any subject-specific effect here.

I agree. I got this wrong at first (and deleted my comment). It's good to have the possibility to model group specific covariance structures.

[wlandau]

One interesting wrinkle from the vignette is: the contrasts(AVISIT) <- MASS::contr.sdif(...) setting affects the parametrization of the (log linear) model for sigma as well as the model for the mean.

Yikes, I'm glad you pointed that out. This makes me think we should prohibit contrasts in all the important factors we are modeling.

In our package ... should the brms.mmrm archetype chosen by the user affect both the mean parametrization and the variance parametrization? I guess yes?

I think it depends on what we are trying to do, and on the scope of the package. For archetypes, we micromanage fixed effects to make it easy and safe to set informative priors. Would we want to do the same for variances? Would we expect users to set carefully-constructed priors on variance components using historical data / expert elicitation? I think this is relatively uncommon, and it is free from the hidden pitfalls of covariate adjustment, so I think it would be feasible for advanced users to directly call brms::set_prior() for this.

I initially assumed that the variance parameterization would be independent of the fixed effect parameterization and handled in a consistent way across all archetype and non-archetype analyses. I think this is feasible because the archetype still retains all the original columns from the data, so you can still write a formula like sigma ~ 1 + AVISIT + TRT01P + AVISIT:TRT01P

[wlandau]

3. In get_draws_sigma() (part of brm_draws_marginal()), use brms::make_standata(formula, data)$X_sigma derive time-specific variances from whatever distributional parameterization for sigma that the user throws at us. This is similar to how the custom marginal mean estimation works (through brm_transform_marginal()).

This is needed because of the "effect size" calculation in brm_marginal_draws(). "Effect size" in this case is the treatment effect divided by sigma. Depending on how sigma is parameterized, this could turn out to be an extremely deep rabbit hole. For example, if we allow sigma ~ AVISIT*TRT01P + AVISIT*BASE + SEX + RACE, then we would need to replicate or extend all the existing marginal mean estimation machinery and make it work for subgroups and archetypes. This does not seem feasible from a package maintenance perspective, especially given that effect size is relatively unimportant for inference. I would prefer to restrict support for effect size to cases without baseline or nuisance variables.