Assessing correlation/surrogacy between 2 estimated treatment effects

Ahmed_Sayed · January 31, 2024, 2:06pm

Hello everyone,

I’m working on a dataset of N trials containing treatment effects for outcomes X and Y, the standard errors of these treatment effects, and the within-trial correlation of treatment effects on X and Y.

Y is generally considered more important/relevant for patients and clinicians than X but, due to restraints on time & resources, most trials focus on X. I am therefore interested in analyzing to what extent a treatment’s effect on X is a good surrogate/predictor of its effect on Y.

Ideally, I would like my analysis to take into account both the uncertainties involved in estimating X and Y as well as their within-trial correlation. An example of a paper that (I think) does this: https://pubmed.ncbi.nlm.nih.gov/30635226/ (the associated statistical appendix clarifies their methodological approach, although it is behind a paywall).

Now, onto a simulated data example and the approaches I’ve tried to tackle this:

#Load foreach and data.table
require("foreach")
require("data.table")

#Sequence over range of (somewhat) plausible surrogate treatment effects "stes"
stes <- seq(-0.5, 0.5, 0.05)

#Run a loop over stes
tes <- foreach(ste = stes,
        counter_id = 1:length(stes), #Add a loop counter as a trial id
        .combine = "rbind") %do% {
          
          #Sample some amount of error from 0.1 to 1
          error <- sample(seq(0.1, 1, 0.1), 1)
          
          #Simulate a distribution of surrogate treatment effects
          #giving rise to the observed treatment effect in this trial
          surrogate_dist <- rnorm(n = 1000,
                                  mean = ste, #With a given mean
                                  sd = error) #With a given standard deviation
          
          #Sample some amount of within-trial correlation of treatment effects from 0.1 to 1
          rho <- sample(seq(0.1, 0.8, 0.1), 1)
          
          #Simulate a distribution of true treatment effects
          #giving rise to the observed treatment effect in this trial
          true_dist <- (rho * surrogate_dist) + sqrt(1 - rho*rho) * rnorm(n = 1000, mean = 0, sd = error)
          
          #Create a data.table representing this trial
          data.table(trial_id = counter_id,
                     surrogate_te = mean(surrogate_dist),
                     true_te = mean(true_dist),
                     surrogate_sete = sd(surrogate_dist),
                     true_sete = sd(true_dist),
                     intra_cor = rho)
        }

We now have a dataset of 21 trials, each with a treatment effect on X and Y, standard errors for these treatment effects, and a within-trial correlation coefficient. The effects on X and Y vary between the trials (reflecting the different degrees of effectiveness of treatments), as do the associated standard errors (reflecting the different degrees of precision due to variation in trial size).

Based on some of the previous literature (https://10.0.3.233/jamainternmed.2021.5726 among other examples), my first impulse was to just run a weighted linear regression. I set the weights as being equal to the inverse of the sum of the variances of treatment effects (such that we let more precise estimates shape our results).

#Run weighted least squares with weights equal to the precision of the observed
#treatment effect
#Run weighted least squares with weights equal to the precision of the observed
#treatment effect
weighted_model <- lm(data = tes, 
                     formula = true_te ~ surrogate_te, 
                     weights = 1/(surrogate_sete^2 + true_sete^2))

The 2 problems with this are (I think) that it ignores the absolute uncertainty in estimating the treatment effects (the weighting procedure only utilizes the relative precision of each trial’s estimate but not the absolute magnitude of their standard errors) and it also ignores within-trial correlations.

Other approaches I’ve tried include:

1) Incorporating error in the treatment effects using the se and me functions of brms:

require("parallel")
#Model 1
m1 <- brm(data = tes,
          prior = prior(normal(0, 1), class = "b"),
          cores = parallel::detectCores() - 1,
          formula = true_te | se(true_sete) ~ me(surrogate_te, surrogate_sete)
          )

As I understand it, this models the actual treatment effects on X and Y as arising from normal distributions, each with a given (observed) mean and standard deviation (equal to the observed standard error). This looks like it accounts for uncertainty in treatment effects, but it ignores within-trial correlation.

2) Using a bivariate normal and modeling the correlation between the residuals:

#Model 2
m2 <- brm(data = tes,
                prior = prior(normal(0, 1), class = "Intercept"),
                cores = parallel::detectCores() - 1,
                formula = bf(true_te | se(true_sete) ~ 1) +
            bf(surrogate_te | se(surrogate_sete) ~ 1) +
            set_rescor(TRUE)
          )

