31:T4ad,pp = { return md` --- ## Bayesian Approach: Power Prior > A power prior can be expressed as a product of the weighted likelihood of ${tex`\theta`} conditional on the historical data and a prior distribution on the parameters before any data is observed: > > ${tex.block` \pi (\theta | D_0, {\color{red}\alpha}) \propto \left[L (\theta | D_0)\right]^ {{\color{red}\alpha}} \pi_0 (\theta) \\~\\ `} > * ${tex`\theta`}: parameter of interest > * ${tex`L (\theta | D_0)`}: likelihood function of the external data > * ${tex`\pi_0 (\theta) `}: initial prior distribution for ${tex`\theta`} > * ${tex` {\color{red} \alpha}`}: power prior parameter, ${tex` 0 \le \alpha \le 1`} > * ${tex` {\color{red} \alpha}`}: controls how much external data to borrow > >
> The posterior distribution of ${tex`\theta`} is: > ${tex.block` \pi (\theta | D, D_0, {\color{red}\alpha}) \propto L (\theta | D) \pi (\theta | D_0, {\color{red}\alpha}) \\~\\ \quad~~~~~~~~~~~~~~~~~~~~~~\quad \propto L(\theta | D) \ \underbrace{L(\theta | D_0)^{{\color{red}\alpha}} \pi_0(\theta)}_{\pi (\theta | D_0, \alpha)} \\~\\ `} > * ${tex`D`}: current study > * ${tex`D_0`}: external data ` }32:T661,CL = { return md` --- ## Frequentist Approach: Composite Likelihood
> ${tex.block` L (\theta | D, D_0) = \prod_{i \in D} f(y_i | \theta) \prod_{i \in D_0} f(y_i | \theta)^{{\color{red}\lambda_i}} \\~\\ `} > where ${tex`{\color{red}\lambda_i}`} is nonnegative weight to be chosen, and can be used to discount patient information from external data scource (i.e., ${tex`D_0`}). > * ${tex`{\color{red} \lambda_i = 1}`} if the patient ${tex`i`} is from the investigational study. > * ${tex`{\color{red} 0 \le \lambda_i \le 1}`} if the patient ${tex`i`} is from teh external data. > * ${tex`f()`} represents the density function for a single observation in either current study or the external data. > >
> Based on maximum likelihood estimation (MLE) of ${tex`L (\theta | D, D_0)`}, the point estimation for the event rate (in the binomial case) can be written: > ${tex.block` \hat{\theta} = \frac{\sum_{i \in D} y_i + \sum_{i \in D_0} {\color{red} \lambda_i }y_i }{N + \sum_{i \in D_0} {\color{red} \lambda_i }} \\~\\ `} > * ${tex`N`}: sample size of the current study ${tex`D`} > * ${tex`N_0`}: sample size of the external data ${tex`D_0`} >
> Nominal number of subjects borrowed ${tex`\cong \sum_{i \in D_0} {\color{red} \lambda_i} `}
(Bayesian power prior: ${tex`\alpha N_0`}) >
> For example: > * Set all ${tex`\lambda_i = 0`}, ${tex`\large \hat{\theta} = \frac{\sum_{i \in D} y_i}{N}`} (event rate of current study only)

> * Set all ${tex`\lambda_i = 1`}, ${tex`\large \hat{\theta} = \frac{\sum_{i \in D} y_i + \sum_{i \in D_0} y_i}{N + N_0}`} (pooling current study and external data) ` }33:T500,PS = { return md` --- ## Propensity Score > The propensity score ${tex`e(x)`} for a subject with a vector ${tex`x`} of observed covariates is the conditional probability of receiving treatment ${tex`A (Z = 1)`} rather than treatment ${tex`B (Z = 0)`} given ${tex`x`}: > ${tex.block` e(x) = Pr(Z = 1| x) \\~\\ `} > > propensity score ${tex`e(x)`} is a balancing score in the sense that it is a function of the observed covariates x such that the conditional distribution of ${tex`x`} given ${tex`e(x)`} is the same for subjects who had received treatment A and subjects who had received treatment B. > > If treatment assignment is strongly ignorable, that is, the treatment assignment Z and the outcome Y are conditionally independent given the covariates, ${tex`x`}, then the average treatment effect at each value of the propensity score is an unbiased estimate of the true treatment effect at that propensity score, and therefore, matching, subclassification and covariate adjustment on the propensity score can produce unbiased estimates of the treatment effect. > > The propensity score is generated by modeling the probability of treatment group membership, based on the baseline covariates, for example, by a multiple logistic regression or discriminant analysis. ` }34:T10ea,{ return htl.html`

[R Code: Propensity score estimation]

${highlightRender( 'r', ` # Propensity score estimation dta_ps_single <- psrwe_est(ex_dta, v_covs = paste("V", 1:7, sep = ""), v_grp = "Group", cur_grp_level = "current", ps_method = "logistic", nstrata = 5) ########################################################################### # ※ indicates nested function psrwe_est <- function (data, ps_fml = NULL, ps_method = c("logistic", "randomforest"), v_covs = "V1", v_grp = "Group", cur_grp_level = 1, v_arm = NULL, ctl_arm_level = NULL, stra_ctl_only = TRUE, nstrata = 5, ...) { if (!identical("data.frame", class(data))) { warning("data should be a data.frame object") } data <- as.data.frame(data) ps_method <- match.arg(ps_method) if (is.null(ps_fml)) { ps_fml <- as.formula(paste(v_grp, "~", paste(v_covs, collapse = "+"), sep = "")) } dnames <- colnames(data) v_fml <- all.vars(ps_fml) if (!(all(v_fml %in% dnames))) { stop("Group or covariate was not found in data") } v_grp <- v_fml[1] if (!is.null(v_arm)) { if (!(v_arm %in% dnames)) stop("Arm was not found in data") is_rct <- TRUE } else { is_rct <- FALSE } d1_inx <- cur_grp_level == data[[v_grp]] keep_inx <- which(d1_inx) if (0 == length(keep_inx)) stop("No current study subjects found in data") if (!is.null(v_arm) & stra_ctl_only) { d1_inx <- d1_inx & ctl_arm_level == data[[v_arm]] } all_ps <- get_ps(data, ps_fml = ps_fml, ps_method = ps_method, ...) # ※ data[["_ps_"]] <- all_ps grp <- rep(1, nrow(data)) grp[which(data[[v_grp]] != cur_grp_level)] <- 0 data[["_grp_"]] <- grp arm <- rep(0, nrow(data)) if (is_rct) { arm[which(data[[v_arm]] != ctl_arm_level)] <- 1 } data[["_arm_"]] <- arm if (nstrata > 0) { d1_ps <- all_ps[which(d1_inx)] strata <- rwe_cut(d1_ps, all_ps, breaks = nstrata, keep_inx = keep_inx) # ※ data[["_strata_"]] <- factor(strata, levels = c(1:nstrata), labels = paste("Stratum", 1:nstrata)) } data <- data %>% arrange('_grp_', '_arm_', '_strata_') %>% mutate('_id_' = 1:n()) rst <- list(data = data, ps_fml = ps_fml, ps_method = ps_method, is_rct = is_rct, nstrata = nstrata) class(rst) <- get_rwe_class("DWITHPS") return(rst) } ########################################################################### # calculate prepensity score by strata # by logistic model or randomforest model get_ps <- function (dta, ps_method = c("logistic", "randomforest"), ps_fml, ntree = 5000, grp = NULL, ps_cov = NULL, ...) { type <- match.arg(ps_method) if (is.null(ps_fml)) ps_fml <- as.formula(paste(grp, "~", paste(ps_cov, collapse = "+"), sep = "")) grp <- all.vars(ps_fml)[1] switch(type, logistic = { glm_fit <- glm(ps_fml, family = binomial, data = dta, ...) est_ps <- glm_fit$fitted }, randomforest = { dta[[grp]] <- as.factor(dta[[grp]]) rf_fit <- randomForest(ps_fml, data = dta, ntree = ntree, ...) est_ps <- predict(rf_fit, type = "prob")[, 2] }) return(est_ps) } # multiple stratum cut by quantile of current study data's PS rwe_cut <- function (x, y = x, breaks = 5, keep_inx = NULL) { cuts <- quantile(x, seq(0, 1, length = breaks + 1)) cuts[1] <- cuts[1] - 1e-100 rst <- rep(NA, length(y)) for (i in 2:length(cuts)) { inx <- which(y > cuts[i - 1] & y <= cuts[i]) rst[inx] <- i - 1 } if (!is.null(keep_inx)) { inx <- which(y[keep_inx] <= cuts[1]) if (0 < length(inx)) { rst[keep_inx[inx]] <- 1 } inx <- which(y[keep_inx] > cuts[length(cuts)]) if (0 < length(inx)) { rst[keep_inx[inx]] <- length(cuts) - 1 } } rst } ` )}

` }35:T54a,ps_plot_1 = { const psdata = ps_demo1.sort((a,b) => +b._id_ - +a._id_); const cuts = d3.range(0,1.2,0.2).map(d=>d3.quantile(psdata.filter(d=>d.Group == 'current').map(d=>+d._ps_), d)); const xpos = d3.range(1, cuts.length).map(d=>+d3.mean([cuts[d-1], cuts[d]]).toFixed(2)); const strataNum = [... new Set(psdata.filter(d=>d._strata_ !== 'NA').map(d=>d._strata_))] .sort() .map(i => ({'N_rwd': psdata.filter(d=>(d._strata_ == i & d.Group == 'rwd')).length, 'N_current': psdata.filter(d=>(d._strata_ == i & d.Group !== 'rwd')).length,}) ); return addAnimation( Plot.plot({ height: 380, color: { legend: true }, x:{ label: 'Propensity Score →', }, marks: [ Plot.dotX(psdata, Plot.dodgeY({x: d=>+d._ps_, r:2, fill: "Group", fillOpacity: 0, title:d=>(+d._ps_).toFixed(4)})), Plot.ruleX(cuts, {stroke: "#979797"}), Plot.text(xpos, {x: d=>d, text: (d,i) => 'Stratum' + (i+1) + '\n\nFrom: ' + cuts[i].toFixed(4) + '\nTo: ' + cuts[i+1].toFixed(4) + '\n\nrwd: ' + strataNum[i].N_rwd + '\ncur: ' + strataNum[i].N_current, fontSize: '11px', fontWeight: 'bold' }), ] }), { type: "circle", attribute: "fillOpacity", endValue: 0.6, delay: 5, replay: replay } ) }36:T7c1,{ return htl.html`

[R Code: Obtain discounting parameters]

${highlightRender( 'r', ` ps_bor_single <- psrwe_borrow(dta_ps_single, total_borrow = 30) ############################################################################# # ※ indicates nested function psrwe_borrow <- function (dtaps, total_borrow, method = c("distance", "inverse_distance"), ...) { is_ps <- inherits(dtaps, what = get_rwe_class("DWITHPS")) is_ps_match <- inherits(dtaps, what = get_rwe_class("DPSMATCH")) stopifnot(is_ps | is_ps_match) method <- match.arg(method) rst_sum <- summary(dtaps, ...) ns0 <- rst_sum$Summary$N_RWD rs <- rst_sum$Summary$Distance if (dtaps$is_rct) { ns1 <- rst_sum$Summary$N_Cur_CTL } else { ns1 <- rst_sum$Summary$N_Current } if (is_ps_match) { method <- "n_current" } borrow <- get_aborrow(total_borrow, ns0, ns1, rs, m_lambda = method, ...) # ※ dtaps$Total_borrow <- total_borrow dtaps$Borrow <- cbind(rst_sum$Summary, borrow) dtaps$Borrow_method <- method class(dtaps) <- get_rwe_class("PSDIST") return(dtaps) } ############################################################################# # steps to calculate power prior alpha get_aborrow <- function (total_borrow, ns0, ns1, rs, m_lambda = c("distance", "inverse_distance", "n_current"), ...) { m_lambda <- match.arg(m_lambda) proportion <- switch(m_lambda, distance = { proportion <- rs/sum(rs) }, inverse_distance = { mrs <- 1/(1 - rs) proportion <- mrs/sum(mrs) }, n_current = { proportion <- ns1/sum(ns1) }) borrow <- apply(cbind(ns0, total_borrow * proportion), 1, min) alpha <- borrow/ns0 data.frame(Proportion = proportion, N_Borrow = borrow, Alpha = alpha) } ` )}

` }37:T13a5,table2 = { const names = ['', '', '1', '2','3','4','5','Total']; const cs = ['', 'Current Study (n)', '40', '40','40','40','40','200']; const rwd = ['', 'RWD (n)', '729', '156','78','50','13','1026']; const oc = ['⇳', 'Overlapping Coefficient', '0.561', '0.721','0.804','0.810','0.796','']; const soc = ['⇳', 'Standardized OC', '15.2%', '19.5%','21.8%', '21.9%','21.6%','100%']; const pb = ['⇳', 'Approx. patients borrowed', '4.56', '5.86','6.53', '6.58','6.47','30']; const pp = ['⇳', 'Power parameter', '0.006', '0.038','0.084', '0.132','0.497','']; const table = d3.create("table").attr('id', 'table2').style('width', 800); let bool3 = true; let bool4 = true; let bool5 = true; let bool6 = true; table .join("tr") .selectAll("th") .data(names) .join("th") .text(d => d) .style("background-color", "#aaa") .style("color", "#fff") .style("text-align", "center") const tbody = table.append("tbody"); tbody.selectAll("tr") .data([cs, rwd]) .join("tr") .selectAll("td") .data(row => row) .join("td") .text(d => d) .style("text-align", "center") const row3 = tbody.append("tr") .selectAll("td") .data(oc) .join('td') .text(d=>d) .style("text-align", "center") .style("cursor", "pointer") .style("color", "steelblue") .style("font-weight", "bold") .on('click', function(){ bool3 = !bool3 step1 .style('display', bool3 ? '' : 'none') row3 .style('color', bool3 ? 'steelblue' : '#1b1e23') .style('font-weight', bool3 ? 'bold' : '') }); const step1 = tbody.append('tr') .html('↑
Step 1

Overlapping Coefficients are calculated for each stratum
by using kernel density estimation
') const row4 = tbody.append("tr") .selectAll("td") .data(soc) .join('td') .text(d=>d) .style("text-align", "center") .style("cursor", "pointer") .on('click', function(){ bool4 = !bool4 step2 .style('display', bool4 ? 'none' : '') row4 .style('color', bool4 ? '#1b1e23' : 'steelblue' ) .style('font-weight', bool4 ? '' : 'bold' ) }); const step2 = tbody.append('tr') .style('display', 'none') .html('↑
Step 2

Standardized OC for each stratum

(i.e., 15.2% = 0.561 / (0.561 + 0.721 + 0.804 + 0.810 + 0.796) )
') const row5 = tbody.append("tr") .selectAll("td") .data(pb) .join('td') .text(d=>d) .style("text-align", "center") .style("cursor", "pointer") .on('click', function(){ bool5 = !bool5 step3 .style('display', bool5 ? 'none' : '') row5 .style('color', bool5 ? '#1b1e23' : 'steelblue' ) .style('font-weight', bool5 ? '' : 'bold' ) }); const step3 = tbody.append('tr') .style('display', 'none') .html('↑
Step 3

The number of external patients allocated to each PS stratum
is proportional to their standardized overlapping coefficient

(i.e., 30 x 15.2% = 4.56 )
') const row6 = tbody.append("tr") .selectAll("td") .data(pp) .join('td') .text(d=>d) .style("text-align", "center") .style("cursor", "pointer") .on('click', function(){ bool6 = !bool6 step4 .style('display', bool6 ? 'none' : '') row6 .style('color', bool6 ? '#1b1e23' : 'steelblue' ) .style('font-weight', bool6 ? '' : 'bold' ) }); const step4 = tbody.append('tr') .style('display', 'none') .html('↑
Step 4

Determining how much info to leverage: α
The info leveraged from each individual external patient depends on
how many external patients in that PS stratum.

The power prior parameter for each individual external patient, α, is inversely proportional to the sample size of external patients in the PS stratum.

(i.e., 4.56 / 729 = 0.006 )
') return table.node(); }38:T1f56,{ return htl.html`

[R Code: Stratum-specific posterior distributions]

${highlightRender( 'r', ` rst_pp <- psrwe_powerp(ps_bor_single, outcome_type = "binary", v_outcome = "Y_Bin", seed = 1234) ################################################################ # ※ indicates nested function psrwe_powerp <- function (dta_psbor, v_outcome = "Y", outcome_type = c("continuous", "binary"), prior_type = c("fixed", "random"), ..., seed = NULL) { stopifnot(inherits(dta_psbor, what = get_rwe_class("PSDIST"))) type <- match.arg(outcome_type) prior_type <- match.arg(prior_type) stopifnot(v_outcome %in% colnames(dta_psbor$data)) old_seed <- NULL if (!is.null(seed)) { if (exists(".Random.seed", envir = .GlobalEnv)) { old_seed <- get(".Random.seed", envir = .GlobalEnv) } set.seed(seed) } rst_obs <- get_observed(dta_psbor$data, v_outcome) # ※ lst_dta <- get_stan_data(dta_psbor, v_outcome, prior_type) # ※ stan_mdl <- if_else("continuous" == type, "powerps", "powerpsbinary") ctl_post <- rwe_stan(lst_data = lst_dta$ctl, stan_mdl = stan_mdl, ...) # ※ ctl_thetas <- extract(ctl_post, "thetas")$thetas is_rct <- dta_psbor$is_rct trt_post <- NULL trt_thetas <- NULL if (is_rct) { trt_post <- rwe_stan(lst_data = lst_dta$trt, stan_mdl = stan_mdl, ...) trt_thetas <- extract(trt_post, "thetas")$thetas } rst_trt <- NULL rst_effect <- NULL if (is_rct) { rst_trt <- get_post_theta(trt_thetas, dta_psbor$Borrow$N_Cur_TRT) # ※ rst_effect <- get_post_theta(trt_thetas - ctl_thetas, dta_psbor$Borrow$N_Current) n_ctl <- dta_psbor$Borrow$N_Cur_CTL } else { n_ctl <- dta_psbor$Borrow$N_Current } rst_ctl <- get_post_theta(ctl_thetas, n_ctl) if (!is.null(seed)) { if (!is.null(old_seed)) { invisible(assign(".Random.seed", old_seed, envir = .GlobalEnv)) } else { invisible(rm(list = c(".Random.seed"), envir = .GlobalEnv)) } } rst <- list(stan_rst = list(ctl_post = ctl_post, trt_post = trt_post), Observed = rst_obs, Control = rst_ctl, Treatment = rst_trt, Effect = rst_effect, Borrow = dta_psbor$Borrow, Total_borrow = dta_psbor$Total_borrow, Method = "ps_pp", Outcome_type = type, is_rct = is_rct) class(rst) <- get_rwe_class("ANARST") rst } ################################################################ get_observed <- function (data, v_covs) { data <- data %>% mutate(Stratum = '_strata_', Group = '_grp_', Arm = '_arm_', Y = (!!as.name(v_covs))) %>% filter(!is.na(Stratum)) rst1 <- data %>% group_by(Group, Arm, Stratum) %>% summarize(N = n(), Mean = mean(Y), StdErr = sd(Y)) rst2 <- data %>% mutate(Stratum = "Overall") %>% group_by(Group, Arm, Stratum) %>% summarize(N = n(), Mean = mean(Y), StdErr = sd(Y)) rst <- data.frame(rbind(rst1, rst2)) rst } get_stan_data <- function (dta_psbor, v_outcome, prior_type) { f_curd <- function(i, d1, d0 = NULL) { cur_d <- c(N1 = length(d1), YBAR1 = mean(cur_d1), YSUM1 = sum(cur_d1)) if (is.null(d0)) { cur_d <- c(cur_d, N0 = 0, YBAR0 = 0, SD0 = 0) } else { cur_d <- c(cur_d, N0 = length(cur_d0), YBAR0 = mean(cur_d0), SD0 = sd(cur_d0)) } list(stan_d = cur_d, y1 = d1, inx1 = rep(i, length(d1))) } get_cur_d <- function (data, i, v_covs, grp_cur = 1, arm_ctl = 0) { f_d <- function(x) { if (0 == nrow(x)) { rst <- NULL } else { rst <- x[, v_covs, drop = TRUE] } rst } cur_d1 <- data[data[["_strata_"]] == i & data[["_grp_"]] == grp_cur & data[["_arm_"]] == arm_ctl, v_covs, drop = FALSE] cur_d0 <- data[data[["_strata_"]] == i & data[["_grp_"]] != grp_cur & data[["_arm_"]] == arm_ctl, v_covs, drop = FALSE] cur_d1t <- data[data[["_strata_"]] == i & data[["_grp_"]] == grp_cur & data[["_arm_"]] != arm_ctl, v_covs, drop = FALSE] if (0 == nrow(cur_d1) | 0 == nrow(cur_d0)) { warning(paste(i, "contains no subjects from group 1 or 0")) } out <- list(cur_d1 = f_d(cur_d1), cur_d0 = f_d(cur_d0), cur_d1t = f_d(cur_d1t)) return(out) } is_rct <- dta_psbor$is_rct data <- dta_psbor$data data <- data[!is.na(data[["_strata_"]]), ] strata <- levels(data[["_strata_"]]) nstrata <- length(strata) ctl_stan_d <- NULL ctl_y1 <- NULL ctl_inx1 <- NULL trt_stan_d <- NULL trt_y1 <- NULL trt_inx1 <- NULL for (i in seq_len(nstrata)) { cur_01 <- get_cur_d(data, strata[i], v_outcome) cur_d1 <- cur_01$cur_d1 cur_d0 <- cur_01$cur_d0 cur_d1t <- cur_01$cur_d1t ctl_cur <- f_curd(i, cur_d1, cur_d0) ctl_stan_d <- rbind(ctl_stan_d, ctl_cur$stan_d) ctl_y1 <- c(ctl_y1, ctl_cur$y1) ctl_inx1 <- c(ctl_inx1, ctl_cur$inx1) if (is_rct) { trt_cur <- f_curd(i, cur_d1t) trt_stan_d <- rbind(trt_stan_d, trt_cur$stan_d) trt_y1 <- c(trt_y1, trt_cur$y1) trt_inx1 <- c(trt_inx1, trt_cur$inx1) } } ctl_lst_data <- list(S = nstrata, A = dta_psbor$Total_borrow, RS = as.array(dta_psbor$Borrow$Proportion), FIXVS = as.numeric(prior_type == "fixed"), N0 = as.array(ctl_stan_d[, "N0"]), N1 = as.array(ctl_stan_d[, "N1"]), YBAR0 = as.array(ctl_stan_d[, "YBAR0"]), SD0 = as.array(ctl_stan_d[, "SD0"]), TN1 = length(ctl_y1), Y1 = ctl_y1, INX1 = ctl_inx1, YBAR1 = as.array(ctl_stan_d[, "YBAR1"]), YSUM1 = as.array(ctl_stan_d[, "YSUM1"])) trt_lst_data <- NULL if (is_rct) { trt_lst_data <- list(S = nstrata, A = 0, RS = as.array(dta_psbor$Borrow$Proportion), FIXVS = as.numeric("fixed" == "fixed"), N0 = as.array(trt_stan_d[, "N0"]), N1 = as.array(trt_stan_d[, "N1"]), YBAR0 = as.array(trt_stan_d[, "YBAR0"]), SD0 = as.array(trt_stan_d[, "SD0"]), TN1 = length(trt_y1), Y1 = trt_y1, INX1 = trt_inx1, YBAR1 = as.array(trt_stan_d[, "YBAR1"]), YSUM1 = as.array(trt_stan_d[, "YSUM1"])) } list(ctl = ctl_lst_data, trt = trt_lst_data) } rwe_stan <- function (lst_data, stan_mdl = c("powerps", "powerpsbinary", "powerp"), chains = 4, iter = 2000, warmup = 1000, control = list(adapt_delta = 0.95), ...) { stan_mdl <- match.arg(stan_mdl) stan_rst <- rstan::sampling(stanmodels[[stan_mdl]], data = lst_data, chains = chains, iter = iter, warmup = warmup, control = control, ...) # ※ stan_rst } get_post_theta <- function (thetas, weights) { ws <- weights/sum(weights) samples <- t(thetas) overall <- apply(samples, 2, function(x) sum(x * ws)) means <- apply(samples, 1, mean) sds <- apply(samples, 1, sd) mean_overall <- mean(overall) sd_overall <- sd(overall) list(Stratum_Samples = samples, Overall_Samples = overall, Stratum_Estimate = cbind(Mean = means, StdErr = sds), Overall_Estimate = cbind(Mean = mean_overall, StdErr = sd_overall)) } ` )}

` }39:T1083,{ return htl.html`

[Stan Code: stanmodels]

${highlightRender( 'Stan', ` $powerp S4 class stanmodel 'powerp' coded as follows: // // TYPICAL POWER PRIOR FOR CONTINUOUS DATA // Power parameter A is fixed // data { //existing data int N0; real YBAR0; real SD0; //current data int TN1; real Y1[TN1]; //target borrowing real A; } transformed data { real a0; real sn0; if (N0 > 0) { a0 = 1 < A/N0 ? 1 : A/N0; sn0 = SD0/sqrt(N0*1.0); } else { a0 = 0; sn0 = 1; } } parameters { real theta; real tau1; } model { //prior theta ~ normal(0, 1000); tau1 ~ cauchy(0, 2.5); //likelihood if (N0 > 0) { target += normal_lpdf(YBAR0 | theta, sn0) * a0; } Y1 ~ normal(theta, tau1); } $powerps S4 class stanmodel 'powerps' coded as follows: // // STRATIFIED POWER PRIOR FOR CONTINUOUS DATA // Power parameter As follows // data { int S; //existing data int N0[S]; real YBAR0[S]; real SD0[S]; //current data int N1[S]; int TN1; real Y1[TN1]; int INX1[TN1]; //prior of vs vector[S] RS; //fix vs int FIXVS; //target borrowing real A; } transformed data { row_vector[S] WS1; for (i in 1:S) { WS1[i] = N1[i]; WS1[i] = WS1[i]/TN1; } } parameters { simplex[S] vs; vector[S] thetas; real taus[S]; } transformed parameters { real as[S]; real sds[S]; for (i in 1:S) { if (0 == N0[i]) { as[i] = 0; } else { if (0 == FIXVS) { as[i] = 1 < A * vs[i] / N0[i] ? 1 : A * vs[i] / N0[i]; } else { as[i] = 1 < A * RS[i] / N0[i] ? 1 : A * RS[i] / N0[i]; } } sds[i] = 0 == as[i] ? 0 : SD0[i] / sqrt(as[i] * N0[i]); } } model { //prior if (A > 0) { target += normal_lpdf(YBAR0 | thetas, sds); } else { thetas ~ normal(0, 1000); } vs ~ dirichlet(RS); taus ~ cauchy(0, 2.5); //likelihood Y1 ~ normal(thetas[INX1], taus[INX1]); } generated quantities { real theta; theta = WS1 * thetas; } $powerpsbinary S4 class stanmodel 'powerpsbinary' coded as follows: // // STRATIFIED POWER PRIOR FOR BINARY DATA // Power parameter As follows dirichlet // data { int S; //existing data int N0[S]; real YBAR0[S]; //current data int N1[S]; int YSUM1[S]; //prior of vs vector[S] RS; //fix vs int FIXVS; //target borrowing real A; } transformed data { row_vector[S] WS1; int sn1; sn1 = sum(N1); for (i in 1:S) { WS1[i] = N1[i]; WS1[i] = WS1[i]/sn1; } } parameters { simplex[S] vs; vector[S] thetas; } transformed parameters { real as[S]; real alphas[S]; real betas[S]; for (i in 1:S) { if (0 == N0[i]) { as[i] = 0; } else { if (0 == FIXVS) { as[i] = 1 < A*vs[i]/N0[i] ? 1:A*vs[i]/N0[i]; } else { as[i] = 1 < A*RS[i]/N0[i] ? 1:A*RS[i]/N0[i]; } } alphas[i] = as[i] * N0[i] * YBAR0[i] + 1; betas[i] = as[i] * N0[i] * (1-YBAR0[i]) + 1; } } model { //prior if (A > 0) { target += beta_lpdf(thetas | alphas, betas); } else { thetas ~ uniform(0,1); } vs ~ dirichlet(RS); //likelihood YSUM1 ~ binomial(N1, thetas); } generated quantities { real theta; theta = WS1 * thetas; } ` )}

` }3a:T1d73,{ return htl.html`

[R Code: Outcome analysis]

${highlightRender( 'r', ` oa_pp <- psrwe_outana(rst_pp, mu = 0.4) oa_pp ################################################################ # ※ indicates nested function psrwe_outana <- function (dta_psrst, method_ci = c("wald", "wilson"), conf_type = c("log_log", "plain"), conf_int = 0.95, alternative = c("less", "greater"), mu = 0, ...) { stopifnot(inherits(dta_psrst, what = get_rwe_class("ANARST"))) stopifnot(dta_psrst$Method %in% c("ps_pp", "ps_cl", "ps_km")) outcome_type <- dta_psrst$Outcome_type is_rct <- dta_psrst$is_rct if (is_rct) { type <- "Effect" } else { type <- "Control" } is_km <- dta_psrst$Method == "ps_km" dta_psrst <- psrwe_ci(dta_psrst, method_ci = method_ci[1], conf_type = conf_type[1], conf_int = conf_int[1], ...) # ※ dta_psrst <- psrwe_infer(dta_psrst, alternative = alternative[1], mu = mu[1]) # ※ rst_conf <- list(Method = dta_psrst$Method, Outcome_type = dta_psrst$Outcome_type, Study_type = ifelse(is_rct, "RCT", "single-arm")) if (is_km) { rst_conf$pred_tp <- dta_psrst$pred_tp } rst_conf$CI <- dta_psrst$CI[c("Method_ci", "Conf_type", "Conf_int")] rst_conf$INFER <- dta_psrst$INFER[c("Method_infer", "Alternative", "Mu")] dtype <- dta_psrst$Observed if (is_km) { dtype <- data.frame(dtype) %>% filter(dta_psrst$pred_tp == T) dtype <- data.frame(dtype[, colnames(dtype) != "T"]) } rst_obs <- dtype if (!is_km) { colnames(rst_obs)[colnames(rst_obs) == "StdErr"] <- "SD" } rst_obs$Group <- factor(rst_obs$Group, levels = c(0, 1), labels = c("RWD", "Cur")) if (is_rct) { rst_obs$Arm <- factor(rst_obs$Arm, levels = c(0, 1), labels = c("ctl", "trt")) } else { rst_obs$Arm <- NULL } col_est <- c("Mean", "StdErr") if (is_km) { col_est <- c(col_est, "T") } dtype <- rbind(dta_psrst[[type]]$Stratum_Estimate[, col_est], dta_psrst[[type]]$Overall_Estimate[, col_est]) if (is_km) { dtype <- data.frame(dtype) %>% filter(dta_psrst$pred_tp == T) dtype <- data.frame(dtype[, colnames(dtype) != "T"]) } rst_est <- data.frame(Stratum = c(dta_psrst$Borrow$Stratum, "Overall")) rst_est <- cbind(rst_est, dtype) if (is_rct) { dtype <- rbind(dta_psrst$Treatment$Stratum_Estimate[, col_est], dta_psrst$Treatment$Overall_Estimate[, col_est]) if (is_km) { dtype <- data.frame(dtype) %>% filter(dta_psrst$pred_tp == T) dtype <- data.frame(dtype[, colnames(dtype) != "T"]) } rst_est_trt <- data.frame(Stratum = c(dta_psrst$Borrow$Stratum, "Overall")) rst_est_trt <- cbind(rst_est_trt, dtype) dtype <- rbind(dta_psrst$Control$Stratum_Estimate[, col_est], dta_psrst$Control$Overall_Estimate[, col_est]) if (is_km) { dtype <- data.frame(dtype) %>% filter(dta_psrst$pred_tp == T) dtype <- data.frame(dtype[, colnames(dtype) != "T"]) } rst_est_ctl <- data.frame(Stratum = c(dta_psrst$Borrow$Stratum, "Overall")) rst_est_ctl <- cbind(rst_est_ctl, dtype) } dtype <- rbind(dta_psrst$CI[[type]]$Stratum_Estimate, dta_psrst$CI[[type]]$Overall_Estimate) if (is_km) { dtype <- cbind(dtype, T = c(dta_psrst[[type]]$Stratum_Estimate$T, dta_psrst[[type]]$Overall_Estimate$T)) dtype <- data.frame(dtype) %>% filter(dta_psrst$pred_tp == T) dtype <- data.frame(dtype[, colnames(dtype) != "T"]) } rst_est <- cbind(rst_est, dtype) if (is_rct) { dtype <- rbind(dta_psrst$CI$Treatment$Stratum_Estimate, dta_psrst$CI$Treatment$Overall_Estimate) if (is_km) { dtype <- cbind(dtype, T = c(dta_psrst$Treatment$Stratum_Estimate$T, dta_psrst$Treatment$Overall_Estimate$T)) dtype <- data.frame(dtype) %>% filter(dta_psrst$pred_tp == T) dtype <- data.frame(dtype[, colnames(dtype) != "T"]) } rst_est_trt <- cbind(rst_est_trt, dtype) dtype <- rbind(dta_psrst$CI$Control$Stratum_Estimate, dta_psrst$CI$Control$Overall_Estimate) if (is_km) { dtype <- cbind(dtype, T = c(dta_psrst$Control$Stratum_Estimate$T, dta_psrst$Control$Overall_Estimate$T)) dtype <- data.frame(dtype) %>% filter(dta_psrst$pred_tp == T) dtype <- data.frame(dtype[, colnames(dtype) != "T"]) } rst_est_ctl <- cbind(rst_est_ctl, dtype) } dtype <- rbind(dta_psrst$INFER[[type]]$Stratum_InferProb, dta_psrst$INFER[[type]]$Overall_InferProb) if (is_km) { dtype <- cbind(dtype, T = c(dta_psrst[[type]]$Stratum_Estimate$T, dta_psrst[[type]]$Overall_Estimate$T)) dtype <- data.frame(dtype) %>% filter(dta_psrst$pred_tp == T) dtype <- data.frame(dtype[, colnames(dtype) != "T"]) } if (dta_psrst$Method == "ps_pp") { colnames(dtype) <- "PostPr" } else { colnames(dtype) <- "p-value" } rst_est <- cbind(rst_est, dtype) rst <- dta_psrst rst$OUTANA <- list(Analysis_Setup = rst_conf, Observed_Summary = rst_obs, Analysis_Summary = rst_est) if (is_rct) { rst$OUTANA$RCT_Summary <- list(Treatment = rst_est_trt, Control = rst_est_ctl) } class(rst) <- get_rwe_class("OUTANA") rst } ################################################################ psrwe_ci <- function (dta_psrst, method_ci = c("wald", "wilson"), conf_type = c("log_log", "plain"), conf_int = 0.95, ...) { stopifnot(inherits(dta_psrst, what = get_rwe_class("ANARST"))) stopifnot(dta_psrst$Method %in% c("ps_pp", "ps_cl", "ps_km")) method_ci <- match.arg(method_ci) outcome_type <- dta_psrst$Outcome_type stopifnot(!(method_ci == "wilson" && outcome_type != "binary")) conf_type <- ifelse(outcome_type == "tte", match.arg(conf_type), NA) if (dta_psrst$Method == "ps_pp") { rst_psci <- get_psci_bayesian(dta_psrst, conf_int, ...) } else { rst_psci <- get_psci_freq(dta_psrst, method_ci, conf_type, conf_int, ...) } rst <- dta_psrst rst$CI <- rst_psci return(rst) } psrwe_infer <- function (dta_psrst, alternative = c("less", "greater"), mu = 0) { stopifnot(inherits(dta_psrst, what = get_rwe_class("ANARST"))) stopifnot(dta_psrst$Method %in% c("ps_pp", "ps_cl", "ps_km")) alternative <- match.arg(alternative) if (dta_psrst$Method == "ps_pp") { rst_psinfer <- get_psinfer_bayesian(dta_psrst, alternative, mu) } else { rst_psinfer <- get_psinfer_freq(dta_psrst, alternative, mu) } rst <- dta_psrst rst$INFER <- rst_psinfer return(rst) } ` )}

` }3b:T46e,table4 = { const names = ['Stratum', 'Mean', 'StdErr', 'Lower', 'Upper', 'PostPr']; const res = [ ['Stratum 1', 0.4399064, 0.07328425, 0.3037564, 0.5900138, 0.30325], ['Stratum 2', 0.2884723, 0.06620974, 0.1699958, 0.4257060, 0.94750], ['Stratum 3', 0.2677781, 0.06265874, 0.1536272, 0.3951384, 0.97750], ['Stratum 4', 0.2565654, 0.06184560, 0.1433269, 0.3858618, 0.98275], ['Stratum 5', 0.3089939, 0.06656854, 0.1871603, 0.4470515, 0.90650], ['Overall', 0.3123432, 0.02929484, 0.2568923, 0.3728607, 0.99800], ] const table = d3.create("table"); table.append("thead") .join("tr") .selectAll("th") .data(names) .join("th") .text(d => d) .style("background-color", "#aaa") .style("color", "#fff") .style("text-align", "center") .style('border', '1px solid black') table.append("tbody") .selectAll("tr") .data(res) .join("tr") .selectAll("td") .data(row => row) .join("td") .text(d => d) .style("text-align", "center") .style('border', '1px solid black') return table.node(); }3c:T46a,{ return htl.html`

[R Code: Overlapping coefficient]

${highlightRender( 'r', ` # calculate Overlapping Coefficients uniq_strata <- unique(data[!is.na(data[["_strata_"]]),"_strata_"]) s <- list() for (i in uniq_strata){ inx_ps0 <- i == data[["_strata_"]] & 0 == data[["_grp_"]] inx_ps1 <- i == data[["_strata_"]] & 1 == data[["_grp_"]] n0_01 <- length(which(inx_ps0)) n1_01 <- length(which(inx_ps1)) ps0 <- data[which(inx_ps0), "_ps_"] ps1 <- data[which(inx_ps1), "_ps_"] cov <- c(ps0, ps1) if (length(unique(cov)) <= 10) { all_x <- c(rep(0, length(ps0)), rep(1, length(ps1))) pt <- apply(prop.table(table(cov, all_x), 2), 1, min) res <- sum(pt) } mn <- pmax(0, min(cov) - 1e-3) mx <- pmin(1, max(cov) + 1e-3) f1 <- approxfun(density(ps1, from = mn, to = mx, bw = "nrd")) f0 <- approxfun(density(ps0, from = mn, to = mx, bw = "nrd")) cur_dist <- try(integrate(function(x) pmin(f1(x), f0(x)),lower = mn, upper = mx, subdivisions = 500)$value) s <- rbind(s, c(i, n0_01, n1_01, cur_dist)) } ` )}

` }3d:T592,{ let w = 800; let xmin0 = pars1_3.mean0 - 5.2; let xmax0 = pars1_3.mean0 + 5.2; let xmin1 = pars1_3.mean1 - 5.2; let xmax1 = pars1_3.mean1 + 5.2; let xmin = d3.min([xmin0, xmin1]); let xmax = d3.max([xmax0, xmax1]); let xs = d3.range( d3.min([pars1_3.mean0 - 3*pars1_3.sd0, pars1_3.mean1 - 3*pars1_3.sd1,]), d3.max([pars1_3.mean0 + 3*pars1_3.sd0, pars1_3.mean1 + 3*pars1_3.sd1,]), 0.01 ); let f0 = xs.map(x => ({x:x, y: jstat.normal.pdf(x, pars1_3.mean0, pars1_3.sd0 ) })); let f1 = xs.map(x => ({x:x, y: jstat.normal.pdf(x, pars1_3.mean1, pars1_3.sd1 ) })); let minf0f1 = f0.map((d,i) => ({x: d.x, y: d.y <= f1[i].y ? d.y : f1[i].y})); let rs = d3.sum(minf0f1.map(d=>d.y)) / 100; let p = plotter({ xDomain: [xmin, xmax], yDomain: [0.5, -0.05], grid: false, width: w, height: w / 2, xLabel: '', yLabel: '', yTickCount: 5, xTickCount: 5, y: { axis: null } }); p.func(x => jstat.normal.pdf(x, pars1_3.mean0, pars1_3.sd0 ), { stroke: 'orange', width: 2, opacity: 0.7, dash: '2,2' }); p.func(x => jstat.normal.pdf(x, pars1_3.mean1, pars1_3.sd1 ), { stroke: 'lightblue', width: 2, opacity: 0.7, }); p.polygon(minf0f1, { fill: '#cc6633', width: 0.2, opacity: 0.4 }); p.text(`Overlapping Coefficient = ` + rs.toFixed(4), pars1_3.mean1, -0.07, { stroke: '#333' }); return p.node; }3e:T445,refs = { return md` --- ## References > * Lilly Q. Yue (2007). Statistical and regulatory issues with the application of propensity score analysis to nonrandomized medical device clinical studies. Journal of Biopharmaceutical Statistics, 17: 1–13, 2007 > * Wang, C., Li, H., Chen, W. C., Lu, N., Tiwari, R., Xu, Y., & Yue, L. Q. (2019). Propensity score-integrated power prior approach for incorporating real-world evidence in single-arm clinical studies. Journal of biopharmaceutical statistics, 29(5), 731-748. > * Wang, C., Lu, N., Chen, W. C., Li, H., Tiwari, R., Xu, Y., & Yue, L. Q. (2020). Propensity score-integrated composite likelihood approach for incorporating real-world evidence in single-arm clinical studies. Journal of biopharmaceutical statistics, 30(3), 495-507. > * Chen, W.C., Wang, C., Li, H., Lu, N., Tiwari, R., Xu, Y. and Yue, L.Q., 2020. Propensity score-integrated composite likelihood approach for augmenting the control arm of a randomized controlled trial by incorporating real-world data. Journal of Biopharmaceutical Statistics, 30(3), pp.508-520. ` } 3f:[{"id":"2fc3bbd9371ccb50cc73ed447404d0486d40af128b86ead9b32c8ae9b8d7272416c7e99dbe91cffbe1d4c2fa0e094c6a7463f4a69b5e98dcf23998b56e1c27f6","href":"https://static.observableusercontent.com/files/2fc3bbd9371ccb50cc73ed447404d0486d40af128b86ead9b32c8ae9b8d7272416c7e99dbe91cffbe1d4c2fa0e094c6a7463f4a69b5e98dcf23998b56e1c27f6","name":"ps_demo1.csv","type":"text/csv","lastModified":1655509398125.047,"size":146216}] 2f:{"title":"Stat Notes | Leveraging External Evidence for Augmenting clinical Study","cells":[{"id":0,"value":"# Stat Notes | Leveraging External Evidence for Augmenting clinical Study\n\n> Propensity score-integrated power prior and Propensity score-integrated composite likelihood (by Chenguang Wang, Heng Li, Wei-Chen Chen, Nelson Lu, Ram Tiwari, Yunling Xu & Lilly Q. Yue) are very intuitive and easy-to-implement approaches for leveraging real-world data in regulatory clinical studies. ","pinned":false,"hidden":false,"mode":"md"},{"id":10,"value":"---\n> According to FDA website, **Real-world data (RWD)** and **real-world evidence (RWE)** are playing an increasing role in health care decisions.\n> * Real-world data are the data relating to patient health status and/or the delivery of health care routinely collected from a variety of sources.\n> * Real-world evidence is the clinical evidence regarding the usage and potential benefits or risks of a medical product derived from analysis of RWD.\n\n> The 21st Century Cures Act, passed in 2016, placing additional focus on the use of external data to support regulatory decision making of medical products.","pinned":false,"hidden":false,"mode":"md"},{"id":136,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":142,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":205,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":226,"value":"{\n return md`\n---\n## Case Study: Leverage RWE in a single-arm study\n\n> This example is based on the R package \"psrwe\". Propensity Score-integrated methods are included in the package for incorporating external data in clinical studies.\n\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":963,"value":"{\n return htl.html`\n

[R Code: Environment Setup]

\n ${highlightRender(\n 'r',\n `\n install.packages(\"psrwe\") \n library(psrwe)\n\n # \"psrwe\" requires \"rstan\" and \"Rcpp\" packages\n`\n)}\n

\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":967,"value":"{\n return md`\n\n> * This is a **single-arm study**. \n> * A total of **1031** RWD subjects and **200** current study subjects \n\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":236,"value":"table1 = Inputs.table(\naq.from(ps_demo1)\n .select([\"Group\", \"Y_Bin\", \"Y_Con\", \"Y_Surv\", \"Status\", \"V1\", \"V2\",\"V3\",\"V4\",\"V5\",\"V6\",\"V7\"])\n .objects()\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":674,"value":"{\n return md`\n

\n### Study Design - 1st Stage\n
\n\n> * Primary endpoint (adverse event, death rate ...)\n> * Associated hypothesis testing
${ tex`\\small H_0 \\ge 0.4`} vs ${tex`\\small H_a < 0.4`} \n\n\n> * Study success criterion \n> * Posterior probability of ${ tex`\\small \\theta < 0.4`} is greater than 0.95 \n> * or p-value < 0.05 in frequentist setting\n\n\n> * Sample size determination \n> * Assume ${ tex`\\small \\theta`} = 0.32 \n> * Set power = 80% and significance level = 0.05 \n> * Then, N = 230\n\n\n> Based on clinical decision (case-by-case basis)\n\n> * Enroll **200** patients in the current investigational study\n> * Borrow **30** (about 13%) patients from the external data\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":716,"value":"{\n return md`\n

\n### Study Design - 2st Stage\n
\n\n> **Propensity score strata**\n\n> Propensity score is used to select patients in ${tex`D_0`} (e.g., from a registry) that are more relevant to the patients in the current study ${tex`D_1`}\n
\n\n> * Propensity scores are estimated by logistic model
(${tex`\\small Group \\sim V1 + V2 + V3 + V4 + V5 + V6 + V7`}) \n
\n> * **5** strata are created based on the quantiles of the current data's propensity scores.\n> * A total of **5** RWD subjects are trimmed and excluded from the final analysis.\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":982,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":271,"value":"viewof replay = Inputs.button(\"Replay\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":259,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":726,"value":"{\n return md`\n

\n\n> **How to determine the discounting parameter ${tex`\\small {\\color{red} \\alpha}`} **\n\n> The table below shows:\n\n> * Leverage all **1026** registry patients\n> * Borrow partial information from each of 1026 registry patients (${tex`\\small {\\color{red} \\alpha < 1}`}) such that the nominal number of patients to be borrowed: **30**.\n> * **Steps** to determine the discounting parameter ${tex`\\small {\\color{red} \\alpha}`} for each registry patient\n\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":988,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":992,"value":"{\n return md`\n> > #### _Click (toggle) each row with the sign **\"⇳\"** to view the detailed process_\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":699,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":884,"value":"{\n return md`\n

\n\n## Outcome \n\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":995,"value":"$38","pinned":false,"hidden":false,"mode":"ojs"},{"id":1068,"value":"$39","pinned":false,"hidden":false,"mode":"ojs"},{"id":1054,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs"},{"id":997,"value":"{\n return md`\n\n> Observed Data Summary\n\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":887,"value":"table3 = {\n const names = ['Group', 'N', 'Mean','SD'];\n const rwd = ['RWD', '1026', '0.3187135', '0.4662049'];\n const cur = ['Current', '200', '0.3000000', '0.4594075'];\n \n const table = d3.create(\"table\");\n \n table.append(\"thead\")\n .join(\"tr\")\n .selectAll(\"th\")\n .data(names)\n .join(\"th\")\n .text(d => d)\n .style(\"background-color\", \"#aaa\") \n .style(\"color\", \"#fff\") \n .style(\"text-align\", \"center\") \n .style('border', '1px solid black')\n \n \n table.append(\"tbody\")\n .selectAll(\"tr\")\n .data([rwd, cur])\n .join(\"tr\")\n .selectAll(\"td\")\n .data(row => row)\n .join(\"td\")\n .text(d => d)\n .style(\"text-align\", \"center\") \n .style('border', '1px solid black')\n \n return table.node();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":896,"value":"{\n return md`\n\n
\n> Analysis Results\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":902,"value":"$3b","pinned":false,"hidden":false,"mode":"ojs"},{"id":905,"value":"{\n return md`\n
\n> With a total of **30** subject borrowed from the RWD, \n\n> * The point estimate is **0.312** with standard error **0.029**\n> * The posterior probability of ${tex`\\theta < 0.4`} is **99.8%**\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":339,"value":"{\n return md`\n

\n## Supplement _(back to study design)_\n
\n> > **Overlapping coefficient**\n\n> > The overlapping coefficient ${tex`r_s`} can be interpreted as the overlapping area of the two density curves ${tex`f_{s,0}`} and ${tex`f_{s,1}`}:\n> ${tex.block`\nr_s = \\int_{0}^{1} min \\left[f_{s,0}(e), f_{s,1}(e) \\right]de. \\\\~\\\\\n`}\n> > ${tex`r_s`} can be equivalently defined as:\n> > ${tex.block`\nr_s = \\int_{0}^{1} min \\left[1, \\frac{d}{du} \\ F_{s,0}(F_{s,1}^{-1}(u)) \\right]du. \\\\~\\\\\n`}\n> > In practice, we can estimate ${tex`r_s`} based on a kernel density estimation of ${tex`f_{s,0}`} and ${tex`f_{s,1}`}, or a kernel estimation of ${tex`F_{s,0}(F_{s,1}^{-1}(u))`} with its corresponding derivative when such derivative exists. Below is a demo to show the similarity changes with distribution change. \n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1028,"value":"$3c","pinned":false,"hidden":false,"mode":"ojs"},{"id":763,"value":"{\n return md`\n

\n> > **Demo: ${tex`r_s`} estimation**\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":366,"value":"viewof pars1_3 = Inputs.form({\n mean0: Inputs.range([1, 7], {step: 0.1, value: 2, label: tex`f_{s,0} \\ \\text{Mean: }(\\mu_0)`}),\n mean1: Inputs.range([1, 5], {step: 0.1, value: 1, label: tex`f_{s,1} \\ \\text{Mean: }(\\mu_1)`}),\n sd0: Inputs.range([0.1, 4], {step: 0.1, value: 2, label: tex`f_{s,0} \\ \\text{SD: }(\\sigma_0)`}),\n sd1: Inputs.range([0.1, 4], {step: 0.1, value: 1, label: tex`f_{s,1} \\ \\text{SD: }(\\sigma_1)`}),\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":360,"value":"$3d","pinned":false,"hidden":false,"mode":"ojs"},{"id":412,"value":"{\n return md`\n

\n> > **PS distribution ( ■current vs ■RWE) by strata**\n\n> > Overlapping Coefficient measures the similarity between the PS distributions.\n`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":333,"value":"plot_splot_2 = {\n\n const psdata = ps_demo1;\n const data = [\n ...psdata.slice().filter(d=>d._strata_ != 'NA').map(d=>({...d, 'strata':'Overall', 'ps': +d._ps_})),\n ...psdata.slice().filter(d=>d._strata_ != 'NA').map(d=>({...d, 'strata':d._strata_, 'ps': +d._ps_}))];\n\n \n return vl.markArea()\n .data(data)\n .transform(\n vl.density('ps').groupby(['strata', 'Group'])\n )\n .encode(\n vl.x().fieldQ(\"value\").axis({grid:false}).title('Propensity Score'), \n vl.y().fieldQ(\"density\").axis(null),\n vl.row().fieldN(\"strata\"),\n vl.color().fieldN(\"Group\").legend(null),\n vl.stroke().fieldN(\"black\").legend(null),\n vl.opacity().value(0.5),\n )\n .height(100)\n .width(width * 0.7)\n .render()\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":217,"value":"$3e","pinned":false,"hidden":false,"mode":"ojs"},{"id":13,"value":"---\n## Appendix","pinned":false,"hidden":false,"mode":"md"},{"id":244,"value":"ps_demo1 = FileAttachment(\"ps_demo1.csv\").csv()","pinned":true,"hidden":false,"mode":"ojs"},{"id":275,"value":"import {addAnimation} from \"@mkfreeman/plot-animation\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":357,"value":"jstat = require('jstat')","pinned":true,"hidden":false,"mode":"ojs"},{"id":354,"value":"import { plotter } from '50dadfdec01c15a8'","pinned":true,"hidden":false,"mode":"ojs"},{"id":457,"value":"html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":920,"value":"highlightRender = (language, src) => {\n let pre = document.createElement('pre');\n pre.innerHTML = `\n ${src}\n `;\n highlight.highlightBlock(pre);\n return pre;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":918,"value":"hljsCSS = html`\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":916,"value":"highlight = require('https://cdn.jsdelivr.net/gh/highlightjs/cdn-release@9.18.5/build/highlight.min.js')","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W3f","stdlib":"1","titling":"auto"}