|
.impute.srp_model <- function(model, data, nsim, parameter_sample, seed = NULL, ...) { |
|
if (!is.null(seed)) { |
|
set.seed(seed) |
|
} |
|
# TODO: convert data to matrix and process in c++ |
|
stopifnot(isa(parameter_sample, "stanfit")) |
|
# extract subject and group id levels for conversion to and back from integer |
|
subject_id_levels <- unique(as.character(data$subject_id)) |
|
group_id_levels <- attr(model, "group_id") # important to maintain ordering |
|
# extract parameter matrices |
|
p <- rstan::extract(parameter_sample, "p")[[1]] |
|
scale <- rstan::extract(parameter_sample, "scale")[[1]] |
|
shape <- rstan::extract(parameter_sample, "shape")[[1]] |
|
|
|
data <- data %>% |
|
arrange(.data$t_sot, .data$subject_id, (.data$t_min + .data$t_max)/2) %>% |
|
mutate( |
|
subject_id = as.integer(factor(as.character(.data$subject_id), levels = subject_id_levels)), |
|
group_id = as.integer(factor(.data$group_id, levels = group_id_levels)) |
|
) |
|
|
|
res <- tribble( |
|
~subject_id, ~group_id, ~from, ~to, ~t_min, ~t_max, ~t_sot, ~iter |
|
) |
|
|
|
visit_spacing <- attr(model, "visit_spacing") |
|
idx <- sample(1:dim(p)[1], size = nsim, replace = TRUE) |
|
|
|
for (i in 1:nrow(data)) { |
|
if (!is.na(data$to[i])) { # observed transition, nothing to sample |
|
res <- bind_rows(res, tidyr::expand_grid(data[i, ], iter = 1:nsim)) |
|
} else { # a censored transition |
|
for (j in 1:nsim) { |
|
g <- data$group_id[i] |
|
k <- idx[j] |
|
sshape <- shape[k, g, ] |
|
sscale <- scale[k, g, ] |
|
if (data$from[i] == "stable") { |
|
# first sample response/progression |
|
pr_response_raw <- p[k, g] # use survival information! |
|
pr_survival_response <- 1 - stats::pweibull(data$t_min[i] - data$t_sot[i], sshape[1], sscale[1]) |
|
pr_survival_progression <- 1 - stats::pweibull(data$t_min[i] - data$t_sot[i], sshape[2], sscale[2]) |
|
pr_response <- pr_response_raw * pr_survival_response / ( |
|
pr_response_raw * pr_survival_response + |
|
(1 - pr_response_raw) * pr_survival_progression |
|
) |
|
response <- stats::rbinom(1, 1, pr_response) |
|
if (response) { |
|
# sample exact response time |
|
t_response <- rtruncweibull( |
|
sshape[1], scale = sscale[1], data$t_min[i], Inf # t_min since time of SoT is known |
|
) |
|
# apply visit scheme |
|
n_visits_response <- t_response %/% visit_spacing[g] |
|
tmin_response <- data$t_min[i] + visit_spacing[g] * n_visits_response |
|
tmax_response <- data$t_min[i] + visit_spacing[g] * (n_visits_response + 1) |
|
# sample subsequent progression, |
|
dt_progression <- stats::rweibull(1, sshape[3], sscale[3]) |
|
# apply visit scheme |
|
n_visits_progression <- (dt_progression + t_response) %/% visit_spacing[g] |
|
tmin_progression <- data$t_min[i] + visit_spacing[g] * n_visits_progression |
|
tmax_progression <- data$t_min[i] + visit_spacing[g] * (n_visits_progression + 1) |
|
res <- bind_rows( |
|
res, tribble( |
|
~subject_id, ~group_id, ~from, ~to, ~t_min, ~t_max, ~t_sot, ~iter, |
|
data$subject_id[i], g, "stable", "response", tmin_response, tmax_response, data$t_sot[i], j, |
|
data$subject_id[i], g, "response", "progression", tmin_progression, tmax_progression, data$t_sot[i], j |
|
) |
|
) |
|
} else { # sample progression directly |
|
dt_progression <- rtruncweibull( |
|
sshape[2], scale = sscale[2], data$t_min[i], Inf # t_min since time of SoT is known |
|
) |
|
# apply visit scheme |
|
n_visits_progression <- dt_progression %/% visit_spacing[g] |
|
tmin_progression <- data$t_min[i] + visit_spacing[g] * n_visits_progression |
|
tmax_progression <- data$t_min[i] + visit_spacing[g] * (n_visits_progression + 1) |
|
res <- bind_rows( |
|
res, tribble( |
|
~subject_id, ~group_id, ~from, ~to, ~t_min, ~t_max, ~t_sot, ~iter, |
|
data$subject_id[i], g, "stable", "progression", tmin_progression, tmax_progression, data$t_sot[i], j |
|
) |
|
) |
|
} # end stable -> progression |
|
} # end from == 1 |
|
if (data$from[i] == "response") { |
|
if (data$from[i - 1] != "stable" || data$subject_id[i - 1] != data$subject_id[i]) { |
|
stop() |
|
} |
|
# sample exact response time |
|
t_response <- rtruncweibull( |
|
sshape[1], scale = sscale[1], data$t_min[i - 1], data$t_max[i - 1] |
|
) |
|
# sample progression time |
|
dt_progression <- rtruncweibull( |
|
sshape[3], scale = sscale[3], data$t_min[i] - t_response, Inf |
|
) |
|
# apply visit scheme |
|
n_visits_progression <- dt_progression %/% visit_spacing[g] |
|
tmin_progression <- data$t_min[i] + visit_spacing[g] * n_visits_progression |
|
tmax_progression <- data$t_min[i] + visit_spacing[g] * (n_visits_progression + 1) |
|
res <- bind_rows( |
|
res, tribble( |
|
~subject_id, ~group_id, ~from, ~to, ~t_min, ~t_max, ~t_sot, ~iter, |
|
data$subject_id[i], g, "response", "progression", tmin_progression, tmax_progression, data$t_sot[i], j |
|
) |
|
) |
|
} # end from == 2 |
|
} # end iterate of j |
|
} # end if/else |
|
} # end iteration over i |
|
# convert subject and group id back |
|
res <- res %>% |
|
mutate( |
|
subject_id = as.character( |
|
factor(.data$subject_id, levels = seq_along(subject_id_levels), labels = subject_id_levels) |
|
), |
|
group_id = as.character( |
|
factor(.data$group_id, levels = seq_along(group_id_levels), labels = group_id_levels) |
|
) |
|
) |
|
return(res) |
|
} |