Visualizing Uncertainty

Part I: Random variables

In data science, we often deal with data that is affected by chance in some way: the data comes from a random sample, the data is affected by measurement error, or the data measures some outcome that is random in nature. Being able to quantify the uncertainty introduced by randomness is one of the most important jobs of a data analyst. Statistical inference offers a framework, as well as several practical tools, for doing this. The first step is to learn how to mathematically describe random variables.

In this section, we introduce random variables and their properties starting with their application to games of chance. We then describe some of the events surrounding the financial crisis of 2007-2008¹ using probability theory. This financial crisis was in part caused by underestimating the risk of certain securities² sold by financial institutions. Specifically, the risks of mortgage-backed securities (MBS) and collateralized debt obligations (CDO) were grossly underestimated. These assets were sold at prices that assumed most homeowners would make their monthly payments, and the probability of this not occurring was calculated as being low. A combination of factors resulted in many more defaults than were expected, which led to a price crash of these securities. As a consequence, banks lost so much money that they needed government bailouts to avoid closing down completely.

Definition of Random variables

Random variables are numeric outcomes resulting from random processes. We can easily generate random variables using some of the simple examples we have shown. For example, define X to be 1 if a bead is blue and red otherwise:

beads <- rep( c("red", "blue"), times = c(2,3))
X <- ifelse(sample(beads, 1) == "blue", 1, 0)

Here X is a random variable: every time we select a new bead the outcome changes randomly. See below:

ifelse(sample(beads, 1) == "blue", 1, 0)

[1] 1

ifelse(sample(beads, 1) == "blue", 1, 0)

[1] 0

ifelse(sample(beads, 1) == "blue", 1, 0)

[1] 0

Sometimes it’s 1 and sometimes it’s 0.

Sampling models

Many data generation procedures, those that produce the data we study, can be modeled quite well as draws from an urn. For instance, we can model the process of polling likely voters as drawing 0s (Republicans) and 1s (Democrats) from an urn containing the 0 and 1 code for all likely voters. In epidemiological studies, we often assume that the subjects in our study are a random sample from the population of interest. The data related to a specific outcome can be modeled as a random sample from an urn containing the outcome for the entire population of interest. Similarly, in experimental research, we often assume that the individual organisms we are studying, for example worms, flies, or mice, are a random sample from a larger population. Randomized experiments can also be modeled by draws from an urn given the way individuals are assigned into groups: when getting assigned, you draw your group at random. Sampling models are therefore ubiquitous in data science. Casino games offer a plethora of examples of real-world situations in which sampling models are used to answer specific questions. We will therefore start with such examples.

Suppose a very small casino hires you to consult on whether they should set up roulette wheels. To keep the example simple, we will assume that 1,000 people will play and that the only game you can play on the roulette wheel is to bet on red or black. The casino wants you to predict how much money they will make or lose. They want a range of values and, in particular, they want to know what’s the chance of losing money. If this probability is too high, they will pass on installing roulette wheels.

We are going to define a random variable $S$ that will represent the casino’s total winnings. Let’s start by constructing the urn. A roulette wheel has 18 red pockets, 18 black pockets and 2 green ones. So playing a color in one game of roulette is equivalent to drawing from this urn:

color <- rep(c("Black", "Red", "Green"), c(18, 18, 2))

The 1,000 outcomes from 1,000 people playing are independent draws from this urn. If red comes up, the gambler wins and the casino loses a dollar, so we draw a -1. Otherwise, the casino wins a dollar and we draw a 1. To construct our random variable $S$, we can use this code:

n <- 1000
X <- sample(ifelse(color == "Red", -1, 1),  n, replace = TRUE)
X[1:10]

 [1] -1  1  1 -1 -1 -1  1  1  1  1

Because we know the proportions of 1s and -1s, we can generate the draws with one line of code, without defining color:

X <- sample(c(-1,1), n, replace = TRUE, prob=c(9/19, 10/19))

We call this a sampling model since we are modeling the random behavior of roulette with the sampling of draws from an urn. The total winnings $S$ is simply the sum of these 1,000 independent draws:

set.seed(100)
X <- sample(c(-1,1), n, replace = TRUE, prob=c(9/19, 10/19))
S <- sum(X)
S

[1] 16

The probability distribution of a random variable

If you run the code above, you see that $S$ changes every time. This is, of course, because $S$ is a random variable. The probability distribution of a random variable tells us the probability of the observed value falling at any given interval. So, for example, if we want to know the probability that we lose money, we are asking the probability that $S$ is in the interval $S<0$.

Note that if we can define a cumulative distribution function $F(a)=\text{Pr}(S\le a)$, then we will be able to answer any question related to the probability of events defined by our random variable $S$, including the event $S<0$. We call this $F$ the random variable’s distribution function.

We can estimate the distribution function for the random variable $S$ by using a Monte Carlo simulation to generate many realizations of the random variable. With this code, we run the experiment of having 1,000 people play roulette, over and over, specifically $B=10,000$ times:

set.seed(200)
n <- 1000
B <- 10000
roulette_winnings <- function(n){
  X <- sample(c(-1,1), n, replace = TRUE, prob=c(9/19, 10/19))
  sum(X)
}
S <- replicate(B, roulette_winnings(n))

Now we can ask the following: in our simulations, how often did we get sums less than or equal to a?