This also looks like it appropriately models observed treatment effects as arising from an unobserved distribution with a given mean and standard error. But again, this would not account for within-trial correlation.

3) Reframing the problem and viewing the two types of outcomes (surrogate and true) as being random effects nested within trials and using the correlation parameter between them to assess surrogacy:

#Convert tes to long format first
tes_long <- data.table(
  te = tes[, .(c(true_te, surrogate_te)), by = trial_id][[2]],
  se = tes[, .(c(true_sete, surrogate_sete)), by = trial_id][[2]],
  outcome = rep(c("True", "Surrogate"), times = nrow(tes)),
  trial_id = rep(tes$trial_id, each = 2),
  intra_cor =  rep(tes$intra_cor, each = 2)
)

#Create a variance-covariance matrix
cov_matrix <- matrix(nrow = nrow(tes_long),
                     ncol = nrow(tes_long)
)

#Fill up the diagonals with the variance for surrogate and true treatment effects
diag(cov_matrix) <- tes_long$se^2

#Calculate and fill in the covariance of all-cause and HF hospitalizations for each trial
for(i in  1:uniqueN(tes_long$trial_id)) {
  cov_matrix[2*i, (2*i) - 1] <- tes_long[trial_id == i, intra_cor][1]*sqrt(cov_matrix[(2*i) - 1, (2*i) - 1])*sqrt(cov_matrix[2*i, 2*i])
  cov_matrix[(2*i) - 1, 2*i] <- tes_long[trial_id == i, intra_cor][1]*sqrt(cov_matrix[(2*i) - 1, (2*i) - 1])*sqrt(cov_matrix[2*i, 2*i])
}

#Convert NAs to 0
cov_matrix[is.na(cov_matrix)] <- 0

#Model 3
m3 <- brm(data = tes_long,
          data2 = list(cov_matrix = cov_matrix),
          prior =
            c(prior(constant(1), class = "sigma")),
          cores = parallel::detectCores() - 1,
          formula = te ~ 
            0 + #No global intercept
            outcome + #A separate intercept for each of the 2 outcomes
            (0 + outcome|trial_id) + #Nest the 2 outcomes within their trials
            fcor(cov_matrix) #Pass the covariance matrix
)

This last option seems like it takes both the standard errors and within-trial correlations into account since I’m allowed to pass a variance-covariance matrix containing said information.

That said, given my lack of formal training, I would greatly appreciate any help regarding what approach(es) may be best suited for this task (and what the key differences between them might be).

florianstijven · February 7, 2024, 7:44pm

Hello Ahmed,

There is extensive literature on the (evaluation) of surrogate endpoints. Your proposed analysis is part of the so-called meta-analytic framework for evaluating surrogate endpoints. In this framework, we distinguish between surrogacy at two levels:

Trial-level surrogacy: This concerns the statistical association between trial-level treatment effects on X and Y, and is quantified by an R-squared measure.
Individual-level surrogacy: This concerns the statistical association between X and Y after correcting for the trial- and treatment effects. Since you do not have individual patient data (at least I suppose) you cannot evaluate this type of surrogacy.

See for example this early paper that introduced the meta-analytic framework. This more recent book describes how surrogate endpoints can be evaluated within this framework, for almost any possible type of endpoints (continuous, binary, survival, …).

The approach you want to apply was also developed some time ago. This paper contains a practical application of a method to correct for (i) the standard errors in the estimated treatment effects and (ii) the within-trial correlation. The method to correct for these different types of uncertainty is described in more detail in van Houwelingen et al. (2002).
Your third model corresponds to this approach I think, but I am not 100% sure. I should look into your third model in some more detail to confirm this.

Kind regards,
Florian

Ahmed_Sayed · February 10, 2024, 4:42am

Hello Florian,

Thank you for the elaborate reply and helpful references!

