---
title: Testing a parametric regression model
date: 21 February 2019
output: slidy_presentation
---

\[
\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\]

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

## The general problem

- You have a parametric regression model for your data
    + e.g., linear with such-and-such variables
- You are worried that it might be **misspecified**, that the true $\mu(x)$
isn't in the model
- Now that we know nonparametric regression, we can test this


## First case: linear specification is wrong

```{r}
# Make up sample data
x <- sort(runif(300, 0, 3)) # Sorting simplifies plotting later
# Impose true regression function $\log{(x+1)} +$ Gaussian noise:
yg <- log(x+1) + rnorm(length(x), 0, 0.15)
# Bind into a data frame:
gframe <- data.frame(x=x, y=yg)
```

## What this looks like

Plot it, plus the true regression curve
```{r}
plot(y~x, data=gframe, xlab="x", ylab="y", pch=16, cex=0.5)
curve(log(1+x), col="grey", add=TRUE, lwd=4)
```

## What linear fit do we get?

Fit the linear model
```{r}
glinfit <- lm(y~x, data=gframe)
print(summary(glinfit), signif.stars=FALSE, digits=2)
```

(N.B.: everything is crazily significant, but the linear model is wrong.)



## What the mis-specified linear fit looks like {.smaller}

```{r}
plot(x, yg, xlab="x", ylab="y")
curve(log(1+x), col="grey", add=TRUE, lwd=4)
abline(glinfit, lwd=4)
legend("bottomright", legend=c("Truth", "Linear"),
       col=c("grey", "black"), lwd=4)
```

## 

MSE of linear model, in-sample: `r signif(mean(residuals(glinfit)^2),3)`.

We'll need to do that a lot, so make it a function:
```{r}
(mse.residuals <- function(model) { mean(residuals(model)^2) })
```

- OFFLINE: Write comments giving the inputs and outputs of the function `ms.residuals`
- OFFLINE: Check that the function works (how?)

##

Fit the non-parametric alternative; use splines as implemented in `gam()`:
```{r, message=FALSE, warning=FALSE}
library(mgcv)
```

```{r}
gnpr <- gam(y~s(x),data=gframe)
```
Add the fitted values from the spline to the plot:
```{r}
plot(x,yg,xlab="x",ylab="y")
curve(log(1+x),col="grey",add=TRUE,lwd=4)
abline(glinfit, lwd=4)
lines(x,fitted(gnpr),col="blue",lwd=4)
legend("bottomright", legend=c("Truth", "Linear", "Spline"),
       col=c("grey", "black", "blue"), lwd=4)
```

##

Calculate the difference in MSEs:
```{r}
(D.hat <- mse.residuals(glinfit) - mse.residuals(gnpr))
```

- OFFLINE: Why do we believe this code works?

## How do we test the linear model?

- $H_0$: $\mu(x)$ is a linear function vs. $H_A$: $\mu(x)$ is not linear
- Test statistic: Difference in MSEs $=D$
    + If the linear model is right, it will predict _at least_ as well as a non-parametric model
	+ If the linear model is wrong, eventually the non-parametric model will beat it
- Need the sampling distribution of $D$ under the null hypothesis
    + Need to simulate the linear model
	+ And need to calculate $D$ over and over


## Simulate the parametric model

Simulate from the parametric model, assuming Gaussian noise
```{r}
sim.lm <- function(linfit, test.x) {
  n <- length(test.x)
  sim.frame <- data.frame(x=test.x)
  sigma <- sqrt(mse.residuals(linfit))  # There are other ways to get sigma
  y.sim <- predict(linfit,newdata=sim.frame)
  y.sim <- y.sim + rnorm(n,0,sigma)
  sim.frame <- data.frame(sim.frame,y=y.sim)
  return(sim.frame)
}
```

- OFFLINE: Write comments describing input and output of the `sim.lm` function
- OFFLINE: How can we check that `sim.lm` is working?
- OFFLINE: How can we switch to non-Gaussian noise?

## Calculate the test statistic

Calculate difference in MSEs between parametric and nonparametric models on a
data frame:
```{r}
calc.D <- function(df) {
  MSE.p <- mse.residuals(lm(y~x,data=df))
  MSE.np <- mse.residuals(gam(y~s(x),data=df))
  return(MSE.p - MSE.np)
}
```

- OFFLINE: Comments once again (with feeling this time)
- OFFLINE: Check that this function works properly (how?)


## Put the pieces together

Calculate the MSE difference on one simulation run:
```{r}
calc.D(sim.lm(glinfit,x))
```
Calculate the MSE difference on 200 simulation runs, so we get a sample
from the null hypothesis:
```{r}
null.samples.D <- replicate(200,calc.D(sim.lm(glinfit,x)))
```
How often does the simulation produce gaps bigger than what we really saw?
```{r}
mean(null.samples.D >= D.hat)
```

