GLM Review (1)

class: center, middle, inverse, title-slide

.title[
# GLM Review (1)
]
.author[
### Yue Jiang
]
.date[
### STA 440 / Spring 2025
]

---

### Review: GLMs

A broad and flexible class of models consisting of three components:

1. An assumed distribution of the outcome `$Y$` (which follows an .vocab[exponential family])
2. The linear predictor `$\mathbf{X}\boldsymbol\beta$` consisting of observed covariates and unknown parameters of interest
3. An invertible .vocab[link function] `$g$` relating the conditional expectation of `$Y$` with the
linear predictor

`$$E(Y | \mathbf{X}) = g^{-1}\left(\mathbf{X}\boldsymbol\beta\right)$$`

---

### Review: GLMs

Many common distributions are members of the exponential family, for instance
normal, exponential, Poisson, Bernoulli, etc., and are thus suitable
choices for the assumed outcome variable.

`$$f(Y | \theta, \phi) = \exp\left\{\frac{y\theta - b(\theta)}{a(\phi)} \right\}c(y, \phi),$$`
(remember that in a GLM we assume that `$\theta$` depends on the linear predictor `$\mathbf{X}\boldsymbol\beta$` "in some way").

---

### Review: GLMs

Exponential families are nice!

They're well-behaved, especially for Bayesian applications, estimation theory,
and statistical inference (e.g., in one-parameter exponential families and UMP
testing).

Exponential family distributions are guaranteed to have conjugate priors, and so 
their posterior predictive distributions can also be explicitly written.

In GLMs, we always assume that there is some functional relationship between
the conditional expectation of the outcome and the .vocab[linear predictor] -
some linear combination of the unknown parameters and the observed covariates.

This is why they are called generalized **linear** models.

---

### Review: GLMs

Finally, this functional relationship relating the linear predictor and the
conditional expectation of `$Y$` is known as the .vocab[link function]. Often
times we take this to be the .vocab[canonical link] function (i.e., when we
take `$\theta$` to be `$\mathbf{X}\boldsymbol\beta$`).

