Exercice

This exercice aims at summarizing what was shown througout the training. It includes dealing with shapefile and epidemiological data in “almost real” conditions. The data represents the number of cases of an imaginary disease across Cambodia called R infections from 2018 to 2022. This infection spread rapidly and started during a training in Phnom Penh. Symptoms are very specific and includes installing Rstudio, loading data and using R software for spatial analysis and mapping.

In R, it exist many differents implementation solution that lead to the same results. The solution presented here just provides one implementation among thousands of possibilities.

R REMINDERS

  1. Comment your code ! (# important informations on the code)

  2. Check your R objects ! (plot(), print(), View() , …)

  3. Listen to R outputs ! (Errors AND Warnings)

  4. Get help ! (?name_of_function, internet, other users)

  5. Keep calm and take a break !

Create your working environment

  1. Create a R project called “RGeotraining”

  2. Download and unzip the training data into a directory called “data/”.

Download example data

  1. Create a directory called “img/” for image outputs

  2. Download and unzip shapefile of Cambodian provinces in the ‘data/’ directory: Province shapefile. A new province have been recently delimited in Cambodia (the 26th one). for this exercice you can use the old shapefile that contains 25 provinces in the data directory.

  3. Check your working directory using getwd()

Load and visualize the data

Steps

  1. Load R libraries sf and mapsf,

  2. Load province shapefile with st_read() and set the projection st_transform(),

  3. Load population data with read.table() and sum the counts over provinces using aggregate(),

  4. Load the number of cases per province from the csv file.

Solution

#============================================
# 2. Load and visualize the data
#============================================

### 2.1 Load R libraries
#-----------------------

# library(dplyr) # deals with dataframes
library(sf) # spatial objects
library(mapsf) # Plot 

### 2.2 Spatial data
#----------------------

# Load data
province_sf <- st_read("data/khm_admbnda_adm1_gov_20181004.shp", quiet = TRUE)
head(province_sf) # We need to define the projection

# set crs
province_sf <- st_transform(province_sf, crs = 32648)

# Visualize
class(province_sf) # type of R object
dim(province_sf) # dimensions of the object = n columns + 1 geometry
summary(province_sf) # Summarize the information of the object
plot(province_sf[,1]) # Plot the first column to look at geometry 
# (if you plot all of column it is sometimes too long)

### 2.3 Population data
#--------------------------
# Load 
population_df <- read.table(file = "data/Population_district_Cambodia.csv",
                            sep = ',', header = TRUE)

# Visualize
head(population_df)
class(population_df)
dim(population_df)
colnames(population_df)
summary(population_df)

# Because we work on province we need to aggregate the values
# Tips: ?aggregate
pop_by_district_df <- aggregate(population_df$T_POP, 
                                by = list(ADM1_PCODE = population_df$ADM1_PCODE),
                                FUN = sum)

colnames(pop_by_district_df) <- c("ADM1_PCODE", "pop") # rename column for later


### 2.4 Cases
#-------------------------
# Load number of cases
cases_df <- read.table(file = "data/R_infection_monthly_cases.csv", 
                       header = TRUE, sep = ";")

head(cases_df)
dim(cases_df)
colnames(cases_df)
summary(cases_df)
# What are these data ? 
# What represents each row ? each column ? 
# What is the time span ?

Merge data with shapefiles

Steps

  1. Merge population with sf polygon using merge(),

  2. Merge the new sf object with number of cases,

  3. Identify the merging issues and correct the datasets, you can compares columns of the datasets using boolean (TRUE/FALSE) operation as == (is equal to), %in% (appear in at least once), ! (negate, invert TRUE and FALSE values) and extract row of interest using [row_selection, column_selection].

Solution

#=======================================================
# 3. Merge all data
#=======================================================

### 3.1 Merge population with sf polygons to prepare for mapping
#---------------------------------------------------------------

# What is the key column to merge the objects (column with IDs) ?
province_pop_sf <- merge(pop_by_district_df, province_sf, by = "ADM1_PCODE")
class(province_pop_sf) # What a mess !!! We lost the geometry !

province_pop_sf <- merge(province_sf, pop_by_district_df, by = "ADM1_PCODE")
class(province_pop_sf) # Much better !!!

head(province_pop_sf)

# Create a maps to represent population
mf_export(province_pop_sf, filename = 'img/Population_province.png', width = 500, res = 100)
mf_map(province_pop_sf)
mf_map(province_pop_sf,
       var = "pop" , 
       inches = .2,
       type = "prop",
       col = "#000066",
       leg_title = "Population")
mf_layout(title = "Population in cambodian provinces")
dev.off()

### 3.2 Merge polygon with number of cases
#---------------------------------------------------------------

# What is (are) the key column(s) ?
province_pop_cases_sf <- merge(province_pop_sf, cases_df, by.x = 'ADM1_EN', by.y = "Province")
class(province_pop_cases_sf) # Perfect !
dim(province_pop_cases_sf) # Why do I have only 18 province instead of 25 ? 7 missing values 

# detect names from cases_df that do not match names in province_pop_sf
!(cases_df$Province %in% province_pop_sf$ADM1_EN) # We have 6 provinces that do not match the names (FALSE value)!
sum(!(cases_df$Province %in% province_pop_sf$ADM1_EN))
# which ones ?
cases_df$Province[!(cases_df$Province %in% province_pop_sf$ADM1_EN)]

# What are the names of these provinces in the sf file ?
province_pop_sf$ADM1_EN[!(province_pop_sf$ADM1_EN %in% cases_df$Province)]
# Some names do not have the same spelling ...

### 3.4 Correct the province names
#----------------------------------------------------------

## OPTION 1 : Open the csv file and change the names one by one 

## OPTION 2 (advanced R): Change the names using R 

# What are the names in the cases files ?
dput(cases_df$Province[!(cases_df$Province %in% province_pop_sf$ADM1_EN)])
# correct the names
cases_df$Province[!(cases_df$Province %in% province_pop_sf$ADM1_EN)] <- c("Ratanak Kiri", 
                                                                          "Banteay Meanchey",
                                                                          "Siemreap", 
                                                                          "Mondul Kiri", 
                                                                          "Preah Sihanouk",
                                                                          "Tboung Khmum")
# Be careful of the order of the provinces names !!!

# Merge again
province_pop_cases_sf <- merge(province_pop_sf, cases_df, by.x = 'ADM1_EN', by.y = "Province")
class(province_pop_cases_sf) # Perfect !
dim(province_pop_cases_sf) # One province is still missing !  

# Which one ? 
province_pop_sf$ADM1_EN[!(province_pop_sf$ADM1_EN %in% cases_df$Province)]
# This province does not exists in our cases data. 
# We need to keep it as a NA value in the merged file 

# Merge again
province_pop_cases_sf <- merge(province_pop_sf, cases_df, 
                               by.x = 'ADM1_EN', by.y = "Province", 
                               all.x = TRUE ) 
# We can use all.x = TRUE to keep all rows from x object event if it is not in the y object
class(province_pop_cases_sf) # Perfect !
dim(province_pop_cases_sf) #  YEAH !!!! Wonderful !!!

# what happen to Pailin ? Let extract the row to look at it :
province_pop_cases_sf[province_pop_cases_sf$ADM1_EN == "Pailin",] 
# NA values have been set for the missing data

Map the count data

Steps

  1. Map the number of cases from 2018 to 2022 (mf_map()) in the same figure (you can split your plotting windows with par(mfrow=c(number_of_line, number_of_columns))) and save it (png()),

  2. Add NA values on the map with mf_map(type = "symb").

Solution

#=======================================================
# 4. Map the count data per year
#=======================================================

### 4.1 Create a maps to represent the number of cases
#---------------------------------------------------

# Save the map in img/ directory

png(filename = "img/Count_R_infections.png", 
          width = 700, res = 100) # Run the code until "dev.off()" function

par(mfrow = c(2,2)) # create subplots

# Count from 2018
mf_map(province_pop_cases_sf)
mf_map(province_pop_cases_sf,
       var = "X2018" , 
       inches = .2,
       type = "prop",
       col = "#000066",
       leg_title = "Cases")
mf_layout(title = "Number of cases per province (2018)")

# Count from 2019
mf_map(province_pop_cases_sf)
mf_map(province_pop_cases_sf,
       var = "X2019" , 
       inches = .2,
       type = "prop",
       col = "#000066",
       leg_title = "Cases")
mf_layout(title = "Number of cases per province (2019)")

# Count from 2020
mf_map(province_pop_cases_sf)
mf_map(province_pop_cases_sf,
       var = "X2020" , 
       inches = .2,
       type = "prop",
       col = "#000066",
       leg_title = "Cases")
mf_layout(title = "Number of cases per province (2020)")

# Count from 2021
mf_map(province_pop_cases_sf)
mf_map(province_pop_cases_sf,
       var = "X2021" , 
       inches = .2,
       type = "prop",
       col = "#000066",
       leg_title = "Cases")
mf_layout(title = "Number of cases per province (2021)")

dev.off()

# What about Pailin ? It is treated as a zero instead of NA value and this is a big issue ! 
# How can we add it as a NA values ? 

# 4.2 Deals with NA values
#-----------------------------------------------

# There is no easy solution design for it yet ... 
# What do you think ? We can try to do it by hand. Here is just a suggestions :  

# Baseline map
mf_map(province_pop_cases_sf)

# Add symbol for Pailin
pailin_sf <- province_pop_cases_sf[province_pop_cases_sf$ADM1_EN == "Pailin",] # select row
pailin_sf$data_avail <- "No data" # Define no data value for pailin
mf_map(pailin_sf ,
       var = "data_avail",
       leg_title = NULL, 
       col = "black",
       cex = 1.5,
       pch = 22,
       type = "symb")

# Count from 2021
mf_map(province_pop_cases_sf,
       var = "X2021" , 
       inches = .2,
       type = "prop",
       col = "#000066",
       leg_title = "Cases")
mf_layout(title = "Number of cases per province (2021)")

# How can I improve this map ? change colors ? fix the scale for all subplot ?

# These count are not really informative since it depends on the population.
# We can compute incidence instead of cases in a new column

Transform data to incidence

Steps

  1. Compute incidence for each column (number of cases/ population * 100,000)

  2. Compute incidence using apply() and by creating a function(){},

  3. Merge incidences with shapefile,

  4. Map incidence for each year in the same figure (par(mfrow)) using a loop for(variable in vector){}, you can call help with ?for (variable in vector) {}

Solution

#=======================================================
### 5. Transform to incidence
#=======================================================

### 5.1 Compute incidence
#---------------------------------------------------------

# OPTION 1 : compute columns by columns
province_pop_cases_sf$incidence_2017 <- province_pop_cases_sf$X2017/province_pop_cases_sf$pop * 100000
province_pop_cases_sf$incidence_2018 <- province_pop_cases_sf$X2018/province_pop_cases_sf$pop * 100000
province_pop_cases_sf$incidence_2019 <- province_pop_cases_sf$X2019/province_pop_cases_sf$pop * 100000
province_pop_cases_sf$incidence_2018 <- province_pop_cases_sf$X2020/province_pop_cases_sf$pop * 100000

# OPTION 2 : use "apply" function (?apply)
# remove geometry to work on dataframe
attribute_df <- st_drop_geometry(province_pop_cases_sf)
# set rownames
row.names(attribute_df) <- attribute_df$ADM1_PCODE # easier to use the ID
# remove useless column
colnames(attribute_df) # show the column names of our dataframe
attribute_cases_df <-  attribute_df[, 15:20] # select columns 15 to 20 (contains cases values)

# Create your own function to compute incidence
compute_incidence <- function(cases, population){ 
  # cases and population are parameters of the function
  # case is a numerical values of a number of cases (can be a single value or a vector)
  # population is a numerical values of population count (can be a single value or a vector)
  # Both parameters must have the same length
  incidence <- cases/population * 100000
  return(incidence)
}

class(compute_incidence)
print(compute_incidence)

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# example:
incidence_by_hand <- province_pop_cases_sf$X2017/province_pop_cases_sf$pop * 100000
incidence_with_function <- compute_incidence(cases = province_pop_cases_sf$X2017, population = province_pop_cases_sf$pop)
incidence_by_hand == incidence_with_function # It gives the same results 
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Apply the same function to each column of the dataframe
attribute_inc_mat <- apply(attribute_cases_df, 2, compute_incidence, population = attribute_df$pop)
class(attribute_inc_mat) # apply return a matrix (table of numerical values only)

# turn it to dataframe to add a column of character
attribute_inc_df <- as.data.frame(attribute_inc_mat)
class(attribute_inc_df)

# Change the columns names to inform on data
colnames(attribute_inc_df) <- paste0('incidence_', colnames(attribute_inc_df))

# Retrieve ID code from row names
attribute_inc_df$ADM1_PCODE <- row.names(attribute_inc_df)

# Have look at the final object
head(attribute_inc_df)
summary(attribute_inc_df)
dim(attribute_inc_df)

# Might look longer to use this option but reminds that when you have larger dataframe 
# in hand this option is way more convenient than dealing with each line with the 
# risque of producing mistakes

### 5.2 Merge with the spatial object
#---------------------------------------------------------

province_pop_inc_sf <- merge(x = province_pop_sf, y = attribute_inc_df, by = "ADM1_PCODE")

class(province_pop_inc_sf)
head(province_pop_inc_sf)
summary(province_pop_inc_sf)
dim(province_pop_inc_sf)

### 5.3 Map incidence
#---------------------------------------------------------

# OPTION 1: Just like we did earlier with the number of cases, plot by plot

# OPTION 2: Use a loop !

# A loop repeat a part of your code for many values
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Simple examples:
for(i in 1:10){
  # i is a variable that will successively takes the values contains in the vector given after 'in' 
  print(i)
}

j <- "Fixed value ouside the loop"
for(i in c("a", "b", "d", "end")){
  print(i)
  print(j)
}
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

png(filename = "img/Incidence_R_infections.png", 
    width = 700, res = 100) # Run the code until "dev.off()" function

par(mfrow = c(2,2)) # create subplots
for(year in c("2018", "2019", "2020", "2021")) {
  print(paste0("Plot incidence for the year ", year))
  
  name_col <- paste0("incidence_X", year) # Define names of the plotted column
  mf_map(province_pop_inc_sf)
  mf_map(province_pop_inc_sf,
         var = c( "pop", name_col) , 
         inches = .2,
         type = "prop_choro",
         col = "#000066",
         leg_title = "Incidence")
  mf_layout(title = paste0(" Incidence of R infections per province (", year, ")"))
  
}
dev.off() # Don't forget to close the plotting window

Go further in analysis (Advanced R) …

The section gives suggestions to go further into the data description.

Work with averaged incidence

We are now interested in the averaged incidence between 2017 and 2022. In other terms, we want to compute the mean values of incidence for each row.

Let look at hospital distribution

Do I have higher incidence if there is more hospital in the province ? In other term, is there a bias in case detected cause by the access to health care ?