---
title: "AJS outbreak report"
output: 
  word_document:
    keep_md: true
---

# Introduction

This is the final template used in the **sitrep** case study.  

Here we would give an executive summary of the report findings and recommentations. 
     


<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
This comment will not show up when you knit the document.

A comment with a title with slashes indicates a name of a code chunk.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->


## Installing and loading required packages 
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// setup \\\
--------------------------------------------------------------------------------

Several packages are required for different aspects of  analysis with *R*. 
You will need to install these before starting. 

These packages can be quite large and may take a while to download in the
field. If you have access to a USB key with these packages, it makes sense to
copy and paste the packages into your computer's R package library 
(run the command .libPaths() to see the folder path). 

For help installing packages, please visit https://r4epis.netlify.com/welcome
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->


```{r setup, include=FALSE, results='hide', message=FALSE, warning=FALSE}
## hide all code chunks in the output, but show errors
knitr::opts_chunk$set(echo = FALSE,       # hide all code chunks in output
                      error = TRUE,       # show errors if they appear, but don't stop
                      fig.width = 6*1.25, # Figure width
                      fig.height = 6      # Figure height
                     )



## set default NA to - in output, define figure width/height
options(knitr.kable.NA = "-")



## Installing required packages for this template
required_packages <- c("knitr",       # create output docs
                       "here",        # find your files
                       "dplyr",       # clean/shape data
                       "forcats",     # clean/shape data
                       "stringr",     # clean text
                       "rio",         # read in data
                       "ggplot2",     # create plots and charts
                       "sitrep",      # MSF field epi functions
                       "linelist",    # Functions for cleaning/standardising data
                       "incidence",   # create epicurves
                       "aweek",       # define epi weeks
                       "epitools",    # 2x2 tables and other epi goodness
                       "epitrix",     # epi helpers and tricks
                       "sf",          # encode spatial vector data
                       "ggspatial")   # plot maps

for (pkg in required_packages) {
  # install packages if not already present
  if (!pkg %in% rownames(installed.packages())) {
    install.packages(pkg)
  }
  
  # load packages to this current session 
  library(pkg, character.only = TRUE)
}



## Set default options for plots and charts

## set default text size to 16 for plots
## give classic black/white axes for plots
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(),
  panel.grid.major.x = element_line(color = "grey60", linetype = 3),
  panel.grid.major.y = element_line(color = "grey60", linetype = 3))
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// define_current_week \\\
--------------------------------------------------------------------------------

You need to set the week you want to report on. Generally, this is the previous
week. Put it below.

aweek::set_week_start will define the beginning of the week. The standard is
Monday.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r define_current_week}
## set the day that defines the beginning of your epiweek.
aweek::set_week_start("Monday")

## set current week 
reporting_week <- aweek::as.aweek("2017-W16")

## sets the labels in ggplot for the epicurves
epicurve_labels <- labs(x = "Calendar week", 
                        y = "Cases (n)", 
                        title = "Cases by week of onset",
                        subtitle = sprintf("Source: MSF data from %s", 
                                           reporting_week)
                       ) 
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// read_fake_data \\\
--------------------------------------------------------------------------------


To play with this template, you can create fake data. Comment out this chunk
when you are using real data.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r read_fake_data, warning = FALSE, message = FALSE}

#linelist_raw <- gen_data("AJS")

## fake lab data
# lab_results <- linelist_cleaned %>% 
#                 select(case_number) %>% 
#                 mutate(test_result = sample(c("Positive", "Negative"),
#                                             nrow(linelist_cleaned), 
#                                             replace = TRUE)
#                        )
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
You should load your data as a linelist, where each row is one case.

There are two options:
- Your data is from DHIS2: use read_DHIS_excel_data
- Your data is NOT FROM DHIS2 and in Excel, CSV, or Stata format: 
use read_nonDHIS_data
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// read_DHIS_data \\\
--------------------------------------------------------------------------------

This section is for data from DHIS2.
It uses the standardized MSF data dictionary.
If you didn't use the standardized data dictionary, go to read_nonDHIS_data.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->

```{r read_DHIS_data, warning = FALSE, message = FALSE}
### Read in data ---------------------------------------------------------------

## Excel file ------------------------------------------------------------------
## For a specific sheet, use "which"
# linelist_raw <- rio::import(here::here("linelist.xlsx"), which = "Sheet1")



## Excel file with password ----------------------------------------------------
## Use this section if your Excel has a password.

# install.packages(c("excel.link", "askpass"))
# library(excel.link)

# linelist_raw <- xl.read.file(here::here("linelist.xlsx"),
#                              xl.sheet = "Sheet1",
#                              password = askpass::askpass(prompt = "please enter file password"))
``` 



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// prep_DHIS_data \\\
--------------------------------------------------------------------------------

This section is for data from DHIS2 that you loaded in read_DHIS_data.

If you didn't use DHIS2 data, go to read_nonDHIS_data.

This step shows you the data dictionary. The data dictionary has variable names
in the data_element_shortname column. Possible values for each variable are
specified in code and name columns. Code has the shortened value and Name has
the full-text value.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->

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

## MSF AJS Dictionary ----------------------------------------------------------
## get MSF standard dictionary for AJS 
#linelist_dict <- msf_dict("AJS", compact = FALSE) %>%
#  select(option_code, option_name, everything())

## look at the standard dictionary by uncommenting the line below
# View(linelist_dict) 



## Clean column names ----------------------------------------------------------
## make a copy of your orginal dataset and name it linelist_cleaned
#linelist_cleaned <- linelist_raw

## define clean variable names using clean_labels from the epitrix package
#cleaned_colnames <- epitrix::clean_labels(colnames(linelist_raw))

## overwrite variable names with defined clean names
#colnames(linelist_cleaned) <- cleaned_colnames



## Standardising values --------------------------------------------------------
## This step fixes the values so that you can read them.
## values like 1/0 will be recoded as "Yes" / "No" based on the dictionary

#linelist_cleaned <- clean_variable_spelling(linelist_cleaned,
#                                            wordlists = filter(linelist_dict,
#                                                               !is.na(option_code)),
#                                            spelling_vars = "data_element_shortname",
#                                            sort_by = "option_order_in_set")
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// read_nonDHIS_data \\\
--------------------------------------------------------------------------------

This section is for data not from DHIS2.
If you have already loaded data from DHIS2, go to read_population_data.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->

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

## Excel file ------------------------------------------------------------------
## to read in a specific sheet use "which"
# YOU MAY NEED TO CHANGE THE FOLDERS REFERENCED IN THE HERE() COMMAND TO MATCH YOUR R PROJECT FOLDER STRUCTURE
linelist_raw <- rio::import(here::here("static", "public_data", "AJS_AmTiman.xlsx"), which = "linelist", na = ".")



## Excel file -- Specific range ------------------------------------------------
## you can specify a range in an excel sheet.
# linelist_raw <- rio::import(here::here("linelist.xlsx"), range = "B2:J102")



## Excel file with password ----------------------------------------------------
## use this section if your Excel has a password.

# install.packages(c("excel.link", "askpass"))
# library(excel.link)

# linelist_raw <- xl.read.file(here::here("linelist.xlsx"),
#                              xl.sheet = "Sheet1",
#                              password = askpass::askpass(prompt = "please enter file password"))


## CSV file --------------------------------------------------------------------
# linelist_raw <- rio::import(here::here("linelist.csv"))



## Stata data file -------------------------------------------------------------
# linelist_raw <- rio::import(here::here("linelist.dat"))
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// prep_nonDHIS_data \\\
--------------------------------------------------------------------------------

This section is for data not from DHIS2. 
If you have already loaded data from DHIS2, go to read_population_data.

It is more difficult to prepare the nonDHIS data. You can do it! It will just
take a little more work.

Checklist to update this script to match your data:
[ ] Comment out all lines in read_DHIS_data and prep_DHIS_data
[ ] Recode your variable names to match the dictionary
[ ] Recode variable contents to match the dictionary


This step shows you the data dictionary. The data dictionary has variable names
in the data_element_shortname column. Possible values for each variable are
specified in code and name columns. Code has the shortened value and Name has
the full-text value.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->

```{r prep_nonDHIS_data, warning = FALSE, message = FALSE}
## MSF AJS Dictionary ----------------------------------------------------------
## get MSF standard dictionary for AJS 
 linelist_dict <- msf_dict("AJS", compact = FALSE) %>%
   select(option_code, option_name, everything())