mean(S <= a)

This will be a very good approximation of $F(a)$ and we can easily answer the casino’s question: how likely is it that we will lose money? We can see it is quite low:

mean(S<0)

[1] 0.0474

We can visualize the distribution of $S$ by creating a histogram showing the probability $F(b)-F(a)$ for several intervals $(a,b]$:

We see that the distribution appears to be approximately normal. A qq-plot will confirm that the normal approximation is close to a perfect approximation for this distribution. If, in fact, the distribution is normal, then all we need to define the distribution is the average and the standard deviation. Because we have the original values (9/19, 10/19) and $n=1000$ plays from which the distribution is created, we already know (or can calculate) these values. We can also easily compute the sample analogs with mean(S) and sd(S). The blue curve you see added to the histogram above is a normal density with the sample average and standard deviation. The gray dashed line is the expected value computed from the known distribution $n\times (\frac{10}{19}\times 1-\frac{9}{19}\times 1)$. They are very close, but are not exactly the same.

This average and this standard deviation have special names. They are referred to as the expected value and standard error of the random variable $S$. This refers to the underlying random variable $S$, while the sample mean and sample standard deviation are what we can always calculate. We will say more about these in the next section where we will discuss an incredibly useful approximation provided by mathematical theory that applies generally to sums and averages of draws from any urn: the Central Limit Theorem (CLT).

Distributions versus probability distributions

Before we continue, let’s make an important distinction and connection between the distribution of a list of numbers and a probability distribution. In the visualization lectures, we described how any list of numbers $x_1,\dots ,x_n$ has a distribution. The definition is quite straightforward. We define $F(a)$ as the function that tells us what proportion of the list is less than or equal to $a$. Because they are useful summaries when the distribution is approximately normal, we define the average and standard deviation. These are defined with a straightforward operation of the vector containing the list of numbers x:

m <- sum(x)/length(x)
s <- sqrt(sum((x - m)^2) / length(x))

A random variable $X$ has a distribution function. To define this, we do not need a list of numbers. It is a theoretical concept. In this case, we define the distribution as the $F(a)$ that answers the question: what is the probability that $X$ is less than or equal to $a$? There is no list of numbers.

However, if $X$ is defined by drawing from an urn with numbers or colors in it, then there is a list: the list of numbers or colors inside the urn. In this case, the distribution of that list (of numbers or colors inside the urn) is the probability distribution of $X$ and the average and standard deviation of that list are the expected value and standard error of the random variable.

Another way to think about it that does not involve an urn is to run a Monte Carlo simulation and generate a very large list of outcomes of $X$. These outcomes are a list of numbers. The distribution of this list will be a very good approximation of the probability distribution of $X$. The longer the list, the better the approximation. But it will always be an approximation - the actual contents of the underlying distribution function may not be knowable. Still, the average and standard deviation of this list will approximate the expected value and standard error of the random variable.

Often, we’d like to know about the probability distribution (the distribution function), but we only have a sample. This is the fundamental challenge of statistical inference: what can we learn about the true distribution function using only a sample? If we can’t learn the distribution function itself, can we at least reject a hypothetical statement about the distribution function?

Notation for random variables

In statistical textbooks, upper case letters are used to denote random variables and we follow this convention here. Lower case letters are used for observed values. You will see some notation that includes both. For example, you will see events defined as $X\le x$. Here $X$ is a random variable, making it a random event, and $x$ is an arbitrary value and not random. So, for example, $X$ might represent the number on a die roll and $x$ will represent an actual value we see 1, 2, 3, 4, 5, or 6. So in this case, the probability of $X=x$ is 1/6 regardless of the observed value $x$. This notation is a bit strange because, when we ask questions about probability, $X$ is not an observed quantity. Instead, it’s a random quantity that we will see in the future. We can talk about what we expect it to be, what values are probable, but not what it is. But once we have data, we do see a realization of $X$. So data scientists talk of what could have been after we see what actually happened.

The expected value and standard error

We have described sampling models for draws. We will now go over the mathematical theory that lets us approximate the probability distributions for the sum of draws. Once we do this, we will be able to help the casino predict how much money they will make. The same approach we use for the sum of draws will be useful for describing the distribution of averages and proportion which we will need to understand how polls work.

The first important concept to learn is the expected value. In statistics books, it is common to use letter $\text{E}$ like this:

$$\text{E}[X]$$

to denote the expected value of the random variable $X$.

A random variable will vary around its expected value in a way that if you take the average of many, many draws, the average of the draws will approximate the expected value, getting closer and closer the more draws you take.

Theoretical statistics provides techniques that facilitate the calculation of expected values in different circumstances. For example, a useful formula tells us that the expected value of a random variable defined by one draw is the average of the numbers in the urn (which is the probability distribution of the random variable). In the urn used to model betting on red in roulette, we have 20 one dollars and 18 negative one dollars. The expected value is thus:

$$\text{E}[X]=(20+-18)/38$$

which is about 5 cents. It is a bit counterintuitive to say that $X$ varies around 0.05, when the only values it takes is 1 and -1. One way to make sense of the expected value in this context is by realizing that if we play the game over and over, the casino wins, on average, 5 cents per game. A Monte Carlo simulation confirms this:

B <- 10^6
x <- sample(c(-1,1), B, replace = TRUE, prob=c(9/19, 10/19))
mean(x)

