poisson_regression.Rmd

# 广义线性模型 {#poisson-regression}


## 线性回归回顾

从线性模型的数学记号
$$
y_n = \alpha + \beta x_n + \epsilon_n \quad \text{where}\quad
\epsilon_n \sim \operatorname{normal}(0,\sigma).
$$

等价于
$$
y_n - (\alpha + \beta x_n) \sim \operatorname{normal}(0,\sigma),
$$

又可以写为
$$
y_n \sim \operatorname{normal}(\alpha + \beta x_n, \, \sigma).
$$


线性回归需要满足四个前提假设：

1. **Linearity **
    - 因变量和每个自变量都是线性关系

2. **Indpendence **
    - 对于所有的观测值，它们的误差项相互之间是独立的

3. **Normality **
    - 误差项服从正态分布

4. **Equal-variance **  
    - 所有的误差项具有同样方差

这四个假设的首字母，合起来就是**LINE**，这样很好记


把这**四个前提**画在一张图中
```{r out.width = '80%', fig.align='center', echo = FALSE}
knitr::include_graphics(path = "images/LINE.png")
```



## 案例

我们从一个有意思的案例开始。

> 在受污染的岛屿附近，金枪鱼出现次数

```{r out.width = '50%', fig.align='center', echo = FALSE}
knitr::include_graphics(path = "images/fishes.png")
```

```{r}
library(tidyverse)

df <- read_rds("./demo_data/fish.rds")
df
```
我们的问题是，污染如何影响鱼类的数量？具体来说是想：建立不同位置金枪**鱼的数量** 与 这个位置的**污染程度**之间的线性关系。

### 线性模型的局限性

先看看变量之间的关系
```{r}
df %>%
  ggplot(aes(x = pollution_level, y = number_of_fish)) +
  geom_point() +
  geom_smooth(method = lm) +
  labs(
    title = "Number of fish counted under different pollution level",
    x = "Pollution level",
    y = "Number of fish counted"
  )
```

线性关系不明显，而且被解释变量甚至出现了负值。


我们再看看线性模型的结果
```{r}
m0 <- lm(number_of_fish ~ pollution_level, data = df)
summary(m0)
```
线性模型的失灵！ 怎么办呢？




### 泊松分布

我们再看看被解释变量的分布：

- 非负的整数0, 1, 2, 3, 4, ...

```{r}
df %>%
  ggplot(aes(x = number_of_fish)) +
  geom_histogram() +
  labs(
    title = "number of fishes (Poisson distribution)"
  )
```

这是典型的泊松分布。

```{r}
generate_pois <- function(lambda_value) {
  tibble(
    lambda = as.character(lambda_value),
    x = seq(1, 10),
    d = dpois(x = x, lambda = lambda_value)
  )
}

tb <- seq(0.1, 1.8, by = 0.2) %>% map_dfr(generate_pois)

tb %>%
  ggplot(aes(x = x, y = d, color = lambda)) +
  geom_point() +
  geom_line() +
  scale_x_continuous(breaks = seq(1, 10, 1)) +
  theme_bw()
```


事实上，生活中很多场合都会遇到计数、二进制、`yes/no`、等待时间等类型的数据，比如

1. 医院每天急诊次数
2. 每年摩托车死亡人数
3. 城市发生火灾的次数

他们有个共同的特征

- **变量**代表单位时间或者区域事件发生的次数，服从泊松分布。泊松分布有什么特点？

$$
y_i \sim \text{Poisson}(\lambda = \lambda_i)
$$




### 正态分布换成泊松分布就行了？

回到我们目的：建立不同位置 **鱼的数量** 与 这个位置的**污染程度**之间线性关系。



在之前线性模型中讲到，对每一次观测，被解释变量服从正态分布，那么，我们用解释变量的线性组合**模拟**正态分布的均值$\mu_i$，即均值$\mu_i$随$x_i$变化

