hierarchical-normal-models.qmd

# 分层正态模型 {#sec-hierarchical-normal-models}

```{r}
#| echo: false

Sys.setenv(CMDSTANR_NO_VER_CHECK = TRUE)
# https://github.com/r-lib/processx/issues/236
Sys.setenv(PROCESSX_NOTIFY_OLD_SIGCHLD = TRUE)
source("_common.R")
```

> This is a bit like asking how should I tweak my sailboat so I can explore the ocean floor.
>
> --- Roger Koenker [^hierarchical-normal-models-1]

[^hierarchical-normal-models-1]: <https://stat.ethz.ch/pipermail/r-help/2013-May/354311.html>

乔治·博克斯说，所有的模型都是错的，但有些是有用的。在真实的数据面前，尽我们所能，结果发现没有最好的模型，只有更好的模型。总是需要自己去构造符合自己需求的模型及其实现，只有自己能够实现，才能在模型的海洋中畅快地遨游。

介绍分层正态模型的定义、结构、估计，分层正态模型与曲线生长模型的关系，分层正态模型与潜变量模型的关系，分层正态模型与线性混合效应的关系。以 **rstan** 包和 **nlme** 包拟合分层正态模型，说明 **rstan** 包的一些用法，比较贝叶斯和频率派方法拟合的结果，给出结果的解释。再对比 16 个不同的 R包实现，总结一般地使用经验，也体会不同 R 包的独特性。

```{r}
#| message: false

library(StanHeaders)
library(ggplot2)
library(rstan)
# 将编译的 Stan 模型与代码文件放在一起
rstan_options(auto_write = TRUE)
# 如果CPU和内存足够，设置成与马尔科夫链一样多
options(mc.cores = 2)
# 调色板
custom_colors <- c(
  "#4285f4", # GoogleBlue
  "#34A853", # GoogleGreen
  "#FBBC05", # GoogleYellow
  "#EA4335"  # GoogleRed
)
rstan_ggtheme_options(
  panel.background = element_rect(fill = "white"),
  legend.position = "top"
)
rstan_gg_options(
  fill = "#4285f4", color = "white",
  pt_color = "#EA4335", chain_colors = custom_colors
)
library(bayesplot)
```

## rstan 包 {#sec-8schools-rstan}

本节以 8schools 数据为例介绍分层正态模型及 **rstan** 包实现，8schools 数据最早来自 @Rubin1981 ，分层正态模型如下：

$$
\begin{aligned}
y_j &\sim \mathcal{N}(\theta_j,\sigma_j^2) \quad
\theta_j = \mu + \tau \times \eta_j \\
\theta_j &\sim \mathcal{N}(\mu, \tau^2) \quad
\eta_j \sim \mathcal{N}(0,1) \\
\mu &\sim \mathcal{N}(0, 100^2) \quad \tau \sim \mathrm{half\_normal}(0,100^2)
\end{aligned}
$$

其中，$y_j,\sigma_j$ 是已知的观测数据，$\theta_j$ 是模型参数， $\eta_j$ 是服从标准正态分布的潜变量，$\mu,\tau$ 是超参数，分别服从正态分布（将方差设置为很大的数，则变成弱信息先验或无信息均匀先验）和半正态分布（随机变量限制为正值）。

### 拟合模型

用 **rstan** 包来拟合模型，下面采用非中心的参数化表示，降低参数的相关性，减少发散的迭代次数，提高采样效率。

```{r}
# 编译模型
eight_schools_fit <- stan(
  model_name = "eight_schools",
  # file = "code/eight_schools.stan",
  model_code = "
  // saved as eight_schools.stan
  data {
    int<lower=0> J;                // number of schools
    array[J] real y;               // estimated treatment effects
    array[J] real <lower=0> sigma; // standard error of effect estimates
  }
  parameters {
    real mu;                // population treatment effect
    real<lower=0> tau;      // standard deviation in treatment effects
    vector[J] eta;          // unscaled deviation from mu by school
  }
  transformed parameters {
    vector[J] theta = mu + tau * eta;        // school treatment effects
  }
  model {
    target += normal_lpdf(mu | 0, 100); 
    target += normal_lpdf(tau | 0, 100);
    target += normal_lpdf(eta | 0, 1);  // prior log-density
    target += normal_lpdf(y | theta, sigma); // log-likelihood
  }
  ",
  data = list( # 观测数据
    J = 8,
    y = c(28, 8, -3, 7, -1, 1, 18, 12),
    sigma = c(15, 10, 16, 11, 9, 11, 10, 18)
  ),
  warmup = 1000, # 每条链预处理迭代次数
  iter = 2000,   # 每条链总迭代次数
  chains = 2,    # 马尔科夫链的数目
  cores = 2,     # 指定 CPU 核心数，可以给每条链分配一个
  verbose = FALSE, # 不显示迭代的中间过程
  refresh = 0,     # 不显示采样的进度
  seed = 20232023  # 设置随机数种子，不要使用 set.seed() 函数
)
```

### 模型输出

用函数 `print()` 打印输出结果，保留 2 位小数。

```{r}
print(eight_schools_fit, digits = 2)
```

值得一提，数据有限而且规律不明确，数据隐含的信息不是很多，则先验分布的情况将会对参数估计结果产生很大影响。Stan 默认采用无信息的先验分布，当使用非常弱的信息先验时，结果就非常不同了。提取任意一个参数的结果，如查看参数 $\tau$ 的 95% 置信区间。

```{r}
print(eight_schools_fit, pars = "tau", probs = c(0.025, 0.975))
```

从迭代抽样数据获得与 `print(fit)` 一样的结果。以便后续对原始采样数据做任意的进一步分析。**rstan** 包扩展泛型函数 `summary()` 以支持对 stanfit 数据对象汇总，输出各个参数分链条和合并链条的后验分布结果。

### 操作数据

抽取数据对象 `eight_schools_fit` 中的采样数据，合并几条马氏链的结果，返回的结果是一个列表。

```{r}
eight_schools_sim <- extract(eight_schools_fit, permuted = TRUE)
```

返回列表中的每个元素是一个数组，标量参数对应一维数组，向量参数对应二维数组。

```{r}
str(eight_schools_sim)
```

对于列表，适合用函数 `lapply()` 配合算术函数计算 $\mu,\tau$ 等参数的均值。

```{r}
fun_mean <- function(x) {
  if (length(dim(x)) > 1) {
    apply(x, 2, mean)
  } else {
    mean(x)
  }
}
lapply(eight_schools_sim, FUN = fun_mean)
```

类似地，计算 $\mu,\tau$ 等参数的分位点。

```{r}
fun_quantile <- function(x, probs) {
  if (length(dim(x)) > 1) {
    t(apply(x, 2, quantile, probs = probs))
  } else {
    quantile(x, probs = probs)
  }
}
lapply(eight_schools_sim, fun_quantile, probs = c(2.5, 25, 50, 75, 97.5) / 100)
```

同理，可以计算最大值 `max()`、最小值 `min()` 和中位数 `median()` 等。

### 采样诊断

获取马尔科夫链迭代点列数据

```{r}
eight_schools_sim <- extract(eight_schools_fit, permuted = FALSE)
```

`eight_schools_sim` 是一个三维数组，1000（次迭代）\* 2 （条链）\* 19（个参数）。如果 `permuted = TRUE` 则会合并马氏链的迭代结果，变成一个列表。

```{r}
# 数据类型
class(eight_schools_sim)
# 1000（次迭代）* 2 （条链）* 19（个参数）
str(eight_schools_sim)
```

