1.1 What Are We Measuring?

Particulate matter (PM) is a mixture of solid particles and liquid droplets suspended in air. Regulatory agencies classify particles by their aerodynamic diameter:

The EPA Air Quality Index (AQI) categories for PM2.5 (24-h average) are:

1.2 Sensor Network

Students deployed Raspberry Pi Zero units equipped with Enviro+ boards (PMS5003 optical particle sensor + BME280 meteorological sensor + MICS6814 gas sensor) at multiple campus locations. Each unit logs:

• pm1_ug_m3, pm25_ug_m3, pm10_ug_m3 — particulate concentrations (µg/m³)
• temperature_C, humidity_pct, pressure_hPa — meteorological variables
• light_lux — ambient light
• oxidised_kO, reduced_kO, nh3_kO — gas-phase proxies (oxidising gases, reducing gases, ammonia)

1.3 Research Questions

1.4 Hypotheses

2.1 Packages

Run the install block once in the RStudio Console — not inside the Rmd.

pkgs <- c("tidyverse", "lubridate", "patchwork", "scales",
          "slider",    "ggridges",  "moments",
          "broom",     "broom.mixed", "knitr", "kableExtra",
          "lme4",      "lmerTest",  "nlme")
to_install <- pkgs[!pkgs %in% rownames(installed.packages())]
if (length(to_install)) install.packages(to_install, quiet = TRUE)

library(tidyverse); library(lubridate); library(patchwork); library(scales)
library(slider);    library(ggridges);  library(moments)
library(broom);     library(knitr);     library(kableExtra)

theme_set(theme_minimal(base_size = 13))

# Color palette — scales automatically to however many sensors are present
sensor_pal <- c(
  "#0072B2", "#D55E00", "#009E73", "#CC79A7",
  "#E69F00", "#56B4E9", "#F0E442", "#999999",
  "#882255", "#44AA99"
)

# pal is built after df exists; placeholder here, rebuilt in Section 6
pal <- sensor_pal   # overwritten once df is loaded

3.1 Correction Strategy

At deployment, each team member recorded:

1. The actual wall-clock time (PST) when they powered the unit on.
2. The time displayed on the Pi’s LCD at that same moment.

The offset is simply:

\[\text{offset} = t_{\text{real}} - t_{\text{Pi}}\]

This constant is then added to every raw timestamp in that sensor’s CSV.

3.2 Calibration Table

# sensor_calibration.R must live in the same directory as this Rmd.
# It creates the `calibration` tibble with one row per sensor (two for PiZ15).
source("sensor_calibration.R")

## Calibration table:
## 
## 
## sensor   team_member   location                                              offset_label       reboot  notes                                                   
## -------  ------------  ----------------------------------------------------  ----------------  -------  --------------------------------------------------------
## PiZ1     NA            NA                                                    NA                      1  No calibration data recorded                            
## PiZ2     NA            NA                                                    NA                      1  No calibration data recorded                            
## PiZ3     NA            NA                                                    NA                      1  No calibration data recorded                            
## PiZ4     NA            NA                                                    NA                      1  No calibration data recorded                            
## PiZ5     NA            NA                                                    NA                      1  No calibration data recorded                            
## PiZ6     NA            NA                                                    NA                      1  No calibration data recorded                            
## PiZ7     Naomi         Outside Blaisdell dorm, by the back stairs            10d 05h 33m 00s         1  NA                                                      
## PiZ8     Sorren        Outside Wig Lounge patio, side facing Lyon            9d 22h 26m 00s          1  NA                                                      
## PiZ9     NA            NA                                                    NA                      1  No calibration data recorded                            
## PiZ10    NA            NA                                                    NA                      1  No calibration data recorded                            
## PiZ11    NA            NA                                                    NA                      1  No calibration data recorded                            
## PiZ12    Ashley        Mudd-Blaisdell Courtyard                              7d 09h 46m 02s          1  NA                                                      
## PiZ13    Mariana       Outside Harwood by the bicycle racks                  10d 00h 31m 17s         1  NA                                                      
## PiZ14    NA            NA                                                    NA                      1  No calibration data recorded                            
## PiZ15    Dessa         Outside Seaver Theater, second floor                  4d 20h 40m 22s          1  First power-on                                          
## PiZ15    Dessa         Outside Seaver Theater, second floor                  14d 02h 51m 23s         2  Second power-on / reboot                                
## PiZ2W1   Kate          Outside 2F Oldenborg window facing Frary              9d 16h 22m 32s          1  Sheet entry '2/9' interpreted as 3/9 — verify with Kate 
## PiZ2W2   Matteo        Outside bathroom window, first floor (building TBD)   13d 15h 28m 10s         1  NA                                                      
## PiZ2W3   Matthew       Outside ground floor, Mudd-Honnold Library            9d 20h 50m 30s          1  NA                                                      
## PiZ2W4   NA            NA                                                    NA                      1  No calibration data recorded                            
## PiZ2W5   Anna          Pomona Organic Farm (near dome)                       9d 12h 35m 54s          1  NA                                                      
## 
## 20 sensors total | 9 with clock corrections | 11 missing data

