rmarkdown_analysis_template.rmd

---
title: "Meta-analysis of survival models in the DataSHIELD platform"
author: "Soumya Banerjee and Tom Bishop"
date: "30 April 2021"
output:
  pdf_document:
    toc: yes
    number_sections: yes
    fig_caption: yes
  html_document:
    toc: yes
    df_print: paged
urlcolor: blue
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=TRUE)
```

# Summary

This is a document that outlines technical notes and a development plan for implementing survival models
and meta-analyzing hazard ratios in the DataSHIELD platform.

# Survival analysis in DataSHIELD

We outline the development and code for implementing survival models and meta-analysis
of hazard ratios in DataSHIELD.

All code is available here:

   * https://github.com/neelsoumya/dsBaseClient/tree/absolute_newbie_client
   
   * https://github.com/neelsoumya/dsBase/tree/absolute_newbie
   
   * https://github.com/neelsoumya/datashield_testing_basic/blob/master/development_plan.rmd
   
   * https://github.com/neelsoumya/datashield_testing_basic/blob/master/development_plan.pdf
   

The computational steps are outlined below. 
The first step is connecting to the server and loading the survival data.
We assume that the reader is familiar with these details.


```{r, include = FALSE}

####################
# Load library
####################
library(knitr)
library(rmarkdown)
library(tinytex)
library(survival)
library(metafor)
library(ggplot2)
library(survminer)
require('DSI')
require('DSOpal')
require('dsBaseClient')

#######################
# Get data
#######################
builder <- DSI::newDSLoginBuilder()

builder$append(server = "study1", 
               url = "http://192.168.56.100:8080/", 
               user = "administrator", password = "datashield_test&", 
               table = "SURVIVAL.EXPAND_NO_MISSING1", driver = "OpalDriver")
builder$append(server = "study2", 
               url = "http://192.168.56.100:8080/", 
               user = "administrator", password = "datashield_test&", 
               table = "SURVIVAL.EXPAND_NO_MISSING2", driver = "OpalDriver")
builder$append(server = "study3", 
               url = "http://192.168.56.100:8080/", 
               user = "administrator", password = "datashield_test&", 
               table = "SURVIVAL.EXPAND_NO_MISSING3", driver = "OpalDriver")

logindata <- builder$build()

##############
# login
##############

# Log onto the remote Opal training servers
connections <- DSI::datashield.login(logins = logindata, assign = TRUE, symbol = "D") 

```

# Creating server-side variables for survival analysis

We now outline some steps for analysing survival data.

* make sure that the outcome variable is numeric 

```
ds.asNumeric(x.name = "D$cens",
             newobj = "EVENT",
             datasources = connections)

ds.asNumeric(x.name = "D$survtime",
             newobj = "SURVTIME",
             datasources = connections)
```

*  convert time id variable to a factor 

```

ds.asFactor(input.var.name = "D$time.id",
            newobj = "TID",
            datasources = connections)

```

*  create in the server-side the log(survtime) variable

```

ds.log(x = "D$survtime",
       newobj = "log.surv",
       datasources = connections)

```

*  create start time variable

```

ds.asNumeric(x.name = "D$starttime",
             newobj = "STARTTIME",
             datasources = connections)

ds.asNumeric(x.name = "D$endtime",
             newobj = "ENDTIME",
             datasources = connections)

```



```{r, include=FALSE}

# make sure that the outcome is numeric 
ds.asNumeric(x.name = "D$cens",
             newobj = "EVENT",
             datasources = connections)

ds.asNumeric(x.name = "D$survtime",
             newobj = "SURVTIME",
             datasources = connections)

# convert time id variable to a factor 
ds.asFactor(input.var.name = "D$time.id",
            newobj = "TID",
            datasources = connections)

# create in the server-side the log(survtime) variable
ds.log(x = "D$survtime",
       newobj = "log.surv",
       datasources = connections)

# create start time variable
ds.asNumeric(x.name = "D$starttime",
             newobj = "STARTTIME",
             datasources = connections)

# create end time variable
ds.asNumeric(x.name = "D$endtime",
             newobj = "ENDTIME",
             datasources = connections)


# check which variables exist
dsBaseClient::ds.ls()

# call new function modify NAMESPACE and DESCRIPTION
# call coxph server side
# client side function is here:
# https://github.com/neelsoumya/dsBaseClient/blob/absolute_newbie_client/R/ds.coxph.SLMA.R
# server side function is here:
# https://github.com/neelsoumya/dsBase/blob/absolute_newbie/R/coxphSLMADS.R
```


# Create survival object and call coxphSLMA()


* use constructed Surv object in *coxph.SLMA()*

```
dsBaseClient::ds.Surv(time='STARTTIME', time2='ENDTIME', 
                      event = 'EVENT', objectname='surv_object',
                      type='counting')

coxph_model_full <- dsBaseClient::ds.coxph.SLMA(formula = 'surv_object~D$age+D$female')

```

* use direct inline call to *survival::Surv()*

```

ds.coxph.SLMA(formula = 'survival::Surv(time=SURVTIME,event=EVENT)~D$age+D$female', 
                dataName = 'D', 
                datasources = connections)

```

* call with *survival::strata()*

```

