1. 分层正态模型

数据集来自 Gelfand et al (1990) 1,研究 30 只小鼠的体重变化,每周测量其体重,连续测量了5周。用$Y_{ij}$表示第$i$只鼠在时间$x_j$测量的体重。部分数据集见下表,其中 Week 1 表示小鼠一周大 $x_1 = 8$,Week 2 表示小鼠两周大 $x_2 = 15$ ,依次类推。

N <- 30
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 # x 的平均值
# 第一张表格
rownames(y) <- paste("Rat", seq(30))
knitr::kable(y,
  format = "markdown", caption = "Rats 数据集",
  row.names = TRUE,
  col.names = paste("Week", seq(5))
)
Week 1 Week 2 Week 3 Week 4 Week 5
Rat 1 151 199 246 283 320
Rat 2 145 199 249 293 354
Rat 3 147 214 263 312 328
Rat 28 157 205 248 289 316
Rat 29 137 180 219 258 291
Rat 30 153 200 244 286 324

为了比较直观地观察数据集,我们绘制如下图形

# 第一张图
png(filename = "rats.png", res = 200, width = 480*200/72, height = 480*200/72,
    type = "cairo", bg = "transparent")
par(mar = c(4, 4, 0.5, 0.5))
matplot(t(y),
  col = hcl.colors(30), type = "b",
  xlab = "Weeks", ylab = "Weight", lty = 2, pch = 1
)
dev.off()

rats

根据图中的显示,我们建立带随机效应的线性模型

$$ Y_{ij} \sim \mathsf{Normal}(\alpha_i + \beta_i(X_{j} - X_{bar}), \tau_{c})\\ \alpha_i \sim \mathsf{Normal}(\alpha_{c}, \tau_{\alpha}) \\ \beta_j \sim \mathsf{Normal}(\beta_{c},\tau_{\beta}) \\ i = 1, 2, \ldots, 30;\quad j = 1, 2, \ldots, 5. $$

根据的先验选择,我们设置的先验如下

$$ \alpha_{c} \sim \mathsf{Normal}(0, 100) \\ \beta_{c} \sim \mathsf{Normal}(0, 100) \\ \tau^2_{c} \sim \mathsf{inv\_{gamma}}(0.001, 0.001) \\ \tau^2_{\alpha} \sim \mathsf{inv\_{gamma}}(0.001, 0.001) \\ \tau^2_{\beta} \sim \mathsf{inv\_{gamma}}(0.001, 0.001). $$

  1. 为什么选择这些先验?后验分布的具体形式能推导吗?参数的似然函数能推导吗?2 见安德鲁的书籍 BDA 第三版
  2. 它们是如何表示成 Stan 代码中的模型块?见 Stan 文档
  3. Stan 是如何编码计算这些统计模型的,自动微分,变分推断?见 Stan 文档
  4. 分层正态模型的 JAGS nimble Stan 实现,参考 bridgesampling 包的文档,该 R 包基于 Xiao-Li Meng and Wing Hung Wong (1996) Simulating ratios of normallizing constants via A simple identity: a theoretical exploration http://www3.stat.sinica.edu.tw/statistica/j6n4/j6n43/j6n43.htm 提出的 Bridge Sampling 比较贝叶斯模型,还可以计算边际似然、贝叶斯因子、后验模型概率,bfrmsBayesFactor 可用于计算贝叶斯因子,前者适用范围更广

加载 rstan 包后,我们需要设置一些模型运行时和模型结果展示阶段的参数

# BUGS rats example (Vol 1, Example 1)
# http://www.openbugs.net/Examples/Rats.html
library(rstan, quietly = TRUE)
# 将编译后的模型写入磁盘,避免重新编译
rstan_options(auto_write = TRUE)
# 如果CPU和内存足够,设置成与马尔科夫链一样多
options(mc.cores = 2)
# Windows 上可能需要指定
# Sys.setenv(USE_Cxx14 = 1)
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
)

后续以Stan编写的模型中,每个参数迭代的链条数为4,即从后验分布中抽样的次数为4,每次抽样的数量为 1000

chains = 4 # 链条数
iter = 1000 # 样本量

接着设置模型参数初始值,每条链开始迭代的初值值都一样,这里选用 WinBUGS 给的初始值,后续我们还考虑随机选择的初始值与之比较

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)

调用 rstan 包中的关键函数 stan 编译运行 Stan 编写的分层正态模型

