Global Measles Surveillance: Quality vs Disease Burden

Analysis of 151 countries reveals surveillance systems under pressure

TidyTuesday
Data Visualization
R Programming
2025
Multi-dimensional scaling analysis of WHO measles surveillance data reveals global patterns in disease burden and healthcare monitoring capabilities. Using data from 151 countries, this visualization identifies ‘crisis zones’ where high measles incidence coincides with poor laboratory confirmation rates, highlighting critical gaps in public health infrastructure and surveillance systems worldwide.
Author

Steven Ponce

Published

June 24, 2025

Figure 1: Scatter plot showing measles surveillance quality versus disease burden across 151 countries. Countries are positioned by disease burden (x-axis) and outbreak variability (y-axis). Points are colored by lab confirmation rates and sized by average measles incidence. A ‘Crisis Zone’ in the upper left shows countries with high measles burden and poor surveillance, while a ‘Success Zone’ in the lower right shows countries with low measles burden and good surveillance. Most countries cluster in the middle, with a clear pattern indicating that better measles surveillance correlates with a lower disease burden.

Steps to Create this Graphic

1. Load Packages & Setup

Show code 
```{r}
#| label: load
#| warning: false
#| message: false
#| results: "hide"

## 1. LOAD PACKAGES & SETUP ----
suppressPackageStartupMessages({
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  tidyverse,  # Easily Install and Load the 'Tidyverse'
  ggtext,     # Improved Text Rendering Support for 'ggplot2'
  showtext,   # Using Fonts More Easily in R Graphs
  janitor,    # Simple Tools for Examining and Cleaning Dirty Data
  scales,     # Scale Functions for Visualization
  glue,       # Interpreted String Literals
  patchwork,  # The Composer of Plots
  scico,      # Colour Palettes Based on the Scientific Colour-Maps 
  ggrepel     # Automatically Position Non-Overlapping Text Labels with 'ggplot2'
  )})

### |- figure size ----
camcorder::gg_record(
  dir    = here::here("temp_plots"),
  device = "png",
   width  =  10,
  height =  8,
  units  = "in",
  dpi    = 320
)

# Source utility functions
suppressMessages(source(here::here("R/utils/fonts.R")))
source(here::here("R/utils/social_icons.R"))
source(here::here("R/utils/image_utils.R"))
source(here::here("R/themes/base_theme.R"))
```

2. Read in the Data

Show code 
```{r}
#| label: read
#| include: true
#| eval: true
#| warning: false

tt <- tidytuesdayR::tt_load(2025, week = 25)

cases_year <- tt$cases_year |> clean_names()

tidytuesdayR::readme(tt)
rm(tt)

#' Note: The number of cases of measles and rubella officially reported by a
#' WHO Member State is only available by July of each year. If any numbers from
#' this provisional data are quoted, they should be properly sourced with a date
#' (i.e. "provisional data based on monthly data reported to WHO (Geneva) as
#' of June 2025"
```

3. Examine the Data

Show code 
```{r}
#| label: examine
#| include: true
#| eval: true
#| results: 'hide'
#| warning: false

glimpse(cases_year)
```

4. Tidy Data