## look at the standard dictionary by uncommenting the line below
# View(linelist_dict) 

## You will need to recode your variables to match the data dictionary. This is
## addressed below.



## Clean column names ----------------------------------------------------------
## This step fixes the column names so they are easy to use in R.

## make a copy of your orginal dataset and name it linelist_cleaned
linelist_cleaned <- linelist_raw
 
## define clean variable names using clean_labels from the epitrix package
cleaned_colnames <- epitrix::clean_labels(colnames(linelist_raw))
 
## overwrite variable names with defined clean names
colnames(linelist_cleaned) <- cleaned_colnames



## Match column names ---------------------------------------------------------
## This step helps you match your variables to the standard variables.
## This step will require some patience. Courage!

## Use the function msf_dict_rename_helper() to create a template based on the
## AJS dictionary. This will copy a rename command like the one above to your
## clipboard.

msf_dict_rename_helper("AJS")

## Paste the result below and your column names to the matching variable.
## Be careful! You still need to be aware of what each variable means and what
## values it takes.
## If there are any variables that are in the MSF dictionary that are not in
## your data set, then you should comment them out, but be aware that some
## analyses may not run because of this. 



## Here is an EXAMPLE for changing a few specific names. function. In this
## example, we have the columns "gender" and "age" that we want to rename as
## "sex" and "age_years". 
## The formula for this is rename(data, NEW_NAME = OLD_NAME).

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

```




<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// read_population_data \\\
--------------------------------------------------------------------------------

This template uses population data to calculate things like incidence.

There are three options:
- You can read in a spreadsheet with age group and region data.
- You can put in the specific populations into the gen_population function. 
- If you have the total or regional populations, you can estimate the age group
from proportions.

Comment out the options you are not using.


Age group proportions are from the OCBA population denominators tool v1. The
proportions below are for sub-Saharan Africa in 2019. They are only an estimate!
If you have more specific proportions, you can use them below.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->

```{r read_population_data, warning = FALSE, message = FALSE}
## Read data -------------------------------------------------------------------
## This step reads in your population data from Excel.
## You may need to rename your columns.


# YOU MAY NEED TO CHANGE THE FOLDERS REFERENCED IN THE HERE() COMMAND TO MATCH YOUR R PROJECT FOLDER STRUCTURE
population_data_region <- rio::import(here("static", "public_data", "AJS_AmTiman_population_revised_for_template.xlsx"), which = "quartier")

## repeat preparation steps as appropriate



## Enter counts directly -------------------------------------------------------

# population_data_age <- gen_population(
#   groups = c("0-2", "3-14", "15-29", "30-44", "45+"),
#   counts = c(3600, 18110, 13600, 8080, 6600),
#   strata = NULL) %>%
#   rename(age_group = groups,
#     population = n)

# population_data_region <- gen_population(
#   groups = c("Village A", "Village B", "Village C", "Village D"),
#   counts = c(10000, 20000, 15000, 5000),
#   strata = NULL) %>%               # do not stratify by gender
#   rename(patient_origin = groups,  # rename columns (syntax is NEW NAME = OLD NAME)
#          population = n)


## Create counts from proportions ----------------------------------------------
## This step helps you estimate sub-group size with proportions.
## You need to replace the total_pop and proportions. You can change the groups
## to fit your needs.

# 
# # estimate population size by age group in years
# population_data_age <- gen_population(total_pop = 5000,
#   groups = c("0-2", "3-14", "15-29", "30-44", "45+"),
#   proportions = c(0.068, 0.3622, 0.276, 0.1616, 0.1321),
#   strata = NULL) %>%
#   rename(age_group = groups,
#     population = n)
# 
# 
# # estimate population size by age group in months (under 2 years) 
# population_data_age_months <- gen_population(total_pop = 5000,
#   groups = c("0-5", "6-8", "9-11","12-24"),
#   proportions = c(0.0164, 0.088, 0.088, 0.034),
#   strata = NULL) %>%
#   rename(age_group_mon = groups,
#     population = n)


## estimate population size by age categories in months and years 
# population_data_age_categories <- gen_population(total_pop = 5000,
#   groups = c("0-5 months", "6-8 months", "9-11 months","12-24 months",
#              "3-14 years", "15-29 years", "30-44 years", "45+ years"),
#   proportions = c(0.0164, 0.088, 0.088, 0.034, 
#                   0.3622, 0.276, 0.1616, 0.1321),
#   strata = NULL) %>%
#   rename(age_category = groups,
#     population = n)
# 
# 
# ## estimate population size by region proportion
# population_data_region <- gen_population(total_pop = 5000,         # set the total population
#   groups = c("Village A", "Village B", "Village C", "Village D"),  # set the groups
#   proportions = c(0.221, 0.174, 0.355, 0.245),                     # set the proportions for each group
#   strata = NULL) %>%               # do not stratify by gender
#   rename(patient_origin = groups,  # rename columns (syntax is NEW NAME = OLD NAME)
#          population = n)

```




<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// lab_data \\\
--------------------------------------------------------------------------------

If you have lab data, you can read it in here.
You will need to make sure that your identifier (e.g. case_number) matches the
linelist.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r lab_data}
## Read data -------------------------------------------------------------------
## This step reads in your lab data from Excel.
## You may need to rename your columns.

## to read in a specific sheet use "which"
# lab_data <- rio::import(here::here("labresults.xlsx"), which = "Sheet1")

## repeat preparation steps as appropriate
## make sure your unique ID variable name matches!

## Merge lab to linelist -------------------------------------------------------
# linelist_cleaned <- left_join(linelist_cleaned, lab_results, 
#                               by = "case_number")
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// browse_data \\\
--------------------------------------------------------------------------------

You'll want to look at your data. Here are a few ways you can explore.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r browse_data, eval = FALSE}
# ## view the first ten rows of data
# head(linelist_cleaned, n = 10)
# 
# ## view your whole dataset interactively (in an excel style format)
# View(linelist_cleaned)
# 
# ## overview of variable types and contents
# str(linelist_cleaned)
# 
# ## get summary: 
# ## mean, median and max values of numeric variables
# ## counts for categorical variables
# ## also gives number of NAs
# summary(linelist_cleaned)
# 
# ## view unique values contained in variables 
# ## you can run this for any column -- just replace the column name
# unique(linelist_cleaned$patient_origin) 

```