I’ve gone through them and, after reading Chapter 4.2 (“The Meta-Analytic Framework”) in the referenced book, I’ve fit models similar to those in equations 15 and 16 in the paper referenced by Buyse et al. (https://pubmed.ncbi.nlm.nih.gov/12933525/).

I was able to use a binomial distribution because, although I don’t have patient-level data, the outcomes are binary and nested within one another (that is, X is a subset of Y, thereby allowing me to reconstruct the data for each of the included trials). The results were relatively similar to the 3rd model in the original post where I worked on pre-computed effect sizes. Still, I think this is probably a bit more elegant in that it’s a 1-stage approach and I don’t have to bootstrap.

The one difference between the models I used and those proposed by Buyse et al. is that I modeled the baseline risk in each trial as a fixed effect rather than as random intercepts like so:

bform <- bf(s_events | trials(n) #Surrogate outcome
            ~ 
              0 + #No intercept
              trial + #Model trial-specific baseline risk as a fixed effect
              ttt + #Model the average effect of treatment on surrogate outcomes
              (0 + ttt|p|trial)) + #Add random effects to model variation in treatment effects on the surrogate outcome
  bf(t_events | trials(n) #Repeat but for the true outcome
     ~
       0 + trial + ttt + (0 + ttt|p|trial))
#The "p" in the above formula tells brms to model the correlation between treatment effects on the surrogate and true outcomes

I did this because (please correct me if I’m wrong) the use of trial random intercepts as suggested by Buyse and colleagues may introduce some (probably small) amount of bias in estimated treatment effects if A) treatment allocations were unequal and B) baseline risks vary across trials (https://onlinelibrary.wiley.com/doi/10.1002/bimj.200900074).

In case anyone needs to do something similar using brms, I’ve attached an example below
where I simulate some data, fit the model, and compute some quantities that may be of interest in an analysis like this.

#Load data.table
require("data.table")
require("dplyr")
#We will generate 40 two-armed trials.
#First, simulate a range of different trial baseline risks and trial sizes
trial_s_baseline_risk <- runif(40, 0.1, 0.4) #Baseline risk for surrogate outcome
trial_t_baseline_risk <- runif(40, 0.4, 0.7) #Baseline risks for true outcome
trial_n <- sample(1000, 40) #N of patients per trial

#Put everything into a data.table
bdata <- data.table(
  ttt = c(0, 1) %>% rep(each = 40), #Treatment arms
  trial = (1:40) %>% rep(times = 2) %>% factor, #Trial ID
  n = trial_n %>% rep(2), #Assume equal N in control/treatment arm
  s_events = rbinom(80, trial_n %>% rep(2), trial_s_baseline_risk), #N of patients experiencing a surrogate event
  t_events = rbinom(80, trial_n %>% rep(2), trial_t_baseline_risk) #N of patients experiencing a true event
  
)

#Load brms
require("brms")
#Build our formula
bform <- bf(s_events | trials(n) #Surrogate outcome
            ~ 
              0 + #No intercept
              trial + #Model trial-specific baseline risk as a fixed effect
              ttt + #Model the average effect of treatment on surrogate outcomes
              (0 + ttt|p|trial)) + #Add random effects to model variation in treatment effects on the surrogate outcome
  bf(t_events | trials(n) #Repeat but for the true outcome
     ~
       0 + trial + ttt + (0 + ttt|p|trial))
#The "p" in the above formula tells brms to model the correlation between treatment effects on the surrogate and true outcomes


#Fit the binomial model
mb <- brm(data = bdata, #Our data
          family = binomial, #Use of a binomial
          #Set some reasonable priors based on subject matter input 
          prior = c(prior(cauchy(0, 1), class = "sd", resp = "sevents"),
                    prior(cauchy(0, 1), class = "sd", resp = "tevents"),
                    prior(normal(0, 1), class = "b", coef = "ttt", resp = "tevents"),
                    prior(normal(0, 1), class = "b", coef = "ttt", resp = "sevents"),
                    prior(normal(0, 2), class = "b")),
          formula = bform, #Our formula
          #N of cores & threads (adjust as needed)
          cores = 16,
          threads = threading(threads = 4, grainsize = NULL, static = FALSE),
          iter = 3000
)

#Perform posterior predictive checks
pp_check(mb, resp = "sevents")
pp_check(mb, resp = "tevents")
#Check that Rhat values are < 1.05
all(rhat(mb) < 1.05)
#Do other diagnostics as needed (e.g., inspection of traceplots)

#Load tidybayes
require("tidybayes")

#Store the draws from the posterior
pdraws <- tidy_draws(mb)

#Now, let us calculate the slope of the relationship
#Since we modeled our treatment effects on the 2 outcomes as arising from a bivariate normal, 
#we can calculate the slope as the correlation coefficient * the ratio of the std. deviation parameters of the two outcomes:
slope <- pdraws$cor_trial__sevents_ttt__tevents_ttt*(pdraws$sd_trial__tevents_ttt/pdraws$sd_trial__sevents_ttt)

#We will also calculate R-squared (the square of our correlation coeffecient)
rsq <- pdraws$cor_trial__sevents_ttt__tevents_ttt^2

#We can also calculate the RMSE, which will be helpful for making predictions a little later.
#The rmse is calculated by multiplying the std. deviation parameter of the true outcome by the square root of 1 - R2
rmse <- pdraws$sd_trial__tevents_ttt*sqrt(1 - rsq)

#Calculate the intercept (average effect on the true outcome minus the average effect on the surrogate outcome multiplied by the slope)
intercept <- (pdraws$b_tevents_ttt - (slope * pdraws$b_sevents_ttt))

#We can make predictions for a trial with a given surrogate treatment effect.
#First, we need to specify the following:
#1) The observed treatment effect on the surrogate (with an associated standard error representing our uncertainty in this observed effect)
observed_ste <- rnorm(n = 12000, #Same N of draws as the posterior)
                      mean = -1,
                      sd = 0.1)

#2) An element of variation centered on 0 with an SD equal to our RMSE
#In essence, this represents the uncertainty that persists after knowing the effect on the surrogate
#In a sense, this is effect of the treatment on the true outcome that is not explained (?mediated?) by its effect
#on the surrogate outcome
residual_uncertainty <- rnorm(n = 12000, #Same N of draws as the posterior
                     mean = 0, #A mean of zero
                     sd = rmse) #
  
#The predicted effect on the true outcome is calculated as
predicted_tte <- intercept + #Our intercept
  observed_ste * slope + #The slope multiplied by the observed surrogate effect
  residual_uncertainty #The residual uncertainty that remains 

#From this, we could obtain several values of interest
predicted_tte %>% density %>% plot #Plot the predicted distribution of treatment effects on the true outcome given the observed surrogate treatment effect and its uncertainty
predicted_tte %>% quantile(c(0.025, 0.975)) #Obtain a 95% prediction interval
sum(predicted_tte < 0)/length(predicted_tte) #Compute the probability that the predicted effect is protective

#Suppose we had a trial with a given sample size and N of surrogate events,
#which results in a particular standard error for the treatment effect on the surrogate outcome.

#What is the minimum protective treatment effect on the surrogate outcome that would need to be observed
#so that the posterior probability of benefit (protection) on the true outcome is 97.5%?

#First, set a sequence over a range of plausible surrogate treatment effects
stes_range <- seq(-2, 2, 0.01)
#And set a standard error implied by the N of events in the trial
implied_standard_error <- 0.05

#Load the foreach package
require("foreach")

#Setup a foreach loop
preds_df <- foreach(ste = stes_range,
                    .combine = "rbind") %do% {
  
          #Generate treatment effect from hypothetical trial
          observed_ste <- rnorm(n = 12000, #Same N of draws as the posterior)
                                mean = ste,
                                sd = implied_standard_error)
          
          #The predicted effect on the true outcome is calculated as
          predicted_tte <- intercept + #Our intercept
            observed_ste * slope + #The slope multiplied by the observed surrogate effect
            residual_uncertainty #The residual uncertainty that remains 
          
          #Put the surrogate treatment effect, the point estimate of the true effect, and the 2.5th/97.5th percentiles into a data.table
          data.table(surrogate_te = ste,
                     pe = median(predicted_tte),
                     q2.5 = quantile(predicted_tte, c(0.025)),
                     q97.5 = quantile(predicted_tte, c(0.975))
                     )
                    }

require("ggplot2")
#Plot our prediction and observe the location where the 95% predictive interval crosses 0
ggplot(data = preds_df,
       aes(x = surrogate_te,
           y = pe,
           ymax = q97.5,
           ymin = q2.5)) +
  geom_line(color = "maroon", linewidth = 1) +
  geom_ribbon(alpha = 0.5, fill = "lightblue") +
  theme_bw()
#(in this example) No such location exists since event risks were sampled independent of treatment allocation.

Thanks again for your help!

Sincerely,
Ahmed

f2harrell · February 10, 2024, 12:10pm

This is a very rich discussion. One tangential comment. Sometimes I find it handy to think of surrogate outcomes as a multiple imputation question. If we have a large number of persons where we assess both the gold standard outcome and the surrogate outcome we can use multiple imputation to impute the gold standard outcome from the surrogate outcome. If you applied the multiple imputation patterns to a new dataset with no gold standard measurements you would be able to make inferences about treatment effects on the gold standard. But you would find that power of the treatment comparison may actually be less than had you used the gold standard in the data at hand. This will seem confusing but what I mean is for example the gold standard is overall survival that takes > 5 years of follow-up and your study is a 2y study. The limited number of all-cause mortalities in the 2y study will limit the power of the treatment comparisons. But it may still have more power than using the surrogate to impute the real outcomes.

Ahmed_Sayed · February 10, 2024, 1:46pm

This is a really helpful analogy. In this case, would it be correct to say that the distribution of the variable being imputed (the true treatment effect), and all of its statistical properties, is some function of:

1) The uncertainty of the variable we’re using to impute (i.e., standard error of the surrogate treatment effect in the new trial).

