6  Joining data and creating maps

6.1 Joining data

In this section we will combine many of the skills that we have learned so far. We will learn how to join datasets together, and then how to create maps using the joined data. We will also learn how to create a map of Zambia using the sf package, and how to plot the data on top of it.

Let’s start by loading in the packages and datafiles that we will be using in this section. We will use the sf package to create maps, and the ggplot2 package to plot the data on top of the maps. We will also use the dplyr package to join the datasets together.

# Load in the packages we will be using
library(tidyverse)
library(sf)
library(terra)
library(lubridate)

# Reload our scramble function 
scram_num <- function(x, offset = 0.5){
  set.seed(10)
  round(runif(1, x*offset, x*(1+offset)))
}

# Load the malaria cases data 
case_data <- read_csv("data/district-cases-long.csv") %>%
  rowwise() %>%
  mutate(count = scram_num(count))

Question

1. What issues are there with the new population dataset? (Hint: Try making a table of the variables)

Let’s also load a new dataset containing 2020 district population estimates from GRID3. This is an older dataset, but it is still useful for learning how to join datasets together. To make things easy, we have already placed this data in a Google Drive folder, and we can load it directly in R using the URL.

We will use this dataset to create a map of Zambia with the malaria cases data plotted on top of it.

# Load in the population data
dist_pop <- read_csv("https://drive.google.com/uc?export=download&id=1ypf9TfZiZsJqZOLYzQ_1EdcGMXkYgFG4") 

# Check out the population data
names(dist_pop)

[1] "...1"     "provive"  "District" "2020 pop"

table(dist_pop$district, useNA = "always")

<NA> 
   0 

At a glance, we can see that the names of the columns are not conducive to our analysis. We can either reload the dataset with the janitor::clean_names() function piped in, or, as this is a dataset with few columns, we can rename the columns manually using rename or using a shortcut in the select() command. In select, we can use the following formula – select(“new name” = “old name”).

Also, we can see from looking at a table of the district variable, that somebody left a “Test” district in, where the population is NA. Let’s remove that with a filter() command.

dist_pop_clean <- dist_pop %>%
  select(c("province" = "provive",
           "district" = "District",
           "pop_2020" = "2020 pop")) %>%
  filter(district != "Test")

There may still be some issues that we haven’t fixed yet (it is usually a good idea to assume that there are some hidden issues with your data), so let’s be on the lookout for those as we continue.

We will also want to take time at the start of our analysis to prepare our datasets so that they are on a similar level of analysis. In this case, we will be joining our case_data dataset with a district-level 2020 population file. Our case_data is already at the district level, but it isn’t aggregated to Annual statistics, which would make our lives easier later on. Let’s aggregate the data to an annual level now, and filter it so it only retains 2020 case data.

case_data_2020 <- case_data %>% filter(!is.na(district)) %>% 
  mutate(year = lubridate::year(period)) %>%
  group_by(year, district, province, data_type) %>% 
  summarise(count = sum(count, na.rm=TRUE)) %>%
  filter(year == 2020) %>%
  ungroup()

Remember to ungroup your data after using group_by() and summarise()!

This is important because if you don’t, the data will still be grouped when you try to use it in other functions, which can cause errors. You can ungroup your data by using the ungroup() function from the dplyr package.

6.1.1 Joining Datasets Together

When we want to combine data in R, there are many potential techniques we can use. The most common dplyr-based method is to use a command called left_join(). Conceptually, this means that we will be picking one dataset to be our primary dataset (think of holding it in your left hand) and then joining a new dataset to that by one or more matching “key” variables (think of adding something from your right hand into your left hand, fitting it together with a magnet). This left_join ensures that we will only get observations from the secondary or right-side dataset that have a match in our primary or left-side dataset. See the below image for an in-depth examination of left_join() from https://www.guru99.com/r-dplyr-tutorial.html#3

This may not always be what you want to do in the future, you may wish to include all observations from both datasets, or only those that don’t match. For these applications, there are many other types of joins, and I would encourage you to read up on them in this nice tutorial https://www.garrickadenbuie.com/project/tidyexplain/

