In this homework, you will focus on:
- estimating simple regression lines and interpreting results,
- understanding coefficient interpretation and prediction,
- conducting statistical inference,
- evaluating model fit, and
- checking assumptions using residual diagnostics.
You will work with the built-in U.S. state dataset
state.x77, which contains socioeconomic
and demographic indicators for 50 states.
Use the help command to learn more about the dataset.
?state
Focus on the state.x77 matrix. Convert it into a data frame and explore its structure.
# Convert matrix to data frame
df_state <- as.data.frame(state.x77)
# Display structure and first rows
str(df_state)
## 'data.frame': 50 obs. of 8 variables:
## $ Population: num 3615 365 2212 2110 21198 ...
## $ Income : num 3624 6315 4530 3378 5114 ...
## $ Illiteracy: num 2.1 1.5 1.8 1.9 1.1 0.7 1.1 0.9 1.3 2 ...
## $ Life Exp : num 69 69.3 70.5 70.7 71.7 ...
## $ Murder : num 15.1 11.3 7.8 10.1 10.3 6.8 3.1 6.2 10.7 13.9 ...
## $ HS Grad : num 41.3 66.7 58.1 39.9 62.6 63.9 56 54.6 52.6 40.6 ...
## $ Frost : num 20 152 15 65 20 166 139 103 11 60 ...
## $ Area : num 50708 566432 113417 51945 156361 ...
head(df_state)
## Population Income Illiteracy Life Exp Murder HS Grad Frost Area
## Alabama 3615 3624 2.1 69.05 15.1 41.3 20 50708
## Alaska 365 6315 1.5 69.31 11.3 66.7 152 566432
## Arizona 2212 4530 1.8 70.55 7.8 58.1 15 113417
## Arkansas 2110 3378 1.9 70.66 10.1 39.9 65 51945
## California 21198 5114 1.1 71.71 10.3 62.6 20 156361
## Colorado 2541 4884 0.7 72.06 6.8 63.9 166 103766
For each independent variable \(X_j\), write a research question stating your expected sign of \(\beta_j\). Formulate hypotheses to test whether Income is related to each of the other variables in the dataset.
(Write your hypotheses and explanations below this section.)
Population: not clear
Illiteracy: higher illiteracy indicate lower human capital -> lower average income.
Life Exp: with years people are become more developed and experienced -> can lead to higher income.
Murder: higher murder rates could indicate bad quality of life -> lower income
HS Grad: if there are better education system and more qualificated people in the country -> higher income
Frost: not clear
Ares: not clear
Create scatter plots of Income against each independent variable.
Describe what you observe from the graphs — is the relationship positive, negative, or unclear?
(Write your hypotheses and explanations below this section.)
par(mfrow = c(2, 2))
#Income vs Population
plot(df_state$Population, df_state$Income,
main = "Income vs Population",
xlab = "Population", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ Population, data = df_state), col = 6, lwd = 2)
#Income vs Illiteracy
plot(df_state$Illiteracy, df_state$Income,
main = "Income vs Illiteracy",
xlab = "Illiteracy (%)", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ Illiteracy, data = df_state), col = 6, lwd = 2)
#Income vs Life Exp
plot(df_state[["Life Exp"]], df_state$Income,
main = "Income vs Life Exp",
xlab = "Life Exp (years)", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ `Life Exp`, data = df_state), col = 6, lwd = 2)
#Income vs Murder
plot(df_state$Murder, df_state$Income,
main = "Income vs Murder",
xlab = "Murder (rate per 100k)", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ Murder, data = df_state), col = 6, lwd = 2)
par(mfrow = c(2, 2))
#Income vs HS Grad
plot(df_state[["HS Grad"]], df_state$Income,
main = "Income vs HS Grad",
xlab = "HS Grad (%)", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ `HS Grad`, data = df_state), col = 6, lwd = 2)
#Income vs Frost
plot(df_state$Frost, df_state$Income,
main = "Income vs Frost",
xlab = "Frost (mean of days)", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ Frost, data = df_state), col = 6, lwd = 2)
#Income vs Area
plot(df_state$Area, df_state$Income,
main = "Income vs Area",
xlab = "Area (square miles)", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ Area, data = df_state), col = 6, lwd = 2)
Income vs Population - there is a weak correlation
Income vs Illiteracy - there is a negative correlation -> higher illiteracy rates lead to lower incomes.
Income vs Life Exp - the correlation is positive -> states with longer life expectancy have higher incomes.
Income vs Murder - negative correlation -> in regions with higher crime rates, the average income is lower.
Income vs HS Grad - positive correlation -> higher education levels lead to higher incomes.
Income vs Frost - weak correlation
Income vs Area - weak correlation
Compute the correlation matrix for all variables and interpret your results.
cor_matrix <- round(cor(df_state), 3)
cor_matrix
## Population Income Illiteracy Life Exp Murder HS Grad Frost Area
## Population 1.000 0.208 0.108 -0.068 0.344 -0.098 -0.332 0.023
## Income 0.208 1.000 -0.437 0.340 -0.230 0.620 0.226 0.363
## Illiteracy 0.108 -0.437 1.000 -0.588 0.703 -0.657 -0.672 0.077
## Life Exp -0.068 0.340 -0.588 1.000 -0.781 0.582 0.262 -0.107
## Murder 0.344 -0.230 0.703 -0.781 1.000 -0.488 -0.539 0.228
## HS Grad -0.098 0.620 -0.657 0.582 -0.488 1.000 0.367 0.334
## Frost -0.332 0.226 -0.672 0.262 -0.539 0.367 1.000 0.059
## Area 0.023 0.363 0.077 -0.107 0.228 0.334 0.059 1.000
From the correlation matrix we see that income is most closely related to HS Grad(0.62) and Life Exp(0.34) -> education and life expectancy increase along with income.
Also there is a negative correlation with Illiteracy (−0.44) and Murder(−0.23) -> higher levels of illiteracy and crime reduce income.
Population, Frost and Area -> weak positive correlation.
Determine how many observations (states) are included in the dataset.
nrow(df_state)
## [1] 50
50 observations in the dataset
Produce summary statistics for all variables.
summary(df_state)
## Population Income Illiteracy Life Exp
## Min. : 365 Min. :3098 Min. :0.500 Min. :67.96
## 1st Qu.: 1080 1st Qu.:3993 1st Qu.:0.625 1st Qu.:70.12
## Median : 2838 Median :4519 Median :0.950 Median :70.67
## Mean : 4246 Mean :4436 Mean :1.170 Mean :70.88
## 3rd Qu.: 4968 3rd Qu.:4814 3rd Qu.:1.575 3rd Qu.:71.89
## Max. :21198 Max. :6315 Max. :2.800 Max. :73.60
## Murder HS Grad Frost Area
## Min. : 1.400 Min. :37.80 Min. : 0.00 Min. : 1049
## 1st Qu.: 4.350 1st Qu.:48.05 1st Qu.: 66.25 1st Qu.: 36985
## Median : 6.850 Median :53.25 Median :114.50 Median : 54277
## Mean : 7.378 Mean :53.11 Mean :104.46 Mean : 70736
## 3rd Qu.:10.675 3rd Qu.:59.15 3rd Qu.:139.75 3rd Qu.: 81162
## Max. :15.100 Max. :67.30 Max. :188.00 Max. :566432
Produce histograms for the variables.
Write a brief comment on the distributions and variability of the data.
Do you need to use logarithmic form?
plot_hist_grid <- function(data, vars) {
par(mfrow = c(2, 2))
for (v in vars) {
x <- as.numeric(data[[v]])
hist(x,
main = v,
xlab = v,
col = "lightgray",
border = "white",
breaks = 8)
}
par(mfrow = c(1, 1))
}
plot_hist_grid(df_state, c("Population", "Income", "Illiteracy", "Life Exp"))
plot_hist_grid(df_state, c("Murder", "HS Grad", "Frost", "Area"))
Population and Area have long right tails -> need to use a
logarithmic form.
Income, Life Exp, HS Grad and Frost are closer to a normal distribution.
Illiteracy and Murder are quite skewed.
Estimate simple regressions of Income on each of the other variables in the dataset.
For each model:
-Present the regression summary.
-Interpret the coefficients and check their significance.
-Comment on ( \(R^2\) ).
predictors <- c("Population", "Illiteracy", "`Life Exp`",
"Murder", "`HS Grad`", "Frost", "Area")
linear_regression <- lapply(predictors, function(var) {
formula <- as.formula(paste("Income ~", var))
lm(formula, data = df_state)
})
lapply(linear_regression, summary)
## [[1]]
##
## Call:
## lm(formula = formula, data = df_state)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1283.19 -417.70 46.62 278.90 1990.44
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.314e+03 1.191e+02 36.226 <2e-16 ***
## Population 2.866e-02 1.943e-02 1.475 0.147
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 607.2 on 48 degrees of freedom
## Multiple R-squared: 0.04336, Adjusted R-squared: 0.02343
## F-statistic: 2.176 on 1 and 48 DF, p-value: 0.1467
##
##
## [[2]]
##
## Call:
## lm(formula = formula, data = df_state)
##
## Residuals:
## Min 1Q Median 3Q Max
## -948.89 -376.20 -49.77 347.00 2024.60
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4951.3 172.3 28.739 < 2e-16 ***
## Illiteracy -440.6 130.9 -3.367 0.00151 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 558.4 on 48 degrees of freedom
## Multiple R-squared: 0.191, Adjusted R-squared: 0.1742
## F-statistic: 11.34 on 1 and 48 DF, p-value: 0.001505
##
##
## [[3]]
##
## Call:
## lm(formula = formula, data = df_state)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1023.75 -467.39 -34.34 339.73 2123.51
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6603.48 4404.26 -1.499 0.1403
## `Life Exp` 155.75 62.13 2.507 0.0156 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 583.8 on 48 degrees of freedom
## Multiple R-squared: 0.1158, Adjusted R-squared: 0.09735
## F-statistic: 6.285 on 1 and 48 DF, p-value: 0.01562
##
##
## [[4]]
##
## Call:
## lm(formula = formula, data = df_state)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1141.64 -483.71 -19.01 403.31 2029.40
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4718.36 192.51 24.510 <2e-16 ***
## Murder -38.30 23.38 -1.638 0.108
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 604.2 on 48 degrees of freedom
## Multiple R-squared: 0.05294, Adjusted R-squared: 0.03321
## F-statistic: 2.683 on 1 and 48 DF, p-value: 0.108
##
##
## [[5]]
##
## Call:
## lm(formula = formula, data = df_state)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1083.13 -277.41 -34.15 241.46 1238.17
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1931.105 462.739 4.173 0.000125 ***
## `HS Grad` 47.162 8.616 5.474 1.58e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 487.1 on 48 degrees of freedom
## Multiple R-squared: 0.3843, Adjusted R-squared: 0.3715
## F-statistic: 29.96 on 1 and 48 DF, p-value: 1.579e-06
##
##
## [[6]]
##
## Call:
## lm(formula = formula, data = df_state)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1192.12 -483.89 -22.45 383.72 1752.04
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4156.380 193.531 21.477 <2e-16 ***
## Frost 2.675 1.662 1.609 0.114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 604.7 on 48 degrees of freedom
## Multiple R-squared: 0.0512, Adjusted R-squared: 0.03144
## F-statistic: 2.59 on 1 and 48 DF, p-value: 0.1141
##
##
## [[7]]
##
## Call:
## lm(formula = formula, data = df_state)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1276.5 -446.6 77.7 435.0 1084.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.251e+03 1.067e+02 39.841 < 2e-16 ***
## Area 2.616e-03 9.684e-04 2.702 0.00951 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 578.4 on 48 degrees of freedom
## Multiple R-squared: 0.132, Adjusted R-squared: 0.1139
## F-statistic: 7.299 on 1 and 48 DF, p-value: 0.009505
Population: p = 0.15 -> NOT statistically significant. Illiteracy: p = 0.001 -> statistically significant Life Exp: p = 0.016 -> statistically significant. Murder: p = 0.11 -> NOT statistically significant. HS Grad: p < 1.58e-06 -> statistically significant Frost: p = 0.11 -> NOT statistically significant. Area: p = 0.009 -> statistically significant
Looking at R^2, we can say that HS Grad, in other words, education, significantly increases income (R^2 = 0.38).
Check the confidence intervals for each of your regression coefficients. Use the 95% and 99% confidence levels.
ci_model <- function(vars) {
for (v in vars) {
cat("\n", v ,"\n")
fit <- lm(df_state[["Income"]] ~ df_state[[v]])
print(confint(fit, level = 0.95))
print(confint(fit, level = 0.99))
}
}
predictors <- c("Population", "Illiteracy", "Life Exp",
"Murder", "HS Grad", "Frost", "Area")
ci_model(predictors)
##
## Population
## 2.5 % 97.5 %
## (Intercept) 4074.65385861 4553.5465782
## df_state[[v]] -0.01040813 0.0677269
## 0.5 % 99.5 %
## (Intercept) 3994.67671892 4.633524e+03
## df_state[[v]] -0.02345701 8.077578e-02
##
## Illiteracy
## 2.5 % 97.5 %
## (Intercept) 4604.9180 5297.7216
## df_state[[v]] -703.7516 -177.4789
## 0.5 % 99.5 %
## (Intercept) 4489.2168 5413.42277
## df_state[[v]] -791.6414 -89.58906
##
## Life Exp
## 2.5 % 97.5 %
## (Intercept) -15458.84891 2251.8831
## df_state[[v]] 30.83403 280.6643
## 0.5 % 99.5 %
## (Intercept) -18416.6172 5209.651
## df_state[[v]] -10.8887 322.387
##
## Murder
## 2.5 % 97.5 %
## (Intercept) 4331.2977 5105.416384
## df_state[[v]] -85.3078 8.713319
## 0.5 % 99.5 %
## (Intercept) 4202.0165 5234.69753
## df_state[[v]] -101.0097 24.41525
##
## HS Grad
## 2.5 % 97.5 %
## (Intercept) 1000.70545 2861.50393
## df_state[[v]] 29.83854 64.48606
## 0.5 % 99.5 %
## (Intercept) 689.94411 3172.26526
## df_state[[v]] 24.05225 70.27234
##
## Frost
## 2.5 % 97.5 %
## (Intercept) 3767.2595304 4545.500257
## df_state[[v]] -0.6667141 6.016515
## 0.5 % 99.5 %
## (Intercept) 3637.289989 4675.469799
## df_state[[v]] -1.782842 7.132643
##
## Area
## 2.5 % 97.5 %
## (Intercept) 4.036208e+03 4.465252e+03
## df_state[[v]] 6.692661e-04 4.563442e-03
## 0.5 % 99.5 %
## (Intercept) 3.964556e+03 4.536904e+03
## df_state[[v]] 1.892204e-05 5.213786e-03
Make scatter plots with regression lines fitted using OLS.
par(mfrow = c(2, 2))
#Income vs Population
plot(df_state$Population, df_state$Income,
main = "Income vs Population",
xlab = "Population", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ Population, data = df_state), col = 6, lwd = 2)
#Income vs Illiteracy
plot(df_state$Illiteracy, df_state$Income,
main = "Income vs Illiteracy",
xlab = "Illiteracy (%)", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ Illiteracy, data = df_state), col = 6, lwd = 2)
#Income vs Life Exp
plot(df_state[["Life Exp"]], df_state$Income,
main = "Income vs Life Exp",
xlab = "Life Exp (years)", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ `Life Exp`, data = df_state), col = 6, lwd = 2)
#Income vs Murder
plot(df_state$Murder, df_state$Income,
main = "Income vs Murder",
xlab = "Murder (rate per 100k)", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ Murder, data = df_state), col = 6, lwd = 2)
par(mfrow = c(2, 2))
#Income vs HS Grad
plot(df_state[["HS Grad"]], df_state$Income,
main = "Income vs HS Grad",
xlab = "HS Grad (%)", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ `HS Grad`, data = df_state), col = 6, lwd = 2)
#Income vs Frost
plot(df_state$Frost, df_state$Income,
main = "Income vs Frost",
xlab = "Frost (mean of days)", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ Frost, data = df_state), col = 6, lwd = 2)
#Income vs Area
plot(df_state$Area, df_state$Income,
main = "Income vs Area",
xlab = "Area (square miles)", ylab = "Income",
pch = 19, col = "blue")
abline(lm(Income ~ Area, data = df_state), col = 6, lwd = 2)
Conduct residual analysis for each model:
Check for outliers.
Produce diagnostic plots.
models <- list(
lm(Income ~ Population, df_state),
lm(Income ~ Illiteracy, df_state),
lm(Income ~ `Life Exp`, df_state),
lm(Income ~ Murder, df_state),
lm(Income ~ `HS Grad`, df_state),
lm(Income ~ Frost, df_state),
lm(Income ~ Area, df_state)
)
par(mfrow = c(2, 2))
for (i in 1:4) {
plot(fitted(models[[i]]), resid(models[[i]]),
xlab = "Fitted values", ylab = "Residuals",
main = paste("Residuals vs Fitted:", names(df_state)[i]),
pch = 19, col = "gray")
abline(h = 0, col = "red", lwd = 2)
lines(lowess(fitted(models[[i]]), resid(models[[i]])), col = "blue", lwd = 2)
}
par(mfrow = c(2, 2))
for (i in 5:7) {
plot(fitted(models[[i]]), resid(models[[i]]),
xlab = "Fitted values", ylab = "Residuals",
main = paste("Residuals vs Fitted:", names(df_state)[i]),
pch = 19, col = "gray")
abline(h = 0, col = "red", lwd = 2)
lines(lowess(fitted(models[[i]]), resid(models[[i]])), col = "blue", lwd = 2)
}
Test the residuals for normality and interpret the results.
shapiro.test(resid(lm(Income ~ Population, data = df_state)))
##
## Shapiro-Wilk normality test
##
## data: resid(lm(Income ~ Population, data = df_state))
## W = 0.96866, p-value = 0.2041
shapiro.test(resid(lm(Income ~ Illiteracy, data = df_state)))
##
## Shapiro-Wilk normality test
##
## data: resid(lm(Income ~ Illiteracy, data = df_state))
## W = 0.94721, p-value = 0.02624
shapiro.test(resid(lm(Income ~ `Life Exp`, data = df_state)))
##
## Shapiro-Wilk normality test
##
## data: resid(lm(Income ~ `Life Exp`, data = df_state))
## W = 0.94793, p-value = 0.02806
shapiro.test(resid(lm(Income ~ Murder, data = df_state)))
##
## Shapiro-Wilk normality test
##
## data: resid(lm(Income ~ Murder, data = df_state))
## W = 0.96523, p-value = 0.1471
shapiro.test(resid(lm(Income ~ `HS Grad`, data = df_state)))
##
## Shapiro-Wilk normality test
##
## data: resid(lm(Income ~ `HS Grad`, data = df_state))
## W = 0.98322, p-value = 0.6934
shapiro.test(resid(lm(Income ~ Frost, data = df_state)))
##
## Shapiro-Wilk normality test
##
## data: resid(lm(Income ~ Frost, data = df_state))
## W = 0.9811, p-value = 0.5992
shapiro.test(resid(lm(Income ~ Area, data = df_state)))
##
## Shapiro-Wilk normality test
##
## data: resid(lm(Income ~ Area, data = df_state))
## W = 0.97625, p-value = 0.4067
For Illiteracy and Life Exp p-value < 0.05, the H0 is rejected, suggesting deviations from a normal distribution.
Write a short analytical comment (4–5 sentences) summarizing your findings.
Discuss:
Which variables appear most strongly related to income.
Whether the observed relationships align with your expectations.
What possible economic mechanisms could explain these relationships.
(Write your hypotheses and explanations below this section.)
The strongest correlation with income level is observed for the HS Grad variable, which is logical, as higher levels of education among the population usually lead to higher incomes. Life Exp and Area also have a notable positive impact, while Illiteracy has a negative correlation with income, which is also quite understandable.
Economically, this can be explained by the fact that regions with good education, longer life expectancy and more developed economies have higher incomes.
Submit your HW in rendered .html format (Press “Knit”)