Show code 
```{r}
#| label: tidy
#| warning: false

# Data Preparation

country_profiles <- cases_year |>
  group_by(country) |>
  filter(
    n() >= 5, # At least 5 years of data
    sum(measles_total, na.rm = TRUE) >= 10 # Meaningful case numbers
  ) |>
  summarise(
    # Core metrics for MDS
    avg_incidence = mean(measles_incidence_rate_per_1000000_total_population, na.rm = TRUE),
    max_incidence = max(measles_incidence_rate_per_1000000_total_population, na.rm = TRUE),
    cv_incidence = sd(measles_incidence_rate_per_1000000_total_population, na.rm = TRUE) /
      (mean(measles_incidence_rate_per_1000000_total_population, na.rm = TRUE) + 0.001),
    lab_confirmation_rate = mean(measles_lab_confirmed / (measles_total + 0.001), na.rm = TRUE),
    years_with_cases = n(),
    .groups = "drop"
  ) |>
  mutate(
    cv_incidence = pmin(cv_incidence, 5), # Cap CV at reasonable level
    lab_confirmation_rate = pmin(lab_confirmation_rate, 1) # Cap at 100%
  ) |>
  filter(complete.cases(across(avg_incidence:years_with_cases)))


# MDS Analysis

# Prepare MDS matrix with log transformation and scaling
mds_matrix <- country_profiles |>
  transmute(
    log_avg_incidence = log10(avg_incidence + 1),
    log_max_incidence = log10(max_incidence + 1),
    log_cv_incidence = log10(cv_incidence + 1),
    lab_confirmation_rate,
    years_with_cases
  ) |>
  scale() |>
  as.matrix()

rownames(mds_matrix) <- country_profiles$country

# Perform MDS and calculate quality metrics
mds_result <- cmdscale(dist(mds_matrix), k = 2, eig = TRUE)
variance_explained <- sum(mds_result$eig[1:2]) / sum(abs(mds_result$eig))
stress <- sum((dist(mds_matrix) - dist(mds_result$points))^2) / sum(dist(mds_matrix)^2)


# Plot Data
plot_data <- tibble(
  country = country_profiles$country,
  MDS1 = mds_result$points[, 1],
  MDS2 = mds_result$points[, 2],
  avg_incidence = country_profiles$avg_incidence,
  lab_confirmation_rate = country_profiles$lab_confirmation_rate
)

# Outliers Data
outlier_datasets <- list(
  top_left = plot_data |> filter(MDS1 <= -3 & MDS2 > 1),
  bottom_right = plot_data |> filter(MDS1 > 2.8 & MDS2 > -1)
)

# Shared Plot Styling & Annotations
base_label_style <- list(
  size = 3.2,
  fontface = "bold",
  color = "black",
  bg.color = "white",
  bg.r = 0.12,
  box.padding = 0.6,
  point.padding = 0.4,
  segment.color = "grey50",
  segment.size = 0.3,
  segment.alpha = 0.8,
  max.overlaps = 12,
  seed = 42,
  min.segment.length = 0.2
)

# Annotation tibble ----
zone_annotations <- tibble(
  zone = c("crisis", "success"),
  xmin = c(-3.1, 0.9),
  xmax = c(-1.9, 2.1),
  ymin = c(1.8, -1.7),
  ymax = c(3.2, -0.7),
  x_text = c(-2.5, 1.5),
  y_text = c(2.5, -1.2),
  label = c(
    "CRISIS ZONE\nHigh burden +\nPoor surveillance",
    "SUCCESS ZONE\nLow burden +\nGood surveillance"
  ),
  color = c("darkred", "darkgreen")
)
```

5. Visualization Parameters

Show code 
```{r}
#| label: params
#| include: true
#| warning: false

### |-  plot aesthetics ----
colors <- get_theme_colors(
  palette = NULL
)

### |-  titles and caption ----
title_text <- str_glue("Global Measles Surveillance: Quality vs Disease Burden")

subtitle_text <- str_glue(
  "Analysis of {nrow(plot_data)} countries reveals surveillance systems under pressure. WHO Provisional Measles Data (as of June 2025)\n",
  "Countries cluster by outbreak characteristics and surveillance capabilities\n\n",
  "Note: Countries positioned closer together have similar outbreak patterns and surveillance capabilities\n",
  "All data is provisional - official annual figures available July 2025\n",
  "Model quality: {round(variance_explained * 100, 1)}% variance explained, Stress = {round(stress, 3)}"
)

caption_text <- create_social_caption(
  tt_year = 2025,
  tt_week = 25,
  source_text =  "WHO Provisional monthly measles and rubella data"
)

### |-  fonts ----
setup_fonts()
fonts <- get_font_families()

### |-  plot theme ----

# Start with base theme
base_theme <- create_base_theme(colors)

# Add weekly-specific theme elements
weekly_theme <- extend_weekly_theme(
  base_theme,
  theme(
    # Axis elements
    axis.text = element_text(color = colors$text, size = rel(0.7)),
    axis.title.x = element_text(color = colors$text, face = "bold", size = rel(0.8), margin = margin(t = 15)),
    axis.title.y = element_text(color = colors$text, face = "bold", size = rel(0.8), margin = margin(r = 10)),

    # Grid elements
    panel.grid.major.y = element_line(color = "gray50", linewidth = 0.05),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),

    # Legend elements
    legend.position = "plot",
    legend.title = element_text(family = fonts$tsubtitle, color = colors$text, size = rel(0.8), face = "bold"),
    legend.text = element_text(family = fonts$tsubtitle, color = colors$text, size = rel(0.7)),

    # Plot margins
    plot.margin = margin(t = 15, r = 15, b = 15, l = 15),
  )
)

# Set theme
theme_set(weekly_theme)
```