Because joining can sometimes work less than perfectly, it is VERY IMPORTANT to check your newly joined dataset, to ensure that things went as planned. Before you join, you should have a rough idea of what you expect the new dataset to look like (# of Variables, # of Observations, etc.), and then make sure to confirm those assumptions after the join goes through.

Question

What do you expect the new dataset to look like after joining it with the district populations data?

For now, let’s try to do a left_join, using our case_data as the primary (left-side) dataset, and the dist_pop_clean dataset as our secondary (right-side) dataset.

case_data_pop <- case_data_2020 %>%
  left_join(dist_pop_clean)

As you can see from the message in your R console, if you do not supply a key variable for R to join by, it will automatically try to identify any set of variables that match between your two datasets. This is helpful, but dangerous. R does not always do a perfect job of matching, and it is much safer to specify what columns to join by in your left_join command, using a “by” statement, as shown below.

case_data_pop <- case_data_2020 %>%
  left_join(dist_pop_clean, by = c("district", "province"))

Question

Now that we have joined our two datasets, what issues can you find with the newly joined dataset?

table(case_data_pop$pop_2020, useNA = "always")

30704.60156 31593.93555 33673.07031 37006.19922 37020.23047 40941.14063 
          4           2           5           3           5           3 
45870.40625 48874.85938 50516.83203 52814.09375 54163.03906 54581.34375 
          2           3           5           5           2           5 
54907.13672 56695.67188 57503.09766 58679.96094 58700.12109 59849.41797 
          3           3           3           5           5           5 
61512.64453 61824.52344 62841.11328 67026.73438  67056.1875 67375.28906 
          5           3           3           3           3           5 
69251.94531    71401.75 77282.84375 77521.48438   77898.125 79467.95313 
          5           5           5           5           5           3 
81023.49219 81508.78125   82079.375 83644.24219 83905.03906 84133.35156 
          5           5           5           5           5           3 
85814.40625 87761.92188 87921.94531  91112.5625   91721.875 92058.78125 
          5           3           5           4           4           3 
 92482.1875 96325.78125 98767.44531 101113.6953 103601.3125 103981.0078 
          3           5           3           3           5           3 
106124.9453  106499.375 107106.2422 107900.7969 108828.8672 110095.6328 
          3           5           3           5           3           4 
114048.4219 115984.1563 118126.1328 118133.6484 118237.6563 118544.1641 
          5           3           3           5           5           3 
118883.8047 123900.9219 132333.6406 133419.8281  135459.875 135854.7813 
          5           3           5           5           3           3 
137886.0313 138663.4375 141698.2969  142577.625 146013.9531 148235.6563 
          3           3           3           5           3           5 
150378.2031 150965.9219 152579.0781 156905.9531  158816.875 161032.0313 
          5           5           3           5           5           5 
162097.4375 162234.1406 163914.8125 168056.4531 168097.3594 169051.4688 
          5           5           5           5           5           3 
170751.8438  174675.875 177167.3906 177633.4531 179343.9531 181282.9688 
          3           5           5           4           3           5 
183040.1875 189030.5469 190355.9219   211248.75 224233.1563  226902.625 
          5           3           3           3           5           5 
236634.4219 237247.6094 238477.2188 241443.6406 241907.3906 244031.2344 
          5           5           3           5           5           3 
247608.5313  269249.125 271302.5938      292503 302964.5313 337449.4688 
          5           3           3           3           5           3 
     338920 375986.1563 452313.1563 458736.2188    513965.5  2389815.75 
          5           3           3           3           3           5 
       <NA> 
         10 

We only added 1 extra variable and no new observations, which is what we expected! However, it seems there are many NAs in this dataset that weren’t in the dist_pops when we loaded it in. This usually happens with joins where two columns that we think are the same (our “keys”) are not actually exactly the same. There may be one value that appears in our right-key column but not in our left-key column, so now NAs have been introduced because there was nothing to match it with in the left-key column. There are lots of ways to find out exactly where we went wrong, but usually with province and district, it is the case that some spelling error was introduced in one of the datasets, or an alternative spelling was used.

Question

Where are the errors coming from in the joined dataset? (Hint: There are at least two discrepancies, one each in our province and district fields respectively)

We can fix these issues using the ‘ifelse’ function. Let’s go back into the dist_pops dataset, fix those issues, and try the join again.

#First fix the issues in the district populations dataset
dist_pop_clean <- dist_pop %>%
  select(c("province" = "provive",
           "district" = "District",
           "pop_2020" = "2020 pop")) %>%
  filter(district != "Test") %>%
  mutate(province = ifelse(district == "Siavonga", "Southern", province),
        district = ifelse(district == "Kalaba", "Kalabo", district))

#Now Join the two datasets
case_data_pop <- case_data_2020 %>%
  left_join(dist_pop_clean, by = c("district", "province"))

Now that we have redone our join, did we fix the problem? Make sure to check your variables to see if things went as expected, or if there are any other issues that we may have missed.

table(case_data_pop$pop_2020, useNA = "always")

30704.60156 31593.93555 33673.07031 37006.19922 37020.23047 40941.14063 
          4           2           5           3           5           3 
45870.40625 48874.85938 50516.83203 52814.09375 54163.03906 54581.34375 
          2           3           5           5           2           5 
54907.13672 56421.33203 56695.67188 57503.09766 58679.96094 58700.12109 
          3           5           3           3           5           5 
59849.41797 61512.64453 61824.52344 62841.11328 67026.73438  67056.1875 
          5           5           3           3           3           3 
67375.28906 69251.94531    71401.75 77282.84375 77521.48438   77898.125 
          5           5           5           5           5           5 
79467.95313 81023.49219 81508.78125   82079.375 83644.24219 83905.03906 
          3           5           5           5           5           5 
84133.35156 85814.40625 87761.92188 87921.94531  91112.5625   91721.875 
          3           5           3           5           4           4 
92058.78125  92482.1875 96325.78125 98767.44531 101113.6953 103293.1406 
          3           3           5           3           3           5 
103601.3125 103981.0078 106124.9453  106499.375 107106.2422 107900.7969 
          5           3           3           5           3           5 
108828.8672 110095.6328 114048.4219 115984.1563 118126.1328 118133.6484 
          3           4           5           3           3           5 
118237.6563 118544.1641 118883.8047 123900.9219 132333.6406 133419.8281 
          5           3           5           3           5           5 
 135459.875 135854.7813 137886.0313 138663.4375 141698.2969  142577.625 
          3           3           3           3           3           5 
146013.9531 148235.6563 150378.2031 150965.9219 152579.0781 156905.9531 
          3           5           5           5           3           5 
 158816.875 161032.0313 162097.4375 162234.1406 163914.8125 168056.4531 
          5           5           5           5           5           5 
168097.3594 169051.4688 170751.8438  174675.875 177167.3906 177633.4531 
          5           3           3           5           5           4 
179343.9531 181282.9688 183040.1875 189030.5469 190355.9219   211248.75 
          3           5           5           3           3           3 
224233.1563  226902.625 236634.4219 237247.6094 238477.2188 241443.6406 
          5           5           5           5           3           5 
241907.3906 244031.2344 247608.5313  269249.125 271302.5938      292503 
          5           3           5           3           3           3 
302964.5313 337449.4688      338920 375986.1563 452313.1563 458736.2188 
          5           3           5           3           3           3 
   513965.5  2389815.75        <NA> 
          3           5           0 

Challenge Question

What are the 10 highest malaria-burden districts by confirmed cases per 1000 population for 2020 (Hint: Use mutate() to make a new column that = cases divided by population times 1000, and make sure to use only the correct data elements)?

How is this list different than if we were only using the raw case numbers? Bonus points to anyone who shows this information using a plot

6.1.2 Calculating incidence rates

Now that we have joined our two datasets, we can start to use them for analysis! In this section, we will learn how to calculate incidence rates using the joined data by apply some of the concepts from our previous sessions.

First, let’s calculate the total cases in each district for 2020. We can do this by using the group_by() and summarise() functions from the dplyr package. We will also calculate the incidence rate by dividing the total cases by the population and multiplying by 1000.

# Calculate the total cases and incidence rate for each district
case_incidence_2020 <- case_data_pop %>%
  filter(data_type %in% c("Confirmed", "Clinical", "Confirmed_Passive_CHW")) %>%
  group_by(district, province) %>%
  summarise(total_cases = sum(count, na.rm = TRUE),
            pop_2020 = mean(pop_2020, na.rm = TRUE),
            cases_per_1000 = (total_cases / pop_2020) * 1000) %>%
  ungroup()

Let’s plot out the incidence rates for each district in 2020. We can use the ggplot2 package to create a bar plot of the incidence rates. We will also use the coord_flip() function to flip the axes so that the districts are on the y-axis and the incidence rates are on the x-axis. We will also use the facet_wrap() function to create a separate plot for each province.

# Plot the incidence rates for each district in 2020
ggplot(case_incidence_2020, aes(x = reorder(district, -cases_per_1000), y = cases_per_1000)) +
  geom_col() +
  coord_flip() +
  facet_wrap(~province, scales = "free_y") +
  labs(title = "Malaria Incidence Rates by District in 2020",
       x = "District",
       y = "Incidence Rate (cases per 1000 population)") +
  theme_minimal()

We can also calculate the incidence rates for each province by using the group_by() and summarise() functions again. All we need to do is remove “district” from the group_by() step.

# Calculate the total cases and incidence rate for each province
province_incidence_2020 <- case_data_pop %>%
  filter(data_type %in% c("Confirmed", "Clinical", "Confirmed_Passive_CHW")) %>%
  # Get District toatls
  group_by(district, province) %>%
  summarise(total_cases = sum(count, na.rm = TRUE),
            pop_2020 = mean(pop_2020, na.rm = TRUE)) %>%
  # Get Province totals
  group_by(province) %>%
  summarise(total_cases = sum(total_cases, na.rm = TRUE),
            pop_2020 = sum(pop_2020, na.rm = TRUE),
            cases_per_1000 = (total_cases / pop_2020) * 1000) %>%
  # Remove grouping
  ungroup()

# Plot the incidence rates for each province in 2020
ggplot(province_incidence_2020, aes(x = reorder(province, -cases_per_1000), y = cases_per_1000)) +
  geom_col() +
  coord_flip() +
  labs(title = "Malaria Incidence Rates by Province in 2020",
       x = "Province",
       y = "Incidence Rate (cases per 1000 population)") +
  theme_minimal()

6.2 Creating maps

Now that we have joined and cleaned our two datasets, we can start to use them for analysis! In this section, we will learn how to create maps using the sf package and how to plot the data on top of the maps. We will also learn how to create a map of Zambia using the sf package, and how to plot the data on top of it.

Let’s start by loading in the province-level shapefile for Zambia. This should be saved in your data folder as an RDS (R data file) file. We will use the readRDS() function to load it in.

# Load in the shapefile for Zambia
province_shp <- readRDS("data/province_shp.rds")

# Check out the shapefile
province_shp

Simple feature collection with 10 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 21.99927 ymin: -18.07742 xmax: 33.70903 ymax: -8.271976
Geodetic CRS:  WGS 84
   geo_province census_pop_22 pop_grid3_2022  area_km2
1       Central       2252483      2252489.2  94879.42
2    Copperbelt       2757539      2757494.0  31224.70
3       Eastern       2454788      2454669.8  69160.09
4       Luapula       1514011      1513987.1  49948.03
5        Lusaka       3079964      3079970.8  22165.02
6      Southern       2381728      2381709.0  85627.94
7      Muchinga        918296       918241.1  69376.69
8      Northern       1618412      1618431.2  77893.17
9  Northwestern       1270028      1270027.8 125632.63
10      Western       1363520      1363523.4 127807.26
                         geometry
1  POLYGON ((28.34492 -15.1360...
2  POLYGON ((27.51905 -12.2943...
3  POLYGON ((32.6429 -14.20324...
4  POLYGON ((28.76563 -11.9935...
5  POLYGON ((27.96051 -15.5838...
6  POLYGON ((28.9341 -15.94181...
7  POLYGON ((32.56641 -10.4470...
8  POLYGON ((30.61348 -11.3285...
9  POLYGON ((22.003 -13.45797,...
10 POLYGON ((22.00166 -14.4503...

class(province_shp)

[1] "sf"         "data.frame"

names(province_shp)

[1] "geo_province"   "census_pop_22"  "pop_grid3_2022" "area_km2"      
[5] "geometry"      

We can see that this is a simple shapefile with only one geometry column, which is what we want. We can also see that the shapefile has a column called “geo_province” which contains the names of the provinces. We will use this column to join our data to the shapefile.

6.2.1 Joining the data to the shapefile

To join the data to the shapefile, we will use the left_join() function from the dplyr package. We will join the province_incidence_2020 dataset to the province_shp shapefile using the “geo_province” column in the shapefile and the “province” column in the dataset.

Note that in this case we want to join our data sets using columns that have different names. In the shapefile, the column is called “geo_province” and in the dataset it is called “province”. We can use the by argument in the left_join() function to specify which columns to join by. We will use the c() function to create a vector of the column names to join by.

# Join the data to the shapefile
province_shp_inc_2020 <- province_shp %>%
  left_join(province_incidence_2020, by = c("geo_province" = "province"))

6.2.2 Creating the incidence rate map

Now that we have joined the data to the shapefile, we can create a map of Zambia with the incidence rates plotted on top of it. We will use the ggplot2 package to create the map. We will use the geom_sf() function to plot the shapefile and the geom_sf() function again to plot the data on top of it. We will also use the scale_fill_viridis_c() function to create a color scale for the incidence rates. We will use the theme_minimal() function to create a clean theme for the map.

# Create the map of Zambia with the incidence rates plotted on top of it
ggplot(province_shp_inc_2020) +
  geom_sf(aes(fill = cases_per_1000), color = "white") +
  scale_fill_viridis_c(option = "plasma", name = "Incidence Rate (cases per 1000 population)") +
  labs(title = "Malaria Incidence Rates by Province in 2020",
       subtitle = "Zambia") +
  theme_minimal() +
  theme(legend.position = "bottom")

Alternativelty, we may want to create our own color palette, based on set incidence values. For instance, we can use the “Zone” classification from the Zambia National Malaria Elimination Centre (NMEC) to create a color palette for the map. We can use the scale_fill_manual() function to create a custom color scale for the map. We will use the c() function to create a vector of the colors to use for each zone.

# use mutate() to create a new column for the zone classification
province_shp_inc_2020 <- province_shp_inc_2020 %>%
  mutate(zone = case_when(
    cases_per_1000 < 1 ~ "0", 
    cases_per_1000 >= 1 & cases_per_1000 <50 ~ "1 - 49", 
    cases_per_1000 >= 50 & cases_per_1000 < 200 ~ "50 - 199",
    cases_per_1000 >= 200 & cases_per_1000 < 500 ~ "200 - 499", 
    cases_per_1000 >=500 ~ "500 and above"), 
    # Convert to factor
    zone = factor(zone, 
    levels = c("0", "1 - 49", "50 - 199", "200 - 499", "500 and above")))

# Incidence map with custom color palette
ggplot(province_shp_inc_2020) +
  geom_sf(aes(fill = zone), color = "white") +
  scale_fill_manual(values = c("#f2f2f2","#cccccc","#fecc5c","#f74436","#ab0e02")) +
  labs(title = "Malaria Incidence Rates by Province in 2020",
       subtitle = "Zambia") +
  theme_minimal() +
  theme(legend.position = "bottom")

6.2.3 Challenge Questions

1. Create a district-level map of Zambia with the incidence rates plotted on top of it using the sf package and the ggplot2 package. Use the custom color palette that you created in the previous section.
2. Try using the filter() function to generate a map of
   1. Only districts in a specific province (e.g. Lusaka)
   2. Only districts with a specific incidence rate (e.g. > 50 cases per 1000 population)
3. Try using the facet_wrap() function to create a separate map for each province.
4. Try using the geom_sf_label() function to add labels to the map for each district. You can use the label argument to specify which column to use for the labels. For example, you can use the district column to label each district on the map.

6.2.4 Bonus section: Interactive maps

The ‘tmap’ package is a great way to create interactive maps in R. It is similar to the ‘ggplot2’ package, but it is specifically designed for creating maps. You can use the tm_shape() function to specify the shapefile to use for the map, and the tm_fill() function to specify the fill color for the map. You can also use the tm_borders() function to add borders to the map.

# Load in the tmap package (you may need to install it first)
library(tmap)

# Select the "view" mode for interactive maps
tmap_mode("view")

# Create the interactive map of Zambia with the incidence rates plotted on top of it
tm_shape(province_shp_inc_2020) +
  tm_fill("cases_per_1000", title = "Incidence Rate (cases per 1000 population)", palette = "plasma") +
  tm_borders() +
  tm_layout(title = "Malaria Incidence Rates by Province in 2020",
            subtitle = "Zambia")

The tmap package also allows you to create static maps, which are similar to the maps created with the ggplot2 package. You can use the tmap_mode("plot") function to switch to static mode. You can also use the tm_shape() function to specify the shapefile to use for the map, and the tm_fill() function to specify the fill color for the map. You can also use the tm_borders() function to add borders to the map.

# Switch to static mode
tmap_mode("plot")

# Create the static map of Zambia with the incidence rates plotted on top of it
tm_shape(province_shp_inc_2020) +
  tm_fill("zone", title = "Incidence Rate (cases per 1000 population)", palette = "plasma") +
  tm_borders() +
  tm_layout(title = "Malaria Incidence Rates by Province in 2020",
            subtitle = "Zambia")

There are many other options for customizing the maps created with the tmap package. You can use the tm_layout() function to customize the layout of the map, and the tm_legend() function to customize the legend. You can also use the tm_text() function to add text labels to the map. For more information on how to use the tmap package, you can check out the tmap documentation.