#' PS-Integrated Weighted Likelihood Estimation
#'
#' For all stratum. Variance estimated by bootstrap or Jack-Knife. Works for
#' one- or two-arm studies when there is one RWD
#'
#' @param data class DWITHPS data frame
#' @param ... parameters for \code{rweWL}
#' @param bs.n number of bootstraps
#' @param lambdas power parameter without standardization by ns0
#' @param m.var method to get variance: jackknife or bootstrap
#' @param seed random seed
#'
#' @export
#'
rwePsWl <- function(data, lambdas, v.outcome = "Y", m.var = c("jk", "bs"),
bs.n = 1000, seed = NULL, ...) {
stopifnot(v.outcome %in% colnames(data));
m.var <- match.arg(m.var);
if (!is.null(seed))
set.seed(seed);
## prepare data
data <- data[!is.na(data[["_strata_"]]), ]
S <- max(data[["_strata_"]]);
## find mwle
rst.theta <- NULL;
for (i in 1:S) {
cur.d1 <- data[data[["_strata_"]] == i & data[["_grp_"]] == 1, v.outcome];
cur.d0 <- data[data[["_strata_"]] == i & data[["_grp_"]] == 0, v.outcome];
ns1 <- length(cur.d1);
ns0 <- length(cur.d0);
if (0 == ns1) {
stop(paste("Stratum ", i,
" contains no subjects from group 1",
sep = ""))
}
cur.lambda <- lambdas[i];
cur.theta <- rweWL(cur.data = cur.d1,
ext.data = cur.d0,
lambda = cur.lambda, ...)
##bootstrap or jackknife
var.theta <- NULL;
if ("bs" == m.var) {
for (j in 1:bs.n) {
cur.d1.bs <- sample(cur.d1, replace = TRUE);
cur.bs <- rweWL(cur.data = cur.d1.bs, ext.data = cur.d0,
lambda = cur.lambda, ...);
var.theta <- c(var.theta, cur.bs);
}
var.mle <- var(var.theta);
} else if ("jk" == m.var) {
for (j in 1:ns1) {
cur.bs <- rweWL(cur.data = cur.d1[-j], ext.data = cur.d0,
lambda = cur.lambda, ...);
var.theta <- c(var.theta, cur.bs);
}
if (ns0 > 0) {
for (j in 1:ns0) {
cur.bs <- rweWL(cur.data = cur.d1, ext.data = cur.d0[-j],
lambda = cur.lambda, ...);
var.theta <- c(var.theta, cur.bs);
}
}
var.mle <- (ns1+ns0-1)/(ns1+ns0)*sum((var.theta - cur.theta)^2);
} else {
var.mle <- 0
}
rst.theta <- rbind(rst.theta, c(ns1, cur.theta, var.mle, ns0));
}
##mwle
ws <- rst.theta[,1] / sum(rst.theta[,1]);
rst.mwle <- sum(ws * rst.theta[, 2])
##rst.bs <- apply(rst.theta[,-(1:2), drop = FALSE], 2, function(x) sum(ws*x));
rst.var <- sum(ws^2 * rst.theta[, 3])
list(mwle = rst.mwle,
var = rst.var,
mwle.strata = rst.theta[, 2],
var.strata = rst.theta[, 3],
ns1 = rst.theta[, 1],
ns0 = rst.theta[, 4])
}
#' PS-Integrated Composite Likelihood Estimation
#'
#' Return all Jack-Knife values for all stratum. Works for
#' one--arm studies when there is one RWD
#'
#' @param data class DWITHPS data frame
#' @param ... parameters for \code{rweWL}
#' @param lambdas power parameter without standardization by ns0
#' @param seed random seed
#'
#' @export
#'
rwePsJkWl <- function(data, lambdas, v.outcome = "Y", ...) {
stopifnot(v.outcome %in% colnames(data));
## prepare data
data[["_id_"]] <- 1:nrow(data)
data <- data[!is.na(data[["_strata_"]]), ]
S <- max(data[["_strata_"]]);
## find mwle
rst_theta <- NULL;
for (i in 1:S) {
cur.d1 <- data[data[["_strata_"]] == i &
data[["_grp_"]] == 1, ]
cur.d0 <- data[data[["_strata_"]] == i &
data[["_grp_"]] == 0, ]
ns1 <- nrow(cur.d1)
ns0 <- nrow(cur.d0)
if (0 == ns1) {
stop(paste("Stratum ", i, " contains no subjects from group 1",
sep = ""))
}
## overall estimate
cur.lambda <- lambdas[i];
cur.theta <- rweWL(cur.data = cur.d1[, v.outcome],
ext.data = cur.d0[, v.outcome],
lambda = cur.lambda, ...)
rst_theta <- rbind(rst_theta,
c(i, NA, NA, ns1, ns0, cur.theta))
for (j in 1:ns1) {
cur.theta <- rweWL(cur.data = cur.d1[-j, v.outcome],
ext.data = cur.d0[, v.outcome],
lambda = cur.lambda, ...)
rst_theta <- rbind(rst_theta,
c(i, 1, cur.d1[j, "_id_"], ns1 - 1, ns0, cur.theta))
}
if (ns0 > 0) {
for (j in 1:ns0) {
cur.theta <- rweWL(cur.data = cur.d1[, v.outcome],
ext.data = cur.d0[-j, v.outcome],
lambda = cur.lambda, ...)
rst_theta <- rbind(rst_theta,
c(i, 0, cur.d0[j, "_id_"], ns1, ns0-1, cur.theta))
}
}
}
colnames(rst_theta) <- c("Strata", "Group", "Inx", "N1", "N0", "Theta")
data.frame(rst_theta)
}
#' PS-Integrated Weighted Likelihood Estimation
#'
#' For all stratum. Variance estimated by Jack-Knife. Works for
#' multiple external RWDs.
#'
#' @param data class D_GPS data frame
#' @param ... parameters for \code{rweWL_G}
#' @param lambdas matrix of number of patients to be borrowed
#'
#' @export
#'
rweGpsWl <- function(data, lambdas, v.outcome = "Y", ...) {
stopifnot(v.outcome %in% colnames(data));
## prepare data
data <- data[!is.na(data[["_strata_"]]), ]
S <- max(data[["_strata_"]])
nd <- max(data[["_grp_"]])
## find mwle
rst_theta <- NULL;
for (i in 1:S) {
cur_d <- list()
ns_d <- NULL
for (j in 1:nd) {
cur_d[[j]] <- data[data[["_strata_"]] == i & data[["_grp_"]] == j, v.outcome]
cur_n <- length(cur_d[[j]])
if (0 == cur_n) {
stop(paste("Stratum ", i, " contains no subjects from group ", "j", sep = ""))
}
ns_d <- c(ns_d, cur_n)
}
cur_lambda <- c(ns_d[1], lambdas[i, ])
## estimate
est_theta <- rweWL_G(cur_d, lambda = cur_lambda, ...)
## jackknife
jk_theta <- NULL;
for (j in 1:nd) {
for (k in 1:ns_d[j]) {
cur_djk <- cur_d
cur_djk[[j]] <- cur_djk[[j]][-k]
cur_theta <- rweWL_G(cur_djk, lambda = cur_lambda, ...)
jk_theta <- c(jk_theta, cur_theta);
}
}
var_theta <- (length(jk_theta)-1) / (length(jk_theta)) * sum((jk_theta - est_theta)^2)
rst_theta <- rbind(rst_theta,
c(est_theta, var_theta, ns_d))
}
##mwle
ws <- rst_theta[, 3] / sum(rst_theta[, 3])
rst_mwle <- sum(ws * rst_theta[, 1])
rst_var <- sum(ws^2 * rst_theta[, 2])
list(mwle = rst_mwle,
var = rst_var,
mwle.strata = rst_theta[, 1],
var.strata = rst_theta[, 2],
ns1 = rst_theta[, 3],
nsj = rst_theta[, -(1:3)])
}
#' Weighted Likelihood Estimation for one external RWD
#'
#' @param cur.data data from current study
#' @param ext.data data from external study
#' @param type type of outcomes
#' @param lambda power parameter without standardization of ns0
#' @param equal.sd boolean. whether sd is the same between the current and external study
#'
#' @export
#'
rweWL <- function(cur.data, ext.data, lambda,
type = c("continuous", "binary"),
equal.sd = TRUE) {
f.ll <- function(pars) {
theta <- pars[1];
sig2.1 <- pars[2];
sig2.0 <- pars[3];
ll <- - n1 * log(sig2.1) / 2;
ll <- ll - n1 * mean((cur.data - theta)^2)/2/sig2.1;
ll <- ll - lambda * log(sig2.0) / 2;
ll <- ll - lambda * mean((ext.data - theta)^2)/2/sig2.0;
ll
}
f.gradient <- function(pars) {
theta <- pars[1];
sig2.1 <- pars[2];
sig2.0 <- pars[3];
g <- numeric(length(pars));
## d logl / d theta
g[1] <- n1/sig2.1*(mean(cur.data) - theta) +
lambda/sig2.0*(mean(ext.data) - theta)
## d logl / d sig2.1
g[2] <- - n1/2/sig2.1 +
n1 * mean((cur.data - theta)^2)/2/sig2.1/sig2.1;
## d logl / d sig2.0
g[3] <- - lambda/2/sig2.0 +
lambda * mean((cur.data - theta)^2)/2/sig2.0/sig2.0;
return(g)
}
type <- match.arg(type);
n1 <- length(cur.data);
## ignore external data
if (0 == lambda) {
ext.data <- cur.data; ## placeholder
equal.sd <- TRUE;
}
init.theta <- (n1 / (n1 + lambda)) * mean(cur.data) +
(lambda/(n1 + lambda)) * mean(ext.data);
if (("continuous" == type & equal.sd) | "binary" == type) {
rst <- init.theta;
} else {
init.sig2.1 <- mean((cur.data - init.theta)^2);
init.sig2.0 <- mean((ext.data - init.theta)^2);
rst <- optim(c(init.theta, init.sig2.1, init.sig2.0),
method = "L-BFGS-B",
fn = f.ll,
lower = c(-Inf, 1e-6, 1e-6), upper = rep(Inf,3),
control= list(fnscale=-1))$par[1];
}
rst;
}
#' Weighted Likelihood Estimation for general number of external RWDs
#'
#' @param dta list of data from each source including current study
#' @param type type of outcomes
#' @param lambda number of patients to be included or borrowed. For current
#' study, it should be equal to the number of patients
#'
#' @export
#'
rweWL_G <- function(lst_dta, lambda, type = c("continuous", "binary")) {
type <- match.arg(type);
all_mean <- NULL
for (i in seq_len(length(lst_dta))) {
all_mean <- c(all_mean, mean(lst_dta[[i]]))
}
sum(all_mean * lambda / sum(lambda))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.