---
title: "Population diagnostics"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{PopulationDiagnostics}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
NOT_CRAN <- identical(tolower(Sys.getenv("NOT_CRAN")), "true")

knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
eval = NOT_CRAN
)
```

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>", message = FALSE, warning = FALSE,
  fig.width = 7
)

library(CDMConnector)
if (Sys.getenv("EUNOMIA_DATA_FOLDER") == "") Sys.setenv("EUNOMIA_DATA_FOLDER" = tempdir())
if (!dir.exists(Sys.getenv("EUNOMIA_DATA_FOLDER"))) dir.create(Sys.getenv("EUNOMIA_DATA_FOLDER"))
if (!eunomiaIsAvailable()) downloadEunomiaData(datasetName = "synpuf-1k")
```

## Introduction
In this example we're going to just create a cohort of individuals with an ankle sprain using the Eunomia synthetic data. 

```{r, message=FALSE, warning=FALSE}
library(CDMConnector)
library(CohortConstructor)
library(CodelistGenerator)
library(PatientProfiles)
library(IncidencePrevalence)
library(PhenotypeR)


con <- DBI::dbConnect(duckdb::duckdb(), 
                      CDMConnector::eunomiaDir("synpuf-1k", "5.3"))
cdm <- CDMConnector::cdmFromCon(con = con, 
                                cdmName = "Eunomia Synpuf",
                                cdmSchema   = "main",
                                writeSchema = "main", 
                                achillesSchema = "main")

cdm$injuries <- conceptCohort(cdm = cdm,
  conceptSet = list(
    "ankle_sprain" = 81151
  ),
  name = "injuries")
```

We can get the incidence and prevalence of our study cohort using `populationDiagnostics()`:
```{r}
pop_diag <- populationDiagnostics(cdm$injuries, 
                                  populationSample = 10000)
```

This function builds on [IncidencePrevalence](https://darwin-eu.github.io/IncidencePrevalence/index.html) R package to perform the following analyses:

-   **Incidence:** It estimates the incidence of our cohort using [estimateIncidence()](https://darwin-eu.github.io/IncidencePrevalence/reference/estimateIncidence.html).
-   **Prevalence:** It estimates the prevalence of our cohort on a year basis using [estimatePeriodPrevalence()](https://darwin-eu.github.io/IncidencePrevalence/reference/estimatePeriodPrevalence.html).

All analyses are performed for:

-   Overall and stratified by age groups: 0 to 17, 18 to 64, 65 to 150. Age groups cannot be modified.
-   Overall and stratified by sex (Female, Male).
-   Restricting the denominator population to those with 0 and 365 of days of prior observation.

## Visualising the results
We can use  [IncidencePrevalence](https://darwin-eu.github.io/IncidencePrevalence/index.html) package to visualise the results obtained.

### Incidence
```{r}
tableIncidence(pop_diag,     
               groupColumn = c("cdm_name", "outcome_cohort_name"),
               hide = "denominator_cohort_name",
               settingsColumn = c("denominator_age_group",
                         "denominator_sex",
                         "denominator_days_prior_observation",
                         "outcome_cohort_name"))
```
```{r}
results <- pop_diag |> 
  omopgenerics::filterSettings(result_type == "incidence") |>
  visOmopResults::filterAdditional(analysis_interval == "years")
plotIncidence(results,
              colour = "denominator_age_group",
              facet = c("denominator_sex", "denominator_days_prior_observation"))
```


### Prevalence
```{r}
tablePrevalence(pop_diag,     
               groupColumn = c("cdm_name", "outcome_cohort_name"),
               hide = "denominator_cohort_name",
               settingsColumn = c("denominator_age_group",
                         "denominator_sex",
                         "denominator_days_prior_observation",
                         "outcome_cohort_name"))
```

```{r}
results <- pop_diag |> 
  omopgenerics::filterSettings(result_type == "prevalence") |>
  visOmopResults::filterAdditional(analysis_interval == "years")
plotPrevalence(results,
               colour = "denominator_age_group",
               facet = c("denominator_sex", "denominator_days_prior_observation"))
```