clock_offsets <- calibration |>
  filter(reboot == 1) |>
  select(sensor, offset_sec)

calibration |>
  filter(reboot == 1) |>
  select(sensor, team_member, location, offset_label, notes) |>
  kable(caption = "Deployment metadata and clock offsets (reboot 1 only)") |>
  kable_styling(full_width = FALSE) |>
  column_spec(3, width = "22em") |>
  row_spec(which(is.na(calibration$offset_label[calibration$reboot == 1])),
           color = "#aaa")

4.1 Reading the CSVs

We auto-discover every *_ea30_sp26_v05.csv file in data_dir so no manual list is needed — and no trailing-comma errors.

# ── Set to the folder containing your CSV files ──────────────────────────────
data_dir <- "/home/mwl04747/RTricks/A08_Raspberry_Pi_Analysis/Data_S2026"

# ── Helper: read one CSV and apply the clock correction ──────────────────────
read_sensor <- function(path, sensor_id, offsets_tbl) {

  offset_s <- offsets_tbl |>
    filter(sensor == sensor_id) |>
    pull(offset_sec)

  if (length(offset_s) == 0 || is.na(offset_s)) {
    warning(sprintf("No clock offset for %s — timestamps uncorrected.", sensor_id))
    offset_s <- 0
  }

  read_csv(path, show_col_types = FALSE) |>
    mutate(
      sensor      = sensor_id,
      datetime_pi = ymd_hms(timestamp, tz = "America/Los_Angeles"),
      datetime    = datetime_pi + seconds(offset_s),      # ← corrected
      date        = as_date(datetime),
      hour        = hour(datetime),
      wday        = wday(datetime, label = TRUE, abbr = TRUE),
      tod         = cut(hour,
                        breaks = c(-1, 5, 11, 17, 20, 23),
                        labels = c("Night (0–5)", "Morning (6–11)",
                                   "Afternoon (12–17)", "Evening (18–20)",
                                   "Night (21–23)"))
    ) |>
    select(-timestamp, -datetime_pi)
}

# ── Auto-discover CSVs ────────────────────────────────────────────────────────
csv_paths  <- list.files(data_dir, pattern = "_ea30_sp26_v05\\.csv$",
                         full.names = TRUE)
sensor_ids <- sub("_ea30_sp26_v05\\.csv$", "", basename(csv_paths))

if (length(csv_paths) == 0)
  stop("No *_ea30_sp26_v05.csv files found. Set data_dir correctly.")

files <- setNames(as.list(csv_paths), sensor_ids)
message(sprintf("Discovered %d file(s): %s",
                length(files), paste(sensor_ids, collapse = ", ")))

raw <- imap_dfr(files, ~ read_sensor(.x, .y, clock_offsets))

glimpse(raw)

## Rows: 15,569
## Columns: 16
## $ temperature_C <dbl> 9.32, 15.78, 12.18, 18.49, 18.38, 18.09, 17.70, 17.89, 1…
## $ pressure_hPa  <dbl> 707.11, 707.11, 697.76, 697.76, 973.24, 973.36, 973.48, …
## $ humidity_pct  <dbl> 74.60, 74.60, 75.24, 75.24, 40.73, 41.67, 42.57, 42.87, …
## $ light_lux     <dbl> 260.37, 152.57, 3.85, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00…
## $ oxidised_kO   <dbl> 22.17, 116.37, 23.83, 11.53, 2.37, 3.78, 5.11, 6.38, 7.6…
## $ reduced_kO    <dbl> 2.67, 7.70, 1.51, 13.72, 624.66, 614.48, 596.72, 578.23,…
## $ nh3_kO        <dbl> 0.67, 2.91, 0.23, 1.18, 70.73, 129.96, 165.98, 187.78, 2…
## $ pm1_ug_m3     <dbl> 2, 1, 7, 5, 4, 5, 5, 5, 5, 5, 5, 4, 5, 5, 5, 5, 5, 5, 5,…
## $ pm25_ug_m3    <dbl> 2, 1, 8, 5, 5, 6, 7, 6, 6, 6, 6, 5, 6, 6, 7, 6, 6, 6, 6,…
## $ pm10_ug_m3    <dbl> 2, 1, 8, 5, 5, 6, 8, 7, 7, 7, 6, 6, 6, 9, 9, 9, 8, 8, 8,…
## $ sensor        <chr> "PiZ11", "PiZ11", "PiZ12", "PiZ12", "PiZ12", "PiZ12", "P…
## $ datetime      <dttm> 2026-02-27 14:57:43, 2026-02-27 15:02:17, 2026-03-09 01…
## $ date          <date> 2026-02-27, 2026-02-27, 2026-03-09, 2026-03-09, 2026-03…
## $ hour          <int> 14, 15, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, …
## $ wday          <ord> Fri, Fri, Mon, Mon, Mon, Mon, Mon, Mon, Mon, Mon, Mon, M…
## $ tod           <fct> Afternoon (12–17), Afternoon (12–17), Night (0–5), Night…

