Should fixed effect priors be correlated in longitudinal models?

[wlandau]

@bjoernholzhauer, @chstock, @yonicd, and I have been discussing if priors on fixed effects in longitudinal models should account for correlations. This is highly relevant to how we think about MAP priors for MMRMs.

I ran a little experiment which shows that in at least one situation, correlations in the prior do not contribute much if the likelihood already models them. I do not know how this generalizes, however. It would be interesting to compare MMRMs with vs without correlations in the priors on the beta parameters.

A quick preliminary frequentist look

Let’s simulate simple longitudinal data with 2 correlated time points.

library(tidyverse)
library(mvtnorm)
set.seed(1234)
example_wide <- rmvnorm(
  n = 1e4L,
  mean = c(0, 0),
  sigma = matrix(c(1, 0.5, 0.5, 1), nrow=2, ncol=2)
) %>%
  as.data.frame() %>%
  as_tibble() %>%
  mutate(
    TRT01P = 1L * (row_number() > n() / 2),
    USUBJID = as.factor(paste0("PATIENT_", row_number()))
  )
example_wide
#> # A tibble: 10,000 × 4
#>        V1      V2 TRT01P USUBJID   
#>     <dbl>   <dbl>  <int> <fct>     
#>  1 -1.09  -0.0444      0 PATIENT_1 
#>  2  0.440 -1.99        0 PATIENT_2 
#>  3  0.545  0.600       0 PATIENT_3 
#>  4 -0.697 -0.677       0 PATIENT_4 
#>  5 -0.776 -1.01        0 PATIENT_5 
#>  6 -0.719 -1.09        0 PATIENT_6 
#>  7 -0.733 -0.139       0 PATIENT_7 
#>  8  0.898  0.142       0 PATIENT_8 
#>  9 -0.729 -1.01        0 PATIENT_9 
#> 10 -0.183  2.12        0 PATIENT_10
#> # ℹ 9,990 more rows

example_long <- example_wide %>%
  pivot_longer(
    cols = c("V1", "V2"),
    names_to = "AVISIT",
    values_to = "AVAL"
  ) %>%
  mutate(AVISIT = as.factor(AVISIT))
example_long
#> # A tibble: 20,000 × 4
#>    TRT01P USUBJID   AVISIT    AVAL
#>     <int> <fct>     <fct>    <dbl>
#>  1      0 PATIENT_1 V1     -1.09  
#>  2      0 PATIENT_1 V2     -0.0444
#>  3      0 PATIENT_2 V1      0.440 
#>  4      0 PATIENT_2 V2     -1.99  
#>  5      0 PATIENT_3 V1      0.545 
#>  6      0 PATIENT_3 V2      0.600 
#>  7      0 PATIENT_4 V1     -0.697 
#>  8      0 PATIENT_4 V2     -0.677 
#>  9      0 PATIENT_5 V1     -0.776 
#> 10      0 PATIENT_5 V2     -1.01  
#> # ℹ 19,990 more rows

Let’s fit an MMRM with unstructured covariance.

library(mmrm)
fit <- mmrm(
  formula = AVAL ~ 0 + AVISIT + us(AVISIT | USUBJID),
  data = example_long
)
attr(fit, "converged")
#> [1] TRUE

We see the 50% correlation between the fixed effects.

vcov(fit) %>%
  cov2cor()
#>           AVISITV1  AVISITV2
#> AVISITV1 1.0000000 0.4943665
#> AVISITV2 0.4943665 1.0000000

And the same correlation between residuals within patients.

fit$cov
#>           V1        V2
#> V1 0.9896992 0.4893578
#> V2 0.4893578 0.9900374

Bayesian models

Next, we compare several Bayesian models which account for different
forms of correlation.

library(cmdstanr)
library(posterior)
run_model <- function(file) {
  y <- as.matrix(example_wide[, c("V1", "V2")])
  n <- nrow(y)
  s <- matrix(c(1, 0.5, 0.5, 1), nrow = 2)
  data <- list(n = n, y = y, s = s)
  model <- cmdstan_model(file)
  model$sample(
    data = data,
    chains = 4,
    iter_warmup = 2e3,
    iter_sampling = 1e4,
    cores = 4,
    refresh = 0
  )
}
print_results <- function(fit, variables = c("beta", "R")) {
  print(fit$summary(variables = variables))
  print("")
  print(cor(as_draws_df(fit)[, c("beta[1]", "beta[2]")]))
}

