Contabilitatea legăturilor într -un model de pericole proporționale Bayesian

Care este problema?

Problema fundamentală este că probabilitatea din modelul original presupune că timpul evenimentului este unic pentru fiecare individ. Această presupunere este rezonabilă atunci când timpul este măsurat în ore, dar devine problematic atunci când utilizați zile și cu atât mai mult cu săptămâni, mai ales dacă studiul acoperă un interval de timp larg.

Atunci când mai mulți indivizi experimentează un eveniment la același timp înregistrat (legături), provocarea este definirea „setului de risc” corespunzător – grupul încă în pericol la acel moment. Două abordări utilizate frecvent pentru manipularea legăturilor sunt Breslow şi Efron metode.

• Breslow Metoda adoptă o abordare simplă presupunând că toate evenimentele legate au același set de riscuri. Le tratează ca și cum s -ar întâmpla secvențial, dar aplică același risc stabilit pentru fiecare eveniment. Acest lucru poate funcționa bine atunci când legăturile sunt rare, dar pot introduce prejudecăți dacă sunt frecvente.

• Efron Metoda rafinează acest lucru prin ajustarea dinamic a setului de riscuri. În loc să trateze toate evenimentele legate ca având loc cu seturi de risc complet, reduce setul de risc în mod treptat pe măsură ce se întâmplă fiecare eveniment. Acest lucru se apropie mai bine de un scenariu în care evenimentele apar cu adevărat într -o succesiune strânsă, mai degrabă decât simultan.

În termeni practice, Efron Metoda oferă o corecție mai exactă atunci când legăturile sunt comune, în timp ce Breslow rămâne o alegere mai simplă din punct de vedere calcul. O altă opțiune este metoda exactă, care calculează probabilitatea, luând în considerare toate ordinea posibilă a evenimentelor legate. În timp ce această abordare este cea mai precisă, aceasta este adesea imposibilă din punct de vedere calcul pentru seturi de date mari. O mare parte din aceasta este descrisă frumos în Hertz-Picciotto și Rockhill, deși metodele originale sunt detaliate de Efron și Breslow. În cele din urmă, aceste note de prelegere ale lui Patrick Breheny oferă o explicație plăcută a algoritmilor pentru manipularea timpilor de supraviețuire legați.

Implementarea metodei EFRON

De când Efron Metoda oferă, în general, estimări mai bune, am ales să o încorporez în modelul Bayesian. Probabilitatea parțială în cadrul acestei abordări este

