Homework 6 — Linear and Logistic Regression
Hernan Vasquez
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")
head(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:
model1 <- lm(LifeExpectancy ~ GDP, data = countries)
summary(model1)##
## Call:
## lm(formula = LifeExpectancy ~ GDP, data = countries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.352 -3.882 1.550 4.458 9.330
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.842e+01 5.415e-01 126.36 <2e-16 ***
## GDP 2.476e-04 2.141e-05 11.56 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.901 on 177 degrees of freedom
## (38 observations deleted due to missingness)
## Multiple R-squared: 0.4304, Adjusted R-squared: 0.4272
## F-statistic: 133.7 on 1 and 177 DF, p-value: < 2.2e-16Report 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 intercept: 68.42 Slope (GDP) : 0.0002476. The slope mean: For every one-unit increase in GDP, the predicted life expectancy increases by about 0.0002476 years in average
What does the R² value tell you about how well GDP explains life expectancy?
R² = 0.4304. The R² value is 0.4304. This means that aboiut 43.04% of the variation in life expectancy is explained by GDP.
Q2 — Multiple Linear Regression
Fit a multiple regression predicting LifeExpectancy from
GDP, Health, and Internet.
# Your code:
model2 <- lm(LifeExpectancy ~ GDP + Health + Internet, data = countries)
summary(model2)##
## Call:
## lm(formula = LifeExpectancy ~ GDP + Health + Internet, data = countries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.5662 -1.8227 0.4108 2.5422 9.4161
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.908e+01 8.149e-01 72.499 < 2e-16 ***
## GDP 2.367e-05 2.287e-05 1.035 0.302025
## Health 2.479e-01 6.619e-02 3.745 0.000247 ***
## Internet 1.903e-01 1.656e-02 11.490 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.104 on 169 degrees of freedom
## (44 observations deleted due to missingness)
## Multiple R-squared: 0.7213, Adjusted R-squared: 0.7164
## F-statistic: 145.8 on 3 and 169 DF, p-value: < 2.2e-16Interpret the coefficient on Health (controlling
for GDP and Internet).
Coefficient on Health = 0.2479 Holding GDP and Internet constant, a one-unit increase in Health is associated with an increase of about 0.248 years in life expectancy, on average.
How does the adjusted R² compare to the simple model in Q1? What does that suggest about adding predictors?
En Q1 -Adjusted R² = 0.4272
En Q2 -Adjusted R² = 0.7164
The adjusted R² increased from 0.4272 to 0.7164. This suggest that adding Health and Internet greatly improves the model and explains much more of the variation in 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?
Then code your check (residual plot + Q-Q plot of residuals) and reflect on what you see.
# Your code:
plot(model1)
Your reflection: The residual plot show a fairly random scatter around zero, suggesting constant variance (homoscedasticity). The Q-Q plot should show points close to the reference line, indicating that the residuals are approximately normally distributed. Overall, the assumptions appear to be reasonably satisfied. —
Q4 — Diagnosing Fit (RMSE)
For the multiple regression in Q2, calculate the RMSE (root mean squared error).
# Hint: sqrt(mean(residuals(model)^2))
rmse <- sqrt(mean(residuals(model2)^2))
rmse## [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 = 4.056417.
RMSE represents the typical prediction error of the model in years of life expectancy. A smaller RMSE indicates that the model’s predictions are more accurate.
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.
The interpretation of the coefficients: Multicollinearity makes the coefficients unstable and increase their standards errors, making it difficult to interpret the individual effect of each predictor.
The reliability of the model: It reduces the reliability of the model because the coefficient estimates can change noticeably with small changes in the data, even if the overall predictions remain similar.
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")
model3 <- glm(Outcome ~ Glucose + BMI + Age,
data = data,
family = "binomial")
summary(model3)##
## 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)?
Glucose: An increase in glucose raises the odds of diabetes
BMI: An increase in BMI raises the odds of diabetes
Age: An increase in age raises the odds of diabetes
Significant predictors: Glucose, BMI, and Age are all statistically significant because their p-values are less than 0.05.
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, data, type = "response")
# data$pred_class <- ifelse(data$pred_prob > 0.5, 1, 0)
# table(Actual = data$Outcome, Predicted = data$pred_class)
pred_prob <- predict(model3, type = "response")
pred_class <- ifelse(pred_prob > 0.5, 1, 0)
actula <- model3$y
table(Predicted = pred_class,
Actual = model.frame(model3)$Outcome)## Actual
## Predicted 0 1
## 0 429 114
## 1 59 150Report the confusion matrix counts: TP, TN, FP, FN.
Predicted 0 429 114 Predicted 1 59 150
TP = 150 TN = 429 FP = 59 FN = 114
Q8 — Accuracy, Precision, Recall
From your confusion matrix, calculate accuracy, precision, and recall.
# Your code:
TP <- 150
TN <- 429
FP <- 59
FN <- 114
accuracy <- (TP + TN) / (TP + TN + FP + FN)
precision <- TP / (TP + FP)
recall <- TP / (TP + FN)
accuracy## [1] 0.7699468 precision## [1] 0.7177033 recall## [1] 0.5681818Report all three values. In a medical screening context, which is more important — precision or recall? Why?
Accuracy: 0.770 Precision:0.718 Recall: 0.568 In a medical screening context, recall is more important than precision because missing a person who has the disease (a false negative) can have serious consequences.
Q9 — ROC and AUC
Plot the ROC curve and compute the AUC.
# install.packages("pROC") if needed
library(pROC)## Warning: package 'pROC' was built under R version 4.6.1## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, var# Your code:
library(pROC)
roc_obj <- roc(model.frame(model3)$Outcome,pred_prob)## Setting levels: control = 0, case = 1## Setting direction: controls < casesplot(roc_obj)
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.
AUC = 0.828
The model has good predictive performance and is much closer to perfect prediction than random guessing.