<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// remove_unused_data \\\
--------------------------------------------------------------------------------

Your data might have empty rows or columns you want to remove.
You can also use this section to create temporary datasets so you can review
specific variables or rows.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r remove_unused_data}
## Drop unused rows  -----------------------------------------------------------
## This step removes blank rows that don't have both a date of admission and ID
## number.

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


## Drop columns ----------------------------------------------------------------
## OPTIONAL: This step shows you how you can remove certain variables.

# linelist_cleaned <- linelist_cleaned %>%
#   select(-c(msf_involvement, treatment_facility_site))


## OPTIONAL: if you want to inspect certain variables, you can select these by
## name or column number. This example creates a reduced dataset for the first
## three columns, age_years, and sex.

# linelist_reduced <- select(linelist_cleaned, c(1:3, "age_years", "sex"))

```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
This part of the script will create and clean variables in your data.

All your cleaning and variable creation should happen in these chunks.
That way, in case something goes wrong, you can push the small arrow at the top
of the chunk to re-run all the code chunks up to the current one.

The chunks are:
- standardise_dates -- will set up and clean dates.
- create_age_group  -- creates the age group variables from age
- create_vars       -- creates variables based on other variables
- factor_vars       -- helps clean factor variables
- vector_vars       -- creates groups of variables for easy use


You must adapt this section according to your data!
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->


<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// standardise_dates \\\
--------------------------------------------------------------------------------

This chunk will help you set up and clean your date variables.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r standardise_dates, warning = FALSE, message=FALSE}
## DHIS2 standard data ---------------------------------------------------------
## If you got your data from DHIS2, use this portion of the code.
## If not, comment this section out and use the below.
# 
# ## make sure all date variables are formatted as dates 
# DATEVARS <- filter(linelist_dict, data_element_valuetype == "DATE") %>%
#   filter(data_element_shortname %in% names(linelist_cleaned)) %>% 
#   ## filter to match the column names of your data
#   pull(data_element_shortname) 
# 
# ## change to dates 
# linelist_cleaned <- linelist_cleaned %>%
#   mutate_at(DATEVARS, linelist::guess_dates,
#             error_tolerance = 0.5)



## Non-DHIS2 data --------------------------------------------------------------
## Use this section if you did not have DHIS2 data. 

## use the guess_dates() function to make a first pass at date variables.
 linelist_cleaned <- linelist_cleaned %>%
   mutate_at(vars(matches("date|Date")), linelist::guess_dates,
           error_tolerance = 0.5)


# ADDED CODE: Individually convert other date columns
linelist_cleaned$dtdeath <- linelist::guess_dates(linelist_cleaned$dtdeath)

linelist_cleaned$dthospitalisation <- linelist::guess_dates(linelist_cleaned$dthospitalisation)



## once you have run guess_dates(), take a look at your date variables.
## here is an example:
# table(linelist_cleaned$date_of_consultation_admission)


## Some dates will be unrealistic or wrong.
## Here is an example of how to manually fix dates. 
## Look at your data and edit as needed.
linelist_cleaned <- mutate(linelist_cleaned,
                           date_of_onset = case_when(
                                 date_of_onset < as.Date("2016-04-01")  ~ as.Date(NA),
                                 date_of_onset > as.Date("2017-10-31")  ~ as.Date(NA),
                                 TRUE                                   ~ date_of_onset
                                 ))


## Create epiweek variable -----------------------------------------------------
## This step creates an epiweek variable from the date of onset.
## You can use date_of_consultation_admission if you are missing many date_of_onset.

linelist_cleaned$epiweek <- aweek::date2week(linelist_cleaned$date_of_onset, 
                                             floor_day = TRUE, 
                                             factor = TRUE)

```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// create_age_group \\\ 
--------------------------------------------------------------------------------

This chunk will help you set up your age group variable.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r create_age_group}

## Age group variables ----------------------------------------------------------
## This step shows you how to create categorical variables from numeric variables.
## We have some intermediate steps on the way.

#ADDED CODE:
# Convert variable age_years to numeric class
linelist_cleaned$age_years <- as.numeric(linelist_cleaned$age_years)

## OPTIONAL: add under 2 years to the age_years variable
## data dictionary defines that under 2s dont have year filled in (but months/days instead)
# linelist_cleaned <- linelist_cleaned %>% 
#   mutate(age_months = case_when(
#     is.na(age_years) & is.na(age_months) ~ as.numeric(age_days / 30), 
#     TRUE                                 ~ as.numeric(age_months)
#   ),
#     age_years = case_when(
#     is.na(age_years) & is.na(age_months) ~ as.numeric(age_days / 365.25),
#     is.na(age_years)                     ~ as.numeric(age_months / 12),
#     TRUE                                 ~ as.numeric(age_years)
#   ))


## OPTIONAL: change those who are above or below a certain age to NA
# linelist_cleaned <- mutate(linelist_cleaned,
#                            age_years = case_when(
#                              is.na(age_years) ~ NA_integer_,
#                              age_years < 0   ~ NA_integer_,
#                              age_years > 120 ~ NA_integer_,
#                              TRUE            ~ as.integer(age_years) #Preserves other values
#                            ))

## OPTIONAL: create an age_months variable from decimal years variable
 linelist_cleaned <- mutate(linelist_cleaned,
                            age_months = case_when(
                              age_years < 5 ~ age_years * 12
                              )) 


## create age group variable for under 5 years based on months
linelist_cleaned$age_group_mon <- age_categories(linelist_cleaned$age_months, 
                                                 breakers = c(0, 6, 9, 12, 24), 
                                                 ceiling = TRUE)

## create an age group variable by specifying categorical breaks
linelist_cleaned$age_group <- age_categories(linelist_cleaned$age_years, 
                                             breakers = c(0, 3, 15, 30, 45))

## alternatively, create an age group variable by specifying a sequence
## this creates age groups for every ten years between 0 and 100.
# linelist_cleaned$age_group <- age_categories(linelist_cleaned$age,
#                                              lower = 0, 
#                                              upper = 100, 
#                                              by = 10)

## If you already have an age group variable defined, you should manually
## arrange the categories in order.
# linelist_cleaned$age_group <- factor(linelist_cleaned$age_group,
#                                      c("0-4y", "5-14y", "15-29y", "30-44y", "45+y"))


## to combine different age categories use the following function 
## this prioritises the smaller unit, i.e. if given months and years, will return months first
## generally, these delineations are NOT used for AJS.
##
# linelist_cleaned <- group_age_categories(linelist_cleaned, 
#                                          years = age_group, 
#                                          months = age_group_mon)
## drop the 0-2 years age level! (need to fix group_age_categories)
# linelist_cleaned <- linelist_cleaned %>% 
#   mutate(age_category = factor(age_category, 
#                                 levels = c("0-5 months",
#                                            "6-8 months",
#                                            "9-11 months",
#                                            "12-24 months",
#                                            "3-14 years",
#                                            "15-29 years",
#                                            "30-44 years", 
#                                            "45+ years")))
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// create_vars \\\
--------------------------------------------------------------------------------