4.2 Quality Control

We apply three QC filters:

1. Timestamp present — rows with no parseable datetime are removed.
2. PM2.5 present — rows with NA in the target variable are removed.
3. Minimum observations — sensors with fewer than 10 valid PM2.5 readings are excluded from analysis (e.g. a unit that only logged one row).

# Sensors with enough data
valid_sensors <- raw |>
  filter(!is.na(pm25_ug_m3)) |>
  count(sensor) |>
  filter(n >= 10) |>
  pull(sensor)

dropped <- setdiff(sensor_ids, valid_sensors)
if (length(dropped) > 0)
  message("Dropped (< 10 obs): ", paste(dropped, collapse = ", "))

df <- raw |>
  filter(sensor %in% valid_sensors,
         !is.na(pm25_ug_m3),
         !is.na(datetime))

# Attach location label for plot facets / legends
df <- df |>
  left_join(
    calibration |>
      filter(reboot == 1) |>
      select(sensor, location),
    by = "sensor"
  ) |>
  mutate(
    sensor   = factor(sensor),
    # Short label: just the location, truncated
    loc_label = str_trunc(coalesce(location, sensor), 35)
  )

Final dataset: 15559 rows across 10 sensor(s).

# Build named color palette from the sensors actually present in df
n_sensors <- n_distinct(df$sensor)
pal        <- setNames(sensor_pal[seq_len(n_sensors)], levels(df$sensor))

4.3 Coverage Summary

df |>
  group_by(sensor, loc_label) |>
  summarise(
    n_obs      = n(),
    start      = min(date),
    end        = max(date),
    days_span  = as.integer(end - start) + 1L,
    pct_missing_pm25 = mean(is.na(pm25_ug_m3)) * 100,
    .groups = "drop"
  ) |>
  kable(digits = 1,
        caption = "Data coverage per sensor after QC",
        col.names = c("Sensor","Location","N obs","Start","End",
                      "Days spanned","% PM2.5 missing")) |>
  kable_styling(full_width = FALSE)

5.1 Summary Table

desc <- df |>
  group_by(sensor) |>
  summarise(
    N        = n(),
    Mean     = mean(pm25_ug_m3),
    SD       = sd(pm25_ug_m3),
    Median   = median(pm25_ug_m3),
    IQR      = IQR(pm25_ug_m3),
    Skewness = moments::skewness(pm25_ug_m3),
    P5       = quantile(pm25_ug_m3, .05),
    P95      = quantile(pm25_ug_m3, .95),
    Max      = max(pm25_ug_m3),
    `% Good (≤12)`       = mean(pm25_ug_m3 <= 12)  * 100,
    `% Moderate (12–35)` = mean(pm25_ug_m3 > 12 & pm25_ug_m3 <= 35.4) * 100,
    `% USG (>35)`        = mean(pm25_ug_m3 > 35.4) * 100,
    .groups = "drop"
  )

kable(desc, digits = 2,
      caption = "PM2.5 descriptive statistics (µg/m³) with EPA AQI breakdown") |>
  kable_styling(full_width = TRUE) |>
  add_header_above(c(" " = 10, "EPA AQI Category (%)" = 3))

5.2 PM1.0 vs PM2.5 — Are They Telling the Same Story?

