Project: Air Quality Report

In this project you will put the whole course together and analyse a real data set from start to finish, the way an analyst does at work. The data is airquality, which comes with R: daily air measurements taken in New York from May to September 1973. You will inspect it, clean it, summarise it by month, draw charts, test a question with statistics, fit a regression model, flag high-pollution days and finish with a short written report.

Every step runs in the browser. Each code block is self-contained, so you can run any step on its own; in RStudio or Positron you would put all the steps, in order, into one script or a Quarto document.

The question we want to answer: when and why is ozone high? Ozone near the ground is a pollutant that forms when sunlight acts on other pollutants, so we expect it to depend on the weather.

Step 1: meet the data

Always start by looking at the size, the columns and the gaps:

aq <- airquality
dim(aq)
str(aq)
colSums(is.na(aq))

[1] 153   6
'data.frame':	153 obs. of  6 variables:
 $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
 $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
 $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
 $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
 $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
 $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...
  Ozone Solar.R    Wind    Temp   Month     Day 
     37       7       0       0       0       0 

There are 153 days and six columns. The help page (?airquality in your IDE) tells you the units: ozone in parts per billion (ppb), solar radiation in langleys, wind in miles per hour and temperature in degrees Fahrenheit. Ozone is missing on 37 days, and solar radiation on 7. Write that down: it will go into the report.

Step 2: tidy the data

The month is stored as a number and the date is split into two columns. Build a real date, a month factor with readable labels in calendar order, and check where the gaps are:

aq <- airquality
aq$date <- as.Date(paste(1973, aq$Month, aq$Day, sep = "-"))
aq$month <- factor(month.abb[aq$Month], levels = month.abb[5:9])
aq$weekday <- weekdays(aq$date)
head(aq[, c("date", "month", "Ozone", "Temp", "Wind")], 4)
table(aq$month, useNA = "ifany")
round(100 * mean(is.na(aq$Ozone)), 1)
tapply(is.na(aq$Ozone), aq$month, sum)

        date month Ozone Temp Wind
1 1973-05-01   May    41   67  7.4
2 1973-05-02   May    36   72  8.0
3 1973-05-03   May    12   74 12.6
4 1973-05-04   May    18   62 11.5

May Jun Jul Aug Sep 
 31  30  31  31  30 
[1] 24.2
May Jun Jul Aug Sep 
  5  21   5   5   1 

• month.abb[aq$Month] turns 5 into "May", and setting levels keeps May to September in calendar order instead of alphabetical order (see Factors).
• Every month is complete in the calendar: 30 or 31 days.
• But 24.2% of ozone readings are missing, and 21 of the 37 gaps are in June. June results will be less reliable than the others. This is the kind of detail that changes how you word your conclusions.

We will not fill the gaps with invented values. Instead each step uses na.rm = TRUE or lets the function drop incomplete rows, and the report says so (see Missing Values).

Step 3: summarise by month

aq <- airquality
aq$month <- factor(month.abb[aq$Month], levels = month.abb[5:9])
by_month <- aggregate(cbind(Ozone, Temp, Wind) ~ month, data = aq,
                      FUN = function(x) round(mean(x), 1))
by_month
median_ozone <- tapply(aq$Ozone, aq$month, median, na.rm = TRUE)
median_ozone
summary(aq$Ozone)

  month Ozone Temp Wind
1   May  23.6 66.7 11.5
2   Jun  29.4 78.2 12.2
3   Jul  59.1 83.9  8.5
4   Aug  60.0 84.0  8.6
5   Sep  31.4 76.9 10.1
May Jun Jul Aug Sep 
 18  23  60  52  23 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     NAs 
   1.00   18.00   31.50   42.13   63.25  168.00      37 

July and August stand out: mean ozone around 60 ppb, more than double May's, and they are also the hottest and least windy months.

A careful reader would spot a trap here. The formula version of aggregate() drops every row with a missing value in any of the listed columns, so the temperature and wind means are computed only on days that also have an ozone reading. That is why May's mean temperature shows 66.7 here, while the full-month mean (in Step 8) is 65.5. Both numbers are correct; they are means of different sets of days. Know which one you are reporting.

Ozone is also skewed: the mean (42.1) is well above the median (31.5), pulled up by a few very high days, up to 168 ppb. For skewed data the median is the safer "typical" value, as you learned in Descriptive Statistics.

Step 4: look at the data

aq <- airquality
aq$month <- factor(month.abb[aq$Month], levels = month.abb[5:9])
par(mfrow = c(1, 2))
boxplot(Ozone ~ month, data = aq, col = "lightblue",
        main = "Ozone by month", xlab = "", ylab = "Ozone (ppb)")
plot(aq$Temp, aq$Ozone, pch = 19, col = "darkorange",
     main = "Ozone vs temperature", xlab = "Temperature (F)", ylab = "Ozone (ppb)")
