The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
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.2These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.