1.1 Overall summary

log_raw |>
  summarise(
    N          = n(),
    Mean_temp  = round(mean(mean_air_temperature_c), 2),
    SD_temp    = round(sd(mean_air_temperature_c),   2),
    Median_temp = round(median(mean_air_temperature_c), 2),
    IQR_temp   = round(IQR(mean_air_temperature_c),  2),
    Min_temp   = round(min(mean_air_temperature_c),  2),
    Max_temp   = round(max(mean_air_temperature_c),  2),
    Skewness   = round(moments::skewness(mean_air_temperature_c), 3),
    Mean_canopy = round(mean(canopy_cover_percent), 1),
    SD_canopy   = round(sd(canopy_cover_percent),   1)
  ) |>
  kable(caption = "Overall descriptive statistics — logger air temperature and canopy cover") |>
  kable_styling(full_width = FALSE)

1.2 By urban setting and income group

log_raw |>
  group_by(urban_setting, demographic_group) |>
  summarise(
    N    = n(),
    Mean = round(mean(mean_air_temperature_c), 2),
    SD   = round(sd(mean_air_temperature_c),   2),
    Med  = round(median(mean_air_temperature_c), 2),
    Canopy_mean = round(mean(canopy_cover_percent), 1),
    .groups = "drop"
  ) |>
  rename(`Urban setting` = urban_setting,
         `Income group`  = demographic_group,
         `Mean temp (°C)` = Mean, `SD` = SD, `Median` = Med,
         `Mean canopy (%)` = Canopy_mean) |>
  kable(caption = "Mean 24-hour air temperature (°C) by urban setting and income group") |>
  kable_styling(full_width = TRUE) |>
  column_spec(1, bold = TRUE)

1.3 Normality check

set.seed(42)
norm_tbl <- log_raw |>
  group_by(urban_setting) |>
  summarise(
    n  = n(),
    W  = round(shapiro.test(mean_air_temperature_c)$statistic, 4),
    p  = round(shapiro.test(mean_air_temperature_c)$p.value,   4),
    Skew = round(moments::skewness(mean_air_temperature_c), 3),
    .groups = "drop"
  ) |>
  mutate(Normal = ifelse(p > 0.05, "Yes", "No"))

kable(norm_tbl,
      caption = "Shapiro-Wilk normality test per urban setting (logger data)",
      col.names = c("Setting","N","W","p-value","Skewness","Normal?")) |>
  kable_styling(full_width = FALSE)

2.1 Temperature distributions by urban setting

ggplot(log_raw, aes(x = urban_setting, y = mean_air_temperature_c,
                    fill = urban_setting)) +
  geom_boxplot(outlier.alpha = 0.4, outlier.size = 1.5, width = 0.55) +
  geom_jitter(aes(color = demographic_group), width = 0.12,
              alpha = 0.6, size = 2) +
  stat_summary(fun = mean, geom = "point",
               shape = 23, size = 3, fill = "white", color = "#1e3a2f") +
  scale_fill_manual(values  = setting_pal, guide = "none") +
  scale_color_manual(values = income_pal, name = "Income group") +
  labs(title    = "24-Hour Mean Air Temperature by Urban Setting",
       subtitle = "Points colored by income group · ◇ = group mean",
       x = NULL, y = "Mean air temperature (°C)") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

2.2 Temperature by income group across settings

ggplot(log_raw, aes(x = demographic_group, y = mean_air_temperature_c,
                    fill = demographic_group)) +
  geom_boxplot(outlier.alpha = 0.35, width = 0.5) +
  stat_summary(fun = mean, geom = "point",
               shape = 23, size = 2.5, fill = "white") +
  facet_wrap(~ urban_setting, nrow = 1) +
  scale_fill_manual(values = income_pal, guide = "none") +
  labs(title    = "Air Temperature by Income Group within Each Urban Setting",
       subtitle = "◇ = mean",
       x = NULL, y = "Mean air temperature (°C)") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 9))

2.3 Canopy cover vs air temperature

ggplot(log_raw, aes(x = canopy_cover_percent, y = mean_air_temperature_c,
                    color = demographic_group, shape = urban_setting)) +
  geom_point(alpha = 0.75, size = 2.5) +
  geom_smooth(aes(group = 1), method = "lm", se = TRUE,
              color = "#1e3a2f", linewidth = 1, linetype = "dashed") +
  scale_color_manual(values = income_pal, name = "Income group") +
  scale_shape_manual(values = c(16, 17, 15, 18, 8), name = "Urban setting") +
  labs(title    = "Canopy Cover vs Mean Air Temperature (24-hour loggers)",
       subtitle = "Dashed line = overall OLS trend",
       x = "Canopy cover (%)", y = "Mean air temperature (°C)") +
  guides(color = guide_legend(nrow = 1),
         shape = guide_legend(nrow = 1))

