---
title: "Week 8 Tutorial"
subtitle: "Data Analysis"
format:
  html:
    toc: true
    toc-location: left
    toc-depth: 3
    highlight-style: github
    code-copy: true
    embed-resources: true
    df-print: paged
    include-in-header:
      text: |
        <style>
          .cell-output pre { max-height: 250px; overflow-y: auto; }
        </style>
editor: visual
---

First things first: Make a new project in RStudio, *and be sure to place it somewhere where you can easily find it* (e.g., Documents, Desktop). Then download this file (the `.qmd` file) from the syllabus website and place it in the folder associated with that project. Then open this file in RStudio and follow along.

I encourage you to work with/talk to others throughout this tutorial, especially during times where I have you code on your own and use AI. When you code, think about the SPARCK framework:

**S: Set up** (What is your end goal?)

**P: Pseudocode** (What are the steps?)

**A: AI** (What are your AI prompts?)

**R: Run** (Run it!)

**C: Check** (Did your code work?)

**K: Knowledge** (What does your code mean?)

## Agenda

- Descriptive statistics

- Visualizing distributions

- t-tests

- Correlations

- Linear regression

- Practice

In the first few weeks of the course, we learned coding basics: what a dataframe is, how to read in data, and how to clean and visualize it. Then we learned about different kinds of data and different ways of getting data, like maps, text, and web scraping. Now we turn to analyzing data: describing it, and testing whether the patterns we see are meaningful.

This is not a stats class, so we won't go deep into the math behind each test.

## Descriptive statistics

When we get a new dataset, the first thing we want to do, before looking at relationships between variables, is understand what our data looks like. We want to be able to *describe* it before we can make *inferences* about it.

### Getting started

Today's data comes from Opportunity Insights, a research group at Harvard that publishes lots of data about income and opportunity in the US. We'll use data from a study that linked tax records to death records for millions of Americans from 2001 to 2014 (Chetty et al., 2016). It has the life expectancy at age 40 (the age a 40-year-old is expected to live to, on average) for men and women in each state and year, broken down by income quartile. Quartile 1 (Q1) is the poorest 25% of people in the US, and Quartile 4 (Q4) is the richest 25%.

