Before You Start

This week you built regression models, interpreted their output, explored multi-variable graphs, and read about what happens when models go wrong. This challenge asks you to put those things together in one investigation.

Run this code to install the openintro package. It only needs to run once, so I set eval=FALSE here to keep the knit from reinstalling it every time.

install.packages("openintro")
library(tidyverse)
library(openintro)
data(county)
head(county)
## # A tibble: 6 × 15
##   name           state  pop2000 pop2010 pop2017 pop_change poverty homeownership
##   <chr>          <fct>    <dbl>   <dbl>   <int>      <dbl>   <dbl>         <dbl>
## 1 Autauga County Alaba…   43671   54571   55504       1.48    13.7          77.5
## 2 Baldwin County Alaba…  140415  182265  212628       9.19    11.8          76.7
## 3 Barbour County Alaba…   29038   27457   25270      -6.22    27.2          68  
## 4 Bibb County    Alaba…   20826   22915   22668       0.73    15.2          82.9
## 5 Blount County  Alaba…   51024   57322   58013       0.68    15.6          82  
## 6 Bullock County Alaba…   11714   10914   10309      -2.28    28.5          76.9
## # ℹ 7 more variables: multi_unit <dbl>, unemployment_rate <dbl>, metro <fct>,
## #   median_edu <fct>, per_capita_income <dbl>, median_hh_income <int>,
## #   smoking_ban <fct>

The county dataset has demographic and economic information for counties across the United States — poverty rate, unemployment, income, homeownership, education level, and more. The quesiton you are going to explore is: what predicts poverty rate at the county level?

Part 1: Explore First

Before fitting any model, look at the data. Make two exploratory graphs. Use at least one color or fill aesthetic and at least one facet_wrap() across the two graphs.

