Why public spaces?

Southern California cities face an intensifying urban heat island effect. Paved surfaces, low tree canopy, and building geometry trap solar radiation and re-radiate it as heat, raising ambient temperatures by 2–5 °C relative to surrounding rural areas. Public spaces — plazas, bus stops, schoolyards, parks — are where residents spend unshaded time, often during the hottest parts of the day.

Understanding how much shade structures, trees, or water features reduce surface-level air temperature is actionable information: it connects directly to park design, tree planting policy, and climate adaptation planning. This activity gives you the analytical skills to answer that question with data you collected yourself.

The study design

Each team placed many sensors at a single site: some sensors in a shaded location (under a tree canopy, a built shade structure, or a building overhang) and some in a nearby exposed, sun-facing location. Sensors were left in place for at least 48 hours, recording temperature every ~5 minutes. The research group have data from multiple sites, each with shade/sun observations.

This structure — repeated measurements nested within sites — is exactly the kind of hierarchical data that requires a mixed-effects model rather than a simple t-test. The key issue is pseudoreplication: the hundreds of readings from a single sensor are not independent data points. Site conditions (neighborhood, proximity to buildings, baseline vegetation) create correlation among readings from the same location. A mixed model handles this correctly by partitioning variance between the site level and the within-site measurement level.

Sensor data files

Each sensor produces a .csv file named after its hardware ID (e.g., 260191640D.csv). I found that it’s best to clean the data by removing the header data then save as CSV with same name but different folder. The cleaned file has no column headers — just four comma-separated columns:

4/10/2026,16:25:13,79.3,26.3
4/10/2026,16:30:13,78.0,25.6
...

Download the class sensor files using the links below. Place all .csv files in data_raw/ inside your project.

⬇ 260191640D.csv ⬇ 260191640E.csv

Note: As more teams contribute data, additional sensor files will appear in the A14 data folder on GitHub. Download all .csv files you find there into your data_raw/ folder.

Metadata file

The metadata file locations_treatments.csv maps each sensor ID to a site name and a treatment label. Download it and place it in data_meta/.

⬇ locations_treatments.csv

It looks like this:

Recommended folder structure

Before writing any R code, create this directory layout inside an RStudio Project. If you haven’t set up projects yet, see R03 — Project Structure.

Download the helper script

We provide a ready-to-run R script that handles all the file-reading and merging in one go. Download it, place it in scripts/, and read through it before running — every line is commented.

⬇ read_and_merge.R

Run it once from the RStudio console (with your project open):

source("scripts/read_and_merge.R")

It will write data_processed/sensor_data_merged.csv. The sections below walk through exactly what that script does, line by line.

Step 1 — Load packages

library(tidyverse)   # readr + dplyr + ggplot2 + lubridate all at once
library(lubridate)   # friendly date/time parsing

Step 2 — Write a function to read one sensor file

The sensor files have no header row. We supply column names manually and parse the date and time into a proper POSIXct datetime object.

read_sensor_file <- function(file_path) {

  # Extract the sensor ID from the file name (everything before ".csv")
  sensor_id <- tools::file_path_sans_ext(basename(file_path))

  read_csv(
    file_path,
    col_names = c("date_str", "time_str", "temp_F", "temp_C"),
    col_types = cols(
      date_str = col_character(),
      time_str = col_character(),
      temp_F   = col_double(),
      temp_C   = col_double()
    )
  ) |>
    mutate(
      sensor_id = sensor_id,
      # Combine date and time strings, then parse with lubridate
      datetime  = mdy_hms(paste(date_str, time_str), tz = "America/Los_Angeles")
    ) |>
    select(sensor_id, datetime, temp_F, temp_C)
}

Let’s unpack the key lines:

• tools::file_path_sans_ext(basename(file_path)) — strips the directory path and the .csv extension, leaving just 260191640D. This becomes the sensor_id column that links to the metadata.
• col_names = c(...) — because the file has no header, we name the columns ourselves.
• mdy_hms(paste(date_str, time_str)) — lubridate::mdy_hms() parses strings in month/day/year hour:minute:second format. Pasting the two strings first gives it "4/10/2026 16:25:13", which it handles correctly.
• tz = "America/Los_Angeles" — always specify a timezone when working with local sensor data. Daylight saving time shifts can introduce phantom hour gaps or duplicates if you leave it as UTC.

Step 3 — Read all sensor files at once

Instead of calling read_sensor_file() once per file, we use purrr::map_dfr() to apply it to every .csv in data_raw/ and stack the results into a single data frame.

