This exercise will build a simple linear model of student
enrollment. The file enrollmentForecast.csv contains information used to
estimate undergraduate enrollment at the University of New Mexico. Each
row represents a single year (1 = 1961), and the variables include: the
enrollment (ROLL), the unemployment rate (UNEM), the number of high
school graduates (HGRAD) and monthly per capita income in 1961 US$
(INC).
Step one: Read in the data
enrollment <- read.csv("enrollmentForecast.csv")
Step two: Look at the data structure
head(enrollment)
## YEAR ROLL UNEM HGRAD INC
## 1 1 5501 8.1 9552 1923
## 2 2 5945 7.0 9680 1961
## 3 3 6629 7.3 9731 1979
## 4 4 7556 7.5 11666 2030
## 5 5 8716 7.0 14675 2112
## 6 6 9369 6.4 15265 2192
str(enrollment)
## 'data.frame': 29 obs. of 5 variables:
## $ YEAR : int 1 2 3 4 5 6 7 8 9 10 ...
## $ ROLL : int 5501 5945 6629 7556 8716 9369 9920 10167 11084 12504 ...
## $ UNEM : num 8.1 7 7.3 7.5 7 6.4 6.5 6.4 6.3 7.7 ...
## $ HGRAD: int 9552 9680 9731 11666 14675 15265 15484 15723 16501 16890 ...
## $ INC : int 1923 1961 1979 2030 2112 2192 2235 2351 2411 2475 ...
summary(enrollment)
## YEAR ROLL UNEM HGRAD INC
## Min. : 1 Min. : 5501 Min. : 5.700 Min. : 9552 Min. :1923
## 1st Qu.: 8 1st Qu.:10167 1st Qu.: 7.000 1st Qu.:15723 1st Qu.:2351
## Median :15 Median :14395 Median : 7.500 Median :17203 Median :2863
## Mean :15 Mean :12707 Mean : 7.717 Mean :16528 Mean :2729
## 3rd Qu.:22 3rd Qu.:14969 3rd Qu.: 8.200 3rd Qu.:18266 3rd Qu.:3127
## Max. :29 Max. :16081 Max. :10.100 Max. :19800 Max. :3345
Step three: Make scatterplots of ROLL against the other
variables
library(ggplot2)
Step four: Build a linear model using the unemployment rate (UNEM)
and number of spring high school graduates (HGRAD) to predict the fall
enrollment (ROLL)
fit1 <- lm(ROLL ~ UNEM + HGRAD, data = enrollment)
fit1
##
## Call:
## lm(formula = ROLL ~ UNEM + HGRAD, data = enrollment)
##
## Coefficients:
## (Intercept) UNEM HGRAD
## -8255.7511 698.2681 0.9423
Step five: Use the summary() and anova() functions to investigate
the model
summary(fit1)
##
## Call:
## lm(formula = ROLL ~ UNEM + HGRAD, data = enrollment)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2102.2 -861.6 -349.4 374.5 3603.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.256e+03 2.052e+03 -4.023 0.00044 ***
## UNEM 6.983e+02 2.244e+02 3.111 0.00449 **
## HGRAD 9.423e-01 8.613e-02 10.941 3.16e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1313 on 26 degrees of freedom
## Multiple R-squared: 0.8489, Adjusted R-squared: 0.8373
## F-statistic: 73.03 on 2 and 26 DF, p-value: 2.144e-11
- The median of the residuals is not close to zero
- UNEM coefficient = 698.3
- HGRAD coefficient = 0.9423
- r-squared = 0.8373
- F-statistic: 73.03
- p-value: 2.144e-11
anova(fit1)
## Analysis of Variance Table
##
## Response: ROLL
## Df Sum Sq Mean Sq F value Pr(>F)
## UNEM 1 45407767 45407767 26.349 2.366e-05 ***
## HGRAD 1 206279143 206279143 119.701 3.157e-11 ***
## Residuals 26 44805568 1723291
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
UNEM F-statistic: 26.349 and p-value = 2.366e-05
HGRAD F-statistic: 119.701 and p-value = 3.157e-11
The large F-statistic for HGRAD suggests that it might be a better
explanatory variable
Step six: Which variable is the most closely related to
enrollment?
High school graduation (HGRAD) is the most closely related variable
to enrollment, given the large F-statistic and small p-value
Step seven: Make a residual plot and check for any bias in the
model
residuals(fit1)
## 1 2 3 4 5 6
## -899.84932 191.63418 618.09763 -417.86171 -1744.03870 -1228.02117
## 7 8 9 10 11 12
## -953.20661 -861.58397 -607.84855 -531.96963 65.96365 -10.15619
## 13 14 15 16 17 18
## -124.18240 -536.67665 -2102.23858 -349.36915 -185.34696 374.51266
## 19 20 21 22 23 24
## -401.07352 -1004.00333 -256.11293 -889.62142 -474.53468 1544.58138
## 25 26 27 28 29
## 1403.20052 1246.31602 2018.95574 2510.88715 3603.54655
hist(residuals(fit1))

The basic histogram shows a fairly well distribution of
residuals
plot(fit1, which = 1)

However, plotting the residuals against the fitted values shows a
pattern
plot(fit1, which = 4)

Using Cook’s distance, there is high influence in the model, which
may bias the final model
Step eight: Use the predict() function to estimate the expected fall
enrollment, if the current year’s unemployment rate is 9% and the size
of the spring high school graduating class is 25,000 students
new_data <- data.frame(UNEM = 9, HGRAD = 25000)
predicted_enrollment <- predict(fit1, newdata = new_data)
predicted_enrollment
## 1
## 21585.58
The predicted enrollment is 21,585.58 students if UNEM is 9% and
HGRAD is 25,000 students
Step nine: Build a second model which includes per capita income
(INC).
fit2 <- lm(ROLL ~ UNEM + HGRAD + INC, data = enrollment)
fit2
##
## Call:
## lm(formula = ROLL ~ UNEM + HGRAD + INC, data = enrollment)
##
## Coefficients:
## (Intercept) UNEM HGRAD INC
## -9153.2545 450.1245 0.4065 4.2749
summary(fit2)
##
## Call:
## lm(formula = ROLL ~ UNEM + HGRAD + INC, data = enrollment)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1148.84 -489.71 -1.88 387.40 1425.75
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.153e+03 1.053e+03 -8.691 5.02e-09 ***
## UNEM 4.501e+02 1.182e+02 3.809 0.000807 ***
## HGRAD 4.065e-01 7.602e-02 5.347 1.52e-05 ***
## INC 4.275e+00 4.947e-01 8.642 5.59e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 670.4 on 25 degrees of freedom
## Multiple R-squared: 0.9621, Adjusted R-squared: 0.9576
## F-statistic: 211.5 on 3 and 25 DF, p-value: < 2.2e-16
The large F-statistic (211.5) and small p-value (< 2.2e-16)
suggests that the model is statistically significant
Step ten: Compare the two models with anova(). Does including this
variable improve the model?
anova(fit1,fit2)
## Analysis of Variance Table
##
## Model 1: ROLL ~ UNEM + HGRAD
## Model 2: ROLL ~ UNEM + HGRAD + INC
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 26 44805568
## 2 25 11237313 1 33568255 74.68 5.594e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The p-value when comparing the two models is 5.594e-09, which
suggests that including the per capita income variable significantly
improves the model’s ability to explain UNM enrollment