GLM Review (2)

class: center, middle, inverse, title-slide

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

---

### Bike crashes

---

### Bike crashes

```
## # 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 \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*}`

---

### Let's fit a model

``` r
m1 <- glm(crashes ~ traffic_vol + pct_rural, 
          data = bike, family = "poisson")

summary(m1)
```

```
## 
## Call:
## glm(formula = crashes ~ traffic_vol + pct_rural, family = "poisson", 
##     data = bike)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  5.9821805  0.0537489 111.299   <2e-16
## traffic_vol  0.0015406  0.0001663   9.263   <2e-16
## pct_rural   -0.0445581  0.0008751 -50.919   <2e-16
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 14931.0  on 99  degrees of freedom
## Residual deviance:  2360.7  on 97  degrees of freedom
## AIC: 2870.3
## 
## Number of Fisher Scoring iterations: 5
```

.question[
What do you think of this model?
]

---

### A few more models

`\begin{align*}
\log\left(E(y_i | \mathbf{x}_i)\right) &= \beta_0 + \beta_1(pop_i) + \beta_2(traffic_i) + \beta_3(rural_i)\\
\end{align*}`

``` r
m2 <- glm(crashes ~ traffic_vol + pct_rural + pop, 
          data = bike, family = "poisson")
```

`\begin{align*}
\log\left( \frac{E(y_i | \mathbf{x}_i)}{pop_i} \right) &= \beta_0 + \beta_1(traffic_i) + \beta_2(rural_i)
\end{align*}`

``` r
m3 <- glm(crashes ~ traffic_vol + pct_rural, offset = log(pop),
          data = bike, family = "poisson")
```

.question[
What are the differences in the two models above?
]

---

### Let's fit a model

``` r
mean(bike$crashes)
```

```
## [1] 74.46
```

``` r
var(bike$crashes)
```

```
## [1] 23032.11
```

---

### Let's fit a model

``` r
m4 <- glm(crashes ~ traffic_vol + pct_rural + I(pop/1000), 
          data = bike, family = "quasipoisson")

summary(m4)
```

```
## 
## Call:
## glm(formula = crashes ~ traffic_vol + pct_rural + I(pop/1000), 
##     family = "quasipoisson", data = bike)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  5.6557252  0.2272488  24.888  < 2e-16
## traffic_vol -0.0000930  0.0007429  -0.125    0.901
## pct_rural   -0.0377608  0.0036378 -10.380  < 2e-16
## I(pop/1000)  0.0012608  0.0001729   7.291  8.7e-11
## 
## (Dispersion parameter for quasipoisson family taken to be 17.17315)
## 
##     Null deviance: 14931.0  on 99  degrees of freedom
## Residual deviance:  1482.6  on 96  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
```

---

### Let's fit a model

``` r
library(MASS)
m5 <- glm.nb(crashes ~ traffic_vol + pct_rural + I(pop/1000), 
          data = bike)

summary(m5)
```

```
## 
## Call:
## glm.nb(formula = crashes ~ traffic_vol + pct_rural + I(pop/1000), 
##     data = bike, init.theta = 2.845949661, link = log)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  5.563e+00  3.090e-01  18.005  < 2e-16
## traffic_vol  6.209e-05  1.362e-03   0.046  0.96365
## pct_rural   -3.793e-02  3.660e-03 -10.362  < 2e-16
## I(pop/1000)  1.702e-03  5.208e-04   3.269  0.00108
## 
## (Dispersion parameter for Negative Binomial(2.8459) family taken to be 1)
## 
##     Null deviance: 580.33  on 99  degrees of freedom
## Residual deviance: 106.57  on 96  degrees of freedom
## AIC: 861.99
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.846 
##           Std. Err.:  0.458 
## 
##  2 x log-likelihood:  -851.990
```

---

### Some discussion questions

.question[
- We fit three models. What can you conclude from each? Which one, if any, of these models might you want to use?
- If you're not satisfied with these models for *any* reason, why is that the case?
]

---

### Another example...

---

### Another example...

---

### Some discussion questions

.question[
- What variables might you want to include?
- What *model* might make sense here?
]