# Per-capita income vs poverty, colored by whether the county is a metro area
ggplot(data = county, aes(x = per_capita_income, y = poverty, color = metro)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_manual(values = c("no" = "#FC8D62", "yes" = "#4C9BE0"),
                     na.value = "grey70") +
  labs(title = "Per-capita income vs. poverty",
       x = "Per-capita income ($)",
       y = "Poverty rate (%)",
       color = "Metro area") +
  theme_minimal()

# Unemployment vs poverty, split into one panel per education level
ggplot(data = filter(county, !is.na(median_edu)),
       aes(x = unemployment_rate, y = poverty)) +
  geom_point(alpha = 0.25, color = "#66C2A5") +
  geom_smooth(method = "lm", se = FALSE, color = "grey30") +
  facet_wrap(~ median_edu) +
  labs(title = "Unemployment vs. poverty, by median education level",
       x = "Unemployment rate (%)",
       y = "Poverty rate (%)") +
  theme_minimal()

ANSWER: What patterns did you notice? Is there a numeric variable that looks like it might predict poverty well? Is there a categorical variable that seems to split the data into meaningfully different groups?

Both graphs point the same way. Per-capita income has a strong, tight negative relationship with poverty: as income climbs, poverty falls, and the metro counties (blue) sit toward the higher-income, lower-poverty end. Unemployment tracks poverty too, but more loosely, with a lot of scatter. The facets by education show the split I care about. Poverty sits lower in the some_college and bachelors panels than in the hs_diploma panel, and the below_hs and bachelors panels barely hold any counties at all. So per_capita_income looks like the best numeric predictor, and median_edu looks like a categorical variable that sorts counties into clearly different poverty levels.

Part 2: Numeric Predictor — Fit and Interpret

Pick the numeric variable from your exploration that you think best predicts poverty and fit a linear model.

county_model <- lm(poverty ~ per_capita_income, data = county)
summary(county_model)
## 
## Call:
## lm(formula = poverty ~ per_capita_income, data = county)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.244  -2.962  -0.485   2.112  33.838 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        3.550e+01  3.443e-01  103.10   <2e-16 ***
## per_capita_income -7.483e-04  1.283e-05  -58.34   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.514 on 3138 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.5203, Adjusted R-squared:  0.5201 
## F-statistic:  3403 on 1 and 3138 DF,  p-value: < 2.2e-16

ANSWER — Interpret the output:

  1. What variable did you choose as your predictor, and why?

  2. Write a sentence interpreting the slope in context — what does a one-unit increase in your predictor actually mean for predicted poverty rate?

  3. Is the intercept meaningful here, or does it involve extrapolation? Explain.

  4. What is the R-squared? What fraction of the variation in county poverty rate does your predictor explain?

  5. Is the slope statistically significant? What does the p-value tell you in plain language?

Answers:

  1. I chose per_capita_income. In the exploration it had one of the strongest negative correlations with poverty of any numeric variable, about -0.72. Median household income edged it out slightly, but the two measure nearly the same thing, and per-capita income is a clean per-person number, so I went with it.

  2. The slope is about -0.00075. Each additional dollar of per-capita income lowers a county’s predicted poverty rate by 0.00075 percentage points. That reads better scaled up: every extra $1,000 in per-capita income predicts a drop of roughly 0.75 percentage points in the poverty rate.

  3. The intercept is about 35.5, the predicted poverty rate when per-capita income is $0. No county comes anywhere near $0 income (the lowest is about $10,500), so the intercept is pure extrapolation. It anchors the line but has no real-world meaning on its own.

  4. The R-squared is about 0.52. Per-capita income explains roughly 52% of the variation in county poverty rates. The other half is driven by everything the model leaves out.

  5. Yes. The p-value is below 2e-16, far under 0.05, so the slope is highly significant. In plain language, if income and poverty had no real relationship, we would almost never see a slope this steep just by chance. The negative link is real, not noise.

Now plot the model:

ggplot(data = county, aes(x = per_capita_income, y = poverty)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", se = TRUE)

Part 3: Check the Residuals

# per_capita_income has 2 missing values, so lm() drops those rows. I pull
# the rows the model actually used (county_model$model) so fitted, residuals,
# and the data all line up.
county_aug <- county_model$model
county_aug$fitted    <- fitted(county_model)
county_aug$residuals <- residuals(county_model)

ggplot(data = county_aug, aes(x = fitted, y = residuals)) +
  geom_point(alpha = 0.3) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed")

ANSWER: What pattern (or lack of pattern) do you see in the residual plot? Does a linear model seem like a reasonable fit?

The residuals are not a clean, formless band. There is a shallow curve: residuals sit above zero at the low and high ends of the fitted values and dip below zero in the middle, which hints the real relationship bends a little instead of running perfectly straight. The spread also widens as fitted poverty rises, a funnel shape, so the model’s error is bigger for higher-poverty counties. A handful of counties even get predicted poverty below zero, which is impossible and shows the straight line over-extrapolates for the richest counties. A linear model is a reasonable first pass and captures the main trend, but the curve and the uneven spread say it is not a perfect fit.

Part 4: Exploring a categorical predictor — Boxplots with Jitter

The county dataset includes median_edu — the median education level in each county. It has four levels: below_hs, hs_diploma, some_college, bachelors.

Before fitting a model with this variable, look at it visually first.

ggplot(data = filter(county, !is.na(median_edu)),
       aes(x = median_edu, y = poverty)) +
  geom_boxplot() +
  geom_jitter(alpha = 0.2, width = 0.2)

ANSWER:

  1. Describe the pattern across education levels. As median education increases, what happens to poverty rate? Is the relationship consistent across all four levels?

  2. Is median_edu a nominal or ordinal categorical variable? How do you know, and why does that distinction matter here?

  3. Look at the spread (height of each box and jitter cloud) across groups. Is the spread similar across all four education levels, or does it differ? What might that mean for our model?

  4. Predict before fitting: based purely on this graph, which education level do you expect R to use as the reference category? Which level do you expect to have the largest coefficient?

Answers:

  1. Poverty drops steadily as median education rises. below_hs counties sit highest, then hs_diploma, then some_college, and bachelors counties are lowest. The direction is consistent across all four levels. The step down from below_hs to hs_diploma looks huge, but below_hs holds only a couple of counties, so that particular jump is shaky.

  2. median_edu is ordinal. The four levels have a natural order (below_hs < hs_diploma < some_college < bachelors), not just different names. That matters because “more education” is a real direction, so we expect poverty to move one way as we climb the levels, and we can talk about the ranking, not only which group a county falls in.

  3. The spread differs. The hs_diploma and some_college boxes are tall with wide jitter clouds because they hold most of the counties. below_hs and bachelors have very few counties, so their spread is hard to trust. Unequal spread plus very unequal group sizes means the model’s estimates for the tiny groups will be far less reliable than for the big ones.

  4. R sorts factor levels alphabetically, and below_hs also happens to be the lowest level, so I expect below_hs to be the reference category. Since below_hs has the highest poverty and every other level sits lower, I expect all the coefficients to be negative, and bachelors, the level furthest from below_hs, to have the largest coefficient in size.

Part 5: Categorical Predictor — Fit and Interpret

county_model_edu <- lm(poverty ~ median_edu, data = county)
summary(county_model_edu)
## 
## Call:
## lm(formula = poverty ~ median_edu, data = county)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.758  -3.858  -0.758   3.056  32.742 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              29.900      4.091   7.308 3.41e-13 ***
## median_eduhs_diploma    -10.642      4.094  -2.599  0.00938 ** 
## median_edusome_college  -16.496      4.094  -4.030 5.71e-05 ***
## median_edubachelors     -19.896      4.179  -4.761 2.02e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.786 on 3136 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.2123, Adjusted R-squared:  0.2115 
## F-statistic: 281.7 on 3 and 3136 DF,  p-value: < 2.2e-16

ANSWER:

  1. Which education level is R using as the reference category? How can you tell from the output?

  2. Interpret one non-reference coefficient (slope) in plain language. What does it mean for counties where median education is hs_diploma compared to the reference group?

  3. Compare R-squared between county_model (your numeric predictor) and county_model_edu (education level). Which predictor explains more variation in poverty rate?

  4. Did the output match your predictions from Part 4? What surprised you, if anything?

Answers:

  1. The reference category is below_hs. I can tell because it is the only one of the four levels with no coefficient row in the output. Its value is folded into the intercept (about 29.9), and the other three coefficients are all measured against it.

  2. The hs_diploma coefficient is about -10.6. Counties whose median education is hs_diploma have a predicted poverty rate about 10.6 percentage points lower than below_hs counties, the reference group. So moving from the below_hs group to the hs_diploma group goes with a large drop in predicted poverty.

  3. county_model (per-capita income) has an R-squared of about 0.52. county_model_edu (education level) has an R-squared of about 0.21. Per-capita income explains far more of the variation in poverty, more than double what education level alone does. Income is the stronger single predictor.

  4. It matched. below_hs was the reference, the coefficients were all negative, and bachelors had the largest one (about -19.9), exactly as I predicted. The surprise was the reference group itself. below_hs holds only two counties, so its intercept carries a huge standard error (about ±4.1). That makes every comparison against it wobbly, even though the overall staircase pattern is clear.

Part 6: Ethics Inside the Model

You just built two models that predict county-level poverty rates from demographic variables. Before finishing, think carefully about what these models are actually doing and what they could be used for.

ANSWER:

  1. From The Signal and the Noise: Silver distinguishes between a model that describes a pattern and one that explains it. Your regression describes a relationship between education (or your chosen numeric variable) and poverty. Does it explain it? What would it take to move from description to explanation?

  2. Imagine a bank uses a model like yours to decide which counties are “too risky” for small business loans — counties with predicted poverty above a threshold are automatically denied. Explain why this could be harmful even if the model’s predictions are statistically accurate.

  3. Think back to your ethical model-building contract from this week. Identify one specific clause and evaluate whether the model you built here satisfies it.

Answers:

  1. My regression only describes. It says counties with higher income tend to have lower poverty and puts a number on how tightly the two move together. It does not explain why. Income could lower poverty, poverty could hold income down, or a third factor like the local job market, school funding, or a region’s history could drive both at once. To move from description to explanation I would need a causal design, not a single snapshot: control for confounders, use time to check that income changes come before poverty changes, or find a natural experiment, and pair that with a real mechanism for why one produces the other. A correlation frozen at one moment cannot do that on its own.

  2. Even if the predictions are statistically accurate, wiring the model into an automatic loan denial is harmful. It turns a description into a gate. Poverty is a symptom of missing capital, so pulling loans out of high-poverty counties starves the exact places that need investment and locks in the poverty the model measured, a feedback loop where the prediction helps make itself come true. It also punishes individuals for their county’s average. A creditworthy owner in a poor county gets denied for something they did not do. Accurate on average is not the same as fair to each case, and the model would launder that unfairness as neutral math.

  3. One clause in my contract was that I would not let a model make a high-stakes decision about a person or group without a human review and a real way to appeal. The model I built here fails that clause if it is used the way question 2 describes. On its own it is just a description, which is fine. Bolted onto an automatic deny rule with no human check and no appeal, it breaks exactly the clause I wrote. The model is not unethical by itself. The use is, and my contract was written about use, not about the math.