1  Data, variables, and distributions

Before we can say anything about a population, we have to be able to say something about the data in front of us. This chapter is about that: how data are arranged, what kinds of variables there are, how to look at them, and how to summarize them.

gss <- readRDS(here::here("data", "gss2024.rds")) |>
  haven::zap_labels()

We work mostly with the single-year file. When we look at change over time we will need the other one.

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       5

Tidy format: columns contain variables, each row is an observation.

Nearly every tool in this book assumes it, and a large fraction of real analytic work is getting data into this shape.1

1 The term is Hadley Wickham’s, and the underlying idea is much older. The useful thing about the definition is that it is checkable: point at a row and ask what it is an observation of. If you cannot answer in one word, the data are not tidy.

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   862

That 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, one row per year. We will build exactly 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.

“Usually fine” is a judgment, though, not a fact. 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? Almost certainly not. Treating it as a number implicitly says yes. We do it anyway, constantly, and Chapter 13 shows how to check whether it cost us anything.

WarningR will not stop you

R happily computes 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 (a number with no meaning whatsoever). There is no such thing as 2.77 marital statuses. Knowing what your variables are is your job.

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.

For this section we’ll use country-level data on internet access from the World Development Indicators, which gives us genuinely continuous variables. The GSS, being mostly categorical, does not.

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] 180

What 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")
Figure 1.1: Internet access by country, 2021.
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")
Figure 1.2: Internet access by country, 2021.

Bin width is a real choice. Too wide and you flatten everything interesting; too narrow and you are looking at noise. Try a few.

A density plot does something similar with a smooth curve instead of bars. It is easier on the eye and slightly harder to interpret honestly, since the smoothing can invent structure that is not there.

Show code
plt(~ intpct,
    data = wdi,
    type = "density",
    xlab = "% households with internet",
    ylab = "Density")
Figure 1.3: Internet access, as a density.
Show code
ggplot(wdi, aes(x = intpct)) +
  geom_density() +
  labs(x = "% households with internet", y = "Density")
Figure 1.4: Internet access, as a 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 = "")
Figure 1.5: Internet access by country, sorted.
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())
Figure 1.6: Internet access by country, sorted.

It shows every observation, in order, with nothing smoothed or binned away, so you can see the shape of the whole distribution and pick out individual countries at either end.

And a bar graph, for a categorical variable:

Show code
plt(~ income, data = wdi, type = "bar", xlab = "")
Figure 1.7: Countries by World Bank income group.
Show code
ggplot(wdi, aes(x = income)) +
  geom_bar() +
  labs(x = "", y = "Count")
Figure 1.8: Countries by World Bank income group.

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")
Figure 1.9: Internet access by region.
Show code
ggplot(wdi, aes(x = intpct, y = region)) +
  geom_point(alpha = .4) +
  labs(x = "% households with internet", y = "")
Figure 1.10: Internet access by region.

Points at identical values land on top of each other, which hides how many there are. Jittering adds a little random noise so the density becomes visible. The noise is cosmetic. Never analyze jittered values.

Show code
plt(intpct ~ region, data = wdi, type = "jitter", alpha = .4,
    xlab = "", ylab = "% households with internet")
Figure 1.11: The same plot, jittered.
Show code
ggplot(wdi, aes(x = intpct, y = region)) +
  geom_jitter(height = .15, width = 0, alpha = .4) +
  labs(x = "% households with internet", y = "")
Figure 1.12: The same plot, jittered.

A scatter plot puts two continuous variables against each other. Back in the GSS, years of schooling against vocabulary score:

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")
Figure 1.13: Years of schooling and vocabulary score. Jittered, since both variables are whole numbers.
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")
Figure 1.14: Years of schooling and vocabulary score. Jittered, since both variables are whole numbers.

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   315
Show code
plt(mean_wordsum ~ factor(degree), data = by_degree, type = "bar",
    xlab = "Degree (0 = less than HS, 4 = graduate)",
    ylab = "Mean words correct")
Figure 1.15: Mean vocabulary score by highest degree.
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")
Figure 1.16: Mean vocabulary score by highest degree.

This is the same comparison-of-conditional-summaries move. In Chapter 3 it is conditional probability; later it becomes a regression coefficient.

1.3.3 Time plots

Some questions need more than one year, and this is where the cumulative file is useful.

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  1447

There is the tidy version of the table we sketched earlier: one row per year, with year and p_support as honest columns.

