We consider two while-alive estimands for recurrent events data \[\begin{align*} \frac{E(N(D \wedge t))}{E(D \wedge t)} \end{align*}\] and the mean of the subject specific events per time-unit \[\begin{align*} E( \frac{N(D \wedge t)}{D \wedge t} ) \end{align*}\] for two treatment-groups in the case of an RCT.
data(hfaction_cpx12)
dtable(hfaction_cpx12,~status)
#>
#> status
#> 0 1 2
#> 617 1391 124
dd <- WA_recurrent(Event(entry,time,status)~treatment+cluster(id),hfaction_cpx12,time=2,
death.code=2)
summary(dd)
#> While-Alive summaries:
#>
#> RMST, E(min(D,t))
#> Estimate Std.Err 2.5% 97.5% P-value
#> treatment0 1.859 0.02108 1.817 1.900 0
#> treatment1 1.924 0.01502 1.894 1.953 0
#>
#> Estimate Std.Err 2.5% 97.5% P-value
#> [treatment0] - [treat.... -0.06517 0.02588 -0.1159 -0.01444 0.0118
#>
#> Null Hypothesis:
#> [treatment0] - [treatment1] = 0
#> mean events, E(N(min(D,t))):
#> Estimate Std.Err 2.5% 97.5% P-value
#> treatment0 1.572 0.09573 1.384 1.759 1.375e-60
#> treatment1 1.453 0.10315 1.251 1.656 4.376e-45
#>
#> Estimate Std.Err 2.5% 97.5% P-value
#> [treatment0] - [treat.... 0.1185 0.1407 -0.1574 0.3943 0.4
#>
#> Null Hypothesis:
#> [treatment0] - [treatment1] = 0
#> _______________________________________________________
#> Ratio of means E(N(min(D,t)))/E(min(D,t))
#> Estimate Std.Err 2.5% 97.5% P-value
#> treatment0 0.8457 0.05264 0.7425 0.9488 4.411e-58
#> treatment1 0.7555 0.05433 0.6490 0.8619 5.963e-44
#>
#> Estimate Std.Err 2.5% 97.5% P-value
#> [treatment0] - [treat.... 0.09022 0.07565 -0.05805 0.2385 0.233
#>
#> Null Hypothesis:
#> [treatment0] - [treatment1] = 0
#> _______________________________________________________
#> Mean of Events per time-unit E(N(min(D,t))/min(D,t))
#> Estimate Std.Err 2.5% 97.5% P-value
#> treat0 1.0725 0.1222 0.8331 1.3119 1.645e-18
#> treat1 0.7552 0.0643 0.6291 0.8812 7.508e-32
#>
#> Estimate Std.Err 2.5% 97.5% P-value
#> [treat0] - [treat1] 0.3173 0.1381 0.04675 0.5879 0.02153
#>
#> Null Hypothesis:
#> [treat0] - [treat1] = 0
WA_recurrent <- function(formula,data,time=NULL,cens.code=0,cause=1,death.code=2,
trans=NULL,cens.formula=NULL,augmentR=NULL,augmentC=NULL,type=NULL,...)
{ ## {{{
cl <- match.call() ## {{{
m <- match.call(expand.dots = TRUE)[1:3]
special <- c("strata", "cluster","offset")
Terms <- terms(formula, special, data = data)
m$formula <- Terms
m[[1]] <- as.name("model.frame")
m <- eval(m, parent.frame())
Y <- model.extract(m, "response")
if (!inherits(Y,"Event")) stop("Expected a 'Event'-object")
if (ncol(Y)==2) {
stop("must give start stop formulation \n");
exit <- Y[,1]
entry <- NULL ## rep(0,nrow(Y))
status <- Y[,2]
} else {
### stop("only right censored data, will not work for delayed entry\n");
entry <- Y[,1]
exit <- Y[,2]
status <- Y[,3]
}
id <- strata <- NULL
if (!is.null(attributes(Terms)$specials$cluster)) {
ts <- survival::untangle.specials(Terms, "cluster")
pos.cluster <- ts$terms
Terms <- Terms[-ts$terms]
id <- m[[ts$vars]]
} else pos.cluster <- NULL
if (!is.null(stratapos <- attributes(Terms)$specials$strata)) {
ts <- survival::untangle.specials(Terms, "strata")
pos.strata <- ts$terms
Terms <- Terms[-ts$terms]
strata <- m[[ts$vars]]
strata.name <- ts$vars
} else { strata.name <- NULL; pos.strata <- NULL}
if (!is.null(offsetpos <- attributes(Terms)$specials$offset)) {
ts <- survival::untangle.specials(Terms, "offset")
Terms <- Terms[-ts$terms]
offset <- m[[ts$vars]]
}
X <- model.matrix(Terms, m)
if (ncol(X)==0) X <- matrix(nrow=0,ncol=0)
### possible handling of id to code from 0:(antid-1)
if (!is.null(id)) {
ids <- sort(unique(id))
nid <- length(ids)
if (is.numeric(id)) id <- fast.approx(ids,id)-1 else {
id <- as.integer(factor(id,labels=seq(nid)))-1
}
orig.id <- id
} else { orig.id <- NULL; nid <- nrow(X);
id <- as.integer(seq_along(exit))-1; ids <- NULL
}
### id from call coded as numeric 1 ->
if (is.null(offset)) offset <- rep(0,length(exit))
if (is.null(weights)) weights <- rep(1,length(exit))
data$id__ <- id ## }}}
## use sorted id for all things
cid <- countID(data,"id__",sorted=TRUE)
### cid <- countID(data,"id__")
data$id__ <- cid$indexid
rrR <- subset(data,cid$reverseCountid==1)
## first var on rhs of formula
vars <- all.vars(formula)
treat.name <- vars[4]
treatvar <- data[,treat.name]
if (!is.factor(treatvar)) stop(paste("treatment=",treat.name," must be coded as factor \n",sep=""));
nlev <- nlevels(treatvar)
nlevs <- levels(treatvar)
ntreatvar <- as.numeric(treatvar)-1
treat.formula <- treat.model <- as.formula(paste(treat.name,"~+1",sep=""))
if (is.null(cens.formula)) cens.formula <- as.formula( paste("~strata(",treat.name,")",collapse=""))
formC <- as.formula( paste("Event(",vars[1],",",vars[2],",",vars[3],"==cens.code)~+cluster(id__)",collapse=""))
formD <- as.formula( paste("Event(",vars[2],",",vars[3],"!=cens.code)~-1+",treat.name,"+cluster(id__)",collapse=""))
form1 <- as.formula( paste("Event(",vars[2],",",vars[3],")~-1+",treat.name,"+cluster(id__)",collapse=""))
### varR <- all.vars(augmentR)
### newaugR <- paste(c(varR),sep="",collapse="+")
### varsR <- c(attr(terms(augmentR), "term.labels"))
## take out intercept
formrec <- update(formula, reformulate(c(".", "-1")))
if (!is.null(augmentR)) {
varsR <- c(attr(terms(augmentR), "term.labels"))
form1X <- update(form1, reformulate(c(".", varsR)))
} else form1X <- form1
###print(formD); print(formrec); print(form1); print(form1X); print(cens.formula);
## ratio of means ## {{{
dd <- resmeanIPCW(formD,data=rrR,cause=1,cens.model=cens.formula,time=time,model="l")
ddN <- recregIPCW(formrec,data=data,cause=cause,death.code=death.code,cens.model=cens.formula,times=time,model="l")
cc <- c(ddN$coef,dd$coef)
cciid <- cbind(ddN$iid,dd$iid)
ratio.means <- estimate(coef=cc,vcov=crossprod(cciid),f=function(p) (p[1:2]/p[3:4]))
ratio.means.log <- estimate(coef=cc,vcov=crossprod(cciid),f=function(p) log(p[1:2]/p[3:4]))
RAW <- list(iid=cciid,coef=cc,time=time,rmst=dd,meanN=ddN,ratio.means=ratio.means,ratio.means.log=ratio.means.log)
## }}}
data <- count.history(data,id="id__",types=cause,status=vars[3])
nameCount <- paste("Count",cause,sep="")
formulaCount <- update.formula(formula,.~+1)
cform <- as.formula(paste("~",nameCount,"+cluster(id__)",sep=""))
formulaCount <- update.formula(formulaCount,cform)
## While-Alive mean of events per time-unit
dataDmin <- evalTerminal(formulaCount,data=data,time=time,death.code=death.code)
rrR[,"ratio__"] <- dataDmin[cid$reverseCountid==1,"ratio"]
if (!is.null(trans)) {
rrR[,"ratio__"] <- dataDmin[cid$reverseCountid==1,4]^trans
}
Yr <- rrR[,"ratio__"]
if (is.null(type)) {
if (is.null(augmentC)) type <- "II" else type <- "I"
}
outae <- binregATE(form1X,rrR,cause=death.code,time=time,treat.model=treat.formula,
Ydirect=Yr,outcome="rmst",model="lin",cens.model=cens.formula,type=type[1],...)
ET <- list(riskDR=outae)
#ids <- countID(data,"id__",sorted=TRUE)
### assume id is ordered
data[,"ratio__"] <- outae$Yipcw[cid$indexid+1]
if (!is.null(augmentC)) { ## {{{
dc0 <- dynCensAug(formC,subset(data,ntreatvar==0),augmentC=augmentC,response="ratio__",time=time)
dc1 <- dynCensAug(formC,subset(data,ntreatvar==1),augmentC=augmentC,response="ratio__",time=time)
###nn <- table(rr$treat)
nid <- nrow(rrR)
treats <- function(treatvar) { treatvar <- droplevels(treatvar)# {{{
nlev <- nlevels(treatvar)
nlevs <- levels(treatvar)
ntreatvar <- as.numeric(treatvar)
return(list(nlev = nlev, nlevs = nlevs, ntreatvar = ntreatvar))
}
fittreat <- function(treat.model, data, id, ntreatvar, nlev) {
if (nlev == 2) {
treat.model <- drop.specials(treat.model, "cluster")
treat <- glm(treat.model, data, family = "binomial")
iidalpha <- lava::iid(treat, id =id)
lpa <- treat$linear.predictors
pal <- expit(lpa)
pal <- cbind(1 - pal, pal)
ppp <- (pal/pal[, 1])
spp <- 1/pal[, 1]
}
else {
treat.modelid <- update.formula(treat.model, . ~ . + cluster(id__))
treat <- mlogit(treat.modelid, data)
iidalpha <- lava::iid(treat)
pal <- predictmlogit(treat, data, se = 0, response = FALSE)
ppp <- (pal/pal[, 1])
spp <- 1/pal[, 1]
}
Xtreat <- model.matrix(treat.model, data)
tvg2 <- 1 * (ntreatvar >= 2)
pA <- c(mdi(pal, 1:length(ntreatvar), ntreatvar))
pppy <- c(mdi(ppp, 1:length(ntreatvar), ntreatvar))
Dppy <- (spp * tvg2 - pppy)
Dp <- c()
for (i in seq(nlev - 1)) Dp <- cbind(Dp, Xtreat * ppp[,
i + 1] * Dppy/spp^2)
DPai <- -1 * Dp/pA^2
out <- list(iidalpha = iidalpha, pA = pA, Dp = Dp, pal = pal,
ppp = ppp, spp = spp, id = id, DPai = DPai)
return(out)
} # }}}
expit <- function(x) 1/(1 + exp(-x))
idW <- rrR[,"id__"]
treatsvar <- rrR[,treat.name]
treats <- treats(treatsvar)
fitt <- fittreat(treat.model, rrR, idW, treats$ntreatvar, treats$nlev)
iidalpha0 <- fitt$iidalpha
wPA <- c(fitt$pA)
riskDRC <- outae$riskDR+c(dc0$augment,dc1$augment)/nid
MGC0 <- sumstrata(dc0$MGCiid/c(wPA[dc0$id+1]),dc0$id,nid)*dc0$n
MGC1 <- sumstrata(dc1$MGCiid/c(wPA[dc1$id+1]),dc1$id,nid)*dc1$n
MGC <- cbind(MGC0,MGC1)
ccaugment <- apply(MGC,2,sum)
riskDRC <- outae$riskDR+ccaugment/nid
iidDRC <- outae$riskDR.iid+MGC/nid
varDRC <- crossprod(iidDRC)
se.riskDRC <- diag(varDRC)^.5
riskDRC <- list(riskDRC=riskDRC,iid=iidDRC,var=varDRC,coef=riskDRC,se=se.riskDRC)
ET <- list(riskDRC=riskDRC,riskDR=outae)
} ## }}}
out <- list(time=time,id=cid$indexid,orig.id=orig.id,trans=trans,cause=cause,cens.code=cens.code,death.code,
RAW=RAW,ET=ET,augmentR=augmentR,augmentC=augmentC)
class(out) <- "WA"
return(out)
} ## }}}
##' @export
print.WA <- function(x,type="log",...) {# {{{
print(summary(x),type=type,...)
}# }}}
##' @export
summary.WA <- function(object,type="p",...) {# {{{
rmst <- estimate(object$RAW$rmst)
rmst.test <- estimate(rmst,contrast=rbind(c(1,-1)))
rmst.log <- estimate(rmst,function(p) log(p))
rmst.test.log <- estimate(rmst.log,contrast=rbind(c(1,-1)))
meanNtD <- estimate(object$RAW$meanN)
meanNtD.test <- estimate(meanNtD,contrast=rbind(c(1,-1)))
meanNtD.log <- estimate(meanNtD,function(p) log(p))
meanNtD.test.log <- estimate(meanNtD.log,contrast=rbind(c(1,-1)))
eer <- estimate(object$RAW$ratio.means)
eedr <- estimate(eer,contrast=rbind(c(1,-1)))
ee <- estimate(coef=object$ET$riskDR$riskDR,vcov=object$ET$riskDR$var.riskDR)
eed <- estimate(ee,contrast=rbind(c(1,-1)))
eer.log <- object$RAW$ratio.means.log
eedr.log <- estimate(eer.log,contrast=rbind(c(1,-1)))
eelog <- estimate(ee,function(p) log(p))
eedlog <- estimate(eelog,contrast=rbind(c(1,-1)))
res <- list(rmst=rmst,rmst.test=rmst.test,meanNtD=meanNtD,meanNtD.test=meanNtD.test,
ratio=eer,test.ratio=eedr,meanpt=ee,test.meanpt=eed)
reslog <- list(rmst=rmst.log,rmst.test=rmst.test.log,
meanNtD=meanNtD.log,meanNtD.test=meanNtD.test.log,
ratio=eer.log,test.ratio=eedr.log,meanpt=eelog,test.meanpt=eedlog)
class(res) <- "summary.WA"
class(reslog) <- "summary.WA"
attr(res,"log") <- (type!="p")
attr(reslog,"log") <- (type!="p")
if (type=="p") return(res) else return(reslog)
}# }}}
##' @export
print.summary.WA <- function(x,...)
{# {{{
if (attr(x,"log"))
cat(paste("While-Alive summaries, log-scale:",x$time,"\n\n"))
else
cat(paste("While-Alive summaries:",x$time,"\n\n"))
cat("RMST, E(min(D,t)) \n")
print(x$rmst)
cat(" \n")
print(x$rmst.test)
cat("mean events, E(N(min(D,t))): \n")
print(x$meanNtD)
cat(" \n")
print(x$meanNtD.test)
cat("_______________________________________________________ \n")
cat("Ratio of means E(N(min(D,t)))/E(min(D,t)) \n")
print(x$ratio)
cat(" \n")
print(x$test.ratio)
cat("_______________________________________________________ \n")
cat("Mean of Events per time-unit E(N(min(D,t))/min(D,t)) \n")
print(x$meanpt)
cat(" \n")
print(x$test.meanpt)
} # }}}
dynCensAug <- function(formC,data,augmentC=~+1,response="Yipcw",time=NULL,Z=NULL) { ## {{{
if (is.null(time)) stop("must give time of response \n")
data$Y__ <- data[,response]
varsC <- c("Y__",attr(terms(augmentC), "term.labels"))
formCC <- update(formC, reformulate(c(".", varsC)))
### www <- rep(1,nrow(data))
### if (treat.specific.cens==1) www <-data$WW1__;
cr2 <- phreg(formCC, data = data, no.opt = TRUE, no.var = 1,Z=Z)
xx <- cr2$cox.prep
icoxS0 <- rep(0,length(cr2$S0))
icoxS0[cr2$S0>0] <- 1/cr2$S0[cr2$S0>0]
S0i <- rep(0, length(xx$strata))
S0i[xx$jumps + 1] <- icoxS0
km <- TRUE
if (!km) {
cumhazD <- c(cumsumstratasum(S0i, xx$strata, xx$nstrata)$lagsum)
St <- exp(-cumhazD)
} else St <- c(exp(cumsumstratasum(log(1 - S0i), xx$strata, xx$nstrata)$lagsum))
nterms <- cr2$p-1
dhessian <- cr2$hessianttime
dhessian <- .Call("XXMatFULL",dhessian,cr2$p,PACKAGE="mets")$XXf
### matrix(apply(dhessian,2,sum),3,3)
timeb <- which(cr2$cumhaz[,1]<time)
### take relevant \sum H_i(s,t) (e_i - \bar e)
covts <- dhessian[timeb,1+1:nterms,drop=FALSE]
### construct relevant \sum (e_i - \bar e)^2
Pt <- dhessian[timeb,-c((1:(nterms+1)),(1:(nterms))*(nterms+1)+1),drop=FALSE]
### matrix(apply(dhessian[,c(5,6,8,9)],2,sum),2,2)
gammatt <- .Call("CubeVec",Pt,covts,1,PACKAGE="mets")$XXbeta
S0 <- cr2$S0[timeb]
gammatt[is.na(gammatt)] <- 0
gammatt[gammatt==Inf] <- 0
Gctb <- St[cr2$cox.prep$jumps+1][timeb]
augmentt.times <- apply(gammatt*cr2$U[timeb,1+1:nterms,drop=FALSE],1,sum)
augment.times <- sum(augmentt.times)
if (!is.null(Z)) {
Zj <- cr2$cox.prep$Z[cr2$cox.prep$jumps+1,][timeb]
Xaugment.times <- apply( augmentt.times*Zj,2,sum)
}
p <- 1
#### iid magic for censoring augmentation martingale #
### int_0^infty gamma(t) (e_i - ebar(s)) 1/G_c(s) dM_i^c
xx <- cr2$cox.prep
nid <- max(xx$id)+1
jumpsC <- xx$jumps+1
rr0 <- xx$sign
S0i <- rep(0,length(xx$strata))
S0i[jumpsC] <- c(1/(icoxS0*St[jumpsC]))
S0i[jumpsC] <- icoxS0
xxtime <- 1*c(xx$time<time)
pXXA <- ncol(cr2$E)-1
EA <- cr2$E[timeb,-1,drop=FALSE]
gammasEs <- .Call("CubeMattime",gammatt,EA,pXXA,p,pXXA,1,0,1,0,PACKAGE="mets")$XXX
gammasE <- matrix(0,length(xx$strata),1)
gammattt <- matrix(0,length(xx$strata),pXXA*1)
jumpsCt <- jumpsC[timeb]
gammasE[jumpsCt,] <- gammasEs
gammattt[jumpsCt,] <- gammatt
gammaEsdLam0 <- apply(gammasE*S0i*xxtime,2,cumsumstrata,xx$strata,xx$nstrata)
gammadLam0 <- apply(gammattt*S0i*xxtime,2,cumsumstrata,xx$strata,xx$nstrata)
XgammadLam0 <- .Call("CubeMattime",gammadLam0,xx$X[,-1],pXXA,p,pXXA,1,0,1,0,PACKAGE="mets")$XXX
Ut <- Et <- matrix(0,length(xx$strata),1)
Ut[jumpsCt,] <- augmentt.times
MGCt <- Ut[,drop=FALSE]-(XgammadLam0-gammaEsdLam0)*c(rr0)
MGCiid <- apply(MGCt,2,sumstrata,xx$id,nid)
MGCiid <- MGCiid/nid
#
nid <- max(cr2$id)
ids <- headstrata(cr2$id-1,nid)
ids <- cr2$call.id[ids]
res <- list(MGCiid=MGCiid,gammat=gammatt,augment=augment.times, id=ids,n=nid)
} ## }}}sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: aarch64-apple-darwin24.2.0
#> Running under: macOS Sequoia 15.2
#>
#> Matrix products: default
#> BLAS: /Users/klaus/.asdf/installs/R/4.4.2/lib/R/lib/libRblas.dylib
#> LAPACK: /Users/klaus/.asdf/installs/R/4.4.2/lib/R/lib/libRlapack.dylib; LAPACK version 3.12.0
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: Europe/Copenhagen
#> tzcode source: internal
#>
#> attached base packages:
#> [1] splines stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] prodlim_2024.06.25 ggplot2_3.5.1 cowplot_1.1.3 mets_1.3.5
#> [5] timereg_2.0.6 survival_3.8-3
#>
#> loaded via a namespace (and not attached):
#> [1] sass_0.4.9 future_1.34.0 generics_0.1.3
#> [4] lattice_0.22-6 listenv_0.9.1 digest_0.6.37
#> [7] magrittr_2.0.3 evaluate_1.0.1 grid_4.4.2
#> [10] mvtnorm_1.3-2 fastmap_1.2.0 jsonlite_1.8.9
#> [13] Matrix_1.7-1 scales_1.3.0 isoband_0.2.7
#> [16] codetools_0.2-20 numDeriv_2016.8-1.1 jquerylib_0.1.4
#> [19] lava_1.8.1 cli_3.6.3 rlang_1.1.4
#> [22] parallelly_1.41.0 future.apply_1.11.3 munsell_0.5.1
#> [25] withr_3.0.2 cachem_1.1.0 yaml_2.3.10
#> [28] tools_4.4.2 parallel_4.4.2 ucminf_1.2.2
#> [31] dplyr_1.1.4 colorspace_2.1-1 globals_0.16.3
#> [34] vctrs_0.6.5 R6_2.5.1 lifecycle_1.0.4
#> [37] MASS_7.3-64 pkgconfig_2.0.3 bslib_0.8.0
#> [40] pillar_1.10.1 gtable_0.3.6 data.table_1.16.4
#> [43] glue_1.8.0 Rcpp_1.0.13-1 xfun_0.50
#> [46] tibble_3.2.1 tidyselect_1.2.1 knitr_1.49
#> [49] farver_2.1.2 htmltools_0.5.8.1 labeling_0.4.3
#> [52] rmarkdown_2.29 compiler_4.4.2