Once weights have been obtained for each event time in the data set, we fit a (potentially stratified) Cox proportional hazards model to this weighted data set to obtain an estimate of the hazard ratio. The confidence interval for the hazard ratio can be derived by bootstrapping the weight estimation and the subsequent model-fitting process.
Example for Pooled Logistic Regression Switching ModelFirst we prepare the data.
sim1 <- tsegestsim(
n = 500, allocation1 = 2, allocation2 = 1, pbprog = 0.5,
trtlghr = -0.5, bprogsl = 0.3, shape1 = 1.8,
scale1 = 360, shape2 = 1.7, scale2 = 688,
pmix = 0.5, admin = 5000, pcatnotrtbprog = 0.5,
pcattrtbprog = 0.25, pcatnotrt = 0.2, pcattrt = 0.1,
catmult = 0.5, tdxo = 1, ppoor = 0.1, pgood = 0.04,
ppoormet = 0.4, pgoodmet = 0.2, xomult = 1.4188308,
milestone = 546, outputRawDataset = 1, seed = 2000)Given that the observations are collected at regular intervals, we can use a pooled logistic regression switching model. In this model, the prognostic baseline variable is represented by bprog, while catlag and the interaction between bprog and catlag serve as time-dependent covariates. Additionally, ns1, ns2, and ns3 are the three terms for the natural cubic spline with \(\nu=3\) degrees of freedom for visit-specific intercepts.
fit1 <- ipcw(
sim1$paneldata, id = "id", tstart = "tstart",
tstop = "tstop", event = "died", treat = "trtrand",
swtrt = "xo", swtrt_time = "xotime", base_cov = "bprog",
numerator = "bprog", denominator = "bprog*catlag",
logistic_switching_model = TRUE, ns_df = 3,
swtrt_control_only = TRUE, boot = FALSE)The fits for the denominator and numerator switching models are as follows, utilizing robust sandwich variance estimates to account for correlations among observations within the same subject.
# denominator switching model fit
fit1$fit_switch[[1]]$fit_den$parest[, c("param", "beta", "sebeta", "z")]
#> param beta sebeta z
#> 1 (Intercept) -6.37314396 0.5980953 -10.65573240
#> 2 bprog 1.73154601 0.3131088 5.53017307
#> 3 catlag 1.02546588 0.5368841 1.91003221
#> 4 bprog.catlag -0.02214233 0.6159246 -0.03594975
#> 5 ns1 0.48314846 0.5906830 0.81794876
#> 6 ns2 3.80068738 1.1802393 3.22026836
#> 7 ns3 -0.50255252 0.7789188 -0.64519242
# numerator switching model fit
fit1$fit_switch[[1]]$fit_num$parest[, c("param", "beta", "sebeta", "z")]
#> param beta sebeta z
#> 1 (Intercept) -6.43591583 0.5776690 -11.14118224
#> 2 bprog 1.82417620 0.2723647 6.69755047
#> 3 ns1 1.06152507 0.5469071 1.94096064
#> 4 ns2 4.33810022 1.1474070 3.78078589
#> 5 ns3 -0.04371903 0.7460213 -0.05860293Since treatment switching is not allowed in the experimental arm, the weights are identical to 1 for subjects in the experimental group. The unstabilized and stablized weights for the control group are plotted below.
# unstabilized weights
ggplot(fit1$data_outcome %>% filter(trtrand == 0),
aes(x = unstabilized_weight)) +
geom_density(fill="#77bd89", color="#1f6e34", alpha=0.8) +
scale_x_continuous("unstabilized weights")
# stabilized weights
ggplot(fit1$data_outcome %>% filter(trtrand == 0),
aes(x = stabilized_weight)) +
geom_density(fill="#77bd89", color="#1f6e34", alpha=0.8) +
scale_x_continuous("stabilized weights")Now we fit a weighted outcome Cox model and compare the treatemnt hazard ratio estimate with the reported.
fit1$fit_outcome$parest[, c("param", "beta", "sebeta", "z")]
#> param beta sebeta z
#> 1 treated 0.13954433 0.015115183 9.232063
#> 2 bprog -0.03834935 0.007769929 -4.935612
exp(fit1$fit_outcome$parest[1, c("beta", "lower", "upper")])
#> beta lower upper
#> 1 1.14975 1.116188 1.184321Now we apply the IPCW method using a Cox proportional hazards model with time-dependent covariates as the switching model.
fit2 <- ipcw(
shilong, id = "id", tstart = "tstart", tstop = "tstop",
event = "event", treat = "bras.f", swtrt = "co",
swtrt_time = "dco",
base_cov = c("agerand", "sex.f", "tt_Lnum", "rmh_alea.c",
"pathway.f"),
numerator = c("agerand", "sex.f", "tt_Lnum", "rmh_alea.c",
"pathway.f"),
denominator = c("agerand", "sex.f", "tt_Lnum", "rmh_alea.c",
"pathway.f", "ps", "ttc", "tran"),
swtrt_control_only = FALSE, boot = FALSE)The fits for the denominator and numerator switching models for the control arm are as follows, utilizing robust sandwich variance estimates to account for correlations among observations within the same subject.
# denominator switching model for the control group
fit2$fit_switch[[1]]$fit_den$parest[, c("param", "beta", "sebeta", "z")]
#> param beta sebeta z
#> 1 agerand 0.007642073 0.01009207 0.7572351
#> 2 sex.fFemale -0.364211003 0.26324814 -1.3835273
#> 3 tt_Lnum 0.042406215 0.04639827 0.9139611
#> 4 rmh_alea.c -0.351616244 0.27378303 -1.2842880
#> 5 pathway.fHR -0.041768569 0.39817164 -0.1049009
#> 6 pathway.fPI3K.AKT.mTOR 0.267055320 0.43180962 0.6184562
#> 7 ps 0.103478553 0.17788020 0.5817317
#> 8 ttc -0.480378314 0.32124349 -1.4953713
#> 9 tran 0.295828936 0.36809457 0.8036765
# numerator switching model for the control group
fit2$fit_switch[[1]]$fit_num$parest[, c("param", "beta", "sebeta", "z")]
#> param beta sebeta z
#> 1 agerand 0.01059477 0.00957653 1.10632644
#> 2 sex.fFemale -0.32452112 0.25680716 -1.26367627
#> 3 tt_Lnum 0.05956492 0.04488927 1.32693011
#> 4 rmh_alea.c -0.29758402 0.25884651 -1.14965436
#> 5 pathway.fHR 0.02467202 0.39869959 0.06188123
#> 6 pathway.fPI3K.AKT.mTOR 0.28245603 0.42341002 0.66709812The fits for the denominator and numerator switching models for the experimental arm are as follows:
# denominator switching model for the experimental group
fit2$fit_switch[[2]]$fit_den$parest[, c("param", "beta", "sebeta", "z")]
#> param beta sebeta z
#> 1 agerand -0.002639441 0.01482504 -0.1780394
#> 2 sex.fFemale 0.428028469 0.49814996 0.8592362
#> 3 tt_Lnum -0.152758042 0.09460056 -1.6147689
#> 4 rmh_alea.c 0.210365664 0.53818671 0.3908786
#> 5 pathway.fHR 1.774466817 1.07337446 1.6531666
#> 6 pathway.fPI3K.AKT.mTOR 0.859950234 1.06150371 0.8101246
#> 7 ps 0.261142698 0.28658566 0.9112204
#> 8 ttc -0.408175689 0.57303595 -0.7123038
#> 9 tran 1.053347294 0.50568571 2.0830078
# numerator switching model for the experimental group
fit2$fit_switch[[2]]$fit_num$parest[, c("param", "beta", "sebeta", "z")]
#> param beta sebeta z
#> 1 agerand 0.001625622 0.01543679 0.1053083
#> 2 sex.fFemale 0.476979559 0.47723504 0.9994647
#> 3 tt_Lnum -0.160020273 0.09602402 -1.6664609
#> 4 rmh_alea.c 0.154833260 0.49850922 0.3105926
#> 5 pathway.fHR 1.841535014 1.06349774 1.7315834
#> 6 pathway.fPI3K.AKT.mTOR 1.064637540 1.05309041 1.0109650Below, the unstabilized and stabilized weights are plotted by treatment group: the left panel displays data for the control group, while the right panel shows data for the experimental group.
# unstabilized weights
ggplot(fit2$data_outcome, aes(x = unstabilized_weight)) +
geom_density(fill="#77bd89", color="#1f6e34", alpha=0.8) +
scale_x_continuous("unstabilized weights") +
facet_wrap(~treated)
# stabilized weights
ggplot(fit2$data_outcome, aes(x = stabilized_weight)) +
geom_density(fill="#77bd89", color="#1f6e34", alpha=0.8) +
scale_x_continuous("stabilized weights") +
facet_wrap(~treated)Finally, we fit the weighted outcome Cox model and compare the treatment hazard ratio estimate with the reported.
fit2$fit_outcome$parest[, c("param", "beta", "sebeta", "z")]
#> param beta sebeta z
#> 1 treated 0.356390611 0.25526832 1.3961412
#> 2 agerand -0.006047034 0.01018596 -0.5936636
#> 3 sex.fFemale -0.487409540 0.24458876 -1.9927716
#> 4 tt_Lnum 0.011244574 0.04147245 0.2711336
#> 5 rmh_alea.c 0.941651485 0.25315061 3.7197283
#> 6 pathway.fHR -0.127273307 0.35878557 -0.3547336
#> 7 pathway.fPI3K.AKT.mTOR -0.166035907 0.34997206 -0.4744262
c(fit2$hr, fit2$hr_CI)
#> [1] 1.4281653 0.8659517 2.3553924RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4