generalized-linear-models.qmd

# 广义线性模型 {#sec-generalized-linear-models}

```{r}
#| echo: false

Sys.setenv(CMDSTANR_NO_VER_CHECK = TRUE)
source("_common.R")
```

## 生成模拟数据 {#sec-simulate-poisson-data}

先介绍泊松广义线性模型，包括模拟和计算，并和 Stan 实现的结果比较。

泊松广义线性模型如下：

$$
\begin{aligned}
\log(\lambda) &= \beta_0 + \beta_1 x_1 + \beta_2 x_2 \\
Y &\sim \mathrm{Poisson}(u\lambda)
\end{aligned}
$$

设定参数向量 $\beta = (\beta_0, \beta_1, \beta_2) = (0.5, 0.3, 0.2)$，观测变量 $X_1$ 和 $X_2$ 的均值都为 0，协方差矩阵 $\Sigma$ 为

$$
\left[
 \begin{matrix}
   1.0 & 0.8  \\
   0.8 & 1.0 
 \end{matrix}
\right]
$$

模拟观测到的响应变量值和协变量值，添加漂移项

```{r}
set.seed(2023)
n <- 2500 # 样本量
beta <- c(0.5, 0.3, 0.2)
X <- MASS::mvrnorm(n, mu = rep(0, 2), Sigma = matrix(c(1, 0.8, 0.8, 1), 2))
u <- rep(c(2, 4), each = n / 2)
lambda <- u * exp(cbind(1, X) %*% beta)
y <- rpois(n, lambda = lambda)
```

## 拟合泊松模型 {#sec-poisson-model}

拟合泊松回归模型

```{r}
fit_poisson_glm <- glm(y ~ X, family = poisson(link = "log"), offset = log(u))
summary(fit_poisson_glm)
```

```{r}
# 对数似然函数值
log_poisson_lik <- logLik(fit_poisson_glm)
# 计算 AIC AIC(fit_poisson_glm)
-2 * c(log_poisson_lik) + 2 * attr(log_poisson_lik, "df")
```

下面用 Stan 编码泊松回归模型，模型代码如下：

```{verbatim, file="code/poisson_log_glm.stan", lang="stan"}
```

Stan 代码主要分三部分：

1.  数据部分 `data`：声明模型的输入数据，数据类型、大小、约束。

2.  参数部分 `parameters`：类似数据部分，声明模型的参数，参数类型、大小。

3.  模型部分 `model`：指定模型参数的先验分布。

4.  生成量 `generated quantities`：拟合模型获得参数估计值后，计算一些统计量。

下面准备数据

```{r}
nchains <- 4 # 4 条迭代链
# 给每条链设置不同的参数初始值
inits_data <- lapply(1:nchains, function(i) {
  list(
    alpha = runif(1, 0, 1),
    beta = runif(2, 1, 10)
  )
})

# 准备数据
poisson_d <- list(
  n = 2500, # 观测记录的条数
  k = 2, # 协变量个数
  X = X, # N x 2 矩阵
  y = y, # N 向量
  log_offset = log(u)
)
```

编译模型，抽样获取参数的后验分布

```{r}
#| message: false

# 加载 cmdstanr 包
library(cmdstanr)
# 编译模型
mod_poisson <- cmdstan_model(
  stan_file = "code/poisson_log_glm.stan",
  compile = TRUE,
  cpp_options = list(stan_threads = TRUE)
)
# 采样拟合模型
fit_poisson_stan <- mod_poisson$sample(
  data = poisson_d, # 观测数据
  init = inits_data, # 迭代初值
  iter_warmup = 1000, # 每条链预处理迭代次数
  iter_sampling = 2000, # 每条链总迭代次数
  chains = nchains, # 马尔科夫链的数目
  parallel_chains = 1, # 指定 CPU 核心数，可以给每条链分配一个
  threads_per_chain = 1, # 每条链设置一个线程
  show_messages = FALSE, # 不显示迭代的中间过程
  refresh = 0, # 不显示采样的进度
  seed = 20222022 # 设置随机数种子，不要使用 set.seed() 函数
)
# 迭代诊断
fit_poisson_stan$diagnostic_summary()
# 输出结果
fit_poisson_stan$summary(c("alpha", "beta", "lp__"))
```

