D7_JGreenstein
2026-10-14
PART 1
dadata <- swiss
str(dadata)## 'data.frame': 47 obs. of 6 variables:
## $ Fertility : num 80.2 83.1 92.5 85.8 76.9 76.1 83.8 92.4 82.4 82.9 ...
## $ Agriculture : num 17 45.1 39.7 36.5 43.5 35.3 70.2 67.8 53.3 45.2 ...
## $ Examination : int 15 6 5 12 17 9 16 14 12 16 ...
## $ Education : int 12 9 5 7 15 7 7 8 7 13 ...
## $ Catholic : num 9.96 84.84 93.4 33.77 5.16 ...
## $ Infant.Mortality: num 22.2 22.2 20.2 20.3 20.6 26.6 23.6 24.9 21 24.4 ...I started off like last week with cor() to get an idea of how the variables react to one another.
cor(dadata)## Fertility Agriculture Examination Education Catholic
## Fertility 1.0000000 0.35307918 -0.6458827 -0.66378886 0.4636847
## Agriculture 0.3530792 1.00000000 -0.6865422 -0.63952252 0.4010951
## Examination -0.6458827 -0.68654221 1.0000000 0.69841530 -0.5727418
## Education -0.6637889 -0.63952252 0.6984153 1.00000000 -0.1538589
## Catholic 0.4636847 0.40109505 -0.5727418 -0.15385892 1.0000000
## Infant.Mortality 0.4165560 -0.06085861 -0.1140216 -0.09932185 0.1754959
## Infant.Mortality
## Fertility 0.41655603
## Agriculture -0.06085861
## Examination -0.11402160
## Education -0.09932185
## Catholic 0.17549591
## Infant.Mortality 1.00000000cor(dadata)[,"Fertility"]## Fertility Agriculture Examination Education
## 1.0000000 0.3530792 -0.6458827 -0.6637889
## Catholic Infant.Mortality
## 0.4636847 0.4165560Then I made my first linear regression model. I started off with using all the possible explanatory variables (a full model).
I want to pay particular attention to the Adjusted R^2 value. It tells me how much of the variance I can expected to be accounted for by the model. So, a larger value the better in this case.
The p-values also help signal which variable(s) are statistically significant. Most of the time, I will use 0.05 or the alpha as the benchmark on whether or not I should keep the H0 of that there’s no linear regression between the 2 variables.
Lastly there’s the coefficients, that tell in what direction the relation of each variable to Fertility is, as well as its strength.
In the lm() function, I used the ‘.’ to indicate I want all the available explanatory variables included.
fullmodel <- lm(Fertility ~ ., data = dadata) # the . tells R to use all other variables as explanatory
summary(fullmodel)##
## Call:
## lm(formula = Fertility ~ ., data = dadata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.2743 -5.2617 0.5032 4.1198 15.3213
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 66.91518 10.70604 6.250 1.91e-07 ***
## Agriculture -0.17211 0.07030 -2.448 0.01873 *
## Examination -0.25801 0.25388 -1.016 0.31546
## Education -0.87094 0.18303 -4.758 2.43e-05 ***
## Catholic 0.10412 0.03526 2.953 0.00519 **
## Infant.Mortality 1.07705 0.38172 2.822 0.00734 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.165 on 41 degrees of freedom
## Multiple R-squared: 0.7067, Adjusted R-squared: 0.671
## F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10Since I started off with a full model, I will use backward elimination, or the process of removing 1 variable at a time from the linear regression model. A common way of going about it is to look at the p-values and remove the largest non-significant variable 1st. Here that would be Examination. So in the new model I just added the ‘- examination’ to the same command from the full model.
Now this is the only non-significant p-value and in fact it made the adjusted R^2 value decrease ever so slightly. Normally I want as big of a R^2 value as possible, but since this decrease is extremely small and the fact that it reduces the model’s complexity I will use this model.
fullmodel_minus_examination <- lm(Fertility ~ . - Examination, data = dadata)
summary(fullmodel_minus_examination)##
## Call:
## lm(formula = Fertility ~ . - Examination, data = dadata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.6765 -6.0522 0.7514 3.1664 16.1422
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 62.10131 9.60489 6.466 8.49e-08 ***
## Agriculture -0.15462 0.06819 -2.267 0.02857 *
## Education -0.98026 0.14814 -6.617 5.14e-08 ***
## Catholic 0.12467 0.02889 4.315 9.50e-05 ***
## Infant.Mortality 1.07844 0.38187 2.824 0.00722 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.168 on 42 degrees of freedom
## Multiple R-squared: 0.6993, Adjusted R-squared: 0.6707
## F-statistic: 24.42 on 4 and 42 DF, p-value: 1.717e-10Now that I have my model, I need to plot it to ensure that the linear regression conditions are still reasonably met.
The 4 plots show that overall the model does indeed meet all the conditions. Yes there is some variance and some extreme values, in particular the Q-Q and Scale-Location plots but nothing too drastic that I would have to remodel. The residuals vs leverage does show me that there are some high leverage observations, but they are well fitted to the model. So overall the model works solid for multiple linear regression with some minor variance.
par(mfrow = c(2,2))
plot(fullmodel_minus_examination)
I also used the Akaike Info Criterion (AIC) to run another test to further backup my decision of using the 2nd model rather than the full model. The AIC takes into account a model’s complexity as well as how well it fits the data. Unlike the adjusted R^2 I want smaller AIC values here. And again the difference is extremely small.
AIC(fullmodel, fullmodel_minus_examination)## df AIC
## fullmodel 7 326.0716
## fullmodel_minus_examination 6 325.2408Part 2
I took this task as to predict whether or not an observation’s Fertility value exceeds a value of 70.
So I made a new integer variable to use, where it returns 1 if its value is 71 and up and 0 for 70 and below.
dadata_Fertility_70 <- as.integer(dadata$Fertility > 70)
table(dadata_Fertility_70)## dadata_Fertility_70
## 0 1
## 23 24I then used the glm() function to make the logistic models. Very similar to the same setup as for the linear regression models, but it needed a couple additional arguments. I decided to use the same two set of variables as I did for the linear regression to see how the change of model impacts the returns. I also made sure to have the dependent variable as my new varaible instead of the original one.
logistic_fullmodel <- glm(dadata_Fertility_70 ~ Agriculture + Examination + Education + Catholic + Infant.Mortality,
data = dadata, family = binomial(link = "logit"))
logistic_model_exam <- glm(dadata_Fertility_70 ~ Agriculture + Education + Catholic + Infant.Mortality,
data = dadata, family = binomial(link = "logit"))
summary(logistic_fullmodel) ##
## Call:
## glm(formula = dadata_Fertility_70 ~ Agriculture + Examination +
## Education + Catholic + Infant.Mortality, family = binomial(link = "logit"),
## data = dadata)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.82826 5.25607 0.919 0.3583
## Agriculture -0.09615 0.04011 -2.397 0.0165 *
## Examination -0.32116 0.13844 -2.320 0.0203 *
## Education -0.12078 0.08610 -1.403 0.1607
## Catholic 0.02078 0.01376 1.509 0.1312
## Infant.Mortality 0.29078 0.21051 1.381 0.1672
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 65.135 on 46 degrees of freedom
## Residual deviance: 32.887 on 41 degrees of freedom
## AIC: 44.887
##
## Number of Fisher Scoring iterations: 6summary(logistic_model_exam)##
## Call:
## glm(formula = dadata_Fertility_70 ~ Agriculture + Education +
## Catholic + Infant.Mortality, family = binomial(link = "logit"),
## data = dadata)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.94524 4.13457 -0.229 0.8192
## Agriculture -0.04973 0.02675 -1.859 0.0631 .
## Education -0.21212 0.08825 -2.404 0.0162 *
## Catholic 0.03316 0.01187 2.795 0.0052 **
## Infant.Mortality 0.22203 0.17261 1.286 0.1983
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 65.135 on 46 degrees of freedom
## Residual deviance: 40.634 on 42 degrees of freedom
## AIC: 50.634
##
## Number of Fisher Scoring iterations: 5Since the logistic regression doesn’t have the adjusted R^2 value, I made sure to use the AIC again to compare the models.
Here Examination seems more significant than in the linear regression. AIC values increase when removing examination by a much larger margin than before. In fact I should try removing infant mortality instead and add exam back in since the infant mortality p-value was larger.
AIC(logistic_fullmodel, logistic_model_exam)## df AIC
## logistic_fullmodel 6 44.88671
## logistic_model_exam 5 50.63442The removal of infant mortality instead of exam worked much better, but still increased the AIC value. But like in linear regression, the difference was relatively small so I will try removing the next largest p valued variable (education).
logistic_model_infant <- glm(dadata_Fertility_70 ~ Agriculture + Examination + Education + Catholic,
data = dadata, family = binomial(link = "logit"))
summary(logistic_model_infant)##
## Call:
## glm(formula = dadata_Fertility_70 ~ Agriculture + Examination +
## Education + Catholic, family = binomial(link = "logit"),
## data = dadata)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 10.44664 3.79705 2.751 0.00594 **
## Agriculture -0.09573 0.03687 -2.596 0.00942 **
## Examination -0.29567 0.12821 -2.306 0.02110 *
## Education -0.14430 0.08420 -1.714 0.08657 .
## Catholic 0.02251 0.01278 1.762 0.07802 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 65.135 on 46 degrees of freedom
## Residual deviance: 35.149 on 42 degrees of freedom
## AIC: 45.149
##
## Number of Fisher Scoring iterations: 6AIC(logistic_fullmodel, logistic_model_exam, logistic_model_infant)## df AIC
## logistic_fullmodel 6 44.88671
## logistic_model_exam 5 50.63442
## logistic_model_infant 5 45.14880That clearly didn’t help the model. So in the end the full model seemed to work the best.
logistic_model_infant_edu <- glm(dadata_Fertility_70 ~ Agriculture + Examination + Catholic,
data = dadata, family = binomial(link = "logit"))
summary(logistic_model_infant_edu)##
## Call:
## glm(formula = dadata_Fertility_70 ~ Agriculture + Examination +
## Catholic, family = binomial(link = "logit"), data = dadata)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.00727 3.40544 2.645 0.00817 **
## Agriculture -0.07337 0.03275 -2.240 0.02506 *
## Examination -0.34698 0.11643 -2.980 0.00288 **
## Catholic 0.01365 0.01121 1.217 0.22352
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 65.135 on 46 degrees of freedom
## Residual deviance: 38.957 on 43 degrees of freedom
## AIC: 46.957
##
## Number of Fisher Scoring iterations: 5AIC(logistic_fullmodel, logistic_model_exam, logistic_model_infant, logistic_model_infant_edu)## df AIC
## logistic_fullmodel 6 44.88671
## logistic_model_exam 5 50.63442
## logistic_model_infant 5 45.14880
## logistic_model_infant_edu 4 46.95666Then I plotted the full model.
At first the plots seemed way off but after remembering that logistic is a different process form linear, the changes are mostly expected. Like the curves in the Residuals vs Fitted for the distribution and the Scale-Location. Overall the plots don’t show any major concern over the model. There are a few particular observations (boundry, entremont and moudon) that may warrant some further investigation.
plot(logistic_fullmodel)