This chunk will help you construct new variables from other variables. It
includes numeric, factor, and character vectors.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r create_vars, warning=FALSE, message=FALSE}

## Numeric variables -----------------------------------------------------------
## This step creates a variable for number of days under observation.
## You can adapt this step to create other calculated variables

## create number of days under observation
#linelist_cleaned <- mutate(linelist_cleaned, 
#                           obs_days = as.numeric(date_of_exit - 
#                                                   date_of_consultation_admission))


#ADDED CODE:
# 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"
  ))
     




## Factor (categorical) variables ----------------------------------------------

## This step creates a variable from another character/factor variable
## You can adapt this step to create other calculated variables

## The variable DIED will be binary (TRUE/FALSE) -- if patient died or not
## if your variable is named ("Dead on arrival", "Dead in facility"), use this
linelist_cleaned$DIED <- str_detect(linelist_cleaned$exit_status2, "Dead")
  
## use this version for coded values (e.g. "DOA", "DD")
# linelist_cleaned$DIED <- linelist_cleaned$exit_status %in% c("DD", "DOA")




## Recode character variables -------------------------------------------------
## This step shows how to fix misspellings in the geographic region variable.
## Ideally, you want these values to match your shapefile and population data!

# linelist_cleaned <- linelist_cleaned %>%
#   mutate(patient_origin = case_when(
#     patient_origin == "Valliages D"       ~ "Village D",
#     patient_origin == "VillageD"          ~ "Village D",
#     patient_origin == "Town C"            ~ "Village C"
#   ))



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



## sometimes, coding is inconsistent across variables -- 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 variables 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 variables
# # create list of variables
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 options
  linelist_cleaned <- linelist_cleaned %>%  
       mutate_at(vars(change_test_vars), forcats::fct_recode,   
       Positive = "Pos",
       Positive = "pos",
       Positive = "yes",
       Positive = "Yes",
       Positive = "Weakly pos",
       Positive = "1",
       Negative = "Neg",
       Negative = "neg",
       Negative = "no",
       Negative = "No",
       Negative = "0"
     )  
  



## Create a variable based on rules from other simple character variables
## If you have access to lab results, you can create a case definition variable 
## the tilda (~) is used to assign the new values (Conf, prob, susp, unknown)
## starting from the specific to the general
## TRUE assigns all remaining rows 
## 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_vars \\\
--------------------------------------------------------------------------------

This chunk will help you clean factor variables.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r factor_vars}

# Force missing values to NA
# important for sex to generate age pyramids
linelist_cleaned$sex <- fct_recode(linelist_cleaned$sex, 
                                   NULL = "Unknown/unspecified")

# 
# # change the order of levels in a single categorical variable 
# # This orders the levels -- since in all figures / tables, 0-4 should come
# # before 4-24, etc
# linelist_cleaned <- linelist_cleaned %>% 
#   mutate(time_to_death = factor(time_to_death, 
#                                 levels = c("0-4 hours", 
#                                            ">4-24 hours", 
#                                            ">24-48 hours", 
#                                            ">48 hours")))
# 
# 
# 
# # Change the order of levels of multiple categorical variables at the same time
# linelist_cleaned <- linelist_cleaned %>%
#   mutate_at(vars(
#            # Looks for variables beginning with "test"
#             starts_with("test")),          
#             fct_relevel,
#             # Sets order of levels
#             "Positive", "Negative", "Not done")




# ADDED CODE:
linelist_cleaned$patient_origin <- toupper(linelist_cleaned$patient_origin)

```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// vector_vars \\\
--------------------------------------------------------------------------------

This chunk creates groups of variables that you might want to use together. That
way, if you want to run the same function over these variables, you can run it
all at once. For example, you may want to look at frequency of all symptoms at
the same time.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r vector_vars}


# create a grouping of all symptoms 
SYMPTOMS <- c("history_of_fever", 
              "fever",
              "nausea_anorexia",
              "vomiting",
              "epigastric_pain_heartburn",
              "generalized_itch",
              "headache",
              "joint_pains",
              "diarrhoea",
              "bleeding", 
              "convulsions",
              "mental_state",
              "other_symptoms"  
              )


# 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"
          )
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// report_setup \\\
--------------------------------------------------------------------------------

This chunk removes cases after your reporting week and defines the start and end
of the reporting period.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r report_setup}
# return the last day of the reporting week
obs_end   <- week2date(str_glue(reporting_week, "-7"))

# filter out cases after end of reporting week
linelist_cleaned <- linelist_cleaned %>% 
  filter(date_of_consultation_admission <= obs_end | is.na(date_of_consultation_admission))

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

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



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// save_cleaned_data \\\
--------------------------------------------------------------------------------

You can save your cleaned dataset as an Excel. 
This automatically names your file "linelist_cleaned_DATE", where DATE is the
current date.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r save_cleaned_data}
# rio::export(linelist_cleaned, paste0("linelist_cleaned_", Sys.Date(), ".xlsx"))
```






<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
After adjusting and running the above code, you will have a clean dataset.
This marks the start of the ANALYSIS portion of the template.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->



### Person
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
This section focuses on who is affected - total, by sex, by age.
There is code to include:
- A bar chart of case numbers of incidence by age group.
- Attack rate (AR)
- Number of deaths (in suspected and confirmed cases)
- Mortality rates
- Case fatality ratio (CFR)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->



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 most affected age group was `r tab_linelist(linelist_cleaned, age_group) %>% slice(which.max(n)) %>% pull(value)` years. 



#### Demographics 

Cases by age group and definition 
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// describe_by_age_group_and_def \\\
--------------------------------------------------------------------------------

This chunk will create a table of cases by age group and case definition.
You can only use this chunk if you have lab results and have defined a case
definition variable.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r describe_by_age_group_and_def, warning=FALSE, message=FALSE}

tab_linelist(linelist_cleaned, 
             # rows as age groups, with case_def as strata (columns)
             age_group, strata = case_def, 
             # we want row totals and column totals
             col_total = TRUE, row_total = TRUE) %>%
  
  # removes the column called variable (which just says "age_group" in this case)
  select(-variable) %>%
  
  # renames the column called value to Age group
  rename("Age group" = value) %>%
  
  # there are several columns called "proportion" - renames them to %
  rename_redundant("%" = proportion) %>%
  
  # any time there is a column that says n, we replace that with "cases (n)"
  augment_redundant(" cases (n)" = " n$") %>%
  
  # makes the table easy to read with only 2 decimal places
  kable(digits = 1)
```


  
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// total_props_agegroup_sex \\\
--------------------------------------------------------------------------------

You can also show cases by proportion of the total population by sex and age.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
Cases by age group and sex 
```{r total_props_agegroup_sex, warning=FALSE, message=FALSE}
tab_linelist(linelist_cleaned, 
             age_group, strata = sex,
             col_total = TRUE, row_total = TRUE, prop_total = TRUE) %>% 
  select(-variable) %>%
  rename("Age group" = value) %>%
  rename_redundant("%" = proportion) %>%
  augment_redundant(" cases (n)" = " n$") %>%
  kable(digits = 1)
