I've got a dataset here called airquality. It's a built-in R dataset containing daily air quality measurements in New York from May to September 1973. It's a classic, but it's also a bit "dirty," which makes it perfect for how we actually work in the real world.

Wait, why is my mean NA?

I want to get a quick sense of the average Ozone levels. Seems straightforward enough. Let's just run a mean on that column.

mean(airquality$Ozone)
# [1] NA

Right. Of course. That's the "R experience." If there is a single NA (Not Available) value in your vector, R conservatively returns NA for the whole calculation because it doesn't want to assume how you'd prefer to handle missing data. It's annoying when you just want a quick number, but it's a safety feature.

I need to tell R to just ignore those gaps. I'll use the na.rm argument (which stands for "NA remove").

mean(airquality$Ozone, na.rm = TRUE)
# [1] 42.12931

There we go. About 42 ppb. Now I can actually work with the numbers.

Narrowing down the smoggy days

Average is fine, but I'm more interested in the extremes. I want to see the days where the air quality was actually dangerous—let's say Ozone levels over 100. I'll try to filter the data frame.

smoggy_days <- airquality[airquality$Ozone > 100, ]
head(smoggy_days)

If I look at the output, I see something weird. I have some rows where the values are NA. Wait, if I asked for Ozone > 100, why am I seeing rows where Ozone is missing? It's because in R, NA > 100 doesn't result in FALSE; it results in NA. When you use that in a row index, R keeps the NA rows because it isn't sure if they should be there or not.

I need to be explicit. I only want rows where Ozone is greater than 100 AND where Ozone is not NA. I'll use subset() here because it handles NAs much more gracefully than the bracket notation.

smoggy_days <- subset(airquality, Ozone > 100)
nrow(smoggy_days)
# [1] 37

37 days. That feels more accurate. Much cleaner.

Seeing if heat drives the pollution

I have a hunch that Ozone spikes are tied to high temperatures. I could run a correlation, but I want to see it first. I'll throw a quick scatter plot together.

plot(airquality$Temp, airquality$Ozone)

The plot shows a clear upward trend, but it's a bit cluttered. I want to highlight those "dangerous" days we just found. I'll plot the whole set in one color, then overlay the smoggy days in red.

plot(airquality$Temp, airquality$Ozone, col = "gray", pch = 16, 
     xlab = "Temperature (F)", ylab = "Ozone (ppb)", main = "Temp vs Ozone")
points(smoggy_days$Temp, smoggy_days$Ozone, col = "red", pch = 16)

Now the relationship is obvious. The red dots are almost all clustered in the top right. To put a hard number on this, I'll use cor(), but again, I have to remember my lesson from the start of this session: the NAs will kill my result if I'm not careful.

cor(airquality$Temp, airquality$Ozone, use = "complete.obs")
# [1] 0.701141

A correlation of 0.7 is a strong positive relationship. Using use = "complete.obs" tells R to only look at pairs where both values are present. If I had left that out, I'd be right back at that frustrating NA result.

📋 Practical Task

Analyzing New York's July Heatwave Peaks

Using the built-in airquality dataset, write a script to perform the following analysis:

• Create a new data frame called july_data that contains only the observations from July (where Month == 7).
• From this July subset, calculate the average Wind speed, ensuring that missing values are ignored.
• Filter the july_data to find "Extreme Days," defined as days where the Temp was above 85 degrees AND the Ozone was above 60 ppb.
• Print the number of "Extreme Days" found in July.