“Uncorrelated” model

Consider a Bayesian model which does not account for correlations either
at the level of the data or the model.

writeLines(readLines("uncorrelated.stan"))
#> data {
#>   int<lower=0> n;
#>   array[n] vector[2] y;
#> }
#> parameters {
#>   vector[2] beta;
#> }
#> model {
#>   y ~ multi_normal(beta, identity_matrix(2));
#>   beta ~ multi_normal(rep_vector(0, 2), identity_matrix(2));
#> }

We fit the model and summarize posterior samples.

fit_uncorrelated <- run_model("uncorrelated.stan")

print_results(fit_uncorrelated, "beta")
#> # A tibble: 2 × 10
#>   variable     mean   median     sd     mad      q5     q95  rhat ess_bulk
#>   <chr>       <dbl>    <dbl>  <dbl>   <dbl>   <dbl>   <dbl> <dbl>    <dbl>
#> 1 beta[1]   0.00510  0.00512 0.0100 0.00999 -0.0115 0.0217   1.00   37891.
#> 2 beta[2]  -0.00712 -0.00703 0.0100 0.0101  -0.0236 0.00922  1.00   35689.
#> # ℹ 1 more variable: ess_tail <dbl>
#> [1] ""
#> Warning: Dropping 'draws_df' class as required metadata was removed.
#>              beta[1]      beta[2]
#> beta[1]  1.000000000 -0.002199508
#> beta[2] -0.002199508  1.000000000

Correlations in the likelihood

Next, consider correlations only in the likelihood.

writeLines(readLines("likelihood.stan"))
#> data {
#>   int<lower=0> n;
#>   array[n] vector[2] y;
#> }
#> parameters {
#>   vector[2] beta;
#>   corr_matrix[2] R;
#> }
#> model {
#>   y ~ multi_normal(beta, R);
#>   beta ~ normal(0, 1);
#>   R ~ lkj_corr(1);
#> }

We fit the model and summarize posterior samples.

fit_likelihood <- run_model("likelihood.stan")

print_results(fit_likelihood)
#> # A tibble: 6 × 10
#>   variable     mean   median      sd     mad      q5     q95  rhat ess_bulk
#>   <chr>       <dbl>    <dbl>   <dbl>   <dbl>   <dbl>   <dbl> <dbl>    <dbl>
#> 1 beta[1]   0.00502  0.00494 0.00994 0.0100  -0.0112 0.0214   1.00   25391.
#> 2 beta[2]  -0.00706 -0.00704 0.00998 0.00997 -0.0236 0.00931  1.00   25123.
#> 3 R[1,1]    1        1       0       0        1      1       NA         NA 
#> 4 R[2,1]    0.497    0.497   0.00676 0.00673  0.486  0.508    1.00   28561.
#> 5 R[1,2]    0.497    0.497   0.00676 0.00673  0.486  0.508    1.00   28561.
#> 6 R[2,2]    1        1       0       0        1      1       NA         NA 
#> # ℹ 1 more variable: ess_tail <dbl>
#> [1] ""
#> Warning: Dropping 'draws_df' class as required metadata was removed.
#>           beta[1]   beta[2]
#> beta[1] 1.0000000 0.4963262
#> beta[2] 0.4963262 1.0000000

Correlations in the prior

Next, consider correlations only in the prior.

writeLines(readLines("prior.stan"))
#> data {
#>   int<lower=0> n;
#>   array[n] vector[2] y;
#>   matrix[2, 2] s;
#> }
#> parameters {
#>   vector[2] beta;
#>   corr_matrix[2] R;
#> }
#> model {
#>   y ~ multi_normal(beta, R);
#>   beta ~ multi_normal(rep_vector(0, 2), s);
#>   R ~ lkj_corr(1);
#> }

We fit the model and summarize posterior samples.

fit_prior <- run_model("prior.stan")

