Univariate visualization/summarization of quantitative variables

Research questions: What is the range of highway miles per gallon (hwy)? What are typical highway mpg? Are highway mpg’s variable across vehicles?

As you visualize quantitative variables, keep in mind:

• Center: Where is the center of the distribution? What is a typical value of the variable?
• Variability: How spread out are the values? A lot or a little?
• Shape: How are values distributed along the observed range? Is the distribution symmetric, right-skewed (tail is longer on the right, pulled further out to the right), left-skewed, bimodal, or uniform (flat)?
• Outliers: Are there any outliers, i.e. values that are unusually large/small relative to the bulk of other values?
• Context: In the context of your research, what do you learn from the plot or table? How would you describe your findings to a broad audience?

We will often talk about the distribution of a variable. This is the values the variables take and how often. You should incorporate the answers to the questions above in your description of a distribution. You should also use numerical summaries, which we’ll talk about shortly. Think about trying to describe the distribution to someone who doesn’t have the picture in front of them.

ggplot(data = mpg) +
  geom_histogram(aes(x = hwy))

ggplot(data = mpg) +
  geom_histogram(aes(x = hwy), bins = 20, fill = "lightblue") +
  labs(x = "Highway MPG", y = "Number of vehicles") +
  theme_minimal() 

mpg %>%
  summarize(
    mean_hwy = mean(hwy),
    median_hwy = median(hwy),
    sd_hwy = sd(hwy),
    IQR_hwy = IQR(hwy)
  )

Multivariate visualization/summarization

We noticed that there was variability across vehicles but didn’t investigate why. Are there other variables that might explain some of the variation? Which variables might have stronger or weaker relationships? How can we visualize these relationships?

We’re interested in the relationship between:

Response/dependent/outcome variable: the variable whose variability we would like to explain (education or income)

Predictors/explanatory variables/independent variables/covariates: variables that might explain some of the variability in the response

Our goal is to construct visualizations that allow us to examine the following features of the relationships among these variables:

• Relationship trends (form and direction of the relationship, eg. positive and linear)
  

• Relationship strength (degree of variability from the trend or how tight the points are around the envisioned trend)

• Outliers in the relationship

ggplot(data = mpg) +
  geom_point(aes(x = displ, y = hwy)) +
  labs(x = "Engine Displacement", y = "Highway MPG") +
  theme_minimal()

mpg %>% 
  summarize(corr = cor(displ, hwy))

ggplot(data = mpg) +
  geom_jitter(aes(x = displ, y = hwy)) + #noice the code change here
  labs(x = "Engine Displacement", y = "Highway MPG") +
  theme_minimal()

ggplot(data = mpg) +
  geom_point(aes(x = displ, y = hwy), alpha = .2) +
  labs(x = "Engine Displacement", y = "Highway MPG") +
  theme_minimal()

ggplot(data = mpg) +
  geom_point(aes(x = displ, y = hwy), size = .5) +
  labs(x = "Engine Displacement", y = "Highway MPG") +
  theme_minimal()

ggplot(data = diamonds) +
  geom_jitter(aes(x = carat, y = price), alpha = .2, size = .5) +
  labs(x = "Carats", y = "Price ($)") +
  theme_minimal()

ggplot(data = mpg) +
  geom_jitter(aes(x = displ, y = hwy, color = drv)) +
  labs(x = "Engine Displacement", y = "Highway MPG", color = "Drive") +
  theme_minimal()

ggplot(data = mpg) +
  geom_jitter(aes(x = displ, y = hwy)) +
  facet_wrap(vars(drv)) +
  labs(x = "Engine Displacement", y = "Highway MPG") +
  theme_minimal()

ggplot(data = mpg) +
  geom_boxplot(aes(x = drv, y = hwy)) +
  labs(x = "Drive", y = "Highway MPG") +
  theme_minimal()

ggplot(data = mpg) +
  geom_density_ridges(aes(x = hwy, y = drv)) +
  labs(x = "Highway MPG", y = "Drive") +
  theme_minimal()

Linear regression with one quantitative explanatory variable

Let’s look at an example.

lm_displ <- lm(hwy ~ displ, data=mpg)

lm_displ

## 
## Call:
## lm(formula = hwy ~ displ, data = mpg)
## 
## Coefficients:
## (Intercept)        displ  
##      35.698       -3.531

Notice a couple things:

1. I have saved this model to an object named lm_displ.
   
2. The output gives us the intercept, \(\hat{\beta}_0\) and the slope or coefficient for the displ term, \(\hat{\beta}_1\).

We can obtain more output (which we’ll want when we start doing inference), using the tidy() function.

tidy(lm_displ)