提取参数 $\mu$ 的迭代点列，绘制迭代轨迹。

```{r}
#| label: fig-8schools-mu-base
#| fig-cap: Base R 绘制参数 $\mu$ 的迭代轨迹
#| fig-showtext: true
#| par: true

eight_schools_mu_sim <- eight_schools_sim[, , "mu"]
matplot(
  eight_schools_mu_sim, xlab = "迭代次数", ylab = expression(mu),
  type = "l", lty = "solid", col = custom_colors
)
abline(h = apply(eight_schools_mu_sim, 2, mean), col = custom_colors)
legend(
  "topleft", legend = paste("chain", 1:2), box.col = "white", 
  inset = 0.01, lty = "solid", horiz = TRUE, col = custom_colors
)
```

也可以使用 **rstan** 包提供的函数 `traceplot()` 或者 `stan_trace()` 绘制参数的迭代轨迹图。

```{r}
#| label: fig-8schools-mu-ggplot2
#| fig-cap: rstan 绘制参数 $\mu$ 的迭代轨迹
#| fig-showtext: true

stan_trace(eight_schools_fit, pars = "mu") +
  labs(x = "迭代次数", y = expression(mu))
```

### 后验分布

可以用函数 `stan_hist()` 或 `stan_dens()` 绘制后验分布图。下图分别展示参数 $\mu$、$\tau$ 的直方图，以及二者的散点图，参数 $\mu$ 的后验概率密度分布图。

```{r}
#| label: fig-8schools-rstan-posterior
#| fig-cap: rstan 包绘制后验分布图
#| fig-showtext: true
#| fig-height: 6

p1 <- stan_hist(eight_schools_fit, pars = c("mu","tau"), bins = 30)
p2 <- stan_scat(eight_schools_fit, pars = c("mu","tau"), size = 1) +
  labs(x = expression(mu), y = expression(tau))
p3 <- stan_dens(eight_schools_fit, pars = "mu") + labs(x = expression(mu))
library(patchwork)
p1 / (p2 + p3)
```

相比于 **rstan** 包，**bayesplot** 包可视化能力更强，支持对特定的参数做变换。**bayesplot** 包的函数 `mcmc_pairs()` 以矩阵图展示多个参数的分布，下图展示参数 $\mu$，$\log(\tau)$ 后验分布图。但是，这些函数都固定了一些标题，不能修改。

```{r}
#| label: fig-8schools-bayesplot-posterior
#| fig-cap: bayesplot 包绘制后验分布图
#| fig-showtext: true
#| fig-height: 6

bayesplot::mcmc_pairs(
  eight_schools_fit, pars = c("mu", "tau"), transform = list(tau = "log")
)
```

## 其它 R 包 {#sec-8schools-others}

### nlme

接下来，用 **nlme** 包拟合模型。

```{r}
# 成绩
y <- c(28, 8, -3, 7, -1, 1, 18, 12)
# 标准差
sigma <- c(15, 10, 16, 11, 9, 11, 10, 18)
# 学校编号
g <- 1:8
```

首先，调用 **nlme** 包的函数 `lme()` 拟合模型。

```{r}
library(nlme)
fit_lme <- lme(y ~ 1, random = ~ 1 | g, weights = varFixed(~ sigma^2), method = "REML")
summary(fit_lme)
```

随机效应的标准差 2.917988 ，随机效应部分的估计

```{r}
ranef(fit_lme)
```

类比 Stan 输出结果中的 $\theta$ 向量，每个学校的成绩估计

```{r}
7.785729 + 2.917988 * ranef(fit_lme)
```

### lme4

接着，采用 **lme4** 包拟合模型，发现 **lme4** 包获得与 **nlme** 包一样的结果。

```{r}
control <- lme4::lmerControl(
  check.conv.singular = "ignore",
  check.nobs.vs.nRE = "ignore",
  check.nobs.vs.nlev = "ignore"
)
fit_lme4 <- lme4::lmer(y ~ 1 + (1 | g), weights = 1 / sigma^2, control = control, REML = TRUE)
summary(fit_lme4)
```

### blme

下面使用 **blme** 包 [@Chung2013] ，**blme** 包基于 **lme4** 包，参数估计结果完全一致。

```{r}
# the mode should be at the boundary of the space.

fit_blme <- blme::blmer(
  y ~ 1 + (1 | g), control = control, REML = TRUE, 
  cov.prior = NULL, weights = 1 / sigma^2
)
summary(fit_blme)
```

### MCMCglmm

**MCMCglmm** 包 [@Hadfield2010] 采用 MCMC 算法拟合数据。

```{r}
schools <- data.frame(y = y, sigma = sigma, g = g)
schools$g <- as.factor(schools$g)
# inverse-gamma prior with scale and shape equal to 0.001
prior1 <- list(
  R = list(V = diag(schools$sigma^2), fix = 1),
  G = list(G1 = list(V = 1, nu = 0.002))
)
# 为可重复
set.seed(20232023)
# 拟合模型
fit_mcmc <- MCMCglmm::MCMCglmm(
  y ~ 1, random = ~g, rcov = ~ idh(g):units, 
  data = schools, prior = prior1, verbose = FALSE
)
# 输出结果
summary(fit_mcmc)
```

R-structure 表示残差方差，这是已知的参数。G-structure 表示随机截距的方差，Location effects 表示固定效应的截距。截距和 **nlme** 包的结果很接近。

### cmdstanr

