roll = read.csv("enrollmentForecast.csv")
library(ggplot2)
Look at the data structure
str(roll)
## '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 ...
names(roll)
## [1] "YEAR" "ROLL" "UNEM" "HGRAD" "INC"
Make scatterplots of ROLL against the other variables
ggplot(roll, aes(x=ROLL, y=UNEM)) + geom_point()
ggplot(roll, aes(x=ROLL, y=HGRAD)) + geom_point()
ggplot(roll, aes(x=ROLL, y=INC)) + geom_point()
Build a linear model using the unemployment rate (UNEM) and the number of spring high school graduates (HGRAD) to predict fall enrollment (ROLL)
ex1 = lm(ROLL ~ UNEM + HGRAD, data = roll)
Use the summary() and anova() functions to investigate the model. Which variable is the most closely related to enrollment?
summary(ex1)
##
## Call:
## lm(formula = ROLL ~ UNEM + HGRAD, data = roll)
##
## 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(ex1)
## 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
The variable most closely related to enrollment is high school graduates based on the lower p value.
Make a residual plot and check for any bias in the model
hist(residuals(ex1))
plot(ex1, which=1)
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
ex2 = data.frame(UNEM = 9, HGRAD = 25000)
predict(ex1, ex2)
## 1
## 21585.58
Build a second model which includes per capita income (INC)
ex3 = lm(ROLL ~ UNEM + HGRAD + INC, data = roll)
Compare the two models with anova(). Does including this variable improve the model?
anova(ex1)
## 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
anova(ex3)
## Analysis of Variance Table
##
## Response: ROLL
## Df Sum Sq Mean Sq F value Pr(>F)
## UNEM 1 45407767 45407767 101.02 2.894e-10 ***
## HGRAD 1 206279143 206279143 458.92 < 2.2e-16 ***
## INC 1 33568255 33568255 74.68 5.594e-09 ***
## Residuals 25 11237313 449493
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Including this variable does seem to improve the model, based on the p values.