2) The inherent extent of variation in the variable being imputed (i.e., the heterogeneity parameter of the true outcome effect in meta-analytic models of trials where you had both measurements)

3) The strength of correlation between the variable being used to impute and the variable being imputed (i.e., the correlation parameter between treatment effects in the model).

The combination of 2. (heterogeneity) and 3. (correlation) results in some amount of residual uncertainty (represented by the RMSE in the above code) that is distinct from the variability in 1. (error in estimated surrogate effect) and which would persist even if 1. were reduced to zero.

So it boils down to 2 questions:

A) Is your estimated surrogate treatment effect precise and large enough to outweigh this residual uncertainty?

B) Is the cost of obtaining this much precision of your surrogate treatment effect estimate less than the cost needed to precisely estimate the true treatment effect?

f2harrell · February 10, 2024, 11:02pm

I like you way you’ve summarized this.

florianstijven · February 12, 2024, 4:43pm

Hello Ahmed,

In your latest model formula (repeated below), you assume that the outcomes X and Y are independent within trials. Although you simulated the data assuming independence between X and Y within trials, X and Y are generally associated when evaluating X as a surrogate for Y; this is the individual-level surrogacy mentioned before. This residual association should be modeled to obtain correct inferences.

bform <- bf(s_events | trials(n) #Surrogate outcome
            ~ 
              0 + #No intercept
              trial + #Model trial-specific baseline risk as a fixed effect
              ttt + #Model the average effect of treatment on surrogate outcomes
              (0 + ttt|p|trial)) + #Add random effects to model variation in treatment effects on the surrogate outcome
  bf(t_events | trials(n) #Repeat but for the true outcome
     ~
       0 + trial + ttt + (0 + ttt|p|trial))

