Model de pericole proporționale bayesiene pentru un design cu țesut în trepte

(Acest articol a fost publicat pentru prima dată pe Generarea noastră de dateși a contribuit cu drag la R-Bloggers). (Puteți raporta problema despre conținutul de pe această pagină aici)

Doriți să vă împărtășiți conținutul pe R-Bloggers? Faceți clic aici dacă aveți un blog sau aici dacă nu.

Am ajuns în sfârșit la capătul drumului. Aceasta este a cincea și ultima postare dintr-o serie care se construiește până la un model de pericole proporționale Bayesian pentru analizarea unui proces-randomizat de cluster cu țesut în trepte. Dacă tocmai vă alăturați, poate doriți să începeți la început.

Modelul prezentat aici integrează tendințele de timp neliniare și efectele aleatorii specifice clusterului-elemente pe care le-am explorat anterior izolat. Nu există nimic fundamental nou în această postare; Reunește totul. Având în vedere că lucrările de bază au fost deja puse, voi păstra comentariul pe scurt și mă voi concentra pe furnizarea codului.

Simularea datelor de la un CRT cu pași în trepte

Voi genera un singur set de date pentru 25 de site-uri, fiecare site care înscrie participanți la studiu pe o perioadă de 30 de luni. Site -urile vor trece de la control la intervenție secvențial, cu un nou site începând în fiecare lună. Fiecare site va înscrie 25 de pacienți în fiecare lună.

Rezultatul ( (Y )) este numărul de zile la un eveniment. Tratamentul ((O)) reduce timpul până la eveniment. Timpurile de supraviețuire depind și de luna de înscriere – un efect pe care l -am exagerat pentru ilustrație. În plus, fiecare site (i ) are un efect specific site-ului (b_i sim n ( mu = 0, sigma = 0,5) )care influențează timpul până la eveniment printre participanții săi.

Iată bibliotecile necesare codului prezentat aici:

library(simstudy)
library(data.table)
library(splines)
library(splines2)
library(survival)
library(survminer)
library(coxme)
library(cmdstanr)

Definiții

def <- defData(varname = "b", formula = 0, variance = 0.5^2)

defS <-
  defSurv(
    varname = "eventTime",
    formula = 
      "..int + ..delta_f * A + ..beta_1 * k + ..beta_2 * k^2 + ..beta_3 * k^3 + b",
    shape = 0.30)  |>
  defSurv(varname = "censorTime", formula = -11.3, shape = 0.36)

Parametri

int <- -11.6
delta_f <-  0.80

beta_1 <-  0.05
beta_2 <-  -0.025
beta_3 <- 0.001

Generarea de date

set.seed(28271)

### Site level data

ds <- genData(25, def, id = "site")                 
ds <- addPeriods(ds, 30, "site", perName = "k") 

# Each site has a unique starting point, site 1 starts period 3, site 2 period 4, etc.

ds <- trtStepWedge(ds, "site", nWaves = 25,     
                   lenWaves = 1, startPer = 3, 
                   grpName = "A", perName = "k")

### Individual level data

dd <- genCluster(ds, "timeID", numIndsVar = 25, level1ID = "id") 
dd <- genSurv(dd, defS, timeName = "Y", censorName = "censorTime", digits = 0,
              eventName = "event", typeName = "eventType")

### Final observed data set

dd <- dd(, .(id, site, k, A, Y, event))

Iată un set de loturi Kaplan-Meier pentru fiecare site și perioada de înscriere. Când un site este în condiția de intervenție, curba KM este roșie. Pentru simplitate, cenzurarea nu este arătat, deși aproximativ 20% din cazurile din acest set de date sunt cenzurate.

Estimarea modelului

Acest model are destul de multe componente în raport cu modelele anterioare, dar nimic nu este cu adevărat nou. Există un spline penalizat pentru efectul timpului și un efect aleatoriu pentru fiecare site. Parametrul principal de interes este încă ( beta ).

Pentru completitate, iată specificația modelului:

( log l ( beta) = sum_ {j = 1}^{j} left ( sum_ {i in d_j} left ( beta a_i + sum_ {m = 1}^m gamma_m x_ {m_i} + b_ {s (i) dreapta) – log_ {r = 0 left( sum_{k in R_j} left(beta A_k + sum_{m=1} ^ M gamma_m X_{m_i} + b_{s(k)} right) – r cdot bar{w}_j right) right) – lambda sum_{m=1}^{M} left( Q^{(2)} gamma dreapta) _m^2 \ )

unde

(J ): Numărul timpului de evenimente unic
(M ): numărul de funcții de bază spline
(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) )

(A_i ): indicator binar pentru tratament
(X_ {m_i} ): valoarea (m^{ text {th}} ) funcție de bază spline pentru (i^{ text {th}} ) observare
(Q^{(2)} ): cea de-a doua matrice de diferență a funcției spline

Parametrii modelului sunt

( beta ): coeficient de tratament
( gamma_m ): coeficient spline pentru (m^ text {th} ) Funcția de bază spline
(b_ {s (i)} ): efect aleatoriu specific clusterului, unde (si)) este grupul de pacient (i )
( lambda ): termenul de penalizare; Acest lucru nu va fi estimat, ci furnizat de utilizator

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) \ gamma_m & sim n (0,2) end {aligned ) )

