--- title: "Two-Stage Randomization for Competing risks and Survival outcomes" author: Klaus Holst & Thomas Scheike date: "2026-05-24" output: rmarkdown::html_vignette: fig_caption: no vignette: > %\VignetteIndexEntry{Two-Stage Randomization for Competing risks and Survival outcomes} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- Two-Stage Randomization for Competing risks and Survival outcomes ===================================================================== Under two-stage randomization we can estimate the average treatment effect $E(Y(i,\bar k))$ of treatment regime $(i,\bar k)$. - treatment $A_0=i$ and - for all responses, randomization $A_1 = (k_1)$, so treatment $k_1$ - response$\times A_1 = (k_1, k_2)$, so treatment $k_1$ if response 1, and treatment $k_2$ if response 2. The estimator can be augmented in different ways: using the two randomizations and the dynamic censoring augmentation. Estimating $\mu_{i,\bar k} = P(Y(i,\bar k,\epsilon=v) <= t)$, restricted mean $E( \min(Y(i,\bar k),\tau))$ or years lost $E( I(\epsilon=v) \cdot (\tau - \min(Y(i,\bar k),\tau)))$ using IPCW weighted estimating equations : \\ The solved estimating equation is \begin{align*} \sum_i \frac{I(min(T_i,t) < G_i)}{G_c(min(T_i ,t))} I(T \leq t, \epsilon=1 ) - AUG_0 - AUG_1 + AUG_C - p(i,j)) = 0 \end{align*} using the covariates from augmentR0 to augment with \begin{align*} AUG_0 = \frac{A_0(i) - \pi_0(i)}{ \pi_0(i)} X_0 \gamma_0 \end{align*} and using the covariates from augmentR1 to augment with \begin{align*} AUG_1 = \frac{A_0(i)}{\pi_0(i)} \frac{A_1(j) - \pi_1(j)}{ \pi_1(j)} X_1 \gamma_1 \end{align*} and censoring augmenting with \begin{align*} AUG_C = \int_0^t \gamma_c(s)^T (e(s) - \bar e(s)) \frac{1}{G_c(s) } dM_c(s) \end{align*} where $\gamma_c(s)$ is chosen to minimize the variance given the dynamic covariates specified by augmentC. - Treatments must be given as factors. - Treatment for the 2nd randomization may depend on response. - Treatment probabilities are estimated by default and uncertainty from this is adjusted for. - Randomization augmentation for the 1st and 2nd randomization is possible. - Censoring model can be stratified on observed covariates (at time 0). - Censoring augmentation is done dynamically over time with time-dependent covariates. Standard errors are estimated using the influence functions of all estimators; tests of differences can therefore be computed subsequently. Data must be in start-stop-status survival format with - one code of `status` indicating response, i.e. 2nd randomization - other codes defining the outcome of interest ``` r library(mets) set.seed(100) n <- 200 ddf <- mets:::gsim(n,covs=1,null=0,cens=1,ce=1,betac=c(0.3,1)) true <- apply(ddf$TTt<2,2,mean) true #> [1] 0.740 0.715 0.340 0.350 datat <- ddf$datat ## set-random response on data, only relevant after status==2 response <- rbinom(n,1,0.5) datat$response <- as.factor(response[datat$id]*datat$Count2) datat$A000 <- as.factor(1) datat$A111 <- as.factor(1) bb <- binregTSR(Event(entry,time,status)~+1+cluster(id),datat,time=2,cause=c(1),response.code=2, treat.model0=A0.f~+1, treat.model1=A1.f~A0.f, augmentR1=~X11+X12+TR, augmentR0=~X01+X02, augmentC=~X01+X02+A11t+A12t+X11+X12+TR, cens.model=~strata(A0.f)) bb #> Simple estimator : #> coef #> A0.f=1, response*A1.f=1 0.5723748 0.19721139 #> A0.f=1, response*A1.f=2 0.7354573 0.08755875 #> A0.f=2, response*A1.f=1 0.2834829 0.10028201 #> A0.f=2, response*A1.f=2 0.4790373 0.10007200 #> #> First Randomization Augmentation : #> coef #> A0.f=1, response*A1.f=1 0.5900338 0.21419118 #> A0.f=1, response*A1.f=2 0.7428143 0.08889887 #> A0.f=2, response*A1.f=1 0.2721992 0.10288280 #> A0.f=2, response*A1.f=2 0.4671886 0.10273188 #> #> Second Randomization Augmentation : #> coef #> A0.f=1, response*A1.f=1 0.5564110 0.24533756 #> A0.f=1, response*A1.f=2 0.7185240 0.09841671 #> A0.f=2, response*A1.f=1 0.2778309 0.09021056 #> A0.f=2, response*A1.f=2 0.4685339 0.09839184 #> #> 1st and 2nd Randomization Augmentation : #> coef #> A0.f=1, response*A1.f=1 0.5981472 0.26485855 #> A0.f=1, response*A1.f=2 0.7302633 0.10170090 #> A0.f=2, response*A1.f=1 0.2699247 0.09093482 #> A0.f=2, response*A1.f=2 0.4620365 0.09859532 estimate(coef=bb$riskG$riskG01[,1],vcov=crossprod(bb$riskG.iid$riskG01)) #> Estimate Std.Err 2.5% 97.5% P-value #> A0.f=1, response*A1.f=1 0.5981 0.26486 0.07903 1.1173 2.392e-02 #> A0.f=1, response*A1.f=2 0.7303 0.10170 0.53093 0.9296 6.946e-13 #> A0.f=2, response*A1.f=1 0.2699 0.09093 0.09170 0.4482 2.994e-03 #> A0.f=2, response*A1.f=2 0.4620 0.09860 0.26879 0.6553 2.783e-06 estimate(coef=bb$riskG$riskG01[,1],vcov=crossprod(bb$riskG.iid$riskG01),f=function(p) c(p[1]/p[2],p[3]/p[4])) #> Estimate Std.Err 2.5% 97.5% P-value #> A0.f=1, response*A1.f=1 0.8191 0.3600 0.1135 1.525 0.022884 #> A0.f=2, response*A1.f=1 0.5842 0.2249 0.1433 1.025 0.009399 estimate(coef=bb$riskG$riskG01[,1],vcov=crossprod(bb$riskG.iid$riskG01),f=function(p) c(p[1]-p[2],p[3]-p[4])) #> Estimate Std.Err 2.5% 97.5% P-value #> A0.f=1, response*A1.f=1 -0.1321 0.266 -0.6534 0.38920 0.6194 #> A0.f=2, response*A1.f=1 -0.1921 0.129 -0.4450 0.06076 0.1365 ``` ``` r ## 2 levels for each response , fixed weights datat$response.f <- as.factor(datat$response) bb <- binregTSR(Event(entry,time,status)~+1+cluster(id),datat,time=2,cause=c(1),response.code=2, treat.model0=A0.f~+1, treat.model1=A1.f~A0.f*response.f, augmentR0=~X01+X02, augmentR1=~X11+X12, augmentC=~X01+X02+A11t+A12t+X11+X12+TR, cens.model=~strata(A0.f), estpr=c(0,0),pi0=0.5,pi1=0.5) bb #> Simple estimator : #> coef #> A0.f=1, response.f*A1.f=1,1 0.5195752 0.18171540 #> A0.f=1, response.f*A1.f=2,1 0.5264182 0.18212146 #> A0.f=1, response.f*A1.f=1,2 0.6343225 0.09083112 #> A0.f=1, response.f*A1.f=2,2 0.6411655 0.08810024 #> A0.f=2, response.f*A1.f=1,1 0.3023142 0.11081204 #> A0.f=2, response.f*A1.f=2,1 0.3199096 0.09719853 #> A0.f=2, response.f*A1.f=1,2 0.5294417 0.13399244 #> A0.f=2, response.f*A1.f=2,2 0.5470372 0.13251846 #> #> First Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.f=1,1 0.6359737 0.21658343 #> A0.f=1, response.f*A1.f=2,1 0.6425345 0.21728891 #> A0.f=1, response.f*A1.f=1,2 0.7129831 0.06901176 #> A0.f=1, response.f*A1.f=2,2 0.7195438 0.06513763 #> A0.f=2, response.f*A1.f=1,1 0.2561355 0.12249210 #> A0.f=2, response.f*A1.f=2,1 0.2790365 0.10488801 #> A0.f=2, response.f*A1.f=1,2 0.4557345 0.14361723 #> A0.f=2, response.f*A1.f=2,2 0.4786356 0.14226926 #> #> Second Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.f=1,1 0.4046498 0.19742912 #> A0.f=1, response.f*A1.f=2,1 0.5687969 0.23376677 #> A0.f=1, response.f*A1.f=1,2 0.5910450 0.09523643 #> A0.f=1, response.f*A1.f=2,2 0.6498579 0.08872953 #> A0.f=2, response.f*A1.f=1,1 0.3039991 0.10638060 #> A0.f=2, response.f*A1.f=2,1 0.3311109 0.09138384 #> A0.f=2, response.f*A1.f=1,2 0.5569428 0.11332978 #> A0.f=2, response.f*A1.f=2,2 0.5164913 0.12796725 #> #> 1st and 2nd Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.f=1,1 0.5212352 0.21098387 #> A0.f=1, response.f*A1.f=2,1 0.6994484 0.26830813 #> A0.f=1, response.f*A1.f=1,2 0.6692919 0.07617047 #> A0.f=1, response.f*A1.f=2,2 0.7328127 0.06615399 #> A0.f=2, response.f*A1.f=1,1 0.2586853 0.11531506 #> A0.f=2, response.f*A1.f=2,1 0.2935352 0.09702298 #> A0.f=2, response.f*A1.f=1,2 0.4809089 0.11457277 #> A0.f=2, response.f*A1.f=2,2 0.4602799 0.13349781 ## 2 levels for each response , estimated treat probabilities bb <- binregTSR(Event(entry,time,status)~+1+cluster(id),datat,time=2,cause=c(1),response.code=2, treat.model0=A0.f~+1, treat.model1=A1.f~A0.f*response.f, augmentR0=~X01+X02, augmentR1=~X11+X12, augmentC=~X01+X02+A11t+A12t+X11+X12+TR, cens.model=~strata(A0.f),estpr=c(1,1)) bb #> Simple estimator : #> coef #> A0.f=1, response.f*A1.f=1,1 0.5733301 0.25088070 #> A0.f=1, response.f*A1.f=2,1 0.6022300 0.25101611 #> A0.f=1, response.f*A1.f=1,2 0.6893786 0.07700325 #> A0.f=1, response.f*A1.f=2,2 0.7182785 0.08178121 #> A0.f=2, response.f*A1.f=1,1 0.2851197 0.09997901 #> A0.f=2, response.f*A1.f=2,1 0.2841367 0.07893846 #> A0.f=2, response.f*A1.f=1,2 0.4821020 0.10343782 #> A0.f=2, response.f*A1.f=2,2 0.4811190 0.10002228 #> #> First Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.f=1,1 0.5934059 0.27131225 #> A0.f=1, response.f*A1.f=2,1 0.6224367 0.27204244 #> A0.f=1, response.f*A1.f=1,2 0.6962560 0.07772829 #> A0.f=1, response.f*A1.f=2,2 0.7252868 0.08319005 #> A0.f=2, response.f*A1.f=1,1 0.2738115 0.10244557 #> A0.f=2, response.f*A1.f=2,1 0.2767634 0.08081569 #> A0.f=2, response.f*A1.f=1,2 0.4661741 0.10481333 #> A0.f=2, response.f*A1.f=2,2 0.4691260 0.10238718 #> #> Second Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.f=1,1 0.4372672 0.23706174 #> A0.f=1, response.f*A1.f=2,1 0.5854429 0.26329507 #> A0.f=1, response.f*A1.f=1,2 0.6685369 0.07730893 #> A0.f=1, response.f*A1.f=2,2 0.7320359 0.08145477 #> A0.f=2, response.f*A1.f=1,1 0.2745442 0.10426421 #> A0.f=2, response.f*A1.f=2,1 0.2988823 0.07578047 #> A0.f=2, response.f*A1.f=1,2 0.5019477 0.09984288 #> A0.f=2, response.f*A1.f=2,2 0.4745778 0.10218550 #> #> 1st and 2nd Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.f=1,1 0.4702327 0.25293879 #> A0.f=1, response.f*A1.f=2,1 0.6127080 0.28587900 #> A0.f=1, response.f*A1.f=1,2 0.6759149 0.07701506 #> A0.f=1, response.f*A1.f=2,2 0.7433289 0.08370622 #> A0.f=2, response.f*A1.f=1,1 0.2634297 0.10565890 #> A0.f=2, response.f*A1.f=2,1 0.2948997 0.07656016 #> A0.f=2, response.f*A1.f=1,2 0.4865367 0.09860797 #> A0.f=2, response.f*A1.f=2,2 0.4699980 0.10193373 ## 2 and 3 levels for each response , fixed weights datat$A1.23.f <- as.numeric(datat$A1.f) dtable(datat,~A1.23.f+response) #> #> response 0 1 #> A1.23.f #> 1 120 23 #> 2 120 25 datat <- dtransform(datat,A1.23.f=2+rbinom(nrow(datat),1,0.5), Count2==1 & A1.23.f==2 & response==0) dtable(datat,~A1.23.f+response) #> #> response 0 1 #> A1.23.f #> 1 120 23 #> 2 111 25 #> 3 9 0 datat$A1.23.f <- as.factor(datat$A1.23.f) dtable(datat,~A1.23.f+response|Count2==1) #> #> response 0 1 #> A1.23.f #> 1 21 23 #> 2 10 25 #> 3 9 0 ### bb <- binregTSR(Event(entry,time,status)~+1+cluster(id),datat,time=2,cause=c(1),response.code=2, treat.model0=A0.f~+1,treat.model1=A1.23.f~A0.f*response.f, augmentR0=~X01+X02, augmentR1=~X11+X12, augmentC=~X01+X02+A11t+A12t+X11+X12+TR, cens.model=~strata(A0.f), estpr=c(1,0),pi1=c(0.3,0.5)) bb #> Simple estimator : #> coef #> A0.f=1, response.f*A1.23.f=1,1 0.5928568 0.20464868 #> A0.f=1, response.f*A1.23.f=2,1 0.6055291 0.20241441 #> A0.f=1, response.f*A1.23.f=3,1 0.5539792 0.19865003 #> A0.f=1, response.f*A1.23.f=1,2 0.7203538 0.09728960 #> A0.f=1, response.f*A1.23.f=2,2 0.7330260 0.08618558 #> A0.f=1, response.f*A1.23.f=3,2 0.6814762 0.07905734 #> A0.f=2, response.f*A1.23.f=1,1 0.3487105 0.14756981 #> A0.f=2, response.f*A1.23.f=2,1 0.2724117 0.09182351 #> A0.f=2, response.f*A1.23.f=3,1 0.2669705 0.10319223 #> A0.f=2, response.f*A1.23.f=1,2 0.5551901 0.15399346 #> A0.f=2, response.f*A1.23.f=2,2 0.4788913 0.12263546 #> A0.f=2, response.f*A1.23.f=3,2 0.4734501 0.12853150 #> #> First Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.23.f=1,1 0.6091794 0.22507560 #> A0.f=1, response.f*A1.23.f=2,1 0.6246762 0.22364858 #> A0.f=1, response.f*A1.23.f=3,1 0.5756603 0.21982886 #> A0.f=1, response.f*A1.23.f=1,2 0.7242102 0.10066515 #> A0.f=1, response.f*A1.23.f=2,2 0.7397070 0.08813234 #> A0.f=1, response.f*A1.23.f=3,2 0.6906911 0.07986256 #> A0.f=2, response.f*A1.23.f=1,1 0.3348012 0.15239160 #> A0.f=2, response.f*A1.23.f=2,1 0.2718558 0.09060009 #> A0.f=2, response.f*A1.23.f=3,1 0.2532931 0.10673462 #> A0.f=2, response.f*A1.23.f=1,2 0.5362379 0.15754640 #> A0.f=2, response.f*A1.23.f=2,2 0.4732925 0.12378699 #> A0.f=2, response.f*A1.23.f=3,2 0.4547298 0.13192021 #> #> Second Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.23.f=1,1 0.4707095 0.21804098 #> A0.f=1, response.f*A1.23.f=2,1 0.6287053 0.22907037 #> A0.f=1, response.f*A1.23.f=3,1 0.6778953 0.30311581 #> A0.f=1, response.f*A1.23.f=1,2 0.6499002 0.12488660 #> A0.f=1, response.f*A1.23.f=2,2 0.7656237 0.07140976 #> A0.f=1, response.f*A1.23.f=3,2 0.6875010 0.08312068 #> A0.f=2, response.f*A1.23.f=1,1 0.2656521 0.17396985 #> A0.f=2, response.f*A1.23.f=2,1 0.2706333 0.08715854 #> A0.f=2, response.f*A1.23.f=3,1 0.3235400 0.07994612 #> A0.f=2, response.f*A1.23.f=1,2 0.4759907 0.16071698 #> A0.f=2, response.f*A1.23.f=2,2 0.4671860 0.11714753 #> A0.f=2, response.f*A1.23.f=3,2 0.5050170 0.11067832 #> #> 1st and 2nd Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.23.f=1,1 0.5068915 0.22640669 #> A0.f=1, response.f*A1.23.f=2,1 0.6601925 0.25328437 #> A0.f=1, response.f*A1.23.f=3,1 0.7128045 0.32561624 #> A0.f=1, response.f*A1.23.f=1,2 0.6681178 0.12092323 #> A0.f=1, response.f*A1.23.f=2,2 0.7709766 0.07310528 #> A0.f=1, response.f*A1.23.f=3,2 0.7088981 0.08171051 #> A0.f=2, response.f*A1.23.f=1,1 0.2472033 0.17493179 #> A0.f=2, response.f*A1.23.f=2,1 0.2703258 0.08650158 #> A0.f=2, response.f*A1.23.f=3,1 0.3183141 0.08043082 #> A0.f=2, response.f*A1.23.f=1,2 0.4540746 0.15704473 #> A0.f=2, response.f*A1.23.f=2,2 0.4667557 0.11617870 #> A0.f=2, response.f*A1.23.f=3,2 0.4960792 0.11127903 ## 2 and 3 levels for each response , estimated bb <- binregTSR(Event(entry,time,status)~+1+cluster(id),datat,time=2,cause=c(1),response.code=2, treat.model0=A0.f~+1, treat.model1=A1.23.f~A0.f*response.f, augmentR0=~X01+X02, augmentR1=~X11+X12, augmentC=~X01+X02+A11t+A12t+X11+X12+TR,cens.model=~strata(A0.f), estpr=c(1,1)) bb #> Simple estimator : #> coef #> A0.f=1, response.f*A1.23.f=1,1 0.5733301 0.25088123 #> A0.f=1, response.f*A1.23.f=2,1 0.6486249 0.25863637 #> A0.f=1, response.f*A1.23.f=3,1 0.5558352 0.25058134 #> A0.f=1, response.f*A1.23.f=1,2 0.6893788 0.07700333 #> A0.f=1, response.f*A1.23.f=2,2 0.7646736 0.11893993 #> A0.f=1, response.f*A1.23.f=3,2 0.6718839 0.07549715 #> A0.f=2, response.f*A1.23.f=1,1 0.2851198 0.09997906 #> A0.f=2, response.f*A1.23.f=2,1 0.2795959 0.10086233 #> A0.f=2, response.f*A1.23.f=3,1 0.2893264 0.12751179 #> A0.f=2, response.f*A1.23.f=1,2 0.4821024 0.10343788 #> A0.f=2, response.f*A1.23.f=2,2 0.4765785 0.12129977 #> A0.f=2, response.f*A1.23.f=3,2 0.4863090 0.13928327 #> #> First Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.23.f=1,1 0.5934060 0.27131280 #> A0.f=1, response.f*A1.23.f=2,1 0.6665511 0.27965619 #> A0.f=1, response.f*A1.23.f=3,1 0.5783224 0.27132992 #> A0.f=1, response.f*A1.23.f=1,2 0.6962562 0.07772837 #> A0.f=1, response.f*A1.23.f=2,2 0.7694013 0.12078013 #> A0.f=1, response.f*A1.23.f=3,2 0.6811727 0.07593118 #> A0.f=2, response.f*A1.23.f=1,1 0.2738116 0.10244561 #> A0.f=2, response.f*A1.23.f=2,1 0.2791785 0.10064910 #> A0.f=2, response.f*A1.23.f=3,1 0.2740035 0.13235513 #> A0.f=2, response.f*A1.23.f=1,2 0.4661745 0.10481340 #> A0.f=2, response.f*A1.23.f=2,2 0.4715413 0.12263585 #> A0.f=2, response.f*A1.23.f=3,2 0.4663664 0.14396909 #> #> Second Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.23.f=1,1 0.4372664 0.23706211 #> A0.f=1, response.f*A1.23.f=2,1 0.5867572 0.24940274 #> A0.f=1, response.f*A1.23.f=3,1 0.6818858 0.34853309 #> A0.f=1, response.f*A1.23.f=1,2 0.6685368 0.07730905 #> A0.f=1, response.f*A1.23.f=2,2 0.7678408 0.07768404 #> A0.f=1, response.f*A1.23.f=3,2 0.6794316 0.07967452 #> A0.f=2, response.f*A1.23.f=1,1 0.2745441 0.10426433 #> A0.f=2, response.f*A1.23.f=2,1 0.2695053 0.09394894 #> A0.f=2, response.f*A1.23.f=3,1 0.3345432 0.09931154 #> A0.f=2, response.f*A1.23.f=1,2 0.5019478 0.09984303 #> A0.f=2, response.f*A1.23.f=2,2 0.4603587 0.11962056 #> A0.f=2, response.f*A1.23.f=3,2 0.5047271 0.12508846 #> #> 1st and 2nd Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.23.f=1,1 0.4702319 0.25293915 #> A0.f=1, response.f*A1.23.f=2,1 0.6117617 0.26980849 #> A0.f=1, response.f*A1.23.f=3,1 0.7132649 0.36786766 #> A0.f=1, response.f*A1.23.f=1,2 0.6759149 0.07701519 #> A0.f=1, response.f*A1.23.f=2,2 0.7733904 0.07957475 #> A0.f=1, response.f*A1.23.f=3,2 0.7066309 0.07677686 #> A0.f=2, response.f*A1.23.f=1,1 0.2634296 0.10565903 #> A0.f=2, response.f*A1.23.f=2,1 0.2691166 0.09330578 #> A0.f=2, response.f*A1.23.f=3,1 0.3296791 0.10056744 #> A0.f=2, response.f*A1.23.f=1,2 0.4865368 0.09860813 #> A0.f=2, response.f*A1.23.f=2,2 0.4605674 0.11777635 #> A0.f=2, response.f*A1.23.f=3,2 0.4971161 0.12478451 ## 2 and 1 level for each response datat$A1.21.f <- as.numeric(datat$A1.f) dtable(datat,~A1.21.f+response|Count2==1) #> #> response 0 1 #> A1.21.f #> 1 21 23 #> 2 19 25 datat <- dtransform(datat,A1.21.f=1,Count2==1 & response==1) dtable(datat,~A1.21.f+response|Count2==1) #> #> response 0 1 #> A1.21.f #> 1 21 48 #> 2 19 0 datat$A1.21.f <- as.factor(datat$A1.21.f) dtable(datat,~A1.21.f+response|Count2==1) #> #> response 0 1 #> A1.21.f #> 1 21 48 #> 2 19 0 bb <- binregTSR(Event(entry,time,status)~+1+cluster(id),datat,time=2,cause=c(1),response.code=2, treat.model0=A0.f~+1, treat.model1=A1.21.f~A0.f*response.f, augmentR0=~X01+X02, augmentR1=~X11+X12, augmentC=~X01+X02+A11t+A12t+X11+X12+TR,cens.model=~strata(A0.f), estpr=c(1,1)) bb #> Simple estimator : #> coef #> A0.f=1, response.f*A1.21.f=1,1 0.6352226 0.11813045 #> A0.f=1, response.f*A1.21.f=2,1 0.6641226 0.12002926 #> A0.f=2, response.f*A1.21.f=1,1 0.3865954 0.09118205 #> A0.f=2, response.f*A1.21.f=2,1 0.3856124 0.07700508 #> #> First Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.21.f=1,1 0.6482593 0.12794220 #> A0.f=1, response.f*A1.21.f=2,1 0.6772901 0.13063539 #> A0.f=2, response.f*A1.21.f=1,1 0.3729074 0.09195877 #> A0.f=2, response.f*A1.21.f=2,1 0.3758593 0.07808807 #> #> Second Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.21.f=1,1 0.6435175 0.11802262 #> A0.f=1, response.f*A1.21.f=2,1 0.6882148 0.11332360 #> A0.f=2, response.f*A1.21.f=1,1 0.3954927 0.09148962 #> A0.f=2, response.f*A1.21.f=2,1 0.3782504 0.08450470 #> #> 1st and 2nd Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.21.f=1,1 0.6685224 0.13033444 #> A0.f=1, response.f*A1.21.f=2,1 0.7138566 0.12655030 #> A0.f=2, response.f*A1.21.f=1,1 0.3839378 0.08986106 #> A0.f=2, response.f*A1.21.f=2,1 0.3758166 0.08434474 ## known weights bb <- binregTSR(Event(entry,time,status)~+1+cluster(id),datat,time=2,cause=c(1),response.code=2, treat.model0=A0.f~+1, treat.model1=A1.21.f~A0.f*response.f, augmentR0=~X01+X02, augmentR1=~X11+X12, augmentC=~X01+X02+A11t+A12t+X11+X12+TR, cens.model=~strata(A0.f),estpr=c(1,0),pi1=c(0.5,1)) bb #> Simple estimator : #> coef #> A0.f=1, response.f*A1.21.f=1,1 0.6410543 0.12005600 #> A0.f=1, response.f*A1.21.f=2,1 0.6486576 0.11864403 #> A0.f=2, response.f*A1.21.f=1,1 0.3780709 0.09171015 #> A0.f=2, response.f*A1.21.f=2,1 0.3940667 0.08377301 #> #> First Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.21.f=1,1 0.6532872 0.12997608 #> A0.f=1, response.f*A1.21.f=2,1 0.6625853 0.12899925 #> A0.f=2, response.f*A1.21.f=1,1 0.3647406 0.09331039 #> A0.f=2, response.f*A1.21.f=2,1 0.3842369 0.08422015 #> #> Second Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.21.f=1,1 0.6435110 0.12122751 #> A0.f=1, response.f*A1.21.f=2,1 0.6751763 0.11432293 #> A0.f=2, response.f*A1.21.f=1,1 0.3957743 0.08413935 #> A0.f=2, response.f*A1.21.f=2,1 0.3772145 0.09110553 #> #> 1st and 2nd Randomization Augmentation : #> coef #> A0.f=1, response.f*A1.21.f=1,1 0.6709750 0.13165583 #> A0.f=1, response.f*A1.21.f=2,1 0.7043078 0.12732297 #> A0.f=2, response.f*A1.21.f=1,1 0.3861201 0.08279377 #> A0.f=2, response.f*A1.21.f=2,1 0.3757694 0.09080739 ``` Two-Stage Randomization CALGB-9823 ================================== We illustrate an analysis of one SMART conducted by Cancer and Leukemia Group B Protocol 8923 (CALGB 8923), Stone et al. (2001). 388 patients were randomized to an initial treatment of GM-CSF ($A_1$) or standard chemotherapy ($A_2$). Patients with complete remission and informed consent to the second stage were then re-randomized to cytarabine only ($B_1$) or cytarabine plus mitoxantrone ($B_2$). We first compute the weighted risk-set estimator based on estimated weights \begin{align*} \Lambda_{A1,B1}(t) & = \sum_i \int_0^t \frac{w_i(s)}{Y^w(s)} dN_i(s) \end{align*} where $w_i(s) = I(A0_i=A1) + I(t>T_R) I(A1_i=B1)/\pi_1(X_i)$, that is 1 when starting on treatment $A1$, and then scaled up by the proportion switching to $B1$ at time $T_R$ for those that do so. This is equivalent to the IPTW (inverse probability of treatment weighted) estimator. We estimate the treatment regimes $A1, B1$ and $A2, B1$ by letting $A10$ indicate those consistent with ending on $B1$. $A10$ starts at 1 and becomes 0 if the subject is treated with $B2$, but stays 1 if treated with $B1$. We can then look at the two strata where $A0=0, A10=1$ and $A0=1, A10=1$. Similarly, we define $A11$ to start at 1, remain 1 if $B2$ is taken, and become 0 if the second randomization is $B1$. - the treatment models apply to all time-points, unless the `weight.var` variable is given (1 for treatments, 0 otherwise) to accommodate a general start-stop format - the treatment model may also depend on a response value - standard errors are based on influence functions and are also computed for the baseline We here use the propensity score model $P(A1=B1|A0)$ that uses the observed frequencies on arm $B1$ among those starting out on either $A1$ or $A2$. ``` r data(calgb8923) calgt <- calgb8923 tm=At.f~factor(Count2)+age+sex+wbc tm=At.f~factor(Count2) tm=At.f~factor(Count2)*A0.f head(calgt) #> id V X Z TR R U delta stop age wbc sex race time status start #> 1 1 0 0 0 0.00 0 13.33 1 13.33 64 128.0 1 1 13.338219 1 0.00 #> 2 2 1 1 0 0.00 0 17.80 1 17.80 71 4.3 2 1 17.802995 1 0.00 #> 3 3 1 0 0 0.00 0 1.27 1 1.27 71 43.6 2 1 1.271527 1 0.00 #> 4 4 1 0 1 0.00 0 24.77 1 24.77 63 72.3 2 1 0.730000 2 0.00 #> 5 4 1 0 1 0.73 1 24.77 1 24.77 63 72.3 2 1 24.772515 1 0.73 #> 6 5 0 1 0 0.00 0 10.37 1 10.37 65 1.4 1 1 10.374479 1 0.00 #> A0.f A0 A1 A11 A12 A1.f A10 At.f lbnr__id Count1 Count2 consent trt2 trt1 #> 1 0 0 0 1 0 0 0 0 1 0 0 -1 -1 1 #> 2 1 1 0 1 0 0 0 1 1 0 0 -1 -1 2 #> 3 0 0 0 1 0 0 0 0 1 0 0 -1 -1 1 #> 4 0 0 0 1 0 0 0 0 1 0 0 -1 -1 1 #> 5 0 0 1 1 1 1 1 1 2 0 1 1 1 1 #> 6 1 1 0 1 0 0 0 1 1 0 0 -1 -1 2 ll0 <- phreg_IPTW(Event(start,time,status==1)~strata(A0,A10)+cluster(id),calgt,treat.model=tm) pll0 <- predict(ll0,expand.grid(A0=0:1,A10=0,id=1)) ll1 <- phreg_IPTW(Event(start,time,status==1)~strata(A0,A11)+cluster(id),calgt,treat.model=tm) pll1 <- predict(ll1,expand.grid(A0=0:1,A11=1,id=1)) plot(pll0,se=1,lwd=2,col=1:2,lty=1,xlab="time (months)",xlim=c(0,30)) plot(pll1,add=TRUE,col=3:4,se=1,lwd=2,lty=1,xlim=c(0,30)) abline(h=0.25) legend("topright",c("A1B1","A2B1","A1B2","A2B2"),col=c(1,2,3,4),lty=1) ``` Survival estimates for four treatment regimes. ``` r summary(pll1,times=12) #> Predictions of type 'surv' #> Showing subjects: 1, 2 #> Showing times: 12 #> #> -- Subject 1 -- #> time surv se lower upper #> 12 0.4557 0.0421 0.3801 0.5462 #> #> -- Subject 2 -- #> time surv se lower upper #> 12 0.5029 0.0417 0.4276 0.5916 summary(pll0,times=12) #> Predictions of type 'surv' #> Showing subjects: 1, 2 #> Showing times: 12 #> #> -- Subject 1 -- #> time surv se lower upper #> 12 0.395 0.0426 0.3197 0.4881 #> #> -- Subject 2 -- #> time surv se lower upper #> 12 0.4279 0.0434 0.3507 0.5221 ``` The propensity score model can be extended to use covariates to get increased efficiency. Note also that the propensity scores for $A0$ will cancel out in the different strata. SessionInfo ============ ``` r sessionInfo() #> R version 4.6.0 (2026-04-24) #> Platform: x86_64-pc-linux-gnu #> Running under: Ubuntu 24.04.4 LTS #> #> Matrix products: default #> BLAS: /home/kkzh/.asdf/installs/r/4.6.0/lib/R/lib/libRblas.so #> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0 #> #> locale: #> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C #> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 #> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 #> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C #> [9] LC_ADDRESS=C LC_TELEPHONE=C #> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C #> #> time zone: Europe/Copenhagen #> tzcode source: system (glibc) #> #> attached base packages: #> [1] stats graphics grDevices utils datasets methods base #> #> other attached packages: #> [1] timereg_2.0.7 survival_3.8-6 mets_1.3.10 #> #> loaded via a namespace (and not attached): #> [1] cli_3.6.6 knitr_1.51 rlang_1.2.0 #> [4] xfun_0.57 otel_0.2.0 future.apply_1.20.2 #> [7] listenv_0.10.1 lava_1.9.1 stats4_4.6.0 #> [10] grid_4.6.0 evaluate_1.0.5 yaml_2.3.12 #> [13] mvtnorm_1.3-7 numDeriv_2016.8-1.1 compiler_4.6.0 #> [16] codetools_0.2-20 Rcpp_1.1.1-1.1 ucminf_1.2.3 #> [19] future_1.70.0 lattice_0.22-9 digest_0.6.39 #> [22] parallelly_1.47.0 parallel_4.6.0 splines_4.6.0 #> [25] Matrix_1.7-5 tools_4.6.0 RcppArmadillo_15.2.6-1 #> [28] globals_0.19.1 ```