PSCI 9590A - Introduction to Quantitative Methods

2.1.1 Using the lm() function

Using lm(), we specify a simple linear regression formula: dependent variable ~ independent variable. Here, we are regressing one dependent variable on one independent variable. The lm() function will generate an intercept and a regression coefficient. The intercept represents the value of the dependent variable when the independent variable is equal to zero. The regression coefficient represents the mean change in the dependent variable for one unit of change in the independent variable (in other words, the slope).

Note that lm() automatically performs listwise deletion of rows with a missing value within the vectors you are using. In other words, it only uses pairwise complete observations.

 # Creating our model using lm()
model1 <- lm(Fertility ~ Infant.Mortality, data=swiss)

# Viewing our regression results using summary()
summary(model1)

## 
## Call:
## lm(formula = Fertility ~ Infant.Mortality, data = swiss)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.672  -5.687  -0.381   7.239  28.565 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)       34.5155    11.7113   2.947  0.00507 **
## Infant.Mortality   1.7865     0.5812   3.074  0.00359 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.48 on 45 degrees of freedom
## Multiple R-squared:  0.1735, Adjusted R-squared:  0.1552 
## F-statistic: 9.448 on 1 and 45 DF,  p-value: 0.003585

How do you interpret the coefficients, the p-values and the r-squared here? Note that the adjusted r-squared is adjusted to the number of predictors in the model, and increases when adding a new variable increases the percentage of explained variance more than would be expected by chance alone.

Let’s now visualize the results of our regression model model1 in a graph.

# Plotting the relationship between the IV and the DV using a scatterplot
ggplot(data=swiss, aes(Infant.Mortality,Fertility)) +
  geom_point() + xlab("Infant Mortality (%)") + ylab("Fertility Index") +
  theme_minimal(base_size = 20) + 
  geom_smooth(method='lm', formula= y~x) + # This is the code that adds the lm line
  xlim(0,30)

2.1.2 Polynomial models

### Loess ###

# Plotting the relationship between the IV and the DV using a scatterplot
ggplot(data=swiss, aes(Infant.Mortality,Fertility)) +
  geom_point() + xlab("Infant Mortality (%)") + ylab("Fertility Index") +
  theme_minimal(base_size = 20) + 
  geom_smooth(method='loess') +  # This is the code that adds the loess line
  xlim(0,30)

## `geom_smooth()` using formula 'y ~ x'

### Polynomial Model ###

# Plotting the relationship between the IV and the DV using a scatterplot
ggplot(data=swiss, aes(Infant.Mortality,Fertility)) +
  geom_point() + xlab("Infant Mortality (%)") + ylab("Fertility Index") +
  theme_minimal(base_size = 20) + 
  geom_smooth(method='lm', formula=y~poly(x,2)) +  # This is the code that adds the polynomial line
  xlim(0,30)

Let’s test drive all of this!

4.1 Linearity of the residuals

We want to make sure that our bivariate relationship is linear, which implies that our residuals are also linear. The red line indicates the relationship between the residuals and the fitted values (i.e. your predicted values). Ideally, the line should be horizontal at zero. Why? Because you should have on average 0 residual if all errors are randomly distributed (just like 100 people who sing badly might have a perfect pitch on average when singing together).

# The plot() function will give you up to 6 diagnostics tests on your model
plot(model1,1) # By indicating "1" as an argument, we select this specific test

Here, the line is indeed horizontal at zero, so we can confirm the linearity of the relationship between infant mortality and fertility rates. If the relationship had not been linear, we could have used a non-linear model or a polynomial function—but we won’t cover this in this class.

4.2 Homoscedasticity

We now want to make sure that the variance of our residuals is similar across the range of predicted values. This would mean that our model always misses the prediction by roughly the same range, across our observations.

The red line on the graph below indicates the relationship between the square root of standardized residuals, i.e. the square root of \[ \frac{\text{observed value - predicted value}}{\sqrt{\text{predicted value}}} \] and your fitted values (i.e. your prediction). Why do we use a square root on the standardized residuals? Because high values get compressed and low values become more spread out. This is not critical—in case you were curious about this, just know that we take the square root because we want the data to be less skewed.

