Empirical, simulated, and ideal data

All distributions (binomial, normal etc.) that we cover in this course are ideal distributions. Real data are very rarely distributed in exactly the same way. Still, as long as real data are distributed in almost the same way, it is useful to assume that they are well described by some ideal distribution. This allows us to use the same set of mathematical tools for every set of similarly distributed data.

Here, we will compare three distributions:

  1. We will collect real data where we will measure the number of successes in a fixed number of trials.
  2. We will simulate a random sample of binomially distributed data using R.
  3. We will compare both to an ideal binomial distribution constructed using R.

Collect empirical data of the N(successes)/N(trials) variety

Each of you will get a 2-sided and a 3-sided fidget spinner. Spin each spinner 12 times. Each time, stop it by hand and determine whether pink or other than pink side came on top. Keep a tally (you may use this website). We will refer to each of the 12 spins as trials, and the number of times the pink came on top in 12 trials will be considered a single observation.

Let each of us record one observation of the number of successes in 12 trials. Each person needs to submit 1 observation for the 2-sided spinner and 1 observation for the 3-sided spinner.

# record all observations for the 2-sided spinner
spinner2 <- c(8, 8, 4, 5, 9, 5, 5, 6, 7, 8, 6, 7, 4, 4, 3)

# record all observations for the 3-sided spinner
spinner3 <- c(3, 5, 6, 3, 1, 6, 2, 5, 6, 2, 4, 5, 3, 4)

# arrange the data in a data.frame (long format)
bino_data <- data.frame("method" = "spinner",
                        "theoretical_probability" = c( rep(0.5, length(spinner2) ),
                                                       rep(0.33, length(spinner3) ) ),
                        "values" = c(spinner2, spinner3))

# plot the frequencies
ggplot(data = bino_data,
       mapping = aes(x = values,
                     fill = method)) +
  geom_histogram(position = position_dodge(0.7),
                 breaks = seq(from = -0.5, to = 12.5, by = 1),
                 col = "lightgrey") +
  scale_fill_manual(values = "grey40") +
  scale_x_continuous(breaks = 0:12,
                     labels = 0:12,
                     limits = c(0,12)) +
  facet_wrap(~theoretical_probability,
             ncol = 1)
## Warning: Removed 4 rows containing missing values (geom_bar).

Simulate a random sample of binomially distributed data with the same parameters

R allows you to simulate data randomly sampled from a distribution with given parameters. For this, we use “r” (“random sample”) functions from the “d/p/q/r” family. Since we are working with binomial data, we will use the function rbinom().

# simulate performing 15 observations of the number of successes in 12 trials if the probability of success = 0.5
simulated0.5 <- rbinom(n = 15,
                       size = 12,
                       prob = 0.5)
simulated0.5
##  [1] 5 4 8 7 7 5 6 6 7 5 8 7 6 6 4
# simulate performing 15 observations of the number of successes in 12 trials if the probability of success = 1/3
simulated0.33 <- rbinom(n = 15,
                        size = 12,
                        prob = 1/3)
simulated0.33
##  [1] 5 3 4 4 2 2 5 4 6 2 4 7 1 0 3
# add to bino_data
new_bino_data <- data.frame("method" = "simulated",
                            "theoretical_probability" = c( rep(0.5, length(simulated0.5) ),
                                                           rep(0.33, length(simulated0.33) ) ),
                            "values" = c( simulated0.5, simulated0.33 ))

bino_data <- rbind(bino_data, new_bino_data)

# plot all frequencies
ggplot(data = bino_data,
       mapping = aes(x = values,
                     fill = method)) +
  geom_histogram(position = position_dodge(0.7),
                 breaks = seq(from = -0.5, to = 12.5, by = 1),
                 col = "lightgrey") +
  scale_fill_manual(values = c("cornflowerblue", "grey40")) +
  scale_x_continuous(breaks = 0:12,
                     labels = 0:12,
                     limits = c(0,12)) +
  facet_wrap(~theoretical_probability,
             ncol = 1)
## Warning: Removed 8 rows containing missing values (geom_bar).

Generate the frequencies of an ideal binomial distribution with the same parameters

We can use a sister function of rbinom() from the “d/p/q/r” family called dbinom() to generate a PDF of an ideal binomial distribution with the same parameters. The output will be a vector of probabilities for each observation (e.g. probability to observe 0 successes in 12 trials, 1 success, 2 successes etc.) on a scale of 0 to 1. We will need to convert these observation probabilities into real numbers of observations since we are performing a finite number of observations.

# simulate an ideal binomial distribution of the number of successes in 12 trials if the probability of success = 0.5
n_successes_vector <- 0:12

dbinom0.5 <- dbinom(x = n_successes_vector,
                    size = 12,
                    prob = 0.5)
dbinom0.5
##  [1] 0.0002441406 0.0029296875 0.0161132813 0.0537109375 0.1208496094
##  [6] 0.1933593750 0.2255859375 0.1933593750 0.1208496094 0.0537109375
## [11] 0.0161132813 0.0029296875 0.0002441406
# simulate an ideal binomial distribution of the number of successes in 12 trials if the probability of success = 1/3
dbinom0.33 <- dbinom(x = n_successes_vector,
                     size = 12,
                     prob = 1/3)
