HW9
2026-08-21
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.3 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsAllCountries <- read_csv("AllCountries.csv")## Rows: 217 Columns: 26
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Country, Code
## dbl (24): LandArea, Population, Density, GDP, Rural, CO2, PumpPrice, Militar...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.str(AllCountries)## spc_tbl_ [217 × 26] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ Country : chr [1:217] "Afghanistan" "Albania" "Algeria" "American Samoa" ...
## $ Code : chr [1:217] "AFG" "ALB" "DZA" "ASM" ...
## $ LandArea : num [1:217] 652.86 27.4 2381.74 0.2 0.47 ...
## $ Population : num [1:217] 37.172 2.866 42.228 0.055 0.077 ...
## $ Density : num [1:217] 56.9 104.6 17.7 277.3 163.8 ...
## $ GDP : num [1:217] 521 5254 4279 NA 42030 ...
## $ Rural : num [1:217] 74.5 39.7 27.4 12.8 11.9 34.5 75.4 8.1 36.9 56.6 ...
## $ CO2 : num [1:217] 0.29 1.98 3.74 NA 5.83 1.29 5.74 4.78 1.9 8.41 ...
## $ PumpPrice : num [1:217] 0.7 1.36 0.28 NA NA 0.97 NA 1.1 0.77 NA ...
## $ Military : num [1:217] 3.72 4.08 13.81 NA NA ...
## $ Health : num [1:217] 2.01 9.51 10.73 NA 14.02 ...
## $ ArmedForces : num [1:217] 323 9 317 NA NA 117 0 105 49 NA ...
## $ Internet : num [1:217] 11.4 71.8 47.7 NA 98.9 14.3 76 75.8 69.7 97.2 ...
## $ Cell : num [1:217] 67.4 123.7 111 NA 104.4 ...
## $ HIV : num [1:217] NA 0.1 0.1 NA NA 1.9 NA 0.4 0.2 NA ...
## $ Hunger : num [1:217] 30.3 5.5 4.7 NA NA 23.9 NA 3.8 4.3 NA ...
## $ Diabetes : num [1:217] 9.6 10.1 6.7 NA 8 3.9 13.2 5.5 7.1 11.6 ...
## $ BirthRate : num [1:217] 32.5 11.7 22.3 NA NA 41.3 16.1 17 13.1 11 ...
## $ DeathRate : num [1:217] 6.6 7.5 4.8 NA NA 8.4 5.8 7.6 9.7 8.9 ...
## $ ElderlyPop : num [1:217] 2.6 13.6 6.4 NA NA 2.5 7.2 11.3 11.4 13.6 ...
## $ LifeExpectancy: num [1:217] 64 78.5 76.3 NA NA 61.8 76.5 76.7 74.8 76 ...
## $ FemaleLabor : num [1:217] 50.3 55.9 16.4 NA NA 76.4 NA 57.1 55.8 NA ...
## $ Unemployment : num [1:217] 1.5 13.9 12.1 NA NA 7.3 NA 9.5 17.7 NA ...
## $ Energy : num [1:217] NA 808 1328 NA NA ...
## $ Electricity : num [1:217] NA 2309 1363 NA NA ...
## $ Developed : num [1:217] NA 1 1 NA NA 1 NA 2 1 NA ...
## - attr(*, "spec")=
## .. cols(
## .. Country = col_character(),
## .. Code = col_character(),
## .. LandArea = col_double(),
## .. Population = col_double(),
## .. Density = col_double(),
## .. GDP = col_double(),
## .. Rural = col_double(),
## .. CO2 = col_double(),
## .. PumpPrice = col_double(),
## .. Military = col_double(),
## .. Health = col_double(),
## .. ArmedForces = col_double(),
## .. Internet = col_double(),
## .. Cell = col_double(),
## .. HIV = col_double(),
## .. Hunger = col_double(),
## .. Diabetes = col_double(),
## .. BirthRate = col_double(),
## .. DeathRate = col_double(),
## .. ElderlyPop = col_double(),
## .. LifeExpectancy = col_double(),
## .. FemaleLabor = col_double(),
## .. Unemployment = col_double(),
## .. Energy = col_double(),
## .. Electricity = col_double(),
## .. Developed = col_double()
## .. )
## - attr(*, "problems")=<pointer: 0x5edfbd5f47e0>summary(AllCountries)## Country Code LandArea Population
## Length :217 Length :217 Min. : 0.01 Min. : 0.0120
## N.unique :217 N.unique :217 1st Qu.: 10.83 1st Qu.: 0.7728
## N.blank : 0 N.blank : 0 Median : 94.28 Median : 6.5725
## Min.nchar: 4 Min.nchar: 3 Mean : 608.38 Mean : 35.0335
## Max.nchar: 30 Max.nchar: 3 3rd Qu.: 446.30 3rd Qu.: 25.0113
## Max. :16376.87 Max. :1392.7300
## NAs :8 NAs :1
## Density GDP Rural CO2
## Min. : 0.1 Min. : 275 Min. : 0.00 Min. : 0.0400
## 1st Qu.: 37.5 1st Qu.: 2032 1st Qu.:19.62 1st Qu.: 0.8575
## Median : 92.1 Median : 5950 Median :38.15 Median : 2.7550
## Mean : 361.4 Mean : 14733 Mean :39.10 Mean : 4.9780
## 3rd Qu.: 219.8 3rd Qu.: 17298 3rd Qu.:57.83 3rd Qu.: 6.2525
## Max. :20777.5 Max. :114340 Max. :87.00 Max. :43.8600
## NAs :8 NAs :30 NAs :3 NAs :13
## PumpPrice Military Health ArmedForces
## Min. :0.1100 Min. : 0.000 Min. : 0.000 Min. : 0.0
## 1st Qu.:0.7450 1st Qu.: 3.015 1st Qu.: 6.157 1st Qu.: 12.0
## Median :0.9800 Median : 4.650 Median : 9.605 Median : 31.5
## Mean :0.9851 Mean : 6.178 Mean :10.597 Mean : 162.1
## 3rd Qu.:1.1800 3rd Qu.: 8.445 3rd Qu.:13.713 3rd Qu.: 146.5
## Max. :2.0000 Max. :31.900 Max. :39.460 Max. :3031.0
## NAs :50 NAs :67 NAs :29 NAs :49
## Internet Cell HIV Hunger
## Min. : 1.30 Min. : 13.70 Min. : 0.100 Min. : 2.50
## 1st Qu.:29.18 1st Qu.: 83.83 1st Qu.: 0.175 1st Qu.: 2.50
## Median :58.35 Median :110.00 Median : 0.400 Median : 6.50
## Mean :54.47 Mean :107.05 Mean : 1.941 Mean :11.25
## 3rd Qu.:78.92 3rd Qu.:127.50 3rd Qu.: 1.400 3rd Qu.:14.80
## Max. :98.90 Max. :328.80 Max. :27.400 Max. :61.80
## NAs :13 NAs :15 NAs :81 NAs :52
## Diabetes BirthRate DeathRate ElderlyPop
## Min. : 1.000 Min. : 7.00 Min. : 1.600 Min. : 1.200
## 1st Qu.: 5.350 1st Qu.:11.40 1st Qu.: 5.800 1st Qu.: 3.600
## Median : 7.200 Median :17.85 Median : 7.250 Median : 6.600
## Mean : 8.542 Mean :20.11 Mean : 7.683 Mean : 8.953
## 3rd Qu.:10.750 3rd Qu.:27.65 3rd Qu.: 9.350 3rd Qu.:14.500
## Max. :30.500 Max. :47.80 Max. :15.500 Max. :27.500
## NAs :10 NAs :15 NAs :15 NAs :24
## LifeExpectancy FemaleLabor Unemployment Energy
## Min. :52.20 Min. : 6.20 Min. : 0.100 Min. : 66
## 1st Qu.:66.90 1st Qu.:50.15 1st Qu.: 3.400 1st Qu.: 738
## Median :74.30 Median :60.60 Median : 5.600 Median : 1574
## Mean :72.46 Mean :57.95 Mean : 7.255 Mean : 2664
## 3rd Qu.:77.70 3rd Qu.:69.25 3rd Qu.: 9.400 3rd Qu.: 3060
## Max. :84.70 Max. :85.80 Max. :30.200 Max. :17923
## NAs :18 NAs :30 NAs :30 NAs :82
## Electricity Developed
## Min. : 39 Min. :1.00
## 1st Qu.: 904 1st Qu.:1.00
## Median : 2620 Median :2.00
## Mean : 4270 Mean :1.81
## 3rd Qu.: 5600 3rd Qu.:3.00
## Max. :53832 Max. :3.00
## NAs :76 NAs :75colSums(is.na(AllCountries))## Country Code LandArea Population Density
## 0 0 8 1 8
## GDP Rural CO2 PumpPrice Military
## 30 3 13 50 67
## Health ArmedForces Internet Cell HIV
## 29 49 13 15 81
## Hunger Diabetes BirthRate DeathRate ElderlyPop
## 52 10 15 15 24
## LifeExpectancy FemaleLabor Unemployment Energy Electricity
## 18 30 30 82 76
## Developed
## 75#Simple Linear Regression Using the AllCountries dataset, fit a simple linear regression model to predict LifeExpectancy (average life expectancy in years) based on GDP (gross domestic product per capita in $US). Report the intercept and slope coefficients and interpret their meaning in the context of the dataset. What does the R² value tell you about how well GDP explains variation in life expectancy across countries?
# Fit simple linear regression model
modelcountries <- lm(LifeExpectancy ~ GDP, data = AllCountries)
summary(modelcountries)##
## Call:
## lm(formula = LifeExpectancy ~ GDP, data = AllCountries)
##
## 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-16coef_countries <- coef(modelcountries)
cat("\nIntercept:", round(coef_countries[1], 3), "\n")##
## Intercept: 68.422cat("GDP Slope:", round(coef_countries[2], 6), "\n\n")## GDP Slope: 0.000248r_squared <- summary(modelcountries)$r.squared
cat("R-squared:", round(r_squared, 3), "\n")## R-squared: 0.43If a country had zero income per person (which isn’t realistic) this is what the model predicts their life expectancy would be. It’s mostly just a mathematical starting point for the equation, not something we can really interpret. Then, for every $1,000 more money people earn on average, we need to look at how much longer they live. The number tells you the yearly increase in life expectancy. For our R-squared value, we use that to help us figure out the impact of different factors. What percentage of these variations in lifespacn, can we blame on differences in income? A value like 0.30 means 30% — so income matters, but so do the other factors.
#Multiple Linear Regression Fit a multiple linear regression model to predict LifeExpectancy using GDP, Health (percentage of government expenditures on healthcare), and Internet (percentage of population with internet access) as predictors. Interpret the coefficient for Health, explaining what it means in terms of life expectancy while controlling for GDP and Internet. How does the adjusted R² compare to the simple regression model from Question 1, and what does this suggest about the additional predictors?
# Fit multiple linear regression model
multi_model <- lm(LifeExpectancy ~ GDP + Health + Internet, data = AllCountries)
summary(multi_model)##
## Call:
## lm(formula = LifeExpectancy ~ GDP + Health + Internet, data = AllCountries)
##
## 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-16adj_r_squared_multi <- summary(multi_model)$adj.r.squared
cat("\nAdjusted R-squared (Multiple Model):", round(adj_r_squared_multi, 3), "\n")##
## Adjusted R-squared (Multiple Model): 0.716cat("R-squared (Simple Model):", round(summary(modelcountries)$r.squared, 3), "\n")## R-squared (Simple Model): 0.43cat("Difference:", round(adj_r_squared_multi - r_squared, 3), "\n")## Difference: 0.286If two countries have the same income and internet access, but one spends 1% more of its government budget on healthcare, the life expectancy difference between them should be about the Health coefficient number of years. When we added healthcare spending and internet access, we need to see if our predictions improved. If the adjusted R² went up, which it did, from 0.43 to 0.716, these new variables helped explain more of the differences in how long people live.
#Checking Assumptions For the simple linear regression model from Question 1 (LifeExpectancy ~ GDP), describe how you would check the assumptions of homoscedasticity and normality of residuals. For each assumption, explain what an ideal outcome would look like and what a violation might indicate about the model’s reliability for predicting life expectancy. Afterwords, code your answer and reflect if it matched the ideal outcome.
#Homoscedasticity
plot(modelcountries, which = 1, main = "Residuals vs Fitted")
#Normality
plot(modelcountries, which = 2, main = "Normal Q-Q")
#Scale-Location Plot
plot(modelcountries, which = 3, main = "Scale-Location")
#Histogram of Residuals
hist(modelcountries$residuals, breaks = 30,
main = "Histogram of Residuals",
xlab = "Residuals", col = "blue")
par(mfrow = c(1, 1)) # Reset plotting panelI examined the Residuals vs Fitted plot from plot. Ideally, residuals should scatter randomly around zero with equal spread across all GDP values. In practice, I observed a fan shape where variance increased at higher GDP levels, indicating heteroscedasticity. This means the model predicts life expectancy better for some countries than others. I examined both the Q-Q plot and a histogram of residuals. Ideally, points should follow the diagonal line and the histogram should be bell-shaped. My data in the histogram showed deviations in the tails and slight left-skewness (median residual = 1.55, not 0), indicating non-normal residuals. The range from -16.35 to +9.33 was asymmetric. The QQ plot was mostly following the diagonal line, though. These findings mean we should examine the p-values and standard error and check the fit further to see if we can make this model more accurate. These are not very egregious violations, however. The model is mostly accurate.
#Diagnosing Model Fit (RMSE and Residuals) For the multiple regression model from Question 2 (LifeExpectancy ~ GDP + Health + Internet), calculate the RMSE and explain what it represents in the context of predicting life expectancy. How would large residuals for certain countries (e.g., those with unusually high or low life expectancy) affect your confidence in the model’s predictions, and what might you investigate further.
predictions_multi <- predict(multi_model, newdata = AllCountries)
residuals_multi <- AllCountries$LifeExpectancy - predictions_multi
RMSE <- sqrt(mean(residuals_multi^2, na.rm = TRUE))
cat("RMSE:", round(RMSE, 3), "years\n\n")## RMSE: 4.056 years# Create residual dataframe with country names
residual_df <- data.frame(
Country = AllCountries$Country,
Actual_LifeExpectancy = AllCountries$LifeExpectancy,
Predicted_LifeExpectancy = round(predictions_multi, 1),
Residual = round(residuals_multi, 1)
)
cat("--- Top 5 Highest Positive Residuals (Actual much higher than predicted) ---\n")## --- Top 5 Highest Positive Residuals (Actual much higher than predicted) ---print(head(residual_df[order(-residual_df$Residual), ], 5))## Country Actual_LifeExpectancy Predicted_LifeExpectancy Residual
## 16 Bangladesh 72.8 63.4 9.4
## 114 Libya 72.1 63.4 8.7
## 97 Italy 83.2 74.9 8.3
## 181 Sri Lanka 75.5 67.8 7.7
## 176 Solomon Islands 71.0 63.4 7.6cat("\n--- Top 5 Largest Negative Residuals (Actual much lower than predicted) ---\n")##
## --- Top 5 Largest Negative Residuals (Actual much lower than predicted) ---print(head(residual_df[order(residual_df$Residual), ], 5))## Country Actual_LifeExpectancy Predicted_LifeExpectancy Residual
## 48 Cote d'Ivoire 54.1 68.7 -14.6
## 112 Lesotho 54.6 67.5 -12.9
## 145 Nigeria 53.9 65.6 -11.7
## 171 Sierra Leone 52.2 63.6 -11.4
## 64 Eswatini 58.3 68.7 -10.4An RMSE of 4 years means if you pick any country, the prediction will usually be within about 4 years of the actual life expectancy. That’s pretty reasonable given life spans range from about 50 to 85 years worldwide. People live longer than the model expects in countries like Bangladesh and Libyam and people live less than the model expects in Lesotho and Nigeria. This indicates there are other factors at play in the countries — public health, diet, crime, war — that the model didn’t capture.
#Hypothetical Example (Multicollinearity in Multiple Regression) Suppose you are analyzing the AllCountries dataset and fit a multiple linear regression model to predict CO2 emissions (metric tons per capita) using Energy (kilotons of oil equivalent) and Electricity (kWh per capita) as predictors. You notice that Energy and Electricity are highly correlated. Explain how this multicollinearity might affect the interpretation of the regression coefficients and the reliability of the model.
cor_test <- cor(AllCountries$Energy, AllCountries$Electricity, use = "complete.obs")
cat("Correlation(Energy, Electricity):", round(cor_test, 3), "\n\n")## Correlation(Energy, Electricity): 0.797collinear_model <- lm(CO2 ~ Energy + Electricity, data = AllCountries)
summary(collinear_model)##
## Call:
## lm(formula = CO2 ~ Energy + Electricity, data = AllCountries)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.7559 -1.1406 -0.2020 0.7143 7.3751
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.998e-01 2.655e-01 3.012 0.00311 **
## Energy 3.122e-03 1.066e-04 29.290 < 2e-16 ***
## Electricity -7.044e-04 5.526e-05 -12.747 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.331 on 131 degrees of freedom
## (83 observations deleted due to missingness)
## Multiple R-squared: 0.899, Adjusted R-squared: 0.8974
## F-statistic: 582.8 on 2 and 131 DF, p-value: < 2.2e-16install.packages("car")## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.6'
## (as 'lib' is unspecified)library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## somevif_values <- vif(collinear_model)
cat("\nVariance Inflation Factors (VIF):\n")##
## Variance Inflation Factors (VIF):cat("\nVIF Rule of Thumb:\n")##
## VIF Rule of Thumb:cat("- VIF > 5: Moderate concern\n")## - VIF > 5: Moderate concerncat("- VIF > 10: Severe multicollinearity\n")## - VIF > 10: Severe multicollinearityThe VIF numbers are 2.74, which is below the threshold of 5. This means Energy and Electricity are correlated, but the correlation isn’t strong enough to seriously distort our model. However, Energy and Electricity are still moderately correlated (approximately r ≈ √[(2.74-1)/2.74] ≈ 0.79), which explains why the model needed to take their effect into account. This model is reliable, since both predictors are highly statistically significant (p < 0.001) and The overall R-squared of 0.899 is very explanatory. The low VIF values suggest this whole model is stable.