probability-sensitivity-analysis.Rmd

---
title: "Probability Sensitivity Analysis (PSA)"
author: "Nathan Green"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
bibliography: bibliography.bib
vignette: >
  %\VignetteIndexEntry{Probability Sensitivity Analysis (PSA)}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

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

## Introduction

PSA is a core part of any cost-effectiveness analysis [@Briggs2012].
Here we will carry this out for a simple decision tree.
This involves repeatedly sampling from a distribution for each branch probability and cost and calculating the total expected value for each set of realisations.

## Example

```{r warning=FALSE, message=FALSE}
library(CEdecisiontree, quietly = TRUE)
library(assertthat, quietly = TRUE)
library(tibble, quietly = TRUE)
library(tidyverse, quietly = TRUE)
library(purrr, quietly = TRUE)
```

We first define the decision tree.
The difference to previous trees is that we now use the _list-column_ feature to define distributions rather than point values.

```{r}
tree_dat <-
  list(child = list("1" = c(2, 3),
                    "2" = NULL,
                    "3" = NULL),
       dat = tibble(
         node = 1:3,
         prob =
           list(
             NA_real_,
             list(distn = "unif", params = c(min = 0, max = 1)),
             list(distn = "unif", params = c(min = 0, max = 1))),
         vals =
           list(
             0L,
             list(distn = "unif", params = c(min = 0, max = 1)),
             list(distn = "unif", params = c(min = 0, max = 1)))))

tree_dat
```

We can now loop over this tree and generate samples of values for the given distributions.
We use the `sample_distributions()` function to do this.
This is all wrapped up in the convenience function `create_psa_inputs`.

```{r warning=FALSE}
library(purrr)

tree_dat_sa <- create_psa_inputs(tree_dat, n = 1000)
```

This results in a list of trees.
```{r}
head(tree_dat_sa, 2)
```

Now it is straightforward to map over each of these trees to obtain the total expected values

```{r}
res <- map_dbl(tree_dat_sa, dectree_expected_values)
head(res)


hist(res, breaks = 20)
```


## References