[1] 0.053422

Note that if we didn’t know the probability distribution but had the urn, we could create the sample $x$.

In general, if the urn has two possible outcomes, say $a$ and $b$, with proportions $p$ and $1-p$ respectively, and you know $p$ from the probability distribution, the average is:

$$\text{E}[X]=pa+(1-p)b$$

TRY IT

What is the expectation of a random variable defined as a roll of a (fair, standard) die?

To see this, notice that if there are $n$ beads in the urn, then we have $np$ $a$s and $n(1-p)$ $b$s and because the average is the sum, $n\times a\times p+n\times b\times (1-p)$, divided by the total $n$, we get that the average is $ap+b(1-p)$.

Now the reason we define the expected value is because this mathematical definition turns out to be useful for approximating the probability distributions of the sum, which then is useful for describing the distribution of averages and proportions. The first useful fact is that the expected value of the sum of the draws is $E[nX]=nE[X]$ when $n$ is a constant:

$$\text{number of draws }\times \text{ average of the numbers in the urn}$$

So if 1,000 people play roulette, the casino expects to win, on average, about 1,000 $\times$ $0.05 = $50. But this is an expected value. How different can one observation be from the expected value? The casino really needs to know this. What is the range of possibilities? If negative numbers are too likely, they will not install roulette wheels. Statistical theory once again answers this question. The standard error (SE) gives us an idea of the size of the variation around the expected value. In statistics books, it’s common to use:

$$\text{SE}[X]$$

to denote the standard error of a random variable. This, like expectation, is a theoretical, population concept based on the probability distribution (e.g. “what is in the urn”)

If our draws are independent, then the standard error of the sum is given by the equation:

$$\sqrt{\text{number of draws }}\times \text{ standard deviation of the numbers in the urn}$$

Using the definition of standard deviation, we can derive, with a bit of math, that if an urn contains two values $a$ and $b$ with proportions $p$ and $(1-p)$, respectively, the standard deviation is:

$$\mid b-a\mid \sqrt{p(1-p)}.$$

So in our roulette example, the standard deviation of the values inside the urn is: $\mid 1-(-1)\mid \sqrt{10/19\times 9/19}$ or:

2 * sqrt(90)/19

[1] 0.998614

The standard error tells us the typical difference between a random variable and its expectation. Since one draw is obviously the sum of just one draw, we can use the formula above to calculate that the random variable defined by one draw has an expected value of 0.05 and a standard error of about 1. This makes sense since we either get 1 or -1, with 1 slightly favored over -1.

Using the formula above, the sum of 1,000 people playing has standard error of about $32:

n <- 1000
sqrt(n) * 2 * sqrt(90)/19

[1] 31.57895

As a result, when 1,000 people bet on red, the casino is expected to win $50 with a standard error of $32. It therefore seems like a safe bet. But we still haven’t answered the question: how likely is it to lose money? Here the CLT will help.

Advanced note: Before continuing we should point out that exact probability calculations for the casino winnings can be performed with the binomial distribution. However, here we focus on the CLT, which can be generally applied to sums of random variables in a way that the binomial distribution can’t.

Population SD versus the sample SD

The standard deviation of a list x (below we use heights as an example) is defined as the square root of the average of the squared differences:

library(dslabs)
x <- heights$height
m <- mean(x)
s <- sqrt(mean((x-m)^2))

The SD is the the square root of the sample variance, and the sample variance is the square of the sample SD. Using mathematical notation we write:

$$\mu =\frac{1}{n}\sum _{i=1}^n{x}_i$$ and $$\sigma =\sqrt{\frac{1}{n}\sum _{i=1}^n(x_i-\mu {)}^2}$$

However, be aware that the sd function returns a slightly different result:

identical(s, sd(x))

[1] FALSE

s-sd(x)

[1] -0.001942661

This is because the sd function R does not return the sd of the list, but rather uses a formula that estimates standard deviations of a population from a random sample $X_1,\dots ,X_N$ which, for reasons not discussed here, divide the sum of squares by the $N-1$.

$$\overline{X}=\frac{1}{N}\sum _{i=1}^N{X}_i,s=\sqrt{\frac{1}{N-1}\sum _{i=1}^N(X_i-\overline{X}{)}^2}$$

You can see that this is the case by typing:

n <- length(x)
s-sd(x)*sqrt((n-1) / n)

[1] 0

For all the theory applications discussed here, you need to compute the actual standard deviation as defined:

sqrt(mean((x-m)^2))

So be careful when using the sd function in R. However, keep in mind that throughout these notes we sometimes use the sd function when we really want the actual SD. This is because when the list size is big, these two are practically equivalent since $\sqrt{(N-1)/N}\approx 1$. That sd() calculates sample SD not population SD by default will be useful in the near future.

Central Limit Theorem

We really want to know the probability that a casino will lose money. We know the probability distribution of roulette which let us calculate the expected value and the standard deviation of a random variable that consists of 1,000 people playing roulette, but we don’t know the distribution of the sum of those random variables.

The Central Limit Theorem (CLT) tells us that when the number of plays (or “draws”), also called the sample size, is large, the probability distribution of the sum of the independent draws is approximately normal. Because sampling models are used for so many data generation processes, the CLT is considered one of the most important mathematical insights in history.