( log l ( beta) = sum_ {j = 1}^{j} left ( sum_ {i in d_j} x_i^ top beta- sum_ {r = 0}^{d_j-1} log left ( sum_ {k in r_j} exp (x_k^ top beta r cd 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} ) reprezintă greutatea medie de risc a persoanelor care se confruntă cu un eveniment la (t_j ):

( bar {w} = frac {1} {d_j} sum_ {i in d_j} exp (x_i^ top beta) )

Ideea cheie aici este că, în loc să tratăm toate evenimentele legate, care au loc simultan cu setul de risc complet (ca în metoda lui Breslow), Efron Reduce treptat riscul stabilit pe măsură ce fiecare eveniment legat este luat în considerare. Aceasta oferă o aproximare mai rafinată a adevăratei probabilități, în special atunci când legăturile sunt comune.

În cazul în care un exemplu numeric simplu este util, iată un exemplu de jucărie de 15 persoane, unde trei împărtășesc un timp comun de 14 zile (arătat ca (Y )), iar două împărtășesc un timp de eveniment de 25 de zile. Caracteristica critică a Efron Corecția este că ( bar {w} ) este calculat prin medie ( text {exp} ( theta) ) în grupul care se confruntă cu evenimentul la un moment dat.

Inițial, suma cumulată a ( text {exp} ( theta) ) Pentru fiecare individ din grup este identic. Cu toate acestea, pe măsură ce fiecare eveniment din grupul legat este procesat secvențial, setul de riscuri este actualizat dinamic, iar contribuția persoanelor care au fost deja contabilizate este redusă treptat. Acest lucru se reflectă în termen (Z = uv )unde (U ) reprezintă greutatea inițială a riscului total și (V ) Reprezintă reducerea incrementală pe măsură ce legăturile sunt procesate.

Simularea datelor

Pentru a genera date pentru testarea modelului Bayesian, simulez un RCT cu 1.000 de persoane randomizate 1: 1 la unul dintre cele două grupuri (reprezentat de (O)) Efectul de tratament este definit de raportul de pericol ( delta_f ). Datele includ cenzurarea și legăturile.

Începem cu încărcarea bibliotecilor necesare:

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

Iată simularea – definițiile, parametrii și generarea de date.

#### Data definitions

def <- defData(varname = "A", formula = "1;1", dist = "trtAssign")

defS <-
  defSurv(
    varname = "timeEvent",
    formula = "-11.6 + ..delta_f * A",
    shape = 0.30
  )  |>
  defSurv(varname = "censorTime", formula = -11.3, shape = .35)

#### Parameters

delta_f <- log(1.5)

#### Generate single data set

set.seed(7398)

dd <- genData(1000, def)
dd <- genSurv(dd, defS, timeName = "Y", censorName = "censorTime", digits = 0,
              eventName = "event", typeName = "eventType", keepEvents  = TRUE)

Graficul Kaplan -Meier arată cele două brațe – roșu este intervenția, iar verde este controlul, iar crucile negre sunt timpii de cenzurări.

Din cele 48 de timpi unici de eveniment, 44 sunt împărtășite de mai multe persoane. Această histogramă arată numărul de evenimente la fiecare punct:

Modelul Stan

Codul Stan de mai jos extinde modelul pentru a analiza datele de supraviețuire pe care le -am prezentat mai devreme, atât de mult din model rămâne același. Am eliminat funcția de căutare binară care a apărut în modelul anterior, înlocuindu -l cu un câmp index care este trecut de R. Există cerințe suplimentare de date care trebuie trimise și de la R Pentru a oferi informații despre legături. Modelul nu are parametri suplimentari. Reglarea probabilității pentru legături are loc în blocul „model”.

Acest cod Stan extinde modelul de supraviețuire pe care l -a descris anterior, o mare parte din structură rămânând neschimbată. Funcția de căutare binară din modelul anterior a fost eliminată și înlocuită cu un câmp index trecut de R. Sunt necesare, de asemenea, noi câmpuri de date legate de legături, care trebuie calculate în R și furnizat modelului. Nu există parametri suplimentari, iar ajustarea probabilității pentru legături este gestionată în blocul „model”.

stan_code <-
"
data {

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

  int N;          // Number of total observations
  vector(N) x;             // Covariates for all observations
  
  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 {
  
  real beta;          // Fixed effects for covariates

}

model {
  
  // Prior
  
  beta ~ normal(0, 4);
  
  // 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) = x(i) * beta;  
  }

  // 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));
  }
  
}

generated quantities {
  real exp_beta = exp(beta);
}
"

stan_model <- cmdstan_model(write_stan_file(stan_code))

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 <- dx(, A)

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(
  N_o = N_obs,
  i_o = i_obs,
  N = N_all,
  x = x_all,
  index = dx$index,
  T = nrow(ties),
  G = max(ties$grp),
  t_grp = ties$grp,
  t_index = ties$index,
  t_adj = ties$adj
)

fit_mle <- stan_model$optimize(data=stan_data, jacobian = FALSE)
fit_mle$draws(format="df")(1,c("beta", "exp_beta"))
## # A tibble: 1 × 2
##    beta exp_beta
##       
## 1 0.401     1.49

cox_model <- coxph(Surv(Y, event) ~ A , data = dd, ties = "efron")
## # A tibble: 1 × 5
##   term  estimate std.error statistic      p.value
##                         
## 1 A        0.401    0.0715      5.62 0.0000000194