```

  
Age pyramid 

There were `r fmt_count(linelist_cleaned, is.na(sex)| sex == "Unknown/unspecified")` cases missing information on sex and `r fmt_count(linelist_cleaned, is.na(age_group))` missing age group.



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// age_pyramid \\\
--------------------------------------------------------------------------------

This chunk creates an age/sex pyramid of your cases.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r age_pyramid, warning=FALSE, message=FALSE}
 plot_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
       )
```



Of the patients, `r fmt_count(linelist_cleaned, patient_facility_type == "Outpatient")` were seen as outpatients and `r fmt_count(linelist_cleaned, patient_facility_type == "Inpatient")` were inpatients.

Age pyramid for under 2  

```{r age_pyramid_under2, warning=FALSE, message=FALSE}
filter(linelist_cleaned, age_years <= 2) %>%
 plot_age_pyramid(age_group = "age_group_mon", 
                  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
       )
```



Cases by symptoms
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// describe_by_symptoms \\\
--------------------------------------------------------------------------------

This chunk gives the counts and proportions for all of the variables in SYMPTOMS
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r describe_by_symptoms, warning=FALSE, message=FALSE}

tab_linelist(linelist_cleaned, SYMPTOMS, keep = "Oui") %>% 
  select(-value) %>%
  # 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
  mutate(variable = str_to_sentence(str_replace_all(variable, "_", " "))) %>%
  # rename accordingly
  rename_redundant("%" = proportion) %>%
  augment_redundant(" (n)" = " n$") %>%
  kable(digits = 1)
```



Cases by lab results 
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// describe_by_labs \\\
--------------------------------------------------------------------------------


This chunk gives the counts and proportions for all of the variables in LABS.
If you do not have lab results, comment out this chunk.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r describe_by_labs, warning=FALSE, message=FALSE}

tab_linelist(linelist_cleaned, LABS, 
             transpose = "value") %>% 
  # fix the way lab test names 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
  mutate(variable = str_to_sentence(str_replace_all(variable, "_", " "))) %>%
  # rename accordingly
  rename("Lab test" = variable) %>%
  rename_redundant("%" = proportion) %>%
  augment_redundant(" (n)" = " n$") %>%
  kable(digits = 1)

```




#### Case fatality ratio 

Of `r fmt_count(linelist_cleaned, patient_facility_type == "Inpatient")` inpatients, there have been `r fmt_count(linelist_cleaned, str_detect(exit_status2, "Dead"), patient_facility_type == "Inpatient")` deaths, of which 
`r fmt_count(linelist_cleaned, exit_status == "Dead on arrival")` 
were dead on arrival. 

 
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// describe_time_to_death \\\
--------------------------------------------------------------------------------

If you have the time_to_death available, you can use this chunk to describe when
inpatients died. If not, comment it out.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r describe_time_to_death, warning = FALSE, message = FALSE}
# filter(linelist_cleaned, DIED, patient_facility_type == "Inpatient") %>% 
#   tab_linelist(time_to_death, col_total = TRUE) %>% 
#   select(-variable) %>%
#   rename("Time (hours)" = value, 
#          "Deaths (n)" = n, 
#          "%" = proportion) %>% 
#   kable(digits = 1)
```



The case fatality ratio among inpatients with known outcomes is below.
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// overall_cfr \\\
--------------------------------------------------------------------------------

This chunk gives the case fatality ratio among inpatients with outcomes.
If you have no deaths reported, the table will not be useful.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r overall_cfr, warning = FALSE, message = FALSE}
# use arguments from above to produce overal CFR
linelist_cleaned %>% 
  filter(patient_facility_type == "Inpatient") %>%
  case_fatality_rate_df(deaths = DIED, mergeCI = TRUE) %>%
  rename("Deaths" = deaths,
         "Cases" = population,
         "CFR (%)" = cfr,
         "95%CI" = ci) %>%
  kable(digits = 1)

```


The case fatality ratio by sex among inpatients with known outcomes is below. 
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// cfr_by_sex \\\
--------------------------------------------------------------------------------

This chunk gives the case fatality ratio among inpatients with outcomes, divided
by sex. If you have no deaths reported, the table will not be useful.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r cfr_by_sex, warning = FALSE, message = FALSE}
linelist_cleaned %>%
  filter(patient_facility_type == "Inpatient") %>%
  mutate(sex = forcats::fct_explicit_na(sex, "-")) %>%
  case_fatality_rate_df(deaths = DIED, group = sex, mergeCI = TRUE, add_total = TRUE) %>%
  rename("Sex" = sex, 
         "Deaths" = deaths, 
         "Cases" = population, 
         "CFR (%)" = cfr, 
         "95%CI" = ci) %>% 
  knitr::kable(digits = 1)
```



CFR by age group among inpatients with known outcomes
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// cfr_by_age_group \\\
--------------------------------------------------------------------------------

This chunk gives the case fatality ratio among inpatients with outcomes, divided
by age group. If you have no deaths reported, the table will not be useful.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r cfr_by_age_group, warning = FALSE, message = FALSE}

linelist_cleaned %>%
  filter(patient_facility_type == "Inpatient") %>%
  case_fatality_rate_df(deaths = DIED, group = age_group, mergeCI = TRUE, add_total = TRUE) %>%
  tidyr::complete(age_group, 
                  fill = list(deaths = 0, 
                              population = 0, 
                              cfr = 0, 
                              ci = 0)) %>% # Ensure all levels are represented
  rename("Age group" = age_group, 
         "Deaths" = deaths, 
         "Cases" = population, 
         "CFR (%)" = cfr, 
         "95%CI" = ci) %>% 
  knitr::kable(digits = 1)
```


<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// cfr_by_case_def \\\
--------------------------------------------------------------------------------

This chunk gives the case fatality ratio among inpatients with outcomes, divided
by case definition. You must have lab results or a case definition defined to
use this. If you have no deaths, the table will not be useful.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r cfr_by_case_def, warning = FALSE, message = FALSE}
linelist_cleaned %>%
  filter(patient_facility_type == "Inpatient") %>%
  case_fatality_rate_df(deaths = DIED, group = case_def,
                        mergeCI = TRUE, add_total = TRUE) %>%
  rename("Case definition" = case_def,
         "Deaths" = deaths,
         "Cases" = population,
         "CFR (%)" = cfr,
         "95%CI" = ci) %>%
  knitr::kable(digits = 1)
```



#### Attack rate

The attack rate per 10,000 population is below (based on available population data available for the catchment area/region of interest). 



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// collect_variables \\\
--------------------------------------------------------------------------------

This defines your population by summing up the population by age group.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r collect_variables, warning=FALSE, message=FALSE}
# define population 
population <- sum(population_data_region$population)
```



Below gives the attack rate per 10,000 population (N = `r format(population, big.mark = ",")`)

<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
///attack_rate \\\
--------------------------------------------------------------------------------

This chunk calculates the attack rate and then shows you the confidence
intervals.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r attack_rate, warning=FALSE, message=FALSE}
ar <- attack_rate(nrow(linelist_cleaned), population, multiplier = 10000)