df |>
  select(sensor, pm1_ug_m3, pm25_ug_m3) |>
  pivot_longer(-sensor, names_to = "fraction", values_to = "conc") |>
  mutate(fraction = recode(fraction,
                           pm1_ug_m3  = "PM1.0",
                           pm25_ug_m3 = "PM2.5")) |>
  group_by(sensor, fraction) |>
  summarise(Mean = mean(conc), Median = median(conc),
            SD = sd(conc), .groups = "drop") |>
  pivot_wider(names_from = fraction,
              values_from = c(Mean, Median, SD)) |>
  mutate(`PM1/PM2.5 ratio (mean)` = `Mean_PM1.0` / `Mean_PM2.5`) |>
  kable(digits = 2,
        caption = "PM1.0 vs PM2.5 summary — a ratio near 1 means nearly all fine
                   particles are in the sub-1µm range (combustion dominated)") |>
  kable_styling(full_width = FALSE)

6.1 Time Series with Rolling Mean

# Disabled: requires ggplot2 >= 3.4 with matching scales/vctrs on R 4.2.2
# Upgrade R to 4.4.x to re-enable this chunk

6.2 PM1.0 and PM2.5 Together

# Disabled: same ggplot2/R 4.2.2 version conflict as timeseries chunk

6.3 Distributions

p_dens <- df |>
  ggplot(aes(x = pm25_ug_m3, fill = sensor, color = sensor)) +
  geom_density(alpha = 0.30, linewidth = 0.7) +
  geom_vline(xintercept = c(12, 35.4),
             linetype = "dashed",
             color    = c("forestgreen", "darkorange"),
             linewidth = 0.7) +
  scale_fill_manual(values  = pal) +
  scale_color_manual(values = pal) +
  labs(title = "PM2.5 Density", x = "PM2.5 (µg/m³)",
       y = "Density", fill = "Sensor", color = "Sensor") +
  theme(legend.position = "top")

p_box <- df |>
  ggplot(aes(x = sensor, y = pm25_ug_m3, fill = sensor)) +
  geom_boxplot(outlier.alpha = 0.25, outlier.size = 0.8, width = 0.5) +
  geom_jitter(width = 0.12, alpha = 0.04, size = 0.5) +
  stat_summary(fun = mean, geom = "point",
               shape = 23, fill = "white", size = 2.5) +
  scale_fill_manual(values = pal) +
  labs(title    = "PM2.5 by Sensor",
       subtitle = "◇ = mean; box = IQR",
       x = "Sensor", y = "PM2.5 (µg/m³)") +
  theme(legend.position = "none")

p_dens + p_box

6.4 Diurnal Pattern (Time-of-Day)

df |>
  group_by(sensor, hour) |>
  summarise(
    mean_pm25 = mean(pm25_ug_m3),
    se_pm25   = sd(pm25_ug_m3) / sqrt(n()),
    mean_pm1  = mean(pm1_ug_m3),
    se_pm1    = sd(pm1_ug_m3)  / sqrt(n()),
    .groups   = "drop"
  ) |>
  pivot_longer(c(mean_pm25, mean_pm1),
               names_to = "fraction", values_to = "mean") |>
  mutate(
    se       = if_else(fraction == "mean_pm25", se_pm25, se_pm1),
    fraction = recode(fraction, mean_pm25 = "PM2.5", mean_pm1 = "PM1.0")
  ) |>
  ggplot(aes(x = hour, y = mean, color = sensor, linetype = fraction)) +
  geom_ribbon(aes(ymin = mean - se, ymax = mean + se,
                  fill = sensor, group = interaction(sensor, fraction)),
              alpha = 0.10, color = NA) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~ fraction) +
  scale_color_manual(values = pal) +
  scale_fill_manual(values  = pal) +
  scale_x_continuous(breaks = seq(0, 23, 3),
                     labels = sprintf("%02d:00", seq(0, 23, 3))) +
  labs(title    = "Diurnal Profile of PM1.0 and PM2.5",
       subtitle = "Mean ± 1 SE by hour of day",
       x = "Hour of Day", y = "Concentration (µg/m³)",
       color = "Sensor", fill = "Sensor", linetype = NULL) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "top")

What to look for: Morning peaks (~7–9 AM) or evening peaks (~5–8 PM) would support H2. A midday dip is consistent with stronger atmospheric mixing (higher planetary boundary layer) diluting surface-level emissions.

6.5 Date × Hour Heatmap

# Disabled: scale_y_date affected by same ggplot2/R 4.2.2 version conflict

6.6 Correlation Matrix

These heatmaps show pairwise Pearson correlations between PM2.5 and all co-measured variables. They help identify which predictors to include in the regression model.

cor_vars <- c("pm25_ug_m3", "pm1_ug_m3", "pm10_ug_m3",
              "temperature_C", "humidity_pct", "pressure_hPa",
              "light_lux", "oxidised_kO", "reduced_kO", "nh3_kO")

for (s in levels(df$sensor)) {
  cat(sprintf("\n\n#### Sensor: %s\n\n", s))

  cor_mat <- df |>
    filter(sensor == s) |>
    select(all_of(cor_vars)) |>
    drop_na() |>
    cor(method = "pearson")

  p <- cor_mat |>
    as.data.frame() |>
    rownames_to_column("var1") |>
    pivot_longer(-var1, names_to = "var2", values_to = "r") |>
    mutate(
      var1  = factor(var1, levels = cor_vars),
      var2  = factor(var2, levels = rev(cor_vars)),
      label = sprintf("%.2f", r)
    ) |>
    ggplot(aes(x = var1, y = var2, fill = r)) +
    geom_tile(color = "white", linewidth = 0.5) +
    geom_text(aes(label = label), size = 3, color = "black") +
    scale_fill_gradient2(
      low  = "#D55E00", mid = "white", high = "#0072B2",
      midpoint = 0, limits = c(-1, 1), name = "r"
    ) +
    labs(title = paste("Correlation Matrix —", s), x = NULL, y = NULL) +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 9),
          axis.text.y = element_text(size = 9),
          plot.title  = element_text(size = 11))

  print(p)
}

7.1 Are the Data Normally Distributed?

set.seed(42)
norm_res <- df |>
  group_by(sensor) |>
  summarise(
    n       = n(),
    W_pm25  = shapiro.test(sample(pm25_ug_m3, min(n(), 5000)))$statistic,
    p_pm25  = shapiro.test(sample(pm25_ug_m3, min(n(), 5000)))$p.value,
    W_pm1   = shapiro.test(sample(pm1_ug_m3,  min(n(), 5000)))$statistic,
    p_pm1   = shapiro.test(sample(pm1_ug_m3,  min(n(), 5000)))$p.value,
    .groups = "drop"
  ) |>
  mutate(
    normal_pm25 = if_else(p_pm25 > 0.05, "Yes", "No"),
    normal_pm1  = if_else(p_pm1  > 0.05, "Yes", "No")
  )