6. Plot

Show code 
```{r}
#| label: plot
#| warning: false

# Initial Plot ----
p <- ggplot(plot_data, aes(x = MDS1, y = MDS2)) +

  # Geoms
  geom_hline(yintercept = 0, color = "grey90", linewidth = 0.25) +
  geom_vline(xintercept = 0, color = "grey90", linewidth = 0.25) +
  geom_point(aes(color = lab_confirmation_rate, size = avg_incidence),
    alpha = 0.85, stroke = 0.2
  ) +
  # Scales
  scale_color_scico(
    name = "Lab Confirmation\nRate",
    labels = scales::percent_format(accuracy = 1),
    palette = "vik",
    direction = -1,
    trans = "identity",
    breaks = c(0, 0.25, 0.50, 0.75, 1.0),
    guide = guide_colorbar(
      title.position = "top",
      title.hjust = 0.5,
      barwidth = 10,
      barheight = 1,
      nbin = 100
    )
  ) +
  scale_size_continuous(
    range = c(1.5, 10),
    name = "Average Incidence (per million)",
    trans = "sqrt",
    breaks = c(1, 10, 100, 500, 2000),
    labels = c("1", "10", "100", "500", "2,000+"),
    guide = guide_legend(
      title.position = "top",
      title.hjust = 0.5,
      override.aes = list(alpha = 0.8, shape = 21),
      direction = "horizontal",
      nrow = 1
    )
  ) +
  # Labs
  labs(
    title = title_text,
    subtitle = subtitle_text,
    caption = caption_text,
    x = "MDS Dimension 1 →\n(Higher values = Lower disease burden)",
    y = "MDS Dimension 2 →\n(Higher values = More outbreak variability)",
  ) +
  # Theme
  theme(
    plot.title = element_text(
      size = rel(1.8),
      family = fonts$title,
      face = "bold",
      color = colors$title,
      lineheight = 1.1,
      margin = margin(t = 5, b = 10)
    ),
    plot.subtitle = element_text(
      size = rel(0.72),
      family = fonts$subtitle,
      color = alpha(colors$subtitle, 0.9),
      lineheight = 1.2,
      margin = margin(t = 5, b = 20)
    ),
    plot.caption = element_markdown(
      size = rel(0.55),
      family = fonts$caption,
      color = colors$caption,
      hjust = 0.5,
      margin = margin(t = 10)
    ),
    # Legend styling
    legend.position = "top",
    legend.box = "horizontal",
    legend.margin = margin(b = 10),
    legend.title = element_text(size = 11, face = "bold"),
    legend.text = element_text(size = 10),
    legend.spacing.x = unit(1, "cm"),
    legend.box.spacing = unit(0.5, "cm"),

    # Grid and axes
    panel.grid.minor = element_blank(),
    panel.grid.major = element_line(color = "grey95", linewidth = 0.3),
    axis.title = element_text(size = 12, face = "bold"),
    axis.text = element_text(size = 10),
  )

# Final Plot ----
p <- p +
  # top-left outliers  
  exec(geom_text_repel,
    data = outlier_datasets$top_left,
    mapping = aes(label = country),
    nudge_x = -0.5,
    nudge_y = -0.4,
    direction = "x",
    !!!base_label_style
  ) +

  # bottom-right outliers 
  exec(geom_text_repel,
    data = outlier_datasets$bottom_right,
    mapping = aes(label = country),
    !!!base_label_style
  ) +

  # zone rectangles
  pmap(zone_annotations, ~ annotate("rect",
    xmin = ..2, xmax = ..3,
    ymin = ..4, ymax = ..5,
    fill = colors$background,
    alpha = 0.7
  )) +

  # zone text labels
  pmap(zone_annotations, ~ annotate("text",
    x = ..6, y = ..7,
    label = ..8,
    hjust = 0.5, vjust = 0.5,
    size = 3.5, fontface = "bold",
    color = ..9
  ))
```

