---
title: Logistic Regression to Generalized Additive Models
author: 36-402, Sec. A, Spring 2019
date: 7 March 2019
output: slidy_presentation
---

```{r, echo=FALSE}
library(knitr)
opts_chunk$set(size="small",background="white", cache=TRUE,
               autodep=TRUE, warning=FALSE, error=FALSE,
               tidy=TRUE, tidy.opts=list(comments=FALSE))
```

\[
\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\newcommand{\Prob}[1]{\mathbb{P}\left[ #1 \right]}
\newcommand{\Var}[1]{\mathbb{V}\left[ #1 \right]}
\]

## A cat picture

```{r, fig.retina=NULL, out.width=800, echo=FALSE}
knitr::include_graphics("2019-03-07-cat-1.jpg")
```

## The logistic curve {.smaller}

\[
\frac{e^u}{1+e^u}  = \frac{1}{1+e^{-u}}
\]

```{r, echo=FALSE}
library(faraway) # for logit and ilogit functions
curve(ilogit(x), from=-10, to=10, xlab="log odds ratio", ylab="probability")
```


## The logistic regression model {.smaller}

\begin{eqnarray}
p(x) & = & \frac{e^{\beta_0 + \beta\cdot x}}{1+e^{\beta_0 + \beta\cdot x}}\\
g(p) & \equiv & \log{\frac{p}{1-p}} = \beta_0 + \beta\cdot x
\end{eqnarray}

```{r, echo=FALSE}
curve(ilogit(x), from=-10, to=10, xlab=expression(beta[0]+beta*x), ylab="probability")
```

## {.smaller}

```{r, echo=FALSE}
curve(ilogit(x), from=-10, to=10, xlab=expression(x), ylab="probability")
curve(ilogit(x+0.5), add=TRUE, col="red")
curve(ilogit(2*x), add=TRUE, col="green")
curve(ilogit(2*x+0.5), add=TRUE, col="blue")
legend("topleft", legend=c(expression(paste(beta[0]==0,",", beta==1)),
                           expression(paste(beta[0]==.5,",", beta==1)),
                           expression(paste(beta[0]==0,",", beta==2)),
                           expression(paste(beta[0]==.5,",", beta==2))),
       col=c("black", "red", "green", "blue"),
       lty="solid")
```

## Why the log-odds ratio, of all things?

- Start with the likelihood = probability of $Y_i$ given $X_i$, as a function of parameters:
\begin{eqnarray}
\Prob{Y_1=y_1, \ldots Y_n=y_n|X_1=x_1, \ldots X_n=x_n} & = & \prod_{i=1}^{n}{\Prob{Y_i=y_i|X_i=x_i}}\\
& = & \prod_{i=1}^{n}{p(x_i)^{y_i} (1-p(x_i))^{1-y_i}}
\end{eqnarray}
- Take the log for the usual reasons:
\begin{eqnarray}
\ell & = & \sum_{i=1}^{n}{y_i\log{(p(x_i))} + (1-y_i)\log{(1-p(x_i))}}\\
& = & \sum_{i=1}^{n}{\log{(1-p(x_i))} + y_i\left(\log{(p(x_i))} - \log{(1-p(x_i))}\right)}\\
& = & \sum_{i=1}^{n}{\log{(1-p(x_i))} + y_i\log{\left(\frac{p(x_i)}{1-p(x_i)}\right)}}
\end{eqnarray}
- $\ell$ depends on $y_i$ only through $\log{\left(\frac{p(x_i)}{1-p(x_i)}\right)} =$ log-odds ratio
    + log-odds ratio is the **natural parameter**
- Logistic regression = make the natural parameter linear in $x$:
\[
\ell(\beta_0, \beta) =  \sum_{i=1}^{n}{-\log{\left(1+e^{\beta_0 + x_i \cdot \beta}\right)}} + \sum_{i=1}^{n}{y_i (\beta_0 + x_i \cdot \beta)}
\]

## Interpretation

- $\beta_j$ is how much a one-unit difference in $x_j$ changes the predicted
_odds ratio_ for $Y$
    + **not** the predicted change in $Y$
	+ Predicted change in $Y$ depends on the start point


```{r, echo=FALSE}
curve(ilogit(x), from=-10, to=10, xlab="log odds ratio", ylab="probability")
points(x=-2, y=ilogit(-2), col="red", pch=15)
points(x=-1, y=ilogit(-1), col="blue", pch=15)
points(x=0.1, y=ilogit(0.1), col="red", pch=16)
points(x=1.1, y=ilogit(1.1), col="blue", pch=16)
points(x=2, y=ilogit(2), col="red", pch=17)
points(x=3, y=ilogit(3), col="blue", pch=17)
```

## Example {.smaller}

```{r}
ch <- read.csv("http://www.stat.cmu.edu/~cshalizi/uADA/19/exams/1/ch.csv")
ch <- ch[,-1] # First column is just an index
ch <- na.omit(ch) # Not rec'd for exam but simplifies a few steps
ch.logistic <- glm(start ~ exports+fractionalization*dominance, data=ch,
                   family="binomial")
coefficients(ch.logistic)
```

## Have another cat picture


```{r, fig.retina=NULL, out.width=400, echo=FALSE}
knitr::include_graphics("2019-03-07-cat-3.jpg")
```

## Residuals

- "Response" residuals $\widehat{\epsilon}_{i,\text{response}} \equiv y_i - p(x_i)$
    + Expectation $=0$ if the model is right
	+ Variance isn't constant even if the model is right
- "Pearson" residuals $\widehat{\epsilon}_{i,\text{Pearson}} \equiv \frac{y_i - p(x_i)}{\sqrt{p(x_i)(1-p(x_i))}}$
    + Expectation $=0$ if the model is right
	+ Variance $=1$ if the model is right
- "Deviance" residuals $\widehat{\epsilon}_{i,\text{deviance}} \equiv \sqrt{2\left(y_i\log{p(x_i)} + (1-y_i)\log{(1-p(x_i))}\right)} \mathrm{sgn}(y_i-p(x_i))$
    + Squared deviance residuals add up to twice the negative log-likelihood = "deviance"

## Residuals: Response

```{r}
plot(ch$exports, residuals(ch.logistic, type="response"),
     xlab="Exports", ylab="Residuals", main="Response residuals")
abline(h=0, col="grey")
lines(smooth.spline(x=ch$exports, y=residuals(ch.logistic, type="response")))
```

## Residuals: Pearson

```{r}
plot(ch$exports, residuals(ch.logistic, type="pearson"),
     xlab="Exports", ylab="Residuals", main="Pearson residuals")
abline(h=0, col="grey")
lines(smooth.spline(x=ch$exports, y=residuals(ch.logistic, type="pearson")))
```

## Squared Residuals: Pearson

```{r}
plot(ch$exports, residuals(ch.logistic, type="pearson")^2,
     xlab="Exports", ylab="Squared residuals", main="Squared Pearson residuals")
abline(h=1, col="grey")
lines(smooth.spline(x=ch$exports, y=residuals(ch.logistic, type="pearson")^2))
```

## Classification

- We classify each point as "should be $Y=1$" or "should be $Y=0$"
```{r}
mean(ifelse(fitted(ch.logistic)<0.5,0,1) != ch$start)
```

Is this good or bad?

```{r}
mean(ch$start)
```

EXERCISE: Why does the model have the _same_ error rate as the constant?


## Calibration

- Logistic regression predicts **probabilities**
- Minimal test: When the model says $\Prob{Y=1} = 0.5$, $Y$ should be 1 half the time
- Probabilities are **calibrated** when events with predicted probability $p$ happen with _frequency_ $p$
- If each event has its own probability, group events with _close_ probabilities


## Checking calibration {.smaller}


```{r}
frequency.vs.probability <- function(p.lower,p.upper=p.lower+0.005,
  model, events) {
  fitted.probs <- fitted(model)
  indices <- (fitted.probs >= p.lower) & (fitted.probs < p.upper)
  matching.probs <- fitted.probs[indices]
  ave.prob <- mean(matching.probs)
  frequency <- mean(events[indices])
  # "Law of total variance": Var[Y]=E[Var[Y|X]] + Var[E[Y|X]]
  total.var <- mean(matching.probs*(1-matching.probs)) + var(matching.probs)
  se <- sqrt(total.var/sum(indices))
  return(c(frequency=frequency,ave.prob=ave.prob,se=se))
}
```

(Can you add comments?)

## {.smaller}

Now apply the function a bunch (why these numbers?)
```{r}
f.vs.p <- sapply(seq(from=0.04, to=0.12, by=0.005), frequency.vs.probability,
                 model=ch.logistic, events=ch$start)
```

This is "turned on its side" from a data frame, so let's fix.  (Can you find a more elegant way to do this?)
```{r}
f.vs.p <- data.frame(frequency=f.vs.p["frequency",],
                      ave.prob=f.vs.p["ave.prob",],
                      se=f.vs.p["se",])
```

## {.smaller}

```{r}
plot(frequency~ave.prob,data=f.vs.p,xlim=c(0,1),ylim=c(0,1),
  xlab="Predicted probabilities",ylab="Observed frequencies")
rug(fitted(ch.logistic),col="grey")
abline(0,1,col="grey")
segments(x0=f.vs.p$ave.prob,y0=f.vs.p$ave.prob-1.96*f.vs.p$se,
  y1=f.vs.p$ave.prob+1.96*f.vs.p$se)
```

## Have another cat photo


```{r, fig.retina=NULL, out.width=800, echo=FALSE}
knitr::include_graphics("2019-03-07-cat-2.jpg")
```



## How do we maximize the log-likelihood?

- Take derivatives, set them to 0
\[
\frac{\partial \ell}{\partial \beta_j} =  \sum_{i=1}^{n}{\left(y_i - p(x_i;\beta_0,\beta)\right) x_{ij}}
\]
- _Like_ the "estimating" or "normal" equations for linear regression,
\[
\frac{\partial MSE}{\partial \beta_j} = 2\sum_{i=1}^{n}{\left(y_i - \beta_0 - x_i \cdot \beta\right) x_{ij}}
\]
   + _Those_ equations have a closed-form solution
   + The estimating equations for _logistic_ regression do not

## How do we maximize (cont'd)?

- Try to turn this into a least-squares problem
- $g(p)$ is linear in $x$ --- can we transform somehow?
    + Note: $y=0$ or $y=1$ so $g(y) = \pm \infty$
- Taylor approximation to the rescue:
\begin{eqnarray}
g(y) & \approx & g(p(x)) + (y-p(x))) g^{\prime}(p)\\
& = & (\beta_0 + \beta\cdot x) + (y-p(x)) g^{\prime}(p)\\
& \equiv & z
\end{eqnarray}
- $z = \text{(linear predictor)} + \text{(mean-0 noise)}$
- **Solution**: Linearly regress $z$ on $x$
    + $\Var{Z|X} = p(x)(1-p(x)) \left(g^{\prime}(p)\right)^2$ so use weighted least squares
	+ Weights, and $Z$'s, depend on current guess about $\beta$

## Iterative weighted least squares / Fisher scoring

0. Start with a guess about $\beta$
1. Calculate $p(x_i)$, $g^{\prime}(p(x_i))$ for each $i$
2. Create
\begin{eqnarray}
z_i & = & g(p(x_i)) + (y_i - p(x_i))g^{\prime}(p)\\
w_i & = & p(x_i)(1-p(x_i)) \left(g^{\prime}(p(x_i))\right)^2
\end{eqnarray}
3. Minimize over $(b_0, b)$ to get the new $(\beta_0, \beta)$:
\[
\sum_{i=1}^{n}{\frac{(z_i - (b_0 + b \cdot x_i))^2}{w_i}}
\]
4. Go back to (1) until the predictions $p(x_i)$ stop changing

## Why does IWLS work?

- It has a fixed point when we've got the right parameters
- It's _also_ Newton's method, applied to minimizing the negative log-likelihood
    + $\theta \leftarrow \theta - (\nabla \nabla \ell(\theta))^{-1} \nabla \ell(\theta)$
    + (Some details about exact vs. expected 2nd derivatives here)
- So iterating _this_ weighted least squares problem 
numerically maximizing the log-likelihood

## Generalized additive model

\[
g(p) = \alpha + \sum_{j=1}^{p}{f_j(x_j)}
\]

- $f_j$ are still partial response functions, but for log-odds of $Y=1$, not for $Y$

- We can use IWLS to estimate the $f_j$: just fit an additive model for
the $z_i$'s instead of a linear model

## Example

```{r}
library(mgcv)
ch.gam.1 <- gam(start ~ s(exports)+s(fractionalization) + s(fractionalization, by=dominance),
                   data=ch, family="binomial")
plot(ch.gam.1)
```

## A better model

```{r}
ch.gam.2 <- gam(start ~ s(peace)+s(lnpop), data=ch, family="binomial")
```

```{r, echo=FALSE}
par(mfrow=c(1,2)); plot(ch.gam.2); par(mfrow=c(1,1))
```

- "Better" by cross-validated log-likelihood

## Generalize!

\begin{eqnarray}
\epsilon(x)  & \equiv  & \Expect{Y|X=x}\\
\eta(x)  & \equiv & \beta_0 + x \cdot \beta ~ (\text{or an additive model or whatever})\\
\eta(x) & = & g(\epsilon(x))\\
Z & \equiv &  g(\epsilon(x)) + (Y-\epsilon(x)) g^{\prime}(\epsilon(x))\\
& = & \eta(x) + (Y-\epsilon(x)) g^{\prime}(\epsilon(x))\\
\Expect{Z|X=x} & = & \eta(x)\\
\Var{Z|X=x} & = & \left( g^{\prime}(\epsilon(x)) \right)^2 \Var{Y|X=x}
\end{eqnarray}

We can do IWLS to recover the $\eta$ function, **without** transforming
the response

## Where does $g()$ come from?

- There's some distribution for $Y|X$
    + E.g., binomial for a binary $Y$
	+ Or maybe Poisson for a count-valued $Y$
	+ Or a gamma for positive-valued continuous $Y$ or...
- Look for the **natural parameter** of that distribution
    + For binomial, $\log{p/(1-p)}$ (as we saw)
	+ Needs some math
	+ Look at how that natural parameter relates to the mean

## Example: Poisson regression

- $Y \sim \mathrm{Pois}(\lambda)$ iff
\[
\Prob{Y=y} = \frac{\lambda^y}{y!}e^{-\lambda}
\]
    + $\Expect{Y} = \lambda$, $\Var{Y} = \lambda$
	+ Natural model for counts ("law of rare events")
- Take log of probability
\[
y\log{\lambda} - \log{(y!)} - \lambda
\]
    + $\Rightarrow$ natural parameter is $\log{\lambda}$
	+ Set $g(m) = \log{m}$

## Example: Poisson regression

- Death in Chicago
```{r}
library(gamair); data(chicago)
chicago.gam <- gam(death ~ s(tmpd) + s(so2median) + s(o3median) + s(pm10median),
                   data=chicago,  family="poisson")
```

## Example: Poisson regression

```{r}
plot(chicago.gam)
```

## Example: Poisson regression

```{r}
plot(death ~ time, data=chicago)
lines(chicago$time, predict(chicago.gam, newdata=chicago, type="response",
                            na.action=na.pass),
      type="l", col="red")
```

## Example: Poisson regression

```{r}
plot(residuals(chicago.gam, type="pearson"))
```






