Un model de pericole proporționale bayesiene cu un spline penalizat

Simularea datelor cu un raport de pericol care variază în timp

Iată R Pachete utilizate în postare:

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

Setul de date simulează un studiu controlat randomizat în care pacienții sunt repartizați fie la grupul de tratament ( (A = 1 )) sau grup de control ( (A = 0 )) în a (1: 1 ) raport. Pacienții se înscriu peste nouă sferturi, cu trimestrul de înscriere notat de (M ), (M in {0, dots, 8 } ). Rezultatul timpului de eveniment, (Y )depinde atât de atribuirea tratamentului, cât și de trimestrul de înscriere. Pentru a introduce neliniaritatea, definesc relația folosind o funcție cubică, cu parametrii adevărați specificați după cum urmează:

defI <- 
  defData(varname = "A", formula = "1;1", dist = "trtAssign") |>
  defData(varname = "M", formula = "0;8", dist = "uniformInt")

defS <-
  defSurv(
    varname = "eventTime",
    formula = "..int + ..beta * A + ..alpha_1 * M + ..alpha_2 * M^2 + ..alpha_3 * M^3",
    shape = 0.30)  |>
  defSurv(varname = "censorTime", formula = -11.3, shape = 0.40)

# parameters

int <- -11.6      
beta <-  0.70
alpha_1 <-  0.10   
alpha_2 <-  0.40    
alpha_3 <- -0.05

Am generat un singur set de date de (640 ) participanți la studiu, (320 ) în fiecare braț. Graficul de mai jos arată curbele Kaplan-Meier cu braț pentru fiecare perioadă de înscriere.

set.seed(7368) # 7362

dd <- genData(640, defI)
dd <- genSurv(dd, defS, timeName = "Y", censorName = "censorTime",
  eventName = "event", typeName = "eventType", keepEvents = TRUE)

Modelul Bayesian

Acest model de pericole proporționale Bayesian se bazează direct pe abordarea din postarea mea anterioară. De la efectul (M ) pe (Y ) Urmează un model neliniar, modelez această relație folosind un spline pentru a ține cont de variația temporală a pericolului. Probabilitatea parțială este o funcție a efectului de tratament și a coeficienților funcției de bază spline, date de:

( L(beta,mathbf{gamma}) = prod_{i=1}^{N} left( frac{exp left(beta A_i + sum_{m=1} ^ M gamma_m X_{m_i} right)} {sum_{j in R(t_i)} expleft(beta A_j + sum_ {m = 1} ^ m gamma_m x_ {m_j} dreapta)} dreapta) ^ { delta_i} )
unde:

• (M ): numărul de funcții de bază spline
• (N ): numărul de observații (cenzurat sau nu)
• (A_i ): indicator binar pentru tratament
• (X_ {m_i} ): valoarea (m^{ text {th}} ) funcție de bază spline pentru (i^{ text {th}} ) observare
• ( delta_i ): indicator de eveniment ( ( delta_i = 1 ) Dacă a avut loc evenimentul, ( delta_i = 0 ) dacă este cenzurat)
• ( beta ): coeficient de tratament
• ( gamma_m ): coeficient spline pentru (m^ text {th} ) Funcția de bază spline
• (R (t_i) ): risc stabilit la timp (t_i ) (inclusiv doar persoane cenzurate după (t_i ))

Componenta spline a modelului este adaptată de la un model pe care l -am descris anul trecut. În această formulare, timpul de petrecere este modelat ca funcție a vectorului ( mathbf {x_i} ) mai degrabă decât perioada în sine. Numărul de funcții de bază este determinat de numărul de noduri, cu fiecare segment al curbei estimat folosind funcții de bază B-spline. Pentru a reduce la minimum suprapunerea, includem un termen de penalizare pe baza celui de-al doilea derivat al funcțiilor de bază B-spline. Puterea acestei penalizări este controlată de un parametru de reglare, ( lambda )care este furnizat modelului.

Codul Stan, furnizat integral aici, a fost explicat în postările anterioare. Diferența principală față de postarea anterioară este adăugarea datelor și parametrilor legați de spline, precum și termenul de penalizare în model.:

stan_code <-
"
functions {

  // Binary search optimized to return the last index with the target value

  int binary_search(vector v, real tar_val) {
    int low = 1;
    int high = num_elements(v);
    int result = -1;

    while (low <= high) {
      int mid = (low + high) %/% 2;
      if (v(mid) == tar_val) {
        result = mid; // Store the index
        high = mid - 1; // Look for earlier occurrences
      } else if (v(mid) < tar_val) {
        low = mid + 1;
      } else {
        high = mid - 1;
      }
    }
    return result;
  }
}

data {

  int K;          // Number of covariates
  int N_o;        // Number of uncensored observations
  vector(N_o) t_o;         // Event times (sorted in decreasing order)

  int N;          // Number of total observations
  vector(N) t;             // Individual times (sorted in decreasing order)
  matrix(N, K) x;          // Covariates for all observations

  // Spline-related data
  
  int Q;          // Number of basis functions
  matrix(N, Q) B;          // Spline basis matrix
  matrix(N, Q) D2_spline;  // 2nd derivative for penalization
  real lambda;             // penalization term
}

parameters {
  vector(K) beta;          // Fixed effects for covariates
  vector(Q) gamma;         // Spline coefficients
}

