Code used in the video vignette"

library(EvidenceSynthesis)

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

Simulate data

Simulate 10 sites:

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:

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

normalApproximation <- approximateLikelihood(
  cyclopsFit = cyclopsFit,
  parameter = "x",
  approximation = "normal"
)
normalApproximation

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

Adaptive approximation

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

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):

fixedFxPooled <- computeFixedEffectMetaAnalysis(populations)
fixedFxPooled

Normal approximation:

fixedFxNormal <- computeFixedEffectMetaAnalysis(normalApproximations)
fixedFxNormal

Adaptive grid approximation:

fixedFxAdaptiveGrid <- computeFixedEffectMetaAnalysis(adaptiveGridApproximations)
fixedFxAdaptiveGrid

Visualization

Normal approximation:

plotMetaAnalysisForest(
  data = normalApproximations,
  labels = paste("Site", 1:10),
  estimate = fixedFxNormal,
  xLabel = "Hazard Ratio"
)

Adaptive grid approximation:

plotMetaAnalysisForest(
  data = adaptiveGridApproximations,
  labels = paste("Site", 1:10),
  estimate = fixedFxAdaptiveGrid,
  xLabel = "Hazard Ratio"
)

Random-effects

Gold standard (pooling data):

randomFxPooled <- computeBayesianMetaAnalysis(populations)
exp(randomFxPooled[, 1:3])

Normal approximation:

randomFxNormal <- computeBayesianMetaAnalysis(normalApproximations)
exp(randomFxNormal[, 1:3])

Adaptive grid approximation:

randomFxAdaptiveGrid <- computeBayesianMetaAnalysis(adaptiveGridApproximations)
exp(randomFxAdaptiveGrid[, 1:3])

Visualization

Normal approximation:

plotMetaAnalysisForest(
  data = normalApproximations,
  labels = paste("Site", 1:10),
  estimate = randomFxNormal,
  xLabel = "Hazard Ratio"
)

Adaptive grid approximation:

plotMetaAnalysisForest(
  data = adaptiveGridApproximations,
  labels = paste("Site", 1:10),
  estimate = randomFxAdaptiveGrid,
  xLabel = "Hazard Ratio"
)


Try the EvidenceSynthesis package in your browser

Any scripts or data that you put into this service are public.

EvidenceSynthesis documentation built on May 31, 2023, 9:05 p.m.