file
指定模型文件的路径
model_name
Stan 代码编译成功后,指定的模型名称
model_code
传递一段字符串,字符串的内容即是 Stan 代码编写的模型
data
参数传递模型的观测数据集,注意数据结构是 R 语言的列表 list,其每个元素的属性需和 Stan 代码编写的模型中设定的参数属性保持一致
algorithm
Stan 可选的抽样算法有 c("NUTS", "HMC", "Fixed_param") 默认和推荐的参数值是 "NUTS",这三者的算法细节及使用情况以后介绍。
chains
参数迭代时设置的马氏链条数目
init
模型中各个参数的初始值,也是用列表来存储,列表中的元素属性和 Stan 代码中的参数也要保持一致,参数名称也要一致,默认值 random 表示随机选择
iter
模型每条链迭代的次数,在这里设置了四条链,迭代总计4000次,抽样间隔默认是1,预处理阶段
thin
抽样间隔,即每迭代多少次抽一组迭代值,默认是1,就是迭代多少次就抽取多少次
warmup
预处理阶段迭代的次数,默认设置 floor(iter/2) 就是每条链迭代总次数的一半,向下取整
seed
模型运行时设置的随机数种子,Stan 内置的算法都是基于 MCMC 的,为了复现结果,设置相同的随机数种子有时是必要的
cores
设置模型运行时可使用的系统 CPU 的核心数,在内存和 CPU 资源充足的情况下,越多越好。注意,可用的 CPU 核心,可用 parallel::detectCores() 检测出来,不要超过这个,此外,数目也不要超过参数 chains 设置的数值,超过这个意义也不大。其默认值是 getOption("mc.cores", 1L) 即获取R语言环境中全局参数列表中的 mc.cores,若无设置就取1
verbose
是否显示模型编译和运行过程中的状态,默认是 FALSE 即不显示中间的迭代过程
control
传递一个列表,包含算法的细节调控,调整采样的策略,在采样效率和准确率中寻找一个平衡,此处不详细介绍,留待详细介绍的算法细节的时候介绍
fit
参数传递一个类型为 S4 的 stanfit 对象实例,默认是 NA,如果参数值不是 NA 而是前一次相关模型拟合的结果,那么模型拟合就会重用之前的结果,可以节省下重编译 Stan 代码的时间
boost_lib = NULL, eigen_lib = NULL
Stan 语言大量使用了 boost 和 eigen 两个 C++ 库,模型编译过程允许使用用户指定的版本
save_dso
参数默认值是 TRUE,该参数只有在参数 fit = NA 时才有作用。意义是从 C++ 编译的动态分享对象 DSO (dynamic shared object) 是否需要保留。如果保留,那么从另一个 R 控制台运行模型时,不再重复编译 C++ 代码,而是直接使用保留的 DSO,这在并行的时候是可以节省编译时间的。
...
... 这种参数在 R 语言当中是比较高级的传参方式了,其实际是传递可变长度的参数向量。其中,参数 refresh 控制采样进度的报告频率,默认值是 max(iter/10, 1),取值 0 表示不显示采样的进度;save_warmup 是否保留预处理阶段的采样过程,默认参数值 TRUE 表示保留。还有参数chain_idinit_rtest_grad 等,这里不一一介绍了,其它参数可见帮助文档 ?rstan::stan
rats_fit <- stan(
  # file = "code/bugs/vol1/rats/rats.stan", 
  model_name = "rats",
  model_code = "
  data {
    int<lower=0> N;
    int<lower=0> T;
    real x[T];
    real y[N,T];
    real xbar;
  }
  parameters {
    real alpha[N];
    real beta[N];
  
    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)

分层正态模型 Stan 拟合结果

print(rats_fit, digits = 1)
Inference for Stan model: rats.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.

              mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
alpha[1]     240.0     0.1  2.7  234.7  238.2  240.0  241.8  245.1  2379    1
alpha[2]     247.7     0.1  2.7  242.3  245.9  247.7  249.4  253.2  2293    1
alpha[3]     252.4     0.1  2.8  246.7  250.6  252.4  254.3  257.8  2296    1
alpha[4]     232.6     0.1  2.7  227.4  230.9  232.7  234.5  237.7  2553    1
alpha[5]     231.6     0.1  2.7  226.2  229.8  231.6  233.4  237.0  2443    1
alpha[6]     249.8     0.1  2.6  244.7  248.0  249.7  251.6  254.7  2171    1
alpha[7]     228.7     0.1  2.6  223.5  226.9  228.6  230.4  234.0  1980    1
alpha[8]     248.3     0.1  2.7  242.7  246.6  248.2  250.0  253.7  2162    1
alpha[9]     283.3     0.1  2.7  277.9  281.5  283.4  285.2  288.5  2060    1
alpha[10]    219.2     0.0  2.5  214.3  217.5  219.2  220.9  224.2  2708    1
alpha[11]    258.3     0.0  2.7  253.1  256.5  258.2  260.1  263.7  2995    1
alpha[12]    228.1     0.1  2.6  223.0  226.3  228.1  229.8  233.3  1985    1
alpha[13]    242.5     0.1  2.7  237.4  240.7  242.4  244.2  247.9  1926    1
alpha[14]    268.3     0.1  2.7  262.9  266.5  268.3  270.1  273.4  2332    1
alpha[15]    242.7     0.1  2.8  237.4  240.8  242.7  244.6  248.3  2225    1
alpha[16]    245.3     0.1  2.6  240.3  243.6  245.3  247.0  250.4  2365    1
alpha[17]    232.2     0.1  2.6  226.8  230.5  232.2  233.9  237.2  2353    1
alpha[18]    240.5     0.1  2.7  235.2  238.7  240.5  242.3  245.7  2023    1
alpha[19]    253.8     0.1  2.7  248.5  252.0  253.8  255.6  258.8  1920    1
alpha[20]    241.6     0.1  2.7  236.2  239.8  241.7  243.5  246.9  2394    1
alpha[21]    248.5     0.1  2.6  243.2  246.8  248.6  250.3  253.5  2314    1
alpha[22]    225.3     0.1  2.7  220.1  223.4  225.2  227.2  230.7  1805    1
alpha[23]    228.5     0.1  2.7  223.6  226.7  228.4  230.2  233.7  2334    1
alpha[24]    245.1     0.1  2.7  239.7  243.3  245.0  246.7  250.5  2420    1
alpha[25]    234.5     0.1  2.8  229.1  232.6  234.6  236.3  239.9  1869    1
alpha[26]    253.9     0.1  2.8  248.5  252.0  254.0  255.8  259.4  1879    1
alpha[27]    254.3     0.1  2.6  249.4  252.5  254.3  256.0  259.3  2400    1
alpha[28]    242.9     0.1  2.7  237.8  241.1  242.9  244.7  248.3  2405    1
alpha[29]    217.9     0.1  2.7  212.5  216.0  217.9  219.6  223.2  2335    1
alpha[30]    241.4     0.1  2.8  235.9  239.6  241.4  243.1  246.9  2132    1
beta[1]        6.1     0.0  0.2    5.6    5.9    6.1    6.2    6.5  2285    1
beta[2]        7.1     0.0  0.2    6.6    6.9    7.0    7.2    7.5  1882    1
beta[3]        6.5     0.0  0.3    6.0    6.3    6.5    6.6    7.0  2247    1
beta[4]        5.3     0.0  0.3    4.8    5.2    5.3    5.5    5.9  1729    1
beta[5]        6.6     0.0  0.2    6.1    6.4    6.6    6.7    7.0  2280    1
beta[6]        6.2     0.0  0.2    5.7    6.0    6.2    6.3    6.7  2230    1
beta[7]        6.0     0.0  0.2    5.5    5.8    6.0    6.1    6.4  2666    1
beta[8]        6.4     0.0  0.2    5.9    6.3    6.4    6.6    6.9  1885    1
beta[9]        7.1     0.0  0.3    6.5    6.9    7.0    7.2    7.6  2122    1
beta[10]       5.8     0.0  0.2    5.4    5.7    5.8    6.0    6.3  2389    1
beta[11]       6.8     0.0  0.2    6.3    6.6    6.8    7.0    7.3  2002    1
beta[12]       6.1     0.0  0.2    5.6    6.0    6.1    6.3    6.6  2128    1
beta[13]       6.2     0.0  0.2    5.7    6.0    6.2    6.3    6.6  2698    1
beta[14]       6.7     0.0  0.2    6.2    6.5    6.7    6.9    7.2  1970    1
beta[15]       5.4     0.0  0.3    4.9    5.2    5.4    5.6    5.9  1975    1
beta[16]       5.9     0.0  0.2    5.4    5.8    5.9    6.1    6.4  2187    1
beta[17]       6.3     0.0  0.2    5.8    6.1    6.3    6.4    6.7  2442    1
beta[18]       5.8     0.0  0.2    5.4    5.7    5.8    6.0    6.3  2299    1
beta[19]       6.4     0.0  0.2    5.9    6.2    6.4    6.6    6.9  2126    1
beta[20]       6.1     0.0  0.2    5.6    5.9    6.1    6.2    6.5  2023    1
beta[21]       6.4     0.0  0.2    5.9    6.2    6.4    6.6    6.9  1801    1
beta[22]       5.9     0.0  0.2    5.4    5.7    5.9    6.0    6.3  2136    1
beta[23]       5.7     0.0  0.3    5.3    5.6    5.8    5.9    6.2  1862    1
beta[24]       5.9     0.0  0.2    5.4    5.7    5.9    6.1    6.4  2427    1
beta[25]       6.9     0.0  0.3    6.4    6.7    6.9    7.1    7.4  2053    1
beta[26]       6.5     0.0  0.2    6.1    6.4    6.5    6.7    7.0  2407    1
beta[27]       5.9     0.0  0.2    5.4    5.7    5.9    6.1    6.4  2528    1
beta[28]       5.8     0.0  0.2    5.4    5.7    5.8    6.0    6.3  2228    1
beta[29]       5.7     0.0  0.2    5.2    5.5    5.7    5.8    6.2  2038    1
beta[30]       6.1     0.0  0.2    5.7    6.0    6.1    6.3    6.6  2378    1
alpha_c      242.5     0.1  2.8  236.4  240.6  242.5  244.3  248.0  1936    1
beta_c         6.2     0.0  0.1    6.0    6.1    6.2    6.3    6.4  1565    1
tausq_c       37.3     0.2  5.7   27.8   33.2   36.8   41.0   49.4  1178    1
tausq_alpha  220.2     1.6 66.6  125.7  173.9  207.8  251.5  378.8  1812    1
tausq_beta     0.3     0.0  0.1    0.1    0.2    0.3    0.3    0.5  1207    1
tau_c          6.1     0.0  0.5    5.3    5.8    6.1    6.4    7.0  1171    1
tau_alpha     14.7     0.0  2.1   11.2   13.2   14.4   15.9   19.5  1919    1
tau_beta       0.5     0.0  0.1    0.4    0.5    0.5    0.6    0.7  1161    1
alpha0       106.5     0.1  3.7   98.9  104.1  106.6  109.0  113.6  1540    1
lp__        -438.4     0.3  6.7 -452.3 -442.8 -438.1 -433.7 -426.0   567    1

Samples were drawn using NUTS(diag_e) at Mon May 13 18:53:46 2019.
For each parameter, 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).