2.4 Canopy cover distribution by income group

ggplot(log_raw, aes(x = canopy_cover_percent, fill = demographic_group,
                    color = demographic_group)) +
  geom_density(alpha = 0.30, linewidth = 0.8) +
  scale_fill_manual(values  = income_pal, name = "Income group") +
  scale_color_manual(values = income_pal, name = "Income group") +
  labs(title    = "Canopy Cover Distribution by Income Group — Logger Sites",
       subtitle = "Are low-income sites systematically under-canopied?",
       x = "Canopy cover (%)", y = "Density")

3.1 One-way ANOVA: does temperature differ by urban setting?

H₀: Mean 24-hour air temperature is equal across all five urban settings. H₁: At least one setting has a different mean temperature.

aov_setting <- aov(mean_air_temperature_c ~ urban_setting, data = log_raw)
summary(aov_setting)

##                Df Sum Sq Mean Sq F value Pr(>F)
## urban_setting   4    2.5  0.6173   0.247  0.911
## Residuals     175  436.6  2.4951

tidy(aov_setting) |>
  mutate(across(where(is.numeric), ~ round(.x, 4))) |>
  kable(caption = "ANOVA table: air temperature ~ urban setting") |>
  kable_styling(full_width = FALSE)

3.2 Post-hoc comparisons (Tukey HSD)

em_setting <- emmeans(aov_setting, ~ urban_setting)
pairs_setting <- pairs(em_setting, adjust = "tukey")
summary(pairs_setting) |>
  as.data.frame() |>
  mutate(across(where(is.numeric), ~ round(.x, 4)),
         sig = case_when(p.value < 0.001 ~ "***",
                         p.value < 0.01  ~ "**",
                         p.value < 0.05  ~ "*",
                         TRUE            ~ "")) |>
  kable(caption = "Tukey HSD pairwise comparisons — urban setting") |>
  kable_styling(full_width = FALSE)

as.data.frame(em_setting) |>
  ggplot(aes(x = reorder(urban_setting, emmean), y = emmean,
             ymin = lower.CL, ymax = upper.CL, color = urban_setting)) +
  geom_pointrange(size = 0.9, linewidth = 1) +
  scale_color_manual(values = setting_pal, guide = "none") +
  coord_flip() +
  labs(title    = "Estimated Marginal Means ± 95% CI",
       subtitle = "Urban setting effect on mean air temperature",
       x = NULL, y = "Estimated mean temperature (°C)")

3.3 Two-way ANOVA: setting × income group interaction

aov_2way <- aov(mean_air_temperature_c ~ urban_setting * demographic_group,
                data = log_raw)
tidy(aov_2way) |>
  mutate(across(where(is.numeric), ~ round(.x, 4)),
         sig = case_when(p.value < 0.001 ~ "***",
                         p.value < 0.01  ~ "**",
                         p.value < 0.05  ~ "*",
                         TRUE            ~ "")) |>
  kable(caption = "Two-way ANOVA: temperature ~ setting * income group") |>
  kable_styling(full_width = FALSE)

3.4 ANOVA assumptions

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
plot(aov_2way)

par(mfrow = c(1, 1))

leveneTest(mean_air_temperature_c ~ urban_setting * demographic_group,
           data = log_raw) |>
  tidy() |>
  mutate(across(where(is.numeric), ~ round(.x, 4))) |>
  kable(caption = "Levene's test for homogeneity of variance") |>
  kable_styling(full_width = FALSE)

3.5 Simple linear regression: canopy cover predicts temperature

lm_canopy <- lm(mean_air_temperature_c ~ canopy_cover_percent, data = log_raw)

tidy(lm_canopy) |>
  mutate(across(where(is.numeric), ~ round(.x, 4))) |>
  kable(caption = "OLS regression: air temperature ~ canopy cover (%)") |>
  kable_styling(full_width = FALSE)

