Replies: 3 comments 6 replies
-
|
This is great! Looking forward to discussing the next steps! |
Beta Was this translation helpful? Give feedback.
-
|
This looks great! in the past i have had to construct mean over time points as the endpoint, eg if t timepoints then each gets 1/t weight instead of 1 (or even a nonlin functional if the clinical team is feeling like they want to test a hypothesis) would this be possible in this construction? |
Beta Was this translation helpful? Give feedback.
-
|
Conditional means priors work best when the transformation matrix is full rank, and averaging across time points might be tricky. I think we would have to set one parameter to be the average and the others to be time-specific contrasts on the average. Kind of like "sum contrasts" from |
Beta Was this translation helpful? Give feedback.
-
|
Nice work as usual @wlandau ! Still thinking about @weberse2 's suggestion in #40. IMO setting the contrast attribute of the AVISIT factor is a simpler way to reparametrize the model and specify the prior on some alternative linear function of the model coefficients. In the code above we need to specify a new formula which does not resemble the usual MMRM By contrast (sorry), if you simply set the contrast attribute for AVISIT, you avoid both these inconveniences. With the alternative contrast attribute on AVISIT, brms will treat this as a reparametrized model, even though the formula is unchanged from the typical MMRM specification. Basic R functions like lm behave in a similar way. In brms, the consequence is you scan specify your prior on the forward difference contrasts rather than the cell means, without writing a completely new model. The syntax would be something like this: suppressPackageStartupMessages({
library(brms)
library(coda)
library(emmeans)
library(mmrm)
library(posterior)
library(tidyverse)
library(zoo)
})
data(fev_data, package = "mmrm")
data <- fev_data %>%
mutate(FEV1_CHG = FEV1 - FEV1_BL) %>%
select(-FEV1) %>%
group_by(USUBJID) %>%
complete(
AVISIT,
fill = as.list(.[1L, c("ARMCD", "FEV1_BL", "RACE", "SEX", "WEIGHT")])
) %>%
ungroup() %>%
arrange(USUBJID, AVISIT) %>%
mutate(USUBJID = as.character(USUBJID)) %>%
group_by(USUBJID) %>%
mutate(FEV1_CHG = na.locf(FEV1_CHG, na.rm = FALSE)) %>%
mutate(FEV1_CHG = na.locf(FEV1_CHG, na.rm = FALSE, fromLast = TRUE)) %>%
ungroup() %>%
filter(!is.na(FEV1_CHG))
formula_cell_means <- brmsformula(FEV1_CHG ~ 0 + AVISIT : ARMCD)
# we can just use a typical MMRM formula here
formula_sdif <- brmsformula(FEV1_CHG ~ AVISIT * ARMCD)
# form forward-difference contrasts
cc <- MASS::contr.sdif(nlevels(data$AVISIT))
cc
# must say I have not wrapped my head around the output from this
# MASS function. Its Moore-Penrose G-inverse makes some sense, in that it
# maps visits to visit differences.
zapsmall(MASS::ginv(cc))
# but this seems to work:
cc_alt <- diag(1, nlevels(data$AVISIT))
colnames(cc_alt) <- c("1", colnames(cc))
rownames(cc_alt) <- levels(data$AVISIT)
cc_alt[lower.tri(cc_alt)] <- 1
cc_alt
solve(cc_alt)
# new data frame in which nothing has changed except the contrasts attribute
# of the AVISIT factor
data_sdif <- data
contrasts(data_sdif$AVISIT) <- cc_alt[,-1] # must drop the intercept from the contrast definition
model_matrix_cell_means <- make_standata(
formula = formula_cell_means,
data = data
)$X
model_matrix_sdif <- make_standata(
formula = formula_sdif,
data = data_sdif
)$X
head(model_matrix_cell_means, n = 8)
head(model_matrix_sdif, n = 8)
get_prior(formula_cell_means, data = data)
# prior now specified for the reparametrized coefficients, even though the
# formula is unchanged
get_prior(formula_sdif, data = data_sdif)
# # cell means model
# fit_cell_means <- brm(formula_cell_means, data, sample_prior = "only",
# prior = c(prior(normal(0, 1), class = "b")))
# reparametrized model in terms of successive differences in visits
fit_sdif <- brm(formula_sdif, data_sdif, sample_prior = "only",
prior = c(prior(normal(0, 1), class = "b")))
# prior draws for the coefficients under each prior specification
# beta_cell_means <- as.matrix(fit_cell_means)
beta_sdif <- as.matrix(fit_sdif)
# beta_cell_means <- beta_cell_means[, grep("^b_*", colnames(beta_cell_means))]
beta_sdif <- beta_sdif[, grep("^b_*", colnames(beta_sdif))]
# to do inference on the cell means, you could just do something like this
# (no need for a transformation)
# define the cells
cells <- data_sdif %>%
select(ARMCD, AVISIT) %>%
arrange(ARMCD, AVISIT) %>%
distinct() %>%
mutate(FEV1_CHG = 0,
cell_label = interaction(ARMCD, AVISIT))
# posterior draws for the linear predictor at these levels
draws <- posterior_linpred(fit_sdif, newdata = cells)
|
Beta Was this translation helpful? Give feedback.
-
|
Thanks for posting, @andrew-bean. As we discussed today, a contrast-based implementation looks far cleaner and more achievable than manually constructing the model matrix from scratch. I think we may still need to use model matrices (from |
Beta Was this translation helpful? Give feedback.
-
|
Good news: the G-inverse should let us easily transform posterior samples from sdif to cell means: transform_sdif_to_cells <- MASS::ginv(model_matrix_sdif) %*%
model_matrix_cell_means
transformed <- model_matrix_sdif %*% transform_sdif_to_cells
max(abs(model_matrix_cell_means - transformed))
#> [1] 1.239092e-14A stumbling block: like you mentioned in our meeting yesterday @andrew-bean, the sdif parameterization is treatment-effect-like and not cell-means-like. In other words, for the treatment group, each successive difference is defined relative to the placebo group. As you said, this may make it difficult to set priors in an intuitive way. > head(model_matrix_sdif, 8)
(Intercept) AVISIT1 AVISIT2 AVISIT3 ARMCDTRT AVISIT1:ARMCDTRT AVISIT2:ARMCDTRT AVISIT3:ARMCDTRT
1 1 0 0 0 1 0 0 0
2 1 1 0 0 1 1 0 0
3 1 1 1 0 1 1 1 0
4 1 1 1 1 1 1 1 1
5 1 0 0 0 0 0 0 0
6 1 1 0 0 0 0 0 0
7 1 1 1 0 0 0 0 0
8 1 1 1 1 0 0 0 0I have experimented with many different contrast matrices for (Intercept) AVISIT1 AVISIT2 AVISIT3 ARMCDTRT AVISIT1:ARMCDTRT AVISIT2:ARMCDTRT AVISIT3:ARMCDTRT
1 0 0 0 0 1 0 0 0
2 0 0 0 0 1 1 0 0
3 0 0 0 0 1 1 1 0
4 0 0 0 0 1 1 1 1
5 1 0 0 0 0 0 0 0
6 1 1 0 0 0 0 0 0
7 1 1 1 0 0 0 0 0
8 1 1 1 1 0 0 0 0Unfortunately, unless we can find a contrast matrix that produces this model matrix, the behind-the-scenes mechanism might have to manually construct said model matrix and define a special formula to go with it. |
Beta Was this translation helpful? Give feedback.
-
As it turns out, the Kronecker product makes this easy: matrix_group <- model.matrix(~ 0 + ARMCD, filter(data, AVISIT == AVISIT[1]))
matrix_time <- diag(4) + lower.tri(diag(4))
matrix_sdif <- kronecker(matrix_treatment, matrix_time)
colnames(matrix_sdif) <- c(paste0("PBO|", seq_len(4)), paste0("TRT|", seq_len(4)))
head(matrix_sdif, n = 8)
#> PBO|1 PBO|2 PBO|3 PBO|4 TRT|1 TRT|2 TRT|3 TRT|4
#> [1,] 0 0 0 0 1 0 0 0
#> [2,] 0 0 0 0 1 1 0 0
#> [3,] 0 0 0 0 1 1 1 0
#> [4,] 0 0 0 0 1 1 1 1
#> [5,] 1 0 0 0 0 0 0 0
#> [6,] 1 1 0 0 0 0 0 0
#> [7,] 1 1 1 0 0 0 0 0
#> [8,] 1 1 1 1 0 0 0 0 |
Beta Was this translation helpful? Give feedback.
-
We might actually want treatment-effect-like sdif in certain situations. It depends on whether we want to model the treatment groups as independent in the prior, or if we want them correlated. This could be weighed against intuition on a case-by-case basis. |
Beta Was this translation helpful? Give feedback.
-
Implemented in #100 via On the original thread, I think we landed well in #96 and #100 with informative prior archetypes. |
Beta Was this translation helpful? Give feedback.
-
This post is a prototype of #40. Borrowing from the methodology of conditional means priors with an identity link function, we put a prior on time differences and get samples of cell means.
Goals
Data
suppressPackageStartupMessages({ library(brms) library(coda) library(emmeans) library(mmrm) library(posterior) library(tidyverse) library(zoo) })We begin with the
FEVdataset from themmrmpackage and informallyimpute missing values (LOCF then NOCB).
Throughout, we highlight patients 1 (treatment) and 2 (placebo). For
each patient, the 4 study visits are shown in chronological order.
Cell means model
We define a simple cell means fixed effect structure. This is the
parameterization where we want to do inference.
formula_cell_meansWe obtain the
brmsmodel matrix for the cell means parameterization.colnames(model_matrix_cell_means)Below, the first four rows are for patient 1 (treatment), and the next
four rows are for patient 2 (placebo).
Reparameterize to differences
Now we create an 8x8 transformation matrix to map the cell means model
matrix to a time difference model matrix. We will put an informative
prior on the time difference parameterization.
transformThe model matrix for time differences is:
Within each treatment group, each visit is modeled as a cumulative sum
of consecutive time differences. The row names show factors from the
data, and the column names are model terms.
We can recover the original model matrix by inverting the
transformation. This will be useful later when we derive posterior
samples of cell means from posterior samples of time differences.
Model the differences
We define an alternative dataset with the new parameterization and
subject visit indicators.
We define an MMRM model formula using the difference terms as fixed
effects and subject visit indicators in the correlation structure.
formula_differencesWe manually set a prior on each of the treatment differences. In
practice, this will be informed by historical data on the specific
treatment differences and not simply standard normal.
prior_differencesWe fit an MMRM using this special formula, data, and prior.
Samples of cell means
Extract the posterior draws of the differences from the model.
draws_differencesWe invert the transformation to derive posterior samples of cell means.
draws_cell_meansIf we had covariates and/or a non-cell-means parameterization to begin with, we could apply the custom marginal means technique
brms.mmrmalready implements for #53.Challenges
It is easy to work with the cell means parameterization, but what if the user begins with a different base model? And with baseline covariates and interactions with baseline, how do we avoid implicitly conditioning on the a reference level at a subset of the data? Figuring out these questions is the next step of this prototype. I hope a combination of workarounds, guardrails, and compromises can give us something powerful and useful for
brms.mmrm.Beta Was this translation helpful? Give feedback.
All reactions