print_results(fit_prior)
#> # A tibble: 6 × 10
#>   variable     mean   median      sd     mad      q5     q95  rhat ess_bulk
#>   <chr>       <dbl>    <dbl>   <dbl>   <dbl>   <dbl>   <dbl> <dbl>    <dbl>
#> 1 beta[1]   0.00504  0.00501 0.00993 0.00994 -0.0112 0.0214   1.00   27589.
#> 2 beta[2]  -0.00712 -0.00712 0.00997 0.00993 -0.0236 0.00928  1.00   27296.
#> 3 R[1,1]    1        1       0       0        1      1       NA         NA 
#> 4 R[2,1]    0.497    0.497   0.00679 0.00674  0.486  0.508    1.00   30903.
#> 5 R[1,2]    0.497    0.497   0.00679 0.00674  0.486  0.508    1.00   30903.
#> 6 R[2,2]    1        1       0       0        1      1       NA         NA 
#> # ℹ 1 more variable: ess_tail <dbl>
#> [1] ""
#> Warning: Dropping 'draws_df' class as required metadata was removed.
#>           beta[1]   beta[2]
#> beta[1] 1.0000000 0.4939576
#> beta[2] 0.4939576 1.0000000

Correlations in both

Next, consider correlations in both the likelihood and the prior.

writeLines(readLines("both.stan"))
#> data {
#>   int<lower=0> n;
#>   array[n] vector[2] y;
#>   matrix[2, 2] s;
#> }
#> parameters {
#>   vector[2] beta;
#>   corr_matrix[2] R;
#> }
#> model {
#>   y ~ multi_normal(beta, R);
#>   beta ~ multi_normal(rep_vector(0, 2), s);
#> }

We fit the model and summarize posterior samples.

fit_both <- run_model("both.stan")

print_results(fit_both)
#> # A tibble: 6 × 10
#>   variable     mean   median      sd     mad      q5     q95  rhat ess_bulk
#>   <chr>       <dbl>    <dbl>   <dbl>   <dbl>   <dbl>   <dbl> <dbl>    <dbl>
#> 1 beta[1]   0.00511  0.00509 0.00998 0.00997 -0.0113 0.0216   1.00   26176.
#> 2 beta[2]  -0.00700 -0.00700 0.00998 0.0100  -0.0234 0.00953  1.00   25484.
#> 3 R[1,1]    1        1       0       0        1      1       NA         NA 
#> 4 R[2,1]    0.497    0.497   0.00672 0.00676  0.486  0.508    1.00   30453.
#> 5 R[1,2]    0.497    0.497   0.00672 0.00676  0.486  0.508    1.00   30453.
#> 6 R[2,2]    1        1       0       0        1      1       NA         NA 
#> # ℹ 1 more variable: ess_tail <dbl>
#> [1] ""
#> Warning: Dropping 'draws_df' class as required metadata was removed.
#>          beta[1]  beta[2]
#> beta[1] 1.000000 0.490407
#> beta[2] 0.490407 1.000000

[chstock]

Thank you, @wlandau, for this case example. Interesting observation ...and food for thought and further discussions. It seems necessary to get our head around how the correlations would contribute, if at all. I don't have an immediate explanation.

[bjoernholzhauer]

To repost what I had sent by email for public discussion. To add to this, I suspect this matters the most when there is a lot of consistent prior information and very little new data. In that case, ignoring prior correlations between parameters is unlikely to be fixed by having a lot of data and then things could go wrong (potentially badly). And I suspect this is exactly the kind of scenario where we're the most likely to do something Bayesian (e.g. early stage proof-of-concept type of study in a well-studied indication).

Original thoughts on this:
I’m not 100% sure, but I’d be worried that doing independent marginal priors on each coefficient would miss something important about the posterior from the previous analysis.

The reason I think so is that we know that a-posteriori (or simply in the variance-covariance matrix of maximum likelihood estimates) there’s kind of mathematically guaranteed to be a correlation between samples for different model terms. E.g. if you have an intercept and a slope, you can always move the intercept higher a bit and then the slope needs to get steeper, or you move it down a little and the slope gets shallower. That’s kind of independent on whether one works in a frequentist or Bayesian setting (I’ll illustrate using frequentist approaches). It’s easy to illustrate this effect with a simple linear model. Here, I simulated something super-basic (just one post-baseline visit + a baseline) and got following correlation of coefficient from a frequentist linear regression:

            (Intercept) TRT01P  BASE