Go to the Opportunity Insights [data library](https://opportunityinsights.org/data/), and find the paper called "The Association Between Income and Life Expectancy in the United States, 2001-2014." Under that paper, find the dataset called "State Life Expectancy Estimates by Year, Gender, and Income Quartile." Download the CSV version, and put it in your project folder. Also open the README next to it, which explains what each column means.

Then load `tidyverse`, and read in the data.

```{r}
library(tidyverse)
le_raw <- read_csv('health_ineq_online_table_5.csv')
```

Take a look at the data, and compare the column names to the README. The column names follow a pattern. For example, `le_agg_q1_F` is the life expectancy (`le`) for the poorest quartile (`q1`) of women (`F`). There are also columns that start with `sd_` (the standard error of each estimate), `count_` (how many people each estimate is based on), and `le_raceadj` (estimates adjusted for differences in racial makeup between states). We'll only use the unadjusted `le_agg` columns today.

Create a new dataframe with only the state name and abbreviation, the year, and the `le_agg` columns.

```{r}
le <- le_raw |>
  select(stateabbrv, statename, year, starts_with('le_agg'))
```

### Basic summaries

There are a couple of functions we can use right off the bat to make sense of our data. Use `str()` to see the structure of `le`. It's not the prettiest output, but it shows every column, the type of data in each column, and the first few values.

```{r}
str(le)
```

Another useful function is `summary()`, which gives you the minimum, maximum, mean, median, and quartiles of every numerical column. Use it on `le`.

```{r}
summary(le)
```

You can also calculate each of these statistics on its own. Calculate the mean, median, standard deviation, and range (the minimum and maximum) of life expectancy for the poorest women.

```{r}
mean(le$le_agg_q1_F)
median(le$le_agg_q1_F)
sd(le$le_agg_q1_F)
range(le$le_agg_q1_F)
```

### Grouped summaries

`summary()` and `mean()` describe the whole dataset at once. Usually, it's more useful to describe groups. Use `group_by()` and `summarize()` to create a new dataframe, with the mean and SD of life expectancy for the poorest women and the richest women in each year.

```{r}
le_by_year <- le |>
  group_by(year) |>
  summarize(mean_q1_F = mean(le_agg_q1_F),
            sd_q1_F = sd(le_agg_q1_F),
            mean_q4_F = mean(le_agg_q4_F),
            sd_q4_F = sd(le_agg_q4_F))

le_by_year
```

Summary tables are great for graphing. Make a line graph with year on the x axis, and two lines: one for the mean of the poorest women, and one for the mean of the richest women.

```{r}
ggplot(data = le_by_year, aes(x = year)) +
  geom_line(aes(y = mean_q1_F), color = 'red') +
  geom_line(aes(y = mean_q4_F), color = 'blue') +
  labs(x = 'Year', y = 'Mean life expectancy at age 40')
```

Now create a dataframe with the mean life expectancy for the poorest men in each state, across all years. Then sort it from highest to lowest.

```{r}
le_by_state <- le |>
  group_by(statename) |>
  summarize(mean_q1_M = mean(le_agg_q1_M)) |>
  arrange(desc(mean_q1_M))

le_by_state
```

## Visualizing distributions

Before we run any statistical tests, it's useful to visualize the spread of our data with a histogram.

Use `ggplot` to make histograms of life expectancy for the poorest men and the richest men, on the same graph. Since they're in different columns, you'll need to add two `geom_histogram()` layers, each with its own `aes(x = ...)`. Make each bar 0.5 years wide with the `binwidth` argument, give each layer a different fill color, and make them partially transparent (using `alpha =`) so you can see where they overlap.

```{r}
ggplot(data = le) +
  geom_histogram(aes(x = le_agg_q1_M), binwidth = 0.5, fill = 'red', alpha = .5) +
  geom_histogram(aes(x = le_agg_q4_M), binwidth = 0.5, fill = 'blue', alpha = .5)
```

## T-tests

The t-test is one of the simplest tools we have to make *inferences* about our data. **A t-test tells us whether the difference between two means is meaningful, or "significant."** All t-tests in R use the same function: `t.test()`.

Essentially, a t-test tells us how likely it is that a difference between two groups is simply due to random noise, or is real. This probability is called the p-value. Lower p-values mean it's more likely that your result is "real." In psychology, we use p = 0.05 as a cutoff for significance.

### Independent samples t-tests

Did life expectancy for the poorest men change between 2001 and 2012? Create two new dataframes, with the data from each year.

```{r}
le01 <- le |>
  filter(year == 2001)

le12 <- le |>
  filter(year == 2012)
```

To compare two columns, give `t.test()` each column as its own argument. Run a t-test that compares life expectancy for the poorest men in 2012 to 2001.

```{r}
t.test(le12$le_agg_q1_M, le01$le_agg_q1_M)
```

Let's read the output from the top:

- t is the t value, and df is the degrees of freedom (which is related to how big our sample is).

- p-value is the p-value.

- The 95 percent confidence interval is a range of plausible values for the *real* difference between the means. If it includes 0, then "no difference" is still plausible.

- The sample estimates are the means of the two groups.

Is the difference significant?

**Answer:** No. The mean is 0.72 years higher in 2012, but the p-value is .056, which is above .05. The confidence interval (-0.02 to 1.46) includes 0.

### Paired samples t-tests

What we just ran is called an **independent samples t-test**. It treats the two columns as if they came from two completely separate groups, like in a between-subjects design (e.g., one group of people gets a treatment, and a different group doesn't).

But that's not our data. The *same* 51 states were measured in both years, so this is more like a within-subjects design. In a within-subjects design, the right test is a **paired samples t-test**, which compares each state *to itself*. Run the same t-test, but add the argument `paired = TRUE`.

```{r}
t.test(le12$le_agg_q1_M, le01$le_agg_q1_M, paired = TRUE)
```

Is the difference significant now?

**Answer:** Yes. The t value went from 1.93 to 4.28, and the p-value went from .056 to .00008.

### One-sample t-tests

The simplest kind of t-test is a one-sample t-test. Instead of comparing two groups, it asks if the mean of one column is different from 0.

Our data has one row for each state in each year, so California shows up 14 times. To keep things simple, we'll use just one year. Create a new dataframe with just the data from 2014. Then add a column that represents the difference in life expectancy between the richest and poorest women in each state. Then run a one-sample t-test using this new column.

```{r}
le14 <- le |>
  filter(year == 2014) |>
  mutate(gap_F = le_agg_q4_F - le_agg_q1_F)

t.test(le14$gap_F)
```

Now run a paired t-test that compares the richest women to the poorest women. Compare its output to the one-sample t-test.

```{r}
t.test(le14$le_agg_q4_F, le14$le_agg_q1_F, paired = TRUE)
```

They're exactly the same - a paired t-test is really just a one-sample t-test on the differences between the pairs.

### T-tests with formulas

What if we wanted to compare life expectancy for the poorest men in New York and in Texas?

This sounds simple, but look at how our data is set up. Our earlier t-tests compared two columns. But here, we don't have a New York column and a Texas column. Instead, all of the life expectancy values are in one column, and another column (`stateabbrv`) tells us which group each row belongs to.

For data set up this way, `t.test()` uses formula notation, which looks like `y ~ x`. This literally asks: How does X *affect* Y?

Create a new dataframe with only the rows for New York and Texas, from every year. Then make a box plot that compares life expectancy for the poorest men in the two states.

```{r}
ny_tx <- le |>
  filter(stateabbrv %in% c('NY', 'TX'))

ggplot(data = ny_tx, aes(x = stateabbrv, y = le_agg_q1_M)) +
  geom_boxplot()
```

Now run a t-test with formula notation, where life expectancy for the poorest men is the Y variable and state is the X variable. (Hint: When you use formula notation, you also need the `data =` argument.)

```{r}
t.test(le_agg_q1_M ~ stateabbrv, data = ny_tx)
```

This is an independent samples t-test, because New York and Texas are two separate groups. (The title says "Welch Two Sample t-test." Welch's t-test is the version of the independent samples t-test that R uses by default.)

## Correlations

A t-test compares a variable across two groups. A **correlation** tells us how strongly two continuous variables are related. The correlation coefficient, *r*, goes from -1 to 1:

- *r* = 1 means a perfect positive relationship (when one goes up, so does the other).

- *r* = 0 means no relationship.

- *r* = -1 means a perfect negative relationship (when one goes up, the other goes down).

Use `cor.test()` to calculate the correlation between life expectancy for the poorest women and the poorest men in 2014, and to find out if the correlation is significant.

```{r}
cor(le14$le_agg_q1_F, le14$le_agg_q1_M)
cor.test(le14$le_agg_q1_F, le14$le_agg_q1_M)
```

It often helps to visualize a correlation. Make a scatterplot with a line of best fit.

```{r}
ggplot(data = le14, aes(x = le_agg_q1_F, y = le_agg_q1_M)) +
  geom_point() +
  geom_smooth(method = 'lm')
```

You can also use `cor()` to create a correlation matrix, where you calculate the correlation between every pair of columns. Make a correlation matrix of all of the `le_agg` columns in `le14`.

```{r}
le14 |>
  select(starts_with('le_agg')) |>
  cor()
```

## Linear regression

Here's a secret. Almost all of the statistical tests you learn about (t-tests, correlations, ANOVAs) can be done with linear regression.

A linear regression asks: **As my X value changes, how much does my Y value change?** It draws the line of best fit between X and Y (like the one in our scatterplot), and tells us how steep it is. Steeper slopes mean bigger effects. We'll learn more about linear regressions in our next class.

In regression, X is often called the predictor (or independent variable), and Y is called the outcome (or dependent variable). A line has two parts:

- The intercept: the value of Y when X is 0 (where the line crosses the y axis).

- The slope: how much Y changes when X goes up by 1.

A linear regression estimates both of these, and then asks whether the slope is significantly different from 0 (a flat line, meaning X has no effect on Y).

### Simple linear regression

To run a linear regression, we use the `lm()` function (which stands for linear model). `lm()` uses formula notation, just like our last `t.test()`. On its own, `lm()` just shows us how the model is set up. To see the statistics, we have to wrap the whole thing in `summary()`. Copy your New York vs. Texas t-test, change `t.test` to `lm`, and wrap it in `summary()`.

```{r}
summary(lm(le_agg_q1_M ~ stateabbrv, data = ny_tx))
```

The most important part of the output is the Coefficients table. It has one row for the intercept, and one row for each predictor.

When X is a category, R picks one group to be the reference group (by default, whichever comes first alphabetically, so here it's NY). It codes the reference group as 0 and the other group as 1. So:

- The `(Intercept)` row is the mean of the reference group (New York): 78.30 years.

- The `stateabbrvTX` row is the slope: how much Y changes when X goes from 0 (New York) to 1 (Texas). R names the row after the column (`stateabbrv`) plus the group being compared to the reference group (`TX`).

Each row has four columns:

- The `Estimate` is the value of the intercept or slope. Here, the slope tells us that the life expectancy of the poorest men is 1.25 years lower in Texas than in New York (78.30 - 1.25 = 77.04, the mean for Texas).

- The `Std. Error` tells us how uncertain the estimate is. Smaller standard errors mean more certain estimates.

- The `t value` is the estimate divided by its standard error, just like in a t-test.

- The `Pr(>|t|)` column is the p-value. For the slope, it tells us whether the slope is significantly different from 0.

Look back at the t-test you ran with formula notation. The difference between the means is the same as the slope here. A regression with one categorical predictor is basically a t-test!

What's nice about linear regression is that it works the same way when X is continuous. Run a linear regression that asks whether life expectancy for the poorest women changed over time, using all of the states and years in your data.

```{r}
summary(lm(le_agg_q1_F ~ year, data = le))
```

Now that X is continuous, "one unit" means one year. So the slope (the Estimate in the `year` row) tells us that, on average, the life expectancy of the poorest women went up by about 0.07 years (less than a month) each year.

### Multiple regression

There are usually *many* things that affect Y. With linear regression, we can put multiple predictors in the same model, to see the effect of each one. There are two ways to add predictors:

- With `+` (e.g., `y ~ x1 + x2`), the model estimates the effect of each predictor while holding the other one constant. This is how researchers "control for" other variables.

- With `*` (e.g., `y ~ x1 * x2`), the model also checks whether the predictors interact: whether the effect of one predictor *depends on* the value of the other. For example, does the change over time depend on which state you're in?

Run a linear regression that asks whether the change over time in life expectancy for the poorest men is different in New York and Texas. Put both `year` and `stateabbrv` in the model as interacting predictors, using `*` between them.

```{r}
summary(lm(le_agg_q1_M ~ year * stateabbrv, data = ny_tx))
```

Ignoring the Intercept, the Coefficients table now has 3 rows, one for each effect: two main effects (`year` and `stateabbrvTX`), and their interaction (`year:stateabbrvTX`).

Interactions can be tricky to read. When a model has an interaction, each main effect is the effect when the *other* predictor is 0:

- The `year` row is the slope over time for the reference group (New York, where `stateabbrv` is 0).

- The `stateabbrvTX` row is the difference between Texas and New York in the year 0. That's way outside our data, which is why the number looks so strange. (Again, we'll learn how to fix this in Week 10.)

- The `year:stateabbrvTX` row is the interaction. It tells us how *different* the slope over time is for Texas than for New York. If it's significant, the two states changed at different rates.

What is the slope for New York? What is the slope for Texas?

**Answer:** The slope for New York is 0.21 years per year. The slope for Texas is 0.21 - 0.21, which is about 0 years per year.

Interactions are much easier to understand when you graph them. Let's do that using SPARCK.

**S:** Make a graph that shows the interaction between year and state.

**P:** Make a scatterplot with year on the x axis and life expectancy for the poorest men on the y axis, using the `ny_tx` data. Give each state its own color, and add a straight line of best fit for each state.

**A:** I have a dataframe in R called `ny_tx`, with columns called `year`, `stateabbrv` (either NY or TX), and `le_agg_q1_M` (life expectancy). Using ggplot, make a scatterplot with `year` on the x axis and `le_agg_q1_M` on the y axis, with the points colored by `stateabbrv`. Add a straight line of best fit for each state.

**R:**

```{r}
ggplot(data = ny_tx, aes(x = year, y = le_agg_q1_M, color = stateabbrv)) +
  geom_point() +
  geom_smooth(method = 'lm')
```

**C:** There are 14 points for each state, one for each year from 2001 to 2014. The New York line goes up, and the Texas line is flat, which matches the slopes in our regression output.

**K:** `color = stateabbrv` inside `aes()` gave each state its own color. Because the color is set in `ggplot()`, every layer uses it, so `geom_smooth(method = 'lm')` drew a separate straight line of best fit for each state. The slopes of these lines are the slopes from our regression.

## Practice

Now it's your turn. Use SPARCK to answer a question about this data. Your analysis must:

- Use at least one t-test, correlation, or linear regression that we **didn't** run in this tutorial (e.g., different columns, groups, states, or years).

- Include a graph.

For example, you could ask whether the gap between the richest and poorest women changed over time, compare two states you're interested in, or test whether the middle quartiles (Q2 and Q3) are different from each other.

**S:**

**P:**

**A:**

**R:**

```{r}

```

**C:**

**K:**

## References

Chetty, R., Stepner, M., Abraham, S., Lin, S., Scuderi, B., Turner, N., Bergeron, A., & Cutler, D. (2016). The association between income and life expectancy in the United States, 2001–2014. *JAMA, 315*(16), 1750–1766. <https://doi.org/10.1001/jama.2016.4226>

Wickham, H., Çetinkaya-Rundel, M., & Grolemund, G. (2023). *R for data science* (2nd ed.). O'Reilly. <https://r4ds.hadley.nz/>