coxph_model_strata <- dsBaseClient::ds.coxph.SLMA(formula = 'surv_object~D$age + 
                          survival::strata(D$female)')

summary(coxph_model_strata)

```


```{r, include=FALSE}

# 1. use constructed surv object in coxph
dsBaseClient::ds.Surv(time='STARTTIME', time2='ENDTIME', event = 'EVENT', objectname='surv_object', type='counting')
coxph_model_full <- dsBaseClient::ds.coxph.SLMA(formula = 'surv_object~D$age+D$female')

# 2. use direct inline call
ds.coxph.SLMA(formula = 'survival::Surv(time=SURVTIME,event=EVENT)~D$age+D$female', dataName = 'D', datasources = connections)

# 3. Call with strata()
coxph_model_strata <- dsBaseClient::ds.coxph.SLMA(formula = 'surv_object~D$age+survival::strata(D$female)')
summary(coxph_model_strata)

```

# Summary of survival objects

We can also summarize a server-side object of type *survival::Surv()*
using a call to *ds.summary()*. This will provide a non-disclosive
summary of the server-side object. An example call is shown below:


```

dsBaseClient::ds.summary(x = 'surv_object')

```



```{r, include=FALSE}

#################################
# summary of coxphSLMA
#################################
dsBaseClient::ds.summary(x = 'surv_object')
dsBaseClient::ds.class(x = 'surv_object')
dsBaseClient::ds.mean(x='surv_object')


#################################
# TODO: Plot survival curves
#################################
# fit <- survival::survfit(formula = 'surv_object~D$age+D$female', data = 'D')
# need ds.survfit() and survfitDS()
# fit_model <- ds.survfit(coxph_model[1])
# plot(fit_model)

# TODO:
# plot(survfit_km, fun="cloglog")

# TODO: 
# ggplot like functionality see other functions

# In dsBaseClient::
# ds.survfit()
#         datashield.aggregate("survfitDS", ....)
#          return (the fit model)

# In dsBase::
# survfitDS(coxph_model)  
#               fit_model <- survival::survfit(coxph_model, newdata = 'D')
#               return (fit_model)

dsBaseClient::ds.survfit(formula='surv_object~1', objectname='survfit_object')
# verify that object has been created
dsBaseClient::ds.ls()

# ERROR summary of survfit not allowed
# dsBaseClient::ds.summary(x='survfit_object')
# dsBaseClient::ds.mean(x='survfit_object')

```


# Diagnostics for Cox proportional hazards models

We have also created functions to test for the assumptions of 
Cox proportional hazards models.

```
dsBaseClient::ds.coxphSLMAassign(formula = 'surv_object~D$age+D$female',
                            objectname = 'coxph_serverside')

dsBaseClient::ds.cox.zphSLMA(fit = 'coxph_serverside')

dsBaseClient::ds.coxphSummary(x = 'coxph_serverside')

```

A diagnostic summary is shown below. 

```{r, echo=FALSE}

dsBaseClient::ds.coxphSLMAassign(formula = 'surv_object~D$age+D$female', objectname = 'coxph_serverside')

dsBaseClient::ds.cox.zphSLMA(fit = 'coxph_serverside')

dsBaseClient::ds.coxphSummary(x = 'coxph_serverside')

```



# Meta-analyze hazard ratios

We now outline how the hazard ratios from the survival models are meta-analyzed. 
We use the *metafor* package for meta-analysis. We show the summary of an 
example meta-analysis and a forest plot below. The forest plot shows a basic
example of meta-analyzed hazard ratios from a survival model (analyzed in
DataSHIELD).

The log-hazard ratios and their standard erros from each study can be
found after running *ds.coxphSLMA()*

The hazard ratios can then be meta-analyzed:

```

metafor::rma(log_hazard_ratio, sei = se_hazard_ratio, method = 'REML')

```

A summary of this meta-analyzed model is shown below.


```{r, echo=FALSE}

# TODO: for each study
for (i_temp_counter in c(1:length(coxph_model_full)))
{
      
}

# list of hazard ratios for first parameter (age) over 3 studies 
input_logHR = c(coxph_model_full$study1$coefficients[1,2], coxph_model_full$study2$coefficients[1,2], coxph_model_full$study3$coefficients[1,2])
input_se    = c(coxph_model_full$study1$coefficients[1,3], coxph_model_full$study2$coefficients[1,3], coxph_model_full$study3$coefficients[1,3])

meta_model <- metafor::rma(input_logHR, sei = input_se, method = 'REML')
summary(meta_model)

#######################################################
# forest plots of final meta-analyzed hazard ratios
#######################################################
```

We now show a forest plot with the meta-analyzed hazard ratios. The hazard ratios come
from the DataSHIELD function *ds.coxphSLMA()*. The hazard ratios are meta-analyzed
using the *metafor* package.

```{r fig.cap="Example forest plot of meta-analyzed hazard ratios.", echo=FALSE}

metafor::forest.rma(x = meta_model)

```    

# Plotting of privacy-preserving survival curves

We also plot privacy preserving survival curves.

```


dsBaseClient::ds.survfit(formula='surv_object~1', objectname='survfit_object')

dsBaseClient::ds.plotsurvfit(formula = 'survfit_object')


```


```{r, echo = FALSE}

dsBaseClient::ds.plotsurvfit(formula = 'survfit_object')

```

\newpage
   
```{r, include=FALSE}

#############################################
# disconnect
#############################################
DSI::datashield.logout(conns = connections)

```

# Dealing with preserving privacy and disclosure checks

Disclosure checks are an integral part of DataSHIELD. We outline some of the checks we have built in to ensure privacy of individual level data. We also outline some issues that need to be discussed.   



# Acknowledgements

We acknowledge the help and support of the DataSHIELD technical team.
We are especially grateful to Paul Burton, Demetris Avraam, Stuart Wheater and
Patricia Ryser-Welch for fruitful discussions and feedback.


# References

  * https://github.com/datashield
  
  * http://www.metafor-project.org
  
  * https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6015946