Multiple regression

Materials for class on Thursday, October 18, 2018

Contents

Slides

Download the slides from today’s lecture:

• PDF (best option)
• PowerPoint

Predicting temperature

At the end of class today, we looked at the factors that predict the daily high temperature in Provo (specifically 40.248752, -111.649216) from January to May 2017. Here’s a cleaner version of that analysis.

Load and wrangle data

First, we load the libraries we’ll be using:

library(tidyverse)
library(lubridate)
library(moderndive)
library(scales)  # For nicer axis labels

Then we load the data. I used the darksky R package to download this historical data from Dark Sky, and then I saved the CSV file online.

# This reads the CSV file directly from the internet
weather_provo_raw <- read_csv("https://andhs.co/provoweather")

Then we clean and wrangle the data, or make adjustments to it so that it’s easier to use and model and plot. I include comments explaining what each line does below:

weather_provo_2017 <- weather_provo_raw %>% 
         # Extract the month name from the date
  mutate(month = month(date, label = TRUE, abbr = FALSE),
         # Extract the month number from the date 
         month_number = month(date, label = FALSE),
         # Extract the day of the week from the date
         weekday = wday(date, label = TRUE, abbr = FALSE),
         # Extract the number of the day of the week (1 = Sunday)
         weekday_number = wday(date, label = FALSE)) %>% 
  # Change all the missing values in the precipType column to "none"
  mutate(precipType = ifelse(is.na(precipType), "none", precipType)) %>% 
  # Only select these columns
  select(date, month, month_number, weekday, weekday_number,
         sunriseTime, sunsetTime, moonPhase, precipProbability, precipType,
         temperatureHigh, temperatureLow, dewPoint, humidity, pressure, 
         windSpeed, cloudCover, visibility, uvIndex)

# Make a subset of the clean data that only includes January through May
winter_spring <- weather_provo_2017 %>% 
  filter(month_number <= 5) %>% 
  # After we filter, R still thinks that the month column has 12 possible
  # categories. Here, we tell it to be a factor again (it already is one) so
  # that it only includes January through May
  mutate(month = factor(month, ordered = FALSE)) %>% 
  # Finally, we multiply some of the percentage variables by 100 so the
  # regression coefficients are easier to interpret
  mutate(humidity = humidity * 100,
         cloudCover = cloudCover * 100,
         precipProbability = precipProbability * 100)

Explore data

First we can look at how the high temperature changed over time. Not surprisingly, it went up.

ggplot(winter_spring, aes(x = date, y = temperatureHigh)) +
  geom_line() +
  geom_smooth(method = "lm", se = FALSE) +
  scale_y_continuous(labels = degree_format()) +
  labs(x = NULL, y = "High temperature") +
  theme_minimal()

Next we look at how humidity has changed over time. It went down:

ggplot(winter_spring, aes(x = date, y = humidity)) +
  geom_line() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(x = NULL, y = "Humidity") +
  theme_minimal()

Neat.

Humidity and temperature

Our main research question is exploring the relationship between humidity and temperature. Are they correlated, and if so, how much does temperature change as humidity increases or decreases?

Our outcome variable (or y or response or dependent variable) is thus temperatureHigh, while our main explanatory variable (or x or independent variable) is humidity.

First we check to see how correlated they are, both visually and with a correlation coefficient:

ggplot(winter_spring, aes(x = humidity, y = temperatureHigh)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  scale_y_continuous(labels = degree_format()) +
  labs(x = "Humidity", y = "High temperature") +
  theme_minimal()

The graph shows that the two variables are definitely negatively correlated (as humidity increases, temperature decreases). The correlation between the two is…

winter_spring %>% 
  get_correlation(temperatureHigh ~ humidity)

## # A tibble: 1 x 1
##   correlation
##         <dbl>
## 1      -0.821

…-0.821. That’s close to 1, so we can say that temperature and humidity are strongly negatively correlated (r = -0.821).

Next, we want to see how much temperature changes as humidity changes. We use a linear model to find the slope and intercept of this line:

model1 <- lm(temperatureHigh ~ humidity, data = winter_spring)
model1 %>% get_regression_table()

## # A tibble: 2 x 7
##   term      estimate std_error statistic p_value lower_ci upper_ci
##   <chr>        <dbl>     <dbl>     <dbl>   <dbl>    <dbl>    <dbl>
## 1 intercept   97.0       2.62       37.0       0   91.9    102.   
## 2 humidity    -0.856     0.049     -17.6       0   -0.952   -0.759

First, we see if the intercept (\(\beta_0\)) provides helpful information. The intercept is 97, which means that on average, during the period of January to May 2017, if humidity was 0, the predicted temperature would be 97˚. That seems… really really high. There were 0 days with 0 humidity during this time period (the lowest humidity level was 25), so it’s a stretch to try to interpret this for non-existent extra-humid days. We can safely disregard the intercept.

The \(\beta_1\) coefficient, or slope of humidity, shows that for every 1 point increase in humidity, there is an associated drop in temperature of 0.856˚, on average. As the weather becomes more humid, the high daily temperature becomes colder. This is probably because humidity in the winter is a sign of snow—I’m guessing this relationship is reversed in summer months, where humidity leads to muggier weather.

Months and temperature

Next, we want to see what the effect of time is on temperature. While humidity explains some of the variation in temperature, from previous experience in being alive, we know that the weather naturally warms up as we approach summer, regardless of humidity or any other factor. We can see this if we look at the average temperature per month:

winter_spring %>% 
  group_by(month) %>% 
  summarize(avg_high = mean(temperatureHigh))

## # A tibble: 5 x 2
##   month    avg_high
##   <fct>       <dbl>
## 1 January      33.9
## 2 February     45.7
## 3 March        56.6
## 4 April        56.8
## 5 May          69.3

Not surprisingly, January was really cold on average, while May was not. The other months fall somewhere in between.

If we include month in the regression model, we’ll get a coefficient for each month except the base case, which is the first category (January in this case):

model2 <- lm(temperatureHigh ~ month, data = winter_spring)
model2 %>% get_regression_table()

## # A tibble: 5 x 7
##   term          estimate std_error statistic p_value lower_ci upper_ci
##   <chr>            <dbl>     <dbl>     <dbl>   <dbl>    <dbl>    <dbl>
## 1 intercept         33.9      1.65     20.5        0    30.6      37.2
## 2 monthFebruary     11.8      2.40      4.94       0     7.11     16.6
## 3 monthMarch        22.8      2.34      9.74       0    18.1      27.4
## 4 monthApril        22.9      2.35      9.75       0    18.3      27.6
## 5 monthMay          35.4      2.34     15.2        0    30.8      40.0

These coefficients are not slopes. They are shifts in the intercept. Notice how the intercept is 33.9—identical to January’s average high temperature. All other month-based coefficients show the adjustment in average temperature per month using January as the base case. February is 11.8˚ warmer than January, on average. March is 22.8˚ warmer than January, on average. April is 22.9˚ warmer than January, on average. May is 35.4˚ warmer than January, on average.

Humidity, months, and temperature

Both humidity and month explain variation in temperatures. This is obvious in this scatterplot, which colors the points by month. Points for May are clustered in the top left corner, while points in January are clustered in the bottom right.

ggplot(winter_spring, aes(x = humidity, y = temperatureHigh)) +
  geom_point(aes(color = month)) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_y_continuous(labels = degree_format()) +
  labs(x = "Humidity", y = "High temperature") +
  theme_minimal()

If we include both humidity and month in the model, we can get an estimate of the effect of humidity controlling for or taking into account time (or month).

model3 <- lm(temperatureHigh ~ humidity + month, data = winter_spring)
model3 %>% get_regression_table()

## # A tibble: 6 x 7
##   term          estimate std_error statistic p_value lower_ci upper_ci
##   <chr>            <dbl>     <dbl>     <dbl>   <dbl>    <dbl>    <dbl>
## 1 intercept        71.8      3.70      19.4    0       64.5     79.1  
## 2 humidity         -0.56     0.052    -10.9    0       -0.662   -0.458
## 3 monthFebruary     6.44     1.85       3.48   0.001    2.78    10.1  
## 4 monthMarch       10.9      2.05       5.30   0        6.83    15.0  
## 5 monthApril       11.2      2.06       5.43   0        7.11    15.3  
## 6 monthMay         20.7      2.20       9.39   0       16.3     25.1

Now we have a slope for humidity, as well as several shifts in intercept based on month. We can start by interpreting the \(\beta_0\) coefficient. This shows the average temperature in January when humidity is at 0. 71˚ still seems really high for January, but recall that there are no days with 0 humidity, so it’s a stretch to try to get anything meaningful out of the intercept.

We can make the following statements about this model:

• Controlling for month, a one point increase in humidity is associated with a 0.56˚ drop in temperature, on average.
• Controlling for humidity, the high temperature in February is 6.4˚ warmer than in January, on average. Note that this is a lot different than the \(\beta\) from the previous model (there it was 11.8). That’s because humidity is now explaining some of the variation in temperature.
• Controlling for humidity, the high temperature in March is 10.9˚ warmer than in January, on average.
• Taking humidity into account, the high temperature in April is 11.2˚ warmer than in January, on average.
• On average, if we take humidity into account, the high temperature in May is 20.7˚ higher than in January.

Humidity, wind, and clouds

For fun, we can see what the effect of other weather-related variables is on high temperature. Here we build a model using humidity, wind speed, and cloud cover to explain variation in temperature.

model4 <- lm(temperatureHigh ~ humidity + windSpeed + cloudCover,
             data = winter_spring)
model4 %>% get_regression_table()

## # A tibble: 4 x 7
##   term       estimate std_error statistic p_value lower_ci upper_ci
##   <chr>         <dbl>     <dbl>     <dbl>   <dbl>    <dbl>    <dbl>
## 1 intercept   103.        3.67      28.0    0       95.5    110.   
## 2 humidity     -0.992     0.08     -12.4    0       -1.15    -0.833
## 3 windSpeed    -0.959     0.533     -1.80   0.074   -2.01     0.093
## 4 cloudCover    0.074     0.036      2.04   0.043    0.002    0.146

We can make the following statements:

• Controlling for wind speed and cloud cover, a one unit increase in humidity is associated with a 0.992˚ lower temperature, on average.
• Taking humidity and cloud cover into account, a 1 MPH increase in wind speed is associated with a lower temperature of 0.959˚, on average.
• On average, controlling for humidity and wind speed, a one percent increase in cloud cover is associated wig a higher temperature of 0.074˚.

Humidity, wind, clouds, and month

We can also account for month when looking at humidity, wind speed, and cloud cover:

model5 <- lm(temperatureHigh ~ humidity + windSpeed + cloudCover + month,
             data = winter_spring)
model5 %>% get_regression_table()

## # A tibble: 8 x 7
##   term          estimate std_error statistic p_value lower_ci upper_ci
##   <chr>            <dbl>     <dbl>     <dbl>   <dbl>    <dbl>    <dbl>
## 1 intercept       74.5       4.41     16.9     0       65.7     83.2  
## 2 humidity        -0.658     0.074    -8.85    0       -0.805   -0.511
## 3 windSpeed       -0.101     0.438    -0.231   0.817   -0.967    0.764
## 4 cloudCover       0.06      0.029     2.07    0.04     0.003    0.117
## 5 monthFebruary    6.14      1.85      3.33    0.001    2.49     9.79 
## 6 monthMarch      10.7       2.04      5.23    0        6.63    14.7  
## 7 monthApril      11.1       2.04      5.43    0        7.04    15.1  
## 8 monthMay        20.8       2.21      9.40    0       16.4     25.2

Phew. Get ready for all these interpretations:

• Controlling for wind speed, cloud cover, and month, a one unit increase in humidity is associated with a 0.658 drop in temperature, on average.
• Taking into account all the other variables in the model, a 1 MPH increase in wind speed is associated with a 0.1˚ lower temperature.
• Taking humidity, wind speed, and month into account, on average a 1 percent increase in cloud cover is associated with a 0.06˚ higher temperature.
• Controlling for humidity, wind speed, and cloud cover, February is—on average—6.14˚ warmer than January.
• Controlling for humidity, wind speed, and cloud cover, the average high temperature in March is 10.7˚ higher than in January
• Taking humidity, wind speed, and cloud cover into account, April is 11.1˚ warmer than January, on average.
• After controlling for humidity, wind speed, and cloud cover, we can see that high temperatures in May are an average of 20.8˚ higher than in January.

The kitchen sink

For fun, here’s a model that explains variation in temperature with humidity, wind speed, cloud cover, the probability of precipitation, visibility, and month. Oh baby.

model6 <- lm(temperatureHigh ~ humidity + windSpeed + cloudCover + 
               precipProbability + visibility + month,
             data = winter_spring)
model6 %>% get_regression_table()

## # A tibble: 10 x 7
##    term              estimate std_error statistic p_value lower_ci upper_ci
##    <chr>                <dbl>     <dbl>     <dbl>   <dbl>    <dbl>    <dbl>
##  1 intercept           88.2       8.50      10.4    0       71.4    105.   
##  2 humidity            -0.764     0.081     -9.48   0       -0.924   -0.605
##  3 windSpeed           -0.039     0.429     -0.09   0.928   -0.887    0.809
##  4 cloudCover           0.019     0.031      0.61   0.543   -0.043    0.081
##  5 precipProbability    0.112     0.041      2.77   0.006    0.032    0.192
##  6 visibility          -0.8       0.74      -1.08   0.282   -2.26     0.664
##  7 monthFebruary        5.90      1.93       3.06   0.003    2.08     9.71 
##  8 monthMarch           9.62      2.20       4.36   0        5.26    14.0  
##  9 monthApril           9.32      2.30       4.06   0        4.78    13.9  
## 10 monthMay            19.4       2.38       8.17   0       14.7     24.1

SO MANY COEFFICIENTS. I’m not going to interpret them all here. I’m sick of typing. Here are few:

• Controlling for all the other variables in the model, a one unit increase in humidity is associated with a 0.764˚ lower temperature, on average.
• Taking into account humidity, wind speed, cloud cover, visibility, and month, a 1% greater likelihood of precipitation is associated with a 0.112˚ higher temperature, on average.
• Controlling for all the variables in the model, February is 5.9˚ warmer than January, on average.

All coefficients from all models

For fun, we can compare the coefficients from all of these models simultaneously, similar to what we did in class with the regression that looked at the effect of small class sizes on test scores. If you install the huxtable package and use the huxreg() function, you can make the same sort of vertical, one-column-per-model table:

library(huxtable)

huxreg(model1, model2, model3, model4, model5, model6)

If the main goal of this whole exercise was to measure the effect of humidity on temperature, we’ve done a fairly good job. In every single model (except #2, where we didn’t include it), the coefficient for humidity is consistently negative and ranges between 0.56 and 0.92.

In general, during the winter in Utah, if you see humidity rising, get your coat.

Clearest and muddiest things

Go to this form and answer these three questions:

1. What was the muddiest thing from class today? What are you still wondering about?
2. What was the clearest thing from class today?
3. What was the most exciting thing you learned?

I’ll compile the questions and send out answers after class.