--- title: The Bootstrap author: 36-402 date: 8 February 2024 output: slidy_presentation: math_method: engine: mathjax url: https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js --- ```{r, include=FALSE} # General set up # Load packages we may use later library(knitr) # Set knitr options for knitting code into the report: # - Save results so that code blocks aren't re-run unless code changes (cache), # _or_ a relevant earlier code block changed (autodep), but don't re-run if the # only thing that changed was the comments (cache.comments) # - Don't clutter R output with messages or warnings (message, warning) # This _will_ leave error messages showing up in the knitted report opts_chunk$set(cache=TRUE, autodep=TRUE, cache.comments=FALSE, message=FALSE, warning=FALSE) set.seed(2024-02-08) ``` ## Last Week Ugly photo omitted, but this slide preserved to keep numbering the same ## The Big Picture \[ \newcommand{\Expect}[1]{\mathbf{E}\left[ #1 \right]} \newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]} \newcommand{\Prob}[1]{\mathrm{Pr}\left( #1 \right)} \newcommand{\Probwrt}[2]{\mathrm{Pr}_{#2}\left( #1 \right)} \] 1. Knowing the sampling distribution of a statistic tells us about statistical uncertainty + standard error + bias + confidence sets + $p$ values 2. The bootstrap principle: _approximate_ the sampling distribution by _simulating_ from a good model of the data, and treating the simulated data just like the real data 3. Sometimes we simulate from the model we're estimating + (model-based or "parametric" bootstrap) 4. Sometimes we simulate by re-sampling the original data + (resampling or "nonparametric" bootstrap) 5. Stronger assumptions $\Rightarrow$ less uncertainty _if we're right_ ## Statistical Uncertainty > - Re-run the experiment (survey, census, ...) and get different data > - $\therefore$ everything we calculate from data would change from run to run > + estimates > + test statistics > + $p$-values > + predicted values > + policy recommendations > - This variability = **statistical uncertainty** > - Quantifying statistical uncertainty = honesty about what we actually know ## Measures of Uncertainty - Standard error = standard deviation of an estimator + (could instead use median absolute deviation, etc.) - $p$-value = Probability we'd see a signal this big if there was just noise - Confidence region = All the parameter values we can't reject at low error rates + **Either** the true parameter is in the confidence region + **or** we are very unlucky + (**or** our model is wrong) - etc., etc. ## Statistical Uncertainty Comes from the Sampling Distribution - Data $X \sim P_X$ for some unknown true distribution $P_X$ - We calculate a statistic $T = \tau(X)$ so it has distribution $P_{T}$, its **sampling distribution** - If we knew $P_{T}$, we could calculate + $\Var{T}$ (and so standard error $=\sqrt{\Var{T}}$) + $\Expect{T}$ (and so bias) + quantiles (and so confidence intervals or $p$-values), etc. ## The Difficulties > - Difficulty 1: Most of the time, $P_{X}$ is very complicated > - Difficulty 2: Most of the time, $\tau$ is a very complicated function > - $\therefore$ We couldn't solve for $P_T$ > - Difficulty 3: Actually, we don't know $P_X$ > - Upshot: We _really_ don't know $P_{T}$ and can't use it to calculate anything ## The Solutions - Classically ($\approx 1900$--$\approx 1975$): Restrict the model and the statistic until you can calculate the sampling distribution, at least for very large $n$ + Paradigm: $t$ distribution for linear models with nice Gaussian noise - Modern ($\approx 1975$--): Use complex models and statistics, but simulate calculating the statistic on the model ## The Monte Carlo Principle > - Generate / simulate $\tilde{X}$ from $P_X$ > - Set $\tilde{T} = \tau(\tilde{X})$ > - Repeat many times > - Use the simulated distribution of the $\tilde{T}$ to approximate $P_{T}$ > + (As a general method, invented by Enrico Fermi in the 1930s, spread through the Manhattan Project) > - Still needs $P_X$ > + Works in HW 3 because we're testing a fixed model ## The Bootstrap Principle 1. Find a good estimate $\hat{P}$ for $P_{X}$ 2. Simulate $\tilde{X}$ from $\hat{P}$, set $\tilde{T} = \tau(\tilde{X})$ 3. Use the simulated distribution of the $\tilde{T}$ to approximate $P_{T}$ + "Pull yourself up by your bootstraps": use $\hat{P}$ to get at uncertainty in itself + Invented by Bradley Efron in the 1970s - First step: find a good estimate $\hat{P}$ for $P_{X}$ ## Model-based Bootstrap If we are using a model, our best guess at $P_{X}$ is $P_{X,\hat{\theta}}$, with our best estimate $\hat{\theta}$ of the parameters #### The Model-based Bootstrap - Get data $X$, estimate $\hat{\theta}$ from $X$ - Repeat $b$ times: + Simulate $\tilde{X}$ from $P_{X,\hat{\theta}}$ (simulate data of same size/"shape" as real data) + Calculate $\tilde{T} = \tau(\tilde{X}$) (treat simulated data the same as real data) - Use simulated distribution of $\tilde{T}$ as $P_{T}$ ## Example: Is Karakedi overweight? ```{r, fig.retina=NULL, out.width=600, echo=FALSE} knitr::include_graphics("kara-as-snookums.png") ``` Concretely: Is she over the 95th percentile of body mass for adult cats? ## Example (cont'd.) ```{r} library(MASS); data(cats); summary(cats) (q95.gaussian <- qnorm(0.95, mean=mean(cats$Bwt), sd=sd(cats$Bwt))) ``` ## Example (cont'd.) Simulate from fitted Gaussian; bundle up estimating 95th percentile into a function ```{r} rcats.gaussian <- function() { rnorm(n=nrow(cats), mean=mean(cats$Bwt), sd=sd(cats$Bwt)) } est.q95.gaussian <- function(x) { m <- mean(x) s <- sd(x) return(qnorm(0.95,mean=m,sd=s)) } ``` ## Example (cont'd.) Simulate, plot the sampling distribution from the simulations ```{r} sampling.dist.gaussian <- replicate(1000, est.q95.gaussian(rcats.gaussian())) plot(hist(sampling.dist.gaussian,breaks=50, plot=FALSE), freq=FALSE) lines(density(sampling.dist.gaussian),lwd=2) abline(v=q95.gaussian, lty="dashed", lwd=4) ``` ## Example (cont'd.) Find standard error and a (crude) confidence interval ```{r} sd(sampling.dist.gaussian) quantile(sampling.dist.gaussian,c(0.025,0.975)) ``` ## Model Checking - As always, if the model isn't right, relying on the model is asking for trouble - Is the Gaussian a good model for cats' weights? ## Example (re-cont'd.) Compare histogram to fitted Gaussian density and to a smooth density estimate ```{r, fig.align="center", echo=FALSE} plot(hist(cats$Bwt,plot=FALSE),freq=FALSE, main="", xlab="Body weight (kg)", xlim=c(0, 1.05*max(cats$Bwt))) curve(dnorm(x,mean=mean(cats$Bwt),sd=sd(cats$Bwt)), add=TRUE,col="purple") lines(density(cats$Bwt),lty="dashed") legend("topleft", legend=c("Histogram", "Smoothed density", "Gaussian"), lty=c("solid", "dashed", "solid"), col=c("black", "black", "purple")) ``` ## Resampling Difficulty: We might not have a trust-worthy model Resource: We do have data, which tells us a lot about the distribution Solution: **Resampling**, treat the sample like a whole population ## The Resampling ("Nonparameteric") Bootstrap - Get data $X$, containing $n$ samples, find $T=\tau(X)$ - Repeat $b$ times: + Generate $\tilde{X}$ by drawing $n$ samples from $X$ _with replacement_ (resample the data) + Calculate $\tilde{T} = \tau(\tilde{X})$ (treat simulated data the same as real data) - Use simulated distribution of $\tilde{T}$ as $P_{T}$ (See backup for $\hat{P}$ implicit here) ## Example, Take 2 Model-free estimate of the 95th percentile = 95th percentile of the data ```{r} (q95.np <- quantile(cats$Bwt,0.95)) ``` How precise is that? ## Example, Take 2 Resampling, re-estimating, and finding sampling distribution ```{r} resample <- function(x) { sample(x,size=length(x),replace=TRUE) } est.q95.np <- function(x) { quantile(x,0.95) } ``` ## Example, Take 2 ```{r, fig.align="center"} sampling.dist.np <- replicate(1000, est.q95.np(resample(cats$Bwt))) plot(density(sampling.dist.np), main="", xlab="Body weight (kg)") abline(v=q95.np,lty=2) ``` ## Example, Take 2 (cont'd) standard error, bias, CI ```{r} sd(sampling.dist.np) mean(sampling.dist.np - q95.np) quantile(sampling.dist.np,c(0.025,0.975)) ``` ## Bootstrapping Regressions - A regression is a model for $Y$ conditional on $X$ \[ Y= \mu(X) + \epsilon, ~ \Expect{\epsilon|X} = 0 \] - Silent on the distribution of $X$, so how do we simulate? - Options, putting less and less trust in the model: + Hold $x_i$ fixed, set $\tilde{y}_i = \hat{\mu}(x_i) + \tilde{\epsilon}_i$ from model's estimated noise distribution (e.g., Gaussian or Poisson) + Hold $x_i$ fixed, set $\tilde{y}_i = \hat{\mu}(x_i) + \tilde{\epsilon}_i$, noise _resampled from the residuals_ + Resample $(x_i, y_i)$ pairs (resample data-points or resample cases) ## Different Regression Bootstraps - Resampling residuals works as long as the noise is IID + Noise could be Gaussian... + Or noise could be very non-Gaussian + Noise does need have same distribution everywhere + Dubious unless the regression model is right - Resampling whole cases works as long as observations are IID + noise needn't be independent of $X$ + needn't be Gaussian + regression model needn't be right ## Cats' Hearts `cats` has weights for cats' hearts, as well as bodies ![](cat-heart.jpg) [(Much cuter than any photo of real cats' hearts)](http://yaleheartstudy.org/site/wp-content/uploads/2012/03/cat-heart1.jpg) How does heart weight relate to body weight? (Useful when Kara's vet needed to know how much heart medicine to prescribe) ## Cats' Hearts (cont'd) ```{r, fig.align="center"} plot(Hwt~Bwt, data=cats, xlab="Body weight (kg)", ylab="Heart weight (g)") cats.lm <- lm(Hwt ~ Bwt, data=cats) abline(cats.lm) ``` ## Cats' Hearts (cont'd) - Coefficients and standard confidence intervals: ```{r} coefficients(cats.lm) confint(cats.lm) ``` - These confidence intervals assume IID Gaussian noise ## Cats' Hearts (cont'd) The residuals don't look that Gaussian: ```{r, echo=FALSE, fig.align="center"} par(mfrow=c(1,2)) plot(cats$Bwt,residuals(cats.lm), xlab="Body weight (kg)", ylab="Residuals for heart weight (g)") plot(density(residuals(cats.lm)), main="", xlab="Residuals for heart weight (g)") par(mfrow=c(1,1)) ``` ## Cats' Hearts (cont'd) Resample residuals: ```{r} sim.cats.resids <- function() { new.cats <- cats noise <- resample(residuals(cats.lm)) new.cats$Hwt <- fitted(cats.lm) + noise return(new.cats) } ``` Re-estimate on new data: ```{r} coefs.cats.lm <- function(df) { fit <- lm(Hwt~Bwt, data=df) return(coefficients(fit)) } ``` ## Cats' Hearts (cont'd) Re-sample to get CIs: ```{r} cats.lm.samp.dist.resids <- replicate(1000, coefs.cats.lm(sim.cats.resids())) t(apply(cats.lm.samp.dist.resids, 1, quantile, c(0.025,0.975))) ``` ## Cats' Hearts (cont'd) Try resampling whole rows: ```{r} resample.data.frame <- function(df) { return(df[resample(1:nrow(df)),]) } ``` ```{r} cats.lm.samp.dist.cases <- replicate(1000, coefs.cats.lm(resample.data.frame(cats))) t(apply(cats.lm.samp.dist.cases,1,quantile,c(0.025,0.975))) ``` ## Comparison - We now have three sets of confidence intervals 1. "Conventional"/Gaussian: ```{r, echo=FALSE} confint(cats.lm) ``` 2. By resampling residuals: ```{r, echo=FALSE} t(apply(cats.lm.samp.dist.resids,1,quantile,c(0.025,0.975))) ``` 3. By resampling rows/cases: ```{r, echo=FALSE} t(apply(cats.lm.samp.dist.cases,1,quantile,c(0.025,0.975))) ``` - Why do the intervals keep getting wider? - Which of these is most reliable? ## Sources of Error in Bootstrapping - **Simulation** Using only $b$ bootstrap replicates + Make this small by letting $b\rightarrow\infty$ + Costs computing time + Diminishing returns: error is generally $\propto 1/\sqrt{b}$ + There are tricks for speeding up simulations, re-using parts of them, etc. - **Approximation** Using $\hat{P}$ instead of $P_{X}$ + Make this small by careful statistical modeling - **Estimation** Only a finite number of samples + Make this small by being careful about what we simulate (example in the backup slides) - For fixed $n$, resampling generally shows more uncertainty than model-based bootstrap - Resampling is less vulnerable to modeling mistakes ## Summing Up - We need to know the sampling distribution of our statistic $T = \tau(X)$ - If we could simulate $\tilde{X}$ from the true distribution of $X$, we'd use the distribution of $\tau(\tilde{X})$ as the sampling distribution ("Monte Carlo") - We simulate from an _estimate_ of the distribution of $X$ instead ("bootstrap") + The simulation could be based on a model + Or on re-sampling the observations + Or some combination ## Backup: What's $\hat{P}$ in Resampling? (After discussion in class) $\hat{P}$ says: - Put probability $1/n$ on each observed value of $X$ - _and_ new values are independent and identically distributed (IID) Symbolically: \[ \hat{P} = \frac{1}{n}\sum_{i=1}^{n}{\delta(x-x_i) \] where the **delta distribution** or **delta function** puts probability 1 on $0$ (and probability 0 on every other value) (Think of an infinitely tall, infinitely narrow Gaussian centered at $0$) ## Backup: Improving on the Crude Confidence Interval Crude CI use distribution of $\tilde{\theta}$ under $\hat{\theta}$ But really we want the distribution of $\hat{\theta}$ under $\theta$ Mathematical observation: Generally speaking, (distributions of) $\tilde{\theta} - \hat{\theta}$ is closer to $\hat{\theta}-\theta_0$ than $\tilde{\theta}$ is to $\hat{\theta}$ ("centering" or "pivoting") ## Backup: The Basic, Pivotal CI quantiles of $\tilde{\theta}$ = $q_{\alpha/2}, q_{1-\alpha/2}$ \[ \begin{eqnarray*} 1-\alpha & = & \Probwrt{q_{\alpha/2} \leq \tilde{\theta} \leq q_{1-\alpha/2}}{\hat{\theta}} \\ & = & \Probwrt{q_{\alpha/2} - \hat{\theta} \leq \tilde{\theta} - \hat{\theta} \leq q_{1-\alpha/2} - \hat{\theta}}{\hat{\theta}} \\ & \approx & \Probwrt{q_{\alpha/2} - \hat{\theta} \leq \hat{\theta} - \theta_0 \leq q_{1-\alpha/2} - \hat{\theta}}{\theta_0}\\ & = & \Probwrt{q_{\alpha/2} - 2\hat{\theta} \leq -\theta_0 \leq q_{1-\alpha/2} - 2\hat{\theta}}{\theta_0}\\ & = & \Probwrt{2\hat{\theta} - q_{1-\alpha/2} \leq \theta_0 \leq 2\hat{\theta}-q_{\alpha/2}}{\theta_0} \end{eqnarray*} \] Basically: re-center the simulations around the empirical data