用图表示参数向量 $\alpha,\beta$ 各个分量的后验分布

rats-alpha rats-beta

# 保存图形
png(filename = "rats-alpha.png", res = 200, width = 480*200/72, height = 480*200/72,
    type = "cairo", bg = "transparent")
plot(rats_fit, pars="alpha")
dev.off()
png(filename = "rats-beta.png", res = 200, width = 480*200/72, height = 480*200/72,
    type = "cairo", bg = "transparent")
plot(rats_fit, pars="beta")
dev.off()

各个参数的迭代轨迹图全部画出,实在太多,但是作为重要的环节,一定要逐一验证,每个参数的迭代过程都是平稳的,只有基于平稳的马氏链才可以做出合理的后验估计。下面以参数 $\alpha_1$ 为例

rats-traceplot-alpha1

png(filename = "rats-traceplot-alpha1.png", res = 200, width = 480*200/72, height = 480*200/72,
    type = "cairo", bg = "transparent")
traceplot(rats_fit, pars = c("alpha[1]"), inc_warmup = FALSE) +
  theme_minimal() +
  labs(x = "Iteration", y = expression(alpha[1]))
dev.off()

考虑使用随机生成的初始值,为了节省版面,常常将 Stan 代码保存到磁盘文件,这里就不妨将前面的 model_code 参数值保存到文件 rats.stan