## 参数后验分布 {#sec-posterior-distribution}

加载 **bayesplot** 包，bayesplot 包提供一系列描述数据分布的绘图函数，比如绘制散点图 `mcmc_scatter()` 。$\beta_1$ 和 $\beta_2$ 的联合分布

```{r}
#| label: fig-stan-scatter
#| fig-cap: $\beta_1$ 和 $\beta_2$ 的联合分布
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 4
#| message: false

library(ggplot2)
library(bayesplot)
mcmc_scatter(fit_poisson_stan$draws(c("beta[1]", "beta[2]")), size = 1) +
  theme_classic() +
  labs(x = expression(beta[1]), y = expression(beta[2]))
```

如果提取采样的数据，也可使用 ggplot2 包绘图，不局限于 bayesplot 设定的风格。

```{r}
#| label: fig-density-filled
#| fig-cap: $\beta_1$ 和 $\beta_2$ 的联合分布
#| fig-showtext: true
#| fig-width: 6
#| fig-height: 5

beta_df <- fit_poisson_stan$draws(c("beta[1]", "beta[2]"), format = "draws_df")
ggplot(data = beta_df, aes(x = `beta[1]`, y = `beta[2]`)) +
  geom_density_2d_filled() +
  facet_wrap(~.chain, ncol = 2) +
  theme_classic() +
  labs(x = expression(beta[1]), y = expression(beta[2]))
```

$\beta_1$ 和 $\beta_2$ 的热力图

```{r}
#| label: fig-stan-hex
#| fig-cap: $\beta_1$ 和 $\beta_2$ 的热力图
#| fig-showtext: true
#| fig-width: 5.5
#| fig-height: 4

mcmc_hex(fit_poisson_stan$draws(c("beta[1]", "beta[2]"))) +
  theme_classic() +
  labs(x = expression(beta[1]), y = expression(beta[2]))
```

各个参数的轨迹图

```{r}
#| label: fig-stan-trace
#| fig-cap: 各个参数的轨迹图
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 5

mcmc_trace(fit_poisson_stan$draws(c("beta[1]", "beta[2]")),
  facet_args = list(
    labeller = ggplot2::label_parsed, strip.position = "top", ncol = 1
  )
) +
  theme_classic()
```

可以将模型参数的后验分布图展示出来

```{r}
#| label: fig-stan-dens
#| fig-cap: 各个参数的分布图（密度图）
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 5

mcmc_dens(fit_poisson_stan$draws(c("beta[1]", "beta[2]")),
  facet_args = list(
    labeller = ggplot2::label_parsed, strip.position = "top", ncol = 1
  )
) +
  theme_classic()
```

后验分布的中位数、80% 区间

```{r}
#| label: fig-stan-areas
#| fig-cap: 各个参数的分布图（岭线图）
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 4
#| message: false

mcmc_areas(fit_poisson_stan$draws(c("beta[1]", "beta[2]")), prob = 0.8) +
  scale_y_discrete(labels = scales::parse_format()) +
  theme_classic()
```

岭线图就是将各个参数的后验分布图放在一起。

```{r}
#| label: fig-stan-ridges
#| fig-cap: 各个参数的分布图（岭线图）
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 4
#| message: false

mcmc_areas_ridges(x = fit_poisson_stan$draws(), pars = c("beta[1]", "beta[2]")) +
  scale_y_discrete(labels = scales::parse_format()) +
  theme_classic()
```

参数的 $\hat{R}$ 潜在尺度收缩因子

```{r}
bayesplot::rhat(fit_poisson_stan, pars = "alpha")
```

后验预测诊断的想法是检查根据拟合模型生成的随机数 $y^{rep}$ 与真实观测数据 $y$ 的接近程度。为直观起见，可以用一系列描述数据分布的图来可视化检验。

