---
title: "Code used in the video vignette"
author: "Martijn Schuemie"
date: "`r Sys.Date()`"
output:
  pdf_document: default
  html_document: default
subtitle: A short demonstration of the EvidenceSynthesis package
vignette: >
  %\VignetteIndexEntry{Code used in the video vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
library(EvidenceSynthesis)
```

This vignette contains the code used in a short video on the EvidenceSynthesis package: [https://youtu.be/dho7E97vpgQ](https://youtu.be/dho7E97vpgQ).

# Simulate data 

Simulate 10 sites:

```{r}
simulationSettings <- createSimulationSettings(
  nSites = 10,
  n = 10000,
  treatedFraction = 0.8,
  nStrata = 5,
  hazardRatio = 2,
  randomEffectSd = 0.5
)
set.seed(1)
populations <- simulatePopulations(simulationSettings)

head(populations[[1]])
table(populations[[1]][, c("x", "y")])
```

# Fit a model locally

Assume we are at site 1:

```{r message = FALSE}
library(Cyclops)

population <- populations[[1]]

cyclopsData <- createCyclopsData(Surv(time, y) ~ x + strata(stratumId),
  data = population,
  modelType = "cox"
)
cyclopsFit <- fitCyclopsModel(cyclopsData)

# Hazard ratio:
exp(coef(cyclopsFit))

# 95% confidence interval:
exp(confint(cyclopsFit, parm = "x")[2:3])
```

# Approximate the likelihood function at one site

## Normal approximation

```{r}
normalApproximation <- approximateLikelihood(
  cyclopsFit = cyclopsFit,
  parameter = "x",
  approximation = "normal"
)
normalApproximation

plotLikelihoodFit(
  approximation = normalApproximation,
  cyclopsFit = cyclopsFit,
  parameter = "x"
)
```

## Adaptive approximation

```{r}
approximation <- approximateLikelihood(
  cyclopsFit = cyclopsFit,
  parameter = "x",
  approximation = "adaptive grid",
  bounds = c(log(0.1), log(10))
)
head(approximation)

plotLikelihoodFit(
  approximation = approximation,
  cyclopsFit = cyclopsFit,
  parameter = "x"
)
```

# Approximate at all sites

```{r}
fitModelInDatabase <- function(population, approximation) {
  cyclopsData <- createCyclopsData(Surv(time, y) ~ x + strata(stratumId),
    data = population,
    modelType = "cox"
  )
  cyclopsFit <- fitCyclopsModel(cyclopsData)
  approximation <- approximateLikelihood(cyclopsFit,
    parameter = "x",
    approximation = approximation
  )
  return(approximation)
}
adaptiveGridApproximations <- lapply(
  X = populations,
  FUN = fitModelInDatabase,
  approximation = "adaptive grid"
)
normalApproximations <- lapply(
  X = populations,
  FUN = fitModelInDatabase,
  approximation = "normal"
)
normalApproximations <- do.call(rbind, (normalApproximations))
```

# Synthesize evidence

## Fixed-effects

Gold standard (pooling data):

```{r message = FALSE,cache = TRUE}
fixedFxPooled <- computeFixedEffectMetaAnalysis(populations)
fixedFxPooled
```

Normal approximation: 

```{r message = FALSE}
fixedFxNormal <- computeFixedEffectMetaAnalysis(normalApproximations)
fixedFxNormal
```

Adaptive grid approximation:

```{r message = FALSE}
fixedFxAdaptiveGrid <- computeFixedEffectMetaAnalysis(adaptiveGridApproximations)
fixedFxAdaptiveGrid
```

### Visualization

Normal approximation: 

```{r message = FALSE, warning = FALSE, fig.width = 9, fig.height = 5}
plotMetaAnalysisForest(
  data = normalApproximations,
  labels = paste("Site", 1:10),
  estimate = fixedFxNormal,
  xLabel = "Hazard Ratio"
)
```

Adaptive grid approximation:

```{r message = FALSE, warning = FALSE, fig.width = 9, fig.height = 5}
plotMetaAnalysisForest(
  data = adaptiveGridApproximations,
  labels = paste("Site", 1:10),
  estimate = fixedFxAdaptiveGrid,
  xLabel = "Hazard Ratio"
)
```

## Random-effects

Gold standard (pooling data):

```{r cache = TRUE, message = FALSE,cache=TRUE}
randomFxPooled <- computeBayesianMetaAnalysis(populations)
exp(randomFxPooled[, 1:3])
```

Normal approximation: 

```{r message = FALSE}
randomFxNormal <- computeBayesianMetaAnalysis(normalApproximations)
exp(randomFxNormal[, 1:3])
```

Adaptive grid approximation:

```{r message = FALSE}
randomFxAdaptiveGrid <- computeBayesianMetaAnalysis(adaptiveGridApproximations)
exp(randomFxAdaptiveGrid[, 1:3])
```

### Visualization

Normal approximation: 

```{r message = FALSE, warning = FALSE, fig.width = 8, fig.height = 5}
plotMetaAnalysisForest(
  data = normalApproximations,
  labels = paste("Site", 1:10),
  estimate = randomFxNormal,
  xLabel = "Hazard Ratio"
)
```

Adaptive grid approximation:

```{r message = FALSE, warning = FALSE, fig.width = 8, fig.height = 5}
plotMetaAnalysisForest(
  data = adaptiveGridApproximations,
  labels = paste("Site", 1:10),
  estimate = randomFxAdaptiveGrid,
  xLabel = "Hazard Ratio"
)
```