Și iată implementarea modelului în Stan:

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;      // Event times (sorted in decreasing order)

  int N;          // Number of total observations
  matrix(N, K) x;          // Covariates for all observations
  array(N) int s;          // Cluster
  
  // Spline-related data
  
  int Q;          // Number of basis functions
  matrix(N, Q) B;          // Spline basis matrix
  matrix(N, Q) Q2_spline;  // 2nd derivative for penalization
  real lambda;    // penalization term
  
  array(N) int index;

  int T;            // Number of records as ties
  int J;            // 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;   // SD of random effect
  
  vector(Q) gamma;         // Spline coefficients
  
}

model {
  
  // Priors
  
  beta ~ normal(0, 1);
  
  // Random effects
  
  b ~ normal(0, sigma_b);
  sigma_b ~ normal(0, 0.5);

  
  // Spline coefficients prior
  
  gamma ~ normal(0, 2);
  
  // Penalization term for spline second derivative
  
  target += -lambda * sum(square(Q2_spline * gamma));
  
  // Compute cumulative sum of exp(theta) in log space (more efficient)
  
  vector(N) theta;
  vector(N) log_sum_exp_theta;
  vector(J) exp_theta_grp = rep_vector(0, J);
  
  int first_in_grp;
  
  // Calculate theta for each observation
  
  for (i in 1:N) {
    theta(i) = dot_product(x(i), beta) + dot_product(B(i), gamma) + b(s(i));
  }
  
  // Compute cumulative sum of log(exp(theta)) from last to first observation
  
  log_sum_exp_theta = rep_vector(0.0, N);
  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));
  }
}
"

Compilarea modelului:

stan_model <- cmdstan_model(write_stan_file(stan_code))

Obținerea datelor de la R la Stan:

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)))

# Spline-related info

n_knots <- 5
spline_degree <- 3
knot_dist <- 1/(n_knots + 1)
probs <- seq(knot_dist, 1 - knot_dist, by = knot_dist)
knots <- quantile(dx$k, probs = probs)
spline_basis <- bs(dx$k, knots = knots, degree = spline_degree, intercept = TRUE)
B <- as.matrix(spline_basis)

Q2 <- dbs(dx$k, knots = knots, degree = spline_degree, derivs = 2, intercept = TRUE)
Q2_spline <- as.matrix(Q2)

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,
  Q = ncol(B),
  B = B,
  Q2_spline = Q2_spline,
  lambda = 0.15,
  index = dx$index,
  T = nrow(ties),
  J = max(ties$grp),
  t_grp = ties$grp,
  t_index = ties$index,
  t_adj = ties$adj
)

Acum probăm din posterior – puteți vedea că este nevoie de destul de mult timp pentru a rula, cel puțin pe 2020 MacBook Pro M1 cu 8 GB RAM:

fit_mcmc <- 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 4 finished in 1847.8 seconds.
## Chain 1 finished in 2202.8 seconds.
## Chain 3 finished in 2311.8 seconds.
## Chain 2 finished in 2414.9 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 2194.3 seconds.
## Total execution time: 2415.3 seconds.
fit_mcmc$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.815  0.815 0.0298 0.0298 0.767 0.865  1.00    3513.    5077.
## 2 sigma_b  0.543  0.535 0.0775 0.0739 0.432 0.683  1.00    3146.    5110.

Estimarea unui model de efecte aleatorii „frecventiste”

La urma urmei, se dovedește că puteți potrivi doar un model de fragilitate cu efecte aleatorii pentru site și un spline pentru perioada de timp (k ) folosind coxmme pachet. În mod evident, acest lucru este mult mai simplu decât tot ce am prezentat aici.

frailty_model <- coxme(Surv(Y, event) ~ A + ns(k, df = 3) + (1 | site), data = dd)
summary(frailty_model)
## Mixed effects coxme model
##  Formula: Surv(Y, event) ~ A + ns(k, df = 3) + (1 | site) 
##     Data: dd 
## 
##   events, n = 14989, 18750
## 
## Random effects:
##   group  variable        sd  variance
## 1  site Intercept 0.5306841 0.2816256
##                   Chisq    df p   AIC   BIC
## Integrated loglik 18038  5.00 0 18028 17990
##  Penalized loglik 18185 27.85 0 18129 17917
## 
## Fixed effects:
##                    coef exp(coef) se(coef)      z      p
## A               0.80966   2.24714  0.02959  27.36 <2e-16
## ns(k, df = 3)1 -2.71392   0.06628  0.04428 -61.29 <2e-16
## ns(k, df = 3)2  1.04004   2.82933  0.07851  13.25 <2e-16
## ns(k, df = 3)3  4.48430  88.61492  0.04729  94.83 <2e-16

Cu toate acestea, avantajul modelului Bayesian este flexibilitatea acestuia. De exemplu, dacă doriți să includeți curbele spline specifice site-ului-efecte de timp specifice site-ului-puteți extinde abordarea Bayesiană pentru a face acest lucru. Modelul actual Bayesian implementează un spline de timp la nivel de studiu, dar încorporarea splinelor specifice sitului ar fi o extensie naturală. Inițial am sperat să implementez spline specifice site-ului folosind mgcv Pachet, dar modelele nu au converg. Sunt destul de încrezător că o extensie Bayesiană ar fi, probabil, ar fi nevoie de resurse substanțiale de calcul. Dacă cineva vrea să încerc asta, cu siguranță aș putea, dar deocamdată, mă voi opri aici.