In biomedical research, biomarker is an important tool to characterize subpopulations with better/worse disease progression, or different response to a treatment/intervention. When comparing two different treatments, to estimate the treatment benefit within a biomarker subpopulation, ideally the biomarker should be collected in two comparative cohorts with the same procedure - for example, collecting biomarker as a baseline stratification factor in a randomized trial.
However, there are a number of cases where it is impossible to collect biomarker samples from both cohorts. For example, the biomarker positive/negative status is only available in the treatment cohort but not observable in the control cohort. In this case, if the treatment benefit in the biomarker positive group is of interest, one needs to compare the biomarker positive group from the treatment cohort against the “unobserved” biomarker positive group from the control cohort. The “unobserved” biomarker positive group is the subjects who are likely to be biomarker positive if the biomarker is measurable. To predict the “unobserved” biomarker positive group, one may predict the hidden status based on other observed associated features.
This package implements multiple approaches for this type of modeling. In addition to the modeling (session 1), the package also provides handy functions on model fitting diagnostics (session 3 and 4), visualization (session 2, 5), and summarization (session 6).
PropensitySub
OS
, OS.event
)Arm
)STRATUM
)library(PropensitySub)
library(dplyr)
head(biomarker, 3)
Patient.ID Sample.ID Arm Weight ECOG Sex Age OS OS.event
1 PID040 SID040 Control 40.0 1 F 65 5.059548 0
2 PID112 SID112 Experimental 45.8 2 F 54 2.398357 1
3 PID079 SID079 Experimental 46.0 1 M 41 10.151951 0
STRATUM STRATUM.next
1 <NA> <NA>
2 Positive Negative
3 Negative Negative
# Control arm subjects have no STRATUM measurments
# Experimental arm subjects partially miss STRATUM measurement.
table(biomarker$Arm, biomarker$STRATUM, useNA = "ifany")
Negative Positive <NA>
Control 0 0 77
Experimental 115 54 6
STRATUM
.Below shows example of 1) imputing missing as negative (assuming missing not at random, and evaluating an extreme case where all missing data are supposed to be negative); 2) treating missing as a separate class: also assuming missing not at random.
indicator_1
: numeric version of STRATUM
where missing in Experimental Arm is imputed as “Negative”.
indicator_2
: numeric version of STRATUM
where missing is kept as a seperate stratum
<- biomarker %>%
biomarker mutate(
indicator_1 = case_when(
== "Positive" ~ 1,
STRATUM == "Negative" ~ 0,
STRATUM # impute missing as "Negative" in Experimental Arm
== "Experimental" & is.na(STRATUM) ~ 0,
Arm is.na(STRATUM) & Arm == "Control" ~ -1
),indicator_2 = case_when(
== "Positive" ~ 1,
STRATUM == "Negative" ~ 0,
STRATUM # keep missing as a seperate stratum in Experimental Arm
== "Experimental" & is.na(STRATUM) ~ 2,
Arm is.na(STRATUM) & Arm == "Control" ~ -1
),# treatment group needs to be factor
Arm = factor(Arm, levels = c("Control", "Experimental"))
)
For subjects in the control arm, to predict a subject’s likelihood of being in different biomarker strata, the first step is calculating P(A | X), where A indicates the biomarker strata and X indicates observed associated features. The model may be learned using the treatment data where both biomarker strata and associated observed features were measured, and be applied in control cohort data to predict the biomarker strata assignments.
Once P(A|X) is calculated, predicted biomarker group of interest in the control cohort can be obtained by inverse probability weighting (IPW) or propensity score matching (PSM). Both IPW and PSM are implemented in the package.
In IPW, the predicted probability of being stratum A will be used as weights when estimating treatment difference of two arms (Hazard ratio for survival endpoint; response rate difference for binary endpoint).
In PSM, the predicted probability of being stratum A will be used for propensity score matching. The matching results will then be used to estimate treatment difference of two arms (Hazard ratio for survival endpoint; response rate difference for binary endpoint).
Below shows different modeling approaches implemented in this package.
plain
: binomial/multinomial glm
# `plain` model with IPW method and aggressive imputation where missing is imputated as negative
<- ipw_strata(
ipw_plain_str2 data.in = biomarker, formula = indicator_1 ~ ECOG + Sex + Age, model = "plain",
indicator.var = "indicator_1", tte = "OS", event = "OS.event", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0)
)# get weighted HRs
$stat ipw_plain_str2
Coefficient Variance HR CI.Lower CI.Upper
Positive -0.4165353 0.3816826 0.6593273 0.2708684 1.6048840
Negative -1.3268812 0.2363171 0.2653034 0.1086968 0.6475431
# check model converge
$converged ipw_plain_str2
[1] TRUE
# get weights
$data %>%
ipw_plain_str2::select(Patient.ID, Arm, STRATUM, ECOG, Sex, Age, indicator_1, pred1, pred0) %>%
dplyrhead()
Patient.ID Arm STRATUM ECOG Sex Age indicator_1 pred1 pred0
1 PID040 Control <NA> 1 F 65 -1 0.3312946 0.6687054
2 PID388 Control <NA> 1 M 41 -1 0.2199562 0.7800438
3 PID087 Control <NA> 0 M 37 -1 0.3005256 0.6994744
4 PID134 Control <NA> 0 F 59 -1 0.4266195 0.5733805
5 PID315 Control <NA> 0 M 57 -1 0.3317387 0.6682613
6 PID062 Control <NA> 1 M 42 -1 0.2211980 0.7788020
# `plain` model with IPW method and missing is kept as another stratum
<- ipw_strata(
ipw_plain_str3 data.in = biomarker, formula = indicator_2 ~ ECOG + Sex + Age, model = "plain",
indicator.var = "indicator_2", tte = "OS", event = "OS.event", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0, "Unknown" = 2)
)# get weighted HRs
$stat ipw_plain_str3
Coefficient Variance HR CI.Lower CI.Upper
Positive -0.4195829 0.3808727 0.6573210 0.27004073 1.6000210
Negative -1.4210601 0.2660317 0.2414579 0.09444277 0.6173254
Unknown -0.4629739 2.6612698 0.6294090 0.09003407 4.4000645
# `plain` model with PSM method and aggressive imputation
<- ps_match_strata(
ps_plain_str2 data.in = biomarker, formula = indicator_1 ~ ECOG + Sex + Age, model = "plain",
indicator.var = "indicator_1", tte = "OS", event = "OS.event", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0)
)# get weighted HRs
$stat ps_plain_str2
Coefficient Variance HR CI.Lower CI.Upper
Positive -0.5140953 0.2599806 0.5980414 0.2406342 1.4862952
Negative -1.1861074 0.1915910 0.3054078 0.1347383 0.6922601
dwc
: doubly weighted control# dwc model with IPW method
<- ipw_strata(
ipw_dwc data.in = biomarker, formula = indicator_2 ~ ECOG + Sex + Age, model = "dwc",
indicator.var = "indicator_2", tte = "OS", event = "OS.event", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0)
)# get weighted HRs
$stat ipw_dwc
Coefficient Variance HR CI.Lower CI.Upper
Positive -0.4221217 0.3804891 0.6556542 0.26943319 1.5955067
Negative -1.4201333 0.2660816 0.2416818 0.09452841 0.6179104
# dwc model with PSM method
<- ps_match_strata(
ps_dwc data.in = biomarker, formula = indicator_2 ~ ECOG + Sex + Age, model = "dwc",
indicator.var = "indicator_2", tte = "OS", event = "OS.event", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0, "Missing" = 2)
)# get weighted HRs
$stat ps_dwc
Coefficient Variance HR CI.Lower CI.Upper
Positive -0.3724204 0.2479238 0.6890645 0.28461633 1.6682454
Negative -1.3035235 0.2154364 0.2715732 0.11870059 0.6213281
Missing 0.0000000 2.0000000 1.0000000 0.07255041 13.7835189
wri
: weight regression imputation, currently only available with IPW method# data process: create a numeric version of next STRATUM to learn from
<- biomarker %>%
biomarker mutate(
indicator_next_2 = case_when(
== "Positive" ~ 1,
STRATUM.next == "Negative" ~ 0,
STRATUM.next == "Experimental" & is.na(STRATUM.next) ~ 2,
Arm is.na(STRATUM.next) & Arm == "Control" ~ -1
)
)<- ipw_strata(
ipw_wri data.in = biomarker, formula = indicator_2 ~ ECOG + Sex + Age, model = "wri",
indicator.var = "indicator_2", indicator.next = "indicator_next_2",
tte = "OS", event = "OS.event", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0)
)# get weighted HRs
$stat ipw_wri
Coefficient Variance HR CI.Lower CI.Upper
Positive -0.4361789 0.2539934 0.6465020 0.2764101 1.5121186
Negative -1.3953601 0.3072783 0.2477438 0.1011774 0.6066276
The treatment benefit within biomarker positive stratum can be calculated as the hazard ratio between the observed biomarker positive group in the treatment cohort against the propensity score matched biomarker positive group in the control cohort. An associated KM can be generated as well. When IPW is used, P(A|X) are used as weights for the HR calculation and KM curves. Binary endpoint is also supported.
# for ipw_plain model results
km_plot_weight(
data.in = ipw_plain_str2$data,
indicator.var = "indicator_1", tte = "OS", event = "OS.event", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0)
)
[[1]]
[[2]]
# to get weights from model result for further usage
$data %>%
ipw_plain_str2::select(Patient.ID, Arm, STRATUM, ECOG, Sex, Age, indicator_1, pred1, pred0) %>%
dplyrhead()
Patient.ID Arm STRATUM ECOG Sex Age indicator_1 pred1 pred0
1 PID040 Control <NA> 1 F 65 -1 0.3312946 0.6687054
2 PID388 Control <NA> 1 M 41 -1 0.2199562 0.7800438
3 PID087 Control <NA> 0 M 37 -1 0.3005256 0.6994744
4 PID134 Control <NA> 0 F 59 -1 0.4266195 0.5733805
5 PID315 Control <NA> 0 M 57 -1 0.3317387 0.6682613
6 PID062 Control <NA> 1 M 42 -1 0.2211980 0.7788020
# for ipw_wri model results
km_plot_weight(
data.in = ipw_wri$data,
indicator.var = "indicator_2", tte = "OS", event = "OS.event", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0)
)
[[1]]
[[2]]
# for ps_dwc model results
km_plot_weight(
data.in = ps_dwc$data,
indicator.var = "indicator_2", tte = "OS", event = "OS.event", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0, "Missing" = 2)
)
[[1]]
[[2]]
[[3]]
The package also implemented the absolute standardized mean difference (ASMD) calculation to assess balance of the observed associated features (Xs)(Austin and Stuart, 2015). With a good model fitting, after IPW or PSM adjustment, one would expect a small number of features’ or no feature’s ASMD exceed certain thresholds (such as 0.1 and 0.25). One would expect to see smaller values comparing after-adjustment ASMDs to pre-adjustment ASMDs.
Theoretical calculation is also implemented. The theoretical calculation provides the expected number of features with ASMD greater than a threshold. The calculation depends on sample size and the number of features in X. One may compare the empirical number of features whose ASMD beyond a threshold to the theoretical calculation to assess model performance.
ipw
method with plain
model<- std_diff(
ipw_plain_diff data.in = ipw_plain_str2$data, vars = c("ECOG", "Sex", "Age"),
indicator.var = "indicator_1", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0),
usubjid.var = "Patient.ID"
)$Positive ipw_plain_diff
var absolute_std_diff type
1 ECOG. 0.12302468 adjusted difference
2 Sex.M 0.04787949 adjusted difference
3 Age. 0.15835312 adjusted difference
4 ECOG. 0.01674508 unadjusted difference
5 Sex.M 0.16075410 unadjusted difference
6 Age. 0.15407885 unadjusted difference
$Negative ipw_plain_diff
var absolute_std_diff type
1 ECOG. 0.18839170 adjusted difference
2 Sex.M 0.02334638 adjusted difference
3 Age. 0.12044671 adjusted difference
4 ECOG. 0.25462615 unadjusted difference
5 Sex.M 0.02819748 unadjusted difference
6 Age. 0.12228381 unadjusted difference
# Visualize differences
std_diff_plot(ipw_plain_diff, legend.pos = "bottom")
$Positive
$Negative
psm
method with dwc
model<- std_diff(
ps_dwc_diff data.in = ps_dwc$data, vars = c("ECOG", "Sex", "Age"),
indicator.var = "indicator_2", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0, "Missing" = 2),
usubjid.var = "Patient.ID"
)$Missing ps_dwc_diff
var absolute_std_diff type
1 ECOG. 0.0000000 adjusted difference
2 Sex.M 0.0000000 adjusted difference
3 Age. 0.2412968 adjusted difference
4 ECOG. 0.0000000 unadjusted difference
5 Sex.M 0.0000000 unadjusted difference
6 Age. 0.2412968 unadjusted difference
# Visualize differences
std_diff_plot(ps_dwc_diff)
$Positive
$Negative
$Missing
<- c("ECOG", "Sex", "Age")
vars <- c(0.15, 0.2)
thresholds <- list("Positive" = 1, "Negative" = 0)
class_int_list <- 2 # randomization ratio
rand_ratio <- lapply(class_int_list, function(x) nrow(biomarker %>% filter(indicator_1 %in% x)))
n_arms <- sapply(n_arms, function(x) {
exp_diff expected_feature_diff(
n.feature = length(vars),
n.arm1 = x,
n.arm2 = x / rand_ratio,
threshold = thresholds
)%>% t()
})
# Expected imbalanced features
exp_diff
Expected # features > 0.15 Expected # features > 0.2
Positive 1.579074 1.1961231
Negative 1.026179 0.6170036
# Calculate the observed imbalanced features in model ipw_plain_diff
<- sapply(thresholds, function(th){
obs_diff_cnt sapply(ipw_plain_diff, function(gp){
<- as.character(subset(gp, type=="adjusted difference" & absolute_std_diff>=th)$var)
ft length(unique(sapply(ft, function(ff)strsplit(ff, split="\\.")[[1]][1])))
})
})colnames(obs_diff_cnt) <- paste0("Observed # features > ", thresholds)
rownames(obs_diff_cnt) <- names(ipw_plain_diff)
# the number of observed imbalanced features for each threshold
# Compare expected to observed # of imbalanced features to check model fit
obs_diff_cnt
Observed # features > 0.15 Observed # features > 0.2
Positive 1 0
Negative 1 0
<- sapply(thresholds, function(th) {
obs_diff_fac sapply(ipw_plain_diff, function(gp) {
<- as.character(subset(gp, type=="adjusted difference" & absolute_std_diff>=th)$var)
ft paste(ft, collapse=",")
})
})colnames(obs_diff_fac) <- paste0("features > ", thresholds)
rownames(obs_diff_fac) <- names(ipw_plain_diff)
# the observed individual features that are imbalanced for each threshold
obs_diff_fac
features > 0.15 features > 0.2
Positive "Age." ""
Negative "ECOG." ""
ipw
method with plain
model<- bootstrap_propen(
boot_ipw_plain data.in = biomarker, formula = indicator_1 ~ ECOG + Sex + Age,
indicator.var = "indicator_1", tte = "OS", event = "OS.event", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0),
estimate.res = ipw_plain_str2, method = "ipw", n.boot = 100
)# get bootstrap CI
$est.ci.mat boot_ipw_plain
Estimate Bootstrap CI low Bootstrap CI high
Positive -0.4165353 -1.283337 0.4502662
Negative -1.3268812 -2.271757 -0.3820054
HR ratio: strata over Positive
Positive 1.000000
Negative 0.402385
# summary statistics from bootstraps
$boot.out.est boot_ipw_plain
Mean Var Median Min
Positive;Coefficient -0.3540341 0.195588381 -0.3698299 -1.683188353
Positive;Variance 0.4310293 0.013133483 0.4027040 0.243995859
Positive;HR 0.7741164 0.134388387 0.6908575 0.185780696
Positive;CI.Lower 0.3123257 0.022428198 0.2830601 0.048395835
Positive;CI.Upper 1.9439027 0.911859777 1.6892788 0.693508873
Negative;Coefficient -1.3246478 0.232409179 -1.2919873 -3.835899987
Negative;Variance 0.2752179 0.012857275 0.2478241 0.145512877
Negative;HR 0.2926887 0.014533128 0.2747243 0.021581906
Negative;CI.Lower 0.1186285 0.002994016 0.1128095 0.002791294
Negative;CI.Upper 0.7330693 0.077532519 0.6865619 0.166868400
Positive over Positive 1.0000000 0.000000000 1.0000000 1.000000000
Negative over Positive 0.4372193 0.061767413 0.4234440 0.042909461
Max VarLog
Positive;Coefficient 0.9218116 NA
Positive;Variance 0.8232510 0.06302566
Positive;HR 2.5138404 0.19558838
Positive;CI.Lower 0.9368164 0.23989078
Positive;CI.Upper 6.7456048 0.18402106
Negative;Coefficient -0.3512042 NA
Negative;Variance 1.0681811 0.08751756
Negative;HR 0.7038400 0.23240918
Negative;CI.Lower 0.2927817 0.38615749
Negative;CI.Upper 1.6920143 0.14463584
Positive over Positive 1.0000000 0.00000000
Negative over Positive 1.9730948 0.32326143
# error status and convergence status
$error.est boot_ipw_plain
[1] 0
$conv.est boot_ipw_plain
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
ipw
method with wri
model<- bootstrap_propen(
boot_ipw_wri data.in = biomarker, formula = indicator_2 ~ ECOG + Sex + Age,
indicator.var = "indicator_2", indicator.next = "indicator_next_2",
tte = "OS", event = "OS.event", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0),
estimate.res = ipw_wri, method = "ipw", n.boot = 100
)# get bootstrap CI
$est.ci.mat boot_ipw_wri
Estimate Bootstrap CI low Bootstrap CI high
Positive -0.4361789 -1.332213 0.4598551
Negative -1.3953601 -2.492321 -0.2983993
HR ratio: strata over Positive
Positive 1.0000000
Negative 0.3832065
psm
method with dwc
model<- bootstrap_propen(
boot_ps_dwc data.in = biomarker, formula = indicator_2 ~ ECOG + Sex + Age,
indicator.var = "indicator_2",
tte = "OS", event = "OS.event", trt = "Arm",
class.of.int = list("Positive" = 1, "Negative" = 0),
estimate.res = ps_dwc, method = "ps", n.boot = 100
)# get bootstrap CI
$est.ci.mat boot_ps_dwc
Estimate Bootstrap CI low Bootstrap CI high
Positive -0.3724204 -1.715192 0.970351514
Negative -1.3035235 -2.600825 -0.006221976
HR ratio: strata over Positive
Positive 1.0000000
Negative 0.3941187
<- list(
boot_models ipw_plain = boot_ipw_plain,
ipw_wri = boot_ipw_wri,
ps_dwc = boot_ps_dwc
)<- c("Estimate", "Bootstrap CI low", "Bootstrap CI high")
cols # get HRs and bootstrap CIs
<- lapply(boot_models, function(x){
boots_dp <- x$est.ci.mat[ , cols, drop = FALSE] %>%
cis exp() %>% round(2) %>% as.data.frame()
colnames(cis) <- c("HR", "LOWER", "UPPER")
%>% mutate(
cis Group = rownames(cis)
)
})<- do.call(`rbind`, boots_dp) %>%
boots_dp mutate(
Methods = rep(c("IPW plain", "IPW wri", "PS dwc"), 2),
Methods_Group = paste(Methods, Group),
Group = factor(Group, levels = c("Positive", "Negative")),
Methods = factor(Methods, levels = c("IPW plain", "IPW wri", "PS dwc"))
%>%
) arrange(Methods, Group)
forest_bygroup(
data = boots_dp, summarystat = "HR", upperci = "UPPER", lowerci = "LOWER",
population.name = "Methods_Group", group.name = "Methods",
color.group.name = "Group",
stat.label = "Hazard Ratio",
stat.type.hr = TRUE, log.scale = FALSE,
endpoint.name = "OS", study.name = "Example Study", draw = TRUE
)
Below are the main functions in this package :
ipw_strata()
: IPW, allows for 2 or multiple classes. It also allows for predicting both patients with unknown class label in control arm and patients with unknown class label in trt arm. Supports both survival endpoint (HR as effect size estimator) and response endpoint (ORR difference as effect size estimator)
ps_match_strata()
: PSM, similar to ipw_strata()
km_plot_weight()
: weighted km, can take outputs from either ipw_strata()
or ps_match_strata()
bootstrap_propen()
: bootstrap function; can take outputs from both ipw_strata()
and ps_match_strata()
; outputs bootstrap CI as well as summary statistics (median, mean, etc) of bootstrap HRs (ORR deltas)
forest_bygroup()
: forestplot function; additional functionality to impose extra summary statistics on forestplot (e.g. median bootstrap HR on top of HR point estimate)
std_diff()
: Calculate feature’s standardized difference across two arms. This function can be used to check covariate balance between the (matched) strata from two arms.
expected_feature_diff()
: Calculate expected number of features showing standardized difference > threshold. This function can be used to assess whether the feature balance in the matched/adjusted data set is comparable to an RCT.