```{r}
#| label: fig-stan-nuts
#| fig-cap: NUTS 能量诊断图
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 4
#| message: false

# mcmc_scatter(fit_poisson_stan$draws(),
#   pars = c("beta[1]", "beta[2]"),
#   np = nuts_params(fit_poisson_stan)
# )

mcmc_nuts_energy(x = nuts_params(fit_poisson_stan), binwidth = 1) +
  ggtitle(label = "NUTS Energy Diagnostic")
```

y 是真实数据，yrep 是根据贝叶斯拟合模型生成的数据。下图是真实数据的密度图和50组生成数据的密度图。

```{r}
#| label: fig-stan-ppcheck-dens
#| fig-cap: 后验预测诊断图（密度图）
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 4

# 抽取 yrep 数据
yrep <- fit_poisson_stan$draws(variables = "y_rep", format = "draws_matrix")
pp_check(y, yrep = yrep[1:50, ], fun = ppc_dens_overlay) +
  theme_classic()
```

观察后验预测区间与真实数据的覆盖情况，不妨取前 50 次观测的数据，即 `y[1:50]` 与第 2 个自变量 `X[1:50, 2]` ，基于后验分布的 500 次采样数据绘制 50% 后验置信区间。

```{r}
#| label: fig-stan-ppcheck-intervals
#| fig-cap: 后验预测诊断图（区间图）
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 4

ppc_intervals(y[1:50], yrep = yrep[1:1000, 1:50], x = X[1:50, 2], prob = 0.5)
```

## 模型评估指标 {#sec-model-evaluation}

**loo** 包可以计算 WAIC

```{r}
fit_poisson_waic <- loo::waic(fit_poisson_stan$draws(variables = "log_lik"))
print(fit_poisson_waic)
```

**loo** 包推荐使用 LOO-CV ，它还提供诊断信息、有效样本量和蒙特卡罗估计。

```{r}
fit_poisson_loo <- fit_poisson_stan$loo(variables = "log_lik", cores = 2)
print(fit_poisson_loo)
```

## 可选替代实现 {#sec-bayesian-brms}

对于常见的统计模型，rstanarm 和 **brms** 包都内置了预编译的 Stan 程序，下面用 **brms** 包的函数 `brm()` 拟合带上述漂移项的泊松广义线性模型，参数估计结果和 Base R 函数 `glm()` 的几乎一致，因编译和抽样的过程比较花费时间，速度不及 Base R。

``` r
# brms
dat <- data.frame(y = y, X = X, u = u)
colnames(dat) <- c("y", "x1", "x2", "u")
fit_poisson_brm <- brms::brm(y ~ x1 + x2 + offset(log(u)),
  data = dat, family = poisson(link = "log")
)
fit_poisson_brm
```

```         
 Family: poisson 
  Links: mu = log 
Formula: y ~ x1 + x2 + offset(log(u)) 
   Data: dat (Number of observations: 2500) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     0.49      0.01     0.47     0.51 1.00     2509     2171
x1            0.29      0.01     0.26     0.32 1.00     1771     1645
x2            0.21      0.01     0.19     0.24 1.00     1727     1847

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
```

调用函数 `brm()` 拟合模型后返回一个 brmsfit 对象 `fit_poisson_brm`，**brms** 包提供很多函数处理该数据对象，比如 `brms::loo()` 计算 LOO-CV

``` r
brms::loo(fit_poisson_brm)
```

```         
Computed from 4000 by 2500 log-likelihood matrix

         Estimate   SE
elpd_loo  -5386.3 37.8
p_loo         2.9  0.1
looic     10772.6 75.5
------
Monte Carlo SE of elpd_loo is 0.0.

All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.
```

输出结果中， LOO IC 信息准则 Loo information criterion，looic 指标的作用类似频率派模型中的 AIC 指标，所以也几乎相同的。

