Homework 6 — Linear and Logistic Regression
CGilbert A
This is the last homework. Part 1 uses linear regression on country-level data. Part 2 uses logistic regression on a medical dataset.
Part 1 — Linear Regression: AllCountries dataset
Download AllCountries.csv from the Datasets folder on
Blackboard. The dataset has 217 countries with variables including GDP,
LifeExpectancy, Health, Internet, CO2, Energy, Electricity, and
more.
countries <- read.csv("AllCountries.csv")
dim(countries) # 217 countries and 26 variables## [1] 217 26head(countries)## Country Code LandArea Population Density GDP Rural CO2 PumpPrice
## 1 Afghanistan AFG 652.86 37.172 56.9 521 74.5 0.29 0.70
## 2 Albania ALB 27.40 2.866 104.6 5254 39.7 1.98 1.36
## 3 Algeria DZA 2381.74 42.228 17.7 4279 27.4 3.74 0.28
## 4 American Samoa ASM 0.20 0.055 277.3 NA 12.8 NA NA
## 5 Andorra AND 0.47 0.077 163.8 42030 11.9 5.83 NA
## 6 Angola AGO 1246.70 30.810 24.7 3432 34.5 1.29 0.97
## Military Health ArmedForces Internet Cell HIV Hunger Diabetes BirthRate
## 1 3.72 2.01 323 11.4 67.4 NA 30.3 9.6 32.5
## 2 4.08 9.51 9 71.8 123.7 0.1 5.5 10.1 11.7
## 3 13.81 10.73 317 47.7 111.0 0.1 4.7 6.7 22.3
## 4 NA NA NA NA NA NA NA NA NA
## 5 NA 14.02 NA 98.9 104.4 NA NA 8.0 NA
## 6 9.40 5.43 117 14.3 44.7 1.9 23.9 3.9 41.3
## DeathRate ElderlyPop LifeExpectancy FemaleLabor Unemployment Energy
## 1 6.6 2.6 64.0 50.3 1.5 NA
## 2 7.5 13.6 78.5 55.9 13.9 808
## 3 4.8 6.4 76.3 16.4 12.1 1328
## 4 NA NA NA NA NA NA
## 5 NA NA NA NA NA NA
## 6 8.4 2.5 61.8 76.4 7.3 545
## Electricity Developed
## 1 NA NA
## 2 2309 1
## 3 1363 1
## 4 NA NA
## 5 NA NA
## 6 312 1Q1 — Simple Linear Regression
Fit a simple linear regression model predicting
LifeExpectancy from GDP.
# Your code:
model <- lm(LifeExpectancy ~ GDP, data = countries)
coef(model) # intercept and slope## (Intercept) GDP
## 6.842208e+01 2.476441e-04summary(model)$r.squared # R-squared## [1] 0.4303976summary(model)$adj.r.squared # adjusted R-squared## [1] 0.4271795summary(model)$coefficients # includes the extra p-value column (Pr(>|t|))## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.842208e+01 5.414692e-01 126.36376 2.062763e-175
## GDP 2.476441e-04 2.141373e-05 11.56473 2.115234e-23Report the intercept and slope. What does the slope mean in plain English (e.g., “for every X increase in GDP, life expectancy increases by Y”)?
The slope is 2.476441e-04 , approximately 0.00024764. For every 1 dollar increase in GDP, the predicted life expectancy increases by about 0.00024764 years or 2.17 hours. The scale on GDP is much larger than dollars, so this slope could be impactful / meaningful.
The y-intercept is 6.842208e+01, approximately 68.422. The y-intercept is the predicted life expectancy for a country with a GDP of x = 0 dollars. That is, a country with GDP of 0 dollars has a predicted life expectancy of approximately 68.4 years. The y-intercept of this model may not be realistic-perhaps extrapolation outside of data domain.
# years(365 days/ 1 year)(24 hours/ 1 day) = hours
years_to_hours <- 0.00024764*(365)*(24)
years_to_hours## [1] 2.169326# so about 2.17 hoursWhat does the R² value tell you about how well GDP explains life expectancy?
The r-squared value is 0.4303976. About 43% of the variation in life expectancy across countries is explained by GDP alone. The remaining 57% of variation could be from other factors that we were not included in this model.
Q2 — Multiple Linear Regression
Fit a multiple regression predicting LifeExpectancy from
GDP, Health, and Internet.
# Your code:
model_3 <- lm(LifeExpectancy ~ GDP + Health + Internet, data = countries)
coef(model_3) # intercept and slope## (Intercept) GDP Health Internet
## 5.908027e+01 2.367169e-05 2.478764e-01 1.903116e-01summary(model_3)$r.squared # R-squared## [1] 0.7213439summary(model_3)$coefficients # includes the extra p-value column (Pr(>|t|))## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.908027e+01 8.149073e-01 72.499370 3.109710e-129
## GDP 2.367169e-05 2.286536e-05 1.035264 3.020246e-01
## Health 2.478764e-01 6.618957e-02 3.744946 2.471736e-04
## Internet 1.903116e-01 1.656372e-02 11.489669 6.014217e-23model_3##
## Call:
## lm(formula = LifeExpectancy ~ GDP + Health + Internet, data = countries)
##
## Coefficients:
## (Intercept) GDP Health Internet
## 5.908e+01 2.367e-05 2.479e-01 1.903e-01Interpret the coefficient on Health (controlling
for GDP and Internet).
# years(365 days/ 1 year) = days
years_to_days_2 <- 0.248*(365)
years_to_days_2## [1] 90.52Controlling for GDP and Internet access, a 1 unit increase in Health is associated with an increase of about 0.248 years (or 90 days) in predicted life expectancy.
How does the adjusted R² compare to the simple model in Q1? What does that suggest about adding predictors?
For this multiple linear regression model, the adjusted R² value is 0.7213. The adjusted R² value for the simple linear model was 0.4272. The simple linear model was able to explain about 43% of the variation in life expectancy. The multiple linear regression model was able to explain about 72% of the variation in life expectancy.
This suggests that adding Health and Internet as predictors substantially improved the model’s adjusted r-squared value.
The p-value for Health coefficient was p = 0.000247. The p-value for Internet coefficient was p = 6.01e-23. These coefficients appeared significant in the multlinear model. While the p-value for GDP went down to p = 0.302. Perhaps Health and Internet have significant association with life expectancy.
–
Q3 — Checking Assumptions
For the simple model in Q1 (LifeExpectancy ~ GDP):
- Briefly describe what you would CHECK to evaluate homoscedasticity and normality of residuals. What would an ideal outcome look like?
“The residual plot shows random scatter around 0, which suggests the homoscedasticity (equal variance) assumption is reasonably satisfied.”
To evaluate homoscedasticity, we would evaluate the residual plot for residuals that are ideally randomly scattered above and below the zero for the whole range of values. Any funnel, curve, or cone shaped residuals would suggest that the homoscedasticity assumption was not satisfied.
“The points on the Q-Q plot closely follow the diagonal line, without much deviation, suggesting residuals are roughly normal.”
To check for normality of residuals, we would evaluate the Q-Q plot of residuals for points falling closely along the diagonal reference line, especially in the tails. If we see curving at the extremes then the residuals may have thicker tails than a normal distribution would.
- Then code your check (residual plot + Q-Q plot of residuals) and reflect on what you see.
# Your code:
plot(model, which = 1) # residuals vs fitted -- homoscedasticity
plot(model, which = 2) # Q-Q plot -- normality
Your reflection:
Looking at the residual plot, we see the residuals forming a curved pattern. “The residual plot shows a curved pattern, which suggests the homoscedasticity assumption is not reasonably satisfied.”
At around GDP of 67 to 70, where a concentration of GDP’s are located, there is a range of residuals from less than -15 to almost 7 or 8 . From 71 or 73, the residuals start to spread out and fall as GDP goes out to 90.
Looking at the Q-Q Plot of Residuals, the middle of the diagonal is closely followed by the residuals. But in the lower left part of the diagonal, the left tail of the data bends under the diagonal and in the upper diagonal the right tail goes above the diagonal. The residuals seem skewed.
“The points deviate from the diagonal line, with noticeable curving at the tails, suggesting residuals are not very normal.”
Q4 — Diagnosing Fit (RMSE)
For the multiple regression in Q2, calculate the RMSE (root mean squared error).
# Hint: sqrt(mean(residuals(model)^2))
rmse_model_3 <- sqrt(mean(residuals(model_3)^2))
rmse_model_3## [1] 4.056417What does the RMSE represent in the context of predicting life expectancy? How would large residuals for certain countries affect your confidence in the model?
The RMSE is 4.056417. In the context of predicting life expectancy, this means that, on average, our model’s predicted life expectancy differs from a country’s actual life expectancy by about 4.1 years. (Using the template)
I am confident in the model for those countries that it fits well.
Countries with unusually large residuals, like rows 171, 38, and 39 that showed up repeatedly in the diagnostic plots, have errors far bigger than the average RMSE of about 4.1 years.
If there was a way to see the similarities in these countries that are associated with larger than average errors (region, income level, data quality issues, or some other unaccounted-for factor), then we could build better models for those countries and then use the current model for the countries that it has approximated well.
Q5 — Multicollinearity (no code)
Suppose you fit a regression predicting CO2 using both
Energy and Electricity. These two predictors
are highly correlated.
Explain in 2-3 sentences how this multicollinearity could affect (a) the interpretation of the coefficients and (b) the reliability of the model.
Since Energy and Electricity are highly correlated ( electricty even being a type of energy , almost a subset of energy), they don’t vary independently of each other in the data. This makes it difficult to isolate “the effect of Energy, holding Electricity constant,” because in practice, when one changes, the other tends to change too. This means the model can’t cleanly attribute CO2’s variation to one predictor versus the other, so the individual coefficients can become sensitive to changes to the data or model and therefore hard to act upon.
This instability also makes the model less reliable because standard errors can get inflated/deflated and individual p-values can change.
It becomes tougher to see which variables are changing / influencing the response variable of CO2.
The model’s overall predictive accuracy (like R²) may still look fine, but the individual coefficients themselves are tough to trust for interpretation.
Part 2 — Logistic Regression: Pima Indians Diabetes
This part uses the Pima Indians Diabetes dataset (768 patients, binary outcome: 0 = no diabetes, 1 = diabetes).
Don’t change this chunk — it loads and cleans the data:
url <- "https://raw.githubusercontent.com/jbrownlee/Datasets/master/pima-indians-diabetes.data.csv"
data <- read.csv(url, header = FALSE)
colnames(data) <- c("Pregnancies", "Glucose", "BloodPressure", "SkinThickness",
"Insulin", "BMI", "DiabetesPedigreeFunction", "Age", "Outcome")
data$Outcome <- as.factor(data$Outcome)
# Replace impossible 0 values with NA
data$Glucose[data$Glucose == 0] <- NA
data$BloodPressure[data$BloodPressure == 0] <- NA
data$BMI[data$BMI == 0] <- NA
colSums(is.na(data))## Pregnancies Glucose BloodPressure
## 0 5 35
## SkinThickness Insulin BMI
## 0 0 11
## DiabetesPedigreeFunction Age Outcome
## 0 0 0Q6 — Fit the Logistic Model
Fit a logistic regression predicting Outcome from
Glucose, BMI, and Age.
# Hint: glm(Outcome ~ Glucose + BMI + Age, data = data, family = "binomial")
model_q6 <- glm(Outcome ~ Glucose + BMI + Age, data = data, family = "binomial")
summary(model_q6)##
## Call:
## glm(formula = Outcome ~ Glucose + BMI + Age, family = "binomial",
## data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.032377 0.711037 -12.703 < 2e-16 ***
## Glucose 0.035548 0.003481 10.212 < 2e-16 ***
## BMI 0.089753 0.014377 6.243 4.3e-10 ***
## Age 0.028699 0.007809 3.675 0.000238 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 974.75 on 751 degrees of freedom
## Residual deviance: 724.96 on 748 degrees of freedom
## (16 observations deleted due to missingness)
## AIC: 732.96
##
## Number of Fisher Scoring iterations: 4Get the summary of the model. For each predictor, does an increase RAISE or LOWER the odds of diabetes? Which predictors are significant (p < 0.05)?
All three coefficients are positive. Therefore higher Glucose higher BMI, and older Age each raise the of diabetes.
Each of the three coefficients are significant and well below p < 0.05 level. Age has the largest p-value at p = 0.000238. The other p-values are much smaller. As stated, all three coefficients are significant.
Q7 — Confusion Matrix
Use threshold 0.5 to convert predicted probabilities into 0/1 predictions, then build a confusion matrix.
# Hint:
# data$pred_prob <- predict(model_q6, data, type = "response")
# data$pred_class <- ifelse(data$pred_prob > 0.5, 1, 0)
# table(Actual = data$Outcome, Predicted = data$pred_class)
data$pred_prob <- predict(model_q6, data, type = "response")
data$pred_class <- ifelse(data$pred_prob > 0.5, 1, 0)
table(Actual = data$Outcome, Predicted = data$pred_class)## Predicted
## Actual 0 1
## 0 429 59
## 1 114 150Report the confusion matrix counts: TP, TN, FP, FN.
The true positive, TP, count is 150. The true negative, TN, count is 429. The false positive, FP, count is 59. The false negative , FN, count is 114.
Q8 — Accuracy, Precision, Recall
From your confusion matrix, calculate accuracy, precision, and recall.
# Your code:
confusion_matrix <- table(Actual = data$Outcome, Predicted = data$pred_class)
TP <- confusion_matrix["1", "1"] #150
TN <- confusion_matrix["0", "0"] #429
FP <- confusion_matrix["0", "1"] #59
FN <- confusion_matrix["1", "0"] #114
TP## [1] 150TN## [1] 429FP## [1] 59FN## [1] 114accuracy <- (TP + TN) / sum(confusion_matrix)
precision <- TP / (TP + FP)
recall <- TP / (TP + FN)
accuracy## [1] 0.7699468precision## [1] 0.7177033recall## [1] 0.5681818Report all three values. In a medical screening context, which is more important — precision or recall? Why?
Accuracy = 0.770 (77.0%), Precision = 0.718 (71.8%), Recall = 0.568 (56.8%)
Recall measures how many of the actual positive cases we correctly found. In a medical screening context, recall matters more because we want to get as many of the actual positive cases as possible. Missing a real case is the more harmful error.
Precision gives a measure of how many of the people who test positive, how many actually have the sickness. Medically, if someone tests positive and further testing shows the person does not have the sickness, then it is not fun for the person, but it is safer to do further testing on someone who actually does not have a positive case.
Q9 — ROC and AUC
Plot the ROC curve and compute the AUC.
# install.packages("pROC") if needed
library(pROC)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, var# Your code:
roc_obj <- roc(data$Outcome, data$pred_prob)## Setting levels: control = 0, case = 1## Setting direction: controls < casesplot(roc_obj, main = "ROC Curve — Pima Diabetes Model", col = "blue")
auc(roc_obj)## Area under the curve: 0.828Report the AUC. Is your model closer to random guessing (AUC = 0.5) or perfect (AUC = 1)? Describe its overall performance in one sentence.
“The AUC is 0.828, which is much closer to 1.0 (perfect) than 0.5 (random guessing), indicating the model does a good job of distinguishing patients with diabetes from those without.”(Using the template)
The recall at the default 0.5 threshold could be improved. (data$pred_prob > 0.5 threshold, 1, 0)