Previously, we discussed that if we know that the distribution of a list of numbers is approximated by the normal distribution, all we need to describe the list are the average and standard deviation. We also know that the same applies to probability distributions. If a random variable has a probability distribution that is approximated with the normal distribution, then all we need to describe the probability distribution are the average and standard deviation, referred to as the expected value and standard error.

We previously ran this Monte Carlo simulation:

n <- 1000
B <- 10000
roulette_winnings <- function(n){
  X <- sample(c(-1,1), n, replace = TRUE, prob=c(9/19, 10/19))
  sum(X)
}
S <- replicate(B, roulette_winnings(n))

The Central Limit Theorem (CLT) tells us that the sum $S$ is approximated by a normal distribution. Amazing!

Using the formulas above, we know that the expected value and standard error are:

n * (20-18)/38

[1] 52.63158

sqrt(n) * 2 * sqrt(90)/19

[1] 31.57895

The theoretical values above match those obtained with the Monte Carlo simulation:

mean(S)

[1] 52.3038

sd(S)

[1] 31.20746

Using the CLT, we can skip the Monte Carlo simulation and instead compute the probability of the casino losing money using this approximation, which is only possible because the CLT tells us we can use the normal pnorm():

mu <- n * (20-18)/38
se <-  sqrt(n) * 2 * sqrt(90)/19
pnorm(0, mu, se)

[1] 0.04779035

which is also in very good agreement with our Monte Carlo result:

mean(S < 0)

[1] 0.0446

It’s pretty powerful to have a function representing the outcome’s distribution, instead of having a (possibly huge) list of sampled outcomes.

How large is large in the Central Limit Theorem?

The CLT works when the number of draws is large. But large is a relative term. In many circumstances as few as 30 draws is enough to make the CLT useful. In some specific instances, as few as 10 is enough. However, these should not be considered general rules. Note, for example, that when the probability of success is very small, we need much larger sample sizes.

By way of illustration, let’s consider the lottery. In the lottery, the chances of winning are less than 1 in a million. Thousands of people play so the number of draws is very large. Yet the number of winners, the sum of the draws, range between 0 and 4. This sum is certainly not well approximated by a normal distribution, so the CLT does not apply, even with the very large sample size. This is generally true when the probability of a success is very low. In these cases, the Poisson distribution is more appropriate.

You can examine the properties of the Poisson distribution using dpois and ppois. You can generate random variables following this distribution with rpois. However, we do not cover the theory here. You can learn about the Poisson distribution in any probability textbook and even Wikipedia³

Statistical properties of averages

There are several useful mathematical results that we used above and often employ when working with data. We list them below.

1. The expected value of the sum of random variables is the sum of each random variable’s expected value. We can write it like this:

$$\text{E}[X_1+X_2+\cdots +X_n]=\text{E}[X_1]+\text{E}[X_2]+\cdots +\text{E}[X_n]$$

If the $X$ are independent draws from the urn, then they all have the same expected value. Let’s call it $\mu$ and thus:

$$\text{E}[X_1+X_2+\cdots +X_n]=n\mu$$

which is another way of writing the result we show above for the sum of draws.

2. The expected value of a non-random constant times a random variable is the non-random constant times the expected value of a random variable. This is easier to explain with symbols:

$$\text{E}[aX]=a\times \text{E}[X]$$

To see why this is intuitive, consider change of units. If we change the units of a random variable, say from dollars to cents, the expectation should change in the same way. A consequence of the above two facts is that the expected value of the average of independent draws from the same urn is the expected value of the urn, call it $\mu$ again:

$$\text{E}[(X_1+X_2+\cdots +X_n)/n]=\text{E}[X_1+X_2+\cdots +X_n]/n=n\mu /n=\mu$$

3. The square of the standard error of the sum of independent random variables is the sum of the square of the standard error of each random variable. This one is easier to understand in math form:

$$\text{SE}[X_1+X_2+\cdots +X_n]=\sqrt{\text{SE}[X_1{]}^2+\text{SE}[X_2{]}^2+\cdots +\text{SE}[X_n{]}^2}$$

The square of the standard error is referred to as the variance in statistical textbooks. Note that this particular property is not as intuitive as the previous three and more in depth explanations can be found in statistics textbooks.

4. The standard error of a non-random constant times a random variable is the non-random constant times the random variable’s standard error. As with the expectation: $$\text{SE}[aX]=a\times \text{SE}[X]$$

To see why this is intuitive, again think of units.

A consequence of 3 and 4 is that the standard error of the average of independent draws from the same urn is the standard deviation of the urn divided by the square root of $n$ (the number of draws), call it $\sigma$:

$$\begin{array}{rl}\text{SE}[(X_1+X_2+\cdots +X_n)/n]& =\text{SE}[X_1+X_2+\cdots +X_n]/n\\ & =\sqrt{\text{SE}[X_1{]}^2+\text{SE}[X_2{]}^2+\cdots +\text{SE}[X_n{]}^2}/n\\ & =\sqrt{{\sigma }^2+{\sigma }^2+\cdots +{\sigma }^2}/n\\ & =\sqrt{n{\sigma }^2}/n\\ & =\sigma /\sqrt{n}\end{array}$$

5. If $X$ is a normally distributed random variable, then if $a$ and $b$ are non-random constants, $aX+b$ is also a normally distributed random variable. All we are doing is changing the units of the random variable by multiplying by $a$, then shifting the center by $b$.