Use of the canonical link function has some nice properties (e.g., guaranteeing
existence of MLEs due to concavity of log-likelihood over a convex set, or
guaranteeing that the observed and Fisher information are the same, etc., which
has really nice properties for fitting GLMs as we'll see later today).

---

### Review: GLMs

.question[
Why would you ever use a non-canonical link function?
]

Sometimes, we might have *good reason* to use a non-canonical link depending on
our assumptions of the data or the particular interpretations we want our
regression estimates to have.

---

### Binary regression

For binary outcomes, we generally choose a .vocab[Bernoulli] distribution for
`$Y$`, such that `$E(Y)$` can be interpreted as the probability of success. As 
you've seen before, the *logit* link is the canonical link for this 
distribution.

What are some other potential choices of link function?

---

### Binary regression

Since `$E(Y)$` lives in `$[0, 1]$`, any function with that same range might do. In
particular, inverses of any continuous CDF. Suppose

`\begin{align*}
Y^\star &= \mathbf{X}\boldsymbol\beta + \epsilon^\star
\end{align*}`

---

### Binary regression

`\begin{align*}
P(Y = 1 | \mathbf{X}) &= P(Y^\star > 0 | \mathbf{X})\\
&= P(\mathbf{X}\boldsymbol\beta + \epsilon^\star > 0 | \mathbf{X})\\
&= P(\epsilon^\star > -\mathbf{X}\boldsymbol\beta | \mathbf{X})\\
&= 1 - F(-\mathbf{X}\boldsymbol\beta)\\
&= F(\mathbf{X}\boldsymbol\beta)
\end{align*}`

(last equality if distribution is symmetric).

.question[
What does this mean?
]

---

### Binary regression

Logistic regression assumes `$$F(\mathbf{X}\boldsymbol\beta) = (1 + \exp(-\mathbf{X}\boldsymbol\beta))^{-1}$$` implying `$$\mathbf{X}\boldsymbol\beta = \log\left(\frac{P(Y = 1 | \mathbf{X})}{1 - P(Y = 1|\mathbf{X})} \right)$$`.

Probit regression assumes `$$F(\mathbf{X}\boldsymbol\beta) = \Phi(\mathbf{X}\boldsymbol\beta)$$` implying `$$\mathbf{X}\boldsymbol\beta = \Phi^{-1}(P(Y = 1 | \mathbf{X}))$$`

---

### Binary regression

Coefficients in such binary regression models are equivalent to covariate effects on a latent variable that following a linear model with errors with distribution function `$F$`.

.question[
How does this relate to our previous discussion of model selection?
]

---

### Estimating bike crashes in NC counties

---

### Data

```
## # A tibble: 6 × 6
##   county       pop med_hh_income traffic_vol pct_rural crashes
##   <chr>      <dbl>         <dbl>       <dbl>     <dbl>   <dbl>
## 1 Alamance  166436          50.5         182        29      77
## 2 Alexander  37353          49.1          13        73       1
## 3 Alleghany  11161          39.7          28       100       1
## 4 Anson      24877          38            79        79       7
## 5 Ashe       27109          41.9          18        85       4
## 6 Avery      17505          41.7          35        89       5
```
Suppose we thought these crashes came from a Poisson distribution, conditional on the covariates (related through a linear predictor). 
.question[
How might we estimate the parameter of that Poisson distribution, given our
observed data? How might we link it to regression coefficients?
]

---

### Poisson regression

`\begin{align*}
\log\large(\underbrace{E(Y | \mathbf{X}}_{\lambda})\large) = \beta_0 + \mathbf{X}^T\boldsymbol\beta
\end{align*}`

.vocab[Generalized linear model] often used for count (or rate) data
- Assumes outcome has Poisson distribution
- Canonical link: log of conditional expectation of response has linear 
relationship with predictors

.question[
Can we differentiate the (log) likelihood function, set it equal to zero, and
solve for the MLEs for `$\boldsymbol\beta = (\beta_0, \beta_1, \cdots, \beta_p)$`
as before?
]

---

### Poisson regression

`\begin{align*}
\log \mathcal{L}&= \sum_{i = 1}^n \left( y_i\log \lambda - \lambda - \log y_i! \right)\\
&= \sum_{i = 1}^n y_i\mathbf{X}_i\boldsymbol\beta - e^{\mathbf{X}_i\boldsymbol\beta} - \log y_i!
\end{align*}`

We would like to solve the equations

`\begin{align*}
\left(\frac{\partial \log \mathcal{L}}{\partial \beta_j}\right) \stackrel{set}{=} \mathbf{0},
\end{align*}`

but there is no closed-form solution, as this is a transcendental equation in 
the parameters of interest.

.question[
How might we solve these equations numerically?
]

---

### Newton's method

.vocab[Score vector] and .vocab[Hessian] for `$\log \mathcal{L}(\boldsymbol\theta | \mathbf{X})$` with `$\boldsymbol\theta = (\theta_1, \cdots, \theta_p)^T$`:

`$$\nabla \log \mathcal{L} = \begin{pmatrix}
\frac{\partial \log \mathcal{L}}{\partial \boldsymbol\theta_1}\\
\vdots\\
\frac{\partial \log \mathcal{L}}{\partial \boldsymbol\theta_p}
\end{pmatrix}$$`
`$$\nabla^2 \log \mathcal{L} = \begin{pmatrix}
\frac{\partial^2 \log\mathcal{L}}{\partial \theta_1^2} & \frac{\partial^2 \log\mathcal{L}}{\partial \theta_1 \theta_2} & \cdots & \frac{\partial^2 \log\mathcal{L}}{\partial \theta_1\theta_p}\\
\frac{\partial^2 \log\mathcal{L}}{\partial \theta_2\theta_1} & \frac{\partial^2 \log\mathcal{L}}{\partial \theta_2^2} & \cdots & \frac{\partial^2 \log\mathcal{L}}{\partial \theta_2\theta_p} \\
\vdots & \vdots & \ddots & \vdots \\
\frac{\partial^2 \log\mathcal{L}}{\partial \theta_p\theta_1} & \frac{\partial^2 \log\mathcal{L}}{\partial \theta_p\theta_2} & \cdots & \frac{\partial^2 \log\mathcal{L}}{\partial \theta_p^2}
\end{pmatrix}$$`
 
---

### Newton's method

We can modify the Newton-Raphson algorithm for higher dimensions:

- Start with initial guess `$\boldsymbol\theta ^{(0)}$`
- Iterate `$\boldsymbol\theta^{(t + 1)} = \boldsymbol\theta^{(t)} - \left(\nabla^2\log\mathcal{L}(\boldsymbol\theta^{(t)} | \mathbf{X}) \right)^{-1} \left( \nabla \log\mathcal{L}(\boldsymbol\theta^{(t)} | \mathbf{X}) \right)$`
- Stop when convergence criterion is satisfied

Under certain conditions, a global maximum exists; this again is guaranteed for 
many common applications.

Computing the Hessian can be computationally demanding (and annoying), but there 
are ways around it in practice.

---

### Poisson regression

`\begin{align*}
\log \mathcal{L}&= \sum_{i = 1}^n y_i\mathbf{X}_i\boldsymbol\beta - e^{\mathbf{X}_i\boldsymbol\beta} - \log y_i!
\end{align*}`

.question[
What are the score vector and Hessian corresponding to the Poisson regression
log-likelihood? What would the Newton-Raphson update steps be?
]

---

### Poisson regression

`\begin{align*}
\log \mathcal{L}&= \sum_{i = 1}^n y_i\mathbf{X}_i\boldsymbol\beta - e^{\mathbf{X}_i\boldsymbol\beta} - \log y_i!\\
\nabla \log \mathcal{L}&= \sum_{i = 1}^n \left(y_i - e^{\mathbf{X}_i\boldsymbol\beta}\right)\mathbf{X}_i^T\\
\nabla^2 \log \mathcal{L} &= -\sum_{i = 1}^n e^{\mathbf{X}_i\boldsymbol\beta}\mathbf{X}_i\mathbf{X}_i^T
\end{align*}`

Newton-Raphson update steps for Poisson regression:

`\begin{align*}
\boldsymbol\beta^{(t+1)} &= \boldsymbol\beta^{(t)} - \left(-\sum_{i = 1}^n e^{\mathbf{X}_i\boldsymbol\beta^{(t)}}\mathbf{X}_i\mathbf{X}_i^T \right)^{-1}\left(\sum_{i = 1}^n \left(y_i - e^{\mathbf{X}_i\boldsymbol\beta^{(t)}}\right)\mathbf{X}_i^T \right)
\end{align*}`

---

### Some discussion questions

.question[
- Do you think Poisson regression was a good idea in this case? Why or why not? 
- What about a non-canonical link? Suppose you used an identity link instead. What might the likelihood equations look like, and how would you solve them?
- How about the logistic regression example from earlier? How about likelihood equations for non-canonical link functions in that setting?
]

---

### Another example...

---

### Another example...