---
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 can fit an exponential/linear link model with 
IPCW adjustment for a specific time-point for the resticted mean survival or years lost for 
competing risks. 
The function can be used for large data and is completely scalable, that is, linear in the data. 
The function can thus be used for large data. 
A nice feature is that influcence functions are computed and are available.

In addition and to summarize

  - the censoring weights can be strata dependent
  - predictions can be computed with standard errors 
  - computation time is linear in data
     - including standard errors
  - clusters can be given and then cluster corrected standard errors are computed



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 fully saturated model with full censoring model this is equal to the integrals of the Kaplan-Meier estimators as
illustrated below. 

We can also compute the integral of the Kaplan-Meier or Cox based  survival estimator to
get the RMST (with 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 IPCW estimator since  
\[ E( I(\epsilon=1) ( t - T \wedge t ) | X) = \int_0^t F_1(s) ds. \] 
For fully saturated model with full 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 
=========================

Computes average treatment effect for restricted mean survival and years lost in competing risks situation

```{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()
```