Note that statistical textbooks use the Greek letters $\mu$ and $\sigma$ to denote the expected value and standard error, respectively. This is because $\mu$ is the Greek letter for $m$, the first letter of mean, which is another term used for expected value. Similarly, $\sigma$ is the Greek letter for $s$, the first letter of standard error.

Law of large numbers

An important implication of the final result is that the standard error of the average becomes smaller and smaller as $n$ grows larger. When $n$ is very large, then the standard error is practically 0 and the average of the draws converges to the average of the urn. This is known in statistical textbooks as the law of large numbers or the law of averages.

Misinterpreting law of averages

The law of averages is sometimes misinterpreted. For example, if you toss a coin 5 times and see a head each time, you might hear someone argue that the next toss is probably a tail because of the law of averages: on average we should see 50% heads and 50% tails. A similar argument would be to say that red “is due” on the roulette wheel after seeing black come up five times in a row. These events are independent so the chance of a coin landing heads is 50% regardless of the previous 5. This is also the case for the roulette outcome. The law of averages applies only when the number of draws is very large and not in small samples. After a million tosses, you will definitely see about 50% heads regardless of the outcome of the first five tosses.

Another funny misuse of the law of averages is in sports when TV sportscasters predict a player is about to succeed because they have failed a few times in a row.

Probabilistic thinking is central in the human experience. How we describe that thinking is mixed, but most of the time we use (rather imprecise) language. With only a few moments of searching, one can find thousands of articles that use probabilistic words to describe events. Here are some examples:

“‘Highly unlikely’ State of the Union will happen amid shutdown” – The Hill

“Tiger Woods makes Masters 15th and most improbable major” – Fox

“Trump predicts ‘very good chance’ of China trade deal” – CNN

Yet people don’t have a good sense of what these things mean. Uncertainty is key to data analytics: if we were certain of things, A study in the 1960s explored the perception of probabilistic words like these among NATO officers. A more modern replication of this found the following basic pattern:

A deep, basic fact about humans is that we struggle to understand probabilities. But visualizing things can help. The graphic above shows the uncertainty about uncertainty (very meta). We can convey all manner of uncertainty with clever graphics. Today’s practical example works through some of the myriad of ways to visualize variable data. We’ll cover some territory that we’ve already hit, in hopes of locking in some key concepts.

Part 2: Statistical Inference and Polls

In this Example we will describe, in some detail, how poll aggregators such as FiveThirtyEight use data to predict election outcomes. To understand how they do this, we first need to learn the basics of Statistical Inference, the part of statistics that helps distinguish patterns arising from signal from those arising from chance. Statistical inference is a broad topic and here we go over the very basics using polls as a motivating example. To describe the concepts, we complement the mathematical formulas with Monte Carlo simulations and R code.

Polls

Opinion polling has been conducted since the 19th century. The general goal is to describe the opinions held by a specific population on a given set of topics. In recent times, these polls have been pervasive during presidential elections. Polls are useful when interviewing every member of a particular population is logistically impossible. The general strategy is to interview a smaller group, chosen at random, and then infer the opinions of the entire population from the opinions of the smaller group. Statistical theory is used to justify the process. This theory is referred to as inference and it is the main topic of this chapter.

Perhaps the best known opinion polls are those conducted to determine which candidate is preferred by voters in a given election. Political strategists make extensive use of polls to decide, among other things, how to invest resources. For example, they may want to know in which geographical locations to focus their “get out the vote” efforts.

Elections are a particularly interesting case of opinion polls because the actual opinion of the entire population is revealed on election day. Of course, it costs millions of dollars to run an actual election which makes polling a cost effective strategy for those that want to forecast the results.

Although typically the results of these polls are kept private, similar polls are conducted by news organizations because results tend to be of interest to the general public and made public. We will eventually be looking at such data.

Real Clear Politics⁴ is an example of a news aggregator that organizes and publishes poll results. For example, they present the following poll results reporting estimates of the popular vote for the 2016 presidential election⁵:

Although in the United States the popular vote does not determine the result of the presidential election, we will use it as an illustrative and simple example of how well polls work. Forecasting the election is a more complex process since it involves combining results from 50 states and DC and we will go into some detail on this later.

Let’s make some observations about the table above. First, note that different polls, all taken days before the election, report a different spread: the estimated difference between support for the two candidates. Notice also that the reported spreads hover around what ended up being the actual result: Clinton won the popular vote by 2.1%. We also see a column titled MoE which stands for margin of error.

In this example, we will show how the probability concepts we learned in the previous content can be applied to develop the statistical approaches that make polls an effective tool. We will learn the statistical concepts necessary to define estimates and margins of errors, and show how we can use these to forecast final results relatively well and also provide an estimate of the precision of our forecast. Once we learn this, we will be able to understand two concepts that are ubiquitous in data science: confidence intervals and p-values. Finally, to understand probabilistic statements about the probability of a candidate winning, we will have to learn about Bayesian modeling. In the final sections, we put it all together to recreate the simplified version of the FiveThirtyEight model and apply it to the 2016 election.

We start by connecting probability theory to the task of using polls to learn about a population.

The sampling model for polls

