Question

Review: Make the model above, and save it as an object named m! - Using get_regression_table(), as the percentage of Black residents increases by 1%, holding all else constant, how much is food_env_index projected to increase?

• As the percentage of Hispanic residents increases by 1%, how much does food_env_index increase? Are these changes statistically significant?

• Using get_regression_summaries(), how much variation does this model explain?

Answer

m <- counties %>%
  lm(formula = food_env_index ~ pop_black + pop_hisplat +
       pop + median_income + democrat_2016) 

# Get coefficients
m %>% get_regression_table()

# Get r-squared
m %>% get_regression_summaries()

Yes! It looks like the Food Environment Index decreases by -0.027 points for every additional percent of Black residents in that county, and it decreases by -0.006 points for every additional percent of Hispanic residents in that county. Both effects are extremely statistically significant (p < 0.001). This model explains 43.7% of the variation in Food Environment Index scores, despite including only 5 variables. Not bad!

Note: Other factors certainly matter here, but for our purposes, this is a descriptive study, meant to show how much key demographic trend with this outcome.

Question

Please recalculate the expected Food Environment Index scores for a city, using setx() and sim(), but this time, please simulate for an average city with 10% Black residents, not 9%, saving your results as zsim1, zresults1, and zev1. Please fill in the blanks below!

My model simulates that an average county where 9% of residents are Black will see a median expected Food Environment Index of 7.47 points, where 95% of those expected values lie between 7.44 and 7.5 points. But, an average county where 10% of residents are Black will see a median expected Food Environment Index of [MEDIAN], where 95% of those expected values lie between [LOWER] and [UPPER]. In other words, the median expected index increases from 7.47 to [MEDIAN], a difference of [MEDIAN - 7.47].

Answer

First, we simulate for an average city with 10% Black residents.

zsim1 <- z %>%
  setx(
    # changed just pop_black from 9 to to 10
    pop_black = 10) %>%
  sim()

Then, we extract the expected values into a data.frame, zresults1.

zresults1 <- zsim1 %>%
  zelig_qi_to_df() %>%
  select(pop_black, ev = expected_value, pv = predicted_value)

We summarize() them to get our median expected value and confidence interval.

zev1 <- zresults1 %>% 
  summarize(lower = quantile(ev, probs = 0.025),
            median = quantile(ev, probs = 0.50),
            upper = quantile(ev, probs = 0.975))

zev1