# 使用随机生成的初始值
rats <- stan(
  #file = "code/bugs/vol1/rats/rats.stan",
  model_code = "
  data {
    int<lower=0> N;
    int<lower=0> T;
    real x[T];
    real y[N,T];
    real xbar;
  }
  parameters {
    real alpha[N];
    real beta[N];
  
    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, iter = iter, init = "random"
)

Stan 模型在如前所述的设置中,运行上面的代码块是非常快的,因为之前编译的模型没有变,只是改变了算法迭代的初始值,可以直接使用

print(rats, digits = 1)
Inference for Stan model: rats.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.

              mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
alpha[1]     239.9     0.1  2.7  234.7  238.1  239.8  241.6  245.3  1628    1
alpha[2]     247.7     0.1  2.6  242.6  245.9  247.7  249.4  252.6  2227    1
alpha[3]     252.4     0.1  2.7  247.0  250.6  252.3  254.2  257.7  2812    1
alpha[4]     232.6     0.1  2.7  227.5  230.7  232.5  234.3  237.9  2284    1
alpha[5]     231.7     0.1  2.6  226.5  229.8  231.7  233.5  236.8  2229    1
alpha[6]     249.7     0.1  2.7  244.6  247.9  249.7  251.5  255.0  2098    1
alpha[7]     228.7     0.1  2.7  223.1  226.9  228.7  230.5  234.0  1874    1
alpha[8]     248.4     0.1  2.7  243.2  246.6  248.4  250.2  253.6  2736    1
alpha[9]     283.3     0.1  2.7  278.0  281.6  283.4  285.2  288.3  2257    1
alpha[10]    219.3     0.1  2.7  214.1  217.5  219.3  221.1  224.4  2117    1
alpha[11]    258.2     0.1  2.8  252.6  256.4  258.1  260.0  263.8  1893    1
alpha[12]    228.1     0.1  2.7  222.8  226.4  228.1  229.9  233.6  1975    1
alpha[13]    242.4     0.1  2.6  237.4  240.8  242.4  244.1  247.6  2328    1
alpha[14]    268.2     0.1  2.6  263.4  266.5  268.2  270.0  273.2  2181    1
alpha[15]    242.8     0.1  2.7  237.6  240.9  242.8  244.7  248.2  1262    1
alpha[16]    245.3     0.1  2.7  240.1  243.4  245.3  247.1  250.5  2137    1
alpha[17]    232.2     0.1  2.7  227.0  230.4  232.2  234.1  237.4  2480    1
alpha[18]    240.3     0.1  2.6  234.9  238.7  240.3  242.1  245.5  2476    1
alpha[19]    253.8     0.1  2.6  248.6  252.0  253.8  255.5  259.0  2621    1
alpha[20]    241.5     0.1  2.7  236.3  239.7  241.6  243.4  247.0  2086    1
alpha[21]    248.6     0.1  2.7  243.2  246.8  248.5  250.4  253.9  1916    1
alpha[22]    225.2     0.1  2.6  219.9  223.5  225.2  226.9  230.3  2170    1
alpha[23]    228.6     0.1  2.5  223.7  226.9  228.6  230.2  233.6  2365    1
alpha[24]    245.1     0.1  2.6  240.0  243.3  245.1  246.8  249.9  2269    1
alpha[25]    234.5     0.1  2.7  229.4  232.6  234.4  236.3  240.0  2272    1
alpha[26]    253.9     0.1  2.7  248.6  252.0  253.9  255.6  259.4  2608    1
alpha[27]    254.4     0.1  2.7  249.2  252.6  254.3  256.1  259.8  2714    1
alpha[28]    242.9     0.1  2.6  237.8  241.2  242.9  244.7  248.0  2290    1
alpha[29]    217.9     0.1  2.7  212.8  216.1  217.9  219.8  223.4  1949    1
alpha[30]    241.4     0.1  2.8  235.8  239.5  241.4  243.3  246.7  2101    1
beta[1]        6.1     0.0  0.2    5.6    5.9    6.1    6.2    6.5  2438    1
beta[2]        7.1     0.0  0.3    6.6    6.9    7.0    7.2    7.5  1905    1
beta[3]        6.5     0.0  0.2    6.0    6.3    6.5    6.6    6.9  2662    1
beta[4]        5.3     0.0  0.3    4.8    5.2    5.3    5.5    5.8  1703    1
beta[5]        6.6     0.0  0.2    6.1    6.4    6.6    6.7    7.0  2734    1
beta[6]        6.2     0.0  0.2    5.7    6.0    6.2    6.3    6.6  2530    1
beta[7]        6.0     0.0  0.2    5.5    5.8    6.0    6.1    6.5  2470    1
beta[8]        6.4     0.0  0.2    5.9    6.3    6.4    6.6    6.9  2131    1
beta[9]        7.0     0.0  0.3    6.5    6.9    7.0    7.2    7.5  2072    1
beta[10]       5.8     0.0  0.2    5.4    5.7    5.8    6.0    6.3  2436    1
beta[11]       6.8     0.0  0.2    6.3    6.6    6.8    7.0    7.3  2030    1
beta[12]       6.1     0.0  0.2    5.7    6.0    6.1    6.3    6.6  2384    1
beta[13]       6.2     0.0  0.2    5.7    6.0    6.2    6.3    6.6  2385    1
beta[14]       6.7     0.0  0.2    6.2    6.5    6.7    6.9    7.2  2611    1
beta[15]       5.4     0.0  0.3    4.9    5.3    5.4    5.6    5.9  2452    1
beta[16]       5.9     0.0  0.2    5.4    5.8    5.9    6.1    6.4  2055    1
beta[17]       6.3     0.0  0.2    5.8    6.1    6.3    6.4    6.7  2314    1
beta[18]       5.8     0.0  0.2    5.4    5.7    5.8    6.0    6.3  2368    1
beta[19]       6.4     0.0  0.3    5.9    6.2    6.4    6.6    6.9  2431    1
beta[20]       6.1     0.0  0.2    5.6    5.9    6.1    6.2    6.5  2612    1
beta[21]       6.4     0.0  0.2    5.9    6.2    6.4    6.6    6.9  2163    1
beta[22]       5.9     0.0  0.2    5.4    5.7    5.9    6.0    6.4  2000    1
beta[23]       5.7     0.0  0.2    5.3    5.6    5.7    5.9    6.2  2699    1
beta[24]       5.9     0.0  0.2    5.4    5.7    5.9    6.1    6.4  3197    1
beta[25]       6.9     0.0  0.2    6.4    6.7    6.9    7.1    7.4  2084    1
beta[26]       6.5     0.0  0.2    6.1    6.4    6.5    6.7    7.1  2199    1
beta[27]       5.9     0.0  0.2    5.4    5.7    5.9    6.1    6.4  3171    1
beta[28]       5.9     0.0  0.2    5.4    5.7    5.9    6.0    6.3  2536    1
beta[29]       5.7     0.0  0.2    5.2    5.5    5.7    5.8    6.2  2137    1
beta[30]       6.1     0.0  0.2    5.6    6.0    6.1    6.3    6.6  2700    1
alpha_c      242.5     0.1  2.7  237.3  240.8  242.4  244.3  247.7  1826    1
beta_c         6.2     0.0  0.1    6.0    6.1    6.2    6.3    6.4  1860    1
tausq_c       37.1     0.2  5.7   27.6   33.0   36.6   40.5   50.1   955    1
tausq_alpha  218.4     2.1 71.2  122.8  169.9  206.5  251.2  391.7  1166    1
tausq_beta     0.3     0.0  0.1    0.1    0.2    0.3    0.3    0.5  1531    1
tau_c          6.1     0.0  0.5    5.3    5.7    6.0    6.4    7.1   954    1
tau_alpha     14.6     0.1  2.2   11.1   13.0   14.4   15.8   19.8  1389    1
tau_beta       0.5     0.0  0.1    0.4    0.4    0.5    0.6    0.7  1449    1
alpha0       106.3     0.1  3.6   98.9  104.0  106.3  108.6  113.4  1666    1
lp__        -438.0     0.3  7.3 -453.9 -442.7 -437.5 -432.9 -425.4   532    1

Samples were drawn using NUTS(diag_e) at Mon May 13 19:03:38 2019.
For each parameter, 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).

从上面的结果来看,随机选择的初始值对结果没有影响,基于 NUTS 采样器实现的算法对初值不敏感,Stan 也具有非常高的采样效率。

Stan 代码尽可能地使用内置的向量化函数,可以极大地加快编译速度,运行速度,达到更好的优化性能。