glance(lm_canopy) |>
  select(r.squared, adj.r.squared, statistic, p.value, nobs) |>
  mutate(across(where(is.numeric), ~ round(.x, 4))) |>
  kable(caption = "Model fit statistics") |>
  kable_styling(full_width = FALSE)

ggplot(log_raw, aes(x = canopy_cover_percent, y = mean_air_temperature_c)) +
  geom_point(aes(color = urban_setting, shape = demographic_group),
             alpha = 0.7, size = 2.5) +
  geom_smooth(method = "lm", se = TRUE,
              color = "#1e3a2f", linewidth = 1) +
  scale_color_manual(values = setting_pal, name = "Setting") +
  scale_shape_manual(values = c(16, 17, 8), name = "Income") +
  annotate("text",
           x = max(log_raw$canopy_cover_percent) * 0.85,
           y = max(log_raw$mean_air_temperature_c) * 0.99,
           label = paste0("R² = ", round(glance(lm_canopy)$r.squared, 3)),
           color = "#1e3a2f", fontface = "bold", size = 4) +
  labs(title    = "OLS Regression: Canopy Cover Predicts Air Temperature",
       x = "Canopy cover (%)", y = "Mean air temperature (°C)")

3.6 Multiple regression: canopy + setting + income group

lm_multi <- lm(mean_air_temperature_c ~
                 canopy_cover_percent + urban_setting + demographic_group,
               data = log_raw)

tidy(lm_multi) |>
  mutate(
    sig = case_when(p.value < 0.001 ~ "***", p.value < 0.01 ~ "**",
                    p.value < 0.05  ~ "*",   TRUE ~ ""),
    across(where(is.numeric), ~ round(.x, 4))
  ) |>
  kable(caption = "Multiple regression: temperature ~ canopy + setting + income") |>
  kable_styling(full_width = FALSE)

glance(lm_multi) |>
  select(r.squared, adj.r.squared, AIC, BIC, nobs) |>
  mutate(across(where(is.numeric), ~ round(.x, 3))) |>
  kable(caption = "Multiple regression model fit") |>
  kable_styling(full_width = FALSE)

3.7 Non-parametric check: Kruskal-Wallis test

kw_setting <- kruskal.test(mean_air_temperature_c ~ urban_setting, data = log_raw)
kw_income  <- kruskal.test(mean_air_temperature_c ~ demographic_group, data = log_raw)

tibble(
  Test      = c("Urban setting", "Income group"),
  Chi_sq    = round(c(kw_setting$statistic, kw_income$statistic), 3),
  df        = c(kw_setting$parameter, kw_income$parameter),
  p_value   = round(c(kw_setting$p.value, kw_income$p.value), 4),
  Sig       = ifelse(c(kw_setting$p.value, kw_income$p.value) < 0.05, "✓", "—")
) |>
  kable(caption = "Kruskal-Wallis tests (non-parametric alternative to ANOVA)") |>
  kable_styling(full_width = FALSE)

4.1 Overall summary

ir_raw |>
  summarise(
    N            = n(),
    Mean_surf    = round(mean(ir_surface_temperature_c),   2),
    SD_surf      = round(sd(ir_surface_temperature_c),     2),
    Median_surf  = round(median(ir_surface_temperature_c), 2),
    IQR_surf     = round(IQR(ir_surface_temperature_c),    2),
    Min_surf     = round(min(ir_surface_temperature_c),    2),
    Max_surf     = round(max(ir_surface_temperature_c),    2),
    Skewness     = round(moments::skewness(ir_surface_temperature_c), 3),
    Mean_canopy  = round(mean(canopy_cover_percent), 1),
    Hour_range   = paste(min(measurement_hour), "–", max(measurement_hour))
  ) |>
  kable(caption = "Overall descriptive statistics — IR surface temperature") |>
  kable_styling(full_width = FALSE)

4.2 By urban setting and income group

ir_raw |>
  group_by(urban_setting, demographic_group) |>
  summarise(
    N    = n(),
    Mean = round(mean(ir_surface_temperature_c),   2),
    SD   = round(sd(ir_surface_temperature_c),     2),
    Med  = round(median(ir_surface_temperature_c), 2),
    Max  = round(max(ir_surface_temperature_c),    2),
    .groups = "drop"
  ) |>
  rename(`Urban setting` = urban_setting, `Income group` = demographic_group,
         `Mean (°C)` = Mean, `Max (°C)` = Max) |>
  kable(caption = "IR surface temperature summary by urban setting and income group") |>
  kable_styling(full_width = TRUE) |>
  column_spec(1, bold = TRUE)

