brm_formula_sigma() default arguments

[wlandau]

An argument could be made that brm_formula_sigma() should have intercept = TRUE, group = TRUE, and group_time = TRUE by default, e.g. to account for heteroskedasticity and different variances for different study arms. The current defaults are all FALSE in order to align with conventional/frequentist MMRMs that do not support distributional regression, but maybe this is not a good reason. What are your thoughts?

[bjoernholzhauer]

Does a frequentist MMRM with the commonly used choice of separate unstructured covariance matrices for each treatment group not result in effectively distributional regression with terms for group, time and group:time? The current default of time = TRUE (i.e. distributional regression with the only predictor being the time factor variable, so distributional regression is already being used) makes sense to match a single unstructured covariance matrix that is assumed to be identical for all groups but is allowed to have different variances at each visit. Whether to include an intercept is much less clear and seems to me to be more of a choice about how you feel comfortable specifying your prior knowledge.

[wlandau]

Before I began using brms, I had mostly seen timepoint-specific variances in multivariate models. I could be wrong, but distributional regression in its full generality has always seemed to me like something unique to brms.

I originally learned about the "separation strategy"/LKJ from this paper, which defines timepoint-specific sigma parameters.

Another example is the mmrm package, which AFAIK assumes there is one single T x T correlation matrix, where T is the number of time points.

library(mmrm)
fit <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT | USUBJID),
  data = fev_data
)
fit$cov
#>           VIS1      VIS2       VIS3       VIS4
#> VIS1 40.553664 14.396045  4.9747288 13.3866534
#> VIS2 14.396045 26.571483  2.7854661  7.4744790
#> VIS3  4.974729  2.785466 14.8978517  0.9082111
#> VIS4 13.386653  7.474479  0.9082111 95.5568420

fit <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(ARMCD | USUBJID),
  data = fev_data
)
#> Error in h_mmrm_tmb_data(formula_parts, data, weights, reml, singular = if (control$accept_singular) "drop" else "error",  : 
#>   time points have to be unique for each subject, detected following duplicates in data:
#>     USUBJID ARMCD
#> 2       PT1   TRT
#> 3       PT1   TRT
#> ...

fit <- mmrm(
  formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT * ARMCD | USUBJID),
  data = fev_data
)
#> Error in mmrm(formula = FEV1 ~ RACE + SEX + ARMCD * AVISIT + us(AVISIT *  : 
#>   1 assertions failed:
#>  * Non-spatial covariance structures must have a single longitudinal variable

But again, even if most people have done it this way in the past, maybe brms.mmrm could be more flexible.

[bjoernholzhauer]

When I started in industry, I was taught (i.e. may just be the local statistical folklore in the group where I started in industry) that this is "the" default MMRM:

PROC MIXED;
  CLASS TRT01P AVISIT USUBJID;
  MODEL CHG ~ TRT01P AVISIT BASE TRT01P*AVISIT AVISIT*BASE 
    / SOLUTION DDFM=KR ALPHA = 0.05;
  REPEATED AVISIT / TYPE=UN SUBJECT = USUBJID R Rcorr GROUP=TRT01P;
  LSMEANS TRT01P*AVISIT / DIFFS PDIFF CL OM E;
RUN;

Note that GROUP=TRT01P means we get separate covariance matrices per group. While brms does not currently let you specify two separate correlation matrices, you can at least allow different variances at each visit for each treatment group by doing distributional regression for the variance with terms for visit, treatment and visit:treatment. The matching code for the frequentist mmrm package would be to use ... + us(AVISIT | TRT01P / USUBJID) in the model formula.

Arguably, this is more robust vs. different variances by group, but people frequently fit linear regression for a single visit assuming equal variance for all treatment groups. The MMRM above ships avoids that assumption (it's of course hard to say in general how much one should be worried about potentially only mild deviations from that assumption).

[wlandau]

You make compelling points, @bjoernholzhauer. @chstock, @yonicd, @andrew-bean, @kkmann, what do you think?