# 使用随机生成的初始值
rats_vec <- stan(
  # file = "code/bugs/vol1/rats/rats_vec.stan", 
  model_code = "
  data {
    int<lower=0> N;
    int<lower=0> T;
    real x[T];
    real y[N,T];
    real xbar;
  }
  transformed data {
    real x_minus_xbar[T];
    real y_linear[N*T];
  
    for (t in 1:T)
      x_minus_xbar[t] = x[t] - xbar;
  
    for (n in 1:N)
      for (t in 1:T)
        y_linear[(n-1)*T + t] = y[n, t];
  }
  parameters {
    real alpha[N];
    real beta[N];
  
    real alpha_c;
    real beta_c;
  
    real<lower=0> tausq_c;
    real<lower=0> tausq_alpha;
    real<lower=0> tausq_beta;
  }
  transformed parameters {
    real<lower=0> tau_c;
    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 {
    real pred[N*T];
  
    for (n in 1:N)
      for (t in 1:T)
        pred[(n-1)*T + t] = fma(beta[n], x_minus_xbar[t], alpha[n]);
  
    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
    y_linear ~ normal(pred, tau_c);  // vectorized
  }",  
  data = list(N = N, T = T, y = y, x = x, xbar = xbar),
  chains = chains, iter = iter, init = "random"
)
print(rats_vec, digits = 1)
Inference for Stan model: 83c0b6cb8d216e4723b1edc83099db35.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.

              mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
alpha[1]     239.9     0.1  2.7  234.8  238.0  239.9  241.8  245.2  1278  1.0
alpha[2]     247.6     0.1  2.7  242.1  245.8  247.6  249.5  252.8   361  1.0
alpha[3]     252.2     0.3  3.1  244.3  250.4  252.4  254.2  258.0   124  1.0
alpha[4]     232.8     0.3  3.1  227.2  230.7  232.6  234.6  240.4   116  1.0
alpha[5]     232.0     0.3  3.1  226.6  230.0  231.7  233.7  240.7    95  1.0
alpha[6]     249.5     0.2  3.0  243.2  247.6  249.6  251.5  255.3   239  1.0
alpha[7]     229.2     0.4  3.4  223.7  227.1  228.9  230.9  239.8    64  1.0
alpha[8]     248.2     0.1  2.9  242.4  246.4  248.3  250.1  254.2   385  1.0
alpha[9]     282.0     1.3  6.9  246.0  281.1  283.1  284.8  288.1    27  1.1
alpha[10]    219.9     0.7  4.6  214.0  217.5  219.3  221.4  239.2    38  1.1
alpha[11]    257.8     0.5  3.7  245.2  256.1  258.1  260.1  263.5    63  1.1
alpha[12]    228.5     0.5  3.4  223.0  226.4  228.2  230.0  239.7    58  1.1
alpha[13]    242.4     0.1  2.8  236.9  240.5  242.4  244.3  247.6  1890  1.0
alpha[14]    267.6     0.9  4.9  246.1  266.1  268.2  270.1  273.6    32  1.1
alpha[15]    242.8     0.1  2.6  237.6  241.0  242.8  244.5  247.9  2598  1.0
alpha[16]    245.2     0.1  2.6  240.1  243.4  245.2  247.0  250.5  1877  1.0
alpha[17]    232.6     0.3  3.1  227.1  230.6  232.4  234.3  240.8    91  1.0
alpha[18]    240.4     0.1  2.8  234.9  238.6  240.4  242.3  245.8  1618  1.0
alpha[19]    253.6     0.3  3.3  245.2  251.8  253.8  255.6  259.3   118  1.0
alpha[20]    241.6     0.1  2.8  236.2  239.7  241.5  243.4  247.1  2047  1.0
alpha[21]    248.4     0.2  2.9  242.7  246.6  248.5  250.4  254.0   355  1.0
alpha[22]    225.7     0.5  3.8  220.0  223.6  225.4  227.1  240.3    60  1.1
alpha[23]    228.8     0.4  3.4  223.4  226.8  228.5  230.3  240.1    71  1.0
alpha[24]    245.1     0.1  2.7  239.7  243.3  245.1  246.9  250.4  1982  1.0
alpha[25]    234.7     0.3  3.1  229.1  232.6  234.5  236.5  241.7   143  1.0
alpha[26]    253.7     0.3  3.2  244.4  252.0  253.9  255.7  259.2    98  1.0
alpha[27]    254.0     0.3  3.3  244.3  252.5  254.2  256.0  259.7   106  1.0
alpha[28]    242.9     0.1  2.7  237.7  241.3  242.9  244.8  248.1  2227  1.0
alpha[29]    218.6     0.7  4.8  212.6  216.1  218.1  220.0  240.0    41  1.1
alpha[30]    241.5     0.1  2.7  236.0  239.6  241.5  243.3  246.8  1777  1.0
beta[1]        6.1     0.0  0.2    5.6    5.9    6.1    6.2    6.5   785  1.0
beta[2]        7.0     0.0  0.3    6.4    6.8    7.0    7.2    7.6    91  1.0
beta[3]        6.5     0.0  0.2    6.0    6.3    6.5    6.6    6.9  1648  1.0
beta[4]        5.4     0.0  0.3    4.8    5.2    5.4    5.5    6.2    73  1.0
beta[5]        6.6     0.0  0.2    6.1    6.4    6.6    6.7    7.0  1072  1.0
beta[6]        6.2     0.0  0.2    5.7    6.0    6.2    6.3    6.6  2051  1.0
beta[7]        6.0     0.0  0.2    5.5    5.8    6.0    6.1    6.5   466  1.0
beta[8]        6.4     0.0  0.2    6.0    6.2    6.4    6.6    6.9  1865  1.0
beta[9]        7.0     0.0  0.3    6.3    6.8    7.0    7.2    7.6   108  1.0
beta[10]       5.9     0.0  0.3    5.4    5.7    5.9    6.0    6.4   237  1.0
beta[11]       6.8     0.0  0.3    6.2    6.6    6.8    7.0    7.3   185  1.0
beta[12]       6.1     0.0  0.2    5.7    6.0    6.1    6.3    6.6  1571  1.0
beta[13]       6.2     0.0  0.2    5.7    6.0    6.2    6.3    6.7  2226  1.0
beta[14]       6.7     0.0  0.2    6.2    6.5    6.7    6.8    7.2   337  1.0
beta[15]       5.4     0.0  0.3    4.9    5.2    5.4    5.6    6.2    56  1.1
beta[16]       5.9     0.0  0.2    5.5    5.8    5.9    6.1    6.4   314  1.0
beta[17]       6.3     0.0  0.2    5.8    6.1    6.3    6.4    6.7  2943  1.0
beta[18]       5.9     0.0  0.2    5.4    5.7    5.9    6.0    6.4   236  1.0
beta[19]       6.4     0.0  0.2    5.9    6.2    6.4    6.6    6.9  2360  1.0
beta[20]       6.1     0.0  0.2    5.6    5.9    6.1    6.2    6.6   711  1.0
beta[21]       6.4     0.0  0.2    5.9    6.2    6.4    6.6    6.9  2454  1.0
beta[22]       5.9     0.0  0.3    5.4    5.7    5.9    6.0    6.4   253  1.0
beta[23]       5.8     0.0  0.3    5.3    5.6    5.8    5.9    6.3   183  1.0
beta[24]       5.9     0.0  0.3    5.4    5.7    5.9    6.1    6.4   318  1.0
beta[25]       6.9     0.0  0.3    6.3    6.7    6.9    7.1    7.4   133  1.0
beta[26]       6.5     0.0  0.2    6.1    6.4    6.5    6.7    7.0  1081  1.0
beta[27]       5.9     0.0  0.3    5.4    5.7    5.9    6.1    6.4   288  1.0
beta[28]       5.9     0.0  0.2    5.4    5.7    5.9    6.0    6.3   254  1.0
beta[29]       5.7     0.0  0.3    5.2    5.5    5.7    5.9    6.3   116  1.0
beta[30]       6.1     0.0  0.2    5.7    6.0    6.1    6.3    6.6  1706  1.0
alpha_c      242.4     0.1  2.6  237.0  240.8  242.4  244.1  247.7  2163  1.0
beta_c         6.2     0.0  0.1    6.0    6.1    6.2    6.3    6.4   555  1.0
tausq_c       43.8     7.6 37.2   27.7   33.5   37.1   41.6  198.1    24  1.1
tausq_alpha  209.0     6.9 68.9    5.5  170.5  203.1  245.2  353.5   100  1.0
tausq_beta     0.3     0.0  0.1    0.0    0.2    0.3    0.3    0.5   111  1.0
tau_c          6.4     0.3  1.7    5.3    5.8    6.1    6.4   14.1    25  1.1
tau_alpha     14.2     0.4  2.9    2.3   13.1   14.3   15.7   18.8    43  1.1
tau_beta       0.5     0.0  0.1    0.0    0.4    0.5    0.6    0.7    47  1.1
lp__        -437.9     0.6  8.0 -453.3 -442.5 -437.8 -433.5 -424.1   178  1.0

Samples were drawn using NUTS(diag_e) at Mon May 13 19:38:13 2019.
For each parameter, 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).