4.3 By time of day

ir_raw |>
  group_by(time_of_day, urban_setting) |>
  summarise(
    N    = n(),
    Mean = round(mean(ir_surface_temperature_c), 2),
    SD   = round(sd(ir_surface_temperature_c),   2),
    .groups = "drop"
  ) |>
  pivot_wider(names_from = time_of_day,
              values_from = c(Mean, SD),
              names_glue  = "{time_of_day}_{.value}") |>
  kable(caption = "Mean IR surface temperature (°C) by time of day and urban setting") |>
  kable_styling(full_width = TRUE)

4.4 Normality check

ir_raw |>
  group_by(urban_setting) |>
  summarise(
    n    = n(),
    W    = round(shapiro.test(sample(ir_surface_temperature_c,
                                     min(n(), 5000)))$statistic, 4),
    p    = round(shapiro.test(sample(ir_surface_temperature_c,
                                     min(n(), 5000)))$p.value, 4),
    Skew = round(moments::skewness(ir_surface_temperature_c), 3),
    .groups = "drop"
  ) |>
  mutate(Normal = ifelse(p > 0.05, "Yes", "No")) |>
  kable(caption = "Shapiro-Wilk normality test per urban setting (IR data)") |>
  kable_styling(full_width = FALSE)

5.1 Surface temperature distributions

p_box2 <- ggplot(ir_raw,
                 aes(x = urban_setting, y = ir_surface_temperature_c,
                     fill = urban_setting)) +
  geom_boxplot(outlier.alpha = 0.2, outlier.size = 0.8, width = 0.55) +
  geom_jitter(aes(color = demographic_group), width = 0.1,
              alpha = 0.25, size = 1) +
  stat_summary(fun = mean, geom = "point",
               shape = 23, size = 2.5, fill = "white", color = "#1e3a2f") +
  scale_fill_manual(values  = setting_pal, guide = "none") +
  scale_color_manual(values = income_pal,  name  = "Income group") +
  labs(title = "IR Surface Temperature by Urban Setting",
       subtitle = "◇ = mean; points colored by income group",
       x = NULL, y = "Surface temperature (°C)") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

p_dens2 <- ggplot(ir_raw,
                  aes(x = ir_surface_temperature_c, fill = urban_setting,
                      color = urban_setting)) +
  geom_density(alpha = 0.25, linewidth = 0.7) +
  geom_vline(xintercept = 35, linetype = "dashed",
             color = "firebrick", linewidth = 0.8) +
  scale_fill_manual(values  = setting_pal, name = "Setting") +
  scale_color_manual(values = setting_pal, name = "Setting") +
  labs(title    = "Density of Surface Temperatures",
       subtitle = "Dashed = 35°C heat stress threshold",
       x = "Surface temperature (°C)", y = "Density")

p_box2 + p_dens2

5.2 Diurnal pattern by urban setting

ir_raw |>
  group_by(urban_setting, measurement_hour) |>
  summarise(
    mean_temp = mean(ir_surface_temperature_c),
    se_temp   = sd(ir_surface_temperature_c) / sqrt(n()),
    .groups   = "drop"
  ) |>
  ggplot(aes(x = measurement_hour, y = mean_temp,
             color = urban_setting, fill = urban_setting)) +
  geom_ribbon(aes(ymin = mean_temp - se_temp, ymax = mean_temp + se_temp),
              alpha = 0.15, color = NA) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 2.5) +
  scale_color_manual(values = setting_pal, name = "Urban setting") +
  scale_fill_manual(values  = setting_pal, guide = "none") +
  scale_x_continuous(breaks = 9:17,
                     labels = paste0(9:17, ":00")) +
  labs(title    = "Diurnal Profile of IR Surface Temperature",
       subtitle = "Mean ± 1 SE by hour of day",
       x = "Hour of day", y = "Surface temperature (°C)") +
  theme(axis.text.x = element_text(angle = 40, hjust = 1))

5.3 Diurnal pattern by income group