To help us understand the connection between polls and what we have learned, let’s construct a similar situation to the one pollsters face. To mimic the challenge real pollsters face in terms of competing with other pollsters for media attention, we will use an urn full of beads to represent voters and pretend we are competing for a $25 dollar prize. The challenge is to guess the spread between the proportion of blue and red beads in this hypothetical urn.

Before making a prediction, you can take a sample (with replacement) from the urn. To mimic the fact that running polls is expensive, it costs you 10 cents per each bead you sample. Therefore, if your sample size is 250, and you win, you will break even since you will pay $25 to collect your $25 prize. Your entry into the competition can be an interval. If the interval you submit contains the true proportion, you get half what you paid and pass to the second phase of the competition. In the second phase, the entry with the smallest interval is selected as the winner.

The dslabs package includes a function that shows a random draw from this urn:

library(tidyverse)
library(dslabs)
take_poll(25)

Think about how you would construct your interval based on the data shown above.

We have just described a simple sampling model for opinion polls. The beads inside the urn represent the individuals that will vote on election day. Those that will vote for the Republican candidate are represented with red beads and the Democrats with the blue beads. For simplicity, assume there are no other colors. That is, that there are just two parties: Republican and Democratic.

Populations, samples, parameters, and estimates

We want to predict the proportion of blue beads in the urn. Let’s call this quantity $p$, which then tells us the proportion of red beads $1-p$, and the spread $p-(1-p)$, which simplifies to $2p-1$.

In statistical textbooks, the beads in the urn are called the population. The proportion of blue beads in the population $p$ is called a parameter. The 25 beads we see in the previous plot are called a sample. The task of statistical inference is to predict the parameter $p$ using the observed data in the sample.

This is a different approach

Above, we knew what was in the urn (the probability distribution) and calculated things from that. Now, we’re in a world where we don’t know what’s in the urn, but we can take samples. The two will meet in the middle – we can caluclate sample versions of things, and see what we can infer about population values.

Here we go!

Can we do this with the 25 observations above? It is certainly informative. For example, given that we see 13 red and 12 blue beads, it is unlikely that $p$ > .9 or $p$ < .1. But are we ready to predict with certainty that there are more red beads than blue in the jar?

We want to construct an estimate of $p$ using only the information we observe. An estimate should be thought of as a summary of the observed data that we think is informative about the parameter of interest. It seems intuitive to think that the proportion of blue beads in the sample $0.48$ must be at least related to the actual proportion $p$. But do we simply predict $p$ to be 0.48? First, remember that the sample proportion is a random variable. If we run the command take_poll(25) four times, we get a different answer each time, since the sample proportion is a random variable.

Note that in the four random samples shown above, the sample proportions range from 0.44 to 0.60. By describing the distribution of this random variable, we will be able to gain insights into how good this estimate is and how we can make it better.

The sample average

Conducting an opinion poll is being modeled as taking a random sample from an urn. We are proposing the use of the proportion of blue beads in our sample as an estimate of the parameter $p$. Once we have this estimate, we can easily report an estimate for the spread $2p-1$, but for simplicity we will illustrate the concepts for estimating $p$. We will use our knowledge of probability to defend our use of the sample proportion and quantify how close we think it is from the population proportion $p$.

We start by defining the random variable $X$ as: $X=1$ if we pick a blue bead at random and $X=0$ if it is red. This implies that the population is a list of 0s and 1s. If we sample $N$ beads, then the average of the draws $X_1,\dots ,X_N$ is equivalent to the proportion of blue beads in our sample. This is because adding the $X$s is equivalent to counting the blue beads and dividing this count by the total $N$ is equivalent to computing a proportion. We use the symbol $\overline{X}$ to represent this average. In general, in statistics textbooks a bar on top of a symbol means the average. The theory we just learned about the sum of draws becomes useful because the average is a sum of draws multiplied by the constant $1/N$:

$$\overline{X}=1/N\times \sum _{i=1}^N{X}_i$$

For simplicity, let’s assume that the draws are independent: after we see each sampled bead, we return it to the urn. In this case, what do we know about the distribution of the sum of draws? First, we know that the expected value of the sum of draws is $N$ times the average of the values in the urn. We know that the average of the 0s and 1s in the urn must be $p$, the proportion of blue beads.

Here we encounter an important difference with what we did in the Probability chapter: we don’t know what is in the urn. We know there are blue and red beads, but we don’t know how many of each. This is what we want to find out: we are trying to estimate $p$.

Parameters

Just like we use variables to define unknowns in systems of equations, in statistical inference we define parameters to define unknown parts of our models. In the urn model which we are using to mimic an opinion poll, we do not know the proportion of blue beads in the urn. We define the parameters $p$ to represent this quantity. $p$ is the average of the urn because if we take the average of the 1s (blue) and 0s (red), we get the proportion of blue beads. Since our main goal is figuring out what is $p$, we are going to estimate this parameter.

The ideas presented here on how we estimate parameters, and provide insights into how good these estimates are, extrapolate to many data science tasks. For example, we may want to determine the difference in health improvement between patients receiving treatment and a control group. We may ask, what are the health effects of smoking on a population? What are the differences in racial groups of fatal shootings by police? What is the rate of change in life expectancy in the US during the last 10 years? All these questions can be framed as a task of estimating a parameter from a sample.

Polling versus forecasting

