---
title: "Calculating prevalence"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{a04_Calculating_prevalence}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  eval = Sys.getenv("$RUNNER_OS") != "macOS"
)
```

```{r, message= FALSE, warning=FALSE, echo=FALSE}
library(here)
```

# Introduction
Prevalence is the total number of people with an ongoing health-related event, such as a medical condition or medication use, at a particular time or during a given period divided by the population at risk. In the previous vignettes we have seen how we can identify a denominator population and define and instantiate an outcome cohort. Prevalence then can be calculated to describe the proportion of people in the denominator population who are in the outcome cohort at a specified time point (point prevalence) or over a given time interval (period prevalence).

In the first plot below, we can We can see at time t+2 that 2 out of 5 people were in an outcome cohort, giving a point prevalence of 40%. In the second figure, period prevalence between t+2 and t+3 was also 40%. However for period prevalence between t and t+1, what do we do with those people who only contributed some time during the period? If we included them we´ll have a period prevalence of 20%, whereas if we require that everyone is observed for the full period to contribute then we´ll have a period prevalence of 33%.  

```{r,echo=FALSE, out.width="80%"}
knitr::include_graphics(here("vignettes/point_prev.png"))
```

```{r,echo=FALSE, out.width="80%"}
knitr::include_graphics(here("vignettes/period_prev.png"))
```

# Outcome definition
Outcome cohorts are defined externally. When creating outcome cohorts for estimating prevalence, important considerations relate to whether to restrict events to the first occurrence in an individuals history or not and how cohort exit is defined. These decisions will necessarily be based on the nature of the proposed outcome (e.g., whether it is an acute or chronic condition) and the research question being investigated. 

In addition, it is typically not recommended to include exclusion requirements when creating outcome cohorts. Restrictions on patient characteristics can be specified when identifying the denominator cohort using `generateDenominatorCohortSet()` or `generateTargetDenominatorCohortSet()`

# Using estimatePointPrevalence() and estimatePeriodPrevalence()
`estimatePointPrevalence()` and `estimatePeriodPrevalence()` are the functions we use to estimate prevalence. To demonstrate its use, let´s load the IncidencePrevalence package (along with a couple of packages to help for subsequent plots) and generate 50,000 example patients using the `mockIncidencePrevalence()` function from whom we´ll create a denominator population without adding any restrictions other than a study period.

```{r, message=FALSE, warning=FALSE}
library(IncidencePrevalence)
library(dplyr)
library(tidyr)

cdm <- mockIncidencePrevalence(
  sampleSize = 20000,
  earliestObservationStartDate = as.Date("1960-01-01"), 
  minOutcomeDays = 365,
  outPre = 0.3
)

cdm <- generateDenominatorCohortSet(
  cdm = cdm, name = "denominator",
  cohortDateRange = c(as.Date("1990-01-01"), as.Date("2009-12-31")),
  ageGroup = list(c(0, 150)),
  sex = "Both",
  daysPriorObservation = 0
)

cdm$denominator %>%
  glimpse()
```

## Using estimatePointPrevalence()
Let´s first calculate point prevalence on a yearly basis. 
```{r, message= FALSE, warning=FALSE}
prev <- estimatePointPrevalence(
  cdm = cdm,
  denominatorTable = "denominator",
  outcomeTable = "outcome",
  interval = "Years"
)

prev %>%
  glimpse()

plotPrevalence(prev, ylim = c(0, NA))
```

We can also calculate point prevalence by calendar month.
```{r, message= FALSE, warning=FALSE}
prev <- estimatePointPrevalence(
  cdm = cdm,
  denominatorTable = "denominator",
  outcomeTable = "outcome",
  interval = "Months"
)

prev %>%
  glimpse()

plotPrevalence(prev, ylim = c(0, NA))
```

By using the estimatePointPrevalence() function, we can further specify where to compute point prevalence in each time interval (start, middle, end). By default, this parameter is set to start. But we can use middle instead like so:

```{r, message= FALSE, warning=FALSE}
prev <- estimatePointPrevalence(
  cdm = cdm,
  denominatorTable = "denominator",
  outcomeTable = "outcome",
  interval = "Years",
  timePoint = "middle"
)

prev %>%
  glimpse()

plotPrevalence(prev, ylim = c(0, NA), options = list(line = FALSE))
```

## Using estimatePeriodPrevalence()
To calculate period prevalence by year (i.e. each period is a calendar year)
```{r, message= FALSE, warning=FALSE}
prev <- estimatePeriodPrevalence(
  cdm = cdm,
  denominatorTable = "denominator",
  outcomeTable = "outcome",
  interval = "Years"
)