For now, we are just interested in the term and estimate columns. The term column is the term/variable in the model equation. The column titled estimate gives the estimated coefficients of the variables/terms (ie. the \(\hat{\beta}\)s.

We could write out the equation above as:

\[ \hat{hwy} = 35.70 - 3.53 displ \]

aug_displ <- augment(lm_displ)

aug_displ %>% 
  slice(1:10)

augment(lm_displ, data = mpg) %>% 
  slice(1:10)

augment(lm_displ, 
        newdata = tibble(displ = c(4.5)))

#example of more than one prediction
augment(lm_displ, 
        newdata = tibble(displ = c(3, 3.5, 4.0, 4.5, 5.0, 5.5)))

#The pipe, %>%, allows us to pass data to the next function. It is read, "and then."
augment(lm_displ, data = mpg) %>% 
  ggplot() +
  geom_jitter(aes(x = displ, y = hwy)) + 
  geom_line(aes(x = displ, y = .fitted), color = "blue") + #the fitted line
  labs(x = "Engine Displacement", y = "Highway MPG") +
  theme_minimal()

Linear regression with one categorical (factor) variable

Now, let’s examine a model that uses a categorical (factor) explanatory variable.

lm_drv <- lm(hwy ~ drv,
                 data=mpg)
tidy(lm_drv)

It did something. WHAT did it do? Any ideas? First, let’s run the following code. How is what we see below related to what we see above?

mpg %>% 
  group_by(drv) %>% 
  summarize(mean_hwy = mean(hwy))

Since R cannot use drv directly, since its values are words not numbers, it is creating 2 new variables:

drvf is 1 if drv is f and 0 otherwise, and drvr is 1 if drv is r and 0 otherwise.

Let’s write out the model equation:

\[ \hat{hwy} = 19.18 + 8.99drvf + 1.83drvr \]

QUESTIONS

1. When will both drvf and drvr take a value of 0?
   
2. What is the predicted value for a front wheel drive vehicle?
   
3. What is the predicted value for a rear wheel drive vehicle?
4. What is the predicted value for a 4 wheel drive vehicle?
5. Interpret each of the coefficients.

In general, in a simple linear regression with one categorical explanatory variable, the intercept is the average response when the explanatory variable is at the “baseline” or “reference” level. That is the level of the categorical variable that does not have a dummy variable created for it. Notice, that is also the average response when all other variables are zero, same as when we had one quantitative explanatory variable.

The other coefficients are the difference in average response between the level of the indicator variable and the reference level. This is similar to “increasing the explanatory variable by one,” but it doesn’t really make sense to talk about it that way.

Explanation

Let’s return to the model output.

tidy(lm_twovars)

Given what we just discussed, try interpreting each of the estimated coefficients from the model on your own.

In general, in a multiple linear regression model, the intercept is the average response when the explanatory variables are all equal to zero. If there are categorical variables in the model, the intercept then reflects when all of those are at the baseline/reference level.

Coefficients of quantitative explanatory variables are the average change in the response variable for a one unit increase in the explanatory variable, with all other variables held fixed (this indicates that we’ve accounted for those variables being in the model).

Coefficients of categorical explanatory variables are the difference in average response between the level of the indicator variable and the reference level, holding all other variables constant.

Prediction

Just as with simple linear regression, we can use the augment() function to help us find predicted values.

To find the predicted responses for all the data used to fit this model:

augment(lm_twovars) %>% 
  slice(1:10)

This now includes all the variables from the orginal dataset, plus the “augmented” variables:

augment(lm_twovars, data = mpg) %>% 
  slice(1:10)

And, we can predict new or future values. REMINDER: se sure the variables are named EXACTLY as the variable name appears in the original dataset, and categorical/factor variables need to have their values in quotes.

augment(lm_twovars, 
        newdata = tibble(displ = c(3, 3.5, 4.0, 4.5, 5.0, 5.5),
                         drv = c("4", "f", "r", "r", "4", "f")))

We can also plot the fitted “line” on top of the original data. You can’t always plot more compled models because at some point there are too many dimensions. What is special about this model?

augment(lm_twovars, data = mpg) %>% 
  ggplot() +
  geom_jitter(aes(x = displ, y = hwy, color = drv)) + 
  geom_line(aes(x = displ, y = .fitted, color = drv)) + #the fitted lines
  labs(x = "Engine Displacement", y = "Highway MPG") +
  theme_minimal()

Constructing Hypotheses

A good scientist/researcher is skeptical of claims until they see solid evidence. If a person claims that two variables (such as engine displacement and highway gas mileage) have a relationship, a good researcher will need to see evidence before accepting a claim as true.

This is the basic framework for hypothesis testing. The hypothesis testing framework generally has two hypotheses: a null hypothesis (usually denoted \(H_0\)) and an alternative hypothesis (usually denoted \(H_A\)). We assume the null hypothesis until/unless we have sufficient evidence to reject it in favor of the alternative hypothesis.

In the context of the simple linear regression model, we assume that the two variables have no linear relationship. We then use our data to measure the evidence of the linear relationship. If we have enough evidence, we can say that the two variables have a significant linear relationship.

Writing the null and alternative hypotheses in sentences is a bit verbose, so researchers usually write the hypotheses as equations. To do this, we need to define some notation. Let’s call the predictor \(x\) and the response variable \(y\). Define \(\beta_0\) to be the intercept and \(\beta_1\) to be the slope. Then our predicted value of y is \[ \hat{y} = {\beta}_0 + {\beta}_1 x , \]

With this setup, we can now translate the null hypothesis (“there is no linear relationship between x and y”) to \(H_0: \beta_1 = 0\). This means that our regression equation is \(\hat{y} = {\beta}_0 + 0 x = \beta_1\); plug in any x value, and the regression equation will produce the same prediction.

The alternative hypothesis is that x and y do have a linear relationship, so the alternative hypothesis is \(H_A: \beta_1 \neq 0\).

Putting these hypotheses together, we have:

\[H_0: \beta_1 = 0\]

\[H_A: \beta_1 \neq 0\]

Let’s return to the model for predicting a car’s highway mileage based on the car’s engine displacement. Our model is:

\[ \widehat{mpg} = \beta_0 + \beta_1 displ. \]

Our null hypothesis is that the car’s highway mileage does not have a linear relationship with the car’s engine displacement; that is \(H_0: \beta_1 = 0\). The alternative hypothesis is that these that the car’s highway mileage does have linear relationship with the car’s engine displacement; that is \(H_A: \beta_1 \neq 0\). Until we have evidence that shows a linear relationship between the two variables, we assume that there is no linear relationship (that is, we assume that \(\beta_1 = 0\).)

Gaining Intuition

To better understand the hypotheses and to get a feel for whether the slope is likely to be nonzero, we revisit our plots:

augment(lm_displ, data = mpg) %>% 
  ggplot() +
  geom_jitter(aes(x = displ, y = hwy)) + 
  geom_line(aes(x = displ, y = .fitted), color = "blue") + 
  labs(x = "Engine Displacement", y = "Highway MPG") +
  theme_minimal()

This step is important. In a couple minutes, you will learn to crunch numbers. However, you are more than a number cruncher. You are a human and you need to know whether the results the computer gives you make sense.

Measuring Evidence Against the Null Hypothesis

Now that have a feel for our data, we are prepared to gather evidence against the null hypothesis. We use p-values to measure evidence against the null hypothesis. A smaller p-value provides more evidence against the null hypothesis.

Once you have constructed your model, calculating a p-value is easy: simply use the **tidy* command, as the relevant p-values are displayed in the last column of our table.

tidy(lm_displ)

Our tidy table has two rows: the first row contains information for \(\beta_0\) and the second row contains information for \(\beta_1\). Because we are testing whether \(\beta_1=0\), we select the p-value from the second row: 2.03897410^{-46} (that is 0.0000000000000000000000000000000000000000000002038974).

This is quite a small number, so we have a lot of evidence against the null hypothesis (that there is no linear relationship between highway gas mileage and engine displacement). In a moment, we will decide whether to reject the null hypothesis or not.

Deciding Whether to Reject

We said that smaller p-values provide more evidence against the null hypothesis. To decide how much evidence is enough to reject the null hypothesis, we need a threshold. Before crunching any numbers, an ethical researcher will select their threshold for how small is small enough to reject the null hypothesis. This threshold is called the significance level. Common thresholds include \(.05, .01,\) and \(.1\). For today, let’s use a significance level of \(.05\). In the future, you should check within your discipline and consult a statistician to select a smart significance level for your study.

If your p-value is smaller than your significance level, reject the null hypothesis in favor of the alternative: you have found a significant linear relationship between the two variables. If your p-value is larger than your significance level, do not reject the null hypothesis: there is not a significant linear relationship between the two variables.

Since our p-value is much smaller than today’s significance level of \(.05\), we reject the null hypothesis in favor of the alternative. There is a significance linear relationship between a car’s engine displacement and its highway gas mileage.

Always conclude in the context of your study. Rather than saying “We reject the null hypothesis” or “\(\beta_1 \neq 0\), write a sentence in plain English (or your language of choice), such as”There is a significance linear relationship between a car’s engine displacement and its highway gas mileage."

Tis better to have regressed and lost than never to have regressed at all.

Your Turn

1. Test whether the highway gas mileage has a linear relationship with the year of manufacture. More specifically, use highway gas mileage as the predictor and use year as the response.

2. Test whether the city gas mileage has a linear relationship with the engine displacement. Use city gas mileage as the response.

Interpreting Tests

Turning Your Skills Up

lm_alt <- lm(hwy ~ displ, data = mpg)

lm_null <- lm(hwy ~ 1, data = mpg)

anova(lm_null, lm_alt)

lm_drv <- lm(hwy ~ drv, data=mpg)
tidy(lm_drv)

anova(lm_drv, lm_null)

Your Turn Again

1. Using the anova function, test whether the highway gas mileage has a linear relationship with the year of manufacture. More specifically, use highway gas mileage as the predictor and use year as the response.

2. Using the anova function, test whether the city gas mileage differs between front-wheel drive, rear-wheel drive, and four-wheel drive cars.

Bigger is Not Necessarily Better

You now have the power to add additional predictors to your model. However, with great power comes great responsibility. Be wise and judiciously add more predictors to your model.

Remember this lesson from Jurassic Park.

To understand why we need to judiciously add predictors, let’s step away from regression and think about a remote control for a television. Imagine you’re a guest in someone’s home. Would you rather use a remote control with 100 buttons such as this remote:

Brain overload.

Or would you rather see a simple remote such as this one:

I would rather use this remote.

Why do I prefer the simple remote? I can do what I need to do. The other remotes are too complicated for me. I don’t know how to use them.

We see a similar lesson here: do not overcomplicate your model. If you add a predictor to your model and your model barely improves, what is the utility of your bulked up model? You should reconsider and pare down your model. Do not include all the predictors that you have on hand. Everything-bagels are delicious, but everything-regression-models are less tempting.

This is my face when you overfit the model.

What is the issue with incorporating many predictors? Actually, there are two main issues.

First, it’s a lot easier to talk about and work with simpler models. Think back to the section on interpreting models. The more predictors you have in your model, the longer and more complicated your interpretation is. If a simpler model will suffice, then keep it simple.

Second, there is a danger of “overfitting.” When too many predictors have been incorporated into a regression equation, we say that the model has been “overfit.” This means that the model has been tailored to resemble the data-set. What’s wrong with this? We don’t know that the model represents reality anymore.

Lizzo absolutely agrees.

How to Keep it Simple

Now that you understand the importance of controlling model complexity, we can discuss how to adhere to the KISS principle in practice. If two models work equally well, but one is simpler, then use the simpler model. To put this another way, in order to justify adding a predictor to your model, you must show the bulked-up model works significantly better. How can you make this decision? Use the p-value from the anova command.

Let’s use a concrete example. Recall your model for highway gas mileage based on whether a car has front-wheel, rear-wheel, or four-wheel drive. We may want to improve the predictive capability of this model by taking engine displacement into account. Then, the two models we are comparing are:

smallmod <- lm(hwy ~ drv, data = mpg)
bigmod <- lm(hwy ~ drv + displ, data = mpg)

Unless the bigmod is noticeably better than the smallmod, we should use the smallmod in order to keep it simple. Let’s use anova to see whether the bigmod is better or whether we should stay with the smallmod.

anova(bigmod, smallmod)

Our p-value is quite small, so the bigmod works noticeably better. We can use the model with both drive and displacement, as it is better than the model with drive only.

Your Turn

1. Is a model with both drive and displacement better than a model with just an intercept \(\beta_0\)? Use highway gas mileage as the response variable.

2. If we already have drive and displacement in the model, should we add the type of transmission? (Hint: the two models to compare will differ by a single variable.)

3. Consider adding to the model that “won” in question 2. Should we stick with this model, or could we justify adding the “type” of car (class)?

Another Approach

Using anova (as you just practiced) works well when you can use it. I didn’t tell you this before, but you can use the anova command to decide whether to add/remove predictors. That is, the predictors in one model must ALL be used in the other model. (We call such models “nested models.”)

You can’t (shouldn’t) use anova to compare, for example, the following two models:

mod1 <- lm(hwy ~ class, data = mpg)
mod2 <- lm(hwy ~ trans, data = mpg)

Why? The first model does not use all the predictors from the second model. Similarly, the second model does not use all the predictors from the first model.

To compare non-nested models, we can use AIC (Akaike Information Criteria). AIC rewards model fit and penalizes for model complexity. For AIC, smaller is better. If two models have AICs within 5 or 10 of each other, than you can make a call on which model you’d like to use. Often the context of your research will help you make this decision.

Let’s find the AIC for each model.

AIC(mod1)

## [1] 1241.573

AIC(mod2)

## [1] 1482.693

The AIC on mod1 is much smaller, which means that we should use mod1. Rather than using a model with a car’s transmission, it is better to use a model with the car’s class.

Your Turn

Create models with highway gas mileage as the response.

1. Use AIC to decide which model to use: a model with transmission only or a model with year of manufacture only.

2. Use AIC to decide which model to use: a model with transmission only or a model with transmission, class, fuel type, and city miles per gallon.

3. Create the best model possible, with the restriction that you can use only a single predictor.