毫无疑问,向量化只是为了提高计算效率,是编程技巧,不会显著影响结果的准确度

2. 简单线性回归模型

小鼠的体重设为响应变量 Y,小鼠的年龄设为 X,不考虑个体之间的差异,也不按年龄分组,

$$ Y = X\beta + \epsilon $$

其中,$Y = (y_1,\ldots, y_{150})^\top, X = (\mathbf{1}, \mathbf{x}), \mathbf{1} = (1,1,\ldots,1)^\top, \mathbf{x} = (x_1, \ldots, x_{150})^\top,\beta = (\beta_1,\beta_2)^\top$,则

$$ \beta = (X^{T}X)^{-1}X^{T}Y \\ \epsilon = Y - \hat{Y} = Y - X(X^{T}X)^{-1}X^{T}Y = (I - X(X^{T}X)^{-1}X^{T})Y $$

# 将矩阵 y 拉直
rats_dat <- data.frame(weight = as.vector(y), day = rep(c(8, 15, 22, 29, 36), each = 30))
# 散点图
plot(data = rats_dat, weight ~ day)
# 响应变量
Y <- rats_dat$weight
# 系数矩阵
X <- cbind(rep(1, 150), rats_dat$day)

计算系数矩阵的条件数

base::rcond(t(X) %*% X)
[1] 0.0002686286
base::kappa(t(X) %*% X)
[1] 3593.169

注意到 t(X) %*% X 是可逆的,计算 $\beta$ 的最小二乘估计

