Nothing
# Helper functions
#' @import dplyr
#' @import rjags
#' @import mvtnorm
#' @import MatchIt
#' @import ggplot2
#' @import survival
#' @import futile.logger
#' @import doParallel
#' @import foreach
#' @importFrom data.table data.table :=
#' @importFrom methods new slot
#' @importFrom stats as.formula binomial glm qnorm rbinom rexp rnorm rweibull time update var
#' @importFrom grDevices png dev.off
`%notin%` <- Negate(`%in%`)
# The piecewise exponential distribution useful for conditions where failure rates change or simulations with a delayed or changing treatment effect
rpwexp <- function(n, rate=1, intervals=NULL, cumulative=FALSE){
if(is.null(intervals)){
if (cumulative){return(cumsum(rexp(n,rate[1])))}else
return(rexp(n,rate[1]))}
k <- length(rate)
if (k==1){
if(cumulative){return(cumsum(rexp(n,rate)))}else
return(rexp(n,rate))
}
if (length(intervals) < k-1) stop("length(intervals) must be at least length(rate) - 1")
tx <- 0
j <- 1
times <- array(0,n)
timex <- cumsum(intervals)
indx <- array(TRUE,n)
for(i in 1:k){
nindx <- sum(indx)
if (nindx==0) break
increment <- rexp(nindx,rate[i])
if (cumulative) times[indx] <- tx + cumsum(increment)
else times[indx] <- tx + increment
if (i<k){
tx <- timex[i]
indx <- (times > timex[i])
}
}
return(times)
}
# Specify sample size
s_trt = function(ssC, ssE, ssExt){
ext = cbind("ext" = rep(1, ssExt), "trt" = rep(0, ssExt))
# int = cbind("ext" = rep(0, ssC + ssE), "trt" = sample(c(rep(1, ssE), c(rep(0, ssC)))))
int = cbind("ext" = rep(0, ssC + ssE), "trt" = c(rep(1, ssE), c(rep(0, ssC))))
rbind(ext, int)
}
# Simulate covariates
s_cov = function(ext, dt, n_cat, n_cont, mu, var, cov, prob){
n = sum(dt[,'ext'] == ext)
flog.debug(paste("Ext =", ext, "trial has sample size =", n, "\n"))
cov_st = sum(grepl("cov", colnames(dt))) + 1
covariate <- NULL
n_cov = sum(n_cont, n_cat)
flog.debug(paste("[s_cov] Total number of covariates to add this time =", n_cov, "\n"))
if (n_cat == 0 & n_cont == 1){
covariate <- as.matrix(rnorm(n, mean = mu, sd = sqrt(var)), ncol = 1)
flog.debug(paste("User chose to simulate one continuous variable. \n"))
} else if (n_cont == 0 & n_cat > 0 & sum(cov) == 0){ # if user want to simulate independent binary variables only
covariate = sapply(1:n_cat, function(i){
flog.debug(paste("User chose to simulate independent binary variable only. \n"))
rbinom(n, 1, prob[i])
})
} else {
sig <- diag(var)
sig[lower.tri(sig)] <- cov
sig[upper.tri(sig)] = t(sig)[upper.tri(sig)] # create a symmetrical matrix
covariate <- rmvnorm(n, mean = mu, sigma = sig)
# flog.debug(paste("User chose to simulate continuous covariate. \n"))
if (n_cat > 0) {
covariate[,1:n_cat] = sapply(1:n_cat, function(i){
cut_val = qnorm(prob[i], mean = mu[i], sd = sqrt(var[i]))
flog.debug(cat("User chose to simulate", n_cov, "correlated binary covariate(s) with cut_val =", cut_val, "\n"))
(covariate[, i] < cut_val)
})
}
}
colnames(covariate) = paste0("cov",cov_st:(cov_st + n_cov - 1))
covariate
}
# raise covariates to power
# dt <- sample_cov
# change = list(c("cov1", "^", "2"),c("cov1", "*", "cov3"), c("cov1", "*", "cov3"), c("cov1", "*", "cov4"))
# keep = "cov1^2"
c_cov = function(dt, change, keep){
dt2 <- dt
new_cov_name = NULL
for (i in 1:length(change)){
if (change[[i]][2] == "^") {
new_cov <- dt[, change[[i]][1]] ^ as.numeric(change[[i]][3])
flog.debug("^")
} else if (change[[i]][2] == "*") {
new_cov <- dt[, change[[i]][1]] * dt[, change[[i]][3]]
flog.debug("*")
} else if (change[[i]][2] == "+") {
new_cov <- dt[, change[[i]][1]] + dt[, change[[i]][3]]
flog.debug("+")
}
dt2 <- cbind(new_cov, dt2)
colnames(dt2)[1]<- paste(change[[i]],collapse="")
new_cov_name = c(new_cov_name, paste(change[[i]],collapse=""))
}
flog.debug(cat("[c_cov] Original covariates to keep include", keep, "new covariates to keep include", new_cov_name, "\n"))
flog.debug(cat("[c_cov] Dataset has variables", colnames(dt2), "\n"))
dt3 <- dt2[, c(keep, new_cov_name, "trt"), drop=FALSE]
flog.debug(cat("[c_cov] Simulate failure time based on", colnames(dt3), "\n"))
dt3
}
# match individuals from the external dataset to the internal tx group
m_cov <- function(dt, match.formula){
dt2 <- data.frame(dt)
dt2$int = 1 - dt2$ext
m.out <- matchit(as.formula(match.formula), method = "nearest", ratio = 1, caliper = 0.2, data = dt2)
#matchit(ext ~ cov1 + cov2, method = "nearest", ratio = 1, caliper = 0.2, data = sample_cov.df) # match on all covariates
dt.match <- match.data(m.out)
ps_message(nrow(dt2[dt2$ext == 1,]) - nrow(dt.match[dt.match$ext == 1,]), " patients from the external arm are excluded. ",
nrow(dt2[dt2$int == 1,]) - nrow(dt.match[dt.match$int == 1,]), " patients from the internal trial are excluded.")
dt.match_ext <- dt.match[dt.match$ext == 1, colnames(dt.match) %notin% c("distance", "weights", "subclass")]
dt_int = dt2[dt2$ext == 0,]
match.all <- as.matrix(rbind(dt.match_ext, dt_int) %>% select(-"int"))
ps_message("After matching, we have ", nrow(match.all), " patients in total.")
return(match.all)
}
# Simulate time-to-events following weibull distribution
s_wb_t = function(dt, lambda, HR, beta, shape){
hr = unique(HR)
flog.debug(paste("[s_wb_t] HR =", hr, "\n"))
n = nrow(dt)
n_cov = sum(grepl("cov", colnames(dt))) #check if the data include any covariate
cov = dt[, colnames(dt) %notin% c("trt", "ext"), drop = FALSE]
flog.debug(cat("[s_wb_t] Simulate time using weibull distribution with sample size =", n, "number of covariates =", n_cov, ":", colnames(cov), "\n"))
flog.debug(cat("[s_wb_t] Dimension of covariates", dim(cov), "length of beta", length(beta), "\n"))
if (n_cov == 0) {
log_k = - (log(lambda) + log(hr) * dt[, "trt"])/shape
} else {
log_k = - (log(lambda) + log(hr) * dt[, "trt"] + cov %*% beta)/shape #log(k), k is the scale
}
# print(paste("length of log of scale for the weibull distribution =", length(log_k)))
sapply(log_k, function(i){rweibull(1, shape, exp(i))})
}
# Simulate time-to-events following piecewise-exponential distribution
s_exp_t = function(dt, lambda, HR, beta, t_itv){
hr = unique(HR)
flog.debug(paste("[s_wb_t] HR =", hr, "\n"))
flog.debug(cat("[s_exp_t] Simulate time using piece-wise exponential distribution with baseline hazard rate of control arm =", lambda, "\n"))
n = nrow(dt)
n_cov = sum(grepl("cov", colnames(dt))) #check if the data include any covariate
if (n_cov == 0) lambda_x = exp(log(hr) * dt[, "trt"]) %*% t(lambda)
else lambda_x = exp(log(hr) * dt[, "trt"] + dt[, colnames(dt) %notin% c("driftHR", "HR", "trt", "ext")] %*% beta) %*% t(lambda)
flog.debug(cat("Baseline hazard rate has dimension of", dim(lambda_x), "\n"))
apply(lambda_x, 1, function(i) rpwexp(1, rate = i, intervals = t_itv, cumulative=FALSE))
}
# Simulate enrollment time following piecewise-exponential distribution
s_exp_e = function(n, gamma, e_itv){
aR <- e_itv
if (n-sum(gamma[1:length(gamma)-1]*e_itv[1:length(gamma)-1])<=0) {
stop ("[s_exp_e] User requested to keep the enrollment rate fixed but the total enrollment is already greater than the specified sample size prior to the last interval of R.")
} else{
aR[length(gamma)] <- (n - sum(gamma[1:length(gamma)-1]*e_itv[1:length(gamma)-1]))/gamma[length(gamma)]
flog.debug(cat("[s_exp_e] Updated enrollment intervals", aR, "\n"))
# Using the cumulative=TRUE option, enrollment times that piecewise constant over time can be generated.
sample(rpwexp(n, rate=gamma, intervals = e_itv[1:length(gamma)-1], cumulative=TRUE))
}
}
# Simulate censored survival times
g_one_t <- function(ext, dt,
gamma, e_itv, # enroll
event, lambda, HR, beta, shape, t_itv, # event
etaC, etaE, d_itv, # drop out
CCOD, CCOD_t, # CCOD
change, keep){
# set new variables to NULL to avoid CRAN warnings
time0 <- time1 <- cnsr <- cnsr1 <- ct <- enterT <- ltfuT <- NULL
# print("lambda");print(lambda)
d <- dt[dt[, "ext"] == ext,] # matrix
flog.debug(cat("[g_one_t] We are simulating covariates for ext =", ext, "All covariates include", colnames(d), "\n"))
x <- data.table(d)
n = nrow(x)
ssC = sum(x$trt == 0); ssE = sum(x$trt == 1)
flog.debug(paste("[g_one_t] Simulate time with sample size =", n, "number of control =", ssC, "number of treatment =", ssE, "\n"))
x[, enterT := s_exp_e(n, gamma, e_itv)] #duration
# update matrix
d_change = d[, colnames(d) %notin% c("driftHR", "HR", "ext"), drop = FALSE]
if (!is.null(change)) d_change = c_cov(d[, colnames(d) %notin% c("driftHR", "HR", "ext"), drop = FALSE], change, keep) #update the matrix
flog.debug(cat("[g_one_t] Covariates for simulating failure times include", colnames(d_change), "\n"))
flog.debug(paste("[g_one_t] HR =", HR, "\n"))
# survival time
if (event == "weibull") x[, time0 := s_wb_t(d_change, lambda, HR, beta, shape)]
else x[, time0 := s_exp_t(d_change, lambda, HR, beta, t_itv)]
# adding LTFU
x[x$"trt" == 0, ltfuT := rpwexp(ssC, etaC, intervals = d_itv, cumulative=FALSE)]
x[x$"trt" == 1, ltfuT := rpwexp(ssE, etaE, intervals = d_itv, cumulative=FALSE)]
x[, time1:= ifelse(time0 > ltfuT, ltfuT, time0)]
x[, cnsr1 := ifelse(time0 > ltfuT, 1, 0)]
x[, ct := enterT + time1] # ct: calendar time a subject had evt/censoring
# add analysis start time
if (CCOD == "fixed-first") { #pre-fixed, calculate from enrollment of first patients
x[, time := ifelse(ct > CCOD_t, CCOD_t - enterT, x$time1)]
x[, cnsr := ifelse(ct > CCOD_t, 1, x$cnsr1)]
ps_message(paste(sum(x$time < 0), "patients enrolled after", CCOD_t,
"days after the study started were excluded from the analysis."))
x = x[x$time >= 0, ]
} else if (CCOD == "fixed-last") {
st = max(x$enterT) + CCOD_t
x[, time := ifelse(ct > st, st - enterT, x$time1)]
x[, cnsr := ifelse(ct > st, 1, x$cnsr1)]
x[, cnsr := ifelse(ct > CCOD_t, 1, x$cnsr1)]
ps_message(paste(sum(x$time < 0), "patients enrolled after", CCOD_t,
"days after the enrollment of last patient were excluded from the analysis."))
}
else { # analysis starts when st events have been observed
st = sort(x$ct[x$cnsr1 == 0])[pmin(sum(x$cnsr1 == 0), CCOD_t)] # compete leaving with event
if (identical(st, numeric(0))) stop ("Please decrease the number of events required to start the analysis (CCOD_t)")
x[, time := ifelse(ct > st, st - enterT, x$time1)]
x[, cnsr := ifelse(ct > st, 1, x$cnsr1)]
ps_message("User chooses to let number of events drive the CCOD", sum(x$time < 0),
"patients were enrolled in the trial after ", st,
"days and were excluded from the analysis")
x = x[x$time >= 0, ]
}
x_fin <- x[, -c("enterT", "time0", "ltfuT", "time1", "cnsr1", "ct")]
as.matrix(sapply(x_fin, as.numeric))
}
# Calculate propensity score
# Function r_post below calls this
c_ps <- function(dt, cov_name){
flog.debug(cat("[c_ps] cov_name =", cov_name, "\n"))
dat = as.data.frame(dt,stringsAsFactors=FALSE)
dat$int = 1 - dat$ext
var_col = c(cov_name, "int")
dat2 = dat[, var_col]
mod <- glm(int ~ ., family = binomial(), data = dat2)
cbind(dt, "wgt" = mod$fitted.values)
}
# Generate posterior sample
r_post <- function(dt,
prior, pred,
r0, alpha, sigma,
n.chains, n.adapt, n.burn, n.iter, seed){
flog.debug(cat("[r_post] dt =", colnames(dt), "\n"))
cov_name = colnames(dt)[colnames(dt) %notin% c("driftHR", "HR", "trt", "ext", "time", "cnsr", "wgt")]
flog.debug(cat("[r_post] cov_name =", cov_name, "\n"))
# pred = c("ps", "cov1")
# dt is a matrix
if (sum(grepl("none", pred)) > 0) {
x = dt[,c("trt"),drop=FALSE]
n_cov = 0
flog.debug("No covariates are used as predictors in the weibull distribution. Any other input is disregarded.\n")
} else if (sum(grepl("all", pred)) > 0) {
x = dt[,c("trt", cov_name)]
n_cov = length(cov_name)
flog.debug(cat("All covariates", cov_name, "are used as predictors in the weibull distribution. Any other input is disregarded. \n"))
} else if (sum(grepl("ps", pred)) == 0) {
x = dt[,c("trt", pred)]
n_cov = length(pred)
flog.debug(cat("Selected covariate(s)", pred, "are used as predictors in the weibull distribution. Any other input is disregarded. \n"))
} else if (sum(grepl("ps", pred)) > 0 & ("wgt" %in% colnames(dt))) { # column wgt already existed
x = dt[,c("trt", "wgt")]
n_cov = 1
flog.debug("Weight already provided in the data. It is used as predictor directly. \n")
} else if (sum(grepl("ps", pred)) > 0 & length(pred) == 1){ # column wgt does not exist
dt_ps <- c_ps(dt, cov_name)
x = dt_ps[,c("trt", "wgt")]
n_cov = 1
flog.debug(cat("Propensity score is calculated based on all coariates", cov_name, "are used as predictors in the weibull distribution.\n"))
} else if (sum(grepl("ps", pred)) > 0 & length(pred) > 1){
dt_ps <- c_ps(dt, pred[pred != "ps"])
x = dt_ps[,c("trt", "wgt")]
n_cov = 1
flog.debug(cat("Propensity score calculated based selected covariate(s)", pred[pred != "ps"],
"ise used as predictors in the weibull distribution.\n"))
}
data_list = list(N = nrow(dt), timev = dt[,'time'], event = 1-dt[,'cnsr'], ext = dt[,'ext'] +1, n_cov = n_cov, x = x)
inits_list = list(.RNG.name="base::Super-Duper", .RNG.seed = seed, r0 = r0, alpha = alpha, beta = rep(0, n_cov + 1)) #.RNG.seed = 1,
out_list = c("alpha", "beta", "r0", "HR_trt_cc", "HR_cc_hc")
if (prior == "full_ext") {
prior = "no_ext"
data_list = list(N = nrow(dt), timev = dt[,'time'], event = 1-dt[,'cnsr'], n_cov = n_cov, x = x)
} else if (prior == "no_ext") { # no borrow
dt2 = dt[dt[, "ext"] == 0,] #remove external trial data
if (sum(grepl("none", pred)) > 0) {
x = dt2[,c("trt"),drop=FALSE]
} else if (sum(grepl("all", pred)) > 0) {
x = dt2[,c("trt", cov_name)]
} else {
x = dt2[,c("trt", pred)]
}
data_list = list(N = nrow(dt2), timev = dt2[,'time'], event = 1-dt2[,'cnsr'], n_cov = n_cov, x = x)
} else if (prior == "cauchy" | prior == "unif") {
inits_list = list(.RNG.name = "base::Super-Duper", .RNG.seed = seed, r0 = r0, alpha = alpha, beta = rep(0, n_cov + 1), sigma = sigma) #.RNG.seed = 1,
out_list = c("alpha", "beta", "r0", "prec", "HR_trt_cc", "HR_cc_hc")
} else if (prior == "gamma") {
out_list = c("alpha", "beta", "r0", "tau", "HR_trt_cc", "HR_cc_hc")
}
prior_txt = system.file("model", paste0(prior, ".txt"), package = "psborrow")
# prior_txt = paste0("inst/model/", prior, ".txt")
flog.debug(cat("[r_post] N =", data_list[['N']], "event =", sum(data_list[['event']] == 1),
"ext =", ifelse(is.null(data_list[['ext']]), NA, sum(data_list[['ext']] == 2)),
"n_cov =", data_list[['n_cov']], "x =", colnames(data_list[['x']]), "\n"
))
flog.debug(cat("[r_post] r0 =", inits_list[['r0']],
"alpha =", inits_list[['alpha']], "beta =", inits_list[['beta']],
"sigma =", ifelse(is.null(inits_list[['sigma']]), NA, inits_list[['sigma']]),
".RNG.seed = ", seed, "\n"
))
jags.out <- jags.model(
prior_txt,
data = data_list,
inits = inits_list,
n.chains = n.chains,
n.adapt = n.adapt,
quiet = TRUE
)
update(jags.out, n.burn = n.burn)
# If user has asked for psborrow to be quiet then we set progress bar to be
# "none" instead of using the default JAGS argument specified by "jags.pb"
progress_bar <- ifelse(
options("psborrow.quiet")[[1]],
"none",
options("jags.pb")[[1]]
)
jags.sample <- coda.samples(
jags.out,
out_list,
n.iter = n.iter,
progress.bar = progress_bar
)
jags.sample
}
#' Conditional Message
#'
#' Simple wrapper function around [message()] that will supress
#' printing messages if the option `psborrow.quiet` is set to `TRUE`
#' i.e.
#' ```
#' options("psborrow.quiet" = TRUE)
#' ```
#'
#' @param ... Values passed onto [message()]
ps_message <- function(...){
if (!getOption("psborrow.quiet")) {
message(...)
}
}
#' Check if user is in psborrow development environment
#'
#' Simple function which leverages the DESCRIPTION file
#' to check if the user is in a development environment
#' for psborrow.
#'
#' @return TRUE/FALSE flag (TRUE = in development environment)
is_psborrow_dev <- function() {
desc <- file.path("DESCRIPTION")
if (file.exists(desc)) {
content <- read.dcf(desc)
pkg <- content[1, ][["Package"]]
return(pkg == "psborrow")
}
return(FALSE)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.