To my knowledge, it is impossible to take this residual association between X and Y into account using brms. Indeed, multivariate logistic regression models cannot be fitted using brms. In a multivariate logistic regression model, we could model the association between X and Y (within trials) through the odds ratio, which also quantifies the individual-level surrogacy.

Still, we could model this residual association in brms using shared patient-level random intercepts, although I find this solution less satisfactory. There may also be other, and better, approaches to take this residual association into account. Below, I simulate some patient-level data and fit your model and a model with patient-level random intercepts. Only the latter model takes the residual association into account. The credible intervals for the correlation between the treatment effects differ between these two models. I expect that this difference will increase with an increasing within-trial association and decreasing trial sizes.

library(tidyverse)
library(brms)
set.seed(1)
# Number of trials.
N_trial = 20
# Number of subjects in each trial in each arm.
N_subjects = 30
# Conditional log odds ratio between X and Y within trials.
log_odds_ratio = 2
expit = function(x) 1 / ( 1 + exp(-x))

data = tibble(
  # Generate trial-level logistic regression parameters.
  surr_intercept_trial = rnorm(n = N_trial, mean = -1),
  surr_effect_trial = rnorm(n = N_trial, mean = 2),
  true_intercept_trial = rnorm(n = N_trial, mean = -2),
  # The treatment effects on the surrogate and true endpoint are correlated.
  true_effect_trial = surr_effect_trial + rnorm(n = N_trial, mean = 0, sd = 0.5),
  trial = 1:N_trial
) %>%
  # Compute the trial-level marginal probabilities for X and conditional
  # probabilities for Y. 
  mutate(
    pi_surr_marg_control = expit(surr_intercept_trial),
    pi_true_cond_surr1_control = expit(true_intercept_trial + log_odds_ratio),
    pi_true_cond_surr0_control = expit(true_intercept_trial),
    pi_surr_marg_treated = expit(surr_intercept_trial + surr_effect_trial),
    pi_true_cond_surr1_treated = expit(true_intercept_trial + log_odds_ratio + true_effect_trial),
    pi_true_cond_surr0_treated = expit(true_intercept_trial + true_effect_trial)
  ) %>%
  # Given the trial-level parameters, we generate patient-level data. 
  rowwise(everything()) %>%
  reframe(bind_rows(
    tibble(
      # Generate surrogate outcomes.
      surr = rbinom(n = N_subjects, size = 1, prob = pi_surr_marg_control),
      # Generate true outcomes, given the surrogate outcome.
      true = rbinom(
        n = N_subjects,
        size = 1,
        prob = ifelse(
          surr == 1,
          pi_true_cond_surr1_control,
          pi_true_cond_surr0_control
        )
      ),
      # Control Group.
      treatment = 0,
      subject_id = trial * 2 * N_subjects + 1:N_subjects
    ),
    # Repeat the data generation for the treatment group.
    tibble(
      surr = rbinom(n = N_subjects, size = 1, prob = pi_surr_marg_treated),
      true = rbinom(
        n = N_subjects,
        size = 1,
        prob = ifelse(
          surr == 1,
          pi_true_cond_surr1_treated,
          pi_true_cond_surr0_treated
        )
      ),
      # Treatment group.
      treatment = 1,
      subject_id = trial * 2 * N_subjects + 1:N_subjects + N_subjects
    )
  )) 

