--- 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")) ```