``` r
# 后验预测检查
brms::pp_check(fit_poisson_brm)
```

## 案例：吸烟喝酒和食道癌的关系 {#sec-esoph}

<!-- 存在有序分类变量数据  -->

本例数据集 esoph 来自 Base R 内置的 datasets 包，是法国伊勒-维莱讷食道癌研究数据，研究吸烟、喝酒与食道癌的关系，量化酒精、烟草、酒精和烟草的交互作用。部分数据集见 @tbl-esoph ，年龄组 agegp、酒精量 alcgp 和烟草量 tobgp 为有序的分类变量，正常来说，年龄越大，吸烟、喝酒对食道癌影响越大。

```{r}
#| label: tbl-esoph
#| tbl-cap: "食道癌研究数据（部分）"
#| echo: false

knitr::kable(head(esoph), col.names = c("年龄组", "酒精量", "烟草量", "实验组", "控制组"))
```

### 描述分析

先来简单统计一下各年龄组、酒精量组的食道癌发病人数

```{r}
xtabs(data = esoph, cbind(ncases, ncontrols) ~ agegp + alcgp)
```

@fig-esoph 描述食道癌发病率与年龄组、酒精量的关系

```{r}
#| label: fig-esoph
#| fig-cap: "食道癌发病率与年龄组、酒精量的关系"
#| fig-showtext: true
#| fig-width: 5
#| fig-height: 3.5

library(ggplot2)
aggregate(cbind(ncases, ncontrols) ~ agegp + alcgp, data = esoph, sum) |>
  ggplot(aes(x = agegp, y = alcgp, fill = ncases / (ncases + ncontrols))) +
  scale_fill_viridis_c(labels = scales::percent_format()) +
  geom_tile() +
  labs(x = "年龄组", y = "酒精量", fill = "发病率")
```

### 拟合模型

响应变量服从二项分布，自变量包含年龄分组 agegp、酒精量 alcgp、烟草量 tobgp 和 酒精量与烟草量的交互作用，建立广义线性模型。

```{r}
fit_glm_esoph <- glm(cbind(ncases, ncontrols) ~ agegp + tobgp * alcgp,
  data = esoph, family = binomial(link = "logit")
)
```

模型输出

```{r}
summary(fit_glm_esoph)
```

整理模型输出后，见 @tbl-glm-esoph

```{r}
#| label: tbl-glm-esoph
#| tbl-cap: "广义线性模型各个参数的估计结果"
#| echo: false

knitr::kable(broom::tidy(fit_glm_esoph), align = "lrrrr")
```

### 与 brms 比较

下面从贝叶斯的视角分析和建模，使用 **brms** 包对该数据拟合，同样是广义线性模型。

``` r
fit_brm_esoph <- brm(ncases | trials(ncases + ncontrols) ~ agegp + tobgp * alcgp, 
                     data = esoph, family = binomial(link = "logit"))
```

```         
 Family: binomial 
  Links: mu = logit 
Formula: ncases | trials(ncases + ncontrols) ~ agegp + tobgp * alcgp 
   Data: esoph (Number of observations: 88) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept          -1.91      0.25    -2.49    -1.51        735 1.01
agegp.L             3.39      0.86     2.13     5.45        674 1.01
agegp.Q            -1.68      0.78    -3.58    -0.50        658 1.01
agegp.C             0.31      0.57    -0.59     1.63        709 1.00
agegpE4            -0.01      0.36    -0.80     0.65        907 1.01
agegpE5            -0.20      0.21    -0.59     0.22       1970 1.00
tobgp.L             0.63      0.20     0.24     1.03       4654 1.00
tobgp.Q             0.03      0.20    -0.38     0.42       3469 1.00
tobgp.C             0.17      0.20    -0.21     0.57       3892 1.00
alcgp.L             1.41      0.22     0.99     1.84       4067 1.00
alcgp.Q            -0.16      0.20    -0.56     0.24       3335 1.00
alcgp.C             0.25      0.19    -0.12     0.62       3870 1.00
tobgp.L:alcgp.L    -0.69      0.42    -1.51     0.16       3878 1.00
tobgp.Q:alcgp.L     0.13      0.43    -0.75     0.97       4249 1.00
tobgp.C:alcgp.L    -0.30      0.44    -1.15     0.58       5149 1.00
tobgp.L:alcgp.Q     0.13      0.41    -0.67     0.94       3127 1.00
tobgp.Q:alcgp.Q    -0.46      0.41    -1.24     0.34       4037 1.00
tobgp.C:alcgp.Q    -0.05      0.40    -0.82     0.74       4490 1.00
tobgp.L:alcgp.C    -0.15      0.38    -0.89     0.58       3507 1.00
tobgp.Q:alcgp.C     0.04      0.37    -0.69     0.75       3274 1.00
tobgp.C:alcgp.C    -0.17      0.36    -0.88     0.54       3773 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).
```