$$
\begin{align*}
y_i \sim & \operatorname{normal}(\mu_i, \, \sigma)\\
    &\operatorname{normal}(\mu_i = \beta_0 + \beta_1 x_i, \, \sigma)
\end{align*}
$$


我们也想如法炮制，正态分布换成**泊松分布**
$$
\begin{align*}
y_i \sim & \text{Poisson}(\lambda_i)\\
         & \text{Poisson}(\lambda_i = \beta_0 + \beta_1 x_i)
\end{align*}
$$

但很遗憾，出现了问题

-  泊松分布的$\lambda_i$ 要求大于等于0，然而$\beta_0 + \beta_1 X_i$势必会出现负数

<!-- -  线性模型要求，所有的误差项具有同样方差；但是泊松分布$E(Y)=Var(Y)=\lambda$，均值越大方差也就越大。线性回归无法满足这个要求。 -->
<!-- $$ -->
<!-- \begin{align*} -->
<!-- \text{Poisson}(\lambda_i) \sim\epsilon_i = y_i - \lambda_i = y_i - (\alpha + \beta x_i) -->
<!-- \end{align*} -->
<!-- $$ -->
<!-- $$ -->
<!-- \begin{align*} -->
<!-- \text{Poisson}(\lambda_i) \sim\epsilon_i = y_i - \lambda_i = y_i - \exp(\alpha + \beta x_i) -->
<!-- \end{align*} -->
<!-- $$ -->

## 泊松回归模型

### 解决办法-连接函数之美

尽管不能使用直接线性模型，但可以间接使用，统计学家想到用log($\lambda_i$) 代替 $\lambda_i$，然后让解释变量的线性组合模拟log($\lambda_i$)，即 
$$
\begin{align*}
y_i \sim & \text{Poisson}(\lambda_i)\\
 \color{red}\log(\lambda_i) = & \beta_0 + \beta_1 x_i
\end{align*}
$$

现在问题迎刃而解：
- log($\lambda_i$) 值域范围是($-\infty$ to $\infty$)，这样既保证了$\lambda_i$ 是正值，又保证了$\beta_0 + \beta_1 x_i$ 可能出现的负值

- 这里的`log()`函数就是**连接函数**， 连接了$x_i$和$\lambda_i$








所以**泊松回归模型**为

$$
\begin{align*}
\log(\lambda_i) = & \beta_0 + \beta_1 x_i
\end{align*}
$$

注意到，这里没有线性回归模型中的误差项。通过极大似然估计（Maximum Likelihood Estimation）计算系数$(\beta)$


```{block, type="danger"}
什么叫最大相似性估计？通俗点讲，给定y的分布以及独立变量(x)的值，那么最有可能的$\beta$系数多是多少？

- step 1  $y$ 服从泊松分布
$$
\operatorname{Pr}\left(y_{i} | \lambda\right)=\frac{\lambda^{y_{i}} e^{-\lambda}}{y_{i} !}, \quad y=0,1,2, \ldots \quad ; \quad \lambda>0
$$

- step 2 回归模型
$$
E\left[y_{i} | x_{i}\right]=\lambda_{i}=\exp \left(x_{i}^{\prime} \beta\right)
$$


- step 3  $\lambda_{i}$ 代入上式，考虑观测值彼此独立，可以将所有$y_i$观测值相乘，

$$
\operatorname{Pr}\left(y_{1}, \ldots, y_{N} | x_{1}, \ldots, x_{N}\right)=\prod_{i=1}^{N} \frac{e^{y_{i} x_{i}^{\prime} \beta} e^{-e^{x_{i}^{\prime}} \beta}}{y_{i} !}
$$


- step 4 然后取对数，得到(joint log-likelihood function)

$$
l\left(\beta | y_{i}, X_{i}\right)=\sum_{i=1}^{n} y_{i} X_{i}^{\prime} \beta-\exp X_{i}^{\prime} \beta-\log y_{i} !
$$


- step 5 求偏导
$$
\frac{\partial l}{\partial \beta}=\sum_{i=1}^{n}\left(y_{i}-\exp X_{i}^{\prime} \beta\right) X_{i}
$$

- step 6 令等式等于0，通过样本值，可以求出系数。
``` 
  
  
  
  
  