abline(lm(Ozone ~ Temp, data = aq), lwd = 2)
par(mfrow = c(1, 1))
round(cor(aq$Temp, aq$Ozone, use = "complete.obs"), 2)
round(cor(aq$Wind, aq$Ozone, use = "complete.obs"), 2)

[1] 0.7
[1] -0.6

Two charts appear side by side. On the left, five light blue boxes show ozone by month: the July and August boxes sit much higher and are taller (more variable) than May, June and September, and a few outlier circles appear above the boxes. On the right, orange points rise from the bottom left to the top right with a black straight line through them: hotter days have more ozone. The points fan out at high temperatures, and at the hot end several sit well above the line, a hint that the relationship curves upward.

The correlations put numbers on what we see: ozone rises with temperature (r = 0.70) and falls with wind (r = −0.60). Wind spreads pollution out; heat and sunshine help ozone form.

Step 5: test the question

Is ozone really higher in the peak months, July and August, than in the rest of the season? The ozone values are skewed and contain outliers, so we use the rank-based Wilcoxon test instead of a t-test (see Hypothesis Tests):

aq <- airquality
aq$season <- ifelse(aq$Month %in% c(7, 8), "Peak (Jul-Aug)", "Other (May, Jun, Sep)")
tapply(aq$Ozone, aq$season, median, na.rm = TRUE)
wilcox.test(Ozone ~ season, data = aq, exact = FALSE)

Other (May, Jun, Sep)        Peak (Jul-Aug) 
                   21                    59 

	Wilcoxon rank sum test with continuity correction

data:  Ozone by season
W = 724.5, p-value = 1.847e-07
alternative hypothesis: true location shift is not equal to 0

The median ozone is 59 ppb in July and August against 21 ppb in the other months, and the p-value (about 0.0000002) is tiny. The difference is far too large to be chance. exact = FALSE tells R to use the normal approximation, because the data has tied values and an exact p-value cannot be computed.

Step 6: model ozone with the weather

A regression model estimates the effect of temperature and wind together (see Linear Regression):

aq <- airquality
model <- lm(Ozone ~ Temp + Wind, data = aq)
round(summary(model)$coefficients, 3)
round(summary(model)$adj.r.squared, 3)
nobs(model)
log_model <- lm(log(Ozone) ~ Temp + Wind, data = aq)
round(coef(log_model), 3)
round(summary(log_model)$adj.r.squared, 3)
predict(model, data.frame(Temp = c(70, 90), Wind = c(10, 5)))

            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -71.033     23.578  -3.013    0.003
Temp           1.840      0.250   7.362    0.000
Wind          -3.055      0.663  -4.607    0.000
[1] 0.561
[1] 116
(Intercept)        Temp        Wind 
     -0.532       0.057      -0.053 
[1] 0.575
       1        2 
27.22439 79.30542 

• Holding wind constant, each extra degree goes with about 1.8 ppb more ozone; holding temperature constant, each extra mph of wind goes with about 3.1 ppb less. Both p-values are far below 0.05 (they round to 0.000 at three decimals).
• The model uses 116 days, because days without an ozone reading are dropped. Adjusted R-squared is 0.56: the two weather variables explain a little over half of the day-to-day variation.
• Because ozone is skewed and the scatter plot fanned out, we also tried modelling log(Ozone). Its fit is slightly better (0.575). On the log scale, each extra degree multiplies ozone by about e^0.057 ≈ 1.06, so ozone rises by roughly 6% per degree. The log model also never predicts negative ozone, which the straight-line model can do on cold, windy days.
• The predictions show the size of the weather effect: a mild, breezy day (70 °F, 10 mph) gives about 27 ppb; a hot, still day (90 °F, 5 mph) about 79 ppb.

Step 7: flag high-ozone days

For this practice report we call a day high-ozone when the reading is above 80 ppb. This cut-off is our own choice for the exercise, not an official air-quality standard:

aq <- airquality
aq$date <- as.Date(paste(1973, aq$Month, aq$Day, sep = "-"))
high <- which(aq$Ozone > 80)
length(high)
round(mean(aq$Temp[high]), 1)
round(mean(aq$Temp[-high]), 1)
head(aq[high, c("date", "Ozone", "Temp", "Wind")][order(-aq$Ozone[high]), ], 5)

[1] 16
[1] 89.2
[1] 76.6
          date Ozone Temp Wind
117 1973-08-25   168   81  3.4
62  1973-07-01   135   84  4.1
99  1973-08-07   122   89  4.0
121 1973-08-29   118   94  2.3
30  1973-05-30   115   79  5.7

