--- title: "LAB 4" --- # Instructions - You will do your work in this `Quarto` document but will submit all your answers on Canvas - *very important*: As you move through the document, make sure to run all code chunks (grey bits) that already have code in them - Write your own code in the empty code chunks - Useful shortcuts: - to run all the code in a specific code chunk, press the green right-facing triangle at the top right of the code chunk - to run all *prior* code chunks, press the downward-facing gray triangle at the top right of the code chunk ```{r load-libraries, warning=FALSE, message=FALSE} library(tidyverse) library(broom) library(AER) library(juanr) library(moderndive) set.seed(1990) ``` # Task 1: GDP Simulate data that shows corruption hurts a country's economic development, with the following characteristics: - You have n = 200 countries - Each country has a corruption level of 50, plus or minus 10 (0 = no corruption, 100 = high corruption) - Each country has a GDP per capita of about 50,000, plus or minus 8,000 - corruption causes a decrease in GDP per capita of 1,000 per unit of corruption ```{r} set.seed(1990) gdp_data = tibble( corruption = rnorm(n = 200, mean =50, sd = 10), gdp = rnorm(n = 200, mean = 50000, sd = 8000) - 1000*corruption, ) ``` Fit a model that estimates the effect of corruption on GDP per capita: ```{r} lm(gdp ~ corruption, data = gdp_data) |> tidy() ``` **Q1: What is the estimated coefficient for the intercept in this model?** Simulate data that shows that: terrorist attacks increase the public's sense of anxiety; and that this increase in anxiety, leads to higher support for more extreme political candidates. Here are the characteristics: - it's a country that has N = 800 cities - there's a lot of terrorist attacks happening. the cities experience, on average, 8 terrorist attacks, plus or minus 1 attack - people in this country are anxious. On a 100-point "anxiety scale", the average city is at 20, +/- 2. - terrorism causes *increases* in anxiety: for every terrorist attack, anxiety increases by 2 points. - Extreme political candidates are not so popular. The average vote share of an extreme political candidate is roughly 15, +/- 3. - as cities get more anxious, support for these candidates *increases*. Say that for every point of anxiety, support for extreme candidates goes up by 0.5 points. ```{r} set.seed(1990) terrorism_data = tibble( terrorist_attack = rnorm(n = 800, mean = 8, sd = 1), anxiety = rnorm(n = 800, mean = 20, sd = 2) + 2*terrorist_attack, ext_cand = rnorm(n = 800, mean = 15, sd = 3) + 0.5*anxiety ) ``` Now, make a scatterplot of the relationship between the number of terror attacks a city has and the level of support for extreme candidates in that city. Color the points by `anxiety` level. Include an OLS line. ```{r} ggplot(terrorism_data, aes(x=terrorist_attack, y=ext_cand, color = anxiety)) + geom_point() + geom_smooth(method= "lm") ``` **Q2: right click the graph --\> "save image as" --\> and save it somewhere on your computer. You will need to submit the plot!** Now, using a model, what is the effect of terrorism on support for extreme political candidates? ```{r} lm(ext_cand ~ terrorist_attack, data = terrorism_data) |> tidy() ``` **Q3: what is the estimated effect of terrorism on support for extreme candidates?** **Q4: would it be a good idea, in this case, to control for anxiety? Why or why not?** Simulate data that shows that the relationship between education and ideology is confounded by class. Here are the characteristics: - You have n = 500 people - Most people have 18 years of education, plus or minus 2 (education) - People are either high income, or low income, equally likely to be one or the other (class) - Most people are perfectly centrist on a (0-100 ideological scale, 0 = conservative, 100 liberal), plus or minus 5 (ideology) - Every additional year of education makes someone more liberal, by about 1 point - People who are high income get, on average, about 5 years more of education - People who are high income are, on average, about 10 points more liberal ```{r} set.seed(1990) income_data = tibble( income = sample(c(0, 1), size = 500, replace = TRUE), education = rnorm(n = 500, mean = 18, sd = 2) + 5 * income, ideology = rnorm(n = 500, mean = 50, sd = 5) + 1*education + 10*income ) ``` Fit a model that estimates the effect of education on ideology. ```{r} lm(ideology ~ education, data = income_data) ``` **Q5: about how much larger is your estimated effect than the true effect? To answer, calculate = (estimated effect / true effect) X 100** 2 / 1 \* 100 200 times more # Task 2: Affairs For this part, you'll look at data on extra-marital affairs, taken from a Psychology Today study in 1969, from the AER package. The affairs variable is coded weirdly, but just pretend it counts the number of affairs the respondent has had (or admits to having). ```{r load-data} data("Affairs") ?Affairs # read about the data and variables head(Affairs, n = 10) ``` First, a blast from the past: find the percent of respondents who had no affairs, one affair, two affairs, and so on: ```{r} affairs_update = Affairs |> group_by(affairs) |> tally() |> mutate(percent = n/sum(n) * 100) affairs_update ``` **Q6: Roughly what percent of respondents have had at least one affair?** Now the key question: does having a child make someone more or less likely to have an extra-marital affair? Estimate the naive model, where you control for nothing: ```{r} lm(affairs ~ children, data = Affairs) |> tidy() ``` **Q7: about how many more affairs do people with children have, on average?** There's good reason to think how long a couple has been married is a **fork** here, since people who have been married longer are more likely to have kids AND more likely to have affairs. Fit a model where you properly adjust for the **fork**: ```{r} lm(affairs ~ children + yearsmarried, data = Affairs) |> tidy() ``` **Q9: how does the estimated effect of children on affairs change?** **Q8: Explain why length of marriage is a fork here?** Finally, describe why a couple's overall happiness might be a **pipe** that we should avoid controlling for. **Q10: How could happiness be a pipe?** **Q11: How could therapy be a collider?** # Task 3: Feelings Look at the feeling thermometer data. Just like in class, we'll pretend this data is the whole of the United States. ```{r} therm ?therm ``` What is the average feeling thermometer for the police in the whole dataset (`therm`)? Calculate the average below: ```{r} therm_police <- therm |> summarize(avg_police = mean(ft_police, na.rm = TRUE)) therm_police ``` **Q12: what is the average level of support for the police?** Now, let's look at how much our estimates might vary across samples. Take 1,000 samples each of size 30 and calculate the average support for police in each. ```{r} set.seed(1990) sample_police = therm |> rep_sample_n(size = 30, reps = 1000, replace = TRUE) |> summarize(avg_police = mean(ft_police, na.rm = TRUE)) min(sample_police$avg_police) ``` **Q13: what is the lowest average level of support for the police that you estimated in your samples?** Now, bias the samples in a form of your choosing so that the sample estimates do not converge on the population parameter. Plot the resulting distribution. ```{r} sample_police = therm |> filter(race != "White") |> rep_sample_n(size = 30, reps = 1000, replace = TRUE) |> summarize(avg_police = mean(ft_police, na.rm = TRUE)) ggplot(sample_police, aes(x = avg_police)) + geom_histogram() ``` **Q14: right click the graph --\> "save image as" --\> and save it somewhere on your computer. You will need to submit the plot!** # Task 4: Bootstrapping Now, let's go back to treating `therm` as what it is: a sample. How much uncertainty should we have in our estimate of average police support in this sample? Generate 1,000 bootstrapped samples of the dataset. Between what two values do 95% of observations fall? ```{r} set.seed(1990) therm |> rep_sample_n(size = nrow(therm), reps = 1000, replace = TRUE) |> summarise(avg_police = mean(ft_police, na.rm = TRUE)) |> summarise(low = quantile(avg_police, .025), # middle 95% means lower bound is .025 mean = mean(avg_police, na.rm = TRUE), high = quantile(avg_police, .975)) # middle 95% means upper bound is .975 ``` **Q15: Between what two values do 95% of averages fall?** Repeat the same procedure above, except this time only for respondents who identify as "Middle Eastern". ```{r} set.seed(1990) therm |> filter(race == "Middle Eastern") |> rep_sample_n(size = nrow(therm), reps = 1000, replace = TRUE) |> summarise(avg_police = mean(ft_police, na.rm = TRUE)) |> summarise(low = quantile(avg_police, .025), # middle 95% means lower bound is .025 mean = mean(avg_police, na.rm = TRUE), high = quantile(avg_police, .975)) # middle 95% means upper bound is .975 ``` **Q16: Between what two numbers do 95% of averages fall?** **Q17: Why is the range of values in Q15 different than in Q16?** Look at the output from the following regression: ```{r} lm(ft_police ~ sex, data = therm) %>% broom::tidy() ``` **Q18 Roughly speaking, what does the value in `std.error` for `sexFemale` tell us?**