gss <- readRDS(here::here("data", "gss2024.rds")) |>
haven::zap_labels()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.
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.
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.
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")
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")
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")
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())
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 = "")
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 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")
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 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")
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 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")
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")
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))
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. 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:
- 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
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")
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.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}} \]
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.
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 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)))
Show code
ggplot() +
xlim(-4, 4) +
geom_function(fun = pnorm) +
labs(x = "SD diff. from mean", y = "" ~ Phi(x) ~ "")
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")
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. 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")
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"))
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")
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 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