data_long = data %>%
  pivot_longer(cols = c("surr", "true"), names_to = "outcome", values_to = "response")

# Model formula that assumes indepndence between X and Y within subjects.
bform_0 <-
  bf(surr ~ 0 + as.factor(trial) + treatment + (0 + treatment | p | trial)) +
  bf(true ~ 0 + as.factor(trial) + treatment + (0 + treatment | p | trial))

# Model formula that allows for within-subject association between X and Y.
bform_1 <-
  bf(
    response ~ 0 + as.factor(trial):outcome + treatment:outcome +
      (0 + as.factor(outcome):treatment | trial) + (1 | subject_id)
  )
  
# Model that assumes indepndence between X and Y within subjects.
mb_0 <- brm(
  data = data,
  family = bernoulli(),
  prior = c(prior(normal(2, 10), class = "sd", resp = "true", group = "trial"),
            prior(normal(2, 10), class = "sd", resp = "surr", group = "trial"),
            prior(normal(0, 5), class = "b")),
  formula = bform_0, 
  cores = 8)
summary(mb_0)

# Model formula that allows for within-subject association between X and Y.
mb_1 <- brm(
  data = data_long,
  family = bernoulli(),
  prior = c(prior(normal(2, 10), class = "sd", group = "trial"),
            prior(normal(2, 10), class = "sd", group = "subject_id"),
            prior(normal(0, 5), class = "b")), 
  formula = bform_1,
  cores = 8
)
summary(mb_1)

We should be careful with these alternative models, however, because the interpretation of the trial-level treatment effects changes when we introduce patient-level random intercepts. Indeed, the trial-level treatment effects are no longer marginal log odds ratios. In light of this, I would prefer the third model in your original post, because it models the “original” treatment effects, e.g., marginal odds ratios, risk differences, etc. Especially when you’re only interested in trial-level surrogacy.

Kind regards,
Florian

Ahmed_Sayed · February 15, 2024, 12:34pm

Hello Florian,

Thanks for getting back, pointing out the oversight on my part, and providing the code examples

Looking back on it, I think it might be possible to fit correlations between X and Y by replacing the fixed trial baseline risks with a single fixed intercept in addition to trial-level random intercepts. On brms, I believe this would be done in the following manner:

bform <- bf(s_events | trials(n) #Surrogate outcome
            ~ 
              1 + #Fixed intercept
              (1|r|trial) + #Model trial-level variation in baseline risks (with r as a correlation parameter)
              ttt + #Model the average effect of treatment on surrogate outcomes
              (0 + ttt|p|trial)) + #Add random effects to model variation in treatment effects on the surrogate outcome
  bf(t_events | trials(n) #Repeat but for the true outcome
     ~
       1 +  (1|r|trial) + ttt + (0 + ttt|p|trial))