### 代码实现

使用`glm()`函数:

```{r, eval = FALSE}
glm(y ~ 1 + x, family = familytype(link = linkfunction), data = )
```

- **formula**: 被解释变量 ~ 解释变量
- **family** : 误差分布（和连接函数），`family = poisson(link="log")`
- **data**   : 数据框




```{r}
m <- glm(number_of_fish ~ pollution_level,
  family = poisson(link = "log"),
  data = df
)
m
```





```{r}
summary(m)
```

```{r}
confint(m)
```




```{r}
broom::tidy(m)
```


## 模型解释

我们建立的模型是
$$
\begin{align*}
y_i \sim & \text{Poisson}(\lambda_i)\\
\log(\lambda_i) = & \beta_0 + \beta_1 x_i
\end{align*}
$$


我个人比较喜欢这样写
$$
\begin{align*}
y_i \sim & \;\text{Poisson}(\lambda_i = \exp(\beta_0 + \beta_1 x_i))
\end{align*}
$$


### 系数

$$
\begin{align*}
\frac{\lambda_1}{\lambda_0} & =
\frac{\exp(\beta_1 + \beta_1*x_1)}{\exp(\beta_0 + \beta_1*x_0)} \\
& = \exp(\beta_1(x_1 - x_0))
\end{align*}
$$

计算当$x_i$增加一个单位时，事件平均发生次数将会是原来的$\exp(\beta_1)$倍。具体系数为：

```{r}
coef(m)
exp(coef(m)[2])
exp(coef(m))
```

- 污染系数为0， `r exp(coef(m)[1])`
- 污染系数从0变到0.5, 引起 `(1/exp(-3.1077*0.5) = 4.7)`倍数的鱼数量下降.
- 污染系数从0变到1,   引起 `r 1/exp(coef(m)[2])` 倍数的鱼数量下降.




### 边际效应(Marginal effects)

即在其他条件不变的情况下，$x_i$增加一个单位，事件的平均发生次数会增加 $100\beta_1 \%$：
类似求偏导$\frac{\partial{\lambda}}{\partial{x}}$


```{r}
margins::margins(m, type = "link")
```


模型是非线性的，所以我们常用更直观的方式评估边际效应，即，自变量直接对因变量的贡献，可以令`type = "response"`，类似求偏导$\frac{\partial{y}}{\partial{x}}$



```{r}
margins::marginal_effects(m, type = "response", se = TRUE) %>%
  as.data.frame() %>%
  dplyr::mutate(pollution_level = df$pollution_level) %>%
  ggplot(aes(x = pollution_level, y = dydx_pollution_level)) +
  geom_point()
```

