Un model de pericole proporționale bayesiene pentru un studiu randomizat în cluster

Simularea datelor cu un efect aleatoriu specific clusterului

Iată R Pachete utilizate în postare:

library(simstudy)
library(ggplot2)
library(data.table)
library(survival)
library(survminer)
library(cmdstanr)

Setul de date simulează un studiu randomizat în cluster în care site -urile sunt randomizate într -un (1: 1 ) raport la grupul de tratament ( (A = 1 )) sau grup de control ( (A = 0 )) Pacienții sunt afiliați cu site-uri și primesc intervenția pe baza randomizării la nivel de sit. Rezultatul timpului de eveniment, (Y )este măsurat la nivel de pacient și depinde atât de atribuirea de tratament a site -ului, cât și de efectele sitului nemăsurat:

defC <- 
  defData(varname = "b", formula = 0, variance = "..s2_b", dist = "normal") |>
  defData(varname = "A", formula = "1;1", dist = "trtAssign")

defS <-
  defSurv(
    varname = "timeEvent",
    formula = "-11.6 + ..delta_f * A + b",
    shape = 0.30
  )  |>
  defSurv(varname = "censorTime", formula = -11.3, shape = .40)
delta_f <- 0.7
s2_b <- 0.4^2

Am generat un singur set de date de 50 de site -uri, cu 25 în fiecare braț. Fiecare site include 100 de pacienți. Graficul de mai jos arată curbele Kaplan-Meier specifice site-ului pentru fiecare braț de intervenție.

set.seed(1821)

dc <- genData(50, defC, id = "site")
dd <- genCluster(dc, "site", numIndsVar = 100, level1ID = "id")
dd <- genSurv(dd, defS, timeName = "Y", censorName = "censorTime",
             eventName = "event", typeName = "eventType", keepEvents = TRUE)

Modelul Bayesian

Întrucât aceasta este a patra iterație a modelului de pericole proporționale Bayesian la care am lucrat, se bazează în mod natural direct pe abordarea din cele trei postări anterioare (aici. Aici și aici). Acum, probabilitatea parțială a jurnalului este o funcție a efectului de tratament și a efectelor aleatorii specifice clusterului, date de:

( log L(beta) = sum_{j=1}^{J} left( sum_{i in D_j} left(beta A_i + b_{s(i)} right) – sum_{r=0}^{d_j-1} log left( sum_{k in R_j} left(beta A_k + b_ {s (k)} dreapta) – r cdot bar {w} _j dreapta) dreapta) \ )

• (D_j ) este setul de persoane care experimentează un eveniment la moment (t_j ).
• (R_j ) este riscul stabilit la momentul (t_j )inclusiv toate persoanele care sunt încă în pericol la acel moment.
• (d_j ) este numărul de evenimente la timp (t_j ).
• (r ) variază de la 0 la (d_j – 1 )iteraând evenimentele legate.
• ( bar {w} _j ) reprezintă greutatea medie de risc a persoanelor care se confruntă cu un eveniment la (t_j ):

( bar {w} _j = frac {1} {d_j} sum_ {i in d_j} left ( beta a_i + b_ {s (i)} dreapta) )

unde:

• (N ): numărul de observații (cenzurat sau nu)
• (A_i ): indicator binar pentru tratament
• ( delta_i ): indicator de eveniment ( ( delta_i = 1 ) Dacă a avut loc evenimentul, ( delta_i = 0 ) dacă este cenzurat)
• ( beta ): coeficient de tratament
• (b_ {s (i)} ): efect aleatoriu specific clusterului, unde (si)) este grupul de pacient (i )
• (R (t_i) ): risc stabilit la timp (t_i ) (inclusiv doar persoane cenzurate după (t_i ))

Distribuțiile anterioare asumate pentru ( beta ) Și efectele aleatorii sunt:

( begin {aligned} beta & sim n (0,4) \ b_i & sim n (0, sigma_b) \ sigma_b & sim t _ { text {student}} (df = 3, mu = 0, sigma = 2) end {aligned} )

stan_code <- 
"
data {
  
  int S;                   // Number of clusters

  int K;                   // Number of covariates
  int N_o;                 // Number of uncensored observations
  array(N_o) int i_o;               // Index in data set

  int N;                   // Number of total observations
  matrix(N, K) x;                   // Covariates for all observations
  array(N) int s;  // Cluster for each observation
  
  array(N) int index;
  
  int T;            // Number of records as ties
  int G;            // Number of groups of ties
  array(T) int t_grp;        // Indicating tie group
  array(T) int t_index;      // Index in data set
  vector(T) t_adj;           // Adjustment for ties (efron)

}

