Pi-Change: Google Mobility Application

This vignette uses processed Google workplace mobility summaries for the United States, New York City, the United Kingdom, and London. Policy and pandemic event dates provide soft prior support for plausible mobility changes. The full raw Google mobility report is too large for a CRAN package, so the package includes the four-region daily summaries needed for this application.

library(PiChange)
library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 4.5.2

data("google_mobility_workplaces")

run_region <- function(region, event_dates) {
  region_data <- google_mobility_workplaces[
    google_mobility_workplaces$region == region,
  ]
  y <- region_data$workplaces_percent_change_from_baseline
  penalty <- construct_penalty(
    time = region_data$date,
    centers = as.Date(event_dates),
    width = 30,
    method = "mbic",
    family = "normal"
  )
  fit <- pi_change(y, penalty, min_seg_len = 30)
  list(
    data = region_data,
    fit = fit,
    centers = penalty$prior$center_indices,
    change_indices = changepoints(fit)
  )
}

us_events <- c("2020-03-15", "2020-05-01", "2020-12-11",
               "2021-03-11", "2022-08-16")
uk_events <- c("2020-03-23", "2020-06-08", "2020-10-31",
               "2020-12-02", "2020-12-08", "2020-12-26",
               "2021-03-29", "2022-02-24")

fits <- list(
  "United States" = run_region("United States", us_events),
  "New York City" = run_region("New York City", us_events),
  "United Kingdom" = run_region("United Kingdom", uk_events),
  "London" = run_region("London", uk_events)
)

data.frame(
  region = rep(names(fits), lengths(lapply(fits, `[[`, "change_indices"))),
  date = do.call(c, lapply(fits, function(x) x$data$date[x$change_indices]))
)
#>                         region       date
#> United States1   United States 2020-03-16
#> United States2   United States 2020-05-08
#> United States3   United States 2020-09-10
#> United States4   United States 2020-11-23
#> United States5   United States 2021-01-01
#> United States6   United States 2021-09-09
#> United States7   United States 2021-11-22
#> United States8   United States 2022-01-01
#> United States9   United States 2022-02-04
#> United States10  United States 2022-04-18
#> United States11  United States 2022-05-23
#> United States12  United States 2022-09-05
#> New York City1   New York City 2020-03-15
#> New York City2   New York City 2020-05-29
#> New York City3   New York City 2021-09-17
#> United Kingdom1 United Kingdom 2020-03-22
#> United Kingdom2 United Kingdom 2020-05-29
#> United Kingdom3 United Kingdom 2020-12-20
#> United Kingdom4 United Kingdom 2021-04-09
#> London1                 London 2020-03-16
#> London2                 London 2020-05-29
#> London3                 London 2020-12-20
#> London4                 London 2021-04-16

Detailed regional result

The standard fit plot shows the prior-informed penalty above the observed series and detected changes. It uses the same visual encoding as the other package examples.

plot(
  fits[["United States"]]$fit,
  ylab = "Workplace mobility percent change",
  title = "United States workplace mobility"
)

Four-region overview

The faceted overview uses the package palette and keeps the observed-series and segment-mean lines at half the width of the prior and detected-change lines.

plot_data <- do.call(rbind, lapply(names(fits), function(region) {
  x <- fits[[region]]
  data.frame(
    date = x$data$date,
    region = region,
    workplaces_percent_change_from_baseline =
      x$data$workplaces_percent_change_from_baseline,
    segment_mean = fitted(x$fit),
    pi_change = seq_len(nrow(x$data)) %in% x$change_indices,
    prior_center = seq_len(nrow(x$data)) %in% x$centers
  )
}))

ggplot(plot_data, aes(x = date, y = workplaces_percent_change_from_baseline)) +
  geom_line(
    aes(color = "Observed series", linetype = "Observed series"),
    linewidth = 0.4
  ) +
  geom_line(
    aes(y = segment_mean, color = "Segment mean", linetype = "Segment mean"),
    linewidth = 0.4
  ) +
  geom_vline(
    data = plot_data[plot_data$pi_change, ],
    aes(
      xintercept = date,
      color = "Detected change",
      linetype = "Detected change"
    ),
    linewidth = 0.8
  ) +
  geom_vline(
    data = plot_data[plot_data$prior_center, ],
    aes(
      xintercept = date,
      color = "Prior center",
      linetype = "Prior center"
    ),
    linewidth = 0.8
  ) +
  facet_wrap(~ region, ncol = 2) +
  scale_color_manual(
    values = c(
      "Observed series" = "#76BDE8",
      "Segment mean" = "#F06292",
      "Prior center" = "#D62728",
      "Detected change" = "#00C853"
    )
  ) +
  scale_linetype_manual(
    values = c(
      "Observed series" = "solid",
      "Segment mean" = "solid",
      "Prior center" = "dashed",
      "Detected change" = "solid"
    )
  ) +
  scale_x_date(date_breaks = "6 months", date_labels = "%b %Y") +
  labs(
    x = "Date",
    y = "Workplace mobility percent change",
    title = "Workplace mobility with PI-Change estimates",
    color = NULL,
    linetype = NULL
  ) +
  theme_minimal(base_size = 11) +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    legend.position = "top",
    legend.justification = "left",
    axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)
  )

The 30-position width represents approximately one month of daily observations, and the minimum segment length prevents changes less than 30 observations apart. Both choices should be varied in a sensitivity analysis when the substantive timing of prior events is uncertain.

Reference

Jacobs, J. and Chen, S. (2026). Pi-Change: A Prior-Informed Multiple Change Point Detection Algorithm. https://doi.org/10.48550/arXiv.2605.01003.