• which() is used on purpose: aq$Ozone > 80 contains NA for the missing days, and which() ignores them, avoiding the junk-row trap from Indexing & Subsetting.
• There are 16 high-ozone days. Their average temperature is 89.2 °F, compared with 76.6 °F on the other days.
• The five worst days all had very light wind (2.3 to 5.7 mph). The worst, 25 August, was not especially hot at 81 °F but almost still. Wind matters as much as heat.

Step 8: write the report

Finally, produce a compact monthly table in text and save the summary as a CSV file for others:

aq <- airquality
aq$month <- factor(month.abb[aq$Month], levels = month.abb[5:9])
report_line <- function(m) {
  d <- aq[aq$month == m, ]
  sprintf("%-4s %3d days | ozone median %5.1f ppb (%2d missing) | mean temp %4.1f F | high-ozone days %2d",
          m, nrow(d), median(d$Ozone, na.rm = TRUE), sum(is.na(d$Ozone)),
          mean(d$Temp), sum(d$Ozone > 80, na.rm = TRUE))
}
writeLines(vapply(levels(aq$month), report_line, character(1)))
summary_table <- data.frame(
  month = levels(aq$month),
  median_ozone = as.numeric(tapply(aq$Ozone, aq$month, median, na.rm = TRUE)),
  high_days = as.numeric(tapply(aq$Ozone > 80, aq$month, sum, na.rm = TRUE))
)
path <- tempfile(fileext = ".csv")
write.csv(summary_table, path, row.names = FALSE)
read.csv(path)

May   31 days | ozone median  18.0 ppb ( 5 missing) | mean temp 65.5 F | high-ozone days  1
Jun   30 days | ozone median  23.0 ppb (21 missing) | mean temp 79.1 F | high-ozone days  0
Jul   31 days | ozone median  60.0 ppb ( 5 missing) | mean temp 83.9 F | high-ozone days  6
Aug   31 days | ozone median  52.0 ppb ( 5 missing) | mean temp 84.0 F | high-ozone days  7
Sep   30 days | ozone median  23.0 ppb ( 1 missing) | mean temp 76.9 F | high-ozone days  2
  month median_ozone high_days
1   May           18         1
2   Jun           23         0
3   Jul           60         6
4   Aug           52         7
5   Sep           23         2

report_line() builds one line per month with sprintf() widths so the columns line up (see Strings & Regular Expressions), and vapply() runs it for every month (see The apply Family). On your computer, write the CSV to a real file name in your project folder instead of tempfile().

Here is the written summary that goes with the numbers:

Ozone in New York, May to September 1973 (R built-in data set airquality)

1. Ozone was highest in July and August (median 60 and 52 ppb) and lowest in
   May (18 ppb). The peak-month difference is statistically significant
   (Wilcoxon test, p < 0.001).
2. Ozone rises with temperature (r = 0.70) and falls with wind (r = -0.60).
   A regression on both explains about 56% of the daily variation: about
   +1.8 ppb per degree F and -3.1 ppb per mph of wind.
3. 16 days had ozone above 80 ppb (a cut-off chosen for this report). They
   averaged 89 F, and the five worst days all had wind below 6 mph.
4. Limitations: 37 of 153 ozone readings (24%) are missing, 21 of them in
   June, so June figures are uncertain. Missing days were left out, not
   filled. This is observational data from one city and one summer: it shows
   associations, not proof of cause.

The limitations paragraph is not optional. Stating what is missing and what the analysis cannot show is what makes a report trustworthy.

Taking it further

• Put Steps 1 to 8 in a Quarto document with #| echo: false chunks and inline numbers, as in R Markdown & Quarto, so the report rebuilds itself from the data.
• Redo Steps 3 and 7 with dplyr verbs (group_by(month) |> summarise(...)) and the charts with ggplot2 (facet_wrap(~ month)), as in The tidyverse & dplyr and ggplot2 Basics.
• Add Solar.R to the regression and see whether adjusted R-squared improves.

Key points

• Start every analysis by checking size, types and missing values, and record the gaps.
• Tidy the data first: real dates, factors in a sensible order, clear column names.
• Summaries, charts, tests and models should tell one consistent story; here, heat and still air go with high ozone.
• Choose methods that fit the data: medians and the Wilcoxon test for skewed values, which() around conditions with NA, a log model for a fanning scatter.
• A good report states its findings with numbers, and its limitations just as clearly.

aq <- airquality
sum(is.na(aq$Solar.R))
round(cor(aq$Solar.R, aq$Ozone, use = "complete.obs"), 2)
m3 <- lm(Ozone ~ Temp + Wind + Solar.R, data = aq)
nobs(m3)
round(summary(m3)$adj.r.squared, 3)
round(coef(m3)["Solar.R"], 3)

[1] 7
[1] 0.35
[1] 111
[1] 0.595
Solar.R 
   0.06