kable(norm_res, digits = 4,
      caption = "Shapiro-Wilk normality test (sampled to ≤ 5000 obs)") |>
  kable_styling(full_width = FALSE)

7.2 Do Sensors Differ? — Wilcoxon Rank-Sum Test

The Wilcoxon rank-sum test (Mann-Whitney U) is the non-parametric equivalent of the two-sample t-test. It tests whether the distributions of two groups have the same location (median), without assuming normality.

Hypotheses: H₀: PM2.5 at sensor A = PM2.5 at sensor B (same distribution); H₁: they differ.

sensors_v <- levels(df$sensor)
combs     <- combn(sensors_v, 2, simplify = FALSE)

wilcox_res <- map_dfr(combs, function(pair) {
  x  <- df |> filter(sensor == pair[1]) |> pull(pm25_ug_m3)
  y  <- df |> filter(sensor == pair[2]) |> pull(pm25_ug_m3)
  wt <- wilcox.test(x, y, conf.int = TRUE, exact = FALSE)
  n1 <- length(x); n2 <- length(y)
  r  <- 1 - 2 * wt$statistic / (n1 * n2)   # rank-biserial effect size
  tibble(
    Comparison      = paste(pair[1], "vs", pair[2]),
    W               = wt$statistic,
    p_value         = wt$p.value,
    Location_shift  = wt$estimate,
    CI_lower        = wt$conf.int[1],
    CI_upper        = wt$conf.int[2],
    Effect_r        = r,
    Significant     = if_else(p_value < 0.05, "✓", "—")
  )
})

wilcox_res |>
  mutate(Significant = cell_spec(
    Significant, bold = TRUE,
    color = if_else(Significant == "✓", "forestgreen", "grey60")
  )) |>
  kable(digits = 4, escape = FALSE,
        caption = "Pairwise Wilcoxon tests for PM2.5; effect size r = rank-biserial correlation") |>
  kable_styling(full_width = FALSE)

7.3 Does PM2.5 Vary by Time of Day? — Kruskal-Wallis Test

kw_res <- df |>
  group_by(sensor) |>
  summarise(
    chi_sq  = kruskal.test(pm25_ug_m3 ~ tod)$statistic,
    df_kw   = kruskal.test(pm25_ug_m3 ~ tod)$parameter,
    p_value = kruskal.test(pm25_ug_m3 ~ tod)$p.value,
    .groups = "drop"
  ) |>
  mutate(significant = if_else(p_value < 0.05, "✓", "—"),
         # Epsilon-squared effect size for KW
         eta_sq = (chi_sq - df_kw + 1) / (n_distinct(df$sensor) - df_kw))

kable(kw_res, digits = 4,
      caption = "Kruskal-Wallis test: does PM2.5 differ across time-of-day groups?") |>
  kable_styling(full_width = FALSE)

df |>
  ggplot(aes(x = tod, y = pm25_ug_m3, fill = tod)) +
  geom_boxplot(outlier.alpha = 0.15, outlier.size = 0.6, width = 0.6) +
  stat_summary(fun = mean, geom = "point",
               shape = 23, fill = "white", size = 2) +
  facet_wrap(~ sensor, nrow = 1) +
  scale_fill_brewer(palette = "RdYlGn", direction = -1) +
  labs(title    = "PM2.5 by Time of Day",
       subtitle = "Kruskal-Wallis test p-values shown above",
       x = NULL, y = "PM2.5 (µg/m³)") +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 35, hjust = 1))

8.1 What Predicts PM2.5? — Multiple Linear Regression

Even though PM2.5 is not normally distributed, OLS regression can still be used to identify predictors, because the assumptions concern the residuals, not the raw response. We will check residuals below.

Model:
\[\text{PM}_{2.5} = \beta_0 + \beta_1 T + \beta_2 \text{RH} + \beta_3 P + \beta_4 \text{light} + \beta_5 \text{ox} + \beta_6 \text{red} + \beta_7 \text{NH}_3 + \varepsilon\]

lm_list <- df |>
  group_by(sensor) |>
  group_map(function(d, key) {
    mod <- lm(pm25_ug_m3 ~ temperature_C + humidity_pct + pressure_hPa +
                light_lux + oxidised_kO + reduced_kO + nh3_kO,
              data = d)
    list(sensor = key$sensor, model = mod,
         tidy   = tidy(mod)   |> mutate(sensor = key$sensor, .before=1),
         glance = glance(mod) |> mutate(sensor = key$sensor, .before=1))
  })

# Coefficient table
map_dfr(lm_list, "tidy") |>
  mutate(sig = case_when(p.value < 0.001 ~ "***",
                         p.value < 0.01  ~ "**",
                         p.value < 0.05  ~ "*",
                         TRUE            ~ "")) |>
  kable(digits = 4,
        caption = "OLS regression coefficients (*** p<.001, ** p<.01, * p<.05)") |>
  kable_styling(full_width = TRUE)

map_dfr(lm_list, "glance") |>
  select(sensor, r.squared, adj.r.squared,
         statistic, p.value, nobs) |>
  kable(digits = 4, caption = "Model fit statistics") |>
  kable_styling(full_width = FALSE)

8.2 Residual Diagnostics

Check the four key OLS assumptions visually:

for (item in lm_list) {
  cat(sprintf("\n\n#### Sensor: %s\n\n", item$sensor))
  par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
  plot(item$model, main = item$sensor)
  par(mfrow = c(1, 1))
}

8.3 Log-Transformed Response

PM concentrations are strictly positive and right-skewed — a natural candidate for log-transformation, which often improves residual normality.

lm_log_list <- df |>
  filter(pm25_ug_m3 > 0) |>    # log(0) = -Inf
  group_by(sensor) |>
  group_map(function(d, key) {
    mod <- lm(log(pm25_ug_m3) ~ temperature_C + humidity_pct + pressure_hPa +
                light_lux + oxidised_kO + reduced_kO + nh3_kO,
              data = d)
    list(sensor = key$sensor, model = mod,
         glance = glance(mod) |> mutate(sensor = key$sensor, .before=1))
  })

map_dfr(lm_log_list, "glance") |>
  select(sensor, r.squared, adj.r.squared, nobs) |>
  kable(digits = 4, caption = "Log-model R² per sensor") |>
  kable_styling(full_width = FALSE)

8.4 Temporal Autocorrelation — ACF / PACF

Air quality measurements taken minutes apart are correlated — a high reading now predicts a high reading in 5 minutes. This autocorrelation violates the independence assumption of OLS and inflates false-positive rates in standard tests. The ACF (autocorrelation function) and PACF (partial ACF) quantify this.

for (s in levels(df$sensor)) {
  cat(sprintf("\n\n#### Sensor: %s\n\n", s))
  pm_vec <- df |> filter(sensor == s) |>
    arrange(datetime) |> pull(pm25_ug_m3)
  par(mfrow = c(1, 2))
  acf(pm_vec,  lag.max = 48,
      main = paste("ACF — PM2.5 |", s),
      col  = "#0072B2", lwd = 1.5)
  pacf(pm_vec, lag.max = 48,
       main = paste("PACF — PM2.5 |", s),
       col  = "#D55E00", lwd = 1.5)
  par(mfrow = c(1, 1))
}

8.5 Mixed Effects Models

# Run once in the Console if not yet installed
install.packages(c("lme4", "nlme", "lmerTest"))

library(lme4)
library(lmerTest)   # adds p-values to lmer() via Satterthwaite df
library(nlme)       # needed for AR(1) correlation structure

m_ri <- lmer(
  pm25_ug_m3 ~ temperature_C + humidity_pct + pressure_hPa +
               light_lux + oxidised_kO + reduced_kO + nh3_kO +
               (1 | sensor),
  data   = df,
  REML   = TRUE
)

summary(m_ri)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## pm25_ug_m3 ~ temperature_C + humidity_pct + pressure_hPa + light_lux +  
##     oxidised_kO + reduced_kO + nh3_kO + (1 | sensor)
##    Data: df
## 
## REML criterion at convergence: 98234.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -9.5190 -0.5043 -0.1934  0.2916 15.6956 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  sensor   (Intercept)  9.276   3.046   
##  Residual             32.046   5.661   
## Number of obs: 15559, groups:  sensor, 10
## 
## Fixed effects:
##                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    2.614e+01  3.556e+00  1.508e+03   7.352 3.20e-13 ***
## temperature_C -1.286e-01  1.748e-02  1.555e+04  -7.356 1.99e-13 ***
## humidity_pct  -4.880e-02  1.179e-02  1.555e+04  -4.138 3.52e-05 ***
## pressure_hPa  -2.119e-02  3.202e-03  1.555e+04  -6.618 3.75e-11 ***
## light_lux     -7.524e-06  6.132e-05  1.555e+04  -0.123   0.9023    
## oxidised_kO   -2.522e-03  1.773e-04  1.555e+04 -14.229  < 2e-16 ***
## reduced_kO     2.233e-03  8.692e-04  1.549e+04   2.569   0.0102 *  
## nh3_kO         1.837e-02  6.923e-04  1.550e+04  26.527  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmpr_C hmdty_ prss_P lght_l oxds_O rdcd_O
## temperatr_C -0.483                                          
## humidty_pct -0.589  0.868                                   
## pressure_hP -0.937  0.304  0.443                            
## light_lux    0.003  0.000 -0.004 -0.005                     
## oxidised_kO -0.008 -0.175 -0.050  0.040  0.003              
## reduced_kO   0.264 -0.049 -0.253 -0.319  0.014 -0.089       
## nh3_kO      -0.254  0.388  0.480  0.167 -0.001 -0.112 -0.520

