The question

  • Ground-level ozone irritates the lungs. It forms when heat and sunlight drive reactions between pollutants in the air.
  • Question: on hotter days, is there more ozone? If so, how much more per degree?
  • Data: R’s built-in airquality set, 153 days of readings in New York City, May to September 1973.
    • Ozone: mean ozone (parts per billion), 1 to 3 PM, Roosevelt Island
    • Temp: maximum daily temperature (°F), LaGuardia Airport
    • Wind: average morning wind speed (mph), LaGuardia Airport

The data

aq <- na.omit(airquality[, c("Ozone", "Temp", "Wind")])
nrow(aq)      # days that have an ozone reading
## [1] 116
summary(aq)
##      Ozone             Temp            Wind       
##  Min.   :  1.00   Min.   :57.00   Min.   : 2.300  
##  1st Qu.: 18.00   1st Qu.:71.00   1st Qu.: 7.400  
##  Median : 31.50   Median :79.00   Median : 9.700  
##  Mean   : 42.13   Mean   :77.87   Mean   : 9.862  
##  3rd Qu.: 63.25   3rd Qu.:85.00   3rd Qu.:11.500  
##  Max.   :168.00   Max.   :97.00   Max.   :20.700

37 days have no ozone reading. I drop them, which leaves n = 116 days.

First look: ozone rises with temperature

The trend is clearly upward: the hottest days have the most ozone. Below about 75°F the points sit low and close together, and above 85°F they shoot up and spread out. A straight line is a reasonable first model, but keep the curve in mind.

The statistical idea: simple linear regression

We model each day’s ozone as a straight-line function of temperature plus random noise:

\[Y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \qquad \varepsilon_i \sim N(0, \sigma^2), \qquad i = 1, \dots, n\]

  • \(Y_i\) is the ozone on day \(i\), and \(x_i\) is that day’s max temperature
  • \(\beta_1\) (slope) is the change in average ozone for each extra 1°F. This is the number the question asks for.
  • \(\beta_0\) (intercept) is the average ozone at 0°F. We have no days anywhere near that, so it has no real-world meaning here.
  • \(\varepsilon_i\) is everything else that moves ozone that day (sunlight, wind, traffic)

We never see the true \(\beta_0\) and \(\beta_1\). We estimate them from the 116 days.

Mathematics: estimating the line and testing it

Least squares picks the line with the smallest \(\sum_{i=1}^{n} (y_i - \hat y_i)^2\):

\[\hat\beta_1 = \frac{\sum_{i=1}^{n}(x_i - \bar x)(y_i - \bar y)}{\sum_{i=1}^{n}(x_i - \bar x)^2}, \qquad \hat\beta_0 = \bar y - \hat\beta_1 \bar x\]

95% confidence interval for the slope, and the t-test of \(H_0: \beta_1 = 0\) (temperature has no linear effect):

\[\hat\beta_1 \pm t_{0.975,\,n-2}\;\mathrm{SE}(\hat\beta_1), \qquad t = \frac{\hat\beta_1 - 0}{\mathrm{SE}(\hat\beta_1)} \sim t_{n-2} \text{ under } H_0\]

How much it explains: \(R^2 = 1 - \dfrac{\sum (y_i - \hat y_i)^2}{\sum (y_i - \bar y)^2}\). With one predictor, \(R^2 = r^2\).

Fitting the model in R

fit <- lm(Ozone ~ Temp, data = aq)
round(summary(fit)$coefficients, 4)
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -146.9955    18.2872 -8.0382        0
## Temp           2.4287     0.2331 10.4177        0
confint(fit, "Temp")
##         2.5 %   97.5 %
## Temp 1.966871 2.890536

Fitted line: predicted ozone = -147 + 2.43 × Temp.

Is the slope real? Inference

  • Estimate: each extra 1°F goes with about 2.43 ppb more ozone on average.
  • 95% CI: (1.97, 2.89) ppb per °F. The whole interval is above 0.
  • Hypothesis test: t = 10.4 on 114 degrees of freedom, so the p-value is far below 0.001. We reject \(H_0: \beta_1 = 0\).
  • Strength: r = 0.7 and \(R^2\) = 0.49. Temperature alone explains about 49% of the day-to-day variation in ozone.
  • In everyday terms: a day 10°F hotter has about 24 ppb more ozone on average.

Checking the model: residuals

It isn’t. The gold curve bends, and the spread grows as the fitted values grow. So the straight line under-predicts the worst ozone days, and the equal-variance assumption is shaky. A model for log(Ozone) would fit better.

Adding wind: an interactive 3D view

Drag to rotate. Model: Ozone ~ Temp + Wind. Holding wind fixed, each 1°F adds about 1.84 ppb, and each extra mph of wind removes about 3.06 ppb. \(R^2\) goes up to 0.57.

What did we learn?

  • Yes, hotter days have more ozone. About 2.4 ppb more per °F, and the 95% CI is well above zero.
  • Temperature is only part of the story. It explains about half the variation. Wind explains more, because windy days blow the ozone away.
  • The relationship isn’t a straight line. The hottest days are worse than the line predicts, so a log scale would model them better.
  • This shows association, not cause. Hot days are also sunny, calm days, and sunlight drives ozone too. The data come from one city in one summer, 1973.
  • Why it matters: forecasters use temperature and wind to issue ozone alerts on hot, still days.

Sources

  • Data: R datasets::airquality. Ozone data from the New York State Department of Conservation and weather data from the National Weather Service, published in Chambers, Cleveland, Kleiner & Tukey (1983), Graphical Methods for Data Analysis, Wadsworth. See ?airquality.
  • Logo and slide style: the ASU logo and the CSS theme come from the HW3 Canvas page.
  • Software: R, ggplot2 and plotly, written in R Markdown ioslides.