Lab 10
Keaton Morehead (u1412276)
2024-06-05
Statistical Modeling
In this lab, we will be doing basic statistical modeling with data from the University of New Mexico to predict undergraduate enrollment.
Reading in the Data and Understanding it
enroll = read.csv("enrollmentForecast.csv")
str(enroll)## '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(enroll)## 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. :3345head(enroll)## 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 2192Scatterplots
Now, we will create scatterplots of the enrollment (ROLL) variable against all other variables.
Enrollment VS Year
library(ggplot2)
ggplot(enroll, aes(x = YEAR, y = ROLL)) +
geom_point(pch = 16, col = "blue") +
xlab("Year") +
ylab("Enrollment") +
theme_minimal()
Enrollment VS Unemployment
ggplot(enroll, aes(x = UNEM, y = ROLL)) +
geom_point(pch = 16, col = "pink") +
xlab("Unemployment") +
ylab("Enrollment") +
theme_minimal()
Enrollment VS Highschool Grads
ggplot(enroll, aes(x = HGRAD, y = ROLL)) +
geom_point(pch = 16, col = "green") +
xlab("Highschool Grads") +
ylab("Enrollment") +
theme_minimal()
Enrollment VS Income
ggplot(enroll, aes(x = INC, y = ROLL)) +
geom_point(pch = 16, col = "purple") +
xlab("Income") +
ylab("Enrollment") +
theme_minimal()
Building a Linear Model
We will now create a linear model, predicting enrollment as a function of unemployment and highschool graduates.
lm1 = lm(ROLL ~ UNEM + HGRAD, data = enroll)This creates a simple linear model, which we will now investigate further using the summary() and ANOVA() (analysis of variance) functions.
summary(lm1)##
## Call:
## lm(formula = ROLL ~ UNEM + HGRAD, data = enroll)
##
## 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-11anova(lm1)## 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 ' ' 1As can be seen above, the very low P value on both the UNEMP and HGRAD variables in the summary statistics suggests that these are statistically significantly correlated with enrollment, and very unlikely to be this highly correlated purely by chance.
Additinally, the very low P value on the F test for the anova() function suggests that this model does a great job of explaining the variance of enrollment based on the inputs unemplyment and highschool graduates.
The variable that is most closely related to enrollment is HGRAD due to a T value and F value greater than that of the UNEM variable.
Making a residual plot
hist(residuals(lm1))
The residuals are normally distributed, with a mean centered near zero, so there does not appear to be any immediate bias.
The predict() Function
We will now predict the enrollment at this university for an unemployment rate of 9% and a graduating highschool class of 25,000 students.
newUNEM = 9
newHGRAD = 25000
newData = data.frame(UNEM = newUNEM, HGRAD = newHGRAD)
predict(lm1, newData)## 1
## 21585.58As you can see, the predict() function using our model has predicted that an unemployment rate of 9% and a graduating class of 25,000 will have 21,585 new enrollees.
Second Model
We will now create a second model that includes per capita income.
lm2 = lm(ROLL ~ UNEM + HGRAD + INC, data = enroll)
summary(lm2)##
## Call:
## lm(formula = ROLL ~ UNEM + HGRAD + INC, data = enroll)
##
## 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-16anova(lm1, lm2)## 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 ' ' 1As can be seen by the anova(lm1) and anova(lm2) functions, the sum of squared errors on model 2 is much lower, indicating a better fitting model. We can say then, that the addition of the income variable helped the model.