ar %>%
  merge_ci_df(e = 3) %>% # merge the lower and upper CI into one column
  rename("Cases (n)" = cases, 
         "Population" = population, 
         "AR (per 10,000)" = ar, 
         "95%CI" = ci) %>% 
  select(-Population) %>% # drop the population column as it is not changing
  knitr::kable(digits = 1, align = "r")
```

Here, we can see that the attack rate for a population of `r format(population, big.mark = ",")` was `r fmt_ci_df(ar, percent = FALSE)`.


<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// attack_rate_by_agegroup \\\
--------------------------------------------------------------------------------

This chunk calculates the attack rate by age group and then gives a table of
attack rate by group.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r attack_rate_by_agegroup, warning=FALSE, message=FALSE}

# cases <- count(linelist_cleaned, age_group) %>%    # cases for each age_group
#   left_join(population_data_age, by = "age_group") # merge population data 
# 
# 
# # attack rate for each group
# attack_rate(cases$n, cases$population, multiplier = 10000, mergeCI = TRUE) %>% 
#   # add the epiweek column to table
#   bind_cols(select(cases, age_group), .) %>% 
#   rename("Age group" = age_group, 
#          "Cases (n)" = cases, 
#          "Population" = population, 
#          "AR (per 10,000)" = ar, 
#          "95%CI" = ci) %>% 
#   kable(digits = 1, format.args = list(big.mark = ",")) # set thousands separator
```




<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#### Mortality attributable to AJS
This section can only be used if you are in a closed population (eg refugee camp).
The assumptions don't hold in an open/community setting.

This chunk calculates the attack rate by age group and then gives a table of
attack rate by group.

This section gives mortality rates attributable to AJS in a closed population. 
It does not calculate all-cause mortality. It assumes that all AJS deaths are 
among inpatients.

This demonstrates three ways of calculating mortality rate based on 
catchment population (twice) and based on hospital population. 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// collect_variables_rates \\\
--------------------------------------------------------------------------------

This chunk calculates key variables to do mortality rate.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r collect_variables_rates, warning=FALSE, message=FALSE}
# 
# # count number of deaths 
# deaths <- sum(linelist_cleaned$DIED)
# 
# # outbreak duration in days 
# obs_time <- as.numeric(obs_end - obs_start)
# 
# # patient observation time 
# pat_obs_time <- linelist_cleaned %>% 
#   filter(!is.na(exit_status)) %>% 
#   summarise(days = sum(obs_days)) %>%
#   pull(days)
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// mortality_rate \\\
--------------------------------------------------------------------------------

To produce a mortality rate (attribuatble to AJS) per 10,000 people use the
following code chunk. This assumes that you are capturing every AJS death in
your population.
Mortality rate attributable to AJS per 10,000 population
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->

```{r mortality_rate, warning=FALSE, message=FALSE}
# 
# mortality_rate(deaths, population, multiplier = 10000, mergeCI = TRUE) %>%
#   rename("Deaths" = deaths, 
#          "Population" = population, 
#          "Mortality (per 10,000)" = `mortality per 10 000`, 
#          "95%CI" = ci) %>% 
#   kable(digits = 1)
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// mortality_rate_CMR \\\
--------------------------------------------------------------------------------

To produce a crude mortality rate attributable to AJS per 10,000 people per 
day, use the folowing code chunk. 

This assumes that you are capturing every AJS death in your population and 
that your population remains stable over the time period of interest. 

In this situation the time period of interest is from the beginning of the 
epiweek your first case occured in, until the last day of the epiweek you are 
currently reporting on. (see this presentation (https://www.odi.org/sites/odi.org.uk/files/odi-assets/events-presentations/776.pdf) for more detail)

Crude mortality rate attributable to AJS per 10,000 population per day
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->


```{r mortality_rate_CMR, warning=FALSE, message=FALSE}
# 
# mortality_rate(deaths, population*obs_time, multiplier = 10000, mergeCI = TRUE) %>%
#   rename("Deaths" = deaths, 
#          "Person-days" = population, 
#          "Mortality (per 10,000/day)" = `mortality per 10 000`, 
#          "95%CI" = ci) %>% 
#   kable(digits = 1)
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// mortality_rate_patients \\\
--------------------------------------------------------------------------------

Alternatively, if you are unsure whether your hospital deaths are
representative of the wider population, use the following code chunk. 
This uses the person days of cases in your linelist with a known outcome. 
However, this will give you an unreasonably high mortality rate, as 
those in hospital will only be the most severely affected. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
Mortality rate attributable to AJS per 10,000 patients per day

```{r mortality_rate_patients, warning=FALSE, message=FALSE}
# mortality_rate(deaths, pat_obs_time, multiplier = 10000, mergeCI = TRUE) %>%
#   rename("Deaths" = deaths, 
#          "Person-days" = population, 
#          "Mortality (per 10,000/day)" = `mortality per 10 000`, 
#          "95%CI" = ci) %>% 
#   kable(digits = 1)
```






### Time
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
This section focuses on the when of the outbreak:
- When did cases fall ill?
- Are numbers increasing or decreasing?

There is code to include an epi curve.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->

There were `r fmt_count(linelist_cleaned, is.na(date_of_onset))` cases missing dates of onset.

<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// epicurve \\\
--------------------------------------------------------------------------------

# This chunk will calculate weekly incidence and plot your epicurve.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r epicurve, warning=FALSE, message=FALSE}

# count cases for each week of outbreak
# you can change the start of your week (eg "Sunday week", "Saturday week")
inc_week_7 <- incidence(linelist_cleaned$date_of_onset, interval = "Monday week")


# plot your epicurve
basic_curve <- plot(inc_week_7, show_cases = TRUE, border = "black", n_breaks = nrow(inc_week_7)) + 
  scale_y_continuous(expand = c(0, 0)) +  # set origin for axes
  # add labels to axes and below chart
  epicurve_labels +
  # change visuals of dates and remove legend title
  epicurve_theme

# show your plot (stored for later use) 
basic_curve + scale_x_incidence(inc_week_7, n_breaks = 6)

## if the outbreak has been going on for a while, your x-axis might look messy.
## to reduce the number of labels, uncomment the below.
## you can customize the number of breaks by changing n_breaks
## you can also use months for your break (here in three-month intervals)

# basic_curve + scale_x_incidence(inc_week_7, n_breaks = 6)
# basic_curve + scale_x_date(breaks = "3 months", date_labels = "%b %Y")


```

The peak of the outbreak was in `r date2week(find_peak(inc_week_7), floor_day = TRUE)`



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// biweekly_epicurve \\\
--------------------------------------------------------------------------------

If needed, you can plot a biweekly epicurve. You can also do this by month,
quarter, or year. You can change the start of your week as needed.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->

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

# create case counts for two week groups of your outbreak overall
inc_week_14 <- incidence(linelist_cleaned$date_of_onset, interval = "2 Monday weeks")


# plot your epicurve
plot(inc_week_14, show_cases = TRUE, border = "black", n_breaks = nrow(inc_week_14)) + 
  scale_y_continuous(expand = c(0,0)) +  # set origin for axes
  # add labels to axes and below chart
  epicurve_labels +
  # change visuals of dates and remove legend title
  epicurve_theme

```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// incidence_by_gender \\\
--------------------------------------------------------------------------------