broom.mixed::tidy(m_ri, effects = "fixed", conf.int = TRUE) |>
  mutate(sig = case_when(p.value < 0.001 ~ "***",
                         p.value < 0.01  ~ "**",
                         p.value < 0.05  ~ "*",
                         TRUE            ~ "")) |>
  kable(digits = 4, format = "html",
        caption = "Random-intercept model: fixed effects (sensor as random intercept)") |>
  kable_styling(full_width = FALSE)

# Intraclass correlation: proportion of variance due to between-sensor differences
vc   <- as.data.frame(VarCorr(m_ri))
icc  <- vc$vcov[vc$grp == "sensor"] /
        sum(vc$vcov)
cat(sprintf("ICC = %.3f — %.1f%% of PM2.5 variance is between-sensor\n",
            icc, icc * 100))

## ICC = 0.224 — 22.4% of PM2.5 variance is between-sensor

m_rs <- lmer(
  pm25_ug_m3 ~ temperature_C + humidity_pct + pressure_hPa +
               light_lux + oxidised_kO + reduced_kO + nh3_kO +
               (1 + humidity_pct | sensor),
  data    = df,
  REML    = TRUE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(m_rs)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## pm25_ug_m3 ~ temperature_C + humidity_pct + pressure_hPa + light_lux +  
##     oxidised_kO + reduced_kO + nh3_kO + (1 + humidity_pct | sensor)
##    Data: df
## Control: lmerControl(optimizer = "bobyqa")
## 
## REML criterion at convergence: 97224.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.2445 -0.4913 -0.1983  0.3147 16.2538 
## 
## Random effects:
##  Groups   Name         Variance Std.Dev. Corr 
##  sensor   (Intercept)   8.88640 2.9810        
##           humidity_pct  0.03657 0.1912   -0.70
##  Residual              29.96047 5.4736        
## Number of obs: 15559, groups:  sensor, 10
## 
## Fixed effects:
##                 Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)    7.964e+00  3.555e+00  1.512e+03   2.240 0.025217 *  
## temperature_C -9.806e-02  1.733e-02  1.546e+04  -5.660 1.54e-08 ***
## humidity_pct   4.176e-02  6.176e-02  9.563e+00   0.676 0.514902    
## pressure_hPa  -7.238e-03  3.199e-03  1.547e+04  -2.263 0.023651 *  
## light_lux      5.142e-07  5.931e-05  1.554e+04   0.009 0.993083    
## oxidised_kO   -2.201e-03  1.733e-04  1.551e+04 -12.703  < 2e-16 ***
## reduced_kO     3.301e-03  8.556e-04  1.485e+04   3.858 0.000115 ***
## nh3_kO         2.275e-02  7.095e-04  1.424e+04  32.069  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmpr_C hmdty_ prss_P lght_l oxds_O rdcd_O
## temperatr_C -0.465                                          
## humidty_pct -0.306  0.163                                   
## pressure_hP -0.939  0.286  0.095                            
## light_lux    0.004 -0.002 -0.002 -0.006                     
## oxidised_kO -0.011 -0.178 -0.008  0.042  0.005              
## reduced_kO   0.238 -0.025 -0.039 -0.297  0.011 -0.101       
## nh3_kO      -0.239  0.362  0.090  0.150  0.002 -0.069 -0.491

# Compare random-intercept vs random-slope via likelihood ratio test
# (refit with ML for comparison)
m_ri_ml <- update(m_ri, REML = FALSE)
m_rs_ml <- update(m_rs, REML = FALSE)

anova(m_ri_ml, m_rs_ml) |>
  as.data.frame() |>
  kable(digits = 4, format = "html",
        caption = "LRT: random intercept vs random intercept + slope") |>
  kable_styling(full_width = FALSE)

# nlme requires data sorted by group and time
df_nlme <- df |>
  arrange(sensor, datetime) |>
  mutate(time_num = as.numeric(datetime))   # seconds since epoch

m_ar1 <- lme(
  pm25_ug_m3 ~ temperature_C + humidity_pct + pressure_hPa +
               light_lux + oxidised_kO + reduced_kO + nh3_kO,
  random     = ~ 1 | sensor,
  correlation = corCAR1(form = ~ time_num | sensor),  # continuous AR(1)
  data       = df_nlme,
  method     = "ML"
)

summary(m_ar1)

## Linear mixed-effects model fit by maximum likelihood
##   Data: df_nlme 
##        AIC      BIC    logLik
##   98174.56 98258.74 -49076.28
## 
## Random effects:
##  Formula: ~1 | sensor
##         (Intercept) Residual
## StdDev:    2.887226 5.659679
## 
## Correlation Structure: Continuous AR(1)
##  Formula: ~time_num | sensor 
##  Parameter estimate(s):
## Phi 
## 0.2 
## Fixed effects:  pm25_ug_m3 ~ temperature_C + humidity_pct + pressure_hPa + light_lux +      oxidised_kO + reduced_kO + nh3_kO 
##                   Value Std.Error    DF    t-value p-value
## (Intercept)   26.133132  3.542421 15542   7.377196  0.0000
## temperature_C -0.128587  0.017476 15542  -7.357976  0.0000
## humidity_pct  -0.048800  0.011792 15542  -4.138326  0.0000
## pressure_hPa  -0.021186  0.003202 15542  -6.615340  0.0000
## light_lux     -0.000008  0.000061 15542  -0.125470  0.9002
## oxidised_kO   -0.002522  0.000177 15542 -14.228192  0.0000
## reduced_kO     0.002227  0.000869 15542   2.561746  0.0104
## nh3_kO         0.018366  0.000692 15542  26.528493  0.0000
##  Correlation: 
##               (Intr) tmpr_C hmdty_ prss_P lght_l oxds_O rdcd_O
## temperature_C -0.485                                          
## humidity_pct  -0.591  0.868                                   
## pressure_hPa  -0.941  0.304  0.443                            
## light_lux      0.003  0.000 -0.004 -0.005                     
## oxidised_kO   -0.008 -0.175 -0.050  0.040  0.003              
## reduced_kO     0.265 -0.049 -0.253 -0.319  0.014 -0.089       
## nh3_kO        -0.255  0.388  0.480  0.167 -0.001 -0.112 -0.520
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -9.5209525 -0.5044998 -0.1934267  0.2915889 15.6991402 
## 
## Number of Observations: 15559
## Number of Groups: 10

# Extract the estimated AR(1) parameter (phi)
phi <- coef(m_ar1$modelStruct$corStruct, unconstrained = FALSE)
cat(sprintf("Estimated AR(1) phi = %.3f\n", phi))

## Estimated AR(1) phi = 0.200

cat("(phi near 1 = strong persistence; near 0 = little autocorrelation)\n")

## (phi near 1 = strong persistence; near 0 = little autocorrelation)

intervals(m_ar1, which = "fixed")$fixed |>
  as.data.frame() |>
  rownames_to_column("term") |>
  kable(digits = 4, format = "html",
        caption = "AR(1) mixed model: fixed effects with 95% CI") |>
  kable_styling(full_width = FALSE)

# AIC comparison across all three models
# (refit m_ri and m_rs with ML to compare fairly with m_ar1)
tibble(
  Model       = c("OLS (per sensor, no RE)",
                  "Random intercept (lmer)",
                  "Random intercept + slope (lmer)",
                  "Random intercept + AR(1) (nlme)"),
  AIC         = c(
    mean(map_dbl(lm_list, ~ AIC(.x$model))),   # mean per-sensor OLS AIC
    AIC(m_ri_ml),
    AIC(m_rs_ml),
    AIC(m_ar1)
  ),
  Notes = c(
    "Ignores pseudoreplication and autocorrelation",
    "Accounts for between-sensor variance",
    "Also allows sensor-specific humidity slope",
    "Accounts for temporal autocorrelation"
  )
) |>
  kable(digits = 1, format = "html", caption = "AIC comparison across model families") |>
  kable_styling(full_width = FALSE)

# Check whether the random-intercept model reduced autocorrelation
# relative to OLS (compare ACF of OLS residuals vs lmer residuals)
ols_resid  <- resid(lm_list[[1]]$model)
lmer_resid <- resid(m_ri)[ df$sensor == levels(df$sensor)[1] ]

par(mfrow = c(1, 2))
acf(ols_resid,  lag.max = 48, main = "ACF: OLS residuals (sensor 1)",
    col = "#D55E00", lwd = 1.5)
acf(lmer_resid, lag.max = 48, main = "ACF: lmer residuals (sensor 1)",
    col = "#0072B2", lwd = 1.5)

par(mfrow = c(1, 1))

8.6 Strongest Relationships: PM2.5 vs Top Predictors

Using the regression results, we plot the two strongest predictors (by |t-statistic|) for each sensor.

top_preds <- map_dfr(lm_list, "tidy") |>
  filter(term != "(Intercept)") |>
  group_by(sensor) |>
  slice_max(abs(statistic), n = 2) |>
  ungroup()

# One figure per sensor — two side-by-side scatter plots
for (s in levels(df$sensor)) {
  cat(sprintf("\n\n#### Sensor: %s\n\n", s))

  d_top <- top_preds |> filter(sensor == s)
  preds <- d_top$term

  plots <- map(preds, function(pred) {
    df |>
      filter(sensor == s) |>
      ggplot(aes(x = .data[[pred]], y = pm25_ug_m3)) +
      geom_point(alpha = 0.08, size = 0.8, color = pal[[s]]) +
      geom_smooth(method = "lm", color = "black",
                  linewidth = 0.8, se = TRUE) +
      labs(title    = pred,
           subtitle = sprintf("t = %.2f, p = %.4f",
                              d_top$statistic[d_top$term == pred],
                              d_top$p.value[d_top$term == pred]),
           x = pred, y = "PM2.5 (µg/m³)") +
      theme(plot.title    = element_text(size = 11),
            plot.subtitle = element_text(size = 9, color = "grey40"))
  })

  print(wrap_plots(plots, ncol = 2))
}