7. Save

Show code 
```{r}
#| label: save
#| warning: false

### |-  plot image ----  
save_plot(
  plot = p, 
  type = "tidytuesday", 
  year = 2025, 
  week = 25, 
  width  = 10,
  height = 8
)
```

8. Session Info

R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 22631)

Matrix products: default


locale:
[1] LC_COLLATE=English_United States.utf8 
[2] LC_CTYPE=English_United States.utf8   
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.utf8    

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base     

other attached packages:
 [1] here_1.0.1      ggrepel_0.9.6   scico_1.5.0     patchwork_1.3.0
 [5] glue_1.8.0      scales_1.3.0    janitor_2.2.0   showtext_0.9-7 
 [9] showtextdb_3.0  sysfonts_0.8.9  ggtext_0.1.2    lubridate_1.9.3
[13] forcats_1.0.0   stringr_1.5.1   dplyr_1.1.4     purrr_1.0.2    
[17] readr_2.1.5     tidyr_1.3.1     tibble_3.2.1    ggplot2_3.5.1  
[21] tidyverse_2.0.0 pacman_0.5.1   

loaded via a namespace (and not attached):
 [1] gtable_0.3.6       httr2_1.0.6        xfun_0.49          htmlwidgets_1.6.4 
 [5] gh_1.4.1           tzdb_0.5.0         vctrs_0.6.5        tools_4.4.0       
 [9] generics_0.1.3     parallel_4.4.0     curl_6.0.0         gifski_1.32.0-1   
[13] fansi_1.0.6        pkgconfig_2.0.3    lifecycle_1.0.4    compiler_4.4.0    
[17] farver_2.1.2       textshaping_0.4.0  munsell_0.5.1      codetools_0.2-20  
[21] snakecase_0.11.1   htmltools_0.5.8.1  yaml_2.3.10        crayon_1.5.3      
[25] pillar_1.9.0       camcorder_0.1.0    magick_2.8.5       commonmark_1.9.2  
[29] tidyselect_1.2.1   digest_0.6.37      stringi_1.8.4      labeling_0.4.3    
[33] rsvg_2.6.1         rprojroot_2.0.4    fastmap_1.2.0      grid_4.4.0        
[37] colorspace_2.1-1   cli_3.6.4          magrittr_2.0.3     utf8_1.2.4        
[41] withr_3.0.2        rappdirs_0.3.3     bit64_4.5.2        timechange_0.3.0  
[45] rmarkdown_2.29     tidytuesdayR_1.1.2 gitcreds_0.1.2     bit_4.5.0         
[49] ragg_1.3.3         hms_1.1.3          evaluate_1.0.1     knitr_1.49        
[53] markdown_1.13      rlang_1.1.6        gridtext_0.1.5     Rcpp_1.0.13-1     
[57] xml2_1.3.6         renv_1.0.3         vroom_1.6.5        svglite_2.1.3     
[61] rstudioapi_0.17.1  jsonlite_1.8.9     R6_2.5.1           systemfonts_1.1.0

9. GitHub Repository

The complete code for this analysis is available in tt_2025_25.qmd.

For the full repository, click here.

10. References

  1. Data Sources:
  • TidyTuesday 2025 Week 25: [Measles cases across the world](https://github.com/rfordatascience/tidytuesday/blob/main/data/2025/2025-06-24
Back to top