title: “Report”
subtitle: “Ebola outbreak in Sierra Leone”
author: “Naga Deepthi Vempati”
output:
word_document: default
date: “r Sys.Date()”
params:
district:
input: select
value: “West II” # can give any other site location
choices: [“West I”, “West II”, “West III”, “Mountain Rural”]
publish_date: !r lubridate::ymd(“2014-12-01”)
knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE)
# This chunk loads packages
# This chunk loads packages
pacman::p_load(
rio, # for importing data
here, # for locating files
janitor, # for data cleaning
lubridate, # for date cleaning
epikit, # for easy inline code
scales, # for percents and date labels
officer, # border lines
flextable, # make HTML tables
gtsummary, # for nice tables
apyramid, # for age/sex pyramids
RColorBrewer, # for colour palettes
gghighlight, # highlighting plot parts
ggExtra, # special plotting functions
viridis, # for color-blind friendly color scales
tsibble, # for epiweeks and other time series analyses
tidyverse # for data management and visualization
)
# Import data -------------------------------------------------------------
here("data", "raw", "surveillance_linelist_20141201.csv")
# surveillance dataset
surv_raw <- import(here("data", "raw", "surveillance_linelist_20141201.csv"))
# hospital datasets
hosp_central <- import(here("data", "raw", "hospitals", "20141201_hosp_central.csv"))
hosp_military <- import(here("data", "raw", "hospitals", "20141201_hosp_military.csv"))
hosp_other <- import(here("data", "raw", "hospitals", "20141201_hosp_other.csv"))
hosp_port <- import(here("data", "raw", "hospitals", "20141201_hosp_port.csv"))
hosp_smmh <- import(here("data", "raw", "hospitals", "20141201_hosp_smmh.csv"))
hosp_missing <- import(here("data", "raw", "hospitals", "20141201_hosp_missing.csv"))
# Import lab data
# (place this in the importing chunk of your script)
lab <- import(here("data", "raw", "lab_results_20141201.xlsx")) %>%
clean_names()
# In the import data code chunk
# Import lab data
# (add to the import code chunk of your R script)
investigations <- import(here("data", "raw", "case_investigations_20141201.xlsx")) %>%
# remove unnecessary columns
select(-c(age, age_unit, sex))
# Exploratory analysis -----------------------------------------
names(surv_raw)
# tabulate <<<<<<<<<<<<<
tabyl(surv_raw, sex)
# tabulate fever in the surv_raw data frame
tabyl(surv_raw, fever)
# a histogram of the case ages
ggplot(data = surv_raw, mapping = aes(x = age))+
geom_histogram()
# to clean the data wjicj has space in the variable
tabyl(surv_raw, `age unit`)
#print class------------------
#class(surv_raw)
#print class, age variable------------------
class(surv_raw$age)
#print class , onset date variable ------------------
class(surv_raw$`onset date`)
#print class, date of report -----------------
class(surv_raw$`date of report`)
#histogram for numeric analysis. wt variable ......
ggplot(data = surv_raw, mapping = aes(x = `wt (kg)`))+
geom_histogram()
#histogram for numeric analysis. ht variable ......
ggplot(data = surv_raw, mapping = aes(x = `ht (cm)`))+
geom_histogram()
# print cross tabulation of two district columns with spaces in the names .........
tabyl(surv_raw, adm3_name_res, adm3_name_det)
# print class of surv_raw data frame
class(surv_raw)
# print class of age column
class(surv_raw$age)
# print class of `onset date` column
class(surv_raw$`onset date`)
# print class of `date of report` column
class(surv_raw$`date of report`)
# plot histogram of `wt kg` column
ggplot(data = surv_raw, mapping = aes(x = `wt (kg)`))+
geom_histogram()
# plot histogram of `ht cm` column
ggplot(data = surv_raw, mapping = aes(x = `ht (cm)`))+
geom_histogram()
# print cross tabulation of two district columns
tabyl(surv_raw, adm3_name_res, adm3_name_det)
# print the class of a column
class(surv$sex)
# print the class of a column
class(surv$date_onset)
# dim() ...................
dim(surv)
dim(hosp)
# duplicates ..........
surv %>%
count(case_id) %>%
filter(n > 1)
hosp %>%
count(ID) %>%
filter(n > 1)
# for each dataset, sort and print the first 10 identifier values
# for the surveillance dataset
surv %>%
arrange(case_id) %>%
pull(case_id) %>%
head(10)
# for the hospital dataset
hosp %>%
arrange(ID) %>%
pull(ID) %>%
head(10)
# TO KNOW THE VARIABLE NAMES .
surv %>% names()
hosp %>% names()
# dim ''''''''''''
dim(test_join)
# cross tabulate the two columns in the joined dataset
test_join %>%
tabyl(sex.x, sex.y)
dim(combined)
dim(lab)
# In the testing area code chunk
# Compare data frame dimensions
# (write and run directly in the Console)
dim(investigations)
dim(combined)
# Compare case_id (our matching variable)
# (write and run directly in the Console)
# for the investigations dataset
investigations %>% arrange(case_id) %>% pull(case_id) %>% head(10)
# for the combined dataset
combined %>% arrange(case_id) %>% pull(case_id) %>% head(10)
timelines
class(timelines_long$date_type)
levels(timelines_long$date_type)
# A vector of dates
seq.Date(from = ymd("2014-05-06"),
to = ymd("2014-11-28"),
by = "week")
# Integrating into yhe epicurve .......
ggplot(data = combined, aes(x = date_onset))+
geom_histogram(
breaks = seq.Date(
from = ymd("2014-05-06"),
to = ymd("2014-11-28"),
by = "week")) +
labs(x = "Date of onset", y = "Incidence") +
theme_bw()
# Monday BEFORE the earliest case
floor_date(ymd("2014-05-06"), unit = "week", week_start = 1)
# Monday AFTER the last case
ceiling_date(ymd("2014-11-28"), unit = "week", week_start = 1)
# Sequence of Mondays from before earliest case, to after latest case
seq.Date(from = floor_date(ymd("2014-05-06"), unit = "week", week_start = 1),
to = ceiling_date(ymd("2014-11-28"), unit = "week", week_start = 1),
by = "week")
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
# Cleaning surveillance data -----------------------------------------
# using the pipe operator for cleaning ... creating a dataset ....
#surv <- surv_raw # make a copy of the raw data
# Cleaning command
# surv <- surv_raw %>% # the raw dataset
# clean_names() # standardize the column names
surv <- surv_raw %>%
# automatically clean column names
clean_names() %>%
# manually clean column names
rename(
date_onset = onset_date,
date_report = date_of_report,
district_res = adm3_name_res,
district_det = adm3_name_det) %>%
# remove unnecessary column
select(-row_num) %>%
# de-duplicate rows
distinct() %>%
# convert date_onset to date class
mutate(date_onset = mdy(date_onset)) %>%
mutate(date_report = mdy(date_report)) %>%
# convert age to numeric class
mutate(age = as.numeric(age)) %>%
# convert "Unknown" sex to NA
mutate(sex = na_if(sex, "Unknown")) %>%
# properly record missing values in many character columns
mutate(across(.cols = where(is.character), .fns = ~na_if(.x, ""))) %>%
# re-code hospital column
mutate(hospital = recode(hospital,
# for reference: OLD = NEW
"Mitilary Hospital" = "Military Hospital",
"Port" = "Port Hospital",
"Port Hopital" = "Port Hospital",
"St. Mark's Maternity Hospital (SMMH)" = "SMMH")) %>%
# recode sex
mutate(sex = recode(sex,
"m" = "male",
"f" = "female")) %>%
# convert negative weight values to NA
mutate(wt_kg = ifelse(wt_kg < 0, NA, wt_kg)) %>%
# create case definition
mutate(case_def = case_when( # create new column case_def
lab_confirmed == TRUE ~ "Confirmed", # Confirmed case if lab positive
epilink == "yes" & fever == "yes" ~ "Suspect", # else, Suspect case if epilink and fever
TRUE ~ "To investigate")) %>% # else, if none of previous are true, assign as "To investigate"
# create age-in-years
mutate(age_years = case_when(
age_unit == "months" ~ age/12, # if age is given in months
age_unit == "years" ~ age, # if age is given in years
is.na(age_unit) ~ age)) %>% # if unit missing assume years, else NA
# create age category column
mutate(age_cat = age_categories( # create new column
age_years, # numeric column to make groups from
lower = 0,
upper = 70,
by = 10)) %>%
# Make date-difference column
mutate(diff = date_report - date_onset) %>%
# New column will be TRUE if the two values are different
mutate(moved = district_res != district_det) %>%
# create new column that prioritizes district of detection
mutate(district = coalesce(district_det, district_res)) %>%
# remove cases “to investigate”
filter(case_def %in% c("Confirmed", "Suspect")) %>%
# re-arrange columns
select(case_id, starts_with("date"), diff, sex, age, age_unit, age_years, age_cat, hospital, district, district_res, district_det, moved, everything())
# Joining data
hosp <- bind_rows(hosp_central, hosp_port, hosp_military,
hosp_smmh, hosp_other, hosp_missing) %>%
# select specific columns from hosp, and re-name ID as case_ID
select(
case_id = ID, # select and rename
date_hospitalisation, # select
time_admission, # select
date_outcome, # select
outcome) # select
# ANTIJOIN WITH HOSP AS THE BASE DATAFRAME ..........
#anti_surv <- anti_join(hosp, surv, by = c("ID" = "case_id")) # Which rows in hosp are NOT present in surv?
# ANTIJOIN , USING SURV AS BASE LINE ......
#anti_hosp <- anti_join(surv, hosp, by = c("case_id" = "ID")) # Which rows in surv are not present in hosp?
# USING LEFT JOIN ..............
#test_join <- surv %>%
#left_join(hosp, by = c("case_id" = "ID"))
# Join the two data frames with a left join .
combined <- left_join(surv, hosp, by = "case_id")
# Join the two data frames with a left-join
combined <- left_join(combined, lab, by = "case_id")
# Join the two data frames with a left-join
# (add to the joining code chunk of your R script)
combined <- left_join(combined, investigations, by = "case_id")
# Clean the new columns that have been joined to 'combined'
combined <- combined %>%
# convert all column names to lower case and remove spaces
clean_names() %>%
# covert new columns to class date
mutate(date_hospitalisation = mdy(date_hospitalisation),
date_outcome = mdy(date_outcome),
date_infection = ymd(date_infection)) %>%
# clean outcome and hospital missings
mutate(outcome = na_if(outcome, ""),
hospital = na_if(hospital, ""))
combined <- combined %>%
mutate(week_onset = yearweek(date_onset, week_start = 1), ## create week of onset variable
week_onset_date = as.Date(week_onset)) ## create a date version
# Patient timelines - pivoting exercise
# Pivoting - patient timelines ----------------------------------------
timelines <- combined %>%
arrange(date_onset) %>% # sort dataset so that earliest are at the top
head(5) %>% # keep only the top 5 rows
select(case_id, starts_with("date")) # keep only certain columns
# Pivot dates longer
timelines_long <- timelines %>%
pivot_longer(
cols = starts_with("date"),
names_to = "date_type",
values_to = "date"
) %>%
mutate(date_type = fct_relevel(date_type, "date_infection", "date_onset", "date_report", "date_hospitalisation", "date_outcome"))
# create plot of patient timelines
ggplot(data = timelines_long, # use the long dataset
mapping = aes(
x = date, # dates of all types displayed along the x-axis
y = case_id, # case_id are discrete, character values
color = date_type, # color of the points
shape = date_type, # shape of the points
group = case_id))+ # this makes the lines appear by color
geom_point(size = 4)+ # show points
geom_line()+ # show lines
theme_minimal()
# USING FCT_LUMP ........
ggplot(data = combined,
mapping = aes(
x = date_onset,
fill = fct_lump_n(district, 3)))+
geom_histogram(binwidth = 7)+
labs(fill = "District")
# Spotlight on district defined at top of script in YAML
surv %>%
filter(district == params$district) %>% # filter data frame to the district in YAML
tabyl(hospital, sex) %>% # begin the cross-tabulation
adorn_totals("both") %>% # add totals for both rows and columns
qflextable()
under_5_count <- surv %>%
filter(age_years < 5) %>%
nrow()
This is a demo situation report on a hypothetical outbreak of Ebola in Sierra Leone from 2014.
This report was produced on CURRENT_DATE
# extract the hospital column, reduce to the unique values, then return the length of that vector (i.e. the number of unique hospital names)
hosp_num <- surv %>% # begin with the surv dataset
pull(hospital) %>% # extract only the hospital column
unique() %>% # use the unique values
na.omit() %>% # remove NA
length() # return only the number
# names of the hospitals
hosp_names <- surv %>%
pull(hospital) %>%
unique() %>%
na.omit() # remove NA from the vector
Key numbers
- In total, there have been 539 cases reported.
- The first case was reported on 2014-06-16.
- The last case was reported on 2014-11-27.
About this report
This report was produced on r Sys.Date().
- In total, there have been
r nrow(surv)cases reported, includingr fmt_count(surv, age_years < 5)children under 5. - The first case was reported on
r min(surv$date_report, na.rm=T). - The last case was reported on
r max(surv$date_report, na.rm=T). - There are
r hosp_numhospitals that have reported cases:r paste(hosp_names, collapse = ", ").
# SUMMARY TABLES ..........................
# Summary tables ---------------------------------------------->>>>>>>>>>
## tables with janitor ----------------------------------------
surv %>%
tabyl(district, hospital)
surv %>%
tabyl(district) %>%
arrange(desc(n)) %>%
adorn_totals() %>%
adorn_pct_formatting()
## tables with dplyr ------------------------------------------
hospital_info <- surv %>%
group_by(hospital) %>% # get statistics for each hospital
summarise(
n_cases = n(), # number of rows (cases)
max_onset = max(date_onset, na.rm = T), # latest onset date
under5 = sum(age_years < 5, na.rm = T), # number of children under 5
vomit_n = sum(vomit == "yes", na.rm=T), # number vomiting
vomit_pct = percent(vomit_n / n_cases), # percent vomiting
max_wt_male = max(wt_kg[sex == "male"], na.rm = T) # max weight among men
)
# GEOM_VIOLIN ............
ggplot(data = surv, mapping = aes(x = district, y = age_years))+
geom_jitter(color = "blue", size = 5, shape = 1)+
geom_violin(alpha = 0.5)
ggplot(data = surv, mapping = aes(x = sex, y = ht_cm))+
geom_violin(alpha = 0.5, size = 2)
#BOXPLOT FILL .....................
ggplot(data = surv, mapping = aes(x = district, y = age_years)) +
geom_jitter() +
geom_boxplot(fill = "blue", alpha = 0.5)
ggplot(data = surv, mapping = aes(x = district, y = age_years, fill = district))+
geom_jitter()+
geom_boxplot(alpha = 0.5)
ggplot(data = surv, mapping = aes(x = district, y = age_years, fill = district, color = district))+
geom_jitter()+
geom_boxplot(alpha = 0)
# Bar plots from count dataset ............
case_counts_district <- surv %>%
count(district) # creating a new dataset
ggplot(data = case_counts_district,
mapping = aes(x = district))+
geom_bar()
# GEOM_SMOOTH .......
ggplot(
data = surv,
mapping = aes(
x = date_report,
y = diff))+
geom_point()+
geom_smooth()
# Histogram breaks .....
ggplot(data = surv, mapping = aes(x = age_years, fill = sex))+
geom_histogram(bins = 100)
ggplot(data = surv, mapping = aes(x = age_years, fill = sex))+
geom_histogram(bins = 5)
ggplot(data = surv, mapping = aes(x = age_years, fill = sex))+
geom_histogram(binwidth = 5)
ggplot(
data = surv,
mapping = aes(x = date_onset))+
geom_histogram()+
labs(
caption = str_glue(
"Data include {num_cases} cases and are current to {current_date}.\n{n_missing_onset} cases are missing date of onset and not shown",
num_cases = nrow(surv),
current_date = format(Sys.Date(), "%d %b %Y"),
n_missing_onset = fmt_count(surv, is.na(date_onset))
)
)
Including Plots
You can also embed plots, for example:
plot(pressure)
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.
# Heat plots
# count cases by hospital-week
hospital_weeks_raw <- combined %>%
count(hospital, week_onset_date)
# create longer dataset of all possible hospital-weeks
expanded <- hospital_weeks_raw %>%
select(hospital, week_onset_date) %>%
tidyr::expand(., week_onset_date, hospital)
# merge so that all hospital-weeks are represented in the data
hospital_weeks <- hospital_weeks_raw %>%
right_join(expanded) %>%
mutate(n = replace_na(n, 0)) %>%
filter(week_onset_date >= ymd("2014-06-15"))
# Adding heatplot
ggplot(data = hospital_weeks,
mapping = aes(x = week_onset_date, y = hospital, fill = n)) +
geom_tile()
# Modifying hospital_weeks to create a new data frame: hospital_weeks_adjusted .................
# Edit the dataset
hospital_weeks_adjusted <- hospital_weeks %>%
mutate(hospital = fct_na_value_to_level(hospital, level = "Missing"),
hospital = fct_relevel(
hospital,
"Central Hospital", "Military Hospital", "Port Hospital", "SMMH", "Other", "Missing"))
# changing the color scheme to a diverging color scale.....................
ggplot(
data = hospital_weeks_adjusted,
mapping = aes(
x = week_onset_date,
y = hospital,
fill = n)) +
geom_tile()+
scale_fill_gradient(low = "skyblue", high = "tomato")
# Changing x asis to show every 2 weeks rathen than month .................
ggplot(
data = hospital_weeks_adjusted,
mapping = aes(
x = week_onset_date,
y = hospital, fill = n)) +
geom_tile()+
scale_fill_gradient(low = "skyblue", high = "tomato")+
scale_x_date(
date_breaks = "2 weeks",
labels = scales::label_date_short()
)
# Making the dataset more presentable ....
# Plot with adjustments
ggplot(data = hospital_weeks_adjusted) +
geom_tile(mapping = aes(
x = week_onset_date,
y = hospital,
fill = n))+
scale_fill_gradient(low = "skyblue", high = "tomato")+
scale_x_date(
date_breaks = "2 weeks",
labels = scales::label_date_short())+
labs(
x = "Week of symptom onset",
y = "Hospital",
fill = "Number of\nweekly cases")+
theme_minimal()