-
Notifications
You must be signed in to change notification settings - Fork 0
/
who-web-scrape.Rmd
286 lines (235 loc) · 12 KB
/
who-web-scrape.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
---
title: "who-web-scrape"
output: github_document
always_allow_html: yes
editor_options:
chunk_output_type: console
---
```{r setup, include=FALSE}
options(width = 100)
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
out.width = "100%"
)
suppressPackageStartupMessages(library(tidyverse))
library(countrycode)
library(emo)
library(gganimate)
library(ggthemes)
library(janitor)
library(leaflet)
library(maps)
library(readxl)
library(rnaturalearth)
library(rvest)
library(here)
```
# Disease incidence and immunization coverage: An exercise in data scraping and visualization `r emo::ji('world_map')` `r emo::ji('woman_health_worker')` `r emo::ji('mag')`
### Motivation
WHO makes a plethora of immunization-related datasets available [online](http://www.who.int/immunization/monitoring_surveillance/data/en/); however, as typically encountered when working with public datasets from international organizations, data are presented as tricky html tables or machine-unfriendly spreadsheets `r emo::ji('robot')` `r emo::ji('boom')`
Here we will use `rvest`, `readxl`, and `tidyverse` to scrape WHO measles data, clean and compile a tidy dataset, and create an interactive map.
### Web scrape
As you see below, WHO presents incidence data by country from 1980 to 2017 for 11 different diseases in html tables or an Excel spreadsheet.
![](images/who_measles_incidence.png)
By inspecting the webpage in the browser ([see here for a brief tutorial](https://blog.rstudio.com/2014/11/24/rvest-easy-web-scraping-with-r/)), we see that the third table element contains the data. `rvest` functions `html_nodes` and `html_table` allow us to scrape the data from that table element and pull them into a dataframe.
```{r}
# parse the html
base_url <- "http://apps.who.int/immunization_monitoring/globalsummary/timeseries/tsincidence"
html <- read_html(str_c(base_url, "measles"))
# grab the third table element and extract its contents
measles_df <- html %>%
html_nodes("table") %>%
.[[3]] %>%
html_table(fill = TRUE, trim = TRUE, header = TRUE) %>%
as_data_frame()
# inspect the output
measles_df
```
#### Clean up the dataframe
Now that we have a dataframe of the reported measles cases, we want to do some cleaning. First, there are quite a lot of empty strings in the table so we want to ensure those are explicit `NA` values. Next, we want to remove the thousands separator which is curiously an apostrophe. Lastly, we want to convert the data type for cases and years from `<chr>` to `<int>`.
```{r}
measles_cases_df <- measles_df %>%
na_if("") %>% # coerce blanks to explicit NAs
rename(country = 1) %>% # rename first column
drop_na(country) %>% # drop any errant rows with NAs for country
gather(year, cases, -country) %>% # reshape df
mutate_at(vars(-country), funs(str_replace_all(., "'", ""))) %>% # remove apostrophe thousands separator
mutate_at(vars(-country), funs(as.integer)) # convert cases to integer
# Let's wrap these steps into a function for later use
clean_web_table <- function(df) {
df %>%
as_data_frame() %>%
na_if("") %>%
rename(country = 1) %>%
drop_na(country) %>%
gather(year, cases, -country) %>%
mutate_at(vars(-country), funs(str_replace_all(., "'", ""))) %>%
mutate_at(vars(-country), funs(as.integer))
}
```
#### Scrape, clean, and compile all disease incidence tables
Let's now adapt our approach so that we can scrape all the available incidence data. To start, we build a vector containing each of the available disease tables. Using `purrr`, we parse each page's html content and table as before and use `map_df` and `clean_web_table` to pull the output into a tidy dataframe.
```{r}
diseases <- c("diphtheria","japenc","measles","mumps","pertussis","rubella","crs","ntetanus","ttetanus","yfever")
diseases %>%
set_names() %>% # set disease names for iteration and .id variable
map_df(
~str_c(base_url, .x) %>% # build urls
read_html() %>% # parse html per page
html_nodes("table") %>% # specify table elements
.[[3]] %>% # grab third table
html_table(fill = TRUE, header = TRUE) %>% # process table
clean_web_table(), # clean up dfs
.id = "disease" # name the id column
)
```
### Process spreadsheet data
Spreadsheets encountered in the wild are often filled with non-tabular data containing row headers, multi-level headers, formatting as data, and/or multiple tables represented across the same sheet
([see this article for spreadsheet best practices](https://amstat.tandfonline.com/doi/full/10.1080/00031305.2017.1375989#.W3mcApNKjeQ)). And, dealing with public datasets from large international organizations is no exception.
Let's grab national immunization coverage estimates to complement the measles incidence data we already scraped. Fortunately, the workbook for these coverage estimates is straightforward to load and clean, as it uses consistent wide data structures and lacks many of the usual spreadsheet quirks `r emo::ji('clinking_glasses')`
Note: I adhere to the [project-oriented workflow](https://www.tidyverse.org/articles/2017/12/workflow-vs-script/) lifestyle, thus the liberal use of `here::here()` from here onwards.
```{r}
# grab WHO immunization coverage estimates
download.file("http://www.who.int/entity/immunization/monitoring_surveillance/data/coverage_estimates_series.xls?ua=1",
here("data", "coverage_estimates_series.xls"))
# peek at the sheet names for measles (MCV1) and note any extraneous sheets (Readme; Regional Global estimates)
here("data", "coverage_estimates_series.xls") %>%
excel_sheets()
read_excel(here("data", "coverage_estimates_series.xls"), sheet = "MCV1")
# let's reshape, clean the table, and create a dataframe for later use
measles_cvg_df <- read_excel(here("data", "coverage_estimates_series.xls"), sheet = "MCV1") %>%
gather(year, coverage, -c(1:4)) %>%
na_if("") %>%
rename(country = Cname) %>%
clean_names() %>% # handy function from janitor() to harmonize variable name format
mutate_at(vars(year, coverage), as.integer) %>%
select(-region)
measles_cvg_df
```
Now that we are familiar with the workbook, let's iterate through each sheet to extract all data from the workbook. We can also pipe all of our reshape and cleaning steps to compile a tidy dataset in one pass.
```{r, warning=FALSE}
excel_path <- here("data", "coverage_estimates_series.xls")
excel_path %>%
excel_sheets() %>%
.[-c(1,length(.))] %>% # drop first and last sheets
set_names() %>%
map_df(
~read_excel(path = excel_path, sheet = .x, trim_ws = TRUE) %>%
gather(year, coverage, -c(1:4)) %>%
na_if("") %>%
rename(country = Cname) %>%
clean_names() %>%
mutate_at(vars(year, coverage), as.integer) %>%
select(-region)
)
```
### Visualize
Ok, now we have reported cases of measles as well as estimated vaccination coverage for measles. Let's use `rnaturalearth` to grab the global mapping data and then employ `leaflet` to build an interactive choropleth map to visualize 2017 reported cases and immunization coverage for measles.
```{r, message=FALSE, warning=FALSE, eval=FALSE}
# prep cases and coverage data (2017) for map
measles_cases_map_dta <- measles_cases_df %>%
mutate(iso_code = countrycode(country, "country.name", "iso3c"), # key for naturalearth join
iso_code = replace(iso_code, country == "Eswatini", "SWZ"), # previously Swaziland
# set bubble size
bubble_radius = findInterval(cases, c(0,25,50,100,500,1000,5000,10000,20000,40000))*2) %>%
filter(year == 2017) %>%
select(-country, -year)
measles_cvg_map_dta <- measles_cvg_df %>%
filter(year == 2017) %>%
select(-year)
# download dataset of latitude/longitude by country for bubble placements
download.file("https://opendata.socrata.com/api/views/mnkm-8ram/rows.csv?accessType=DOWNLOAD",
here("data", "lat_lon.csv"))
lat_lon <- read_csv(here("data", "lat_lon.csv")) %>%
select(iso_a3 = `Alpha-3 code`, lat = `Latitude (average)`, lon = `Longitude (average)`) %>%
add_row(iso_a3 = "SSD", lat = 6.9, lon = 31.3) # South Sudan was missing
# load naturalearth data
map_dta <- ne_countries(scale = 50, returnclass = "sf") %>%
# address some missing iso3 codes
mutate(iso_a3 = replace(iso_a3, sovereignt == "France", "FRA"),
iso_a3 = replace(iso_a3, sovereignt == "Denmark", "DNK"),
iso_a3 = replace(iso_a3, sovereignt == "Norway", "NOR")) %>%
select(iso_a3) %>%
# join reported cases, coverage, and lat/lon data to naturalearth data
left_join(., measles_cases_map_dta, by = c("iso_a3" = "iso_code")) %>%
left_join(., measles_cvg_map_dta, by = c("iso_a3" = "iso_code")) %>%
left_join(., lat_lon, by = "iso_a3")
# set palette and tooltip formats
cvg_pal <- colorBin(palette = "RdYlGn",
domain = map_dta$coverage,
bins = c(0,49,79,89,100),
reverse = FALSE)
cases_lbl <- str_glue("<b>{map_dta$country}</b><br>
Reported cases of measles (2017): {map_dta$cases}")
cvg_lbl <- str_glue("<b>{map_dta$country}</b><br>
MCV1 coverage (2017): {map_dta$coverage}%")
# build map
leaflet(data = map_dta) %>%
setView(lat = 25, lng = 30, zoom = 2.5) %>%
addProviderTiles("CartoDB.Positron", group = "Simple") %>%
# immunization coverage
addPolygons(
fillColor = ~cvg_pal(coverage),
weight = 1,
color = "white",
fillOpacity = 0.8,
popup = cvg_lbl
) %>%
# reported cases
addCircleMarkers(
lng = ~lon, lat = ~lat,
radius = ~bubble_radius,
fillOpacity = 0.7,
color = "purple",
stroke = FALSE,
popup = cases_lbl
) %>%
addLegend(
pal = cvg_pal,
values = ~map_dta$coverage,
title = "Measles (MCV1) coverage, 2017",
position = "topright",
opacity = 0.7,
na.label = "Not available",
labFormat = labelFormat(
suffix = "%",
transform = function(coverage) sort(coverage, decreasing = FALSE)
)
)
```
As the rmarkdown `github_document` output format does not support interactive `leaflet` maps, here is a static image of the map generated from the above code:
![](images/who_measles_leaflet_static.png)
The hmtl file containing the interactive `leaflet` map can be found [here](who-web-scrape.html).
Let's use the excellent `gganimate` package to create an animated map of annually reported measles cases since 2000.
```{r, message=FALSE, warning=FALSE, eval=FALSE}
# set up measles cases dataframe for mapping
lat_lon <- read_csv(here("data", "lat_lon.csv")) %>%
select(iso_a3 = `Alpha-3 code`, lat = `Latitude (average)`, lon = `Longitude (average)`) %>%
add_row(iso_a3 = "SSD", lat = 6.9, lon = 31.3) # South Sudan was missing
map_anim_df <- measles_cases_df %>%
mutate(
iso_code = countrycode(country, "country.name", "iso3c"), # key for lat/lon join
iso_code = replace(iso_code, country == "Eswatini", "SWZ") # previously Swaziland
) %>%
left_join(., lat_lon, by = c("iso_code" = "iso_a3")) %>%
na_if(0) %>% # suppress zeros as we only want cases
filter(year >= 2000)
# build animated map
map_anim <- ggplot(data = map_anim_df) +
borders(database = "world", colour = "gray85", fill = "gray80") +
geom_point(aes(x = lon, y = lat, size = cases), colour = "red", alpha = 0.4) +
scale_size_continuous(labels = scales::comma,
breaks = c(10,1000,10000,50000,100000,200000),
range = c(1, 10)) +
scale_y_continuous(limits = c(-60, 90)) + # exclude Antarctica
coord_quickmap() +
labs(title = "Measles reported cases, {frame_time}",
caption = "Source: WHO vaccine-preventable diseases: monitoring system 2018 global summary",
size = "") +
theme_map(base_size = 16) +
transition_time(year)
animate(map_anim, nframes = 18, fps = 1, detail = 2, length = 36, width = 1000)
anim_save(here("images", "measles_anim.gif"))
```
![](images/measles_anim.gif)