输出结果和 `glm()` 有不少差别的。

## 案例：哥本哈根住房状况调查

<!-- 响应变量是分类有序的变量 -->

数据集 housing 哥本哈根住房状况调查中的次数分布表，`Sat` 住户对目前居住环境的满意程度，是一个有序的因子变量，`Infl` 住户对物业管理的感知影响程度，`Type` 租赁住宿类型，如塔楼、中庭、公寓、露台，`Cont` 联系居民可与其他居民联系(低、高)，`Freq` 每个类中的居民人数，调查的人数。

```{r}
data("housing", package = "MASS")
str(housing)
```

响应变量是居民对居住环境满意度 Sat ，分三个等级，且存在强弱，等级，大小之分。

```{r}
# 因子变量的处理
options(contrasts = c("contr.treatment", "contr.poly"))
# 有序逻辑回归
housing_mass <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE)
summary(housing_mass)
```

计算置信区间

```{r}
# 剖面
confint(profile(housing_mass), level = 0.95)
```

## 习题 {#sec-bayesian-exercises}

1.  分析挑战者号航天飞机 O 型环数据。**DAAG** 包的 orings 数据集记录美国挑战者号航天飞机 O 型环在不同温度下发生 Erosion 腐蚀和 Blowby 串气的失效数量。 @fig-cdplot-orings 展示航天飞机 O 型环在不同温度下失效的分布图（条件密度图）：随着温度升高，O 型环越来越不容易失效。请分别用 Base R 函数 `glm()` 和 **cmdstanr** 包建模分析 O 型环数据。

    ```{r}
    #| label: fig-cdplot-orings
    #| fig-cap: 航天飞机 O 型环在不同温度下失效的条件密度图
    #| fig-showtext: true
    #| fig-width: 6
    #| fig-height: 4
    #| code-fold: true
    #| echo: !expr knitr::is_html_output()

    # data(orings, package = "DAAG")
    orings <- readRDS(file = "data/orings.rds")
    ggplot(orings, aes(x = Temperature, y = after_stat(count))) +
      geom_density(aes(fill = Total > 0), position = "fill", bw = 2) +
      scale_y_continuous(labels = scales::label_percent()) +
      scale_fill_grey(labels = c("TRUE" = "是", "FALSE" = "否")) +
      theme_classic() +
      labs(x = "温度", y = "比例", fill = "失效")
    ```

2.  基于数据集 infert 分析自然流产和人工流产后的不育情况，

    ```{r}
    #| eval: false
    #| code-fold: true
    #| echo: !expr knitr::is_html_output()

    infert_glm <- glm(
      case ~ age + parity + education + spontaneous + induced,
      data = infert, family = binomial()
    )
    summary(infert_glm)

    # conditional logistic regression
    library(survival)
    infert_survival <- clogit(
      case ~ age + parity + education + spontaneous + induced + strata(stratum), data = infert
    )
    summary(infert_survival)
    ```

3.  根据 @sec-nuclear-pollution-concentration 的数据，建立贝叶斯空间广义线性混合模型，用 Stan 预测核辐射强度的分布。