--- title: Additive Models author: 36-402, Section A, Spring 2019 date: 19 February 2019 output: slidy_presentation --- \[ \newcommand{\TrueRegFunc}{\mu} \newcommand{\EstRegFunc}{\widehat{\mu}} \newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]} \newcommand{\Var}[1]{\mathbb{V}\left[ #1 \right]} \DeclareMathOperator*{\argmin}{argmin} % thanks, wikipedia! \] ```{r, include=FALSE} 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) ``` ## Previously: {.smaller} - We want to predict $Y$ from $X$ with small MSE, $\Expect{(Y-m(X))^2}$ - Optimum would be $\TrueRegFunc(x) = \Expect{Y|X=x}$; then $Y=\TrueRegFunc(X)+\epsilon$ with $\Expect{\epsilon|X=x}=0$ for all $x$ + $\sigma^2 \equiv \Expect{(Y-\mu(X))^2} > 0$ (generally) - We try to learn/estimate $\TrueRegFunc$ from data $(x_1, y_1), \ldots (x_n, y_n)$, getting an approximation $\EstRegFunc$ - Linear smoothers: \[ \EstRegFunc(x_0) = \sum_{i=1}^{n}{y_i w(x_0, x_i)} \] ## Previously (cont'd): {.smaller} - Kernel smoothing, 1D: \[ \EstRegFunc(x_0) = \sum_{i=1}^{n}{y_i \frac{K(\frac{x_0-x_i}{h})}{\sum_{k=1}^{n}{K(\frac{x_0-x_k}{h})}}} \] - Kernel smoothing, $p$-dimensional: \[ \EstRegFunc(x_0) = \sum_{i=1}^{n}{y_i \frac{\prod_{j=1}^{p}{K(\frac{x_{0j}-x_{ij}}{h_j})}}{\sum_{k=1}^{n}{\prod_{j=1}^{p}{K(\frac{x_{0j}-x_{kj}}{h_j})}}}} \] - Bias of kernel smoothing is $O(h^2)$ in any dimension - Variance of kernel smoothing is $O(n^{-1} h^{-p})$ in $p$ dimensions ## EXERCISE - Bias is $O(h^2)$, variance is $O(n^{-1} h^{-p})$ - MSE of kernel smoothing is \[ MSE(h) = \sigma^2 + b h^4 + v n^{-1}h^{-p} \] + (plus negligible remainders) 1. What's the best bandwidth? How does it change with $n$ and $p$? 2. Does the bias at the best bandwidth $\rightarrow 0$ as $n\rightarrow\infty$? 3. Does the variance at the best bandwidth $\rightarrow 0$ as $n\rightarrow\infty$? 4. Does $\EstRegFunc \rightarrow \TrueRegFunc$ as $n\rightarrow\infty$? ## SOLUTION #### Part 1: \begin{eqnarray} \left.\frac{dMSE}{dh}\right|_{h=h_{opt}} & = & 0 = 0 + 3bh_{opt}^3 - \frac{vp}{n}h_{opt}^{-p-1}\\ h_{opt}^{4+p} & = & \frac{vp}{3b}n^{-1}\\ h_{opt} & = & \left(\frac{vp}{3b}\right)^{1/(4+p)} n^{-\frac{1}{4+p}}\\ h_{opt} & = & O(n^{-\frac{1}{4+p}}) \end{eqnarray} #### Part 2: \begin{eqnarray} \mathrm{bias} & = & O(h_{opt}^2)\\ & = & O(n^{-\frac{2}{4+p}}) \rightarrow 0 \end{eqnarray} #### Part 3: \begin{eqnarray} \mathrm{variance} & = & O(n^{-1} h_{opt}^{-p})\\ & = & O(n^{-1} n^{-\frac{p}{4+p}})\\ & = & O(n^{-\frac{4}{4+p}}) \rightarrow 0 \end{eqnarray} #### Part 4: Yes, if bias and variance both $\rightarrow 0$, then the estimate must converge on the truth ## Kernel regression is universally consistent - In any dimension, as $n\rightarrow \infty$, the bias of $\EstRegFunc$ $\rightarrow 0$ and the variance $\rightarrow 0$, so \[ \EstRegFunc \rightarrow \TrueRegFunc \] - Kernel regression is **universally consistent** + At least if $\TrueRegFunc$ is smooth enough for a 2nd-order Taylor series ## Kernel regression can converge very slowly > - The MSE at the best bandwidth is \[ MSE = \sigma^2+O(h_{opt}^4) = \sigma^2 + O(n^{-\frac{4}{4+p}}) \] > + So MSE approaches that of $\TrueRegFunc$ > - BUT we need exponentially more data as the dimension grows > + With $n={10}^{12}$ and $p=1$, $n^{-4/5} = `r signif((1e12)^(-4/5),2)`$ > + With $n={10}^{12}$ and $p=100$, $n^{-4/(104)} = `r signif((1e12)^(-1/26), 2)`$ ## Kernel regression can converge very slowly ```{r, echo=FALSE} curve(x^(-1),from=1,to=1e15,log="x",xlab="n",ylab="Excess MSE") curve(x^(-4/5),add=TRUE,lty="dashed") curve(x^(-1/26),add=TRUE,lty="dotted") legend("topright",legend=c(expression(n^{-1}), expression(n^{-4/5}),expression(n^{-1/26})), lty=c("solid","dashed","dotted")) ``` ## The problem isn't kernel regression - $k$-nearest neighbors is also universally consistent, and also has the same rate of convergence + So do splines, trees, random forests, deep neural networks, ... - In fact, this is the best possible rate for _any_ universally-consistent method + You can converge faster for _some_ $\TrueRegFunc$, if you don't converge at all for others - This --- exponentially slow convergence as $p$ grows --- is the **curse of dimensionality** ## Lifting the curse > - Impose special structure on $\TrueRegFunc$, either because: > + We believe it's true, or > + We're willing to make that approximation > - Extreme example: Only use linear functions > - If $\TrueRegFunc$ is linear, then $MSE = \sigma^2 + O(p/n)$ > - If $\TrueRegFunc$ isn't linear, then $MSE = \sigma^2 + (\mathrm{bias\ of\ best\ linear\ model})^2 + O(p/n) = \sigma^2 + O(1)$ > - Can we do something less extreme than linear models? ## Additive models {.smaller} - Linear model: \[ Y = \beta_0 + \sum_{j=1}^{p}{\beta_j x_j} + \epsilon \] - **Additive** model: \[ Y = \alpha + \sum_{j=1}^{p}{f_j(x_j)} + \epsilon \] + Every linear model is additive, $f_j(x_j) = \beta_j x_j$ + Not every additive model is linear, e.g, $x_1 + \sin{(x_2)} + \log{(1+x_3)}$ - The $f_j$ are called **partial response** functions + To get total response (= prediction for $Y$), add up all the partial responses, plus $\alpha$ + All else being equal, if you change $x_j$ to $x_j^{\prime}$, the predicted $Y$ changes by $f_j(x_j^{\prime}) - f_j(x_j)$ + Like transforming your predictors in a linear model ## So how do we estimate an additive model? - Re-arrange terms \[ f_j(x_j) + \epsilon = Y - (\alpha + \sum_{k \neq j}{f_k(x_k)}) \] - If we knew all the other functions, we could estimate $f_j$ by 1D smoothing #### **Backfitting** (or **Gauss-Seidel**) 0. Start with $\alpha = \overline{y}$, and $f_1 = f_2 \ldots = f_p = 0$ 1. For $j \in 1:p$: a. $r_i \leftarrow y_i - (\alpha + \sum_{k\neq j}{f_k(x_{ik})})$ (**partial residual**) b. Use your favorite smoother to regress $r$'s on $x_j$ (1D smoothing), get function $s_j$ c. $f_j \leftarrow s_j - \overline{s_j}$ (keep functions centered at 0) 2. Repeat (1) until things stop changing ## In R: `mgcv` {.smaller} - Several packages available, older `gam` and newer, nicer `mgcv` - Main function is called `gam` (Generalized Additive Model) + We'll get to the "generalized" in a few weeks - Syntax is very like `lm` _but_ can tell it to smooth variables: ```{r, warning=FALSE, message=FALSE} library(mgcv) x <- matrix(runif(1000, min=-pi, max=pi), ncol=2) y <- 0.5*x[,1]+sin(x[,2])+rnorm(n=nrow(x), sd=0.1) df <- data.frame(y=y, x1=x[,1], x2=x[,2]) babys.first.am <- gam(y ~ s(x1)+s(x2), data=df) par(mfrow=c(1,2)); plot(babys.first.am) ``` ## `s()` {.smaller} - People often read `s()` as "smooth", but it really means "spline" - A **smoothing spline** is a curve balancing fit to data against wiggliness: \[ s_{\lambda} = \argmin_{m}{\frac{1}{n}\sum_{i=1}^{n}{(y_i - m(x_i))^2} + \lambda\int{\left(\frac{d^2 m}{dx^2}(x)\right)^2 dx}} \] + Solution for each $\lambda$ is a smooth, piecewise-cubic polynomial + $s_{\lambda}$ is linear in $\mathbf{y}$ (surprisingly) + Polynomials are _fast_ to calculate and estimate + $\lambda\rightarrow 0$ gives interpolation; $\lambda\rightarrow\infty$ gives a linear model + (There are also other kinds of spline) - `gam()` knows about _lots_ of other smoother techniques + `te()`, `ti()`, etc., etc., etc. + Mostly about how to handle interactions + See textbook ## Methods for `mgcv::gam` - `predict`: Works just like with `lm` (or `npreg` for that matter), see `help(predict.gam)` for details - `fitted`, `residuals`: Work just like with `lm` - `plot`: We just saw this; many options, some described in text, see `help(plot.gam)` for details - `coefficients`: _Mostly_ coefficients of the spline polynomials --- NOT INTERPRETABLE - Can include linear terms, e.g., `gam(y ~ x1+s(x2))` - Can include interaction terms, e.g., `gam(y ~ s(x1, x2))` ## A small real-data example ```{r} library(gamair) data(chicago) chicago.am <- gam(death ~ s(tmpd) + s(pm10median), data=chicago) ``` ## A small real-data example (cont'd) ```{r} par(mfrow=c(1,2)); plot(chicago.am) ``` ## A small real-data example (cont'd) ```{r} # Handling missing values by padding in NAs: am.fits <- predict(chicago.am, newdata=chicago, na.action=na.pass) plot(death ~ time, data=chicago) points(chicago$time, am.fits, pch="x", col="red") ``` ## But what about the curse? - If $\TrueRegFunc$ is additive, then \[ MSE = \sigma^2 + O(p n^{-4/5}) \] which is pretty good - If $\TrueRegFunc$ is not additive, then \[ MSE = \sigma^2 + (\mathrm{bias\ of\ best\ additive\ model})^2 + O(pn^{-4/5}) \] + This is still better, as $n\rightarrow\infty$, than linear models ## Summing up - Estimating regression functions in high dimensions is **intrinsically hard** - We need either a _lot_ of data OR special assumptions about the function - Additive models are an easy step beyond linear models with a lot of advantages - Get used to typing `gam` rather than `lm`