# 求解线性方程组
beta <- solve(t(X) %*% X, t(X) %*%Y)
# 调用 lm 函数
fit <- lm(weight ~ day, data = rats_dat)
summary(fit)
Call:
lm(formula = weight ~ day, data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-38.253 -11.278   0.197   7.647  64.047 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 106.5676     3.2099   33.20   <2e-16 ***
day           6.1857     0.1331   46.49   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.13 on 148 degrees of freedom
Multiple R-squared:  0.9359,    Adjusted R-squared:  0.9355 
F-statistic:  2161 on 1 and 148 DF,  p-value: < 2.2e-16

下面,我们将线性模型拟合的结果呈现到图上

png(filename = "rats-lm.png", res = 200, width = 480*200/72, height = 480*200/72,
    type = "cairo", bg = "transparent")
par(mar = c(4, 4, 0.5, 0.5))	
plot(data = rats_dat, weight ~ day, col = "#4285f4")
abline(fit, col = "#EA4335")
dev.off()

rats-lm

如果矩阵 $X^{T}X$ 病态,奇异,条件数会很大,那么只有以广义逆替代 $$ \beta = (X^{T}X)^{-}X^{T}Y \\ \epsilon = Y - \hat{Y} = Y - X(X^{T}X)^{-}X^{T}Y = (I - X(X^{T}X)^{-}X^{T})Y $$

矩阵的广义逆
$A$$n\times p$ 阶矩阵,若存在一个 $p\times n$ 矩阵 $X$,使得 $AXA=A$,则称 $X$$A$ 的广义逆,记作 $X=A^{-}$
  1. 对任意矩阵 $A$,一定存在(不一定唯一)广义逆
  2. $A$ 非退化,即可逆,则 $A^-$ 唯一,且 $A^- = A^{-1}$
  3. $A(A^{\top}A)^{-}A^{\top}$ 为投影阵,也叫对称幂等阵,且与 $(A^{\top}A)^{-}$ 的取法无关

矩阵的广义逆可由函数 MASS::ginv 计算,矩阵 $X$ 的广义逆

MASS::ginv(X)
           [,1]       [,2]       [,3]       [,4]       [,5]       [,6]       [,7]
[1,]  0.0276190  0.0276190  0.0276190  0.0276190  0.0276190  0.0276190  0.0276190
[2,] -0.0009524 -0.0009524 -0.0009524 -0.0009524 -0.0009524 -0.0009524 -0.0009524
           [,8]       [,9]      [,10]      [,11]      [,12]      [,13]      [,14]
[1,]  0.0276190  0.0276190  0.0276190  0.0276190  0.0276190  0.0276190  0.0276190
[2,] -0.0009524 -0.0009524 -0.0009524 -0.0009524 -0.0009524 -0.0009524 -0.0009524
          [,15]      [,16]      [,17]      [,18]      [,19]      [,20]      [,21]
[1,]  0.0276190  0.0276190  0.0276190  0.0276190  0.0276190  0.0276190  0.0276190
[2,] -0.0009524 -0.0009524 -0.0009524 -0.0009524 -0.0009524 -0.0009524 -0.0009524
          [,22]      [,23]      [,24]      [,25]      [,26]      [,27]      [,28]
[1,]  0.0276190  0.0276190  0.0276190  0.0276190  0.0276190  0.0276190  0.0276190
[2,] -0.0009524 -0.0009524 -0.0009524 -0.0009524 -0.0009524 -0.0009524 -0.0009524
          [,29]      [,30]      [,31]      [,32]      [,33]      [,34]      [,35]
[1,]  0.0276190  0.0276190  0.0171429  0.0171429  0.0171429  0.0171429  0.0171429
[2,] -0.0009524 -0.0009524 -0.0004762 -0.0004762 -0.0004762 -0.0004762 -0.0004762
          [,36]      [,37]      [,38]      [,39]      [,40]      [,41]      [,42]
[1,]  0.0171429  0.0171429  0.0171429  0.0171429  0.0171429  0.0171429  0.0171429
[2,] -0.0004762 -0.0004762 -0.0004762 -0.0004762 -0.0004762 -0.0004762 -0.0004762
          [,43]      [,44]      [,45]      [,46]      [,47]      [,48]      [,49]
[1,]  0.0171429  0.0171429  0.0171429  0.0171429  0.0171429  0.0171429  0.0171429
[2,] -0.0004762 -0.0004762 -0.0004762 -0.0004762 -0.0004762 -0.0004762 -0.0004762
          [,50]      [,51]      [,52]      [,53]      [,54]      [,55]      [,56]
[1,]  0.0171429  0.0171429  0.0171429  0.0171429  0.0171429  0.0171429  0.0171429
[2,] -0.0004762 -0.0004762 -0.0004762 -0.0004762 -0.0004762 -0.0004762 -0.0004762
          [,57]      [,58]      [,59]      [,60]      [,61]      [,62]      [,63]
[1,]  0.0171429  0.0171429  0.0171429  0.0171429  6.667e-03  6.667e-03  6.667e-03
[2,] -0.0004762 -0.0004762 -0.0004762 -0.0004762 -5.421e-19 -5.421e-19 -5.421e-19
          [,64]      [,65]      [,66]      [,67]      [,68]      [,69]      [,70]
[1,]  6.667e-03  6.667e-03  6.667e-03  6.667e-03  6.667e-03  6.667e-03  6.667e-03
[2,] -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19
          [,71]      [,72]      [,73]      [,74]      [,75]      [,76]      [,77]
[1,]  6.667e-03  6.667e-03  6.667e-03  6.667e-03  6.667e-03  6.667e-03  6.667e-03
[2,] -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19
          [,78]      [,79]      [,80]      [,81]      [,82]      [,83]      [,84]
[1,]  6.667e-03  6.667e-03  6.667e-03  6.667e-03  6.667e-03  6.667e-03  6.667e-03
[2,] -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19
          [,85]      [,86]      [,87]      [,88]      [,89]      [,90]      [,91]
[1,]  6.667e-03  6.667e-03  6.667e-03  6.667e-03  6.667e-03  6.667e-03 -0.0038095
[2,] -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19 -5.421e-19  0.0004762
          [,92]      [,93]      [,94]      [,95]      [,96]      [,97]      [,98]
[1,] -0.0038095 -0.0038095 -0.0038095 -0.0038095 -0.0038095 -0.0038095 -0.0038095
[2,]  0.0004762  0.0004762  0.0004762  0.0004762  0.0004762  0.0004762  0.0004762
          [,99]     [,100]     [,101]     [,102]     [,103]     [,104]     [,105]
[1,] -0.0038095 -0.0038095 -0.0038095 -0.0038095 -0.0038095 -0.0038095 -0.0038095
[2,]  0.0004762  0.0004762  0.0004762  0.0004762  0.0004762  0.0004762  0.0004762
         [,106]     [,107]     [,108]     [,109]     [,110]     [,111]     [,112]
[1,] -0.0038095 -0.0038095 -0.0038095 -0.0038095 -0.0038095 -0.0038095 -0.0038095
[2,]  0.0004762  0.0004762  0.0004762  0.0004762  0.0004762  0.0004762  0.0004762
         [,113]     [,114]     [,115]     [,116]     [,117]     [,118]     [,119]
[1,] -0.0038095 -0.0038095 -0.0038095 -0.0038095 -0.0038095 -0.0038095 -0.0038095
[2,]  0.0004762  0.0004762  0.0004762  0.0004762  0.0004762  0.0004762  0.0004762
         [,120]     [,121]     [,122]     [,123]     [,124]     [,125]     [,126]
[1,] -0.0038095 -0.0142857 -0.0142857 -0.0142857 -0.0142857 -0.0142857 -0.0142857
[2,]  0.0004762  0.0009524  0.0009524  0.0009524  0.0009524  0.0009524  0.0009524
         [,127]     [,128]     [,129]     [,130]     [,131]     [,132]     [,133]
[1,] -0.0142857 -0.0142857 -0.0142857 -0.0142857 -0.0142857 -0.0142857 -0.0142857
[2,]  0.0009524  0.0009524  0.0009524  0.0009524  0.0009524  0.0009524  0.0009524
         [,134]     [,135]     [,136]     [,137]     [,138]     [,139]     [,140]
[1,] -0.0142857 -0.0142857 -0.0142857 -0.0142857 -0.0142857 -0.0142857 -0.0142857
[2,]  0.0009524  0.0009524  0.0009524  0.0009524  0.0009524  0.0009524  0.0009524
         [,141]     [,142]     [,143]     [,144]     [,145]     [,146]     [,147]
[1,] -0.0142857 -0.0142857 -0.0142857 -0.0142857 -0.0142857 -0.0142857 -0.0142857
[2,]  0.0009524  0.0009524  0.0009524  0.0009524  0.0009524  0.0009524  0.0009524
         [,148]     [,149]     [,150]
[1,] -0.0142857 -0.0142857 -0.0142857
[2,]  0.0009524  0.0009524  0.0009524

结果是 $X$

X %*% MASS::ginv(X) %*% X

简单线性回归模型残差标准差为 16.13,比之前分层正态模型 14.2,14.6,14.7 都要大

$\hat{\epsilon}$ 作为残差向量,满足 $\mathrm{E}(\hat{\epsilon}) = 0,\mathrm{Cov}(\hat{\epsilon}) = \sigma^2(I - X(X^{T}X)^{-}X^{T})$,则

  1. $\sigma^2$ 的无偏估计,也是最小二乘估计为 $\hat{\sigma}^2=\frac{\|Y - X\hat{\beta}\|^2}{n-r}, r = \mathrm{rank}(X)$
  2. $\sigma^2$ 的最大似然估计为 $\frac{n-r}{n}\hat{\sigma}^2$,且 $\frac{(n-r)\hat{\sigma}^2}{\sigma^2} \sim \chi^2_{n-r}$

计算矩阵 $X$ 的秩

qr(X)$rank
[1] 2

$\hat{\sigma}^2$

norm((diag(rep(1, 150)) - X %*% MASS::ginv(t(X) %*% X) %*% t(X)) %*% Y, type = "2") / (150 - 2)

此处等同于

norm((diag(rep(1, 150)) - X %*% solve(t(X) %*% X) %*% t(X)) %*% Y, type = "2") / (150 - 2)
[1] 1.326064

顺便一提,在 R语言中,比 t(X) %*% X 更加高效的计算方式是 crossprod(X),所以 $\beta$ 的计算可以变为

beta <- solve(crossprod(X), crossprod(X, Y))
           [,1]
[1,] 106.567619
[2,]   6.185714

更多关于线性回归的矩阵计算,见 COS海外线上沙龙第4期 邱怡轩的演讲 — 统计计算的背后 — 公式和代码之间隔了几本线代。楼主看后,写了一份 笔记,以飨读者。

3. 生长曲线模型

测量数据的时间变化添加进简单线性模型,就是生长曲线模型 Growth Curve Model,是多元线性模型的一种特例

小鼠体重 $Y$ 随时间 $t$ 的变化,以多项式的形状生长。假定不同白鼠的观测值是不相关的,而同一只白鼠的$p$次观测却是相关的。3

$$Y = \beta_0 + \beta_{1}t + \cdots + \beta_{k-1}t^{k-1},k = 1, 2, \cdots$$

其数据形式如下

$$ \begin{pmatrix} y_{11} & y_{12} & \cdots & y_{1p}\\ y_{21} & y_{22} & \cdots & y_{2p} \\ \vdots & \vdots & & \vdots \\ y_{n1} & y_{n2} & \cdots & y_{np} \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{pmatrix} (\beta_0,\beta_1,\ldots,\beta_{k-1}) \begin{pmatrix} 1 & 1 & \cdots & 1 \\ t_1 & t_2 & \cdots & t_p \\ \vdots & \vdots & & \vdots \\ t_{1}^{k-1} & t_{2}^{k-1} & \cdots & t_{p}^{k-1} \end{pmatrix} + \begin{pmatrix} \epsilon_{ij} \end{pmatrix} $$

其中 $n = 1,2,\cdots, N = 30, p = 1,2,\cdots,5$,求解上述模型,就是在已知$n\times p$的矩阵 $Y$$p$ 维向量 $t = (t_1,t_2,\cdots,t_p)$ 的情况下,确定 $k$$k$$\beta$ 向量。方便起见,我们将上述方程改写为矩阵形式,并和之前的线性模型记号保持一致

$$ Y = \mathbf{1}\beta^{\top}T^{\top} + \epsilon $$

library(lavaan)
colnames(y) <- c("t1","t2","t3","t4","t5")
model <- ' i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4 + 1*t5
           s =~ 0*t1 + 1*t2 + 2*t3 + 3*t4 + 4*t5 '
fit <- growth(model, data=y)
summary(fit)
lavaan 0.6-4 ended normally after 181 iterations

  Optimization method                           NLMINB
  Number of free parameters                         10

  Number of observations                            30

  Estimator                                         ML
  Model Fit Test Statistic                      68.986
  Degrees of freedom                                10
  P-value (Chi-square)                           0.000

Parameter Estimates:

  Information                                 Expected
  Information saturated (h1) model          Structured
  Standard Errors                             Standard

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)
  i =~                                                
    t1                1.000                           
    t2                1.000                           
    t3                1.000                           
    t4                1.000                           
    t5                1.000                           
  s =~                                                
    t1                0.000                           
    t2                1.000                           
    t3                2.000                           
    t4                3.000                           
    t5                4.000                           

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)
  i ~~                                                
    s                10.625    9.304    1.142    0.253

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)
   .t1                0.000                           
   .t2                0.000                           
   .t3                0.000                           
   .t4                0.000                           
   .t5                0.000                           
    i               155.759    2.065   75.437    0.000
    s                44.525    0.808   55.080    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .t1               32.761   12.914    2.537    0.011
   .t2               13.005    5.219    2.492    0.013
   .t3                5.181    3.498    1.481    0.139
   .t4               17.105    8.426    2.030    0.042
   .t5              157.548   45.272    3.480    0.001
    i               112.957   33.179    3.404    0.001
    s                15.909    5.124    3.105    0.002

