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.


1. Explore the dataset

Use the help command to learn more about the dataset.

?state

2. Variables and structure

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

3. Research questions

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

4. Visualize the data

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

5. Correlation analysis

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.

6. Number of observations

Determine how many observations (states) are included in the dataset.

nrow(df_state)
## [1] 50

50 observations in the dataset

7. Summary statistics

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

8. Distributions of the data

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.

9. Simple linear regressions

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).

10. Confidence intervals

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

11. Visualization with fitted regression lines

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)

12. Residual analysis

Conduct residual analysis for each model:

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)
}

13. Normality of residuals

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.

14. Final comment

Write a short analytical comment (4–5 sentences) summarizing your findings.

Discuss:

(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.

Deliverable

Submit your HW in rendered .html format (Press “Knit”)