This would give us the correlation (r) between the trial-level random intercepts of the true and surrogate outcomes (thus modeling the usually positive correlation between X and Y within trials).

This does force us to forsake the N trial-specific baseline-risk fixed effects though. Which has the potential to introduce small amounts of bias (in the treatment effect estimates) if there are imbalances in treatment allocation and baseline risks across trials.

Having looked at the Buyse et al. paper in a bit more detail, I think the last model I fit closely corresponds with the “Reduced random-effects model” they describe in the paper (just below Eq. 26) and also in Section 4.3.2 of the book. Their model is described by the two equations (4.12 in the book):

Tij = µT + (β + bi) Zij + εTi for the true treatment effect (T)
Sij = µS + (α + ai) Zij + εSij for the surrogate treatment effect (S)

And a variance-covariance matrix D (4.13 in the book):

D = \begin{pmatrix}daa & dab &\\ & dbb \end{pmatrix}

Where daa is the variance of the true treatment effect, dbb is the variance of the surrogate treatment effect, and dab is their covariance.

The two differences are:

My use of N trial-level fixed effect parameters rather than a single global intercept (uT and uS in the above models)
the lack of a residual correlation in my models.

For 1., I think the choice to use N parameters is probably justifiable in my use-case because of the known variation in baseline risks between trials (that would not be adequately captured by a global intercept).

For 2., I believe it is only necessary for individual-level surrogacy but does not bias treatment effect estimates. Also, since we are fitting logistic regression models, I am not quite sure that a correlation between residuals can be computed, since there is no real error term (at least not explicitly) in our model.

We definitely agree regarding the difficulty of interpreting models with patient-level random intercepts. Thanks a lot for your help with this!

Sincerely,
Ahmed

florianstijven · February 16, 2024, 10:44am

Hello Ahmed,

Thank you for your response! This discussion certainly also helps me improving my own understanding of these methods.

I want to clarify the residual association. Let (X, Y)' | \boldsymbol{u}, A be the conditional distiribution of the outcomes given the trial-level random effects, \boldsymbol{u}, and the treatment indicator, A. When we only include trial-level random effects and do not account for residual association, then we assume that X \perp Y | \boldsymbol{u}, A, i.e., within each trial and treatment group, X and Y are independent. Note that this residual association can be perfectly captured be the following odds ratio: \frac{P(X = 1, Y = 1 | \boldsymbol{u}, A) P(X = 0, Y = 0 | \boldsymbol{u}, A)}{P(X = 1, Y = 0 | \boldsymbol{u}, A) P(X = 0, Y = 1 | \boldsymbol{u}, A)}.

As you stated correctly, it would not make sense to compute the correlation between residuals when X and Y are binary. However, there are other approaches to modeling the residual association. First, we could introduce patient-level random effects as I outlined in my previous post. Second, we could model the odds ratio introduced in the previous paragraph; this would be a multivariate logistic regression model. As I stated in my previous post, I do not think multivariate logistic regression is implemented in brms.

I partly support your argument that the residual association is unimportant for trial-level surrogacy. However, I think this argument only holds in specific settings.
Look at equation 23 in this paper. This equation highlights the bias in estimating the trial-level surrogacy when we ignore the correlation between the trial-level treatment effects that is due to sampling variability (i.e., given that we repeat the same trial multiple times, what is the correlation between the estimated treatment effects?). This latter correlation is a consequence of the residual association. Hence, this bias is an indirect consequence of the residual association. Equation 23 shows that the bias is large when 2 conditions are simultaneously satisfied:

Strong correlation between trial-level treatment effects that is due to sampling variability. This correlation is strong when the residual association is strong.
The standard errors for the estimated treatment effects are large as compared to the standard deviation of the “true” treatment effects. This is likely to hold when the trial-level sample sizes are small.

Kind regards,
Florian

Ahmed_Sayed · February 17, 2024, 4:27pm

Hello Florian,

Thank you for the explanation and the useful reference. I was indeed mistaken in that the binomial model I fit didn’t take into account the within-trial correlation of treatment effects. As you said, such information can only be (I think) supplied to brm via the fcor argument, as in the 3rd example in the original post.

The differences are not huge in my case (likely due to the presence of a number of large trials with relatively small standard errors). As you say, I imagine in cases where small sample trials dominate, the bias may turn out to be more sizeable.

This has been a really informative discussion!

Sincerely,
Ahmed