parameters {
  
  vector(K) beta;          // Fixed effects for covariates
  vector(S) b;             // Random effects
  real sigma_b;   // Variance of random effect
  
}

model {
  
  // Prior
  
  beta ~ normal(0, 4);
  b ~ normal(0, sigma_b);
  sigma_b ~ student_t(3, 0, 2);
  
  // Calculate theta for each observation to be used in likelihood
  
  vector(N) theta;
  vector(N) log_sum_exp_theta;
  vector(G) exp_theta_grp = rep_vector(0, G);
  
  int first_in_grp;

  
  for (i in 1:N) {
    theta(i) = dot_product(x(i), beta) + b(s(i));  
  }
  
  // Computing cumulative sum of log(exp(theta)) from last to first observation
  
  log_sum_exp_theta(N) = theta(N);
  
  for (i in tail(sort_indices_desc(index), N-1)) {
    log_sum_exp_theta(i) = log_sum_exp(theta(i), log_sum_exp_theta(i + 1));
  }

  // Efron algorithm - adjusting cumulative sum for ties
  
  for (i in 1:T) {
    exp_theta_grp(t_grp(i)) += exp(theta(t_index(i)));
  }

  for (i in 1:T) {
  
    if (t_adj(i) == 0) {
      first_in_grp = t_index(i);
    }

    log_sum_exp_theta(t_index(i)) =
      log( exp(log_sum_exp_theta(first_in_grp)) - t_adj(i) * exp_theta_grp(t_grp(i)));
  }
  
  // Likelihood for uncensored observations

  for (n_o in 1:N_o) {
    target += theta(i_o(n_o)) - log_sum_exp_theta(i_o(n_o));
  }
}
"

Pregătirea datelor pentru a trece la Stancompilând Stan cod și eșantionare din model folosind cmdstanr:

dx <- copy(dd)
setorder(dx, Y)
dx(, index := .I)

dx.obs <- dx(event == 1)
N_obs <- dx.obs(, .N)
i_obs <- dx.obs(, index)

N_all <- dx(, .N)
x_all <- data.frame(dx(, .(A)))
s_all <- dx(, site)

K <- ncol(x_all)                 # num covariates - in this case just A
S <- dx(, length(unique(site)))

ties <- dx(, .N, keyby = Y)(N>1, .(grp = .I, Y))
ties <- merge(ties, dx, by = "Y")
ties <- ties(, order := 1:.N, keyby = grp)(, .(grp, index))
ties(, adj := 0:(.N-1)/.N, keyby = grp)

stan_data <- list(
  S = S,
  K = K,
  N_o = N_obs,
  i_o = i_obs,
  N = N_all,
  x = x_all,
  s = s_all,
  index = dx$index,
  T = nrow(ties),
  G = max(ties$grp),
  t_grp = ties$grp,
  t_index = ties$index,
  t_adj = ties$adj
)

# compiling code

stan_model <- cmdstan_model(write_stan_file(stan_code))

# sampling from model

fit <- stan_model$sample(
  data = stan_data,
  seed = 1234, 
  iter_warmup = 1000,
  iter_sampling = 4000,
  chains = 4,
  parallel_chains = 4,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## Chain 2 finished in 40.0 seconds.
## Chain 3 finished in 40.0 seconds.
## Chain 1 finished in 40.3 seconds.
## Chain 4 finished in 40.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 40.2 seconds.
## Total execution time: 40.6 seconds.

Media posterioară pentru ( beta )efectul tratamentului este destul de aproape de valoarea „adevărată” de 0,70, la fel și estimarea abaterii standard a efectului aleatoriu (am folosit (SD = 0,4 ) În procesul de generare a datelor):

fit$summary(variables = c("beta", "sigma_b"))
## # A tibble: 2 × 10
##   variable  mean median     sd    mad    q5   q95  rhat ess_bulk ess_tail
##                        
## 1 beta(1)  0.659  0.660 0.124  0.124  0.455 0.863  1.00    1529.    3265.
## 2 sigma_b  0.417  0.413 0.0472 0.0458 0.347 0.501  1.00   15493.   11171.

Postarea finală din această serie va include cod pentru a simula datele dintr-un proces-randomizat cu cluster în trepte, cu un rezultat de timp pentru eveniment. Acest model va integra atât componentele spline, cât și cele aleatorii. Sunt curios să văd cât de bine funcționează, deoarece resursele de calcul necesare ar putea fi substanțiale.