You can plot weekly incidence by gender using this chunk.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->

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

# get counts by gender
inc_week_7_gender <- incidence(linelist_cleaned$date_of_onset, 
                               interval = "Monday week", 
                               groups = linelist_cleaned$sex)

plot(inc_week_7_gender, 
     show_cases = FALSE, # this removes black borders
     border = "black", 
     n_breaks = nrow(inc_week_7_gender)) + 
  scale_y_continuous(expand = c(0,0)) +  # set origin for axes
  epicurve_labels +
  epicurve_theme + scale_x_incidence(inc_week_7, n_breaks = 6)

```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// incidence_by_case_def \\\
--------------------------------------------------------------------------------

This chunk shows the weekly incidence by case definition. You must have a
defined case definition variable.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r incidence_by_case_def, warning=FALSE, message=FALSE}
inc_week_7_casedef <- incidence(linelist_cleaned$date_of_onset, 
                                interval = "Monday week", 
                                groups = linelist_cleaned$case_def)
plot(inc_week_7_casedef, show_cases = TRUE, border = "black", n_breaks = nrow(inc_week_7_casedef)) + 
  scale_y_continuous(expand = c(0, 0)) +  # set origin for axes
  # add labels to axes and below chart
  epicurve_labels +
  # change visuals of dates, remove legend title and move legend to bottom
  epicurve_theme + scale_x_incidence(inc_week_7, n_breaks = 6)
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// incidence_by_water_source \\\
--------------------------------------------------------------------------------

You can similarly stratify your incidence by water source (or any other
categorical variable!).
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->

```{r incidence_by_water_source, warning=FALSE, message=FALSE}
#inc_week_7_water <- incidence(linelist_cleaned$date_of_onset, 
#                              interval = "Monday week", 
#                              groups = linelist_cleaned$water_source)


#plot(inc_week_7_water, show_cases = TRUE, border = "black", 
#     n_breaks = nrow(inc_week_7_water)) + 
#  scale_y_continuous(expand = c(0,0)) +  # set origin for axes
#  epicurve_labels +
#  epicurve_theme + scale_x_incidence(inc_week_7, n_breaks = 6)
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// incidence_by_sex_facility \\\
--------------------------------------------------------------------------------

This chunk limits the data to only inpatients, then stratifies by sex.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
Cases by week of onset among inpatients by sex 

```{r incidence_by_sex_facility, warning=FALSE, message=FALSE}
inc_week_7_sex_fac <- linelist_cleaned %>%
  filter(patient_facility_type == "Inpatient") %>%
  with(incidence(date_of_onset, interval = "Monday week", groups = sex))



plot(inc_week_7_sex_fac, show_cases = TRUE, border = "black", n_breaks = nrow(inc_week_7_sex_fac)) + 
  scale_y_continuous(expand = c(0, 0)) +  # set origin for axes
  epicurve_labels +
  epicurve_theme + scale_x_incidence(inc_week_7, n_breaks = 6)

```




<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// attack_rate_per_week \\\
--------------------------------------------------------------------------------

This chunk creates a table of attack rate by week
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
Attack rate per 10,000 population by week 

```{r attack_rate_per_week, warning = FALSE, message = FALSE}
# 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) %>% 
  knitr::kable(digits = 1, align = "r")
```
 


<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// cumulative_attack_rate_per_week \\\
--------------------------------------------------------------------------------

This chunk calculates the cumulative attack rate per week.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
Cumulative attack rate per 10,000 population per week

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

attack_rate(cases$cumulative, population, multiplier = 10000) %>% 
  bind_cols(select(cases, epiweek), .) %>% # add the epiweek column to table
  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) %>% 
  knitr::kable(digits = 1, align = "r")
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// cumulative_attack_rate_per_week \\\
--------------------------------------------------------------------------------

This chunk calculates the CFR among inpatients as a proportion per week.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
Case fatality ratio as a proportion among inpatients by week

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

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


cfr %>%
  merge_ci_df(e = 4) %>% # merge the lower and upper CI into one column
  rename("Epiweek" = epiweek, 
         "Deaths" = deaths, 
         "Cases" = population, 
         "CFR (%)" = cfr, 
         "95%CI" = ci) %>% 
  knitr::kable(digits = 1, align = "r")
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
You could plot the AR (in the population) and CFR (among inpatients only)
together with the epicurve by epiweek.

We will do this in three steps: 
- creating the AR graph (ar_line_graph)
- creating the CFR graph (cfr_line_graph)
- combining and plotting 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->


<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// ar_line_graph \\\
--------------------------------------------------------------------------------

This chunk sets up the AR graph. (It does not print the graph automatically.)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r ar_line_graph, warning=FALSE, message=FALSE}

ar_plot <- ggplot(ar, aes(x = week2date(epiweek) + (7 * 0.5), 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, 0)) +
  # scale the x axis the same as the incidence curve. Expand forces it to align. 
  incidence::scale_x_incidence(inc_week_7, 
                               n_breaks = nrow(inc_week_7), 
                               expand = c(0, 7 * 1.5)) +
  # 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
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// cfr_line_graph \\\
--------------------------------------------------------------------------------

This chunk sets up the CFR graph. (It does not print the graph automatically.)
If you do not have any deaths in your data set, this will not show anything.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r cfr_line_graph, warning=FALSE, message=FALSE}

cfr_plot <- ggplot(cfr, aes(x = week2date(epiweek) + (7 * 0.5), 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, 0)) +  
  # scale the x axis the same as the incidence curve. Expand forces it to align. 
  incidence::scale_x_incidence(inc_week_7, 
                               n_breaks = nrow(inc_week_7), 
                               expand = c(0, 7 * 1.5)) +
  # 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 
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// epicurve_ar_cfr \\\
--------------------------------------------------------------------------------

This chunk prints the AR and CFR with the epicurve.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r epicurve_ar_cfr, warning = FALSE, message = FALSE, fig.height = 10}

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

cowplot::plot_grid(
  basic_curve + nofx,
  ar_plot + nofx,
  cfr_plot,
  align = "v", # align plots vertically
  axis = "lr", # only by their left and right margins
  ncol = 1     # allow only one column
)
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// describe_admissions_by_epiweek \\\
--------------------------------------------------------------------------------

This chunk creates a table showing inpatient admissions by  case definition 
and epiweek.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
Inpatient admissions by case definition and week  

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

linelist_cleaned %>% 
  filter(patient_facility_type == "Inpatient") %>%
  tab_linelist(epiweek, strata = case_def,
             col_total = TRUE,
             row_total = TRUE) %>% 
  select(-variable) %>%
  rename_redundant("%" = proportion) %>%
  augment_redundant(" (n)" = " n$") %>%
  kable(digits = 1)
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// describe_admissions_by_epiweek \\\
--------------------------------------------------------------------------------

This chunk creates a table showing inpatient exits by type per epiweek.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
Inpatient discharges by reason for exit and week 

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