ir_raw |>
  group_by(demographic_group, measurement_hour) |>
  summarise(
    mean_temp = mean(ir_surface_temperature_c),
    se_temp   = sd(ir_surface_temperature_c) / sqrt(n()),
    .groups   = "drop"
  ) |>
  ggplot(aes(x = measurement_hour, y = mean_temp,
             color = demographic_group, fill = demographic_group)) +
  geom_ribbon(aes(ymin = mean_temp - se_temp, ymax = mean_temp + se_temp),
              alpha = 0.18, color = NA) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 2.5) +
  scale_color_manual(values = income_pal, name = "Income group") +
  scale_fill_manual(values  = income_pal, guide = "none") +
  scale_x_continuous(breaks = 9:17, labels = paste0(9:17, ":00")) +
  labs(title    = "Diurnal Profile by Income Group",
       subtitle = "Do low-income sites stay hotter throughout the day?",
       x = "Hour of day", y = "Surface temperature (°C)") +
  theme(axis.text.x = element_text(angle = 40, hjust = 1))

5.4 Canopy cover vs surface temperature by time of day

ggplot(ir_raw,
       aes(x = canopy_cover_percent, y = ir_surface_temperature_c,
           color = demographic_group)) +
  geom_point(alpha = 0.4, size = 1.5) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 1.1) +
  facet_wrap(~ time_of_day) +
  scale_color_manual(values = income_pal, name = "Income group") +
  labs(title    = "Canopy Cover vs Surface Temperature by Time of Day",
       subtitle = "Separate regression lines per income group",
       x = "Canopy cover (%)", y = "Surface temperature (°C)") +
  guides(color = guide_legend(nrow = 1))

5.5 Income group × urban setting heatmap

ir_raw |>
  group_by(urban_setting, demographic_group) |>
  summarise(mean_temp = mean(ir_surface_temperature_c), .groups = "drop") |>
  ggplot(aes(x = demographic_group, y = fct_rev(urban_setting),
             fill = mean_temp)) +
  geom_tile(color = "white", linewidth = 0.8) +
  geom_text(aes(label = round(mean_temp, 1)),
            color = "white", fontface = "bold", size = 4) +
  scale_fill_gradient2(low = "#2e86ab", mid = "#f6ae2d", high = "#f26419",
                       midpoint = 42, name = "°C") +
  labs(title = "Mean IR Surface Temperature (°C)",
       subtitle = "Urban setting × income group",
       x = NULL, y = NULL) +
  theme(legend.position = "right")

6.1 One-way ANOVA: surface temperature by urban setting

aov2_setting <- aov(ir_surface_temperature_c ~ urban_setting, data = ir_raw)
tidy(aov2_setting) |>
  mutate(across(where(is.numeric), ~ round(.x, 4)),
         sig = case_when(p.value < 0.001 ~ "***", p.value < 0.01 ~ "**",
                         p.value < 0.05 ~ "*", TRUE ~ "")) |>
  kable(caption = "ANOVA: IR surface temperature ~ urban setting") |>
  kable_styling(full_width = FALSE)

6.2 One-way ANOVA: surface temperature by income group

aov2_income <- aov(ir_surface_temperature_c ~ demographic_group, data = ir_raw)
tidy(aov2_income) |>
  mutate(across(where(is.numeric), ~ round(.x, 4)),
         sig = case_when(p.value < 0.001 ~ "***", p.value < 0.01 ~ "**",
                         p.value < 0.05 ~ "*", TRUE ~ "")) |>
  kable(caption = "ANOVA: IR surface temperature ~ income group") |>
  kable_styling(full_width = FALSE)

6.3 Tukey HSD post-hoc: urban setting

em2 <- emmeans(aov2_setting, ~ urban_setting)
pairs(em2, adjust = "tukey") |>
  as.data.frame() |>
  mutate(across(where(is.numeric), ~ round(.x, 4)),
         sig = case_when(p.value < 0.001 ~ "***", p.value < 0.01 ~ "**",
                         p.value < 0.05 ~ "*", TRUE ~ "")) |>
  kable(caption = "Tukey HSD: urban setting pairwise comparisons (IR data)") |>
  kable_styling(full_width = FALSE)

6.4 Kruskal-Wallis: time of day effect

kw_tod <- kruskal.test(ir_surface_temperature_c ~ time_of_day, data = ir_raw)

tibble(
  Test    = "Time of day (Morning / Midday / Afternoon)",
  Chi_sq  = round(kw_tod$statistic, 3),
  df      = kw_tod$parameter,
  p_value = round(kw_tod$p.value,   4),
  Sig     = ifelse(kw_tod$p.value < 0.05, "✓", "—")
) |>
  kable(caption = "Kruskal-Wallis: does surface temperature vary by time of day?") |>
  kable_styling(full_width = FALSE)