sensor_files <- list.files(
  path    = "data_clean",
  pattern = "\\.csv$",      # only files ending in .csv
  full.names = TRUE          # return full relative paths
)

sensor_data <- map_dfr(sensor_files, read_sensor_file)

glimpse(sensor_data)

result <- data.frame()
for (f in sensor_files) {
  result <- rbind(result, read_sensor_file(f))
}

Step 4 — Load the metadata and merge

locations <- read_csv(
  "data_meta/locations_treatments.csv",
  col_types = cols(
    Sensor_id = col_character(),
    Site      = col_character(),
    Treatment = col_character()
  )
) |>
  rename(sensor_id = Sensor_id)   # match the column name in sensor_data

sensor_merged <- sensor_data |>
  left_join(locations, by = "sensor_id")

# Integrity check: are any sensors unmatched?
missing <- sensor_merged |>
  filter(is.na(Treatment)) |>
  distinct(sensor_id)

if (nrow(missing) > 0) {
  warning("These sensor IDs have no metadata row: ",
          paste(missing$sensor_id, collapse = ", "))
}

A left_join() keeps every row in sensor_data and attaches the matching Site and Treatment columns from locations. Rows without a match get NA — hence the integrity check.

Step 5 — Save the merged file

write_csv(sensor_merged, "data_processed/sensor_data_merged.csv")
message("Saved: data_processed/sensor_data_merged.csv")

From this point on, your analysis script loads the processed file rather than re-reading all the raw sensor files. This keeps analysis fast and reproducible.

Time series: diel temperature cycles

The most important first plot is the raw time series, coloured by treatment. You should see a characteristic diel cycle — rising temperatures in the morning, a mid-afternoon peak, then cooling through the evening.

ggplot(sim_list, aes(x = datetime, y = temp_C,
                     colour = treatment, group = sensor_id)) +
  geom_line(linewidth = 0.6, alpha = 0.85) +
  facet_wrap(~ site, ncol = 1, scales = "free_y") +
  scale_colour_manual(
    values = c(Sun = "#c8832a", Shade = "#2d5a3d"),
    name   = "Treatment"
  ) +
  labs(
    x = "Date / Time",
    y = "Air temperature (°C)",
    title = "Diel temperature cycles — public space sensors",
    subtitle = "Each panel = one site; one line per sensor"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    strip.background = element_rect(fill = "#1e3a2f", colour = NA),
    strip.text       = element_text(colour = "white", face = "bold", size = 10),
    legend.position  = "top",
    plot.title       = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
    plot.subtitle    = element_text(colour = "#7a7a72", size = 10)
  )

Simulated 4-day diel temperature records for shaded and sun-exposed sensors across three sites. Shaded sensors (green) run consistently cooler and show a dampened daily amplitude compared with sun-exposed sensors (amber).

Boxplot: treatment comparison

A boxplot strips away the time dimension and summarises the temperature distribution for each treatment across all readings — a useful complement to the time series.

ggplot(sim_list, aes(x = treatment, y = temp_C, fill = treatment)) +
  geom_boxplot(notch = TRUE, outlier.size = 0.8, alpha = 0.85) +
  facet_wrap(~ site, nrow = 1) +
  scale_fill_manual(values = c(Sun = "#e8a84a", Shade = "#4a7c5e"),
                    guide = "none") +
  labs(
    x = NULL, y = "Air temperature (°C)",
    title    = "Temperature distributions by treatment",
    subtitle = "Notches indicate 95% CI around the median"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    strip.background = element_rect(fill = "#1e3a2f", colour = NA),
    strip.text       = element_text(colour = "white", face = "bold"),
    plot.title       = element_text(family = "serif", face = "bold", colour = "#1e3a2f")
  )

Distribution of all temperature readings by treatment and site. The median, interquartile range, and outliers summarise the full 4-day record for each sensor group.

Summary statistics

Compute group means as a quick sanity check — these should roughly match what you see in the boxplots.

sim_list |>
  group_by(site, treatment) |>
  summarise(
    n           = n(),
    mean_temp   = round(mean(temp_C, na.rm = TRUE), 2),
    sd_temp     = round(sd(temp_C,   na.rm = TRUE), 2),
    min_temp    = round(min(temp_C,  na.rm = TRUE), 1),
    max_temp    = round(max(temp_C,  na.rm = TRUE), 1),
    .groups = "drop"
  ) |>
  knitr::kable(
    col.names = c("Site", "Treatment", "n", "Mean (°C)", "SD (°C)", "Min (°C)", "Max (°C)"),
    caption   = "Descriptive statistics for each sensor group."
  )