dbinom0.33
##  [1] 7.707347e-03 4.624408e-02 1.271712e-01 2.119520e-01 2.384460e-01
##  [6] 1.907568e-01 1.112748e-01 4.768921e-02 1.490288e-02 3.311751e-03
## [11] 4.967626e-04 4.516023e-05 1.881676e-06
# by default, all probabilities are on a 0..1 scale
# calculate the most likely number of observations of each number of successes
# if we perform a total of 15 observations
dbinom0.5 <- 15 * dbinom0.5
dbinom0.5
##  [1] 0.003662109 0.043945313 0.241699219 0.805664062 1.812744141 2.900390625
##  [7] 3.383789062 2.900390625 1.812744141 0.805664062 0.241699219 0.043945313
## [13] 0.003662109
dbinom0.33 <- 15 * dbinom0.33
dbinom0.33
##  [1] 1.156102e-01 6.936612e-01 1.907568e+00 3.179280e+00 3.576691e+00
##  [6] 2.861352e+00 1.669122e+00 7.153381e-01 2.235432e-01 4.967626e-02
## [11] 7.451439e-03 6.774035e-04 2.822515e-05
# the numbers of observations must be integers, so round up the values
dbinom0.5 <- round(dbinom0.5,
                   digits = 0)
dbinom0.5
##  [1] 0 0 0 1 2 3 3 3 2 1 0 0 0
dbinom0.33 <- round(dbinom0.33,
                    digits = 0)
dbinom0.33
##  [1] 0 1 2 3 4 3 2 1 0 0 0 0 0

We are halfway there! There is one more important thing we need to change before we can plot our ideal distributions, though. Both our real (e.g. spinner2) and simulated (e.g. simulated0.5) data are simple collections of individual observations. It means that if we observed 1 success in 12 trials twice, these vectors will contain “1” two times. In contrast, the vectors that store our ideal data (dbinom0.5 and dbinom0.33) contain frequencies of each observation in n_successes_vector. If we observed 1 success in 12 trials twice, this information will be stored like so:

0 1 2 3 4 5 6 7 8 9 10 11 12 <– n_successes_vector

x 2 x x x x x x x x x x x <– dbinom...

Therefore, we need to convert dbinom0.5 and dbinom0.33 from frequency vector to simple sets of all the individual observations. For example, 1 2 2 1 0 0 0 0 0 0 0 0 0 should become 0 1 1 2 2 3.

# use a for loop
ideal0.5 <- c() # create an empty vector

for (n_successes in n_successes_vector) {
  ideal0.5 <- c(ideal0.5,
                rep(n_successes, dbinom0.5[n_successes + 1] ) ) # add 1 because the frequency of having 0 successes is stored in the 1st element of dbinom0.5 etc.
}
ideal0.5
##  [1] 3 4 4 5 5 5 6 6 6 7 7 7 8 8 9
ideal0.33 <- c() # create an empty vector

for (n_successes in n_successes_vector) {
  ideal0.33 <- c(ideal0.33,
                 rep(n_successes, dbinom0.33[n_successes + 1] ) ) # add 1 because... see above
}
ideal0.33
##  [1] 1 2 2 3 3 3 4 4 4 4 5 5 5 6 6 7
# add to bino_data
new_bino_data <- data.frame("method" = "ideal",
                            "theoretical_probability" = c( rep(0.5, length(ideal0.5) ),
                                                           rep(0.33, length(ideal0.33) ) ),
                            "values" = c( ideal0.5, ideal0.33 ))
bino_data <- rbind(bino_data, new_bino_data)

# plot all frequencies
ggplot(data = bino_data,
       mapping = aes(x = values,
                     fill = method)) +
  geom_histogram(position = position_dodge(0.7),
                 breaks = seq(from = -0.5, to = 12.5, by = 1),
                 col = "lightgrey") +
  scale_fill_manual(values = c("black", "cornflowerblue", "grey40")) +
  scale_x_continuous(breaks = 0:12,
                     labels = 0:12,
                     limits = c(0,12)) +
  facet_wrap(~theoretical_probability,
             ncol = 1)
## Warning: Removed 8 rows containing missing values (geom_bar).

How likely is each empirical observation in an ideal binomial distribution with given parameters?

We can determine how likely each of the individual observations is given the ideal distribution. We did this earlier in the polio vaccine example. To do this, we can use sister functions of rbinom() from the “d/p/q/r” family called dbinom() (which we used above) and pbinom(). The difference is that dbinom() returns the Y value from the PDF (what is the chance of observing exactly the given value?) and pbinom() returns the Y value from the CDF (what is the chance of observing the given value or a more extreme value?)

# how likely is it to observe *exactly* 4 successes in 12 trials given the ideal binomial distribution with a 0.5 probability of success?
dbinom(x = 4,
       size = 12,
       prob = 0.5)
## [1] 0.1208496
# *exactly* 8 successes?
dbinom(x = 8,
       size = 12,
       prob = 0.5)
## [1] 0.1208496
# 4 *or fewer* successes?
pbinom(q = 4,
       size = 12,
       prob = 0.5)
## [1] 0.1938477
# is it consistent with the eCDF? (note that the values may not match exactly due to rounding)
plot(ecdf(ideal0.5))

# 8 *or fewer* successes?
pbinom(q = 8,
       size = 12,
       prob = 0.5)
## [1] 0.927002
# 8 *or more* successes? (i.e. more than 7 successes)
pbinom(q = 7,
       size = 12,
       prob = 0.5,
       lower.tail = FALSE)
## [1] 0.1938477