6.5 Simple regression: canopy cover predicts surface temperature

lm2_canopy <- lm(ir_surface_temperature_c ~ canopy_cover_percent, data = ir_raw)

tidy(lm2_canopy) |>
  mutate(across(where(is.numeric), ~ round(.x, 4))) |>
  kable(caption = "OLS regression: surface temperature ~ canopy cover (%)") |>
  kable_styling(full_width = FALSE)

glance(lm2_canopy) |>
  select(r.squared, adj.r.squared, statistic, p.value, nobs) |>
  mutate(across(where(is.numeric), ~ round(.x, 4))) |>
  kable(caption = "Model fit (simple regression)") |>
  kable_styling(full_width = FALSE)

6.6 Multiple regression: canopy + hour + setting + income

lm2_multi <- lm(ir_surface_temperature_c ~
                  canopy_cover_percent +
                  measurement_hour     +
                  urban_setting        +
                  demographic_group,
                data = ir_raw)

tidy(lm2_multi) |>
  mutate(
    sig = case_when(p.value < 0.001 ~ "***", p.value < 0.01 ~ "**",
                    p.value < 0.05  ~ "*",   TRUE ~ ""),
    across(where(is.numeric), ~ round(.x, 4))
  ) |>
  kable(caption = "Multiple regression: surface temp ~ canopy + hour + setting + income") |>
  kable_styling(full_width = FALSE)

glance(lm2_multi) |>
  select(r.squared, adj.r.squared, AIC, BIC, nobs) |>
  mutate(across(where(is.numeric), ~ round(.x, 3))) |>
  kable(caption = "Multiple regression model fit (IR data)") |>
  kable_styling(full_width = FALSE)

6.7 Model comparison: does adding hour improve fit?

lm2_base   <- lm(ir_surface_temperature_c ~
                   canopy_cover_percent + urban_setting + demographic_group,
                 data = ir_raw)
lm2_withhour <- lm2_multi   # already includes measurement_hour

anova(lm2_base, lm2_withhour) |>
  tidy() |>
  mutate(across(where(is.numeric), ~ round(.x, 4))) |>
  kable(caption = "Likelihood ratio test: does measurement hour improve model fit?") |>
  kable_styling(full_width = FALSE)

tibble(
  Model   = c("Base (canopy + setting + income)",
              "Full (+ measurement hour)"),
  AIC     = round(c(AIC(lm2_base), AIC(lm2_multi)), 2),
  adj_R2  = round(c(summary(lm2_base)$adj.r.squared,
                    summary(lm2_multi)$adj.r.squared), 4)
) |>
  kable(caption = "AIC and adjusted R² comparison") |>
  kable_styling(full_width = FALSE)

6.8 Interaction: does canopy effect vary by income group?

lm2_interact <- lm(ir_surface_temperature_c ~
                     canopy_cover_percent * demographic_group +
                     measurement_hour + urban_setting,
                   data = ir_raw)

tidy(lm2_interact) |>
  filter(str_detect(term, "canopy|demographic|:")) |>
  mutate(
    sig = case_when(p.value < 0.001 ~ "***", p.value < 0.01 ~ "**",
                    p.value < 0.05  ~ "*",   TRUE ~ ""),
    across(where(is.numeric), ~ round(.x, 4))
  ) |>
  kable(caption = "Interaction model: canopy × income group coefficients (selected terms)") |>
  kable_styling(full_width = FALSE)

newdata <- expand.grid(
  canopy_cover_percent = seq(0, 95, by = 5),
  demographic_group    = income_lvl,
  measurement_hour     = 13,          # midday
  urban_setting        = "Park"       # reference setting
) |>
  mutate(demographic_group = factor(demographic_group, levels = income_lvl),
         urban_setting     = factor(urban_setting,     levels = setting_lvl))

newdata$pred_temp <- predict(lm2_interact, newdata = newdata)

ggplot(newdata, aes(x = canopy_cover_percent, y = pred_temp,
                   color = demographic_group)) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(values = income_pal, name = "Income group") +
  labs(title    = "Predicted Surface Temperature by Canopy Cover and Income Group",
       subtitle = "From interaction model · Midday · Park setting · 95% CI not shown",
       x = "Canopy cover (%)", y = "Predicted surface temperature (°C)")

6.9 Residual diagnostics

par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
plot(lm2_multi, main = "IR Multiple Regression Diagnostics")

par(mfrow = c(1, 1))