---
title: Significance Testing etc.
author: 36-402, Section A
date: 7 February 2019
output: slidy_presentation
bibliography: locusts.bib
---

\[
\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\newcommand{\StdErr}[1]{\mathrm{se}\left[ #1 \right]}
\newcommand{\EstStdErr}[1]{\widehat{\mathrm{se}}\left[ #1 \right]}
\newcommand{\Prob}[1]{\mathbb{P}\left( #1 \right)}
\newcommand{\Var}[1]{\mathbb{V}\left[ #1 \right]}
\]

```{r, include=FALSE}
set.seed(2019)
```

## Linear Regression and Hypothesis Testing (Once More with Feeling) {.smaller}

```{r}
stocks <- read.csv("http://www.stat.cmu.edu/~cshalizi/uADA/19/hw/03/stock_history.csv")
stocks$MAPE <- with(stocks, Price/Earnings_10MA_back)
returns.on.mape.transforms <- lm(Return_10_fwd ~ I(1/MAPE) + MAPE + I(MAPE^2), data=stocks)
summary(returns.on.mape.transforms)
```

## What Is the Model?  What Are the Hypotheses?

\[
R_i = \beta_0 + \beta_1\frac{1}{M_i} + \beta_2 M_i + \beta_3 M_i^2 + \epsilon_i, ~ \Expect{\epsilon|M} = 0
\]
`lm` also assumes  $\epsilon_i$ IID Gaussian, $~ \mathcal{N}(0, \sigma^2)$

- Sampling distribution of $\hat{\beta}_j$ under _all_ these assumptions:
\[
\hat{\beta}_j \sim \mathcal{N}(\beta_j, \StdErr{\beta_j}^2)
\]


#### What Hypotheses Does `lm` Test?

- For $\beta_0$, $H_{00}$: That model is right, and $\beta_0=0$ exactly vs. That model is right, and $\beta_0 \neq 0$
- For $\beta_1$, $H_{01}$: That model is right, and $\beta_1=0$ exactly vs. That model is right, and $\beta_1 \neq 0$
- For $\beta_2$, $H_{02}$: That model is right, and $\beta_2=0$ exactly vs. That model is right, and $\beta_2 \neq 0$
- For $\beta_3$, $H_{03}$: That  model is right, and $\beta_3=0$ exactly vs. That model is right, and $\beta_3 \neq 0$

## How Does R Test These Hypotheses?

- Need a **test statistic** = a function of the data with different distributions under null and alternative
    + Generally, we pick test statistics which tend to be small under the null, and large under the alternative
- R uses the **Wald test**, where the test statistic is
\[
T_j \equiv \frac{\hat{\beta}_j}{\EstStdErr{\hat{\beta}_j}}
\]
    + Works for testing $\beta_j=0$ --- how to modify to test $\beta_j=b \neq 0$?
- We know the sampling distribution of $T_j$ under these modeling assumptions _and_ the null hypothesis $H_{0j}$:
\[
T_j \sim t_{n-4} \approx N(0,1) ~ \mathrm{for\ large}\ n
\]
- We get a certain value of $T_j$ on the data, say $\hat{T}_j$
\[
P_j=\Prob{|T_j| \geq |\hat{T}_j|; H_{0j}}
\]

## What "Statistically Significant" Means

- We can reliably or confidently detect a difference from null in the direction of the alternative
    + We wouldn't see this big a departure from the null, in the direction of the alternative, _if the null were true_
	+ unless we were really unlucky
	+ or unless the model is wrong
- "Statistically detectable" might have been a better name

## What Tends to Make Things Significant?

- $T_j = \hat{\beta}_j/\EstStdErr{\hat{\beta}_j}$
- As $n\rightarrow \infty$, $\hat{\beta}_j \rightarrow \beta_j$ and $\EstStdErr{\hat{\beta}_j} \rightarrow \StdErr{\hat{\beta}_j}$
- So if $\beta_j$ is large, $T_j$ will tend to be large
- What about $\StdErr{\hat{\beta}_j}$?
    + $\StdErr{\hat{\beta}_j} \propto \sigma$
    + $\StdErr{\hat{\beta}_j} \propto 1/\sqrt{n}$
	+ As $\Var{X_j} \uparrow \infty$, $\StdErr{\hat{\beta}_j} \downarrow 0$
	+ As $X_j$ becomes correlated with the other $X_k$'s, $\StdErr{\hat{\beta}_j} \uparrow \infty$ (worst case: collinearity)
	+ (Really: MSE of a linear regression of $X_j$ on the other $X_k$'s)

## What Tends to Make Things Significant?

> - For fixed $n$, the detectably-different-from-zero coefficients belong to variables with
>    + Large true coefficients $\beta_j$
> 	 + Large variances of that variable, $\Var{X_j}$
>	 + Little correlation with the other variables

## What Tends to Make Things Significant?

```{r}
x <- with(na.omit(stocks), cbind(1/MAPE, MAPE, MAPE^2))
colnames(x) <- c("1/MAPE", "MAPE", "MAPE^2")
var(x)
cor(x)
```

## What Makes Things Significant? (cont'd)

> - As $n$ grows, **every** coefficient which isn't **exactly** zero becomes significant
>    + $T_j \rightarrow \pm \infty$, **unless** $\beta_j=0$
>    + $P_j \rightarrow 0$ exponentially fast
>    + "The $P$-value is a measure of sample size"
>    + More [here](http://bactra.org/weblog/1111.html)
> - If $\beta_j=0$, then $P_j \sim \mathrm{Unif}(0,1)$ for **all** $n$
>   + Assuming the model is right
>   + So about 1/20 of really-0 coefficients will look significant

## Star-Gazing Is a Bad Way to Pick Variables {.smaller}

```{r}
n <- 1e4; p <- 40
predictors <- matrix(rnorm(n*p), nrow=n, ncol=p); response<-rnorm(n) # no relationship at all
summary(lm(response ~ predictors))
```

## What's the Point of Significance Testing?

> - You have a real reason to think $\beta_j=0$ is a "live option"
>     + **And** you've included the right covariates
> 	  + **And** you're using the right kind of model
>    + Good evidence when test has power to detect if $H_0$ is false
> - Similarly: you have reasons why $\beta_j=b$ might be worth considering
> - Or: you want to know whether $H_0$ is _close enough_ to true that you can't detect the difference
>    + Especially: goodness-of-fit tests and diagnostics
>	 + "All the hypotheses that could be true" = confidence set


## Confidence Sets

- A confidence set (interval, rectangle, ellipse, ...) is a bet:
    + **Either** the true parameter is in the set, $\theta \in C$
	+ **or** something really improbable happened (1-confidence level)
	+ (**or** the model is wrong)
- Every confidence set can be a test
    + Reject $\theta=\theta_0$ if $\theta_0 \not \in C$
	+ Retain if $\theta_0 \in C$
- Every test gives a confidence set
    + $\theta_0 \in C$ if test doesn't reject $\theta_0$
	+ Wald test for coefficients $\Rightarrow$ usual confidence intervals


## Simulations

- Nothing changes in the logic, just how we calculate the sampling distribution
- If you want to test a hypothesis, run simulations where that hypothesis is true, and see if they look like the data
    + E.g., we simulate with $\beta_j=0$ but all the other $\beta_k$ set to their optimal values
	+ If we didn't know that $\hat{\beta_j} \sim N(\beta_j, \StdErr{\hat{\beta}_j}^2)$ in the linear-Gaussian model, but we did have `rnorm()`, we could do this simulation, and we'd get $P$-values that match `lm`'s calculations
	+ "Looks like" comes down to the test statistic
- In the homework, we simulated the hypothesis that $R_i = \frac{1}{M_i} + \epsilon_i$, with $\epsilon_i$ following a generalized $t$-distribution
    + Used the $t$ rather than normal to get a bit more flexibility in the shape
	+ Also, $t$-distributed noise is used a lot in financial modeling
	+ Many other possibilities for the noise

## Bootstrapping

- In the bootstrap, we don't simulate a fixed model, we craft a simulation so it's close to the data
    + Either by fitting a model to simulate
	+ Or by resampling
- This is good at getting things like standard errors and confidence regions
    + Confidence region = our estimate would be in this set with high probability
	+ Can always invert a confidence region to test a hypothesis
- Advantage of resampling is that it assumes very little about the true distribution
- Advantage of model-based bootstrap is that it uses the data more efficiently, _if_ the model is good


## Summing Up

- When faced with a hypothesis test, ask:
    + What exactly are we testing, and why?
	+ What is the null hypothesis, and what alternative are we testing against?
	+ What test statistic do we use to discriminate between null and alternative?
	+ How do we find the distribution of the test statistic under the null?
	+ Do we have power to detect that alternative with that test statistic?
- Often, what you really want is a confidence set
    + $=$ all the specific hypotheses the data would let you retain
	+ Still need to worry about model accuracy to calculate probabilities
- Simulation is a way of calculating probabilities
    + Can simulate under specific hypotheses
	+ Or simulate under something close to the real distribution (bootstrap)

## References

- The correct line on hypothesis tests, and when, why and how they are useful: @Mayo-error, @Mayo-Cox-frequentist
- $P$-values tend to zero exponentially fast (unless the null is true): @Bahadur-introduces-Bahadur-efficiency, @Bahadur-limit-theorems, @van-der-Vaart-asymptotic-stats
- Crash course in linear regression modeling, including's actually being tested, diagnostics, model-building: [http://www.stat.cmu.edu/~cshalizi/TALR/]
    + Chapters 7, 12 and 16 are especially relevant to today
- All the linear-Gaussian-model theory you could ever want: @Seber-Lee-linear-reg-2nd

## From discussion: "P-Hacking" {.smaller}

- Manipulating the data and/or analysis until you get the significant results you want
    + Or, sometimes, any significance at all (cf. the [example I tried to describe in class](https://arstechnica.com/science/2017/04/the-peer-reviewed-saga-of-mindless-eating-mindless-research-is-bad-too/) )
- Three main ways it happens:
     + Deciding which data points to include or exclude as outliers
	 + Deciding how to define variables --- [Recent example](http://dx.doi.org/10.1038/s41562-018-0506-1) with different measures of "screen time" and "depression"
	 + Deciding which model specification to use
- Andy Gelman calls this the "garden of forking paths" (after the @Borges-ficciones story about alternate histories)
- Brutal solution to $p$-hacking by specification search: data-splitting
    + Randomly split the data into two parts
	+ Do exploration, try out variations, etc., on one part, and _only_ on one part
	+ _After_ you have fixed on a model, estimate it and look at inferential statistics using second part
	+ Do not use the second part for anything else
    + Cf. Exercise 3.4 in Chapter 3 (computationally challenging)

## From discussion: What about the F test?

- R reports an F test at the bottom of its output for `lm()`

#### What hypothesis does this test?

- $H_{0F}$: The linear-Gaussian model is right, and $\beta_1=\beta_2=\beta_3=0$ vs. The linear-Gaussian model is right, and not all slopes are exactly 0
- This is **not** testing "Is the linear-Gaussian model right?"



## Bibliography