Before we continue, let’s make an important clarification related to the practical problem of forecasting the election. If a poll is conducted four months before the election, it is estimating the $p$ for that moment and not for election day. The $p$ for election night might be different since people’s opinions fluctuate through time. The polls provided the night before the election tend to be the most accurate since opinions don’t change that much in a day. However, forecasters try to build tools that model how opinions vary across time and try to predict the election night results taking into consideration the fact that opinions fluctuate. We will describe some approaches for doing this in a later section.

Properties of our estimate: expected value and standard error

To understand how good our estimate is, we will describe the statistical properties of the random variable defined above: the sample proportion $\overline{X}$. Remember that $\overline{X}$ is the sum of independent draws so the rules we covered in the probability chapter apply.

Using what we have learned, the expected value of the sum $N\overline{X}$ is $N\times$ the average of the urn, $p$. So dividing by the non-random constant $N$ gives us that the expected value of the average $\overline{X}$ is $p$. We can write it using our mathematical notation:

$$\text{E}(\overline{X})=p$$

We can also use what we learned to figure out the standard error: the standard error of the sum is $\sqrt{N}\times$ the standard deviation of the urn. Can we compute the standard error of the urn? We learned a formula that tells us that it is $(1-0)\sqrt{p(1-p)}$ = $\sqrt{p(1-p)}$. Because we are dividing the sum by $N$, we arrive at the following formula for the standard error of the average:

$$\text{SE}(\overline{X})=\sqrt{p(1-p)/N}$$

This result reveals the power of polls. The expected value of the sample proportion $\overline{X}$ is the parameter of interest $p$ and we can make the standard error as small as we want by increasing $N$. The law of large numbers tells us that with a large enough poll, our estimate converges to $p$.

If we take a large enough poll to make our standard error about 1%, we will be quite certain about who will win. But how large does the poll have to be for the standard error to be this small?

One problem is that we do not know $p$, so we can’t compute the standard error. However, for illustrative purposes, let’s assume that $p=0.51$ and make a plot of the standard error versus the sample size $N$:

From the plot we see that we would need a poll of over 10,000 people to get the standard error that low. We rarely see polls of this size due in part to costs. From the Real Clear Politics table, we learn that the sample sizes in opinion polls range from 500-3,500 people. For a sample size of 1,000 and $p=0.51$, the standard error is:

sqrt(p*(1-p))/sqrt(1000)

[1] 0.01580823

or 1.5 percentage points. So even with large polls, for close elections, $\overline{X}$ can lead us astray if we don’t realize it is a random variable. Nonetheless, we can actually say more about how close we get the $p$.

Central Limit Theorem

The Central Limit Theorem (CLT) tells us that the distribution function for a sum of draws is approximately normal. You also may recall that dividing a normally distributed random variable by a constant is also a normally distributed variable. This implies that the distribution of $\overline{X}$ is approximately normal.

In summary, we have that $\overline{X}$ has an approximately normal distribution with expected value $p$ and standard error $\sqrt{p(1-p)/N}$.

Now how does this help us? Suppose we want to know what is the probability that we are within 1% from $p$. We are basically asking what is

$$\text{Pr}(|\overline{X}-p|\le .01)$$ which is the same as:

$$\text{Pr}(\overline{X}\le p+.01)-\text{Pr}(\overline{X}\le p-.01)$$

Can we answer this question? We can use the mathematical trick we learned in the previous lecture. Subtract the expected value and divide by the standard error to get a standard normal random variable, call it $Z$, on the left. Since $p$ is the expected value and $\text{SE}(\overline{X})=\sqrt{p(1-p)/N}$ is the standard error we get:

$$\text{Pr}(Z\le \frac{.01}{\text{SE}(\overline{X})})-\text{Pr}(Z\le -\frac{.01}{\text{SE}(\overline{X})})$$

One problem we have is that since we don’t know $p$, we don’t know $\text{SE}(\overline{X})$. But it turns out that the CLT still works if we estimate the standard error by using $\overline{X}$ in place of $p$. We say that we plug-in the estimate. Our estimate of the standard error is therefore:

$$\widehat{\text{SE}}(\overline{X})=\sqrt{\overline{X}(1-\overline{X})/N}$$ In statistics textbooks, we use a little hat to denote estimates. The estimate can be constructed using the observed data and $N$.

Now we continue with our calculation, but dividing by $\widehat{\text{SE}}(\overline{X})=\sqrt{\overline{X}(1-\overline{X})/N})$ instead. In our first sample we had 12 blue and 13 red so $\overline{X}=0.48$ and our estimate of standard error is:

x_hat <- 0.48
se <- sqrt(x_hat*(1-x_hat)/25)
se

[1] 0.09991997

And now we can answer the question of the probability of being close to $p$. The answer is:

pnorm(0.01/se) - pnorm(-0.01/se)

[1] 0.07971926

Therefore, there is a small chance that we will be close. A poll of only $N=25$ people is not really very useful, at least not for a close election.

Margin of Error

Earlier we mentioned the margin of error. Now we can define it because it is simply 1.96 times the standard error, which we can now estimate. In our case it is:

1.96*se

[1] 0.1958431

Why do we multiply by 1.96? Because if you ask what is the probability that we are within 1.96 standard errors from $p$, we get:

$$\text{Pr}(Z\le 1.96\text{SE}(\overline{X})/\text{SE}(\overline{X}))-\text{Pr}(Z\le -1.96\text{SE}(\overline{X})/\text{SE}(\overline{X}))$$ which is:

$$\text{Pr}(Z\le 1.96)-\text{Pr}(Z\le -1.96)$$

which we know is about 95%:

pnorm(1.96)-pnorm(-1.96)

[1] 0.9500042

Hence, there is a 95% probability that $\overline{X}$ will be within $1.96\times \widehat{SE}(\overline{X})$, in our case within about 0.2, of $p$. Note that 95% is somewhat of an arbitrary choice and sometimes other percentages are used, but it is the most commonly used value to define margin of error. We often round 1.96 up to 2 for simplicity of presentation.

In summary, the CLT tells us that our poll based on a sample size of $25$ is not very useful. We don’t really learn much when the margin of error is this large. All we can really say is that the popular vote will not be won by a large margin. This is why pollsters tend to use larger sample sizes.

From the table above, we see that typical sample sizes range from 700 to 3500. To see how this gives us a much more practical result, notice that if we had obtained a $\overline{X}$=0.48 with a sample size of 2,000, our standard error $\widehat{\text{SE}}(\overline{X})$ would have been 0.0111714. So our result is an estimate of 48% with a margin of error of 2%. In this case, the result is much more informative and would make us think that there are more red balls than blue. Keep in mind, however, that this is hypothetical. We did not take a poll of 2,000 since we don’t want to ruin the competition.

A Monte Carlo simulation

(Optional) Suppose we want to use a Monte Carlo simulation to corroborate the tools we have built using probability theory. To create the simulation, we would write code like this:

B <- 10000
N <- 1000
x_hat <- replicate(B, {
  x <- sample(c(0,1), size = N, replace = TRUE, prob = c(1-p, p))
  mean(x)
})

The problem is, of course, we don’t know p. We could construct an urn like the one pictured above and run an analog (without a computer) simulation. It would take a long time, but you could take 10,000 samples, count the beads and keep track of the proportions of blue. We can use the function take_poll(n=1000) instead of drawing from an actual urn, but it would still take time to count the beads and enter the results.

One thing we therefore do to corroborate theoretical results is to pick one or several values of p and run the simulations. Let’s set p=0.45. We can then simulate a poll:

p <- 0.45
N <- 1000

x <- sample(c(0,1), size = N, replace = TRUE, prob = c(1-p, p))
x_hat <- mean(x)

In this particular sample, our estimate is x_hat. We can use that code to do a Monte Carlo simulation:

B <- 10000
x_hat <- replicate(B, {
  x <- sample(c(0,1), size = N, replace = TRUE, prob = c(1-p, p))
  mean(x)
})

To review, the theory tells us that $\overline{X}$ is approximately normally distributed, has expected value $p=$ 0.45 and standard error $\sqrt{p(1-p)/N}$ = 0.0157321. The simulation confirms this:

mean(x_hat)

[1] 0.4500761

sd(x_hat)

[1] 0.01579523

A histogram and qq-plot confirm that the normal approximation is accurate as well:

Of course, in real life we would never be able to run such an experiment because we don’t know $p$. But we could run it for various values of $p$ and $N$ and see that the theory does indeed work well for most values. You can easily do this by re-running the code above after changing p and N.

The spread

The competition is to predict the spread, not the proportion $p$. However, because we are assuming there are only two parties, we know that the spread is $p-(1-p)=2p-1$. As a result, everything we have done can easily be adapted to an estimate of $2p-1$. Once we have our estimate $\overline{X}$ and $\widehat{\text{SE}}(\overline{X})$, we estimate the spread with $2\overline{X}-1$ and, since we are multiplying by 2, the standard error is $2\widehat{\text{SE}}(\overline{X})$. Note that subtracting 1 does not add any variability so it does not affect the standard error.

For our 25 item sample above, our estimate $p$ is .48 with margin of error .20 and our estimate of the spread is 0.04 with margin of error .40. Again, not a very useful sample size. However, the point is that once we have an estimate and standard error for $p$, we have it for the spread $2p-1$.

Bias: why not run a very large poll?

For realistic values of $p$, say from 0.35 to 0.65, if we run a very large poll with 100,000 people, theory tells us that we would predict the election perfectly since the largest possible margin of error is around 0.3%:

One reason is that running such a poll is very expensive. Another possibly more important reason is that theory has its limitations. Polling is much more complicated than picking beads from an urn. Some people might lie to pollsters and others might not have phones. But perhaps the most important way an actual poll differs from an urn model is that we actually don’t know for sure who is in our population and who is not. How do we know who is going to vote? Are we reaching all possible voters? Hence, even if our margin of error is very small, it might not be exactly right that our expected value is $p$. We call this bias. Historically, we observe that polls are indeed biased, although not by that much. The typical bias appears to be about 1-2%. This makes election forecasting a bit more interesting and we will talk about how to model this in our Assignment for this week.

Footnotes

1. https://en.wikipedia.org/w/index.php?title=Financial_crisis_of_2007%E2%80%932008

2. https://en.wikipedia.org/w/index.php?title=Security_(finance)

3. https://en.wikipedia.org/w/index.php?title=Poisson_distribution

4. http://www.realclearpolitics.com

5. http://www.realclearpolitics.com/epolls/2016/president/us/general_election_trump_vs_clinton-5491.html