Bootstrap method

“Taking repeated samples from the population” is a theoretical concept and hard to grasp since we don’t know the population. But a good substitute for this in practice is taking repeated samples with replacement from the sample. This is called the bootstrap method.

Taking samples with and without replacement. Consider taking a sample of 10 balls from a bingo ball cage out of 100 ping pong balls that are numbered from 1 to 100. To take the sample, select one ball at a time.

• With replacement: You put the ball back in the cage before you select the next ball. (In other words, the ball can be selected multiple times.)

• Without replacement: Once you select a ball, don’t put it back in the cage. (The ball can only be selected at most once.)

Let’s take 100 bootstrap samples from the nminer dataset. (Usually you’d take thousands of samples but let’s keep the computation time low.)

The workhorse here is the sample_n() function: we are taking samples from nminer, the same size of the original dataset, with replacement.

set.seed(24601)
N <- 100
boot_samples <- tibble()
for(i in 1:N){
  boot_samples <- bind_rows(
    boot_samples, 
    tibble(
      n = i, 
      sample_n(nminer, size = nrow(nminer), replace = TRUE)
    )
  ) |> 
    select(n, Minerab, Eucs)
}

Plots of the first 6 bootstrap samples, along with the fitted mean lines. Note that the fitted lines (and the \(\hat{\beta}_0\) and \(\hat{\beta}_1\) values) vary according to which observations are in the sample.

boot_samples |> 
  filter(n <= 6) |> 
  ggplot(
    mapping = aes(x = Eucs, y = Minerab)
  ) + 
  geom_point() + 
  geom_smooth(
    method = "glm", method.args = list(family = poisson(link = "log")), se = FALSE
  ) + 
  facet_wrap(~n)

This fits the model for each of the 100 samples and extracts the \(\hat{\beta}_1\) values.

beta1hat <- rep(0, N)
for(i in 1:N){
  beta1hat[i] <- glm(Minerab ~ Eucs, family = poisson(link="log"), 
                       data = boot_samples |> filter(n == i)
                     ) |> 
    tidy() |> 
    filter(term == "Eucs") |> 
    pull(estimate)
}

Histogram of \(\hat{\beta}_1\) values. This is a representation of the sampling distribution of \(\hat{\beta}_1\).

tibble(beta1hat = beta1hat) |> 
  ggplot(
    mapping = aes( x= beta1hat)
  ) + 
  geom_histogram(bins=30)

Bootstrap estimate of the standard error.

sd(beta1hat)

## [1] 0.02140763

Compare with the tidy() std.error value of 0.0124. These are not all that close but in the same ballpark.

Fisher information and standard errors

The standard errors given by tidy() come from the log-likelihood function – specifically, the negative second derivative of the log-likelihood function, which is called the Fisher information.

This measures the “sharpness” of the peak of the log-likelihood function (See Fig 4.4 in Section 4.5.2).

A sharper peak means a better estimate with more information and, therefore, less variation. The standard error comes from the inverse of this information (1 divided by sqrt(information)). (See Section 4.5.3)

MLEs follow a normal distribution

Before we discuss the test statistic and the p-value of \(\hat{\beta}_1\), it helps to know the distribution of \(\hat{\beta}_1\).

See Section 4.9.2: Properties of maximum likelihood estimates (like \(\hat{\beta}_1\)). #5 says: “MLEs are asymptotically normally distributed” with

• mean = true parameter value (\(\beta_1\))

• standard deviation calculated from Fisher information as described above.

Let’s return to the other quantities given by with tidy() table.

tidy(fit)

## # A tibble: 2 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   -0.876    0.283      -3.10 1.95e- 3
## 2 Eucs           0.114    0.0124      9.17 4.77e-20

Hypothesis test for \(\beta_1\)

The test statistic and p-value are calculated within the hypothesis test between two competing hypotheses:

• Null hypothesis: \(H_0: \beta_1\)= 0 (no relationship between Minerab and Eucs)

• Alternative hypothesis: \(H_a: \beta_1\neq 0\) (some relationship between Minerab and Eucs)

This test measures how compatible the data and the null hypothesis are. (Test statistic and p-value are based on the null hypothesis that \(H_0: \beta_1= 0\))

Test statistic is a z score

Recall the formula for z scores used in normal calculations:

\[z = \frac{value - mean}{sd}\]

Interpretation of z-score: How standard deviations the value is above or below the mean.

Focusing on the (normal) distribution of \(\hat{\beta}_1\) under \(H_0: \beta_1\)=0, the corresponding z-score formula for \(\hat{\beta}_1\) calculates the test statistic.

• value = \(\hat{\beta}_1\)

• mean = \(\beta_1\) (= 0 under \(H_0\))

• sd estimated by the standard error.

\[ z = \frac{\hat{\beta}_1 - 0}{SE_{\hat{\beta}_1}}\]

Our test statistic is z = 9.17

Interpretation of test statistic: The \(\hat{\beta}_1\) is 9.17 standard deviations above what we expect if \(H_0\) is true.

Assuming \(H_0:\beta_1=0\), should we expect this test statistic (based on what we know about the normal distribution)? No! To be 9 standard deviations away from the mean is a ton for the normal distribution (We said earlier that more than 2 or 3 constitutes an outlier).

The test statistic \(z\) follows a normal distribution with a mean 0 and standard deviation 1. (standard normal distribution).