Note that the average of the standardized residuals is 0 and their standard deviation is 1 (like any standard normal distribution, which is a specific type of normal distribution with mean = 0 and sd = 1). Ideally, the red line should be horizontal with equally (randomly) spread points, since this would represent the same amount of standardized residuals regardless of our fitted value.

# The plot() function will give you up to 6 diagnostics tests on your model
plot(model1,3) # By indicating "3" as an argument, we select this specific test

Here, our line is roughly horizontal, but we observe that our variance seems to be slightly higher for provinces with fitted values between 63 and 70 for the fertility rate (although they remain random, just within a broader range, hence the straight red line). What do you think?

4.3 Normality of residuals

We want to make sure that our residuals are distributed normally, which is one of the assumptions of OLS. The y axis is here again the standardized residuals (but not squared, this time). The x axis is the theoretical quantile associated with any standard normal distribution. The one thing that is critical to understand here is this: the diagonal dotted line represents what the standardized residuals should look like if they were truly normally distributed.

# The plot() function will give you up to 6 diagnostics tests on your model
plot(model1,2) # By indicating "2" as an argument, we select this specific test

# We can also visualize this like we did last time... 
# See how the small jump on the line in the previous graph fits
# with the patterns in the graph below? 
plot(density(model1$residuals))

Here, our observations are roughly normally distributed. If they weren’t, we might have wanted to include missing covariates that could explain that our predictor misses out on certain predictions (ex.: if we had disproportionally high residuals for all provinces with high fertility rates, and all these provinces with high fertility rates were located at high altitudes, we might have wanted to control for altitude to remove this bias in case this is what is causing the distortion).

4.4 Outliers and high leverage points (extra material)

Some of your data points might be outliers (i.e. extreme values on the dependent variable) or have high leverage (i.e. extreme values on the independent variable). In the test below, we want to make sure that our outliers are within 3 standard deviations from the mean standardized residuals value (which, by definition, is 0). We also want to see whether standardized residuals have a leverage smaller than \(2(p + 1)/n\) (where p is the number of predictors, in this case 1, and n is the number of observations, in this case 47, so \(2(p + 1)/n\) = \(2(1 + 1)/47\) = \(4/47\) = \(0.08510638\)).

# The plot() function will give you up to 6 diagnostics tests on your model
plot(model1,5) # By indicating "5" as an argument, we select this specific test

Here, none of the observations are more distant than 3 standard deviations from the mean, so we don’t need to remove outliers. For a few of our observations, however, leverage is above 0.08.

# Printing out the rows with a leverage above 0.08
swiss_withLmInfo[swiss_withLmInfo$.hat>0.08,]

## # A tibble: 4 × 9
##   .rownames  Fertility Infant.Mor…¹ .fitted .resid   .hat .sigma .cooksd .std.…²
##   <chr>          <dbl>        <dbl>   <dbl>  <dbl>  <dbl>  <dbl>   <dbl>   <dbl>
## 1 Porrentruy      76.1         26.6    82.0 -5.94  0.135    11.6 2.41e-2 -0.556 
## 2 Glane           92.4         24.9    79.0 13.4   0.0843   11.4 6.84e-2  1.22  
## 3 La Vallee       54.3         10.8    53.8  0.490 0.235    11.6 3.68e-4  0.0489
## 4 Conthey         75.5         15.1    61.5 14.0   0.0814   11.4 7.18e-2  1.27  
## # … with abbreviated variable names ¹​Infant.Mortality, ²​.std.resid

4.5 Influential values (extra material)

This graph gives us the Cook’s distance for each observation, which indicates the influence that each data point has in the model. How is it calculated? Essentially by removing one given data point from the model and recalculating the regression. It then tells you how much the values in the model have changed when this one particular data point was excluded. If the model changes a lot from removing simply one observation, this one data point had a lot of influence!

There isn’t a clear cutoff here. Some people say that a Cook’s distance of 0.5 or even 1 is too much, other say that you should simply remove observations that are 4/n (where n is the total number of observations, in this case 47, so 4/47 = 0.08), or simply those that over the very distant from the rest.

# The plot() function will give you up to 6 diagnostics tests on your model
plot(model1,4) # By indicating "4" as an argument, we select this specific test

# Checking the province for observation #8
rownames(swiss)[8]

## [1] "Glane"

Here, we might want to remove provinces Glane, Conthey, Sierre and Ville de Genève because they are above the 0.08 threshold.