---
title: "Average treatment effect (ATE) for Restricted mean survival and years lost of Competing risks"
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{Average treatment effect (ATE) for Restricted mean survival and years lost of Competing risks}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(mets)
```


The `resmeanATE` or `resmeanIPCW`/`rmstIPCW` functions fit an exponential or linear link model with
IPCW adjustment at a specific time point for the restricted mean survival or years lost in
competing risks. The computation is linear in the data size, making it suitable for large datasets.
Influence functions are computed and available for downstream inference.

Key features:

  - censoring weights can be stratum-dependent
  - predictions can be computed with standard errors
  - computation time is linear in data size, including standard errors
  - cluster-corrected standard errors are available via the `clusters` argument



RMST 
====

Regression for rmst outcome $E(T \wedge t | X) = exp(X^T \beta)$ based on IPCW adjustment for censoring, thus solving
the estimating equation 
\begin{align*}
 & X^T [ (T \wedge t) \frac{I(C > T \wedge t)}{G_C(T \wedge t,X)} - exp(X^T \beta) ] = 0 .
\end{align*}
This is done with the resmeanIPCW function. 
For a fully saturated model with a complete censoring model this is equal to the integral of the Kaplan-Meier estimator, as
illustrated below.

We can also compute the integral of the Kaplan-Meier or Cox-based survival estimator to
get the RMST (using the `resmean_phreg` function):
\[ \int_0^t S(s|X) ds. \]

For competing risks, the years lost can be computed via cumulative incidence functions (`cif_yearslost`)
or as an IPCW estimator, since
\[ E( I(\epsilon=1) ( t - T \wedge t ) | X) = \int_0^t F_1(s) ds. \]
For a fully saturated model with a complete censoring model these estimators are equivalent, as illustrated below.


```{r}
set.seed(101)

     data(bmt); bmt$time <- bmt$time+runif(nrow(bmt))*0.001

     # E( min(T;t) | X ) = exp( a+b X) with IPCW estimation 
     out <- resmeanIPCW(Event(time,cause!=0)~tcell+platelet+age,bmt,
                     time=50,cens.model=~strata(platelet),model="exp")
     summary(out)
     
      ### same as Kaplan-Meier for full censoring model 
     bmt$int <- with(bmt,strata(tcell,platelet))
     out <- resmeanIPCW(Event(time,cause!=0)~-1+int,bmt,time=30,
                                  cens.model=~strata(platelet,tcell),model="lin")
     estimate(out)
     out1 <- phreg(Surv(time,cause!=0)~strata(tcell,platelet),data=bmt)
     rm1 <- resmean_phreg(out1,times=30)
     summary(rm1)
     
     ## competing risks years-lost for cause 1  
     out <- resmeanIPCW(Event(time,cause)~-1+int,bmt,time=30,cause=1,
                                 cens.model=~strata(platelet,tcell),model="lin")
     estimate(out)
     ## same as integrated cumulative incidence 
     rmc1 <- cif_yearslost(Event(time,cause)~strata(tcell,platelet),data=bmt,times=30,cause=1)
     summary(rmc1)

     ## plotting the years lost for different horizon's and the two causes 
     par(mfrow=c(1,3))
     plot(rm1,years.lost=TRUE,se=1)
     ## cause refers to column of cumhaz for the different causes
     plot(rmc1,cause=1,se=1)
     plot(rmc1,cause=2,se=1)
```

Based on the output from the IPCW estimators we can derive any measure of interest

```{r}
estimate(out)

measures <- function(p) {
 ratio1 <- p[1]/p[2]; ratio2 <- p[2]/p[1]; dif1 <- p[4]-p[1]; dif2 <- p[3]-p[1]
 m <- c(dif1,dif2,ratio1,ratio2)
 return(m)
}

 labs <- c("dif4-1","dif3-1","ratio 1/2","ratio 2/1")
 estimate(out,f=measures,labels=labs)
```

Average treatment effect
=========================

This section computes the average treatment effect for restricted mean survival and years lost in the competing risks setting.

```{r}
 dfactor(bmt) <- tcell~tcell
 bmt$event <- (bmt$cause!=0)*1
 out <- resmeanATE(Event(time,event)~tcell+platelet,data=bmt,time=40,treat.model=tcell~platelet)
 summary(out)
 
 out1 <- resmeanATE(Event(time,cause)~tcell+platelet,data=bmt,cause=1,time=40,
		    treat.model=tcell~platelet)
 summary(out1)

 out2 <- resmeanATE(Event(time,cause)~tcell+platelet,data=bmt,cause=2,time=40,
		    treat.model=tcell~platelet)
 summary(out2)
```



SessionInfo
============


```{r}
sessionInfo()
```