linelist_cleaned %>% 
  filter(patient_facility_type == "Inpatient") %>%
  tab_linelist(epiweek, strata = exit_status2, 
               col_total = TRUE, row_total = TRUE) %>% 
  select(-variable) %>%
  rename("Week" = value) %>%
  rename_redundant("%" = proportion) %>%
  augment_redundant(" (n)" = " n$") %>%
  kable(digits = 1)
```






### Place 
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
This section focuses on the where of the outbreak: what area is affected, how
many villages, and so on.

There is code to include maps based on distribution of cases. You must have a
shapefile to create this map.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->


#### Descriptive



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// describe_by_region_facility \\\
--------------------------------------------------------------------------------

This chunk creates a basic table of the number and percent of cases by region
and by facility. If you only have one facility, you can remove the strata.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
Cases by region and facility type

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

tab_linelist(linelist_cleaned, patient_origin, 
             strata = patient_facility_type,
             col_total = TRUE, row_total = TRUE) %>% 
  select(-variable) %>%
  rename("Region" = value) %>%
  rename_redundant("%" = proportion) %>%
  augment_redundant(" (n)" = " n$") %>%
  kable(digits = 1)
```

  

<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// describe_by_region_outcome \\\
--------------------------------------------------------------------------------

This chunk creates a basic table of the outcomes among inpatients by region.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
Inpatient discharges by reason for exit and region 

```{r describe_by_region_outcome, warning = FALSE}

tab_linelist(filter(linelist_cleaned,
                    patient_facility_type == "Inpatient"),
             patient_origin, strata = exit_status2,
             col_total = TRUE, row_total = TRUE) %>% 
  select(-variable) %>%
  rename("Region" = value) %>%
  rename_redundant("%" = proportion) %>%
  augment_redundant(" (n)" = " n$") %>%
  kable(digits = 1)
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// attack_rate_by_region \\\
--------------------------------------------------------------------------------

If you do not have a shapefile, you may want to calculate attack rates by
region.

Consider facet wrapping by a larger unit if you have many regions (eg if you
have patients from 10+ villages, you may want to show your tables by health
zone)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
Attack rage per 10,000 population by region 

```{r attack_rate_by_region, warning = FALSE}

# count cases by region
cases <- count(linelist_cleaned, 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
  kable(digits = 1, align = "r", format.args = list(big.mark = ",")) # set thousands separator
```


<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// bar_attack_rate_by_region \\\
--------------------------------------------------------------------------------

You can then plot the above table on a bar plot. AR is on the y-axis, and it
will show regions in descending order by AR.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r bar_attack_rate_by_region, warning = FALSE}

ar_region <- filter(ar_region, !is.na(`AR (per 10,000)`))

ggplot(ar_region, 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 = paste0("Source: MSF data from ", reporting_week)) + 
  epicurve_theme
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// mortality_rate_region \\\
--------------------------------------------------------------------------------

You could also calculate mortality rate by region (check the mortality code
chunk in Person section for assumptions).
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
Mortality rate per 10,000 population by region 

```{r mortality_rate_region, warning = FALSE}

deaths <- group_by(linelist_cleaned, patient_origin) %>%
  filter(str_detect(exit_status2, "Dead")) %>% 
  summarise(deaths = n()) %>% # count deaths by region
  left_join(population_data_region, by = "patient_origin") # merge population data 

mortality_rate(deaths$deaths, deaths$population, multiplier = 10000) %>%
  # add the region column to table
  bind_cols(select(deaths, patient_origin), .) %>% 
  merge_ci_df(e = 4) %>% # merge the lower and upper CI into one column
  rename("Region" = patient_origin, 
         "Deaths" = deaths, 
         "Population" = population, 
         "Mortality (per 10,000)" = `mortality per 10 000`, 
         "95%CI" = ci) %>% 
  kable(digits = 1)
```




#### Maps 
<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// read_shapefiles \\\
--------------------------------------------------------------------------------

To create maps, you need to have a shapefile of the area. Often, the MSF GIS
unit can provide shapefiles.

Your shapefile can be a polygon or points. Polygons do not need to be contiguous.

The names of the polygons or points MUST match the names in your linelist.

Your coordinate reference system needs to be WGS84.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r read_shapefiles, message=FALSE}
## fake map data - DELETE if you are using real data
#map <- gen_polygon(regions = unique(linelist_cleaned$patient_origin))


## read in shapefile
# map <- read_sf(here::here("mapfolder", "region.shp"))
map <- read_sf(here("content", "training", "Walk-through", "mapfolder", "Quartiersshape.shp"))

## check the coordinate reference system (CRS)
st_crs(map)

## if CRS not WGS84, reset it
map <- st_set_crs(map, value = 4326) # Sets to WGS84

```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// dot_map \\\
--------------------------------------------------------------------------------

If you have the coordinates for your cases, you can create a dot map showing
where the cases are.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r dot_map, message=FALSE}
# ggplot() +
#   # load shapefile
#   geom_sf(data = map, col = "grey50") + 
#   # add in case coordinates
#   geom_point(mapping = aes(x = lat, y = lon), 
#              data = linelist_cleaned)
```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// choropleth_maps \\\
--------------------------------------------------------------------------------

Once you have loaded your shapefile, you can map the case counts or attack rates.

You will want to make your counts or AR categorical. R will not do this
automatically (unlike QGIS or ArcGIS).

This chunk will walk you through several steps:
- Create table with categories of counts or ARs by region.
- Join your table with your shapefile.
- Choose which variable you will use.

Make sure you delete or comment out the section you are not using.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r choropleth_maps, message = FALSE, warning = FALSE}

## Create a categories of ARs by region ----------------------------------------

# define maximum
max_ar    <- max(ar_region$`Upper 95%CI`, na.rm = TRUE) # define your highest AR

# 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 table to shapefile -----------------------------------------------------
mapsub <- left_join(map, ar_map, by = c("Name" = "Region"))



## Plot cases or AR by geography -----------------------------------------------

## 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") + 
  # remove coordinates and axes
  theme_void() 

```



<!-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/// map_for_loop_epiweek \\\
--------------------------------------------------------------------------------

This chunk will create a map and barplot by week.

You need to:
- Decide if you want to show counts, AR, or categories of counts/AR.
- Define categories for your variable.
- Replace `Cases (n)` and `AR (per 10,000)` or "categories", appropriately.


If you have a lot of map regions, you may want to use facet_wrap to show
sub-units.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -->
```{r map_for_loop_epiweek, message = FALSE, warning = FALSE, fig.width = 12}

## 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 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 = 100) # four breaks rounded to nearest 100
             ))


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


# go through each epiweek, fiter and plot the data
for (i in unique(cases$epiweek)) {
  this_ar <- filter(ar, epiweek == i)
  
  # 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 
  
  this_ar <- filter(this_ar, !is.na(`AR (per 10,000)`))

  barplot <- ggplot(this_ar, 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), limits = c(0, max_ar)) +  
    # add labels to axes and below chart
    labs(x = "Region", y = "AR (per 10,000)", 
         captions = paste0("Source: MSF data from ", reporting_week)) + 
    epicurve_theme
  
  
  
  # combine the barplot and map plot into one
  print(
    cowplot::plot_grid(
      barplot + labs(title = paste0("Epiweek: ", i)),
      map_plot,
      nrow = 1,
      align = "h",
      axis = "tb"
    )
  )
}
```

