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, 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.
library(tidyverse)
Warning: package 'ggplot2' was built under R version 4.4.3
Warning: package 'tidyr' was built under R version 4.4.3
Warning: package 'purrr' was built under R version 4.4.3
Warning: package 'dplyr' was built under R version 4.4.3
-- Attaching core tidyverse packages ------------------------ tidyverse 2.0.0 --
v dplyr 1.2.0 v readr 2.1.5
v forcats 1.0.0 v stringr 1.5.1
v ggplot2 4.0.2 v tibble 3.2.1
v lubridate 1.9.4 v tidyr 1.3.2
v purrr 1.2.1
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Rows: 714 Columns: 44
-- Column specification --------------------------------------------------------
Delimiter: ","
chr (2): stateabbrv, statename
dbl (42): st, year, le_agg_q1_F, le_raceadj_q1_F, le_agg_q2_F, le_raceadj_q2...
i Use `spec()` to retrieve the full column specification for this data.
i Specify the column types or set `show_col_types = FALSE` to quiet this message.
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.
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.
Another useful function is summary(), which gives you the minimum, maximum, mean, median, and quartiles of every numerical column. Use it on le.
summary(le)
stateabbrv statename year le_agg_q1_F
Length:714 Length:714 Min. :2001 Min. :75.44
Class :character Class :character 1st Qu.:2004 1st Qu.:81.51
Mode :character Mode :character Median :2008 Median :82.40
Mean :2008 Mean :82.56
3rd Qu.:2011 3rd Qu.:83.63
Max. :2014 Max. :89.06
le_agg_q2_F le_agg_q3_F le_agg_q4_F le_agg_q1_M
Min. :78.37 Min. :81.10 Min. :78.43 Min. :71.79
1st Qu.:84.18 1st Qu.:85.38 1st Qu.:87.02 1st Qu.:75.03
Median :85.19 Median :86.67 Median :88.24 Median :76.27
Mean :85.15 Mean :86.57 Mean :88.09 Mean :76.49
3rd Qu.:86.09 3rd Qu.:87.73 3rd Qu.:89.43 3rd Qu.:77.88
Max. :90.74 Max. :91.55 Max. :92.64 Max. :82.12
le_agg_q2_M le_agg_q3_M le_agg_q4_M
Min. :75.19 Min. :74.80 Min. :75.69
1st Qu.:79.81 1st Qu.:82.62 1st Qu.:84.78
Median :80.81 Median :83.71 Median :86.08
Mean :80.96 Mean :83.75 Mean :85.90
3rd Qu.:82.04 3rd Qu.:84.88 3rd Qu.:87.10
Max. :87.60 Max. :90.21 Max. :91.14
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.
mean(le$le_agg_q1_F)
[1] 82.55999
median(le$le_agg_q1_F)
[1] 82.40423
sd(le$le_agg_q1_F)
[1] 1.678479
range(le$le_agg_q1_F)
[1] 75.44415 89.06374
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.
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.
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.
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.
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.
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.
t.test(le12$le_agg_q1_M, le01$le_agg_q1_M)
Welch Two Sample t-test
data: le12$le_agg_q1_M and le01$le_agg_q1_M
t = 1.9303, df = 99.479, p-value = 0.05642
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.0200778 1.4601107
sample estimates:
mean of x mean of y
76.60453 75.88451
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.
Paired t-test
data: le12$le_agg_q1_M and le01$le_agg_q1_M
t = 4.2784, df = 50, p-value = 8.485e-05
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
0.3819977 1.0580353
sample estimates:
mean difference
0.7200165
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.
le14 <- le |>filter(year ==2014) |>mutate(gap_F = le_agg_q4_F - le_agg_q1_F)t.test(le14$gap_F)
One Sample t-test
data: le14$gap_F
t = 22.244, df = 50, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
5.787118 6.935997
sample estimates:
mean of x
6.361558
Now run a paired t-test that compares the richest women to the poorest women. Compare its output to the one-sample t-test.
Paired t-test
data: le14$le_agg_q4_F and le14$le_agg_q1_F
t = 22.244, df = 50, p-value < 2.2e-16
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
5.787118 6.935997
sample estimates:
mean difference
6.361558
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.
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.)
t.test(le_agg_q1_M ~ stateabbrv, data = ny_tx)
Welch Two Sample t-test
data: le_agg_q1_M by stateabbrv
t = 4.6616, df = 17.966, p-value = 0.000195
alternative hypothesis: true difference in means between group NY and group TX is not equal to 0
95 percent confidence interval:
0.688465 1.818472
sample estimates:
mean in group NY mean in group TX
78.29738 77.04392
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.
cor(le14$le_agg_q1_F, le14$le_agg_q1_M)
[1] 0.6537527
cor.test(le14$le_agg_q1_F, le14$le_agg_q1_M)
Pearson's product-moment correlation
data: le14$le_agg_q1_F and le14$le_agg_q1_M
t = 6.0476, df = 49, p-value = 1.977e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.4612737 0.7874640
sample estimates:
cor
0.6537527
It often helps to visualize a correlation. Make a scatterplot with a line of best fit.
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.
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().
summary(lm(le_agg_q1_M ~ stateabbrv, data = ny_tx))
Call:
lm(formula = le_agg_q1_M ~ stateabbrv, data = ny_tx)
Residuals:
Min 1Q Median 3Q Max
-1.2092 -0.3704 -0.1302 0.4770 1.3635
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 78.2974 0.1901 411.794 < 2e-16 ***
stateabbrvTX -1.2535 0.2689 -4.662 8.21e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7114 on 26 degrees of freedom
Multiple R-squared: 0.4553, Adjusted R-squared: 0.4343
F-statistic: 21.73 on 1 and 26 DF, p-value: 8.213e-05
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.
summary(lm(le_agg_q1_F ~ year, data = le))
Call:
lm(formula = le_agg_q1_F ~ year, data = le)
Residuals:
Min 1Q Median 3Q Max
-7.5132 -1.0521 -0.1263 1.0316 6.2509
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -62.48349 30.82856 -2.027 0.0431 *
year 0.07225 0.01536 4.705 3.05e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.654 on 712 degrees of freedom
Multiple R-squared: 0.03015, Adjusted R-squared: 0.02879
F-statistic: 22.14 on 1 and 712 DF, p-value: 3.053e-06
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.
summary(lm(le_agg_q1_M ~ year * stateabbrv, data = ny_tx))
Call:
lm(formula = le_agg_q1_M ~ year * stateabbrv, data = ny_tx)
Residuals:
Min 1Q Median 3Q Max
-0.6721 -0.2452 -0.0299 0.2376 0.6368
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -338.46566 49.75840 -6.802 4.91e-07 ***
year 0.20760 0.02479 8.376 1.39e-08 ***
stateabbrvTX 415.66925 70.36900 5.907 4.28e-06 ***
year:stateabbrvTX -0.20768 0.03505 -5.925 4.10e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3739 on 24 degrees of freedom
Multiple R-squared: 0.8611, Adjusted R-squared: 0.8438
F-statistic: 49.61 on 3 and 24 DF, p-value: 1.933e-10
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:
ggplot(data = ny_tx, aes(x = year, y = le_agg_q1_M, color = stateabbrv)) +geom_point() +geom_smooth(method ='lm')
`geom_smooth()` using formula = 'y ~ x'
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:
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/
