Post your solutions, interpretation, and code for your peers to see.

MULTIVARIATE REGRESSION

Using the “swiss” dataset, build the best multiple regression model you can for the variable Fertility.

data("swiss")
head(swiss)
##              Fertility Agriculture Examination Education Catholic
## Courtelary        80.2        17.0          15        12     9.96
## Delemont          83.1        45.1           6         9    84.84
## Franches-Mnt      92.5        39.7           5         5    93.40
## Moutier           85.8        36.5          12         7    33.77
## Neuveville        76.9        43.5          17        15     5.16
## Porrentruy        76.1        35.3           9         7    90.57
##              Infant.Mortality
## Courtelary               22.2
## Delemont                 22.2
## Franches-Mnt             20.2
## Moutier                  20.3
## Neuveville               20.6
## Porrentruy               26.6
#Create a multiple regression model
model<-lm(Fertility~Agriculture+Examination+Education+Catholic+Infant.Mortality,data = swiss)
summary(model)
## 
## Call:
## lm(formula = Fertility ~ Agriculture + Examination + Education + 
##     Catholic + Infant.Mortality, data = swiss)
## 
## 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-10
par(mfrow=c(2,2))
plot(model)

#AIC is a statistical metric to compare, rank, and select best-fit statistical model among a set of alternatives. Lower AIC = better fit/less complex.
#Use step AIC function to find a model that minimizes AIC 
AIC_model<-step(model,direction = "both",trace = FALSE)
summary(AIC_model)
## 
## Call:
## lm(formula = Fertility ~ Agriculture + Education + Catholic + 
##     Infant.Mortality, data = swiss)
## 
## 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-10
formula(AIC_model)
## Fertility ~ Agriculture + Education + Catholic + Infant.Mortality
par(mfrow=c(2,2))
plot(AIC_model)

AIC(model)
## [1] 326.0716
AIC(AIC_model)
## [1] 325.2408
summary(model)$r.squared
## [1] 0.706735
summary(AIC_model)$r.squared
## [1] 0.6993476

The stepwise model retained 4/5 predictors; it excluded the “examination” variable. Excluding the examination variable improves the AIC (from 326 to 325). It is not surprising that the step function removed the “examination” variable since that was the only variable with p> 0.05.

After removing the examination variable, the p-values for the remaining predictors decrease for all variables except agriculture. This suggests stronger statistical evidence for these coefficients in the four-variable model.

LOGISTIC REGRESSION

Then build a logistic regression model for predicting Fertility>70.0.

#Create binary outcome variable, 1 if Fertility is >70, 0 if Fertility is less than or equal to 70
swiss$HighFertility<-as.integer(swiss$Fertility>70)

#Number of observations in each category
table(swiss$HighFertility)
## 
##  0  1 
## 23 24

The outcome categories are evenly balanced: There are 24 districts with fertility >70 and 23 districts with fertility.

A log model estimates the log-odds of “high fertility” (>70) as a function of the 5 explanatory variables.

#Fit model
log_model<-glm(HighFertility~Agriculture+Examination+Education+Catholic+Infant.Mortality,data = swiss,family = binomial(link = "logit"))
summary(log_model)
## 
## Call:
## glm(formula = HighFertility ~ Agriculture + Examination + Education + 
##     Catholic + Infant.Mortality, family = binomial(link = "logit"), 
##     data = swiss)
## 
## 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: 6

At 5% significance level, Agriculture and Examination have statistically significant negative coefficients. This means that, holding the other variables constant, higher values of agriculture and examination are associated with lower odds of fertility>70.

#Predicted probabilities
predicted_prob<-predict(log_model,type="response")
#Plots
library(ggplot2)
ggplot(swiss,aes(x=Agriculture,y=predicted_prob))+geom_point( color="green")+geom_smooth(method = "glm",method.args=list(family="binomial"),color="darkgreen",se=FALSE)+labs(title = "Logistic Regression: Agriculture", x= "% of males in agriculture",y="Probability Fertility > 70")
## `geom_smooth()` using formula = 'y ~ x'
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

ggplot(swiss,aes(x=Examination,y=predicted_prob))+geom_point( color="blue")+geom_smooth(method = "glm",method.args=list(family="binomial"),color="darkblue",se=FALSE)+labs(title = "Logistic Regression: Examination", x= "% draftees receiving highest marks on army exam",y="Probability Fertility > 70")
## `geom_smooth()` using formula = 'y ~ x'
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

The two models (multiple linear regression and logistic regression) address different questions: one predicts a numerical fertility measure, while the other estimates the probability of belonging to the high-fertility category.