model {
  
  // Prior
  
  beta ~ normal(0, 4);
  
  // Spline coefficients prior
  
  gamma ~ normal(0, 4);
  
  // Penalization term for spline second derivative
  
  target += -lambda * sum(square(D2_spline * gamma));
  
  // Calculate theta for each observation to be used in likelihood
  
  vector(N) theta;
  vector(N) log_sum_exp_theta;
  
  for (i in 1:N) {
    theta(i) = dot_product(x(i), beta) + dot_product(B(i), gamma);  
  }
  
  // Compute 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
    log_sum_exp_theta(i) = log_sum_exp(theta(i), log_sum_exp_theta(i + 1));
  }

  // Likelihood for uncensored observations
  
  for (n_o in 1:N_o) {
    int start_risk = binary_search(t, t_o(n_o)); // Use binary search
    
    real log_denom = log_sum_exp_theta(start_risk);
    target += theta(start_risk) - log_denom;
  }
}
"

Pentru a estima modelul, trebuie să pregătim datele pentru a trece Stancompilează Stan cod, apoi probe din model folosind cmdstanr:

dx <- copy(dd)
setorder(dx, Y)

dx.obs <- dx(event == 1)
N_obs <- dx.obs(, .N)
t_obs <- dx.obs(, Y)

N_all <- dx(, .N)
t_all <- dx(, Y)
x_all <- data.frame(dx(, .(A)))

# 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$M, probs = probs)
spline_basis <- bs(dx$M, knots = knots, degree = spline_degree, intercept = TRUE)
B <- as.matrix(spline_basis)

D2 <- dbs(dx$M, knots = knots, degree = spline_degree, derivs = 2, intercept = TRUE)
D2_spline <- as.matrix(D2)

K <- ncol(x_all)             # num covariates - in this case just A

stan_data <- list(
  K = K,
  N_o = N_obs,
  t_o = t_obs,
  N = N_all,
  t = t_all,
  x = x_all,
  Q = ncol(B),
  B = B,
  D2_spline = D2_spline,
  lambda = 0.10
)

# compiling code

stan_model <- cmdstan_model(write_stan_file(stan_code))

# sampling from model

fit <- stan_model$sample(
  data = stan_data,
  iter_warmup = 1000,
  iter_sampling = 4000,
  chains = 4,
  parallel_chains = 4,
  max_treedepth = 15,
  refresh = 0
)
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 64.1 seconds.
## Chain 3 finished in 64.5 seconds.
## Chain 2 finished in 65.2 seconds.
## Chain 1 finished in 70.6 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 66.1 seconds.
## Total execution time: 70.8 seconds.

Media posterioară (și mediana) pentru ( beta )efectul de tratament, este destul de aproape de valoarea „adevărată” de 0,70:

fit$summary(variables = c("beta", "gamma"))
## # A tibble: 10 × 10
##    variable   mean median     sd    mad     q5   q95  rhat ess_bulk ess_tail
##                           
##  1 beta(1)   0.689  0.689 0.0844 0.0857  0.551 0.828  1.00    3664.    4002.
##  2 gamma(1) -1.75  -1.77  1.33   1.35   -3.91  0.468  1.00    1364.    1586.
##  3 gamma(2) -1.59  -1.60  1.33   1.35   -3.75  0.626  1.00    1360.    1551.
##  4 gamma(3) -1.22  -1.24  1.33   1.35   -3.39  0.978  1.00    1365.    1515.
##  5 gamma(4) -0.115 -0.127 1.33   1.35   -2.28  2.09   1.00    1361.    1576.
##  6 gamma(5)  1.97   1.95  1.34   1.35   -0.206 4.20   1.00    1366.    1581.
##  7 gamma(6)  2.63   2.61  1.33   1.34    0.452 4.84   1.00    1358.    1586.
##  8 gamma(7)  1.08   1.05  1.33   1.34   -1.08  3.28   1.00    1360.    1505.
##  9 gamma(8) -0.238 -0.260 1.33   1.34   -2.40  1.97   1.00    1355.    1543.
## 10 gamma(9) -0.914 -0.935 1.33   1.35   -3.07  1.30   1.00    1356.    1549.

Cifra de mai jos arată splina estimată și intervalul credibil de 95%. Linia verde reprezintă raportul de pericol de jurnal median posterior pentru fiecare perioadă (în raport cu perioada de mijloc, 4), banda umbrită indicând intervalul credibil corespunzător. Punctele violet reprezintă raporturile de pericol de jurnal implicate de procesul de generare a datelor. De exemplu, raportul de pericol de jurnal care compară perioada 1 la perioada 4 pentru ambele brațe este:

( begin{array}{c} (-11.6 + 0.70A +0.10times1 + 0.40 times 1^2 -0.05times1^3) – (-11.6 + 0.70A +0.10times4 + 0.40 times 4^2 -0.05times4^3) = \ (0.10 + 0.40 – 0.05) – (0.10 times 4 + 0.40 times 16 – 0,05 ori 64) = \ 0.45 – 3,60 = -3.15 end {array} )

Se pare că posterior median se aliniază destul de bine cu valorile utilizate în procesul de generare a datelor:

Pentru următoarea postare, voi prezenta un alt scenariu care include efecte aleatorii pentru un studiu randomizat în cluster (dar nu va include spline).