一般地，**rstan** 包使用的 stan 框架版本低于 **cmdstanr** 包，从 **rstan** 包切换到 **cmdstanr** 包，需要注意语法、函数的变化。**rstan** 和 **cmdstanr** 使用的 Stan 版本不同导致参数估计结果不同，结果可重复的条件非常苛刻，详见 [Stan 参考手册](https://mc-stan.org/docs/reference-manual/reproducibility.html)。在都是较新的版本时，Stan 代码不需要做改动，如下：

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

此处，给参数 $\mu,\tau$ 添加了非常弱（模糊）的先验，结果将出现较大不同。

```{r}
#| message: false

eight_schools_dat <- list(
  J = 8,
  y = c(28, 8, -3, 7, -1, 1, 18, 12),
  sigma = c(15, 10, 16, 11, 9, 11, 10, 18)
)
library(cmdstanr)
mod_eight_schools <- cmdstan_model(
  stan_file = "code/eight_schools.stan",
  compile = TRUE, cpp_options = list(stan_threads = TRUE)
)
fit_eight_schools <- mod_eight_schools$sample(
  data = eight_schools_dat, # 数据
  chains = 2,            # 总链条数
  parallel_chains = 2,   # 并行数目
  iter_warmup = 1000,    # 每条链预处理的迭代次数
  iter_sampling = 1000,  # 每条链采样的迭代次数
  threads_per_chain = 2, # 每条链设置 2 个线程
  seed = 20232023,       # 随机数种子
  show_messages = FALSE, # 不显示消息
  refresh = 0 # 不显示采样迭代的进度
)
```

结果保留 3 位有效数字，模型输出如下：

```{r}
fit_eight_schools$summary(.num_args = list(sigfig = 3, notation = "dec"))
```

模型采样过程的诊断结果如下：

```{r}
fit_eight_schools$diagnostic_summary()
```

分层模型的参数 $\mu,\log(\tau)$ 的后验联合分布呈现经典的漏斗状。

```{r}
#| label: fig-8schools-funnels
#| fig-cap: 参数 $\mu,\log(\tau)$ 的联合分布
#| fig-width: 5
#| fig-height: 4
#| dev: 'tikz'
#| fig-process: !expr to_png

bayesplot::mcmc_scatter(
  fit_eight_schools$draws(), pars = c("mu", "tau"), 
  transform = list(tau = "log"), size = 2
) + labs(x = "$\\mu$", y = "$\\log(\\tau)$")
```

```{r}
#| eval: false
#| echo: false

bayesplot::mcmc_pairs(
  fit_eight_schools$draws(), pars = c("mu", "tau"), 
  transform = list(tau = "log")
)

eight_schools_df <- fit_eight_schools$draws(c("mu", "tau"), format = "draws_df")
ggplot(data = eight_schools_df, aes(x = mu, y = log(tau))) +
  geom_point(color = "#4285f4") +
  geom_density_2d(color = "#FBBC05", linewidth = 1) +
  theme_bw() +
  labs(x = expression(mu), y = expression(log(tau)))
```

对于调用 **cmdstanr** 包拟合的模型，适合用 **bayesplot** 包来可视化后验分布和诊断采样。

## 案例：rats 数据 {#sec-thirty-rats}

rats 数据最早来自 @gelfand1990 ，记录 30 只小鼠每隔一周的重量，一共进行了 5 周。第一次记录是小鼠第 8 天的时候，第二次测量记录是第 15 天的时候，一直持续到第 36 天。下面在 R 环境中准备数据。

```{r}
# 总共 30 只老鼠
N <- 30
# 总共进行 5 周
T <- 5
# 小鼠重量
y <- structure(c(
  151, 145, 147, 155, 135, 159, 141, 159, 177, 134,
  160, 143, 154, 171, 163, 160, 142, 156, 157, 152, 154, 139, 146,
  157, 132, 160, 169, 157, 137, 153, 199, 199, 214, 200, 188, 210,
  189, 201, 236, 182, 208, 188, 200, 221, 216, 207, 187, 203, 212,
  203, 205, 190, 191, 211, 185, 207, 216, 205, 180, 200, 246, 249,
  263, 237, 230, 252, 231, 248, 285, 220, 261, 220, 244, 270, 242,
  248, 234, 243, 259, 246, 253, 225, 229, 250, 237, 257, 261, 248,
  219, 244, 283, 293, 312, 272, 280, 298, 275, 297, 350, 260, 313,
  273, 289, 326, 281, 288, 280, 283, 307, 286, 298, 267, 272, 285,
  286, 303, 295, 289, 258, 286, 320, 354, 328, 297, 323, 331, 305,
  338, 376, 296, 352, 314, 325, 358, 312, 324, 316, 317, 336, 321,
  334, 302, 302, 323, 331, 345, 333, 316, 291, 324
), .Dim = c(30, 5))
# 第几天
x <- c(8.0, 15.0, 22.0, 29.0, 36.0)
xbar <- 22.0
```

重复测量的小鼠重量数据 rats 如下 @tbl-rats 所示。

```{r}
#| label: tbl-rats
#| tbl-cap: 小鼠重量数据（部分）
#| echo: false

rownames(y) <- 1:30
knitr::kable(head(y), col.names = paste("第", c(8, 15, 22, 29, 36), "天"), row.names = TRUE)
```

小鼠重量数据的分布和变化情况见下图，由图可以假定 30 只小鼠的重量服从正态分布，而30 只小鼠的重量呈现一种线性增长趋势。

```{r}
#| label: fig-rats
#| fig-cap: 30 只小鼠 5 次测量的数据
#| fig-subcap: 
#| - 小鼠重量的分布
#| - 小鼠重量的变化
#| fig-showtext: true
#| par: true
#| echo: false
#| fig-width: 5
#| fig-height: 4.5
#| layout-ncol: 2

matplot(y, xlab = "小鼠编号", ylab = "小鼠重量")
matplot(t(y), xlab = "测量次数", ylab = "小鼠重量", pch = 1)
```

## 频率派方法 {#sec-rats-frequentist}

### nlme {#sec-rats-nlme}

**nlme** 包适合长格式的数据，因此，先将小鼠数据整理成长格式。

```{r}
rats_data <- data.frame(
  weight = as.vector(y), 
  rats = rep(1:30, times = 5), 
  days = rep(c(8, 15, 22, 29, 36), each = 30)
)
```

将 30 只小鼠的重量变化及回归曲线画出来，发现各只小鼠的回归线的斜率几乎一样，截距略有不同。不同小鼠的出生重量是不同，前面 Stan 采用变截距变斜率的混合效应模型拟合数据。

```{r}
#| label: fig-rats-lm
#| fig-cap: 小鼠重量变化曲线
#| fig-showtext: true
#| fig-width: 6
#| fig-height: 7

ggplot(data = rats_data, aes(x = days, y = weight)) +
  geom_point() +
  geom_smooth(formula = "y ~ x", method = "lm", se = FALSE) +
  theme_bw() +
  facet_wrap(facets = ~rats, labeller = "label_both", ncol = 6) +
  labs(x = "第几天", y = "重量")
```

小鼠的重量随时间增长，不同小鼠的情况又会有所不同。作为一个参照，首先考虑变截距的随机效应模型。

$$
y_{ij} = \beta_0 + \beta_1 * x_j + \alpha_i + \epsilon_{ij}, \quad i = 1,2,\ldots,30. \quad j = 1,2,3,4,5
$$

其中，$y_{ij}$ 表示第 $i$ 只小鼠在第 $j$ 次测量的重量，一共 30 只小鼠，共测量了 5 次。固定效应部分是 $\beta_0$ 和 $\beta_1$ ，分别表示截距和斜率。随机效应部分是 $\alpha_i$ 和 $\epsilon_{ij}$ ，分别服从正态分布$\alpha_i \sim \mathcal{N}(0, \sigma^2_{\alpha})$ 和 $\epsilon_{ij} \sim \mathcal{N}(0, \sigma^2_{\epsilon})$ 。$\sigma^2_{\alpha}$ 和 $\sigma^2_{\epsilon}$ 分别表示组间方差（group level）和组内方差（individual level）。

```{r}
library(nlme)
rats_lme0 <- lme(data = rats_data, fixed = weight ~ days, random = ~ 1 | rats)
summary(rats_lme0)
```

当然，若考虑不同小鼠的生长速度不同（变化不是很大），可用变截距和变斜率的随机效应模型表示生长曲线模型，下面加载 **nlme** 包调用函数 `lme()` 拟合该模型。

```{r}
library(nlme)
rats_lme <- lme(data = rats_data, fixed = weight ~ days, random = ~ days | rats)
summary(rats_lme)
```

模型输出结果中，固定效应中的截距项 `(Intercept)` 对应 106.56762，斜率 `days` 对应 6.18571。Stan 模型中截距参数 `alpha0` 的后验估计是 106.332，斜率参数 `beta_c` 的后验估计是 6.188。对比 Stan 和 **nlme** 包的拟合结果，可以发现贝叶斯和频率方法的结果是非常接近的。截距参数 `alpha0` 可以看作小鼠的初始（出生）重量，斜率参数 `beta_c` 可以看作小鼠的生长率 growth rate。

函数 `lme()` 的输出结果中，随机效应的随机截距标准差 10.7425835，对应 `tau_alpha`，表示每个小鼠的截距偏移量的波动。而随机斜率的标准差为 0.5105447，对应 `tau_beta`，相对随机截距标准差来说很小。残差标准差为 6.0146608，对应 `tau_c`，表示与小鼠无关的剩余量的波动，比如测量误差。总之，和 Stan 的结果有所不同，但相去不远。主要是前面的 Stan 模型没有考虑随机截距和随机斜率之间的相关性，这可以进一步调整 [@sorensen2016] 。

```{r}
# 参数的置信区间
intervals(rats_lme, level = 0.95)
```

Stan 输出中，截距项 alpha、斜率项 beta 参数的标准差分别是 `tau_alpha` 和 `tau_beta` ，残差标准差参数 `tau_c` 的估计为 6.1。简单起见，没有考虑截距项和斜率项的相关性，即不考虑小鼠出生时的重量和生长率的相关性，一般来说，应该是有关系的。函数 `lme()` 的输出结果中给出了截距项和斜率项的相关性为 -0.343，随机截距和随机斜率的相关性为 -0.159。

计算与 Stan 输出中的截距项 `alpha_c` 对应的量，结合函数 `lme()` 的输出，截距、斜率加和之后，如下

```{r}
106.56762 + 6.18571 * 22
```

值得注意，Stan 代码中对时间 days 做了中心化处理，即 $x_t - \bar{x}$，目的是降低采样时参数 $\alpha_i$ 和 $\beta_i$ 之间的相关性，而在拟合函数 `lme()` 中没有做处理，因此，结果无需转化，而且更容易解释。

```{r}
fit_lm <- lm(weight ~ days, data = rats_data)
summary(fit_lm)
```

采用简单线性模型即可获得与 **nlme** 包非常接近的估计结果，主要是小鼠重量的分布比较正态，且随时间的变化非常线性。

### lavaan

**lavaan** 包 [@Rosseel2012] 主要是用来拟合结构方程模型，而生长曲线模型可以放在该框架下。所以，也可以用 **lavaan** 包来拟合，并且，它提供的函数 `growth()` 可以直接拟合生长曲线模型。

```{r}
#| message: false

library(lavaan)
# 设置矩阵 y 的列名
colnames(y) <- c("t1","t2","t3","t4","t5")
rats_growt_model <- " 
  # intercept and slope with fixed coefficients
  intercept =~ 1*t1 + 1*t2 + 1*t3 + 1*t4 + 1*t5
  days =~ 0*t1 + 1*t2 + 2*t3 + 3*t4 + 4*t5 

  # if we fix the variances to be equal, the models are now identical.
  t1 ~~ resvar*t1    
  t2 ~~ resvar*t2
  t3 ~~ resvar*t3
  t4 ~~ resvar*t4
  t5 ~~ resvar*t5
"
```

其中，算子符号 `=~` 定义潜变量，`~~` 定义残差协方差，intercept 表示截距， days 表示斜率。假定 5 次测量的测量误差（组内方差）是相同的。拟合模型的代码如下：

```{r}
rats_growth_fit <- growth(rats_growt_model, data = y)
```

提供函数 `summary()` 获得模型输出，结果如下：

```{r}
summary(rats_growth_fit, fit.measures = TRUE)
```

输出结果显示 **lavaan** 包的函数 `growth()` 采用极大似然估计方法。协方差部分 `Covariances:` 随机效应中斜率和截距的协方差。截距部分 `Intercepts:` 对应于混合效应模型的固定效应部分。方差部分 `Variances:` 对应于混合效应模型的随机效应部分，包括残差方差、斜率和截距的方差。不难看出，这和前面 **nlme** 包的输出结果差别很大。原因是 **lavaan** 包将测量的次序从 0 开始计，0 代表小鼠出生后的第 8 天。也就是说，**lavaan** 采用的是次序标记，而不是实际数据。将测量发生的时间（第几天）换算成次序（第几次），并从 0 开始计，则函数 `lme()` 的输出和函数 `growth()` 就一致了。

```{r}
# 重新组织数据
rats_data2 <- data.frame(
  weight = as.vector(y), 
  rats = rep(1:30, times = 5), 
  days = rep(c(0, 1, 2, 3, 4), each = 30)
)
# ML 方法估计模型参数
rats_lme2 <- lme(data = rats_data2, fixed = weight ~ days, random = ~ days | rats, method = "ML")
summary(rats_lme2)
```

可以看到函数 `growth()` 给出的截距和斜率的协方差估计为 8.444，函数 `lme()` 给出对应截距和斜率的标准差分别是 10.652390 和 3.496588，它们的相关系数为 0.227，则函数 `lme()` 给出的协方差估计为 `10.652390*3.496588*0.227` ，即 8.455，协方差估计比较一致。同理，比较两个输出结果中的其它成分，函数 `growth()` 给出的残差方差估计为 36.176，则残差标准差估计为 6.0146，结合函数 `lme()` 给出的 `Random effects:` 中 `Residual`，结果完全一样。函数 `growth()` 给出的 `Intercepts:` 对应于函数 `lme()` 给出的固定效应部分，结果也是完全一样。

针对模型拟合对象 `rats_growth_fit` ，除了函数 `summary()` 可以汇总结果，**lavaan** 包还提供 `AIC()` 、 `BIC()` 和 `logLik()` 等函数，分别可以提取 AIC、BIC 和对数似然值， `AIC()` 和 `logLik()` 结果与前面的函数 `lme()` 的输出是一样的，而 `BIC()` 不同。

### lme4

当采用 **lme4** 包拟合数据的时候，发现输出结果与 **nlme** 包几乎相同。

```{r}
rats_lme4 <- lme4::lmer(weight ~ days + (days | rats), data = rats_data)
summary(rats_lme4)
```

### glmmTMB

glmmTMB 包基于 Template Model Builder (TMB) ，拟合广义线性混合效应模型，公式语法与 **lme4** 包一致。

```{r}
#| message: false

rats_glmmtmb <- glmmTMB::glmmTMB(weight ~ days + (days | rats), REML = TRUE, data = rats_data)
summary(rats_glmmtmb)
```

结果与 **nlme** 包完全一样。

### MASS

**MASS** 包的结果与前面完全一致。

```{r}
rats_mass <- MASS::glmmPQL(
  fixed = weight ~ days, random = ~ days | rats, 
  data = rats_data, family = gaussian(), verbose = FALSE
)
summary(rats_mass)
```

### spaMM

**spaMM** 包的结果与前面完全一致。

```{r}
#| message: false

rats_spamm <- spaMM::fitme(weight ~ days + (days | rats), data = rats_data)
summary(rats_spamm)
```

``` markdown
 --------------- Random effects ---------------
Family: gaussian( link = identity ) 
         --- Random-coefficients Cov matrices:
 Group        Term   Var.   Corr.
  rats (Intercept)  110.1        
  rats        days 0.2495 -0.1507
# of obs: 150; # of groups: rats, 30 
```

随机效应的截距方差 110.1，斜率方差 0.2495，则标准差分别是 10.49 和 0.499，相关性为 -0.1507。

``` markdown
 -------------- Residual variance  ------------
phi estimate was 36.1755 
```

残差方差为 36.1755，则标准差为 6.0146。

### hglm

**hglm** 包 [@rönnegård2010] 可以拟合分层广义线性模型，线性混合效应模型和广义线性混合效应模型，随机效应和响应变量服从的分布可以很广泛，使用语法与 **lme4** 包一样。

```{r}
rats_hglm <- hglm::hglm2(weight ~ days + (days | rats), data = rats_data)
summary(rats_hglm)
```

固定效应的截距和斜率都是和 **nlme** 包的输出结果一致。值得注意，随机效应和模型残差都是以发散参数（Dispersion parameter）来表示的，模型残差方差为 37.09572，则标准差为 6.0906，随机效应的随机截距和随机斜率的方差分别为 103.4501 和 0.2407，则标准差分别为 10.1710 和 0.4906，这与 **nlme** 包的结果也是一致的。

### mgcv

先考虑一个变截距的混合效应模型

$$
y_{ij} = \beta_0 + \beta_1 * x_j + \alpha_i + \epsilon_{ij}, \quad i = 1,2,\ldots,30. \quad j = 1,2,3,4,5
$$

假设随机效应服从独立同正态分布，等价于在似然函数中添加一个岭惩罚。广义可加模型在一定形式下和上述混合效应模型存在等价关系，在广义可加模型中，可以样条表示随机效应。**mgcv** 包拟合代码如下。

```{r}
#| message: false

library(mgcv)
rats_data$rats <- as.factor(rats_data$rats)
rats_gam <- gam(weight ~ days + s(rats, bs = "re"), data = rats_data)
```

其中，参数取值 `bs = "re"` 指定样条类型，re 是 Random effects 的简写。

```{r}
summary(rats_gam)
```

其中，残差的方差 Scale est. = 67.303 ，则标准差为 $\sigma_{\epsilon} = 8.2038$ 。随机效应的标准差如下

```{r}
gam.vcomp(rats_gam, rescale = TRUE)
```

`rescale = TRUE` 表示恢复至原数据的尺度，标准差 $\sigma_{\alpha} = 14.033$。可以看到，固定效应和随机效应的估计结果与 **nlme** 包等完全一致。若考虑变截距和变斜率的混合效应模型，拟合代码如下：

```{r}
rats_gam1 <- gam(
  weight ~ days + s(rats, bs = "re") + s(rats, by = days, bs = "re"),
  data = rats_data, method = "REML"
)
summary(rats_gam1)
```

输出结果中，固定效应部分的结果和 **nlme** 包完全一样。

```{r}
gam.vcomp(rats_gam1, rescale = TRUE)
```

输出结果中，依次是随机效应的截距、斜率和残差的标准差（标准偏差），和 **nlme** 包给出的结果非常接近。

**mgcv** 包还提供函数 `gamm()`，它将混合效应和固定效应分开，在拟合 LMM 模型时，它类似 **nlme** 包的函数 `lme()`。返回一个含有 lme 和 gam 两个元素的列表，前者包含随机效应的估计，后者是固定效应的估计，固定效应中可以添加样条（或样条表示的简单随机效益，比如本节前面提及的模型）。实际上，函数 `gamm()` 分别调用 **nlme** 包和 **MASS** 包来拟合 LMM 模型和 GLMM 模型。

```{r}
rats_gamm <- gamm(weight ~ days, random = list(rats = ~days), method = "REML", data = rats_data)
# LME
summary(rats_gamm$lme)
# GAM
summary(rats_gamm$gam)
```

## 贝叶斯方法 {#sec-rats-bayesianism}

### rstan {#sec-rats-rstan}

初始化模型参数，设置采样算法的参数。

```{r}
# 迭代链
chains <- 4
# 迭代次数
iter <- 1000
# 初始值
init <- rep(list(list(
  alpha = rep(250, 30), beta = rep(6, 30),
  alpha_c = 150, beta_c = 10,
  tausq_c = 1, tausq_alpha = 1,
  tausq_beta = 1
)), chains)
```

接下来，基于重复测量数据，建立线性生长曲线模型：

$$
\begin{aligned}
\alpha_c &\sim \mathcal{N}(0,100) \quad \beta_c  \sim \mathcal{N}(0,100) \\
\tau^2_{\alpha} &\sim \mathrm{inv\_gamma}(0.001, 0.001) \\
\tau^2_{\beta}  &\sim \mathrm{inv\_gamma}(0.001, 0.001) \\
\tau^2_c &\sim \mathrm{inv\_gamma}(0.001, 0.001) \\
\alpha_n &\sim \mathcal{N}(\alpha_c, \tau_{\alpha})  \quad
\beta_n  \sim \mathcal{N}(\beta_c, \tau_{\beta}) \\
y_{nt} &\sim \mathcal{N}(\alpha_n + \beta_n * (x_t - \bar{x}), \tau_c) \\
& n = 1,2,\ldots,N \quad t = 1,2,\ldots,T
\end{aligned}
$$

其中， $\alpha_c,\beta_c,\tau_c,\tau_{\alpha},\tau_{\beta}$ 为无信息先验，$\bar{x} = 22$ 表示第 22 天，$N = 30$ 和 $T = 5$ 分别表示实验中的小鼠数量和测量次数，下面采用 Stan 编码、编译、采样和拟合模型。

```{r}
rats_fit <- stan(
  model_name = "rats",
  model_code = "
  data {
    int<lower=0> N;
    int<lower=0> T;
    vector[T] x;
    matrix[N,T] y;
    real xbar;
  }
  parameters {
    vector[N] alpha;
    vector[N] beta;

    real alpha_c;
    real beta_c;          // beta.c in original bugs model

    real<lower=0> tausq_c;
    real<lower=0> tausq_alpha;
    real<lower=0> tausq_beta;
  }
  transformed parameters {
    real<lower=0> tau_c;       // sigma in original bugs model
    real<lower=0> tau_alpha;
    real<lower=0> tau_beta;

    tau_c = sqrt(tausq_c);
    tau_alpha = sqrt(tausq_alpha);
    tau_beta = sqrt(tausq_beta);
  }
  model {
    alpha_c ~ normal(0, 100);
    beta_c ~ normal(0, 100);
    tausq_c ~ inv_gamma(0.001, 0.001);
    tausq_alpha ~ inv_gamma(0.001, 0.001);
    tausq_beta ~ inv_gamma(0.001, 0.001);
    alpha ~ normal(alpha_c, tau_alpha); // vectorized
    beta ~ normal(beta_c, tau_beta);  // vectorized
    for (n in 1:N)
      for (t in 1:T)
        y[n,t] ~ normal(alpha[n] + beta[n] * (x[t] - xbar), tau_c);
  }
  generated quantities {
    real alpha0;
    alpha0 = alpha_c - xbar * beta_c;
  }
  ",
  data = list(N = N, T = T, y = y, x = x, xbar = xbar),
  chains = chains, init = init, iter = iter,   
  verbose = FALSE, refresh = 0, seed = 20190425
)
```

模型输出结果如下：

```{r}
print(rats_fit, pars = c("alpha", "beta"), include = FALSE, digits = 1)
```

`alpha_c` 表示小鼠 5 次测量的平均重量，`beta_c` 表示小鼠体重的增长率，$\alpha_i,\beta_i$ 分别表示第 $i$ 只小鼠在第 22 天（第 3 次测量或 $x_t = \bar{x}$ ）的重量和增长率（每日增加的重量）。

对于分量众多的参数向量，比较适合用岭线图展示后验分布，下面调用 **bayesplot** 包绘制参数向量 $\boldsymbol{\alpha},\boldsymbol{\beta}$ 的后验分布。

```{r}
#| label: fig-rats-alpha
#| fig-cap: 参数 $\boldsymbol{\alpha}$ 的后验分布
#| fig-showtext: true
#| fig-width: 6
#| fig-height: 8
#| message: false

# plot(rats_fit, pars = "alpha", show_density = TRUE, ci_level = 0.8, outer_level = 0.95)
bayesplot::mcmc_areas_ridges(rats_fit, pars = paste0("alpha", "[", 1:30, "]")) +
  scale_y_discrete(labels = scales::parse_format()) 
```

参数向量 $\boldsymbol{\alpha}$ 的后验估计可以看作 $x_t = \bar{x}$ 时小鼠的重量，上图即为各个小鼠重量的后验分布。

```{r}
#| label: fig-rats-beta
#| fig-cap: 参数 $\boldsymbol{\beta}$ 的后验分布
#| fig-showtext: true
#| fig-width: 6
#| fig-height: 8
#| message: false

# plot(rats_fit, pars = "beta", ci_level = 0.8, outer_level = 0.95)
bayesplot::mcmc_areas_ridges(rats_fit, pars = paste0("beta", "[", 1:30, "]")) +
  scale_y_discrete(labels = scales::parse_format()) 
```

参数向量 $\boldsymbol{\beta}$ 的后验估计可以看作是小鼠的重量的增长率，上图即为各个小鼠重量的增长率的后验分布。

### cmdstanr

从 rstan 包转 cmdstanr 包是非常容易的，只要语法兼容，模型代码可以原封不动。

```{r}
#| message: false

library(cmdstanr)
mod_rats <- cmdstan_model(
  stan_file = "code/rats.stan",
  compile = TRUE, cpp_options = list(stan_threads = TRUE)
)
fit_rats <- mod_rats$sample(
  data = list(N = N, T = T, y = y, x = x, xbar = xbar), # 数据
  chains = 2,            # 总链条数
  parallel_chains = 2,   # 并行数目
  iter_warmup = 1000,    # 每条链预处理的迭代次数
  iter_sampling = 1000,  # 每条链采样的迭代次数
  threads_per_chain = 2, # 每条链设置 2 个线程
  seed = 20232023,       # 随机数种子
  show_messages = FALSE, # 不显示消息
  adapt_delta = 0.9,     # 接受率
  refresh = 0 # 不显示采样迭代的进度
)
```

模型输出

```{r}
# 显示除了参数 alpha 和 beta 以外的结果
vars <- setdiff(fit_rats$metadata()$stan_variables, c("alpha", "beta"))
fit_rats$summary(variables = vars)
```

诊断信息

```{r}
fit_rats$diagnostic_summary()
```

### brms

**brms** 包是基于 **rstan** 包的，基于 Stan 语言做贝叶斯推断，提供与 lme4 包一致的公式语法，且扩展了模型种类。

```{r}
#| eval: false

rats_brms <- brms::brm(weight ~ days + (days | rats), data = rats_data)
summary(rats_brms)
```

``` markdown
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: weight ~ days + (days | rats) 
   Data: rats_data (Number of observations: 150) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~rats (Number of levels: 30) 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)          11.27      2.23     7.36    16.08 1.00     2172     2939
sd(days)                0.54      0.09     0.37     0.74 1.00     1380     2356
cor(Intercept,days)    -0.11      0.24    -0.53     0.39 1.00      920     1541

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   106.47      2.47   101.61   111.23 1.00     2173     2768
days          6.18      0.11     5.96     6.41 1.00     1617     2177

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     6.15      0.47     5.30     7.14 1.00     1832     3151

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

### rstanarm

**rstanarm** 包与 **brms** 包类似，区别是前者预编译了 Stan 模型，后者根据输入数据和模型编译即时编译，此外，后者支持的模型范围更加广泛。

```{r}
#| eval: false

library(rstanarm)
rats_rstanarm <- stan_lmer(formula = weight ~ days + (days | rats), data = rats_data)
summary(rats_rstanarm)
```

``` markdown
Model Info:
 function:     stan_lmer
 family:       gaussian [identity]
 formula:      weight ~ days + (days | rats)
 algorithm:    sampling
 sample:       4000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 150
 groups:       rats (30)

Estimates:
                                      mean    sd      10%     50%     90%  
(Intercept)                         106.575   2.236 103.789 106.559 109.415
days                                  6.187   0.111   6.048   6.185   6.329
sigma                                 6.219   0.497   5.626   6.183   6.862
Sigma[rats:(Intercept),(Intercept)] 103.927  42.705  57.329  98.128 159.086
Sigma[rats:days,(Intercept)]         -0.545   1.492  -2.361  -0.402   1.162
Sigma[rats:days,days]                 0.304   0.112   0.181   0.285   0.445

MCMC diagnostics
                                    mcse  Rhat  n_eff
(Intercept)                         0.043 1.000 2753 
days                                0.003 1.005 1694 
sigma                               0.015 1.001 1172 
Sigma[rats:(Intercept),(Intercept)] 1.140 1.000 1403 
Sigma[rats:days,(Intercept)]        0.054 1.006  772 
Sigma[rats:days,days]               0.003 1.000 1456 

For each parameter, mcse is Monte Carlo standard error, 
n_eff is a crude measure of effective sample size, 
and Rhat is the potential scale reduction factor 
on split chains (at convergence Rhat=1).
```

固定效应的部分，截距和斜率如下：

``` markdown
Estimates:
                                      mean    sd      10%     50%     90%  
(Intercept)                         106.575   2.236 103.789 106.559 109.415
days                                  6.187   0.111   6.048   6.185   6.329
```

模型残差的标准差 sigma、随机效应 Sigma 的随机截距的方差 103.927 、随机斜率的方差 0.304 及其协方差 -0.545。

``` markdown
sigma                                 6.219   0.497   5.626   6.183   6.862
Sigma[rats:(Intercept),(Intercept)] 103.927  42.705  57.329  98.128 159.086
Sigma[rats:days,(Intercept)]         -0.545   1.492  -2.361  -0.402   1.162
Sigma[rats:days,days]                 0.304   0.112   0.181   0.285   0.445
```

**rstanarm** 和 **brms** 包的结果基本一致的。

### blme

**blme** 包 [@Chung2013] 基于 **lme4** 包 [@Bates2015] 拟合贝叶斯线性混合效应模型。参考前面 **rstan** 小节中关于模型参数的先验设置，下面将残差方差的先验设置为逆伽马分布，随机效应的协方差设置为扁平分布。发现拟合结果和 **nlme** 和 **lme4** 包的几乎一样。

```{r}
rats_blme <- blme::blmer(
  weight ~ days + (days | rats), data = rats_data,
  resid.prior = invgamma, cov.prior = NULL
)
summary(rats_blme)
```

与 **lme4** 包的函数 `lmer()` 所不同的是参数 `resid.prior` 、`fixef.prior` 和 `cov.prior` ，它们设置参数的先验分布，其它参数的含义同 `lme4` 包。`resid.prior = invgamma` 表示残差方差参数使用逆伽马分布，`cov.prior = NULL` 表示随机效应的协方差参数使用扁平先验 flat priors。

### rjags

**rjags** [@rjags] 是 JAGS 软件的 R 语言接口，可以拟合分层正态模型，再借助 **coda 包** [@coda2006] 可以分析 JAGS 返回的各项数据。

JAGS 代码和 Stan 代码有不少相似之处，最大的共同点在于以直观的统计模型的符号表示编码模型，仿照 Stan 代码， JAGS 编码的模型（BUGS 代码）如下：

```{verbatim, file="code/rats.bugs", lang="bugs"}
```

转化主要集中在模型块，注意二者概率分布的名称以及参数含义对应关系，JAGS 使用 precision 而不是 standard deviation or variance，比如正态分布中的方差（标准偏差）被替换为其倒数。JAGS 可以省略类型声明（初始化模型时会补上），最后，JAGS 不支持 Stan 中的向量化操作，这种新特性是独特的。

```{r}
#| message: false

library(rjags)
# 初始值
rats_inits <- list(
  list(".RNG.name" = "base::Marsaglia-Multicarry", 
       ".RNG.seed" = 20222022, 
       "alpha_c" = 100, "beta_c" = 6, "tau_c" = 5, "tau_alpha" = 10, "tau_beta" = 0.5),
  list(".RNG.name" = "base::Marsaglia-Multicarry", 
       ".RNG.seed" = 20232023, 
       "alpha_c" = 200, "beta_c" = 10, "tau_c" = 15, "tau_alpha" = 15, "tau_beta" = 1)
)
# 模型
rats_model <- jags.model(
  file = "code/rats.bugs",
  data = list(x = x, y = y, N = 30, T = 5, xbar = 22.0),
  inits = rats_inits, 
  n.chains = 2, quiet = TRUE
)
# burn-in
update(rats_model, n.iter = 2000)
# 抽样
rats_samples <- coda.samples(rats_model,
  variable.names = c("alpha_c", "beta_c", "sigma_alpha", "sigma_beta", "sigma_c"),
  n.iter = 4000, thin = 1
)
# 参数的后验估计
summary(rats_samples)
```

输出结果与 rstan 十分一致，且采样速度极快。类似地，`alpha0 = alpha_c - xbar * beta_c` 可得 alpha0 = 242.4752 - 22 \* 6.1878 = 106.3436。

### MCMCglmm

同前，先考虑变截距的混合效应模型，**MCMCglmm** 包 [@Hadfield2010] 给出的拟合结果与 **nlme** 包很接近。

```{r}
## 变截距模型
prior1 <- list(
  R = list(V = 1, nu = 0.002),
  G = list(G1 = list(V = 1, nu = 0.002))
)
set.seed(20232023)
rats_mcmc1 <- MCMCglmm::MCMCglmm(
  weight ~ days, random = ~ rats,
  data = rats_data, verbose = FALSE, prior = prior1
)
summary(rats_mcmc1)
```

随机效应的方差（组间方差）为 211.4 ，则标准差为 14.539。残差方差（组内方差）为 68.77，则标准差为 8.293。

再考虑变截距和斜率的混合效应模型。

```{r}
## 变截距、变斜率模型
prior2 <- list(
  R = list(V = 1, nu = 0.002),
  G = list(G1 = list(V = diag(2), nu = 0.002))
)
set.seed(20232023)
rats_mcmc2 <- MCMCglmm::MCMCglmm(weight ~ days,
  random = ~ us(1 + days):rats,
  data = rats_data, verbose = FALSE, prior = prior2
)
summary(rats_mcmc2)
```

G-structure 代表随机效应部分，R-structure 代表残差效应部分，Location effects 代表固定效应部分。**MCMCglmm** 包的这套模型表示术语源自商业软件 [ASReml](https://vsni.co.uk/software/asreml) 。

随机截距的方差为 124.1327，标准差为 11.1415，随机斜率的方差 0.2783，标准差为 0.5275，随机截距和随机斜率的协方差 -0.7457，相关系数为 -0.1268，这与 **nlme** 包结果很接近。

### INLA

同前，先考虑变截距的混合效应模型。

```{r}
#| message: false

library(INLA)
inla.setOption(short.summary = TRUE)
# 数值稳定性考虑
rats_data$weight <- rats_data$weight / 400
# 变截距
rats_inla1 <- inla(weight ~ days + f(rats, model = "iid", n = 30), 
                  family = "gaussian", data = rats_data)
# 输出结果
summary(rats_inla1)
```

再考虑变截距和斜率的混合效应模型。

```{r}
# https://inla.r-inla-download.org/r-inla.org/doc/latent/iid.pdf
# 二维高斯随机效应的先验为 Wishart prior
rats_data$rats <- as.integer(rats_data$rats)
rats_data$slopeid <- 30 + rats_data$rats
# 变截距、变斜率
rats_inla2 <- inla(
  weight ~ 1 + days + f(rats, model = "iid2d", n = 2 * 30) + f(slopeid, days, copy = "rats"),
  data = rats_data, family = "gaussian"
)
# 输出结果
summary(rats_inla2)
```

::: callout-warning
对于变截距和斜率混合效应模型，还未完全弄清楚 INLA 包的输出结果。固定效应部分和残差部分都是和前面一致的，但不清楚随机效应的方差协方差矩阵的估计与 INLA 输出的对应关系。参考[《Bayesian inference with INLA》](https://becarioprecario.bitbucket.io/inla-gitbook/index.html)第 3 章第 3 小节。
:::

## 总结 {#sec-hierarchical-normal-models-summary}

基于 rats 数据建立变截距、变斜率的分层正态模型，也是线性混合效应模型的一种特殊情况，下表给出不同方法对模型各个参数的估计及置信区间。

|                | $\beta_0$ | $\beta_1$ | $\sigma_0$ | $\sigma_1$ | $\rho_{\sigma}$ | $\sigma_{\epsilon}$ |
|-----------------|-----------|-----------|-----------|-----------|-----------|-----------|
| nlme (REML)    | 106.568   | 6.186     | 10.743     | 0.511      | -0.159          | 6.015               |
| lme4 (REML)    | 106.568   | 6.186     | 10.744     | 0.511      | -0.16           | 6.015               |
| glmmTMB (REML) | 106.568   | 6.186     | 10.743     | 0.511      | -0.16           | 6.015               |
| MASS (PQL)     | 106.568   | 6.186     | 10.495     | 0.500      | -0.15           | 6.015               |
| spaMM (ML)     | 106.568   | 6.186     | 10.49      | 0.499      | -0.15           | 6.015               |
| hglm           | 106.568   | 6.186     | 10.171     | 0.491      | \-              | 6.091               |
| mgcv (REML)    | 106.568   | 6.186     | 10.311     | 0.492      | \-              | 6.069               |

: 频率派方法比较 {#tbl-rats-freqentist-compare}

表中给出截距 $\beta_0$ 、斜率 $\beta_1$ 、随机截距 $\sigma_0$、随机斜率 $\sigma_1$、随机截距和斜率的相关系数 $\rho_{\sigma}$、残差 $\sigma_{\epsilon}$ 等参数的估计及 95% 的置信区间，四舍五入保留 3 位小数。固定效应部分的结果完全相同，随机效应部分略有不同。

|                 | $\beta_0$ | $\beta_1$ | $\sigma_0$ | $\sigma_1$ | $\rho_{\sigma}$ | $\sigma_{\epsilon}$ |
|-----------------|-----------|-----------|-----------|-----------|-----------|-----------|
| rstan (NUTS)    | 106.4     | 6.2       | 14.6       | 0.5        | \-              | 6.1                 |
| cmdstanr (NUTS) | 106       | 6.19      | 14.5       | 0.513      | \-              | 6.09                |
| brms (NUTS)     | 106.47    | 6.18      | 11.27      | 0.54       | -0.11           | 6.15                |
| rstanarm (NUTS) | 106.575   | 6.187     | 10.194     | 0.551      | -0.0969         | 6.219               |
| blme (REML)     | 106.568   | 6.186     | 10.787     | 0.512      | -0.160          | 5.949               |
| rjags (Gibbs)   | 106.344   | 6.188     | 14.623     | 0.518      | \-              | 6.073               |
| MCMCglmm (MCMC) | 106.40    | 6.19      | 11.14      | 0.53       | -0.13           | 6.18                |

: 贝叶斯方法比较 {#tbl-rats-bayesian-compare}

其中，**INLA** 结果的转化未完成，表格中暂缺。**rstan** 、 **cmdstanr** 和 **rjags** 未考虑随机截距和随机斜率的相关性，因此，相关系数暂缺。MCMC 是一种随机优化算法，在不同的实现中，可重复性的要求不同，设置随机数种子仅是其中的一个必要条件，故而，每次运行程序结果可能略微不同，但不影响结论。Stan 相关的 R 包输出结果中，**rstan** 保留 1 位小数，**cmdstanr** 保留 3 位有效数字，**brms** 保留 2 位小数，**rstanarm** 小数点后保留 3 位有效数字，各不相同，暂未统一处理。

## 习题 {#sec-hierarchical-models-exercises}

1.  四个组的重复测量数据，如下表所示，建立贝叶斯线性混合效应模型/分层正态模型分析数据，与 nlme 包拟合的结果对比。

    ```{r}
    #| label: tbl-exer
    #| tbl-cap: 实验数据
    #| echo: false

    y <- c(
      62, 60, 63, 59,
      63, 67, 71, 64, 65, 66,
      68, 66, 71, 67, 68, 68,
      56, 62, 60, 61, 63, 64, 63, 59
    )
    group <- c(rep(1, 4), rep(2, 6), rep(3, 6), rep(4, 8))
    id <- c(1:4, 1:6, 1:6, 1:8)
    dat <- data.frame(y = y, group = group, id = id)
    dat2 <- reshape(dat, direction = "wide", timevar = "group", idvar = "id")
    options(knitr.kable.NA = '')
    knitr::kable(dat2, col.names = c("编号", "第1组", "第2组", "第3组", "第4组"), row.names = FALSE)
    ```

    $$
    \begin{aligned}
    y_{ij}   \sim \mathcal{N}(\theta_i, \sigma^2) &\quad
    \theta_i \sim \mathcal{N}(\mu, \tau^2) \\
    (\mu,\log \sigma, \tau) &\sim \mathrm{uniform\ prior} \\
    i = 1,2,3,4 &\quad j = 1,2, \ldots, n_i
    \end{aligned}
    $$

    $y_{ij}$ 表示第 $i$ 组的第 $j$ 个测量值，$\theta_i$ 表示第 $i$ 组的均值，$\mu$ 表示整体的均值，$\sigma^2$ 表示组内的方差，$\tau^2$ 表示组内的方差。

    ```{r}
    library(nlme)
    fit_lme <- lme(data = dat, fixed = y ~ 1, random = ~ 1 | group)
    summary(fit_lme)
    ```

    随机效应（组间标准差）$\tau^2$ 3.419288 、残差效应（组内标准差）$\sigma^2$ 2.366309。截距 $\mu$ 64.01266 代表整体的均值。各组的均值如下：

    ```{r}
    64.01266 + ranef(fit_lme)
    ```

    也可以调用 **rjags** 包连接 JAGS 软件做贝叶斯推理，JAGS 代码如下：

    ```{verbatim, file="code/hnm.bugs", lang="bugs"}
    ```

    完整的运行代码如下：

    ```{r}
    #| message: false

    library(rjags)
    # 参考值
    mu_a <- min(y)
    mu_b <- max(y)
    log_sigma_b <- 2 * log(sd(y))
    tau_b <- 2 * sd(y)

    J <- 4            # 4 个组
    n <- length(y)    # 观察值数量
    N <- 1500         # 总采样数
    nthin <- 1        # 采样间隔
    nchains <- 2      # 2 条链
    ndiscard <- N / 2 # 预处理阶段 warm-up / burn-in

    # 初始值
    jags_inits <- list(
      list(".RNG.name" = "base::Marsaglia-Multicarry", 
           ".RNG.seed" = 20222022, 
           "theta" = rep(3, 4), "mu" = 60, "log_sigma" = 0, "tau" = 1.5),
      list(".RNG.name" = "base::Marsaglia-Multicarry", 
           ".RNG.seed" = 20232023, 
           "theta" = rep(2, 4), "mu" = 60, "log_sigma" = 1, "tau" = 0.375)
    )
    # Call JAGS from R
    jags_model <- jags.model(
      file = "code/hnm.bugs",
      data = list("y" = y, "group" = group, "J" = J, "n" = n),
      inits = jags_inits, n.chains = nchains, quiet = TRUE
    )
    # burn-in
    update(jags_model, n.iter = ndiscard)
    # 抽样
    jags_samples <- coda.samples(jags_model,
      variable.names = c('theta','mu','sigma','tau'), n.iter = N
    )
    # 参数的后验估计
    summary(jags_samples)
    ```

2.  基于 **lme4** 包中学生对老师的评价数据 `InstEval` 建立（广义）线性混合效应模型分析数据。将响应变量（学生评价）视为有序的离散型变量，比较观察两个模型拟合效果（lme4、GLMMadaptive、spaMM 都不支持有序的响应变量，brms 则支持各类有序回归，使用语法与 lme4 完全一样。但是，由于数据规模比较大，计算时间数以天计，可考虑用 Stan 直接编码）。再者，从 Stan 实现的贝叶斯模型来看，感受 Stan 建模的灵活性和扩展性。（**nlme** 包不支持此等交叉随机效应的表达。）

    ```{r}
    data(InstEval, package = "lme4")
    str(InstEval)
    ```

    -   因子型变量 `s` 表示 1-2972 位参与评分的学生。
    -   因子型变量 `d` 表示 1-2160 位上课的讲师。
    -   因子型变量 `dept` 表示课程相关的 1-15 院系。
    -   因子型变量 `service` 表示讲师除了授课外，是否承担其它服务。
    -   数值型变量 `y` 表示学生给课程的评分，1-5 分对应从坏到很好。

    ```{r}
    # 数值型的响应变量
    fit_lme4 <- lme4::lmer(y ~ 1 + service + (1 | s) + (1 | d) + (1 | dept), data = InstEval)
    summary(fit_lme4)
    ```

    **lme4** 包不支持响应变量为有序分类变量的情形，可用 **ordinal** 包，此等规模数据，拟合模型需要 5-10 分钟时间。

    ```{r}
    #| eval: false
    #| echo: true

    # 有序因子型的响应变量
    InstEval$y <- factor(InstEval$y, ordered = TRUE)
    library(ordinal)
    fit_ordinal <- clmm(
      y ~ 1 + service + (1 | s) + (1 | d) + (1 | dept),
      data = InstEval, link = "probit", threshold = "equidistant"
    )
    summary(fit_ordinal)

    ## MCMCglmm
    library(MCMCglmm)
    prior2 <- list(
      R = list(V = 1, nu = 0.002),
      G = list(
        G1 = list(V = 1, nu = 0.002),
        G2 = list(V = 1, nu = 0.002),
        G3 = list(V = 1, nu = 0.002)
      )
    )
    # 响应变量视为数值变量
    fit_mcmc2 <- MCMCglmm(
      y ~ service, random = ~ s + d + dept, family = "gaussian",
      data = InstEval, verbose = FALSE, prior = prior2
    )
    # 响应变量视为有序的分类变量
    fit_mcmc3 <- MCMCglmm(
      y ~ service, random = ~ s + d + dept, family = "ordinal",
      data = InstEval, verbose = FALSE, prior = prior2
    )
    ```

    当数据量较大时，**MCMCglmm** 包拟合模型需要很长时间，放弃，此时，Stan 的相对优势可以体现出来了。Stan 适合大型复杂概率统计模型。