Show code
plt(p_support ~ year, data = grass_trend, type = "b",
    xlab = "Year", ylab = "Proportion supporting legalization",
    ylim = c(0, 1))
Figure 1.17: Support for legalizing marijuana, GSS 1973-2024.
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")
Figure 1.18: Support for legalizing marijuana, GSS 1973-2024.

Support runs from about 19% in 1973 to about 69% in 2024. Whatever you think about the policy, that is a remarkable amount of movement in public opinion, and no single year of data could show it.

1.4 Descriptive statistics

We will distinguish between descriptive statistics for three different variable types:

  1. Continuous (interval, ratio, and some ordinal variables)
  2. Binary
  3. 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. We will not address that in this book, but you should know that dropping rows is a decision and not a neutral default.

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")
Figure 1.19: Distribution of wordsum.
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")
Figure 1.20: 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.339731

The 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.96374

1.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.227945

The 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 unremarkable, while a score of 0 or 10 is genuinely unusual.

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\).

Roman letters for what we compute, Greek for what we want to know. We almost never observe the Greek letters. The rest of this book is about what we can responsibly say about them anyway.

1.5 The normal distribution

The normal distribution is the one we will keep coming back to. When we resample and compute the mean, for example, our results will converge to that shape. That is the subject of 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)))
Figure 1.21: Normal probability density function.
Show code
ggplot() +
  xlim(-4, 4) +
  geom_function(fun = dnorm) +
  labs(title = "Normal probability density function",
       x = "SD diff. from mean",
       y = "" ~ phi(x) ~ "")
Figure 1.22: Normal probability density function.

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")
Figure 1.23: Two normal densities. The narrower one exceeds a density of 1.
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 = "")
Figure 1.24: Two normal densities. The narrower one exceeds a density of 1.

1.5.2 Cumulative density function

The density answers “how concentrated is probability here?” The cumulative density function answers the question we usually want: “what is the probability of a value this large or smaller?”

Show code
plt(pnorm(x) ~ x, type = "l", lwd = 2,
    xlab = "SD diff. from mean", ylab = expression(Phi(x)))
Figure 1.25: The normal cumulative density function.
Show code
ggplot() +
  xlim(-4, 4) +
  geom_function(fun = pnorm) +
  labs(x = "SD diff. from mean", y = "" ~ Phi(x) ~ "")
Figure 1.26: The normal cumulative density function.

Read it off at any point on the x-axis. At one standard deviation above the mean the curve has reached about 0.84:

pnorm(1)
[1] 0.8413447
Show 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")
Figure 1.27: Cumulative probability up to one standard deviation above the mean.
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)))
Figure 1.28: Cumulative probability up to one standard deviation above the mean.

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)
Figure 1.29: Distribution of wordsum with normal dist.
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()`).
Figure 1.30: Distribution of wordsum with normal dist.

Not bad. The curve is symmetric; the data are slightly lumpy. The curve extends past 10 and below 0, which is impossible for a ten-item test.

We can look at the same comparison cumulatively, which is often more revealing because nothing depends 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")
Figure 1.31: ECDF vs. Normal.
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")
Figure 1.32: ECDF vs. Normal.

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"))
Figure 1.33: Two datasets differing in one value. The means differ; the medians do not.
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"))
Figure 1.34: Two datasets differing in one value. The means differ; the medians do not.

The same thing happens in real data. 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. A handful of people reporting 24 hours pull the mean upward; they cannot pull the median, which only cares about position. The interquartile range (the distance from the 25th to the 75th percentile) is the robust counterpart to the SD.

Neither the mean nor the median is more correct; they answer different questions. The mean is what you want if the total matters. The median is what you want if the typical case matters, which is why income is nearly always reported as 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.3939539

So 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")
Figure 1.35: College degree, as a proportion.
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")
Figure 1.36: College degree, as a proportion.

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)
Figure 1.37: Position after 100 coin flips, 1000 people.
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")
Figure 1.38: Position after 100 coin flips, 1000 people.

The bell comes from the adding up. That is why the normal distribution turns up in places with nothing obviously normal about them, and Chapter 3 makes the point properly.

There is a nice physical simulation of this worth watching.

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 but is a judgment R will not make for you
  • look at a variable before summarizing it: histograms, densities, sorted dotplots, and strip plots each reveal something the others hide
  • 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