---
title: "{sitrep} outbreaks template vignette"
subtitle: "A guided example of how to use the R Markdown template to produce a situation report"
output: 
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Outbreaks_Walkthrough}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Overview

This page demonstrate use of the {sitrep} R Markdown template for analysis of an outbreak of Acute Jaundice Syndrome (AJS), using data from an outbreak in Am Timan, Chad. *This case study uses data that is not standardized with the MSF data dictionary, and so reflects a challenging use case.*

### AJS/HEV outbreak in Am Timan, Chad (2016-2017)

We will be using data and examples from a real outbreak of Hepatitis E virus (HEV) infection which occurred in the Chadian town of Am Timan between October 2016 and April 2017

*DISCLAIMER: The real data from the outbreak has been used for the training materials linked to the outbreak templates of the R4epis project. The Chadian Ministry of Health (MoH) has approved the use of this data for training purposes. Please note, that some data has been adapted in order to best achieve training objectives. Also, the GPS coordinates included in the dataset do not correspond to real cases identified during this outbreak. They have been generated exclusively for training purposes.*

In early September 2016, a cluster of severely ill pregnant women with jaundice was detected at the Am Timan Hospital in the maternity ward. Following rapid tests conducted in Am Timan and confirmatory testing in the Netherlands, this was confirmed as due to HEV infection. Thus the MoH and MSF outbreak response was activated.

The response consisted of four components:

1)  active community based surveillance (CBS)\
2)  clinical assessment and management of ill cases at the hospital\
3)  water chlorination activities at most water points in the town\
4)  hygiene promotion

**The outbreak linelist data was combined data from the active CBS data and the clinical data from the hospital.** The CBS functioned with community health workers visiting all households in the town every two weeks and actively searching for people with Acute Jaundice Syndrome (AJS). For this group of persons (suspected cases), only those that were visibly ill or that were pregnant were referred to the hospital for clinical assessment and admission if required. Persons that self-reported to the hospital or that arrived after referral would undergo a clinical assessment and a rapid test for HEV for diagnosis. Thus only for people assessed at the hospital were we able to capture a confirmed case status.

For the duration of the outbreak we detected 1193 suspected cases, 100 confirmed cases and discarded 150 cases with AJS who were not positive for HEV infection.

[Link to Spina and Lenglet *et al.* 2017](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0188240)

## Getting started

### Software installations

To begin using R, ensure you have the following free software installed:

1.  [R](https://cran.r-project.org/bin/windows/base/) (the free statistical software)
2.  [RStudio](https://www.rstudio.com/products/rstudio/download/#download) (the free user interface)
3.  [RTools](https://cran.r-project.org/bin/windows/Rtools/) (only needed for Windows computers)

If you are a new R user, we highly suggest doing the following before beginning to use these templates:

1.  Complete our 5 free online introductory [R tutorials](www.appliedepi.org/tutorial). To begin, create your free Applied Epi account, as explained in the linked page
2.  Review the [R basics chapter](https://epirhandbook.com/en/r-basics.html#r-basics) of our Epidemiologist R Handbook
3.  Review the [Reports with R Markdown chapter](https://epirhandbook.com/en/reports-with-r-markdown.html) of the Epidemiologist R Handbook

Once you have an understanding of how R and RStudio work, know how to run basic commands, have practiced running R code on your computer, and have read about R Markdown reports, then you are ready to use our {sitrep} R Markdown templates to produce situation reports.

### Install the {sitrep} R package

To access the templates, you must install the "sitrep" R package. R packages are often written in text with curly brackets, like this: {sitrep}

1.  Open RStudio
2.  Run the following command in the R Console pane:

```{r, eval=FALSE, echo=TRUE}
# install the {remotes} package, which is used to install the {sitrep} package
install.packages("remotes")
```

3.  Now run this command, also in the R Console pane. The {sitrep} package is currently available from a location on Github, a site to store and share code.

```{r, eval=FALSE, echo=TRUE}
# use {remotes} to install {sitrep} from Github
remotes::install_github("R4EPI/sitrep", dependencies = TRUE)
```

### Download the solution to follow along

We suggest that you [Download the "solution" R Markdown template](/public_data/Final_casestudy_AJS_template.Rmd) used in this vignette, in order to follow-along step-by-step.

### Open a new RStudio project

Working in R is best done from within an "RStudio project" that contains all the datasets, R scripts, and outputs for a single analytical project.

You should review the Epidemiologist R Handbook chapters on [RStudio projects](https://epirhandbook.com/en/r-projects.html). Below are brief instructions on how to create a new RStudio project:

> Select *File -\> New Project* from the top RStudio menu.

> Creating a new R project will create:

> -   A new project directory (folder)\
> -   A R project file (.Rproj) in the project folder that serves as a shortcut to open the project via RStudio

> Note:

> -   You can make sub-folders for data, scripts, etc.\
> -   All files produced by R scripts saved to the project will also save into the project folder\
> -   Navigate and open project files, such as scripts, using the Files tab of the lower-right RStudio pane.\
>     \*For an RMarkdown to produce the final report correctly, you must have read/write permissions for the location of the R project folder. On an MSF computer this is generally the c:/temp folder.

One significant benefit of working within an R project is that all associated files (datasets, etc.) are located within the same folder. Thus, to send the project to someone else, **simply zip the folder and send**. All files referenced in the scripts are included and if you utilized the {here} package, no file paths are unique to your computer.

## Open a {sitrep} R Markdown template

For this vignette, you will need to open the template for "AJS" (Acute Jaundice Syndrome). See the instructions and GIF below.

### Opening a template

1.  Ensure you have installed the {sitrep} package as detailed above\
2.  Go to the File menu, select "New File", and select "R Markdown"\
3.  In the pop-up window, select "From template"
4.  From the list of templates, select "AJS Outbreak Report" - the template will open in RStudio\
5.  Go to File, "Save", and write a name for this R Markdown script; save it into your RStudio Project folder

![This above GIF shows how to open one of the templates, after having installed the {sitrep} R package](images/opening_template.gif){width="549"} \### Folder structure

It is highly recommended that you store your R Markdown script and any data used in the same folder (the RStudio project). You should review the Epidemiologist R Handbook chapters on [RStudio projects](https://epirhandbook.com/en/r-projects.html).

## Become familiar with the template

### Header

The very top of the R Markdown template consists of a header surrounded by `---` lines. This is called the "YAML" header, in which you can edit the title of your document. The other settings in the header define the default document type produced (Microsoft Word) when the RMarkdown is "knit".

![](images/AJS_Intro.png)

### Introduction to the template

This text just below the header is for your reference only. If left alone it will appear at the beginning of your report, so you should replace it with a more relevant introduction. Any text you write here can be formatted with *italics*, **bold**, bullets, numbers, tables, and in-line code. See this [RMarkdown cheatsheet](/images/rmarkdown_cheatsheet_2.0.pdf) for tips on formatting.

### Explanatory text

Green text located between these special characters (`<!--` and `-->`), is for your information only. This text will **not** appear in the final Word document report.

## Installing and loading packages

```{r, eval=TRUE, echo=FALSE}
# Install and load required packages for this template
pacman::p_load(
  knitr,           # create output docs
  here,            # find your files
  rio,             # read in data
  forcats,         # handle ordinal columns
  lubridate,       # handle dates
  dplyr,           # clean/shape data
  tidyr,           # clean/shape data
  stringr,         # clean text
  ggplot2,         # create plots and charts
  slider,          # moving averages
  parsedate,       # guessing dates
  patchwork,       # combine multiple plots
  sitrep,          # MSF field epi functions
  epitabulate,     # helper functions
  janitor,         # clean data
  matchmaker,      # matching data
  gtsummary,       # tables
  flextable,
  tsibble,         # epiweeks
  sf,              # encode spatial vector data
  ggspatial,       # plot maps
  purrr,           # iterate over data
  ggrepel,         # text on maps
  DT               # interactive tables for this walkthrough
)
```


**The first code chunk, `setup`, is very important and must be kept.** When the final report is knit, this chunk will run without producing any output, but it is required for the template to work correctly.

This `setup` code chunk does the following:

-   Set default settings for tables such as size, error handling, and display of missing values

-   Run a command using the {pacman} package to install and load the R packages that are used in this template.

-   Set default text size in graphic plots

-   Establish visual theme settings for epidemic curve plots

### About the packages used in this template

#### Generic packages

| R package   | Use                                                                        |
|-------------|----------------------------------------------------------------------------|
| {knitr}     | to create the output document (pdf, doc, html)                             |
| {here}      | to locate files within the RStudio project                                 |
| {rio}       | to import data from various formats (csv, xlsx, tsv, etc.)                 |
| {dplyr}     | to clean and handle data                                                   |
| {forcats}   | to clean the data and create ordinal columns ("factors")                 |
| {stringr}   | to clean and handle text/characters (known as "strings")                   |
| {ggplot2}   | to visualise your data in plots                                            |
| {lubridate} | to handle dates and times                                                  |
| {gtsummary} | to produce summary tables                                                  |
| {flextable} | to format summary tables for presentation                                  |
| {purrr}     | to iterate processes - repeating the same process many times with small variations |

#### More epidemiology-specific packages

| R package            | Use                                                                |
|----------------------|--------------------------------------------------------------------|
| {sitrep}             | includes the templates and functions useful for field epidemiology |
| {slider}             | to calculate moving averages                                       |
| {parsedate}          | to interpret "messy" date formats                              |
| {matchmaker}         | to clean data using data dictionaries                                 |
| {tsibble}            | to conduct time-series analysis, and produce epidemiological weeks      |
| {sf} and {ggspatial} | to produce maps                                                 |


## Define the report settings  
Typically, near the top of a report you establish settings that will be used throughout the report.  

### Define the reporting date
This chunk continues the set-up of time-related parameters. There are potentially two items to edit in this chunk:  

1. Set the `reporting_date` object to the latest date for data to be used in your report. This value is used in later code chunks, for example when producing epidemic curves. Follow the format in the code: `ymd("YYYY-MM-DD")`

```{r, eval=TRUE, echo=TRUE}
reporting_date <- ymd("2018-04-23")
```

*If you want to use the current date (as per your computer) which will be automatically updated, place `Sys.Date()` within the `ymd()` function, like this: `ymd(Sys.Date())`. Or try, `ymd(Sys.Date()-3)`.*  


### Define the current week

2. Adjust whether your epidemiological weeks begin on Mondays, Sundays, etc. The second command in this code chunk creates the object `reporting_week` using the `report_date` object from the above command.  

*The only part of this command that you may need to change is the number near the end*. Use `1` for Mondays, or `7` for Sundays.  

```{r, eval=TRUE, echo=TRUE}
reporting_week <- tsibble::yearweek(reporting_date, 1)
```


### Set plot defaults  

In this next chunk, you define the labels that will be used in epidemic curve plots. This code should not need any adjustment. 

You may notice that the `subtitle = ` label uses the object `reporting_week` in curly braces - this means that the value of `reporting_week` defined in the previous chunk will show in its place.  

```{r, echo=FALSE, eval=TRUE}
ggplot2::theme_set(theme_classic(base_size = 18))

## sets the theme in ggplot for epicurves
epicurve_theme <- theme(
  axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
  legend.title = element_blank(),
  legend.position = "bottom",
  panel.grid.major.x = element_line(color = "grey60", linetype = 3),
  panel.grid.major.y = element_line(color = "grey60", linetype = 3))
```


```{r echo=TRUE}
## sets the labels in ggplot for the epicurves
epicurve_labels <- labs(x = "Calendar week", 
                        y = "Cases (n)", 
                        title = "Cases by week of onset",
                        subtitle = str_glue("Source: MSF data from {reporting_week}")
                       ) 
```






## Import data

This {sitrep} template offers three (3) dataset options:  

1. Generating a fake linelist for practice  
2. Using standardised DHIS2 data  
3. Using non-standard data

Ensure that you use the code chunks for your situation, and either delete or commend out (using \#) the code for the other scenarios.  

![](images/AJS_data_import.png){width="549"} 

### Import data for this vignette  

**This vignette will demonstrate using *non-standard data***. 

The AJS linelist spreadsheet is named "AJS_Am_Timan.xlsx", and is stored in the {sitrep} R package. We created the function `download_outbreak_linelist()` which will copy this file into a location of your choice on your computer.  

Run the commands below. Then look for a pop-up window and use it to select your RStudio project folder as the file destination (the pop-up may appear *behind* RStudio out-of-sight).  

```{r, echo=TRUE}
# load the sitrep package
library(sitrep)
```

```{r, eval= FALSE, echo=TRUE}
# run the command below as it is (empty parentheses).
# Look for a pop-up window, and select your RStudio project as the file destination
download_outbreak_linelist()


```
  

```{r, echo=FALSE}
# BACK-END load internal data for this vignette
linelist_raw <- AJS_linelist_internal
   
```

For this vignette, **comment out (using \#) or delete the `read_fake_data`, `read_DHIS_data`, and `prep_DHIS_data` chunks**, as they are not necessary and will cause errors.

A successful import should result in a dataframe object named `linelist_raw` in your RStudio Environment pane, consisting of 1447 observations and 122 columns.

![](images/linelist_raw_imported.png){width="549"}


### Import the population data into R

Now that the linelist data are saved in your RStudio project, you can import them into R using `import()`, as described in the Epidemiologist R Handbook's [Import and export](https://epirhandbook.com/en/import-and-export.html) chapter.

To import a specific sheet from an xlsx spreadsheet, use the `which = ` argument as shown below. You can also use the `na = ` argument to specify a value that should be converted to `NA` across the entire sheet (e.g. 99 or "Missing", or "" which represents an empty space). 

Th example command below imports an Excel spreadsheet called "AJS_AmTiman.xlsx" which is not password protected, and is stored in the RStudio project you are using for this analysis. The `which =` argument of `import()` specifies that the sheet "linelist" should be imported. You would edit this command as appropriate for your situation.

Note: If you saved the data into a subfolder of your RStudio project, you would need to also use the `here()` function, as described in the Epidemiologist R Handbook's [Import and export](https://epirhandbook.com/en/import-and-export.html) chapter.

```{r echo=TRUE, eval=FALSE}
# example code for importing an Excel file saved in RStudio project

# Import the linelist from your RStudio project (top folder)
linelist_raw <- import("AJS_AmTiman.xlsx", which = "linelist")

# Import from a "data" subfolder of the RStudio project
#linelist_raw <- import(here("data", "AJS_AmTiman.xlsx"), which = "linelist")
```














## Align columns with data dictionary  

As in this vignette, if your dataset is not initially aligned with the MSF data dictionary, the column names and values must be standardised. This is because later parts of the template expect certain column names and values.  

**This process can take some time** and requires looking back-and-forth between your dataset, RStudio, and the MSF data dictionary. In this vignette, here are the steps taken:  


### Review the MSF data dictionary

Uncomment and run this command, found in the `prep_nonDHIS_data` chunk, to view the MSF data dictionary for the specific disease ("AJS" in this example).

```{r eval = FALSE}
# Creates object linelist_dict using the msf_dict function from the sitrep package
linelist_dict <- msf_dict("AJS", compact = FALSE) %>%
  select(option_code, option_name, everything())
```

The dataframe `linelist_dict` should appear in your Environment pane. You can view the data dictionary by clicking on `linelist_dict` in the Environment pane, or running the command `View(linelist_dict)` (note capital V).  

![](images/AJS_Data_Dictionary.png){width="549"}

### Clean the column names

These steps standardize how your column names are written, such as changing spaces and dots to underscores ("_"). **Uncomment these lines of code in the `prep_nonDHIS_data` chunk**.  

**First, make a copy of the dataframe `linelist_raw` but with a new name: `linelist_cleaned`.**  
Throughout the template you will modify and improve this `linelist_cleaned` dataframe. However, you can always return to the `linelist_raw` version for reference.  

```{r, echo=TRUE}
## make a copy of your orginal dataset and name it linelist_cleaned
linelist_cleaned <- linelist_raw
```

**Second, use the `clean_names()` function** from the {janitor} package to fix any column names with non-standard syntax.

```{r, echo=TRUE, eval=TRUE}
linelist_cleaned <- janitor::clean_names(linelist_cleaned)
```



### Match your dataset's columns to MSF's AJS data dictionary  

Standardized column names are required for this template to work smoothly, and the column names in the Am Timan linelist do not align with the names this template expects.  

We need to change many column names, so we have developed a special function to help map our columns to the expected column names. The template also offers code that you can un-comment if there are only a few column names to change.  
This is done by un-commenting and running a command found in near the end of the chunk `prep_nonDHIS_data`. We run the command `msf_dict_rename_helper()` with "AJS" to get the correct helper for this template.

```{r eval = FALSE, warning = FALSE, message = FALSE}
msf_dict_rename_helper("AJS")
```

Running this command copies a block of code text to our clipboard. Now, *paste* from your clipboard into the template in the bottom of the code chunk. For the AJS template, the code that is pasted will look like the code below. This code uses the function `rename()` to change column names.  

*Be sure to paste the AJS rename helper code into an existing R code chunk (e.g. the bottom of `prep_nonDHIS_data`), not into the white space of the RMarkdown script.*

![](images/Variable_Helper_Paste_Revised.gif){width="549"}

#### Complete the mapping of column names

The columns will be in alphabetical order, sorted by if the appear in the template (REQUIRED) or not (optional).

Some columns---such as lab results or symptoms---are listed as REQUIRED, but are not absolutely necessary for the template to run. If you do not have a column and are unsure if it's necessary, do a search in the template to see if it's marked as OPTIONAL in the code chunk.

To the right of each equals sign, and before each comma, type the *exact* names of the columns from your dataset that correspond to the expected MSF data dictionary columns on the left. The last un-commented line should not have a comma at the end. 

Be sure to comment (#, ctrl + shift + c) out the lines of each data dictionary column in the helper that is not present in your dataset!


![](images/AJS_Variable_Names.gif){width="549"}

If you see this error:  
`Error in is_symbol(expr) : argument "expr" is missing, with no default`  
then you likely forgot to comment (#) the line for a column you did not use.

And now you can see the result:  

![](images/AJS_Rename_Variables_Proof.png){width="549"}


Below is the code for the column assignments used in this case study walk-through. You can copy or reference this when building your template for this vignette.  

```{r, echo=TRUE}
## Below are the column assignments used in this case study example

linelist_cleaned <- rename(linelist_cleaned,
  # age_days                          =   , # INTEGER_POSITIVE  (Dataset only has years)
  # age_months                        =   , # INTEGER_POSITIVE  (Dataset only has years)
  age_years                         = age, # INTEGER_POSITIVE  (REQUIRED)
  bleeding                          = medbleeding, # BOOLEAN           (REQUIRED)
  case_number                       = hevid, # TEXT              (REQUIRED)
  chikungunya_onyongnyong           = chik, # TEXT              (REQUIRED)
  # convulsions                       =   , # BOOLEAN           (Not in dataset)
  date_of_consultation_admission    = dtmedical, # DATE              (REQUIRED)
  # date_of_exit                      =   , # DATE              (Not in dataset)
  date_of_onset                     = dtjaundice, # DATE              (REQUIRED)
  dengue                            = dengue, # TEXT              (REQUIRED)
  # dengue_rdt                        =   , # TEXT              (Not in dataset)
  diarrhoea                         = meddiar, # BOOLEAN           (REQUIRED)
  # ebola_marburg                     =   , # TEXT              (Not in dataset)
  epigastric_pain_heartburn         = medepigastric, # BOOLEAN           (REQUIRED)
  exit_status                       = outcomehp, # TEXT              (REQUIRED)
  fever                             = medfever, # BOOLEAN           (REQUIRED)
  generalized_itch                  = meditching, # BOOLEAN           (REQUIRED)
  headache                          = medheadache, # BOOLEAN           (REQUIRED)
  hep_b_rdt                         = medhepb, # TEXT              (REQUIRED)
  hep_c_rdt                         = medhepc, # TEXT              (REQUIRED)
  hep_e_rdt                         = medhevrdt, # TEXT              (REQUIRED)
  # history_of_fever                  =   , # BOOLEAN           (Not in dataset)
  joint_pains                       = medarthralgia, # BOOLEAN           (REQUIRED)
  # lassa_fever                       =   , # TEXT              (Not in dataset)
  # malaria_blood_film                =   , # TEXT              (Not in dataset)
  malaria_rdt_at_admission          = medmalrdt, # TEXT              (REQUIRED)
  mental_state                      = medmental, # BOOLEAN           (REQUIRED)
  nausea_anorexia                   = mednausea, # BOOLEAN           (REQUIRED)
  other_arthropod_transmitted_virus = arbovpcr, # TEXT              (REQUIRED)
  other_cases_in_hh                 = medothhhajs, # BOOLEAN           (REQUIRED)
  # other_pathogen                    =   , # TEXT              (Not in dataset)
  other_symptoms                    = medother, # BOOLEAN           (REQUIRED)
  patient_facility_type             = hospitalised, # TEXT              (REQUIRED)
  patient_origin                    = quartier, # ORG UNIT          (REQUIRED)
  sex                               = sex, # TEXT              (REQUIRED)
  test_hepatitis_a                  = medhavelisa, # TEXT          (REQUIRED)
  test_hepatitis_b                  = medhbvelisa, # TEXT          (REQUIRED)
  test_hepatitis_c                  = medhcvelisa, # TEXT          (REQUIRED)
  test_hepatitis_e_genotype         = hevgenotype, # TEXT          (REQUIRED)
  # test_hepatitis_e_igg              =   , # TEXT              (In same column as elisa result)
  test_hepatitis_e_igm              = hevrecent, # TEXT          (REQUIRED)
  test_hepatitis_e_virus            = medhevelisa, # TEXT          (REQUIRED)
  # time_to_death                     =   , # TEXT              (Not in dataset)
  # typhoid                           =   , # TEXT              (Not in dataset)
  vomiting                          = medvomit, # BOOLEAN           (REQUIRED)
  # water_source                      =   , # TEXT              !(Split across many columns)
  yellow_fever                      = yf, # TEXT              (REQUIRED)
  # arrival_date_in_area_if_3m        =   , # DATE              (optional)
  date_lab_sample_taken             = medblooddt, # DATE              (optional)
  # date_of_last_vaccination          =   , # DATE              (optional)
  # delivery_event                    =   , # TRUE_ONLY         (optional)
  detected_by                       = referredby, # TEXT              (optional)
  # event_file_type                   =   , # TEXT              (optional)
  # ever_received_hepatitis_e_vaccine =   , # TEXT              (optional)
  # foetus_alive_at_admission         =   , # TEXT              (optional)
  # lab_comments                      =   , # TEXT              (optional)
  # msf_involvement                   =   , # TEXT              (optional)
  patient_origin_free_text          = block, # TEXT              (optional)
  pregnancy_outcome_at_exit         = medppoutcome, # TEXT              (optional)
  pregnant                          = medpreg, # TEXT              (optional)
  # recent_travel                     =   , # BOOLEAN           (optional)
  # residential_status_brief          =   , # TEXT              (optional)
  # traditional_medicine_details      =   , # TEXT              (optional)
  # traditional_medicines             =   , # BOOLEAN           (optional)
  treatment_facility_site           = hpid, # TEXT              (optional)
  # treatment_location                =   , # ORGANISATION_UNIT (optional)
  trimester                         = medpregtri # TEXT              (optional)
)
```













## Provide population counts  

The population estimates for regions and age groups are collected in order to calculate incidence, attack rates, and mortality rates. **The template allows several ways to provide population estimates, but importantly, the data must be stored under specific column names.** In your real-use case, the specific age groupings and region names in your data may differ from the examples, but the column names must follow the expected standard.  

You can do one of the following:  

A) Import a spreadsheet file with the estimates in the correct format  
B) Use the function `gen_population()` to derive the estimates from proportions based on globally accepted standards  
C) Use the function `gen_population()` to directly enter counts  

![](images/AJS_Population_Counts.png){width="549"}

**For this case study, we use a spreadsheet containing population data by quartier (stored inside the {sitrep} R package).**  

### For this vignette

Because the vignette data are stored within the {sitrep} package, you can copy them to your RStudio project using the below commands. Look for the pop-up window (it may appear behind RStudio) and assign the destination as your RStudio project.  

```{r, echo=TRUE, eval=FALSE}
library(sitrep)
download_outbreak_pop()
```


```{r, echo=FALSE}
# BACK-END load internal data for vignette use
population_data_region <- AJS_pop_internal
```

**Note the name change** from the default `population_data`, because the data we are importing is for regions (quartiers).  

Note also that the columns names in the spreadsheet have already been edited to match the expected column names ("patient_origin" and "population").



### Import the population data into R

Now that the population data are saved in your RStudio project, you can import them into R using `import()`, as described in the Epidemiologist R Handbook's [Import and export](https://epirhandbook.com/en/import-and-export.html) chapter.

To import a specific sheet from an xlsx spreadsheet, use the `which = ` argument as shown below. You can also use the `na = ` argument to specify a value that should be converted to `NA` across the entire sheet (e.g. 99 or "Missing", or "" which represents an empty space). 

Th example command below imports an Excel spreadsheet called "AJS_AmTiman_population.xlsx" which is not password protected, and is stored in the RStudio project you are using for this analysis. The `which =` argument of `import()` specifies that the sheet "quartier" should be imported. You would edit this command as appropriate for your situation.

Note: If you saved the data into a subfolder of your RStudio project, you would need to also use the `here()` function, as described in the Epidemiologist R Handbook's [Import and export](https://epirhandbook.com/en/import-and-export.html) chapter.

```{r echo=TRUE, eval=FALSE}
# example code for importing an Excel file saved in RStudio project

# Import dataset from your RStudio project (top folder)
population_data_region <- import("AJS_AmTiman_population.xlsx", which = "quartier")

# Import from a "data" subfolder of the RStudio project
#population_data_region <- import(here("data", "AJS_AmTiman_population.xlsx"), which = "quartier")
```

Before moving on, delete or comment (#) all the code lines in the chunk that create fake data or proportions for `population_data_region`, `population_data_age`, and `population_data_age_months`.  

**If you choose not to complete the population estimates section at all, delete or comment (#) all code in this section, and expect that later sections on attack and mortality rates will not produce output.**




## Provide laboratory testing data  

In this example, laboratory testing data is already included in the linelist dataset. Therefore, **there is no need to import a separate laboratory testing dataset**.  

In your personal use-case, if you have a separate laboratory dataset, use the alterate code provided in the template to import and then join that dataset to the linelist dataset. You can join them based on an ID or case number column that is present in both datasets. For more information about joins, see the Epidemiologist R Handbook chapter on [Joining data](https://epirhandbook.com/en/joining-data.html).  





## Browse data

This code chunk can be used to explore your data, but **remember to eventually comment out (#) these code lines if you do not want these outputs in your final report**.

View the first 10 rows of the data
```{r echo=TRUE, eval=FALSE}
head(linelist_cleaned, n = 10)
```

Wiew your whole dataset interactively. Remember that `View` needs to be written with a capital *V*. 

```{r echo=TRUE, eval=FALSE}
View(linelist_cleaned)
```

Overview of column types and contents
```{r echo=TRUE, eval=FALSE}
str(linelist_cleaned)
```

Gives mean, median and max values of columns, gives counts for categorical columns, and also gives number of NAs.

```{r echo=TRUE, eval=FALSE}
summary(linelist_cleaned)
```

View unique values contained in columns 

```{r echo=TRUE, eval=FALSE}
unique(linelist_cleaned$sex)
```








## Standardize and clean data  

Because the Am Timan linelist is not yet aligned with the MSF data dictionary and what the template expects, there are several data cleaning steps we must complete.  


### Re-code missing values from periods (.) to `NA` at import  

Sometimes, missing values in the dataset are represented by a period (.), or by a number such as 99. This causes problems because R expects missing values to be coded as `NA`. For example, if missing dates are coded as a period (.), then R gives an error because the periods are considered character values that do not fit the expected YYYY-MM-DD date format.  

If you encounter a problem like this, modify the general import command in the `read_nonDHIS_data` chunk of the template to include the argument `na = ` within the `import()` command. Writing `na = "."` specifies a period as the value in the Excel sheet that R will consider to be "missing". Alternatively, your scenario may require you to write `na = 99` (no quotation marks, as 99 is a number). As the data are re-imported, all cells with that value are now changed to `NA`.  

```{r eval=FALSE, echo=TRUE}
# An example of a revised import command (don't forget the comma between arguments!)
linelist_raw <- import("AJS_AmTiman.xlsx", which = "linelist", na = ".")
```




### Drop ineligible observations  

The next code chunk in the script, `remove_unused_data`, drops observations with missing `case_number` or `date_of_consultation_admission`. Because the Am Timan dataset contains observations for patients seen in the community and at the hospital, these criteria may not be relevant.  

We can check, using the code below, and see that of the 1447 observations, there are 0 missing `case_number` or `dateentry`, but 616 observations missing `date_of_consultation_admission`. These 616 community-identified cases are still of interest for our report, so we will not drop them.  


```{r echo=TRUE}
# Print the number of observations missing case_number, dateentry, and date_of_consultation_admission
nrow(linelist_cleaned %>% filter(is.na(case_number))) 

nrow(linelist_cleaned %>% filter(is.na(dateentry))) 

nrow(linelist_cleaned %>% filter(is.na(date_of_consultation_admission))) 
```

Notice that the commands above returns the number of rows missing `date_of_consultation_admission`, but does not overwrite `linelist_cleaned` with anything, as there is not an assignment operator `<-`. 

**After viewing the above results, we must comment out (#) the two lines that `filter()` the `linelist_cleaned` dataset - we do *not* want to remove those 616 observations.**

```{r, echo=TRUE}

# linelist_cleaned <- linelist_cleaned %>% 
#   filter(!is.na(case_number) & !is.na(date_of_consultation_admission)) 

```




### Standardise dates  

In `linelist_cleaned` there are seven columns containing dates (`dateentry`, `epiquestdate`, `date_of_consultation_admission`, `date_of_onset`, `date_lab_sample_taken`, `dthospitalisation`, `dtdeath`).  

Because the Am Timan dataset contains date columns not found in the data dictionary, we will **comment out the first code lines related to DHIS-aligned data** (those commands convert only date columns known to the MSF data dictionary).  

Instead, we uncomment and **use the commands applicable to non-DHIS2 data**. You may see messages warning that a few date entries are not in the correct timeframe - this is ok and these entries will be addressed in the next step.  

The command below uses the function `parse_date()` from the function {parsedate} to "guess" each individual date in the column - this approach can handle very messy data in which some dates are written as "17 October 2016" and others may be written as "2016-10-27". Read more about the `parse_date()` function [here](https://readr.tidyverse.org/reference/parse_datetime.html).  


```{r results='hide', echo=TRUE}
## Non-DHIS2 data --------------------------------------------------------------
## Use this section if you did not have DHIS2 data. 

## use the parse_date() function to make a first pass at date columns.
linelist_cleaned <- linelist_cleaned %>%
  mutate(
    across(.cols = matches("date|Date"),
           .fns  = ~ymd(parsedate::parse_date(.x))))
```

The code above successfully converted many of our date columns (remember you can check a column's class like `class(linelist$date_of_onset`). 

```{r, echo=TRUE}
class(linelist_cleaned$date_of_onset)
```


However, this command did not detect the date columns named with "dt", such as `dtdeath` and `dthospitalization`. We must **add code to convert those columns manually**:  

```{r, echo=TRUE}
# Individually convert column using the parse_date() function
linelist_cleaned <- linelist_cleaned %>% 
  mutate(dtdeath = ymd(parse_date(dtdeath)),
         dthospitalisation = ymd(parse_date(dthospitalisation)))
```


Next we uncomment and **modify code to correct unrealistic dates**. We have browsed our data and know that there are observations with `date_of_onset` outside the reasonable range:


```{r, echo=TRUE}
# Check range of date_of_onset values, ignoring (removing) missing values
range(linelist_cleaned$date_of_onset, na.rm = TRUE)
```

We convert dates outside the expected range (April 2016 to October 2017) to missing using `case_when()`. Note that when making the assignment on the right-hand side (RHS), wrap `NA` in `as.Date()` because all right-hand side values must be the same class ("Date" in this scenario).  

As you modify this chunk for your own situation, remember that **each left-hand side (LHS) of the `~` must be a logical statement** (not just a value), and to include commas at the end of each `case_when()` line (except the last one). In addition, it is best to write `case_when()` lines from most specific at top to most general/catch-all at the bottom. You can read more about `case_when()` [on this tidyverse site](https://dplyr.tidyverse.org/reference/case_when.html).  



```{r, echo=TRUE}
# Convert dates before April 2016 or after October 2017 to missing (NA)
linelist_cleaned <- linelist_cleaned %>% 
  mutate(date_of_onset = case_when(
    date_of_onset < as.Date("2016-04-01")  ~ as.Date(NA),
    date_of_onset > as.Date("2017-10-01")  ~ as.Date(NA),
    TRUE                                   ~ date_of_onset))
                           
```

Check the range of dates again:  


```{r, echo=TRUE}
# Check range of date_of_onset values, ignoring (removing) missing values
range(linelist_cleaned$date_of_onset, na.rm = TRUE)
```

We also must **use the provided code to create a column called `epiweek`**. Although there are already columns that give the epidemiological weeks of various data points, it is safer to calculate a new column ourselves, AND later code chunks rely on the column being named `epiweek`.

```{r, echo=TRUE}
# Create column epiweek
linelist_cleaned <- linelist_cleaned %>% 
  mutate(
    
    ## create an epiweek column 
    epiweek = tsibble::yearweek(
      date_of_onset,
      week_start = 1),
    
    ## create a date version of epiweek 
    epiweek_date = as.Date(epiweek))  
```



### Create age groups  

**To create a column for `age_group` we must first clean the column `age_years`**. If we look at the range of values in `age_years`, we see something strange:  

```{r, echo=TRUE}
# See the range of age_years values, removing (excluding) NA
range(linelist_cleaned$age_years, na.rm = TRUE)
```

We know there are ages older than 9 years. So we check `class(linelist_cleaned$age_years)` and see that R is reading this column as class *character*, not numeric!  

```{r, echo=TRUE}
class(linelist_cleaned$age_years)
```


We must convert it by **adding the following command** to the script:  

```{r, echo=TRUE}
# Convert column age_years to numeric class
linelist_cleaned <- linelist_cleaned %>% 
  mutate(age_years = as.numeric(age_years))
```

If we run the `range()` command again, we can see that the corrected range is 0.08 to 75.  

```{r, echo=TRUE}
# See the range of age_years values, removing (excluding) NA
range(linelist_cleaned$age_years, na.rm = TRUE)
```

Then, the chunk contains three age-group commands marked "OPTIONAL":  

* We do **not** want to use "add under 2 years to the age_years column" because it assumes that we already have a column `age_months`, which we do not. **Comment out (#) this code.**  

* We also do **not** need to use "change those who are above or below a certain age to NA", because we already know our range of ages and do not have any outside an expected range. **Comment out (#) this code.**  

* We **do** want to use "create an age_months column from decimal years column", as we do have decimal years. This command will **create an age_months column** that has a value if a patient is under 5 years. **Uncomment and use this code.**  



```{r, echo=TRUE}
# OPTIONAL: create an age_months column from decimal years column
linelist_cleaned <- linelist_cleaned %>%
  mutate(age_months = case_when(
            age_years < 5 ~ age_years * 12))
```


**Now we can create the columns `age_group_mon` and `age_group`.**

First we use the code to create an `age_group_mon` column for children under 5 years based on `age_months`.  

```{r, echo=TRUE}
## create age group column for under 2 years based on months
linelist_cleaned <- linelist_cleaned %>% 
  mutate(age_group_mon = age_categories(
    age_months, 
    breakers = c(0, 6, 9, 12, 23), 
    ceiling = TRUE))
```

...and we use the similar code below to create `age_group` as groupings of the column `age_years`. The function `age_categories()` is used with its `breakers =` argument set to a vector `c( )` of numbers: 0, 3, 15, 30, and 45. Thus, a 30-year old patient will be included in an age group named "30-44". Read more about `age_categories()` by searching the Help pane.

```{r, echo=TRUE}
## create an age group column by specifying categorical breaks
linelist_cleaned <- linelist_cleaned %>% 
  mutate(age_group = age_categories(
    age_years,
    breakers = c(0, 2, 15, 30, 45)))
```

Finally, ensure that the remaining code in the `create_age_group` chunk is commented out (#).  



### Create and clean outcome columns  

In the next chunk (`create_vars`) we can **comment out/ignore the code that makes a new numeric column `obs_days`**, because we do not have a `date_of_exit` column in our dataset. 

While the template directs us to next create a `DIED` column based on the column `exit_status` containing the characters "Dead", we must first **clean our `exit_status` column, which is currently in French**. 

**Add code that uses `case_when()` to assign new values in a new `exit_status2` column**, as in the code below. We do this by **modifing code from just below in the template.** The code uses the function `case_when()` to re-define the dataframe `linelist_cleaned` as itself, but mutated to create the new column `exit_status2`. The values in `exit_status2` are based on the values in `exit_status`, such that when `exist_status == "Décédé"`, the value in `exit_status2` is "Dead", and so on. 

As you modify this chunk for your own situation, remember that **each left-hand side (LHS) of the `~` must be a logical statement** (not just a value), and to include commas at the end of each `case_when()` line (except the last one). It is best to write `case_when()` lines from most specific at top to most general/catch-all at bottom. You can read more about `case_when()` [at this tidyverse site](https://dplyr.tidyverse.org/reference/case_when.html).  

Another tip is that all the right-hand side (RHS) values must be the same class (either character, numeric, etc.). So, if your other RHS values are character and you want one particular RHS value to be missing, you cannot just write `NA` on the RHS. Instead you must use the special character version of `NA`, `NA_character_`. This may seem like unnecessary trouble, but it ensures consistency among your data in the long term.  

```{r, echo=TRUE}
# Tabulate to see all possible values of exit_status
linelist_cleaned %>% tabyl(exit_status, show_na = FALSE)

# Create exit_status2 from values in exit_status
linelist_cleaned <- linelist_cleaned %>% 
  mutate(exit_status2 = case_when(
    exit_status == "Décédé"          ~ "Dead",
    exit_status == "-"               ~ NA_character_,
    exit_status == "Echappé"         ~ "Left",
    exit_status == "Déchargé/Guéri"  ~ "Discharged"
  ))

# Tabulate the NEW exit_status2 column to check correct assignment
# Tabulate to see all possible values of exit_status
linelist_cleaned %>% tabyl(exit_status2, show_na = FALSE)
```


**Now we can make the `DIED` column**, which is referenced in later code chunks. 

This command creates `DIED` as a logical (`TRUE` or `FALSE`) column, depending on whether each observation meets the criteria to the right of the assignment operator `<-`. Read more about the `%in%` operator on the [R Basics chapter](https://epirhandbook.com/en/r-basics.html#in) of the Epi R Handbook 

However, we must **modify the existing command to detect within the NEW column `exit_status2`, not `exit_status`.**  


```{r, echo=TRUE}
## Note we are directing R to look within the NEW exit_status2 column
linelist_cleaned <- linelist_cleaned %>% 
  mutate(DIED = str_detect(exit_status2, "Dead"))
```


### **Re-code values in `patient_facility_type`**  

When we assigned our columns to match the data dictionary, we used the column `hospitalisation` as the column `patient_facility_type`. **However, the values in that column do not match those expected by the template.** In the data dictionary, `patient_facility_type` should have values of "Inpatient" or "Outpatient." Currently, the values are "Oui" and "Non". In later code chunks, analyses are restricted to observations where `patient_facility_type == "Inpatient"`, thus, we must align the values to match the data dictionary.  

**To clean these values we uncomment and modify code also from the `create_vars` chunk, found under the heading "recode character columns" (the template example is of fixing incorrect spellings of village names).**  

Included below are extra steps before and after the `case_when()` command using `tabyl()` to verify the correct conversion of values. Remove these two commands after verifying, as otherwise their outputs will appear in the report.  


```{r, echo=TRUE}
# View all the values in patient_facility_type
linelist_cleaned %>% tabyl(patient_facility_type)

# Convert the values
linelist_cleaned <- linelist_cleaned %>% 
  mutate(patient_facility_type = case_when(
    patient_facility_type == "Oui" ~ "Inpatient",
    patient_facility_type == "Non" ~ "Outpatient"
  ))

# Re-check that the values converted sucessfully
linelist_cleaned %>% tabyl(patient_facility_type)

```



### Recode values in lab test columns  

For later steps, we will need many of our testing columns to have the values "Positive" and "Negative" instead of "Pos" and "Neg". We can uncomment and use this helper code to make that change. Some notes about this step:  

* You may receive warnings about unknown levels in some columns - this is generally okay but check the changes visually to be sure with `View(linelist_cleaned)`.  
* One column (`test_hepatitis_e_igm`) has values 0 and 1 that we want to change to "Positive" and "Negative". The `fct_recode()` function expects *character* RHS inputs (within quotes) - so put "0" and "1" in quotes, as below, and confirm the change has occurred as expected.  
* One column (`test_hepatitis_b`) has a value recorded as "Weakly pos". For the purposes of this exercise we re-code this to "Positive".  


```{r, echo=TRUE}
## sometimes, coding is inconsistent across columns -- for example, "Yes" / "No"
## may be coded as Y, y, yes, 1 / N, n, no, 0. You can change them all at once!
## Create a list of the columns you want to change, and run the following.
## You may need to edit this code if options are different in your data.

# # create list of columns
change_test_vars <- c("hep_e_rdt", "hep_c_rdt", "hep_b_rdt", "test_hepatitis_a", "test_hepatitis_b", "test_hepatitis_c", "malaria_rdt_at_admission", "test_hepatitis_e_genotype", "test_hepatitis_e_igm", "test_hepatitis_e_virus", "hevpcr", "other_arthropod_transmitted_virus")

# standardize values
linelist_cleaned <- linelist_cleaned %>%
  mutate(
    across(
      .cols = all_of(change_test_vars),
      .fns = ~forcats::fct_recode(
        .x,
        Positive = "Pos",
        Positive = "pos",
        Positive = "yes",
        Positive = "Yes",
        Positive = "Weakly pos",
        Positive = "1",
        Negative = "Neg",
        Negative = "neg",
        Negative = "no",
        Negative = "No",
        Negative = "0")))
``` 

**Lastly for the chunk `create_vars` we must create a case definition column.** In this use of `case_when()`, the last line left-hand side (LHS) is TRUE, which serves as a catch-all for any other possible values that have not yet met the criteria of the earlier `case_when()` lines.  

Note how this code checks the column `hep_e_rdt` for "Positive". The earlier cleaning steps converting `hep_e_rdt` values from "Pos" to "Positive" were necessary for this case definition to properly apply.  

In addition we need to **change the value looked for in `other_cases_in_hh` from "Yes" to "Oui".**

```{r, echo=TRUE}
# You MUST modify this section to match your case definition. The below
# uses positive RDT for Confirmed and epi link only for Probable.

linelist_cleaned <- linelist_cleaned %>%
  mutate(case_def = case_when(
    is.na(hep_e_rdt) & is.na(other_cases_in_hh)           ~ NA_character_,
    hep_e_rdt == "Positive"                               ~ "Confirmed",
    hep_e_rdt != "Positive" & other_cases_in_hh == "Oui"  ~ "Probable",
    TRUE                                                  ~ "Suspected"
  ))
```



### Factor columns  

Factor columns require specialized functions to manipulate (see the Epidemiologist R Handbook chapter [Factors](https://epirhandbook.com/en/factors.html)). The one factor column we should clean now is `sex`, with the following code:  

```{r, echo=TRUE}
linelist_cleaned <- linelist_cleaned %>% 
  mutate(sex = fct_recode(sex, NULL = "Unknown/unspecified"))
```                                   
  
And we can verify the unique values in column `sex` with a quick table:

```{r, echo=TRUE}
linelist_cleaned %>% 
  tabyl(sex)
```

                                   
**Comment out (#) the rest of the `factor_vars` chunk**, such as the code to change the order of levels in categorical columns. The Am Timan dataset does not include a column `time_to_death` and we do not need to change the order of any categorical columns.  



### Cleaning patient origin  

We must **add one additional cleaning step** necessary for this dataset.  

The column `patient_origin` (originally the column *quartier* in the raw dataset) is referenced in the place analyses chunks, for example when calculating attack rates by region. In those steps, the column `patient_origin` in `linelist_cleaned` is *joined to* the column `patient_origin` of the data frame `population_data_region` (which was imported in the population and lab data steps).  

However, the place names in `population_data_region` are ALL CAPITAL LETTERS. This is not true in `linelist_cleaned` - so the join will not register any matches. The easiest way to solve this problem is to **redefine the `linelist_cleaned` column into all capital letters as well**, using the base R function `toupper()`, which means "to upper case":  

```{r, echo=TRUE}
linelist_cleaned <- linelist_cleaned %>% 
  mutate(patient_origin = toupper(patient_origin))
```


### Gather common columns for later analysis

These columns names are stored in vectors that are created using the function `c()`. These vectors of names will be referenced in later code chunks. This way if you want to run the same function over these columns you can simply use the named vector rather than typing out each column individually.  

The default template code creates two vectors - one for symptom columns (`SYMPTOMS`) and one for laboratory testing columns (`LABS`). 

```{r}
# vectors of column names ----------------------------------------------------
# create a grouping of all symptoms 

SYMPTOMS <- c("generalized_itch", 
#             "history_of_fever", 
              "fever",
              "joint_pains",
              "epigastric_pain_heartburn",
              "nausea_anorexia",
              "vomiting", 
              "diarrhoea",
              "bleeding", 
              "headache",
              "mental_state",
#              "convulsions",
              "other_symptoms"  
              )
```


```{r}
# create a grouping of all lab tests 
LABS <- c("hep_b_rdt", 
          "hep_c_rdt",
          "hep_e_rdt",
          "test_hepatitis_a",
          "test_hepatitis_b",
          "test_hepatitis_c",
#          "test_hepatitis_e_igg",
          "test_hepatitis_e_igm" ,
          "test_hepatitis_e_genotype",
          "test_hepatitis_e_virus",
          "malaria_rdt_at_admission",
#          "malaria_blood_film", 
          "dengue",
#          "dengue_rdt", 
          "yellow_fever",
#          "typhoid",
          "chikungunya_onyongnyong", 
#          "ebola_marburg",
#          "lassa_fever",
          "other_arthropod_transmitted_virus" 
#          "other_pathogen"
          )

```







### Final setup for report  

The next code chunk, `report_setup`, defines important date-related parameters used to produce the report, and filters to remove any observations with date of onset after the `reporting_week` (`reporting_week` is an object defined at the beginning of the script).  

The second command in this chunk (shown below) uses `filter()` to remove observations later than the end of the pre-defined `reporting_week`. However, `filter()` also removes observations with missing `date_of_consultation_admission`. To avoid this **we must can add ` | is.na(date_of_onset)` into the filter** (This is read as: "keep observations where `date_of_consultation_admission` is less than or equal to `obs_end`, OR, if `date_of_consultation_admission` is missing").  


```{r echo=FALSE}
# return the last day of the reporting week
obs_end   <- as.Date((reporting_week + 1)) - 1

# filter out cases after end of reporting week
linelist_cleaned <- linelist_cleaned %>% 
  filter(epiweek <= reporting_week)

# define the first week of outbreak (date of first case)
first_week <- min(linelist_cleaned$epiweek)


# outbreak start 
# return the first day in the week of first case 
obs_start <- as.Date(first_week)

## pull the unique weeks which occur
all_weeks <- seq(first_week,
                 reporting_week, by = 1) 

## change weeks to dates
all_weeks_date <-  as.Date(all_weeks) 
```



### Export if desired  

And finally, if desired you can **export the cleaned dataset** for other purposes.  

The command `export()` is also from the package {rio}.  

First, the object that you want to export is named (`linelist_cleaned`).  

Then, the function `str_glue()` creates the name for the new file by combining "AmTiman_linelist_cleaned_" with the current date and finally the extention ".xlsx". Don't be confused - this is equivalent of writing `rio::export(linelist_cleaned, "AmTiman_linelist_cleaned_2019-08-24.xlsx")`, but using `str_glue()` and `Sys.Date()` means the current date will always be used automatically.

```{r save_cleaned_data, eval = FALSE}
## OPTIONAL: save your cleaned dataset! 
## put the current date in the name so you know!
rio::export(linelist_cleaned, str_glue("AmTiman_linelist_cleaned_{Sys.Date()}.xlsx"))
```

















# Analysis

After cleaning the data, there are three sections of the template that produce tabular, graphic, and cartographic (map) analysis outputs.  



## Person

This first section produces descriptive analyses about the patient's demographics and attack rates.  





### Demographic Tables:

The person analysis section begins with a couple sentences that contain **in-line code** - code embedded in normal RMarkdown text (not within an R code chunk). The second sentence **in-line code** inserts the number of males and females by counting the observations with "Male" and "Female" in the variable `sex`. **Because our dataset's variable `sex` contains `M` and `F` instead, we must modify this in-line code (or modify our variable) so the `fmt_count()` function is searching for the correct terms.**  


```
From the start of the outbreak up until `r reporting_week` there were a 
total of `r nrow(linelist_cleaned)` cases. There were
`r fmt_count(linelist_cleaned, sex == "F")` females and
`r fmt_count(linelist_cleaned, sex == "M")` males. 

```

The first demographic table presents patients by their age group (the table's rows) and their relationship with the case definition (the table's columns). This chunk uses pipes to link many different functions together and produce a table (`%>%` is the pipe operator - see the [R Basics chapter](https://epirhandbook.com/en/r-basics.html#piping) of the Epidemiologist R Handbook:  

* `select()` keeps only the columns of interest
* `mutate()` converts any instances of `NA` in the `case_def` column to the word "Missing", so the table's presentation is more clean. Note this change only applies in this table, not in the underlying dataset.  
* `tbl_summary()` creates a summary table by age group with counts and percentages, with columns for each `case_def` value.  
* `add_overall()` is a function from {gtsummary} that adds overall totals to the table  
* `bold_labels()` is from {gtsummary} and makes the variable names bold.  
* `as_flex_table()` is from the {flextable} package, which allows more detailed editing of the table presentation  
* `bold()` makes the header bold  
* `set_table_properties()` aligns the table to the maximum width of the word document report.  




```{r describe_by_age_group_and_def, warning = FALSE, message = FALSE}
linelist_cleaned %>% 
  ## only keep variables of interest
  select(age_group, case_def) %>% 
  ## make NAs show up as "Missing" (so they are included in counts)
  mutate(case_def = fct_explicit_na(case_def, "Missing")) %>%
  ## create a table with counts/percentages by column
  tbl_summary(by = case_def,
              percent = "column",
              label = list(age_group = "Age group")) %>% 
  ## add totals
  add_overall() %>% 
  ## make variable names bold 
  bold_labels() %>% 
  # change to flextable format
  as_flex_table() %>%
  # make header text bold (using {flextable})
  bold(part = "header") %>% 
  # make your table fit to the maximum width of the word document
  set_table_properties(layout = "autofit")
```



You can change the argument `by =` to refer to any number of columns (in this next example, by `sex`).  

```{r total_props_agegroup_sex, warning = FALSE, message = FALSE}

linelist_cleaned %>% 
  select(age_group, sex) %>% 
  mutate(sex = fct_explicit_na(sex, "Missing")) %>%
  tbl_summary(by = sex,
              percent = "column",
              label = list(age_group = "Age group")) %>% 
  add_overall() %>% 
  ## make variable names bold 
  bold_labels() %>% 
  # change to flextable format
  as_flex_table() %>%
  # make header text bold (using {flextable})
  bold(part = "header") %>% 
  # make your table fit to the maximum width of the word document
  set_table_properties(layout = "autofit")
```

Outside of the R code chunk, note the **in-line code** that gives the counts and percents of rows missing sex and age_group.  




### Age pyramids 


To print an age pyramid in your report, use the code below. A few things to note:  

* The **column for `split_by =` should have two non-missing value options** (e.g. Male or Female, Oui or Non, etc. Three will create a messy output.)  
* The variable names work with or without quotation marks  
* The dashed lines in the bars are the midpoint of the un-stratified age group   
* You can adjust the position of the legend by replacing `legend.position = "bottom"` with "top", "left", or "right"  
* Read more by searching "age_pyramid" in the Help tab of the lower-right RStudio pane  

You can make this a pyramid of months by supplying `age_group_mon` to the `age_group =` argument.  



```{r age_pyramid, fig.width=6}
# plot age pyramid by sex 
age_pyramid(linelist_cleaned,
                  age_group = "age_group", 
                  split_by = "sex") + 
  labs(y = "Cases (n)", x = "Age group") + # change axis  labels
  theme(legend.position = "bottom",     # move legend to bottom
        legend.title = element_blank(), # remove title
        text = element_text(size = 18)  # change text size
       )
```

**To have an age pyramid of patients under 2 by month age groups, it is best to add a `filter()` step to the beginning of the code chunk, as shown below**. This selects for `linelist_cleaned` observations that meet the specified critera and passes that subsetted dataset through the "pipes" to `age_pyramid()`. If this `filter()` step is not added, you will see that the largest pyramid bars are of "missing". These are the patients old enough to not have a months age group.  

If you add this filtering step, **you must also modify `age_pyramid()` by removing its first argument `linelist_cleaned,`**. The dataset is already given to the command in the `filter()` and is passed to `age_pyramid()` via piping.  

Note that the filter step does not drop any observations from the `linelist_cleaned` object itself. Because the filter is **not** being assigned (`<-`) to over-write `linelist_cleaned`, this filter is only temporarily applied for the purpose of producing the age pyramid.

```{r age_pyramid_under_two, fig.width = 6}
# plot age pyramid by month groups, for observations under 2 years
linelist_cleaned %>% 
  filter(age_years <= 2) %>%
  age_pyramid(age_group = "age_group_mon",
                   split_by = "sex") +
                  # stack_by = "case_def") +
  labs(y = "Cases (n)", x = "Age group (months)") + # change axis  labels (nb. x/y flip)
  theme(legend.position = "bottom",     # move legend to bottom
        legend.title = element_blank(), # remove title
        text = element_text(size = 18)  # change text size
       )
```

### Inpatient statistics  

The text following the age pyramids uses **in-line code** to describe the distribution of outpatient and inpatient observations, and descriptive statistics of the length of stay for inpatients. **The Am Timan dataset does not have length of stay, so it is best to delete those related sentences related to `obs_days` for the final report**.  

### Symptom and lab descriptive tables  

This next code also uses the `tbl_summary()` function to **create descriptive tables of all the variables that were included in the `SYMPTOMS` variable list**.  

* In the `tbl_summary()` function, `SYMPTOMS` (the value supplied to the second argument) is an object we defined in data cleaning that is a list of variables to tabulate. If this code produces an error about an "Unknown column", ensure that the variables in the object `SYMPTOMS` are *all* present in your dataset (and spelled correctly). 
* Also in `tbl_summary()`, the argument `keep =` must represent the character value to be counted for the the table. As these Am Timan variables are still in French, **we change `keep = "Yes"` to `keep = "Oui"`**.  
* The `mutate()` function is aesthetically changing variable names with underscores to spaces.   

```{r describe_by_symptoms, warning=FALSE, message=FALSE}
# get counts and proportions for all variables named in SYMPTOMS

linelist_cleaned %>% 
  select(all_of(SYMPTOMS)) %>% 
  tbl_summary() %>% 
  # fix the way symptoms are displayed
  # str_replace_all switches underscore for space in the variable column
  # str_to_sentence makes the first letter capital, and all others lowercase
  modify_table_body(
    ~.x %>% 
      mutate(label = str_to_sentence(str_replace_all(label, "_", " ")))
  ) %>% 
  # change to flextable format
  as_flex_table() %>%
  # make header text bold (using {flextable})
  bold(part = "header") %>% 
  # make your table fit to the maximum width of the word document
  set_table_properties(layout = "autofit")

```
These tables may be large and unwieldy at first, until you implement further data cleaning steps. Using `fct_recode()` to convert unexpected or mistakenly typed values as desired. For example, "OuI" should be recoded as "Oui" and "U" as "Unknown".  

The code for making a table with the `LAB` columns is very similar:


```{r describe_by_lab, warning=FALSE, message=FALSE}
# get counts and proportions for all variables named in SYMPTOMS

linelist_cleaned %>% 
  select(all_of(LABS)) %>% 
  tbl_summary() %>% 
  # fix the way symptoms are displayed
  # str_replace_all switches underscore for space in the variable column
  # str_to_sentence makes the first letter capital, and all others lowercase
  modify_table_body(
    ~.x %>% 
      mutate(label = str_to_sentence(str_replace_all(label, "_", " ")))
  ) %>% 
  # change to flextable format
  as_flex_table() %>%
  # make header text bold (using {flextable})
  bold(part = "header") %>% 
  # make your table fit to the maximum width of the word document
  set_table_properties(layout = "autofit")

```

The step in data cleaning where we converted 0, 1, "yes", "pos", "neg", etc. to standardized "Positive" and "Negative" was crucial towards making this table readable.



### Case Fatality Ratio (CFR)

The opening text of this chunk with in-line code must be edited to match our Am Timan data. The second in-line code references the variable `exit_status` - **this must now reference the variable `exit_status2`**. Also, we do not have "Dead on Arrival" recorded in our dataset, so **that part of the sentence should be deleted**.  

Likewise, **the code section on time-to-death does not apply to our dataset and should be deleted.**  

Overall CFR is produced with the code below. Note the following:  

* This code requires the variable `patient_facility_type`, that it has a value "Inpatient", and the `DIED` variable.  
* A filter is applied that restricts this table to Inpatient observations only.  
* The function `add_cfr()` from the package {epitabulate} modifies the table to include the CFR.  

```{r overall_cfr, warning = FALSE, message = FALSE}
# use arguments from above to produce overal CFR
linelist_cleaned %>%
  filter(patient_facility_type == "Inpatient") %>%
  select(DIED) %>%
  gtsummary::tbl_summary(
    statistic = everything() ~ "{N}",
    label = DIED ~ "Inpatients"
    ) %>%
  modify_header(stat_0 = "Cases (N)") %>% 
  # Use wrapper function to calculate cfr
  epitabulate::add_cfr(deaths_var = "DIED") %>% 
  ## make variable names bold 
  bold_labels() %>% 
  # change to flextable format
  as_flex_table() %>%
  # make header text bold (using {flextable})
  bold(part = "header") %>% 
  # make your table fit to the maximum width of the word document
  set_table_properties(layout = "autofit")

```

The next code adds the `include = "sex"` argument to `tbl_summary()` (which stratifies the CFR by sex).  

```{r cfr_by_sex, warning = FALSE, message = FALSE}
linelist_cleaned %>%
  filter(patient_facility_type == "Inpatient") %>%
  select(DIED, sex) %>%
  mutate(sex = forcats::fct_explicit_na(sex, "Missing")) %>%
  gtsummary::tbl_summary(
    include = sex,
    statistic = sex ~ "{n}",
    label = sex ~ "Gender") %>%
  modify_header(stat_0 = "Cases (N)") %>% 
  # Use wrapper function to calculate cfr
  epitabulate::add_cfr(deaths_var = "DIED") %>%
  # make variable names bold 
  bold_labels() %>% 
  # change to flextable format
  as_flex_table() %>%
  # make header text bold (using {flextable})
  bold(part = "header") %>% 
  # make your table fit to the maximum width of the word document
  set_table_properties(layout = "autofit")
```

This table shows the CFR by age group. Note how the `include = ` argument is now assigned to `age_group`.


```{r cfr_by_age_group, warning = FALSE, message = FALSE}

linelist_cleaned %>% 
  filter(patient_facility_type == "Inpatient") %>%
  select(DIED, age_group) %>%
  gtsummary::tbl_summary(
    include = age_group,
    statistic = age_group ~ "{n}",
    missing = "no",
    label = age_group ~ "Age group"
    ) %>%
  modify_header(stat_0 = "Cases (N)") %>%
  # Use wrapper function to calculate cfr
  add_cfr(deaths_var = "DIED") %>%
  # make variable names bold 
  bold_labels() %>% 
  # change to flextable format
  as_flex_table() %>%
  # make header text bold (using {flextable})
  bold(part = "header") %>% 
  # make your table fit to the maximum width of the word document
  set_table_properties(layout = "autofit")
```


The commented code below examines **CFR by case definition**. Note that this is dependent upon our working `case_def` variable.

```{r cfr_by_case_def, warning = FALSE, message = FALSE}
# Use if you have enough confirmed cases for comparative analysis 

linelist_cleaned %>%
  filter(patient_facility_type == "Inpatient") %>%
  select(DIED, case_def) %>%
  gtsummary::tbl_summary(
    include = case_def,
    statistic = case_def ~ "{n}",
    missing = "no",
    label = case_def ~ "Case definition"
    ) %>%
  modify_header(stat_0 = "Cases (N)") %>%
  # Use wrapper function to calculate cfr
  add_cfr(deaths_var = "DIED") %>%
  # make variable names bold 
  bold_labels() %>% 
  # change to flextable format
  as_flex_table() %>%
  # make header text bold (using {flextable})
  bold(part = "header") %>% 
  # make your table fit to the maximum width of the word document
  set_table_properties(layout = "autofit")
```



### Attack Rate
To use the attack rate section, we need to create an object `population` from the sum of population counts in the population figures. **Because we only imported region-based population counts, we must change this command** to reflect that we do not have `population_data_age`, but rather `population_data_region`.

```{r eval=FALSE}
# OLD command from template 
# population <- sum(population_data_age$population)
```

Running the correct command and printing the value of population, we see that the sum population across regions is estimated to be 62336.  

```{r}
# CORRECTED command for Am Timan exercise
population <- sum(population_data_region$population)

population
```

The first line of code below creates a multi-part object `ar` with the number of cases, population, attack rate per 10,000, and lower and upper confidence intervals (you can run just this line to verify yourself). This function comes from our package {sitrep}. The subsequent commands alter the aesthetics of to produce a neat table with appropriate column names. 

```{r attack_rate, warning = FALSE, message = FALSE}
# calculate the attack rate information and store them in object "ar""
ar <- attack_rate(nrow(linelist_cleaned), population, multiplier = 10000)

# print the table
ar %>%
  merge_ci_df(e = 3, digits = 1) %>% # merge the lower and upper CI into one column
  rename("Cases (n)" = cases, 
         "Population" = population, 
         "AR (per 10,000)" = ar, 
         "95%CI" = ci) %>% 
  # produce styled output table with auto-adjusted column widths with {flextable}
  qflextable() %>% 
  # make header text bold (using {flextable})
  bold(part = "header") %>% 
  # make your table fit to the maximum width of the word document
  set_table_properties(layout = "autofit") %>% 
  ## set to only show 1 decimal place 
  colformat_double(digits = 1)
```

We are unable to calculate the attack rate by age group, because we do not have population counts for each age group. **Comment out (#) this code.**  

Mortality attributable to AJS is also not appropriate for this example. **Comment out (#) the 4 related code chunks.**


## Time


The next section of the template produces analyses and outputs related to time, such as epidemic curves.

The first in-line code sentence clearly states how many cases are missing date of onset, as this is important to be able to correctly interpret the epidemic curve.  

### Epidemic curve  

The `basic_curve` epidemic curve is achieved using `ggplot()` and `geom_bar()` on the case linelist.  


```{r epicurve, fig.width = 6}

## make plot
basic_curve <- linelist_cleaned %>%
  filter(!is.na(epiweek_date)) %>%
  
  ## plot those not missing using the date variable 
  ggplot(aes(x = epiweek_date)) +
  ## plot as a bar
  geom_bar(width = 7,
           fill = "red",
           colour = "black") +
  ## align plot with bottom of plot
  scale_y_continuous(expand = c(0, NA)) +
  ## control x-axis breaks and display
  scale_x_date(
    breaks = all_weeks_date,
    labels = all_weeks) + 
  ## use pre-defined labels
  epicurve_labels +
  ## use pre-defined theme (design) settings
  epicurve_theme

# show your plot (stored for later use) 
basic_curve
```

This x-axis is difficult to read with every week displayed. To increase readability we can uncomment and use the helper code at the bottom of that code chunk. This modifies the already-defined `basic_curve` object (with a `+`) to specify that the x-axis shows only a specified number of breaks/dates (in this case 6). Note: A warning that a scale for x is being replaced is okay.  


The data are still shown by week - it is only the labels on the scale that have changed.



```{r, fig.width = 6}
# Modifies basic_curve to show only 6 breaks in the x-axis 
basic_curve + scale_x_date(breaks = "3 months", date_labels = "%b %Y")
```

An in-line code sentence then calculates and prints the week of the peak of the outbreak.  

### Moving average  

The next code chunk uses the {slider} function to calculate a moving average (mean) of number of cases in the two weeks prior, plotted as a line on top of the epidemic curve. You can adjust the window in the top line of this chunk (default set to 14 days).  

```{r epicurve_moving, fig.width = 6}
## define the width of the moving average window here
moving_window_days <- 14

## calculates the moving average using counts by epiweek
moving_avg <- linelist_cleaned %>% 
  ## count cases by epiweek
  count(epiweek_date) %>% 
  ## create column with sliding window for day and window prior
  mutate(indexed_window = slider::slide_index_dbl(
    n,                       # calculate on new_cases
    .i = epiweek_date,       # indexed with epiweek_date 
    .f = ~mean(.x, na.rm = TRUE),     # function is mean() with missing values removed
    .before = days(moving_window_days - 1)))

## create plot
moving_plot <-  basic_curve + 
  geom_line(data = moving_avg,
       mapping = aes(x = epiweek_date, y = indexed_window)) +
  
  # adjust title and subtitle for this particular plot
  labs(title = str_glue("Cases by date of onset"),
       subtitle = str_glue("{moving_window_days}-day moving average
                           Source: MSF data from {reporting_week}"))

# show your plot (stored for later use) 
moving_plot + scale_x_date(breaks = "3 months", date_labels = "%b %Y")
```


### Biweekly epicurve  

A biweekly epidemic curve can be produces using the next code chunk. In the first step, the linelist is aggregated into weekly counts and then a biweekly sum is created using {slider}. Finally, In the `geom_bar()` function, the `width = ` argument is set to equal `14` (days).  

```{r biweekly_epicurve, fig.width = 6}

## calculates the moving sum using counts by epiweek
moving_sum <- linelist_cleaned %>% 
  ## count cases by date of onset
  count(epiweek_date) %>% 
  mutate(biweekly = slide_sum(x = n, 
            after = 1, 
            step = 2)) %>% 
  filter(!is.na(biweekly))

ggplot(moving_sum, aes(x = epiweek_date, y = biweekly)) +
  ## plot as a bar
  geom_bar(
    ## use the count we calculated above
    stat = "identity", 
    ## make each bar two weeks wide
    width = 14,
    fill = "red",
    colour = "black") +
  ## align plot with bottom of plot
  scale_y_continuous(expand = c(0, NA)) +
  ## control x-axis breaks and display
  scale_x_date(
    breaks = all_weeks_date,
    labels = all_weeks) + 
  ## use pre-defined labels
  epicurve_labels +
  ## use pre-defined theme (design) settings
  epicurve_theme +
  scale_x_date(breaks = "3 months", date_labels = "%b %Y")

```


Using the next code chunk, you can group the data in the epidemic curve by the values in `sex`, thus producing a "stacked" epidemic curve. You can also change out `sex` for another column.  


```{r incidence_by_gender, fig.width = 6, fig.height=8}

# get counts by gender
linelist_cleaned %>%
  filter(!is.na(epiweek_date)) %>%
  ## plot those not missing using the date variable 
  ggplot(aes(x = epiweek_date, fill = fct_rev(fct_explicit_na(sex)))) +
  ## plot as a bar
  geom_bar(width = 7,
           colour = "black") +
  # align plot with bottom of plot
  scale_y_continuous(expand = c(0, NA)) +
  # control x-axis breaks and display
  scale_x_date(
    breaks = all_weeks_date,
    labels = all_weeks) + 
  # use pre-defined labels
  epicurve_labels +
  # use pre-defined theme (design) settings
  epicurve_theme + 
  scale_x_date(breaks = "3 months", date_labels = "%b %Y")
```

Such as by, `case_def`:  

```{r incidence_by_case_def, fig.width=6, fig.height=8}

linelist_cleaned %>%
  filter(!is.na(epiweek_date)) %>%
  
  ## plot those not missing using the date variable 
  ggplot(aes(x = epiweek_date, fill = fct_rev(fct_explicit_na(case_def)))) +
  ## plot as a bar
  geom_bar(width = 7,
           colour = "black") +
  # align plot with bottom of plot
  scale_y_continuous(expand = c(0, NA)) +
  # control x-axis breaks and display
  scale_x_date(
    breaks = all_weeks_date,
    labels = all_weeks) +
  # use pre-defined labels
  epicurve_labels +
  # use pre-defined theme (design) settings
  epicurve_theme + 
  scale_x_date(breaks = "3 months", date_labels = "%b %Y")
```

Remember also that you can limit the entire plot to a subset of observations by first filtering the dataset and passing the result to the `ggplot()` function, as is done in the above examples filtering to remove rows with a missing `epiweek_date` value.   


### Analysis of other indicators by time  

The chunk `attack_rate_per_week` does not produce an epidemic curve, but does print a table of the the attack rate by week. This table could be very long if your epidemic is protracted, so you may choose to use a different unit of time than week.   


```{r attack_rate_per_week}

# counts and cumulative counts by week
cases <- linelist_cleaned %>%
  arrange(date_of_onset) %>%        # arrange by date of onset
  count(epiweek, .drop = FALSE) %>% # count all epiweeks and include zero counts
  mutate(cumulative = cumsum(n))    # add a cumulative sum


# attack rate for each week
ar <- attack_rate(cases$n, population, multiplier = 10000) %>% 
  bind_cols(select(cases, epiweek), .) # add the epiweek column to table


ar %>%
  merge_ci_df(e = 4) %>% # merge the lower and upper CI into one column
  rename("Epiweek" = epiweek, 
         "Cases (n)" = cases, 
         "Population" = population, 
         "AR (per 10,000)" = ar, 
         "95%CI" = ci) %>% 
  # produce styled output table with auto-adjusted column widths with {flextable}
  qflextable() %>% 
  # make header text bold (using {flextable})
  bold(part = "header") %>% 
  # make your table fit to the maximum width of the word document
  set_table_properties(layout = "autofit") %>% 
  ## set to only show 1 decimal place 
  colformat_double(j = -1, digits = 1)
```

The next chunks offers a similar calculation, but for cumulative attack rate, and CFR  


### Plotting AR, CFR, and the Epidemic Curve togethr


These code chunks produce plots of the AR (in the population), CFR (among inpatients only), and the AF and CFR as line graphs together with the Epidemic Curve.  

**The chunk, `ar_line_graph`, prepares a plot of attack rate by epiweek, but does not yet print it.** You can print it alone by running a command of its name `ar_plot`.

This `ar_plot` depends on you having already run the `attack_rate_per_week` chunk. If you have not, you may see "Error: epiweek not found"

Each of the plots involves several steps:  

1. Creating the data object e.g. `ar` or `cfr`  
2. Creating the plot e.g. `ar_plot` or `cfr_plot`  

The plots are finally combined with the {patchwork} package into one image.  

The code below uses the `ar` object created above to produce *and store* the AR plot.  


```{r ar_line_graph, fig.width=6, fig.height=8}

ar_plot <- ggplot(ar, aes(x = as.Date(epiweek), group = 1)) +
  # add confidence intervals as a ribbon
  geom_ribbon(aes(ymin = lower, ymax = upper), 
              color = "blue",  fill = "blue", 
              linetype = 2, alpha = 0.2, show.legend = FALSE) +
  # add AR as a line
  geom_line(aes(y = ar), color = "blue", show.legend = FALSE) +
  # set origin for axes
  scale_y_continuous(expand = c(0, NA)) +
  # control x-axis breaks and display
  scale_x_date(
    breaks = all_weeks_date,
    labels = all_weeks) +
  # add labels to axes and below chart
  labs(x = "Calendar week", y = "AR [95% CI]", subtitle = "Attack Rate (per 10,000)") + 
  # change visuals of dates and remove legend title
  epicurve_theme
```


Then code is run to produce *and store* the CFR plot.  

First we make `cfr` from the `linelist_cleaned`:  

```{r cfr_per_week, warning = FALSE, message = FALSE}

# calculate CFR among inpatients by week
cfr <- linelist_cleaned %>%
  filter(patient_facility_type == "Inpatient") %>%
  mutate(epiweek = as.character(epiweek)) %>% 
  case_fatality_rate_df(str_detect(exit_status, "Dead"), group = epiweek)
```

Then we make the plot and save it as `cfr_plot`. Note that it is not yet printed.  

```{r cfr_line_graph, fig.width=6, fig.height=8}

cfr_plot <- ggplot(cfr, aes(x = as.Date(yearweek(as.character(epiweek)))), group = 1) +
  # add confidence intervals as a ribbon
  geom_ribbon(aes(ymin = lower, ymax = upper), 
              color = "red", fill = "red", 
              linetype = 2, alpha = 0.2, show.legend = FALSE) +
  # add CFR as a line
  geom_line(aes(y = cfr), color = "red", show.legend = FALSE) +
  # set origin for axes
  scale_y_continuous(expand = c(0, NA)) +  
  # control x-axis breaks and display
  scale_x_date(
    breaks = all_weeks_date,
    labels = all_weeks) +
  # add labels to axes and below chart
  labs(x = "Calendar week", y = "CFR [95% CI]", 
       subtitle = "Case Fatality Ratio [95% CI] Among Inpatients") + 
  # change visuals of dates and remove legend title
  epicurve_theme + 
  scale_x_date(breaks = "3 months", date_labels = "%b %Y")
```

To create a combination visualisation, we must adjust some final setting for the x-axes.  

```{r epicurve_ar_cfr, warning = FALSE, message = FALSE, fig.height = 10}

no_x <- theme(axis.text.x = element_blank(),
              axis.title.x = element_blank())
```

And finally we plot the Epidemic Curve, the AR plot, and the CFR plot all over each other.  

```{r, fig.height=10, fig.width=8}
## use patchwork to stack the three plots on top of one another
(basic_curve + no_x) / (ar_plot + no_x) / cfr_plot
```

### Admissions and exits

The final two tables print the admissions and exits by epidemic week, using code previously described.  












## Place



The final section of the template produces tabular analysis outputs related to *place*. If you have shapefiles for your situation, you can also produce maps.

### Tabular place analyses (non-maps)

To get a basic **descriptive table of cases by region** (which in our case was originally the variable *quartier* from the raw linelist), we use `tbl_cross()` function from {gtsummary} on the columns `patient_origin` and `patient_facility_type` (inpatient, outpatient, or missing). 

```{r describe_by_region_facility}

linelist_cleaned %>% 
  tbl_cross(
    row = patient_origin, 
    col = patient_facility_type, 
    percent = "column", 
    label = list(patient_origin = "Region", 
                 patient_facility_type = "Facility type")
  ) %>% 
  ## make variable names bold 
  bold_labels() %>% 
  # change to flextable format
  as_flex_table() %>%
  # make header text bold (using {flextable})
  bold(part = "header") %>% 
  # make your table fit to the maximum width of the word document
  set_table_properties(layout = "autofit")

```


A table of attack rates by region may be useful, and is produced in the chunk `attack_rate_by_region`.  

**First**, like when we calculated attack rate by epiweek in the `time_analysis` section, the cases are counted by `patient_origin`. This list of counts is then **joined** to the imported population data for each region **by `patient_origin` (the name of the region/quartier)**. These data are stored in the object `cases`.  

*Recall that this join requires exact matches* between `patient_origin` values in `linelist_cleaned$patient_origin` and `population_data_region$patient_origin` to match a region's cases to its population. In the data cleaning steps we converted the values in `linelist_cleaned` to all capital letters to match those in `population_data_region`.  

```{r attack_rate_by_region, warning = FALSE}

# count cases by region
cases <- linelist_cleaned %>% 
  count(patient_origin) %>%
  # add in population data
  left_join(population_data_region, by = "patient_origin") 


# calculate attack rate for region
ar_region <- attack_rate(cases$n, cases$population, multiplier = 10000) %>% 
  # add the region column to table
  bind_cols(select(cases, patient_origin), .) %>% 
  rename("Region" = patient_origin, 
         "Cases (n)" = cases, 
         "Population" = population, 
         "AR (per 10,000)" = ar, 
         "Lower 95%CI" = lower,
         "Upper 95%CI" = upper) 

ar_region %>% 
  merge_ci_df(e = 4) %>% # merge lower and upper CI in to one column 
  rename("95%CI" = ci) %>%  # rename single 95%CI column
  # produce styled output table with auto-adjusted column widths with {flextable}
  qflextable() %>% 
  # make header text bold (using {flextable})
  bold(part = "header") %>% 
  # make your table fit to the maximum width of the word document
  set_table_properties(layout = "autofit") %>% 
  ## set to only show 1 decimal place 
  colformat_double(digits = 1)
```


You could then also plot this on a bar chart with confidence intervals.  

```{r, fig.width=8, fig.height=8}
ar_region %>% 
  drop_na(`AR (per 10,000)`) %>% 
ggplot(aes(x = reorder(Region, `AR (per 10,000)`),
               y = `AR (per 10,000)`)) +
  # plot as bars (identity = as is)
  geom_bar(stat = "identity", col = "black", fill = "red") + 
  # add CIs
  geom_errorbar(aes(ymin = `Lower 95%CI`, ymax = `Upper 95%CI`), width = 0.2) +
  # set origin for axes
  scale_y_continuous(expand = c(0,0)) +  
  # add labels to axes and below chart
  labs(x = "Region", y = "AR (per 10,000)", 
       captions = str_glue("Source: MSF data from {reporting_week}")) + 
  epicurve_theme
```


Note that there were many quartiers without population denominators that make the graph difficult to read. Thus, **we modified the code by filtering `ar_region` with `drop_na()` to remove these. Now any quartiers with missing denominators will not appear in the graph.  




```{r bar_attack_rate_by_region, fig.width=8, fig.height=8}
# filter ar for non-zero quartiers
ar_region <- filter(ar_region, !is.na(`AR (per 10,000)`))

# plot with the region on the x axis sorted by increasing ar
# ar value on the y axis 
ggplot(ar_region, aes(x = reorder(Region, `AR (per 10,000)`),
               y = `AR (per 10,000)`)) + 
  geom_bar(stat = "identity", col = "black", fill = "red") + # plot as bars (identity = as is)
  geom_errorbar(aes(ymin = `Lower 95%CI`, ymax = `Upper 95%CI`), width = 0.2) + # add CIs
  scale_y_continuous(expand = c(0,0)) +  # set origin for axes
  # add labels to axes and below chart
  labs(x = "Region", y = "AR (per 10,000)", 
       captions = str_glue("Source: MSF data from {reporting_week}")) + 
  epicurve_theme
```

A similar analysis can be done to show mortality rates (be sure to use `exit_status2`!)



### Maps 

Mapping in R can be intimidating at first because it is different than ArcGIS, QGIS, or other mapping software that primarily use a point-and-click interface. However, you will see that using commands to create maps can be much faster, more replicable, and transparent for whomever you share code with.  


*Shapefiles are a format for storing geometric location and attribute data for geographic features - for example **either points (e.g. GPS points of cases or households surveyed), lines (e.g. road network), or polygons (e.g. district boundaries)* 

We do have shapefiles for the Am Timan exercise. The shapefile for the quartiers is accessible from within the {sitrep} R package by running the following command:  

```{r}
# saves shapefile stored in R package as "map"
map <- shape_quartier_internal
```
 

A shapefile is actually made of several files. Viewing the shapefile in your computer's internal folders you will often see 4 or 5 files with the same name but different extensions (e.g. .shp, .prj, etc.). As long as these files are located in the same folder the shapefile will load correctly when you import with the .shp extension.

Because we have shapefiles for this exercise, we **must delete the command that generates the fake shapefile** 

`map <- gen_polygon(regions = unique(linelist_cleaned$patient_origin))`

Note the checklist for using shapefiles:  

* Your shapefile can be a polygon or points (polygons do not have to be contiguous - they do not have to touch each other)  
* Ensure that the region names in the shapefile match *exactly* the region names in `linelist_cleaned` (in most shapefiles there is an ID or NAME variable in the spreadsheet file)  
* Make sure the coordinate reference system (CRS) of the shapefile is WGS84 (this is the standard coordinate system used by most GPS systems, but it is good to ask whomever collected to data to be sure)  
 
When giving the name of a shapefile to `here()`, list the file name with a .shp extension.
   
A generic example of importing a shapefile might look like this:  
  
  
```{r read_shapefiles_display, eval=FALSE}
## reading in a shapefile 
map <- read_sf(here("my_map_folder", "Quartiersshape.shp"))
```


The next code will check the coordinate reference system (CRS) of the shapefile, as it was read-in. In our exercise, we see an EPSG code of 4326, which is the correct code for WGS84 (you can also see below the EPSG code that the datum is WGS84).

```{r}
## check the coordinate reference system (CRS)
st_crs(map)
```

If there is no CRS specified for your shapefile, you can set the CRS to WGS84 with the following command, which uses the EPSG code of 4326.   

```{r, eval=FALSE }
## set the CRS if not present using EPSG value
## this is the most commonly used 
map <- st_set_crs(map, value = 4326) # Sets to WGS84
```

### Attack rate map  

These commands reference the `ar_region` table prepared in previous sections of the report, which contains region names, number of cases, populaton AR, and 95%CI.  

#### Create categories of ARs by region  

**Create a *categorical* variable of the attack rate.** These categories will appear in the map as the legend - and the regions will be colored by these groups.  

1) The first command defines `max_ar` as the highest attack rate value in the object `ar_region`. This is referenced later.  

2) In the next command, `breakers` is defined as a sequence of integers. Zero is included by default because it should be a separate color on the final map. After zero, quartiles are determined up to the `max_ar` using the {sitrep} function `find_breaks()`.  

3) In the last command, `mutate()` is used to make `ar_map` (a copy `ar_region`), and add the `categories` variable. This variable is created using `breakers` to define the breaks in the numeric variable `AR (per 10,000)`. 

```{r}
# save the maximum AR in the ar_region table
max_ar    <- max(ar_region$`Upper 95%CI`, na.rm = TRUE)

# create groups - one group will be 0 only, then up to 4 more groups evenly
# divided to maximum
breakers <- as.integer(c(
                         # include zero as a standalone group
                         0, 
                         # 1 to 4 divisions, snapping the boundaries to the nearest 100
                         find_breaks(n = max_ar, breaks = 4, snap = 100)
             ))

## create a categorical variable using the age_categories function 
## (we aren't using ages - but it functions the same way!)
ar_map <- mutate(ar_region, 
             categories = age_categories(`AR (per 10,000)`, 
                                         breakers = breakers))



```


#### Join the table to the shapefile  

**Join the object `map` (shapefile) to `ar_map` (a data frame)**. The spatial object `map` is joined to the data frame `ar_map` with a **left join**. Both of these objects contain tables (you can click them in the Environment pane to see this). In the command, they are linked by the `Name` variable in `ar` and the `Region` variable in `map`.  

**You may need to modify your code so that the join occurs on `Name` in the object `map`, and not `name`.**

```{r}
## Join table to shapefile -----------------------------------------------------
mapsub <- left_join(map, ar_map, by = c("Name" = "Region"))

```

Make sure you modify your `left_join` code to use "Name" with a capital N, as that is how the variable is spelled in the object 'map'. In your other data these variable names may be completely different and must be changed

The resulting object `mapsub` is a product of this **left join**. A left join will include all values from the Left-hand Side (LHS) and only include values from the Right-hand Side (RHS) if they match to LHS values. Now, the attack rates are joined to the region names and to their GIS geometries. You can verify this by Viewing the object `mapsub` (click on it in the Environment pane).

#### Plot cases or AR by region


**Now, `ggplot()` is used to create the map.** Remember that `ggplot()` subcommands must be connected by plus symbols at the end of each line (`+`) (except the last line). Read more about ggplot2 and `ggplot()` in the Advanced R Basics page.  

**`geom_sf()`** identifies `mapsub` as providing the data for the map, and filling the polygons by values in the variable `categories`.  

**`coord_sf()`** set to `NA` prevents gridlines from being drawn on the map.  

**`annotation_scale()`** adds a scalebar at the bottom of the map.  

**`scale_fill_brewer()`** uses a color palatte from orange to red ("OrRd"), keeps all color levels even if there are no polygons at that level (`drop=FALSE`), and specifies a title for the legend. If you want a different color scale, use this code from the package **RColorBrewer** to see the possible options (hint: after running this code click "Zoom" in your Plot/Viewer pane to see more clearly). Then replace "OrRd" with your desired code. You can read more about colors [here](http://www.sthda.com/english/wiki/print.php?id=85).  

```{r, fig.width=6}
library(RColorBrewer)
display.brewer.all()
```

**`geom_sf_text()`** specifies labels for the polygons, including which variable to use for labels (`Name`), and the text colour. **Again, we need to modify this `ggplot()` reference to the variable "name" to now reference "Name" with a capital N. In your other data the variable may be completely different.**

**`theme_void()`** specifies a theme for the map that does not include coordinates or axes.  

```{r fig.width=8}
## you could also fill by cases using `Cases (n)` in the fill option instead of `AR (per 10,000)`
ggplot() +
  # shapefile as polygon
  geom_sf(data = mapsub, aes(fill = categories), col = "grey50") + 
  # needed to avoid gridlines being drawn
  coord_sf(datum = NA) + 
  # add a scalebar
  annotation_scale() + 
  # color the scale to be perceptually uniform 
  # drop FALSE keeps all levels 
  # name allows you to change the legend title 
  scale_fill_brewer(drop = FALSE, palette = "OrRd", name = "AR (per 10,000)") + 
  # label polygons
  geom_sf_text(data = mapsub, aes(label = Name), colour = "grey50", size = 2) + 
  # remove coordinates and axes
  theme_void() 
```

**You can use this alternative to make the labels not overlap.** 

**`ggrepel()`** arranges point/polygon labels so they do not overlap and are easier to read. You can see more examples of `ggrepel()` [in this vignette](https://mran.microsoft.com/snapshot/2017-08-20/web/packages/ggrepel/vignettes/ggrepel.html)  



The `geom_sf_text()` line has been replaced by ggrepel.



```{r fig.width=8, fig.height=8}
## plotting with polygons ------------------------------------------------------

ggplot() +
  geom_sf(data = mapsub, aes(fill = categories), col = "grey50") + # shapefile as polygon
  coord_sf(datum = NA) + # needed to avoid gridlines being drawn
  annotation_scale() + # add a scalebar
  scale_fill_brewer(drop = FALSE, palette = "OrRd", name = "AR (per 10,000)") + 
  
  # Optional ggrepel section to make labels appear more nicely
  ggrepel::geom_label_repel(data = mapsub, aes(label = Name, geometry = geometry), 
                            stat = "sf_coordinates", 
                            min.segment.length = 0) +

  theme_void()  # remove coordinates and axes
```


### Map iteration  


Lastly, **there is an efficient way to create attack rate maps and bar plots of incidence for each epiweek**. In this case, we want to make the map and bar chart for each region. You can read more about iteration in the Epidemiologist R Handbook chapter [Iteration, Loops, and Lists](https://epirhandbook.com/en/iteration-loops-and-lists.html).  

A few things to note about this chunk:  

* The first half is very similar to the set-up of the previous map.  
* You must also modify the `left_join` command to list the correct variable names (e.g. **modify to join on `Name` instead of `name`**)
* These maps are not labeled with region names.  
* The creation of the barplots depends on the object `epicurve_theme`, which was defined in the Time section of the template. If that section is not created, then this object will be missing and cause errors.  

First we run preparation code that should look very familiar to the steps above.
```{r map_for_loop_epiweek, fig.width=8, fig.height=8}
## Prepare data ---------------------------------------------------------------

# change region variable to a factor so that zero counts can be included
linelist_cleaned$patient_origin <- as.factor(linelist_cleaned$patient_origin)


# case counts
cases <- linelist_cleaned %>% 
  group_by(epiweek) %>%
  count(patient_origin, .drop = FALSE) %>%   # cases for each week by region
  left_join(population_data_region, by = "patient_origin")    # merge population data 


# attack rate for region
ar <- attack_rate(cases$n, cases$population, multiplier = 10000) %>% 
  # add the region column to table
  bind_cols(select(cases, epiweek, patient_origin), .) %>% 
  rename("Region" = patient_origin, 
         "Cases (n)" = cases, 
         "Population" = population, 
         "AR (per 10,000)" = ar, 
         "Lower 95%CI" = lower,
         "Upper 95%CI" = upper)

## define the last 10 weeks available in your dataset
which_weeks <- tail(unique(cases$epiweek), n = 10)

## only keep the last four weeks of data 
ar <- ar %>% 
  filter(epiweek %in% which_weeks)

# define the maximum number of cases / max AR for the color palette 
max_cases <- max(cases$n, na.rm = TRUE) 
max_ar    <- max(ar$`Upper 95%CI`, na.rm = TRUE)


# define breaks for standardising color palette
breakers <- as.integer(c(0, # include zero as a standalone group
             find_breaks(max_ar, breaks = 4, snap = 25) # four breaks rounded to nearest 25
             ))


# add a categorical variable for AR
# (replace with `Cases (n)` if you want to map cases)
ar <- mutate(ar, 
             categories = age_categories(`AR (per 10,000)`, 
                                         breakers = breakers))
```

Now, we can use the {purrr} package to *iterate* (repeat a process many times). For each week in the final 10 weeks of the outbreak, we left join the attack rate and the shapefile. Then we create a map, create a barplot, and combine the two plots together.  


```{r, fig.width=8, fig.height=8}

# go through each epiweek, filter and plot the data
purrr::walk(.x = which_weeks, 
           .f = ~{
  this_ar <- filter(ar, epiweek == yearweek(.x))
  
  ## map 
  mapsub <- left_join(map, this_ar, by = c("Name" = "Region"))
  
  ## choropleth 
  map_plot <- ggplot() +
    ## shapefile as polygon
    geom_sf(data = mapsub, aes(fill = categories), col = "grey50") + 
    ## needed to avoid gridlines being drawn
    coord_sf(datum = NA) + 
    ## add a scalebar
    annotation_scale() + 
    ## color the scale to be perceptually uniform (keep levels)
    scale_fill_brewer(drop = FALSE, 
                    palette = "OrRd", 
                    name = "AR (per 10,000)") +  
    ## remove coordinates and axes
    theme_void() 
  
  
  ## plot with the region on the x axis sorted by increasing ar
  ## ar value on the y axis 
  barplot <- this_ar %>%
    # Filter values in this_ar to only show quartiers with attack rates > 0 in barplot
    filter(`AR (per 10,000)` > 0) %>% 
    # make plot
    ggplot(aes(x = reorder(Region, `AR (per 10,000)`),
                                 y = `AR (per 10,000)`)) + 
    ## plot as bars (identity = as is)
    geom_col(col = "black", fill = "red") + 
     ## add CIs
    geom_errorbar(aes(ymin = `Lower 95%CI`, ymax = `Upper 95%CI`), width = 0.2) +
    ## set origin for axes
    scale_y_continuous(expand = c(0, 0), limits = c(0, max_ar)) +  
    ## add labels to axes and below chart
    labs(x = "Region", y = "AR (per 10,000)") +  
    epicurve_theme
  
  
  
  ## combine the barplot and map plot into one with patchwork and add a title
  ## in the top-left corner and a caption in the bottom-right
  print(
    barplot + map_plot + plot_annotation(
    title   = str_glue("Epiweek: {.x}"),
    caption = str_glue("Source: MSF data from {reporting_week}")
    )
  )
  
})
```







