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)

Enrollment versus unemployment rate - not obviously correlated

ggplot(enrollment, aes(x = ROLL, y = UNEM, col = YEAR)) + 
  geom_point()

Enrollment versus spring high school graduates in NM - maybe correlated

ggplot(enrollment, aes(x = ROLL, y = HGRAD, col = YEAR)) + 
  geom_point()

Enrollment versus per capita income in ABQ - maybe correlated

ggplot(enrollment, aes(x = ROLL, y = INC, col = YEAR)) + 
  geom_point()

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
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 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