## Look at the distribution

Plot histogram of the sampling distribution, and the observed value:
```{r}
hist(null.samples.D,n=31,xlim=c(min(null.samples.D),1.2*D.hat),probability=TRUE,
     main="Distribution of D under the null model",
     xlab="D")
abline(v=D.hat, col="red")
```


## Second case: linear model is properly specified

Deliberately  left uncommented so that you can explore

```{r}
y2 <- 0.2+0.5*x + rnorm(length(x),0,0.15)
  # Why don't I have to make up x values again?
y2.frame <- data.frame(x=x,y=y2)
plot(x,y2,xlab="x",ylab="y",pch=16)
abline(0.2,0.5,col="grey",lwd=4)
```

## 

```{r}
y2.fit <- lm(y~x,data=y2.frame)
null.samples.D.y2 <- replicate(200,calc.D(sim.lm(y2.fit,x)))
(D.hat2 <- calc.D(y2.frame))
hist(null.samples.D.y2,n=31,probability=TRUE, main="Distribution of D when the linear model is right", xlab="D")
abline(v=D.hat2, col="red")
mean(null.samples.D.y2 >= D.hat2)
```

## Extensions/lessons

- Don't need to guess at specific nonlinearities; let the smoother figure it out
- Extends to testing a linear model with multiple predictors easily
    + Choice of whether the alternative then is an additive model vs. a general nonparametric regression
	+ Additive: more power to detect additive nonlinearities, can easily miss nonlinear interactions
	+ General model: sensitive to more mis-specifications but less power over-all (curse of dimensionality again)
- Nothing special about linear models --- the same idea works for testing any parametric regression model

## Another approach, based on residuals

\[
Y = \mu(X) + \epsilon, ~ \Expect{\epsilon|X} = 0
\]

- $\Rightarrow$ if we've got the right model, $\Expect{Y-\mu(X)|X=x} = 0$ for all $x$
- $\Rightarrow$ smooth the residuals and see if we get (close to) 0


```{r}
par(mfrow=c(1,2))
plot(gam(residuals(glinfit) ~ s(x)),residuals=TRUE,
  lwd=4,ylab="Residuals of the linear model", main="Reality nonlinear")
abline(h=0,col="grey")
plot(gam(residuals(y2.fit) ~ s(x)), residuals=TRUE,
     lwd=4, ylab="Residuals of the linear model", main="Reality linear")
abline(h=0,col="grey")
```

## Smoothing residuals, cont'd

- We'd need to calculate how far the smoothed-residuals curve is from 0
    + Once common choice is the $L_2$ norm, $\|r\|^2_2 \equiv \int{r^2(x) dx}$, or $\int{r^2(x) f(x) dx}$ to weight by the pdf $f(x)$
	+ Could approximate the weighted version by $\frac{1}{n}\sum_{i=1}^{n}{r^2(x_i)}$ (why?)
- Then simulate to get a null distribution for that test statistic






## Backup: Interactions in `gam()`

- Can extend the additive model to include interaction terms:
\[
Y = \alpha + \sum_{j=1}^{p}{f_j(x_j)} + \sum_{j=1}^{p}{\sum_{k=j+1}^{p}{f_{jk}(x_j, x_k)}} + \epsilon
\]
- Estimate $f_{jk}$ by taking partial residuals and doing a 2D smooth against $x_j$ _and_ $x_k$
- How do we do this in R?
    + `s(xj, xk)` fits a 2D ("thin plate") spline, good if `xj` and `xk` are on the same scale
	+ `te(xj, xk)` fits a "tensor-product" spline, good if `xj` and `xk` are measured on _different_ scales

- When do you want interactions?
    + When it makes sense (location interacts latitude and longitude)
	+ Diagonstic: scatter-plot residuals vs. `xj` and `xk`, look for patterns

- Reading: textbook, section 8.5


## Convergence when we add interactions to a GAM

- Convergence rate for a $k$-dimensional smoother is $O(n^{-4/(4+k)})$
- MSE of true regression function $\mu$ is $\sigma^2$
- For the additive model (w/ no interactions),
\[
MSE = \sigma^2 + (\text{approximation bias of best additive model})^2 + O(pn^{-4/5})
\]
- For the model with _all_ two-way interactions,
\[
MSE = \sigma^2 + (\text{approximation bias of best two-way model})^2 + O(p^2 n^{-4/6})
\]
    + There are ${p \choose 2} = O(p^2)$ possible interactions, each converging at rate $O(n^{-4/6})$, and these dominate the $O(n^{-4/5})$ of the additive terms
- For the model with all $k$-way interactions,
\[
MSE = \sigma^2 + (\text{approximation bias})^2 + O\left({p \choose k} n^{-4/(4+k)}\right)
\]


