Dr. Blake,
When I initially tried to break up the year “pre and post freak out” I used the code below and changed the initial dates we set to pull the data from Media Cloud. I was concerned on the Distribution of Pairs when I broke down the time because there is a bar in the negative, but as you explained earlier today, it makes sense because the salience may have flipped between the two topics.
My other concern was my “mean” is negative, and in your discussion, that sounded like that shouldn’t happen.
| Descriptive Statistics: Pair Differences | ||||
| count | mean | sd | min | max |
|---|---|---|---|---|
| 36.000 | −1.056 | 27.815 | −112.000 | 54.000 |
So I ran it a second time, not shown below, and flipped my variables for this time and it was still weird. I hope that makes sense.
CODE: start_date <- “2025-10-19” end_date <- “2026-06-21”
############################################################
# STEP 1: Install and Load Required Packages
############################################################
required_packages <- c(
"tidyverse",
"httr2",
"lubridate",
"plotly",
"knitr",
"kableExtra"
)
for (pkg in required_packages) {
if (!requireNamespace(pkg, quietly = TRUE)) {
install.packages(pkg)
}
}
library(tidyverse)
library(httr2)
library(lubridate)
library(plotly)
library(knitr)
library(kableExtra)
############################################################
# STEP 2: Verify Media Cloud API Key
############################################################
if (!nzchar(Sys.getenv("MEDIACLOUD_KEY"))) {
stop(
"No Media Cloud API key was found.\n",
"Run the setup section once,\n",
"then restart R and try again."
)
}
mc_key <- Sys.getenv("MEDIACLOUD_KEY")
############################################################
# STEP 3: Define the Date Range
############################################################
start_date <- "2025-10-19"
end_date <- "2026-06-21"
############################################################
# STEP 4: Define Search Queries
############################################################
vaccine_query <-
"((measles OR mumps OR pertussis OR whooping cough) AND
(outbreak* OR cases OR vaccine* OR infection* OR transmission OR MMR OR measles vaccine))"
food_query <-
"((cyclospora OR salmonella OR listeria OR e.coli)
AND (food recall* OR recalled food OR recall OR recalls OR contamination OR foodborne illness OR food safety OR produce))"
############################################################
# STEP 5: Create Topic Definitions
############################################################
topics <- c(
"Vaccine-preventable disease outbreaks" = vaccine_query,
"Food recall and contamination events" = food_query
)
############################################################
# STEP 6: Create a Function to Retrieve Daily
# Story Counts
############################################################
mc_counts_over_time <- function(
query,
start_date,
end_date,
collection_id = 34412234
) {
resp <-
request(
"https://search.mediacloud.org/api/search/count-over-time"
) |>
req_headers(
Authorization = paste("Token", mc_key)
) |>
req_url_query(
q = query,
start = start_date,
end = end_date,
platform = "onlinenews-mediacloud",
cs = collection_id
) |>
req_perform()
results <- resp |>
resp_body_json()
map_dfr(
results$count_over_time$counts,
as_tibble
) |>
mutate(
date = as.Date(date)
)
}
############################################################
# STEP 7: Download Daily Counts for All Topics
############################################################
topic_data <-
imap_dfr(
topics,
function(search_string, topic_name) {
message(
"Downloading: ",
topic_name
)
results <-
mc_counts_over_time(
query = search_string,
start_date = start_date,
end_date = end_date
)
Sys.sleep(35)
results |>
mutate(
topic = topic_name
)
}
)
############################################################
# STEP 8: Examine the Downloaded Daily Data
############################################################
head(topic_data)
glimpse(topic_data)
############################################################
# STEP 9: Convert Daily Counts to Weekly Counts
############################################################
weekly_counts <-
topic_data |>
mutate(
week = floor_date(
date,
unit = "week"
)
) |>
group_by(
topic,
week
) |>
summarise(
stories = sum(count),
.groups = "drop"
)
############################################################
# STEP 10: Examine Weekly Counts
############################################################
weekly_counts
############################################################
# STEP 11: Summarize Coverage by Topic
############################################################
topic_summary <-
weekly_counts |>
group_by(topic) |>
summarise(
Total_Stories = sum(stories),
Weekly_Minimum = min(stories),
Weekly_Maximum = max(stories),
Weekly_Median = median(stories),
Weekly_Mean = round(
mean(stories),
2
),
Weekly_SD = round(
sd(stories),
2
),
.groups = "drop"
) |>
arrange(
desc(Total_Stories)
)
Table <-
topic_summary |>
kbl(
caption =
"Media Cloud Topic Summary Statistics",
align = "lrrrrrr"
) |>
kable_styling(
full_width = FALSE,
bootstrap_options = c(
"striped",
"hover",
"condensed"
)
)
Table
############################################################
# STEP 12: Create an Interactive Plotly Graph
############################################################
weekly_counts <- weekly_counts |>
mutate(
stories_formatted = formatC(
stories,
format = "d",
big.mark = ","
)
)
Graph <- plot_ly(
data = weekly_counts,
x = ~week,
y = ~stories,
color = ~topic,
colors = "Set1",
type = "scatter",
mode = "lines+markers",
customdata = ~topic,
text = ~stories_formatted,
hovertemplate =
paste(
"%{customdata}
",
"Week: %{x|%Y-%m-%d}
",
"Stories: %{text}
",
""
)
) |>
layout(
title =
"Weekly U.S. Coverage Volume",
hovermode = "x unified",
xaxis = list(
title = "Week"
),
yaxis = list(
title = "Weekly Story Count",
separatethousands = TRUE
)
)
Graph
weekly_counts_wide <-
weekly_counts %>%
mutate(
topic = make.names(topic)
) %>%
select(week, topic, stories) %>%
pivot_wider(
names_from = topic,
values_from = stories
)
# ============================================================
# Setup: Install and Load Required Packages
# ============================================================
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("plotly")) install.packages("plotly")
if (!require("gt")) install.packages("gt")
if (!require("gtExtras")) install.packages("gtExtras")
if (!require("broom")) install.packages("broom")
library(tidyverse)
library(plotly)
library(gt)
library(gtExtras)
library(broom)
options(scipen = 999)
# ============================================================
# Data Import
# ============================================================
mydata <- weekly_counts_wide
# Specify the two variables involved
mydata$V1 <- mydata$Vaccine.preventable.disease.outbreaks
mydata$V2 <- mydata$Food.recall.and.contamination.events
# ============================================================
# Compute Pair Differences
# ============================================================
mydata$PairDifferences <- mydata$V2 - mydata$V1
# ============================================================
# Interactive Histogram of Pair Differences
# ============================================================
hist_plot <- plot_ly(
data = mydata,
x = ~PairDifferences,
type = "histogram",
marker = list(color = "#1f78b4", line = list(color = "black", width = 1))
) %>%
layout(
title = "Distribution of Pair Differences",
xaxis = list(title = "Pair Differences"),
yaxis = list(title = "Count"),
shapes = list(
list(
type = "line",
x0 = mean(mydata$PairDifferences, na.rm = TRUE),
x1 = mean(mydata$PairDifferences, na.rm = TRUE),
y0 = 0,
y1 = max(table(mydata$PairDifferences)),
line = list(color = "red", dash = "dash")
)
)
)
# ============================================================
# Descriptive Statistics
# ============================================================
desc_stats <- mydata %>%
summarise(
count = n(),
mean = mean(PairDifferences, na.rm = TRUE),
sd = sd(PairDifferences, na.rm = TRUE),
min = min(PairDifferences, na.rm = TRUE),
max = max(PairDifferences, na.rm = TRUE)
)
desc_table <- desc_stats %>%
gt() %>%
gt_theme_538() %>%
tab_header(title = "Descriptive Statistics: Pair Differences") %>%
fmt_number(columns = where(is.numeric), decimals = 3)
# ============================================================
# Normality Test (Shapiro-Wilk)
# ============================================================
shapiro_res <- shapiro.test(mydata$PairDifferences)
shapiro_table <- tidy(shapiro_res) %>%
select(statistic, p.value, method) %>%
gt() %>%
gt_theme_538() %>%
tab_header(title = "Normality Test (Shapiro-Wilk)") %>%
fmt_number(columns = c(statistic, p.value), decimals = 4) %>%
tab_source_note(
source_note = "If the P.VALUE is 0.05 or less, the number of pairs is fewer than 40, and the distribution of pair differences shows obvious non-normality or outliers, consider using the Wilcoxon Signed Rank Test results instead of the Paired-Samples t-Test results."
)
# ============================================================
# Reshape Data for Repeated-Measures Plot
# ============================================================
df_long <- mydata %>%
pivot_longer(cols = c(V1, V2),
names_to = "Measure",
values_to = "Value")
# ============================================================
# Repeated-Measures Boxplot (Interactive, with Means)
# ============================================================
group_means <- df_long %>%
group_by(Measure) %>%
summarise(mean_value = mean(Value), .groups = "drop")
boxplot_measures <- plot_ly() %>%
add_trace(
data = df_long,
x = ~Measure, y = ~Value,
type = "box",
boxpoints = "outliers",
marker = list(color = "red", size = 4),
line = list(color = "black"),
fillcolor = "royalblue",
name = ""
) %>%
add_trace(
data = group_means,
x = ~Measure, y = ~mean_value,
type = "scatter", mode = "markers",
marker = list(
symbol = "diamond", size = 9,
color = "black", line = list(color = "white", width = 1)
),
text = ~paste0("Mean = ", round(mean_value, 2)),
hoverinfo = "text",
name = "Group Mean"
) %>%
layout(
title = "Boxplot of Repeated Measures (V1 vs V2) with Means",
xaxis = list(title = "Measure"),
yaxis = list(title = "Value"),
showlegend = FALSE
)
# ============================================================
# Parametric Test (Paired-Samples t-Test)
# ============================================================
t_res <- t.test(mydata$V2, mydata$V1, paired = TRUE)
t_table <- tidy(t_res) %>%
select(statistic, parameter, p.value, conf.low, conf.high, method) %>%
gt() %>%
gt_theme_538() %>%
tab_header(title = "Paired-Samples t-Test") %>%
fmt_number(columns = c(statistic, p.value, conf.low, conf.high), decimals = 4)
t_summary <- mydata %>%
select(V1, V2) %>%
summarise_all(list(Mean = mean, SD = sd)) %>%
gt() %>%
gt_theme_538() %>%
tab_header(title = "Group Means and SDs (t-Test)") %>%
fmt_number(columns = everything(), decimals = 3)
# ============================================================
# Nonparametric Test (Wilcoxon Signed Rank)
# ============================================================
wilcox_res <- wilcox.test(mydata$V1, mydata$V2, paired = TRUE,
exact = FALSE)
wilcox_table <- tidy(wilcox_res) %>%
select(statistic, p.value, method) %>%
gt() %>%
gt_theme_538() %>%
tab_header(title = "Wilcoxon Signed Rank Test") %>%
fmt_number(columns = c(statistic, p.value), decimals = 4)
wilcox_summary <- mydata %>%
select(V1, V2) %>%
summarise_all(list(Mean = mean, SD = sd)) %>%
gt() %>%
gt_theme_538() %>%
tab_header(title = "Group Means and SDs (Wilcoxon)") %>%
fmt_number(columns = everything(), decimals = 3)
# ============================================================
# Results Summary (in specified order)
# ============================================================
hist_plot
desc_table
shapiro_table
boxplot_measures
t_table
t_summary
wilcox_table
wilcox_summary