Why not a simple t-test or ANOVA?

You might be tempted to run t.test(temp_C ~ treatment, data = sensor_merged). The problem: that test assumes every row in sensor_merged is an independent observation. But readings from the same sensor 5 minutes apart are highly correlated — they share the same physical location, the same sensor hardware, and the same ambient conditions. Treating them as independent observations would give you thousands of pseudo-replicates and an absurdly small p-value even if the treatment effect were negligible.

The mixed-effects model solves this by explicitly modelling the site-level variance as a random effect. It estimates how much temperatures vary across sites (random intercepts) and how consistent the shade effect is (fixed effect), while accounting for the non-independence of within-sensor readings.

Model formulation

We fit the following model using lme4::lmer():

\[\text{temp}_{ijk} = \underbrace{\beta_0 + \beta_1 \cdot \text{treatment}_i}_{\text{fixed effects}} + \underbrace{b_{0j}}_{\text{random intercept for site } j} + \varepsilon_{ijk}\]

where:

• \(\beta_0\) = population-average temperature for the reference treatment (Sun)
• \(\beta_1\) = fixed effect of shade (the quantity we most care about)
• \(b_{0j} \sim \mathcal{N}(0, \sigma^2_{\text{site}})\) = site-specific deviation from the mean
• \(\varepsilon_{ijk} \sim \mathcal{N}(0, \sigma^2)\) = within-site residual

In lme4 formula notation this is: temp_C ~ treatment + (1 | site)

Fitting the model

library(lme4)
library(lmerTest)   # adds Satterthwaite p-values to lmer output

model1 <- lmer(
  temp_C ~ treatment + (1 | site),
  data = sim_list,
  REML = TRUE         # Restricted Maximum Likelihood — preferred for variance estimation
)

summary(model1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: temp_C ~ treatment + (1 | site)
##    Data: sim_list
## 
## REML criterion at convergence: 30770.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5548 -0.8033 -0.4120  0.8309  2.5930 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  site     (Intercept) 1.348    1.161   
##  Residual             5.007    2.238   
## Number of obs: 6912, groups:  site, 3
## 
## Fixed effects:
##                  Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)      24.74583    0.67143    2.00645   36.85 0.000722 ***
## treatmentShade   -4.03872    0.05383 6908.00000  -75.03  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## treatmntShd -0.040

Interpreting the output

The summary() output has two main sections:

Fixed effects — these are the coefficients you care about most:

• (Intercept) — estimated mean temperature for the Sun group (the reference level), averaged across all sites.
• treatmentShade — the estimated change in temperature when moving from Sun to Shade. A negative value means shade is cooler.

Random effects — the variance components:

• site variance — how much baseline temperature varies across sites, independent of treatment. A large value here means the sites you chose are quite different from each other.
• Residual variance — within-site, within-sensor noise after accounting for the fixed and random effects.

vc <- as.data.frame(VarCorr(model1))
icc <- vc$vcov[1] / sum(vc$vcov)
round(icc, 3)

Model comparison (optional but recommended)

How do we know the random intercept for site is necessary? Compare the model with and without it using a likelihood ratio test:

model_null <- lm(temp_C ~ treatment, data = sim_list)   # no random effect

# AIC comparison — lower is better
AIC(model_null, model1)

##            df      AIC
## model_null  3 31895.40
## model1      4 30778.36

If model1 has a substantially lower AIC (rule of thumb: difference > 4), the random effect improves model fit and should be retained.

Residual diagnostics

Always check that the model assumptions are reasonably met. The two most important plots are:

par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))

# Residuals vs fitted
plot(fitted(model1), residuals(model1),
     pch = 16, cex = 0.5, col = "#4a7c5e88",
     xlab = "Fitted values", ylab = "Residuals",
     main = "Residuals vs. Fitted")
abline(h = 0, lty = 2, col = "#c8832a")

# Normal Q-Q
qqnorm(residuals(model1), pch = 16, cex = 0.5, col = "#4a7c5e88",
       main = "Normal Q–Q")
qqline(residuals(model1), col = "#c8832a", lty = 2)

Left: residuals vs. fitted values, checking for heteroscedasticity. Right: Q–Q plot of residuals, checking for normality. Both look reasonable for this simulated dataset.

par(mfrow = c(1, 1))

What to look for:

• Residuals vs. Fitted: points should scatter randomly around zero with no obvious fan-shape (which would indicate heteroscedasticity).
• Q–Q plot: points should fall close to the diagonal line. Systematic curves suggest non-normality. Slight deviations at the tails are common and rarely a serious problem with large datasets.

For more detail on these diagnostic plots and what to do when assumptions are violated, see Module 11 — Multiple Regression & GLMs.

Effect size and confidence intervals

The model coefficient gives you the point estimate for the shade effect, but you should also report uncertainty. Use confint():

confint(model1, method = "Wald")   # Wald CIs are fast; use "profile" for publication

##                    2.5 %    97.5 %
## .sig01                NA        NA
## .sigma                NA        NA
## (Intercept)    23.429852 26.061815
## treatmentShade -4.144223 -3.933208

Report the shade effect as: “Shaded locations were on average X.X °C cooler than sun-exposed locations (95% CI: [lower, upper]; p = …, linear mixed-effects model with random intercept for site).”

Visualizing the model output

A clean way to communicate the fixed effect is to plot the raw data with the model-estimated group means overlaid:

# Extract fixed-effect predictions with confidence intervals
library(broom.mixed)

fe_ci <- broom.mixed::tidy(model1, effects = "fixed", conf.int = TRUE)

# Reference level is Sun; add Shade prediction manually
intercept <- fe_ci$estimate[fe_ci$term == "(Intercept)"]
shade_est <- fe_ci$estimate[fe_ci$term == "treatmentShade"]
shade_ci_lo <- fe_ci$conf.low[fe_ci$term == "treatmentShade"]
shade_ci_hi <- fe_ci$conf.high[fe_ci$term == "treatmentShade"]

pred_df <- data.frame(
  treatment = c("Sun", "Shade"),
  mean_est  = c(intercept, intercept + shade_est),
  ci_lo     = c(fe_ci$conf.low[1],  intercept + shade_ci_lo),
  ci_hi     = c(fe_ci$conf.high[1], intercept + shade_ci_hi)
) |> mutate(treatment = factor(treatment, levels = c("Sun", "Shade")))

ggplot() +
  geom_jitter(
    data  = sim_list,
    aes(x = treatment, y = temp_C, colour = treatment),
    width = 0.15, alpha = 0.08, size = 0.6
  ) +
  geom_errorbar(
    data = pred_df,
    aes(x = treatment, ymin = ci_lo, ymax = ci_hi),
    width = 0.12, linewidth = 1.1, colour = "#1a1a18"
  ) +
  geom_point(
    data  = pred_df,
    aes(x = treatment, y = mean_est, colour = treatment),
    size  = 4, shape = 18
  ) +
  scale_colour_manual(values = c(Sun = "#c8832a", Shade = "#2d5a3d"),
                      guide  = "none") +
  labs(
    x = NULL, y = "Air temperature (°C)",
    title    = "Estimated shade treatment effect",
    subtitle = "Diamonds = model-estimated means ± 95% CI; dots = raw observations"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title    = element_text(family = "serif", face = "bold", colour = "#1e3a2f"),
    plot.subtitle = element_text(colour = "#7a7a72", size = 10)
  )

Group means (±95% CI) from the mixed-effects model, overlaid on individual sensor readings (jittered). The model-estimated effect of shade is shown by the vertical gap between the two means.

Adding time as a covariate

The model above treats all time points equally. But you likely noticed a strong diel cycle — temperatures are far higher at 2 pm than at 6 am. You can add hour of day as a fixed effect to partial out this variation:

sensor_merged <- sensor_merged |>
  mutate(hour = as.numeric(format(datetime, "%H")))

model_time <- lmer(
  temp_C ~ treatment + hour + (1 | site),
  data = sensor_merged,
  REML = TRUE
)

summary(model_time)

Adding hour as a linear covariate improves the model fit and gives you a cleaner estimate of the shade effect by removing the systematic time-of-day trend from the residuals. For a cyclic covariate like hour of day, you might also consider sin and cos transformations (see Module 14 — Time Series).

Random slopes

The model above assumes the shade effect is identical at every site. To test whether the effect varies across sites, add a random slope:

model_rs <- lmer(
  temp_C ~ treatment + (treatment | site),
  data = sensor_merged,
  REML = TRUE
)

# Compare to random-intercept model
anova(model1, model_rs)

If model_rs fits significantly better, the shade effect is genuinely heterogeneous — some sites benefit more from shade than others.

Mixed-effects model deep dive

This activity has only scratched the surface of mixed-effects models. For a thorough treatment of random intercepts, random slopes, nested versus crossed effects, and model selection strategies, work through: