--- title: "While Alive estimands for Recurrent Events" author: Klaus Holst & Thomas Scheike date: "`r Sys.Date()`" output: rmarkdown::html_vignette: fig_caption: yes fig_width: 7.15 fig_height: 5.5 vignette: > %\VignetteIndexEntry{While Alive estimands for Recurrent Events} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(mets) ``` While Alive estimands for Recurrent Events ========================================== 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. ```{r} data(hfaction_cpx12) dtable(hfaction_cpx12,~status) dd <- WA_recurrent(Event(entry,time,status)~treatment+cluster(id),hfaction_cpx12,time=2, death.code=2) summary(dd) 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]