R 内置的 ChickWeight 数据集,测量50只不同的小鸡体重随时间和喂食方式的变化,函数 SSlogis 求解

生长曲线模型的进一步延伸,就是还可能包含随机效应,那么首先应该检验随机效应的存在性4,然后估计固定效应和随机效应的参数。更多关于生长曲线建模的介绍,请看这里,混合效应模型的关于随机效应的解释 https://www.stata.com/manuals13/meme.pdf

估计

最大似然到限制极大似然估计的突破进展,就是改进方差成分的估计

Bartlett, M. S. (1937) Properties of Sufficiency and Statistical Tests. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 160:901, 268–282. https://www.jstor.org/stable/96803 其收录在 Breakthroughs in Statistics 卷1 Foundations and Basic Theory,同时修正了他论文的一个错误

Dietrich Von Rosen (1991) The growth curve model: a review. Communications in Statistics - Theory and Methods, 20:9, 2791-2822, doi:10.1080/03610929108830668

Verbyla, A. P. and Venables, W. N. (1988). An extension of the growth curve model, Biometrika, 75:1, 129–138. https://doi.org/10.1093/biomet/75.1.129

Potthoff, Richard F. and Roy, S. N. (1964) A generalized multivariate analysis of variance model useful especially for growth curve problems. Biometrika, 51:3-4, 313-326, <10.1093/biomet/51.3-4.313>


  1. Alan E. Gelfand and Susan E. Hills and Amy Racine-Poon and Adrian F. M. Smith (1990) Illustration of Bayesian Inference in Normal Data Models Using Gibbs Sampling. Journal of the American Statistical Association, 85:412, 972-985.
  2. 参考 Xiao Liu 在人大数据挖掘讨论班上的分享 https://wp-contents.netlify.com/talks/2016-xiao-liu-hierarchical-normal-model-with-mcmc.html
  3. 参考王松桂、史建红、尹素菊和吴密霞所著的《线性模型引论》第116页第四章参数估计的例4.8.1
  4. Zaixing Li (2015) Testing for Random Effects in Growth Curve Models. Communications in Statistics - Theory and Methods, 44:3, 564-572. doi:10.1080/03610926.2012.746988