My model simulates that an average county where 9% of residents are Black will see a median expected Food Environment Index of 7.47 points, where 95% of those expected values lie between 7.44 and 7.5 points. But, an average county where 10% of residents are Black will see a median expected Food Environment Index of [7.44 ], where 95% of those expected values lie between [7.41] and [7.47`]. In other words, the median expected index increases from 7.47 to [7.44], a difference of [-0.03].

If we were to visualize our two data.frames, zev and zev1, the confidence intervals would look like this.

Question

Please simulate what happens to the food environment index in an average county given an increasingly large Hispanic/Latinx population (10, 20, 30, 40, 50, and 60%), using zelig_qi_to_df() and use geom_ribbon() to plot your confidence interval. Be sure to summarize() your data to get your confidence interval, using group_by(). Save your data.frame as zev3.

Answer

Generate your confidence interval with the code below:

zev3 <- z %>%
  setx(pop_hisplat = c(10, 20, 30, 40, 50, 60)) %>%
  sim() %>%
  # Extract data.frame
  zelig_qi_to_df() %>%
  # Grab just our variables of interest
  select(pop_hisplat, ev = expected_value) %>%
  # For each value of pop_hisplat we supplied,
  group_by(pop_hisplat) %>%
  # Get the median, lower, and upper confidence intervals
  summarize(lower = quantile(ev, probs = 0.025),
            median = quantile(ev, probs = 0.50),
            upper = quantile(ev, probs = 0.975))

Then visualize it as a ribbon plot below.

zev3 %>%
  ggplot(mapping = aes(x = pop_hisplat, y = median, ymin = lower, ymax = upper)) +
  # Add a transparent ribbon, with a black outline and blue fill
  geom_ribbon(alpha = 0.5, color = "black", fill = "seagreen") +
  # Add a dashed line to plot the median
  geom_line(linetype = "dashed") +
  # Add a theme and labels
  labs(x = "% Hispanic/Latinx Residents", y = "Expected\nFood Environment Index",
       subtitle = "Change in Expected Index\nwith 95% Confidence Intervals",
       caption = "From 1000 simulations for an average county.") +
  theme_bw()

Question

What does mydist represent, and why does it matter whether its coefficients are correlated?

Answer

Our object mydist contains not just one set of model coefficients, but 1000 slightly varied versions of our model coefficients. These values were randomly drawn from a multivariate normal distribution, and that kind of distribution is special because it makes sure that simulated coefficients (like pop_black and democrat_2016) keep the same correlations and patterns from our original model.

Definitions

• Estimation Uncertainty: the uncertainty in our predictions due to not knowing Beta and Alpha perfectly, because we had fewer than an infinite number of observations in our data. Our multivariate normal distribution helps us account for this.

• Fundamental Uncertainty: the uncertainty in our predictions due to innumerable chance events, like weather, illness, people spilling their coffee, banana peels, etc. that may influence our outcome but are not contained within our predictors. Random samples from a normal distribution help us account for this.

Examples

• Estimation Uncertainty: If our coefficients were Marvel superheros, estimation uncertainty is akin to the feeling that I’m pretty sure the superheroes will team up and defeat the bad guy, but I’m not exactly sure how. We know that Thor will go do something hotheaded, and we know that the Hulk will SMASH, but we’re not exactly sure whether that means a little SMASH or a lot of SMASH, for example. So, we are uncertain on exactly what will happen, but we kind of already know what is going to happen when we go to the movies.

• Fundamental Uncertainty: If our coefficients were Marvel superheros, fundamental uncertainty is akin to the butterfly effect or the series What If, where tiny changes alter the outcome in innumerable, unpredictable ways.

Question

In your data.frames myev and mypv, you simulated 1000 predicted values and 1000 expected values, and calculated the 95% confidence intervals for each. Which confidence interval was wider? Now compare that to the table above. Generally, does the confidence interval of predicted values tend to be wider or narrower than the confidence interval of expected values?

Answer

In our example above with zev and zpv, the confidence interval for our predicted values (5.772 to 9.17) was much wider than the confidence interval for our expected values (7.436 to 7.5). This matches the table above, which says predicted values produce wider confidence intervals than expected values.

Coding Estimation Uncertainty

Now, we’re going to make 1000 predictions from our 1000 simulated coefficients! These will be saved into the pv column in the object zsim. We will make them by plugging our predictor data x into each of our 1000 slightly different model equations using our 1000 slightly different simulated coefficients, to get our 1000 slightly different predictions. These values reflect estimation uncertainty.

Model Equation: \(Y_{simulated} = Alpha + Beta_{1} \times X_{1.observed} + Beta_{2} \times X_{2.observed} + ...\)

diysim <- mydist %>%
  mutate(
    # Use our simulated coefficients, 
    # saved in mydist under each variable name
    ysim = intercept + 
      # Multiply each by single value for variables in x
      pop_black * x$pop_black + 
      pop_hisplat * x$pop_hisplat +
      pop * x$pop + 
      median_income * x$median_income + 
      democrat_2016 * x$democrat_2016)

diysim %>% # reorder variables
  select(id, ysim, intercept:democrat_2016) %>%
  head(3)

For example, here are 3 of our simulated outcomes (lines), compared with all 1000 of them (histogram).

Coding Fundamental Uncertainty

Next, we can account for fundamental uncertainty using the rnorm() function.

rnorm() samples 1 or more values from a hypothetical distribution of values. To make that hypothetical distribution, we must give it a center (mean) and a standard deviation to estimate its width (sd). Our distribution will have a mean equal to our simulated value ysim, and a sd equal to sigma, the average variation of our model residuals. This distribution reflects the full range of simulated values we might get due to fundamental uncertainty. Some are larger, some are smaller, but most are still close to our simulated value. We can sample just 1 value from this distribution (n = 1) or many values, eg. 1000 (n = 1000).

• To get the predicted value, we will randomly draw one value from our distribution approximating fundamental uncertainty, giving us our predicted value pv.

• To get the expected value, we will randomly draw 1000 values from our distribution approximating fundamental uncertainty, and calculate the mean of those 1000 values, giving us our expected value ev.

We will save them both in a data.frame called diyresults.

diyresults <- diysim %>% 
  # For each of our 1000 simulated outcomes,
  # which already accounted for estimation uncertainty,
  group_by(id) %>%
  summarize(
    # Record the value we used to calculate everything (9%) 
    pop_black = x$pop_black,
    # Calculate the predicted value by drawing 1 value
    # from our normal distribution
    # to represent fundamental uncertainty
    pv = rnorm(n = 1, mean = ysim, sd = sigma),
    # Calculate the expected value by drawing 1000 values
    # from our normal distribution
    # and average them to eliminate fundamental uncertainty
    ev = rnorm(n = 1000, mean = ysim, sd = sigma) %>% mean())

Compare

Finally, let’s compare our version against Zelig’s version, by stacking our do-it-yourself data.frame of projections diyresults on top of our Zelig-made data.frame, and visualizing the distributions of these simulations. We see that our estimates line up nearly perfectly! Our expected values have a very thin spread, because we averaged to get rid of fundamental uncertainty, while predicted values have very wide spread, to reflect fundamental uncertainty.

compare <- bind_rows(
  # Get our version
  diyresults %>% # labelling it by 'type'
    mutate(type = "Do-It-Yourself"),
  # Get Zelig's version
  zresults %>% # labelling it by 'type'
    mutate(type = "Zelig") )

compare %>%
  pivot_longer(cols = c(ev, pv), 
               names_to = "measure", values_to = "value") %>%
  ggplot(mapping = aes(x = type, y = value, color = type)) +
  geom_jitter(alpha = 0.5, shape = 21, size = 3, fill = "white") +
  geom_violin(alpha = 0.75, color = "black") +
  facet_wrap(~measure) +
  # Add labels
  labs(x = "Type", y = "Projected\nFood Environment Index",
       subtitle = "Equivalent Results",
       caption = "From 1000 simulations of average county\ngiven 9% Black Residents") +
  scale_color_viridis(option = "mako", discrete = TRUE,
                       begin = 0.2, end = 0.8) +
  theme_bw() +
  guides(color = "none")

Question

What does the variable ysim in our data.frame diysim refer to? How did we calculate it?

Answer

• The variable ysim refers to 1000 different simulated outcomes.

• We calculated it by taking 1000 slightly different model equations, each a row in mydist, plugging in the values of our predictors from x (eg. pop_black = 9), and calculating the output of each of these equations. This produced our 1000 different simulated outcomes, ysim.

• Their variation reflects estimation uncertainty, or how our model coefficients could be slightly off due to chance.

Question

When we calculate a single predicted value, how many values do we draw from a normal distribution? Why? How about when we calculate a single expected value? How does this help us account for uncertainty?

Answer

• When we calculate a single predicted value, we draw just 1 value from a normal distribution using rnorm(). This adds some randomness to our predicted values, helping spread them out to demonstrate fundamental uncertainty.

• When we calcualte a single expected value, we draw 1000 value from a normal distribution using rnorm(), and then take their average. This approximates the fundamental uncertainty of our estimate, then lets us get rid of it by finding the most likely, most expected of our predictions, the average.

Question

Reproduce the visual below! Please simulate and visualize the change in expected Food Environment Index scores as the share of Black residents increases from 10 to 30 to 50 percent in an otherwise average county with a median income of 40,000 USD (~25th percentile), 50,000 USD (~50th percentile), and 60,000 USD (~75th percentile). What trends do you find?

Answer

• Results ($40,000): We observe that as the share of Black residents increases from 10% to 50% in an average county with a median income of $40,000, our model simulates a decrease in expected Food Environment Index scores from a median of ~7 to ~5.8 points.

• Results ($50,000): However, as the median income increases by 10,000 USD to 50,000 USD, the expected Food Environment Index score also increases, by roughly 0.40 points.

• Results ($60,000): Given 60,000 USD, an otherwise average county still sees a decrease in expected Food Environment Index scores as the share of Black residents increases from 10% to 50%, but it now decreases from a median of ~7.8 to ~6.6 points.

• Interpretation: This shows that while racial discrimination’s long term effects impact access to healthy foods, income helps mediate this effect, meaning that poorer communities are especially at risk. See code below!

lc8 <- z2 %>%
  setx(pop_black = c(10, 30, 50,  
                     10, 30, 50,
                     10, 30, 50),
       median_income = c(40000, 40000, 40000,
                         50000, 50000, 50000,
                         60000, 60000, 60000)) %>%
  sim() %>%
  zelig_qi_to_df() %>%
  select(pop_black, median_income, ev = expected_value) %>%
  # 
  group_by(pop_black, median_income) %>%
  summarize(lower = quantile(ev, probs = 0.025),
            median = quantile(ev, probs = 0.50),
            upper = quantile(ev, probs = 0.975))

lc8 # View it!

lc8 %>%
  # Recode, changing it to a categorical variable
  mutate(median_income = median_income %>% recode_factor(
    "40000" = "$40,000",
    "50000" = "$50,000",
    "60000" = "$60,000")) %>%
  ggplot(mapping = aes(x = pop_black, y = median, 
                       ymin = lower, ymax = upper, 
                       fill = median_income)) +
  # Add Confidence Intervals
  geom_ribbon(alpha = 0.5, color = "white") +
  geom_linerange(size = 1) +
  # Add median expected value
  geom_point(size = 3) +
  # Split up by region, ordered by median
  facet_wrap(~reorder(median_income, median) ) +
  # Label median expected values
  geom_text(mapping = aes(label = round(median, 1)), 
            size = 5, nudge_x = 5) +
  # Cut fill legend
  guides(fill = "none") +
  # Add labels
  labs(x = "% Black Residents",
       y = "Expected\nFood Environment Index",
       subtitle = "Effects of Racial Discrimination and Income\non Access to Healthy Food",
       caption = "From 1000 simulations of average counties\n(Points = median; ribbons = 95% confidence intervals)") +
  scale_fill_viridis(option = "mako", discrete = TRUE,
                     begin = 0.2, end = 0.8) +
  theme_bw()