gss <- readRDS(here::here("data", "gss2024.rds")) |>
haven::zap_labels()1 Data, variables, and distributions
Before we can say anything about a population, we need to be able to describe the data we have. This chapter is about how data are arranged, what kinds of variables there are, and how to look at and summarize them.
We’ll mostly use the single-year file in this chapter. We’ll need the cumulative file when we look at change over time.
1.1 Data structure
Data, for our purposes, arrive as a rectangle.
gss |>
select(age, educ, degree, marital, wordsum) |>
head(10)# A tibble: 10 x 5
age educ degree marital wordsum
<dbl> <dbl> <dbl> <dbl> <dbl>
1 33 16 3 5 5
2 64 16 4 5 8
3 69 14 2 1 NA
4 19 12 1 5 6
5 70 13 1 3 NA
6 53 14 2 1 7
7 48 13 1 1 6
8 30 14 1 3 NA
9 60 14 2 1 NA
10 25 12 1 5 5Tidy format: columns contain variables, each row is an observation.
Nearly every tool in this book assumes data are in this format. A lot of real analytic work is just getting data into this shape.1
1 The term “tidy data” comes from Hadley Wickham.
1.1.1 Untidy data
Untidy data usually means a table where the columns are not variables but values. Suppose we summarize support for legalizing marijuana by year and lay it out with one column per year:
attitude 1973 1976 2024
1 Support legalization 18.7% 28.7% 68.6%
2 Sample size 1,471 1,447 862That is readable, and for a report it may be exactly what you want. But 1973 is not a variable. It is a value of the variable year. To compute with these data we would want them the other way around, with one row per year. We’ll make that table later in the chapter.
1.2 Types of variables
| Ratio | dollars; points (e.g., basketball) |
| Interval | degrees Celsius |
| Ordinal | clothing sizes; Likert scales |
| Nominal | race; sex; country |
The first two types are continuous or numeric. The second two types are categorical. Ordinal variables are often treated as numeric and this is usually fine.
But it’s still a choice you are making. Consider degree in the GSS, coded 0 for less than high school through 4 for a graduate degree. The categories are clearly ordered, but is the distance from 0 to 1 the same as from 3 to 4? Probably not. Treating degree as a number assumes that it is. We do this all the time anyway, and Chapter 13 shows one way to check whether it makes a difference.
WarningR will not stop you
R will happily compute the mean of a nominal variable if it is stored as a number. The GSS codes marital as 1 through 5, so mean(gss$marital) returns 2.77. This number means nothing at all! R doesn’t know that, so you have to.
NoteThe origins of “statistics”
The word statistics comes from the fact that it was information about the state. We’ll focus on information like this for now rather than thinking about samples of individuals.
1.3 Visualization basics
Consider two types of plots:
- univariate plots
- bivariate plots
These are also types of distributions.
The GSS is mostly categorical, so for this section we’ll use country-level data on internet access from the World Development Indicators.
wdi <- readRDS(here::here("data", "WDI.rds")) |>
filter(region != "Aggregates") |>
select(country,
iso = iso3c,
intpct = IT.NET.USER.ZS,
income,
region) |>
drop_na()
nrow(wdi)[1] 180What kinds of variables are these?
1.3.1 Univariate plots
Here’s a histogram.
Show code
plt(~ intpct,
data = wdi,
type = type_hist(breaks = "Sturges"),
main = "Internet access by country, 2021",
sub = "World Development Indicators data",
xlab = "% households with internet")
Show code
ggplot(wdi, aes(x = intpct)) +
geom_histogram(bins = nclass.Sturges(wdi$intpct),
fill = tableau10[1], color = "white", alpha = 0.8) +
labs(title = "Internet access by country, 2021",
subtitle = "World Development Indicators data",
x = "% households with internet",
y = "Count")
The width of the bins makes a big difference. If they’re too wide you’ll smooth over the interesting stuff, and if they’re too narrow you’ll mostly see noise. It’s a good idea to try a few.
A density plot does something similar with a smooth curve instead of bars. It’s easier on the eyes, but the smoothing can show structure that isn’t really there.
Show code
plt(~ intpct,
data = wdi,
type = "density",
xlab = "% households with internet",
ylab = "Density")
Show code
ggplot(wdi, aes(x = intpct)) +
geom_density() +
labs(x = "% households with internet", y = "Density")
Here’s a dotplot with the countries sorted by rank.
Show code
plt(~ sort(intpct),
data = wdi,
main = "Internet access by country, 2021",
sub = "World Development Indicators data",
ylab = "% households with internet",
xaxt = "n",
xlab = "")
Show code
ggplot(wdi |> arrange(intpct) |> mutate(rank = row_number()),
aes(x = rank, y = intpct)) +
geom_point(color = tableau10[1]) +
labs(title = "Internet access by country, 2021",
subtitle = "World Development Indicators data",
y = "% households with internet",
x = NULL) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank())
This shows every single observation, so nothing is smoothed over or hidden in a bin.
And a bar graph, for a categorical variable:
Show code
plt(~ income, data = wdi, type = "bar", xlab = "")
Show code
ggplot(wdi, aes(x = income)) +
geom_bar() +
labs(x = "", y = "Count")
1.3.2 Bivariate plots
A strip plot puts a continuous variable against a categorical one and shows every case.
Show code
plt(intpct ~ region, data = wdi, type = "p", alpha = .4,
xlab = "", ylab = "% households with internet")
Show code
ggplot(wdi, aes(x = intpct, y = region)) +
geom_point(alpha = .4) +
labs(x = "% households with internet", y = "")
Points at identical values land on top of each other, which hides how many there are. Jittering adds a little random noise so you can see how many there are. The noise is just for display (don’t ever analyze jittered values!).
Show code
plt(intpct ~ region, data = wdi, type = "jitter", alpha = .4,
xlab = "", ylab = "% households with internet")
Show code
ggplot(wdi, aes(x = intpct, y = region)) +
geom_jitter(height = .15, width = 0, alpha = .4) +
labs(x = "% households with internet", y = "")
A scatter plot shows two continuous variables. This one goes back to the GSS and plots vocabulary score against years of schooling:
Show code
ed_word <- gss |>
filter(!is.na(educ), !is.na(wordsum)) |>
select(educ, wordsum)
plt(wordsum ~ educ, data = ed_word, type = "jitter", alpha = .2,
xlab = "Years of schooling", ylab = "Words correct")
Show code
ggplot(ed_word, aes(x = educ, y = wordsum)) +
geom_jitter(alpha = 0.2, width = 0.3, height = 0.3) +
labs(x = "Years of schooling", y = "Words correct")
A bivariate bar graph compares a summary across groups.
by_degree <- gss |>
filter(!is.na(degree), !is.na(wordsum)) |>
group_by(degree) |>
summarize(mean_wordsum = mean(wordsum), n = n())
by_degree# A tibble: 5 x 3
degree mean_wordsum n
<dbl> <dbl> <int>
1 0 4.48 187
2 1 5.91 992
3 2 6.08 185
4 3 7.21 477
5 4 7.63 315Show code
plt(mean_wordsum ~ factor(degree), data = by_degree, type = "bar",
xlab = "Degree (0 = less than HS, 4 = graduate)",
ylab = "Mean words correct")
Show code
ggplot(by_degree, aes(x = factor(degree), y = mean_wordsum)) +
geom_col() +
labs(x = "Degree (0 = less than HS, 4 = graduate)",
y = "Mean words correct")
We will keep coming back to this idea of comparing a summary across groups. In Chapter 3 it shows up as conditional probability, and later it shows up as a regression coefficient.
1.3.3 Time plots
A time plot needs more than one year of data, so here we’ll use the cumulative file.
gss_all <- readRDS(here::here("data", "gss-1972-2024.rds"))
grass_trend <- gss_all |>
select(year, grass) |>
filter(!is.na(grass)) |>
mutate(year = as.numeric(year), support = as.numeric(grass) == 1) |>
group_by(year) |>
summarize(p_support = mean(support), n = n())
head(grass_trend, 3)# A tibble: 3 x 3
year p_support n
<dbl> <dbl> <int>
1 1973 0.187 1471
2 1975 0.213 1414
3 1976 0.287 1447This is the tidy version of the table from earlier, with one row per year.
Show code
plt(p_support ~ year, data = grass_trend, type = "b",
xlab = "Year", ylab = "Proportion supporting legalization",
ylim = c(0, 1))
Show code
ggplot(grass_trend, aes(x = year, y = p_support)) +
geom_line() +
geom_point() +
ylim(0, 1) +
labs(x = "Year", y = "Proportion supporting legalization")
Support runs from about 19% in 1973 to about 69% in 2024. That’s a huge change in public opinion!
1.4 Descriptive statistics
We will distinguish between descriptive statistics for three different variable types:
- Continuous (interval, ratio, and some ordinal variables)
- Binary
- Multinomial or categorical (nominal and some ordinal)
Let’s get a few variables to work with.
d <- gss |>
select(wordsum, # continuous
age, # continuous
educ, # continuous (make binary/ordinal)
marital) |> # nominal
drop_na() |>
mutate(marital_chr = case_when(marital == 1 ~ "married",
marital == 2 ~ "widowed",
marital %in% c(3, 4) ~ "sep. or div.",
marital == 5 ~ "never mar."))
nrow(d)[1] 2084
Note
Deleting cases with any missing data is sometimes OK, but there are often better ways to handle it. That’s a topic for a later course!
1.4.1 Continuous: wordsum
How many of the following words can you correctly define (picking the closest synonym via multiple choice)?
- Adept
- Audible
- Consume
- Coherent
- Emulate
- Erroneous
- Fortitude
- Misnomer
- Reverent
- Stimulus
I’m not 100% sure these are the words. But ChatGPT was pretty confident about it!
Show code
plt(~ wordsum, data = d, type = "hist", breaks = -0.5:10.5,
xlab = "Words Correct", ylab = "Count",
main = "Distribution of wordsum",
sub = "Source: 2024 General Social Survey")
Show code
ggplot(d, aes(x = wordsum, y = after_stat(count * 100 / nrow(d)))) +
geom_histogram(binwidth = 1, color = "white") +
scale_x_continuous(breaks = 0:10) +
labs(x = "Words Correct", y = "% of sample",
caption = "Source: 2024 General Social Survey",
title = "Distribution of wordsum")
1.4.2 Center and spread
We can use numbers to summarize a variable from a sample rather than having to reproduce the entire column of data every time.
1.4.2.1 Center
- mean
- median
- mode
1.4.2.2 Spread
- variance
- standard deviation
- interquartile range
We will focus on the mean, variance, and standard deviation first.
1.4.3 Mean and notation
\(\bar{x}\) is pronounced “x-bar” and is the mean of the variable \(x\) in a particular sample. We often use \(x\) when we are talking about a variable.
\[ \bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i \]
\(\Sigma\) means to sum; \(i\) is an index for each observation; \(n\) is the number of observations in the sample. So we are summing the values of \(x\) for each observation from the first \((i=1)\) to the last \((i=n)\) and then dividing by \(n\).
mean(d$wordsum)[1] 6.339731The mean is 6.34.
1.4.4 Variance
The sample variance tells you how spread out the data points are.
\[ s^2 = \frac{1}{n-1} \sum_{i=1}^{n}(x_i-\bar{x})^2 \]
This is sort of the average squared deviation from the mean. We divide by \(n-1\) for reasons you don’t need to worry about right now. We use squared deviations instead of absolute deviations for many reasons we are also not going to talk about right now!
var(d$wordsum)[1] 4.963741.4.5 Standard deviation
The variance \((s^2)\) has many desirable properties we’re not ready to discuss. Its main disadvantage is that it’s in squared units of the variable. By taking the square root, we get an interpretable value.
\[ s = \sqrt{\frac{1}{n-1} \sum_{i=1}^{n}(x_i-\bar{x})^2} \]
The standard deviation, \(s\), is a “typical deviation” from the mean.
sd(d$wordsum)[1] 2.227945The mean of wordsum is 6.34. The standard deviation is 2.23. We’ll talk more about how to use these values soon. For now, just remember that a deviation from the mean of that size or less would not be unusual. So anything between 4.1 and 8.6 is pretty typical.
1.4.6 Sample and population
So far, we’ve defined and discussed these as sample statistics rather than population parameters. The notation is slightly different for populations (although researchers are not always consistent).
- The sample mean is \(\bar{x}\); the population mean is \(\mu\).
- The sample variance is \(s^2\); the population variance is \(\sigma^2\).
- The sample SD is \(s\); the population SD is \(\sigma\).
In general, we use Roman letters for things we compute from a sample and Greek letters for the population values we want to know about. We almost never get to observe the Greek ones. A lot of this book is about what we can say about them anyway.
1.5 The normal distribution
We are going to see a lot of the normal distribution. When we resample and compute the mean, for example, our results will converge to that shape. We’ll see this in Chapter 3.
\[ f(x) = \frac{1}{\sigma \sqrt{2\pi}} e^{-\frac{(x - \mu)^2}{2\sigma^2}} \]
Warning
This is a probability density function. Don’t freak out about this. The important thing is to see \(\mu\) (the mean) and \(\sigma\) (the standard deviation). This just means that the probability of seeing a particular observation is a function of the mean and SD of the distribution.
Show code
x <- seq(-4, 4, length.out = 400)
plt(dnorm(x) ~ x, type = "l", lwd = 2,
main = "Normal probability density function",
xlab = "SD diff. from mean", ylab = expression(phi(x)))
Show code
ggplot() +
xlim(-4, 4) +
geom_function(fun = dnorm) +
labs(title = "Normal probability density function",
x = "SD diff. from mean",
y = "" ~ phi(x) ~ "")
1.5.1 What is “probability density”?
For a truly continuous variable, the probability that a variable takes on an exact value (say a height of 170.0000… cm) is zero.
This is quite different than, say, the probability that a fair coin comes up heads (.5) or that a person answers “yes” to a question about abortion in a population.
Tip
You could ask the probability that a person’s height is, say, greater than or equal to 169.5 and less than 170.5. As the width of this “window” shrinks to zero, the probability also shrinks to zero. But we can talk about the density of the probability in that area.
One consequence: density can be higher than 1!
Show code
plt(dnorm(x) ~ x, type = "l", lwd = 2, ylim = c(0, 1.4),
main = "Normal probability density function",
xlab = "SD diff. from mean", ylab = expression(phi(x)))
lines(x, dnorm(x, sd = 0.3), lty = 2, lwd = 2)
legend("topright", legend = c(expression(sigma == 1), expression(sigma == 0.3)),
lty = 1:2, lwd = 2, bty = "n")
Show code
ggplot() +
xlim(-4, 4) +
geom_function(fun = dnorm, aes(linetype = "sigma = 1")) +
geom_function(fun = dnorm, args = list(mean = 0, sd = .3),
aes(linetype = "sigma = 0.3")) +
labs(title = "Normal probability density function",
x = "SD diff. from mean",
y = "" ~ phi(x) ~ "", linetype = "")
1.5.2 Cumulative density function
The cumulative density function tells us the probability of seeing a value that large or smaller.
Show code
plt(pnorm(x) ~ x, type = "l", lwd = 2,
xlab = "SD diff. from mean", ylab = expression(Phi(x)))
Show code
ggplot() +
xlim(-4, 4) +
geom_function(fun = pnorm) +
labs(x = "SD diff. from mean", y = "" ~ Phi(x) ~ "")
For example, the probability of a value one standard deviation above the mean or lower is about 0.84:
pnorm(1)[1] 0.8413447Show code
plt(pnorm(x) ~ x, type = "l", lwd = 2,
xlab = "SD diff. from mean", ylab = expression(Phi(x)))
xs <- x[x <= 1]
polygon(c(xs, rev(xs)), c(pnorm(xs), rep(0, length(xs))),
col = adjustcolor(tableau10[1], alpha.f = 0.5), border = NA)
text(-2, 0.5, "Pr(x <= 1) = .84")
Show code
ggplot(data.frame(x = c(-4, 4)), aes(x = x)) +
stat_function(fun = pnorm) +
stat_function(fun = pnorm, geom = "area", xlim = c(-4, 1),
fill = tableau10[1], alpha = 0.6) +
annotate("text", x = -2, y = 0.5, label = "Pr(x <= 1) = .84", size = 4) +
labs(x = "SD diff. from mean", y = expression(Phi(x)))
1.5.3 Normal distribution: wordsum
Based on what we have already computed, we can approximate the distribution of wordsum using a normal distribution with a mean of 6.34 and a SD of 2.23.
We can write this as
\[ \text{wordsum} \sim \mathcal{N}(6.34, 2.23) \]
The first number is the mean and the second is the standard deviation.
How good is this approximation?
Show code
plt(~ wordsum, data = d, type = "hist", breaks = -0.5:10.5, freq = FALSE,
xlab = "Words Correct", ylab = "Density",
main = "Distribution of wordsum with normal dist.")
curve(dnorm(x, mean(d$wordsum), sd(d$wordsum)), add = TRUE, lwd = 2)
Show code
ggplot(d, aes(x = wordsum)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 1, color = "white") +
stat_function(fun = dnorm,
args = list(mean = mean(d$wordsum), sd = sd(d$wordsum)),
linewidth = 1.1) +
scale_x_continuous(breaks = -1:13, limits = c(-1, 13)) +
labs(x = "Words Correct", y = "Density",
caption = "Source: 2024 General Social Survey",
title = "Distribution of wordsum with normal dist.")Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_bar()`).
Not bad! But the data are a little lumpier than the curve, and the curve goes below 0 and past 10, which is impossible on a ten-item test.
We can also make the same comparison cumulatively. This is often more useful because it doesn’t depend on the bin width.
ecdf_data <- d |>
group_by(wordsum) |>
summarize(p = n() / nrow(d)) |>
mutate(cp = cumsum(p))Show code
plt(cp ~ I(wordsum + .5), data = ecdf_data, type = "s", lwd = 2,
xlim = c(-1, 13), ylim = c(0, 1),
xlab = "Words Correct", ylab = "Cumulative Probability")
curve(pnorm(x, mean(d$wordsum), sd(d$wordsum)), add = TRUE, lty = 2, lwd = 2)
legend("topleft", legend = c("observed", "normal"), lty = 1:2, lwd = 2, bty = "n")
Show code
ggplot(ecdf_data) +
stat_ecdf(aes(x = wordsum + .5, # shift to center of implied bar
y = cp),
geom = "step") +
stat_function(fun = pnorm,
args = list(mean = mean(d$wordsum), sd = sd(d$wordsum)),
linetype = "dashed") +
scale_x_continuous(breaks = -1:13, limits = c(-1, 13)) +
labs(x = "Words Correct", y = "Cumulative Probability")
1.6 Robust statistics
In inferential statistics (making inferences from samples to populations), we focus on the mean and standard deviation.
The median is used more as a descriptive statistic. It is called a robust statistic because it is insensitive to outliers. For example, the median age in the 2024 GSS is 49. This would be true even if we took the oldest person and made them 900 years old!
Data prep
median_data <- tibble(x1 = 1:11,
x2 = c(1:10, 20)) |>
pivot_longer(everything())Show code
plt(name ~ value | name, data = median_data, type = "p", cex = 1.5,
xlab = "", ylab = "", legend = FALSE)
text(c(5, 5, 10, 10), c(1.25, 2.25, 1.25, 2.25),
labels = c("mean = 6", "mean = 6.82", "median = 6", "median = 6"))
Show code
ggplot(median_data, aes(x = value, y = name, color = name)) +
geom_point() +
theme(legend.position = "none") +
labs(y = "", x = "") +
annotate("text",
x = c(5, 5, 10, 10),
y = c(1.25, 2.25, 1.25, 2.25),
label = c("mean = 6", "mean = 6.82", "median = 6", "median = 6"))
Here’s a real example using television hours in the GSS:
tv <- gss$tvhours[!is.na(gss$tvhours)]
c(mean = mean(tv), median = median(tv), max = max(tv), IQR = IQR(tv)) mean median max IQR
3.302509 2.000000 24.000000 3.000000 Half of respondents report 2 hours a day or fewer, but the mean is 3.3. The few people who say they watch 24 hours a day pull the mean up, but they can’t move the median. The interquartile range (the distance from the 25th to the 75th percentile) is the robust version of the SD.
Neither the mean nor the median is more correct. They just tell you different things. This is why variables with a long tail, like income, are usually summarized with a median.
1.7 Bernoulli distribution
You’ve seen this before but some statistical distributions have only two options. If we want to describe the proportion of US adults who have a college degree, we can describe this as a Bernoulli distribution.
d <- d |>
mutate(college = if_else(educ >= 16, TRUE, FALSE))
mean(d$college)[1] 0.3939539So college is Bernoulli with \(p = 0.394\).
Show code
plt(~ college, data = d, type = "bar",
xlab = "College degree?", ylab = "Count",
sub = "Source: 2024 General Social Survey")
Show code
ggplot(d, aes(x = college, y = after_stat(count / nrow(d)))) +
geom_bar() +
labs(x = "College degree?", y = "Proportion",
caption = "Source: 2024 General Social Survey")
1.7.1 One- and two-parameter distributions
The normal distribution has two parameters, \(\mu\) and \(\sigma\). This is because the normal distribution is defined by the location of its center and the width of its spread.
The Bernoulli distribution has only one parameter, which is \(p\) (sometimes people use \(\pi\)). This is just the probability of a “yes,” or, as it is often called, a “success.”
But this doesn’t mean that the Bernoulli doesn’t have center and spread.
1.7.2 Spread of the Bernoulli distribution
Variance is a measure of uncertainty about where the data are. Imagine two alternatives: a Bernoulli distribution with \(p = .01\) and one with \(p = .50\). There’s a lot more uncertainty about the latter!
So the spread is also a function of \(p\). In other words, \(p\) determines both center and spread.
For a variable, \(X\), \(\text{Var}[X] = p(1-p)\). Therefore it’s also true that \(\text{SD}[X] = \sqrt{p(1-p)}\).
p <- mean(d$college)
c(variance = p * (1 - p), sd = sqrt(p * (1 - p))) variance sd
0.2387542 0.4886248 1.7.3 From Bernoulli to normal
The normal distribution can be derived as the sum of many Bernoulli trials. For example, imagine we start with 100 people standing on the halfway line of a football field. Each person flips a coin and, if it’s heads, takes a step forward (say one meter). If tails, they take a step backward (one meter). What would things look like after 100 trials?
set.seed(522)
take_walk <- function(steps = 100) {
sum(sample(c(-1, 1), steps, replace = TRUE))
}
walks <- tibble(person = 1:1000) |>
rowwise() |>
mutate(position = take_walk())
c(mean = mean(walks$position), sd = sd(walks$position)) mean sd
0.56200 10.17171 Show code
plt(~ position, data = walks, type = "hist", freq = FALSE,
xlab = "Meters from the halfway line", ylab = "Density")
curve(dnorm(x, mean(walks$position), sd(walks$position)), add = TRUE, lwd = 2)
Show code
ggplot(walks, aes(x = position)) +
geom_histogram(aes(y = after_stat(density)), binwidth = 2, color = "white") +
stat_function(fun = dnorm,
args = list(mean = mean(walks$position), sd = sd(walks$position)),
linewidth = 1) +
labs(x = "Meters from the halfway line", y = "Density")
The normal shape comes from adding up lots of small random steps. This is why the normal distribution shows up in so many places. We’ll come back to this in Chapter 3.
Here’s a nice physical simulation of this.
1.8 Recap
- tidy format: columns contain variables, each row is an observation
- variables are ratio, interval, ordinal, or nominal; the first two are continuous, the last two categorical
- ordinal variables are often treated as numeric, which is usually fine
- histograms, density plots, dotplots, and bar graphs show one variable; strip plots, scatter plots, and time plots show two
- the mean and median answer different questions; the median is robust to outliers and the mean is not
- variance is sort of the average squared deviation from the mean, and the standard deviation is its square root
- sample quantities get Roman letters (\(\bar{x}\), \(s\)), population quantities get Greek (\(\mu\), \(\sigma\))
- the normal distribution has two parameters, \(\mu\) and \(\sigma\), and its \(y\) axis is density rather than probability, which is why it can exceed 1
- the Bernoulli has only \(p\), which fixes both its center and its spread: \(\text{Var}[X] = p(1-p)\)
- the normal distribution can be derived as the sum of many Bernoulli trials