(Intercept)        1.00  -0.58 -0.58
TRT01P            -0.58   1.00  0.00
BASE              -0.58   0.00  1.00

In a MMRM fit on fev_data from the mmrm package (model: FEV1 ~ AVISIT + ARMCD + ARMCD:AVISIT + RACE + SEX + us(AVISIT | USUBJID)), there are more coefficients that are strongly correlated:

1.00  -0.59  -0.72 -0.38 -0.60  -0.41 -0.30 -0.36  0.42    0.50   0.27
-0.59   1.00   0.61  0.35  0.49   0.01 -0.01  0.02 -0.71   -0.42  -0.25
-0.72   0.61   1.00  0.37  0.59   0.02 -0.01  0.02 -0.43   -0.69  -0.26
-0.38   0.35   0.37  1.00  0.31   0.00 -0.01  0.01 -0.25   -0.26  -0.71
-0.60   0.49   0.59  0.31  1.00   0.04 -0.06  0.00 -0.70   -0.82  -0.44
-0.41   0.01   0.02  0.00  0.04   1.00  0.44  0.09  0.01    0.02  -0.01
-0.30  -0.01  -0.01 -0.01 -0.06   0.44  1.00  0.07  0.02    0.01   0.00
-0.36   0.02   0.02  0.01  0.00   0.09  0.07  1.00 -0.03   -0.04  -0.01
0.42  -0.71  -0.43 -0.25 -0.70   0.01  0.02 -0.03  1.00    0.60   0.35
0.50  -0.42  -0.69 -0.26 -0.82   0.02  0.01 -0.04  0.60    1.00   0.36
0.27  -0.25  -0.26 -0.71 -0.44  -0.01  0.00 -0.01  0.35    0.36   1.00

(column/row order: (Intercept) AVISITVIS2 AVISITVIS3 AVISITVIS4 ARMCDTRT RACEBlack or African American RACEWhite SEXFemale AVISITVIS2:ARMCDTRT AVISITVIS3:ARMCDTRT AVISITVIS4:ARMCDTRT).

My suspicion is that these are such strong correlations that they ideally shouldn’t be ignored. I admit though that I’m not 100% sure how much things go wrong if one ignores these correlations.

You are right though, this is a dangerous territory, where getting wrong which indices to intervene on (or messing up the ordering of indices etc.) could result in setting priors on different coefficients than intended. If one finds a reliable way of making that work automatically, then that’s a benefit of having that inside a package that avoids error-prone manual specification.

I’m not entirely sure whether many univariate mixture approximations cause more problems than one multivariate one. The time to do it and the time-cost during MCMC fitting is probably not that big a deal. The univariate ones are easier to check visually (and otherwise), while for the multivariate one that’s harder (although marginal distributions would still be easy to compare).

R code used:

library(tidyverse)
library(mvtnorm)

set.seed(1234)
example <- rmvnorm(n=100000, mean=c(0, 0), sigma=matrix(c(1, 0.5, 0.5, 1), nrow=2, ncol=2)) %>%
  as_tibble() %>%
  rename(BASE=V1, AVAL=V2) %>%
  filter(BASE > -0.25) %>%
  mutate(TRT01P = 1L*(row_number()>n()/2),
         AVAL = AVAL - 0.5*TRT01P,
         CHG = AVAL - BASE)

lmfit1 <- lm(CHG ~ 1 + TRT01P + BASE,
             data=example)

vcov(lmfit1) %>% 
  cov2cor() %>%
  round(2)

library(mmrm)
fit <- mmrm(
  formula = FEV1 ~ AVISIT + ARMCD + ARMCD:AVISIT + RACE + SEX  + us(AVISIT | USUBJID),
  data = fev_data
)

vcov(fit) %>%
  cov2cor() %>%
  round(2)

[wlandau]

@bjoernholzhauer, this makes sense. I added a section on multivariate mixture priors in a new MAP prior vignette: #138. If you're curious and happen to spot any errors, or if you see an easier way to do this, please let me know.