prev %>%
  glimpse()

plotPrevalence(prev)
```


To calculate period prevalence by month (i.e. each period is a calendar month)
```{r, message= FALSE, warning=FALSE}
prev <- estimatePeriodPrevalence(
  cdm = cdm,
  denominatorTable = "denominator",
  outcomeTable = "outcome",
  interval = "Months"
)

prev %>%
  glimpse()

plotPrevalence(prev, ylim = c(0, NA))
```

When using the estimatePeriodPrevalence() function, we can set the fullContribution parameter to decide whether individuals are required to be present in the database throughout the interval of interest in order to be included (fullContribution=TRUE). If not, individuals will only be required to be present for one day of the interval to contribute (fullContribution=FALSE), which would be specified like so:  

```{r, message= FALSE, warning=FALSE}
prev <- estimatePeriodPrevalence(
  cdm = cdm,
  denominatorTable = "denominator",
  outcomeTable = "outcome",
  interval = "Months",
  fullContribution = FALSE
)

prev %>%
  glimpse()

plotPrevalence(prev)
```

## Stratified analyses
If we specified multiple denominator populations results will be returned for each. Here for example we define three age groups for denominator populations and get three sets of estimates back when estimating prevalence.

```{r}
cdm <- generateDenominatorCohortSet(
  cdm = cdm, name = "denominator_age_sex",
  cohortDateRange = c(as.Date("1990-01-01"), as.Date("2009-12-31")),
  ageGroup = list(c(0, 39),
                  c(41, 65),
                  c(66, 150)),
  sex = "Both",
  daysPriorObservation = 0
)
prev <- estimatePeriodPrevalence(
  cdm = cdm,
  denominatorTable = "denominator_age_sex",
  outcomeTable = "outcome"
)

plotPrevalence(prev, facet = "denominator_age_group")
```

While we specify time-varying stratifications when defining our denominator populations, if we have time-invariant stratifications we can include these at the the estimation stage. To do this we first need to add a new column to our denominator cohort with our stratification variable. Here we´ll add an example stratification just to show the idea. Note, as well as getting stratified results we´ll also get overall results.

```{r}
cdm$denominator <- cdm$denominator %>% 
  mutate(group = if_else(as.numeric(subject_id)  < 500, "first", "second")) 

prev <- estimatePeriodPrevalence(
  cdm = cdm,
  denominatorTable = "denominator",
  outcomeTable = "outcome",
  strata = "group"
)

plotPrevalence(prev, 
               colour = c("group"))
```

We can also stratify on multiple variables at the same time.

```{r}
cdm$denominator <- cdm$denominator %>% 
  mutate(group_1 = if_else(as.numeric(subject_id)  < 1500, "first", "second"),
         group_2 = if_else(as.numeric(subject_id)  < 1000, "one", "two"))

prev <- estimatePeriodPrevalence(
  cdm = cdm,
  denominatorTable = "denominator",
  outcomeTable = "outcome",
  strata = list(c("group_1"), # for just group_1
                c("group_2"), # for just group_2
                c("group_1", "group_2"))  # for group_1 and group_2
)

plotPrevalence(prev, 
               facet = c("group_1", "group_2"))
```

## Other parameters

In the examples above, we have used calculated prevalence by months and years, but it can be also calculated by weeks, months or for the entire time period observed (overall). In addition, the user can decide whether to include time intervals that are not fully captured in the database (e.g., having data up to June for the last study year when computing period prevalence rates). By default, incidence will only be estimated for those intervals where the database captures all the interval (completeDatabaseIntervals=TRUE).

Given that we can set `estimatePointPrevalence()` and `estimatePeriorPrevalence()` to exclude individuals based on certain parameters (e.g., fullContribution), it is important to note that the denominator population used to compute prevalence rates might differ from the one calculated with `generateDenominatorCohortSet()`. Along with our central estimate, 95 % confidence intervals are calculated using the Wilson Score method. 

## Attrition
`estimatePointPrevalence()` and `estimatePeriorPrevalence()` will generate a table with point and period prevalence rates for each of the time intervals studied and for each combination of the parameters set, respectively. We can also view attrition associated with performing the analysis.
```{r, message= FALSE, warning=FALSE}
prev <- estimatePeriodPrevalence(
  cdm = cdm,
  denominatorTable = "denominator",
  outcomeTable = "outcome",
  interval = "Years",
  fullContribution = TRUE
)
tablePrevalenceAttrition(prev, settingsColumns = character())
```