纵坐标是事件的平均发生次数（增加）下降的比例，可以看到随着x变大，下降趋缓。更多边际效应的内容，可参考[这里](https://cran.r-project.org/web/packages/margins/vignettes/Introduction.html)

### 拟合

```{r}
fitted(m) %>% 
  head()
```


实质上就是$exp(\beta_0 +  \beta_1 * pollution_level)$

```{r}
intercept <- coef(m)[1]
beta <- coef(m)[2]

df %>%
  dplyr::mutate(theory_pred = fitted(m)) %>%
  dplyr::mutate(
    myguess_pred = exp(intercept + beta * pollution_level)
  )
```



```{r}
df %>%
  dplyr::mutate(theory_pred = fitted(m)) %>%
  ggplot(aes(x = pollution_level, y = theory_pred)) +
  geom_point()
```


```{r}
pred <- predict(m, type = "response", se = TRUE) %>% as.data.frame()
pred %>% 
  head()
```
 


```{r}
df_pred <- df %>%
  dplyr::mutate(
    fit = pred$fit,
    se_fit = pred$se.fit
  )
df_pred
```


```{r}
real_df <-
  tibble(
    x = seq(0, 1, length.out = 100),
    y = 4 * exp(-3.2 * x)
  )

df_pred %>%
  ggplot(aes(x = pollution_level, y = number_of_fish)) +
  geom_point() +
  geom_pointrange(aes(
    y = fit,
    ymax = fit + se_fit,
    ymin = fit - se_fit
  ), color = "red") +
  # geom_point(aes(y = fit + se_fit), color = "red") +
  # geom_point(aes(y = fit - se_fit), color = "red") +
  geom_line(data = real_df, aes(x = x, y = y), color = "black") +
  labs(
    title = "Number of fish counted under different pollution level",
    x = "Pollution level",
    y = "Number of fish counted"
  )
```


### 模型评估

```{r}
knitr::include_graphics(path = "images/model_evaluation.png")
```

- 过度离散，负二项式分布模型
- 零膨胀

```{r, eval=FALSE}
margins::margins(m)

ggeffects::ggpredict(m)

ggeffects::ggpredict(m, terms = c("pollution_level"))

performance::model_performance(m)

performance::check_model(m)

performance::check_overdispersion(m)

performance::check_zeroinflation(m)
```


## 思考

### 这两者为什么是相等的？

```{r}
d <- tibble(
  x = 1:100,
  y = 4 + 2 * x + rnorm(100)
)

lm(y ~ x, data = d)

glm(y ~ x,
  family = gaussian(link = "identity"),
  data = d
)
```


> `lm` 是 `glm`的一种特殊情形。



### `log link` 与 `log transforming`  

在案例中

- log link

```{r, eval=FALSE}
glm(number_of_fish ~ pollution_level,
  family = gaussian(link = "log"),
  data = df
)
```


- 先对 number_of_fish 取对数后，然后线性回归 

```{r, eval=FALSE}
# lm(log(number_of_fish) ~ pollution_level, data = df)
glm(log(number_of_fish) ~ pollution_level,
  family = gaussian(link = "identity"),
  data = df
)
```


**这两者有什么区别？**

比较两者的结果，发现有很大的差别，尤其是斜率，为什么呢？ 因为这是两个不同的模型：

- 第一个模型中 `link="log"`，模型并没有直接变换数据，**只是使用了原始数据的均值**，均值对数计算后，建立线性关系，即

$$
\begin{align*}
\log(\lambda_i)  & = \beta_0 + \beta_1 x_i \\
\lambda_i & = \exp(\beta_0 + \beta_1 x_i)
\end{align*}
$$

注意到这里`family = gaussian`, 因此，误差项满足高斯分布

$$
\begin{align*}
y_i - \lambda_i &\sim \operatorname{normal}(0,\sigma)\\
y_i - \exp(\beta_0 + \beta_1 x_i) &\sim \operatorname{normal}(0,\sigma) \\
y_i &\sim \operatorname{normal}(\exp(\beta_0 + \beta_1 x_i), \, \sigma)
\end{align*}
$$


$$
y_i = \exp(\beta_0 + \beta_1 x_i) + \epsilon_i\quad \epsilon_i \sim \operatorname{normal}(0,\sigma)
$$

因此对不同均值$\lambda_i$，误差项的方差是一样的



- 第二个模型 `log transforming`，是直接转换原始数据，然后建立模型$\log(y_i) = \alpha + \beta x_i$，**原始数据y的均值和方差都改变了**，一旦log(y) 变回 y时，

$$
\log(y_i) =\beta_0 + \beta_1 x_i + \epsilon_i, \quad \epsilon_i \sim \operatorname{normal}(0,\sigma)
$$


$$
y_i = \exp(\beta_0 + \beta_1 x_i) * \exp(\epsilon_i), \quad \epsilon_i \sim \operatorname{normal}(0,\sigma)
$$

误差的方差也会随着均值变化。



### 更多分布
```{r, eval=FALSE}
x <- c(1, 2, 3, 4, 5)
y <- c(1, 2, 4, 2, 6)

regNId <- glm(y ~ x, family = gaussian(link = "identity"))
regNlog <- glm(y ~ x, family = gaussian(link = "log"))
regPId <- glm(y ~ x, family = poisson(link = "identity"))
regPlog <- glm(y ~ x, family = poisson(link = "log"))
regGId <- glm(y ~ x, family = Gamma(link = "identity"))
regGlog <- glm(y ~ x, family = Gamma(link = "log"))
regIGId <- glm(y ~ x, family = inverse.gaussian(link = "identity"))
regIGlog <- glm(y ~ x, family = inverse.gaussian(link = "log"))
```

```{r, eval=FALSE}
dx <- tibble(
  x = c(1, 2, 3, 4, 5),
  y = c(1, 2, 4, 2, 6)
)

dx %>%
  ggplot(aes(x = x, y = y)) +
  geom_point()
```


```{r, eval=FALSE}
regNId <- glm(y ~ x, family = gaussian(link = "identity"), data = dx)
regNId

dx %>%
  mutate(pred = predict(regNId, type = "response")) %>%
  ggplot(aes(x = x, y = y)) +
  geom_point() +
  geom_line(aes(y = pred, group = 1))
```




```{r, eval=FALSE}
regPlog <- glm(y ~ x, family = poisson(link = "log"), data = dx)
regPlog

dx %>%
  mutate(pred = predict(regPlog, type = "response")) %>%
  ggplot(aes(x = x, y = y)) +
  geom_point() +
  geom_line(aes(y = pred, group = 1))
```


```{r, eval=FALSE}
regGId <- glm(y ~ x, family = Gamma(link = "identity"), data = dx)
regGId

dx %>%
  mutate(pred = predict(regGId, type = "response")) %>%
  ggplot(aes(x = x, y = y)) +
  geom_point() +
  geom_line(aes(y = pred, group = 1))
```

### 更复杂的模型
以后再说
```{r, eval=FALSE}
glm(number_of_fish ~ 1 + (1 | pollution_level),
  family = poisson(link = "log"),
  data = df
)
```



## 小结

```{r}
knitr::include_graphics(path = "images/One_Picture.png")
```

第 \@ref(logistic-regression) 章接着讲广义线性模型中的logistic回归模型。


<!-- ## 贝叶斯泊松回归 -->

<!-- ```{r} -->
<!-- m <- glm(number_of_fish ~ pollution_level,  -->
<!--          family = poisson(link = "log"), -->
<!--          data = df) -->
<!-- m -->
<!-- ``` -->


<!-- ```{r} -->
<!-- library(brms) -->
<!-- m_brms <- brm(number_of_fish ~ pollution_level,  -->
<!--               family = poisson(link = "log"), -->
<!--               data = df,  -->
<!--               prior = c(set_prior("normal(0,100)", class="b")) -->
<!--               ) -->
<!-- m_brms -->
<!-- ``` -->


<!-- ```{r} -->
<!-- plot(m_brms) -->
<!-- ``` -->


<!-- ```{r} -->
<!-- df %>% tidybayes::add_fitted_draws(m_brms) -->
<!-- ``` -->


<!-- ```{r} -->
<!-- conditional_effects(m_brms) -->
<!-- ``` -->




```{r, echo = F}
# remove the objects
# rm(list=ls())
rm(beta, d, df, df_pred, generate_pois, intercept, m, m0, pred, real_df, tb)
```

```{r, echo = F, message = F, warning = F, results = "hide"}
pacman::p_unload(pacman::p_loaded(), character.only = TRUE)
```