Closed mattsecrest closed 1 year ago
Using the Poisson trick it's quite easy actually. We can also process the data in Stan for the cuts.
## Stan
pem <- cmdstanr::cmdstan_model(stan_file = write_stan_file(
"
data {
int<lower=0> N;
vector[N] trt;
vector[N] time;
vector[N] cens;
int<lower = 1> M;
vector[M] starts;
vector[M] durations;
}
transformed data {
matrix[N,M] T;
matrix[N,M] D;
for(j in 1:M) {
T[,j] = fmin(fdim(time, starts[j]), durations[j]);
}
for(j in 1:M) {
for(i in 1:N) {
D[i,j] = (starts[j] <= time[i] && time[i] < starts[j] + durations[j]) * (1 - cens[i]);
}
}
}
parameters {
real beta_trt;
vector[M] alpha;
}
transformed parameters {
real HR_trt = exp(beta_trt);
}
model {
vector[N] lp;
beta_trt ~ normal(0, 1e+05);
alpha ~ normal(0, 1e+05);
for (j in 1:M){
lp = alpha[j] + trt * beta_trt;
target += (lp' * D[,j]) - (exp(lp)' * T[,j]);
}
}
"))
pem_samples <- pem$sample(
data = list(
N = nrow(example_matrix),
trt = example_matrix[,"trt"],
cens = example_matrix[,"cnsr"],
time = example_matrix[,"time"],
M = 4,
starts = c(0, 1.83, 6, 23),
durations = diff(c(c(0, 1.83, 6, 23), max(example_matrix[,"time"])))
))
pem_results <- pem_samples$summary(
variables = c("beta_trt", "alpha[1]", "alpha[2]", "alpha[3]", "alpha[4]"),
mean = mean,
median = median
)
pem_results$exp_mean <- exp(pem_results$mean)
pem_results$exp_median <- exp(pem_results$median)
pem_results
## Comparison
em <- as.data.frame(example_matrix)
em$event <- 1 - em$cnsr
res_eha_pchreg <- eha::pchreg(
Surv(time, event) ~ trt,
cuts = c(0, 1.83, 6, 23, 50),
data = em
)
summary(res_eha_pchreg)
res_eha_pchreg$hazards
Closing in favor of newer issue https://github.com/Genentech/psborrow2/issues/223
as.matrix(survSplit(data = data.frame(example_matrix), cut = c(5, 10), end = "time", event = "cnsr"))