17.3 MC-simulated \(p\)-values

Let’s reconsider the 24/7 data set, where we have \(k=7\) observations of ‘heads’ in \(N=24\) tosses of a coin.

# 24/7 data
k_obs <- 7
n_obs <- 24

The question of interest is whether the coin is fair, i.e., whether \(\theta_c = 0.5\). R’s built-in function binom.test calculates a binomial test and produces a \(p\)-value which is calculated precisely (since this is possible and cheap in this case).

binom.test(7,24)

## 
##  Exact binomial test
## 
## data:  7 and 24
## number of successes = 7, number of trials = 24, p-value = 0.06391
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.1261521 0.5109478
## sample estimates:
## probability of success 
##              0.2916667

It is also possible to approximate a \(p\)-value by Monte Carlo simulation. Notice that the definition of a \(p\)-value repeated here from Section 16.2 is just a statement about the probability that a random variable (from which we can take samples with MC simulation) delivers a value below a fixed threshold:

\[ p\left(D_{\text{obs}}\right) = P\left(T^{|H_0} \succeq^{H_{0,a}} t\left(D_{\text{obs}}\right)\right) % = P(\mathcal{D}^{|H_0} \in \{D \mid t(D) \ge t(D_{\text{obs}})\}) \]

So here goes:

# specify how many Monte Carlo samples to take
x_reps <- 500000

# build a vector of likelihoods (= the relevant test statistic)
#   for hypothetical data observations, which are 
#   sampled based on the assumption that H0 is true
lhs <- map_dbl(1:x_reps, function(i) {
  # hypothetical data assuming H0 is true
  k_hyp <- rbinom(1, size = n_obs, prob = 0.5)
  # likelihood of that hypothetical observation
  dbinom(k_hyp, size = n_obs, prob = 0.5)
})

# likelihood (= test statistic) of the observed data
lh_obs = dbinom(k_obs, size = n_obs, prob = 0.5)

# proportion of samples with a lower or equal likelihood than 
#   the observed data 
mean(lhs <= lh_obs) %>% show()

## [1] 0.06414

Monte Carlo sampling for \(p\)-value approximation is always possible, even for cases where we cannot rely on known simplifying assumptions.