在 Stan 中开发非线性增长曲线模型的分层版本
Developing hierarchical version of nonlinear growth curve model in Stan
以下模型是 Preece 和 Baines (1978, Annals of Human Biology) 的模型 1，用于描述人类的成长。
这个模型我的Stan代码如下：
```{stan output.var="test"}
data {
int<lower=1> n;
ordered[n] t; // age
ordered[n] y; // height of human
}
parameters {
positive_ordered[2] h;
real<lower=0, upper=t[n-1]>theta;
positive_ordered[2] s;
real<lower=0>sigma;
}
model {
h[1] ~ uniform(0, y[n]);
h[2] ~ normal(180, 20);
sigma ~ student_t(2, 0, 1);
s[1] ~ normal(5, 5);
s[2] ~ normal(5, 5);
theta ~ normal(10, 5);
y ~ normal(h[2] - (2*(h[2] - h[1]) * inv(exp(s[1]*(t - theta)) + exp(s[2]*(t - theta)))), sigma);
}
```
我现在想创建此模型的分层版本，其参数建模在男孩之间相同（例如测量变异性 sigma）或分层版本（例如 h_1，成人身高）。
至于参数 sigma 和 theta，我想保持这些先验在男孩中与 sigma ~ student_t(2, 0, 1); 和 theta ~ normal(10, 5); 相同。
不幸的是，我几乎没有实施层次建模的经验，而且我一直在努力尝试做贝叶斯教科书中简单的二项式层次模型之外的任何层次示例（参见[=的第 12 章，第 358-359 页） 44=]统计反思 作者：Richard McElreath）。但是，我确实了解分层建模背后的理论，正如安德鲁·盖尔曼 (Andrew Gelman) 的 贝叶斯数据分析 第 5 章和 统计反思 第 12 章所写作者：理查德·麦克埃尔瑞斯。
我很想知道如何在 Stan 中实现这种类型的层次模型。理想情况下，我正在寻找代码旁边的解释，以便我将来可以学习如何独立实现这些类型的层次模型示例。
我的数据样本如下：
age boy01 boy02 boy03 boy04 boy05 boy06 boy07 boy08 boy09 boy10 boy11 boy12 boy13 boy14 boy15 boy16 boy17 boy18
1 1 81.3 76.2 76.8 74.1 74.2 76.8 72.4 73.8 75.4 78.8 76.9 81.6 78 76.4 76.4 76.2 75 79.7
2 1.25 84.2 80.4 79.8 78.4 76.3 79.1 76 78.7 81 83.3 79.9 83.7 81.8 79.4 81.2 79.2 78.4 81.3
3 1.5 86.4 83.2 82.6 82.6 78.3 81.1 79.4 83 84.9 87 84.1 86.3 85 83.4 86 82.3 82 83.3
4 1.75 88.9 85.4 84.7 85.4 80.3 84.4 82 85.8 87.9 89.6 88.5 88.8 86.4 87.6 89.2 85.4 84 86.5
5 2 91.4 87.6 86.7 88.1 82.2 87.4 84.2 88.4 90 91.4 90.6 92.2 87.1 91.4 92.2 88.4 85.9 88.9
6 3 101. 97 94.2 98.6 89.4 94 93.2 97.3 97.3 100. 96.6 99.3 96.2 101. 101. 101 95.6 99.4
7 4 110. 105. 100. 104. 96.9 102. 102. 107. 103. 111 105. 106. 104 106. 110. 107. 102. 104.
8 5 116. 112. 107. 111 104. 109. 109 113. 108. 118. 112 113. 111 113. 117. 115. 109. 112.
9 6 122. 119. 112. 116. 111. 116. 117. 119. 114. 126. 119. 120. 117. 120. 122. 121. 118. 119
10 7 130 125 119. 123. 116. 122. 123. 126. 120. 131. 125. 127. 124. 129. 130. 128 125. 128
我承认小数点不够精确。数据采用 tibble table 的形式，它似乎不响应 R 的常用命令以获得更高的精度。为了保持一致性，最好忽略第 5 行之后的所有行，因为第 1 - 5 行显示原始数据中存在的完整精度。
完整数据中，年龄为
> Children$age
[1] 1.00 1.25 1.50 1.75 2.00 3.00 4.00 5.00 6.00 7.00 8.00 8.50 9.00 9.50 10.00 10.50 11.00 11.50 12.00 12.50
[21] 13.00 13.50 14.00 14.50 15.00 15.50 16.00 16.50 17.00 17.50 18.00
并且有 39 个男孩，以与上述示例相同的宽数据格式列出。
免责声明：首先让我们使用 Stan 拟合一个（非分层的）非线性增长模型。
我们读入了示例数据。
library(tidyverse);
df <- read.table(text = "
age boy01 boy02 boy03 boy04 boy05 boy06 boy07 boy08 boy09 boy10 boy11 boy12 boy13 boy14 boy15 boy16 boy17 boy18
1 1 81.3 76.2 76.8 74.1 74.2 76.8 72.4 73.8 75.4 78.8 76.9 81.6 78 76.4 76.4 76.2 75 79.7
2 1.25 84.2 80.4 79.8 78.4 76.3 79.1 76 78.7 81 83.3 79.9 83.7 81.8 79.4 81.2 79.2 78.4 81.3
3 1.5 86.4 83.2 82.6 82.6 78.3 81.1 79.4 83 84.9 87 84.1 86.3 85 83.4 86 82.3 82 83.3
4 1.75 88.9 85.4 84.7 85.4 80.3 84.4 82 85.8 87.9 89.6 88.5 88.8 86.4 87.6 89.2 85.4 84 86.5
5 2 91.4 87.6 86.7 88.1 82.2 87.4 84.2 88.4 90 91.4 90.6 92.2 87.1 91.4 92.2 88.4 85.9 88.9
6 3 101. 97 94.2 98.6 89.4 94 93.2 97.3 97.3 100. 96.6 99.3 96.2 101. 101. 101 95.6 99.4
7 4 110. 105. 100. 104. 96.9 102. 102. 107. 103. 111 105. 106. 104 106. 110. 107. 102. 104.
8 5 116. 112. 107. 111 104. 109. 109 113. 108. 118. 112 113. 111 113. 117. 115. 109. 112.
9 6 122. 119. 112. 116. 111. 116. 117. 119. 114. 126. 119. 120. 117. 120. 122. 121. 118. 119
10 7 130 125 119. 123. 116. 122. 123. 126. 120. 131. 125. 127. 124. 129. 130. 128 125. 128", header = T, row.names = 1);
df <- df %>%
gather(boy, height, -age);
我们定义了 Stan 模型。
model <- "
data {
int N; // Number of observations
real y[N]; // Height
real t[N]; // Time
}
parameters {
real<lower=0,upper=1> s[2];
real h_theta;
real theta;
real sigma;
}
transformed parameters {
vector[N] mu;
real h1;
h1 = max(y);
for (i in 1:N)
mu[i] = h1 - 2 * (h1 - h_theta) / (exp(s[1] * (t[i] - theta)) + (exp(s[2] * (t[i] - theta))));
}
model {
// Priors
s ~ cauchy(0, 2.5);
y ~ normal(mu, sigma);
}
"
这里我们考虑 s[1] 和 s[2].
根据数据拟合模型。
library(rstan);
options(mc.cores = parallel::detectCores());
rstan_options(auto_write = TRUE);
model <- stan(
model_code = model,
data = list(
N = nrow(df),
y = df$height,
t = df$age),
iter = 4000);
6个模型参数汇总：
summary(model, pars = c("h1", "h_theta", "theta", "s", "sigma"))$summary
# mean se_mean sd 2.5% 25% 50%
#h1 131.0000000 0.000000000 0.0000000 131.0000000 131.0000000 131.0000000
#h_theta 121.6874553 0.118527828 2.7554944 115.4316738 121.1654809 122.2134014
#theta 6.5895553 0.019738319 0.5143429 5.4232740 6.4053479 6.6469534
#s[1] 0.7170836 0.214402086 0.3124318 0.1748077 0.3843143 0.8765256
#s[2] 0.3691174 0.212062373 0.3035039 0.1519308 0.1930381 0.2066811
#sigma 3.1524819 0.003510676 0.1739904 2.8400096 3.0331962 3.1411533
# 75% 97.5% n_eff Rhat
#h1 131.0000000 131.0000000 8000.000000 NaN
#h_theta 123.0556379 124.3928800 540.453594 1.002660
#theta 6.8790801 7.3376348 679.024115 1.002296
#s[1] 0.9516115 0.9955989 2.123501 3.866466
#s[2] 0.2954409 0.9852540 2.048336 6.072550
#sigma 3.2635849 3.5204101 2456.231113 1.001078
那么这是什么意思呢？从 s[1] 和 s[2] 的 Rhat 值可以看出这两个参数的收敛性存在问题。这是因为 s[1] 和 s[2] 是不可区分的：不能同时估计它们。 s[1] 和 s[2] 上更强的正则化先验可能会将 s 参数之一驱动为零。
我不确定我是否理解 s[1] 和 s[2] 的要点。从统计建模的角度来看，您无法在我们正在考虑的简单非线性增长模型中获得两个参数的估计值。
更新
正如这里所承诺的那样，这是一个更新。这变得很长 post，我试图通过添加额外的解释使事情尽可能清楚。
初步评论
使用 positive_ordered 作为 s 的数据类型在解决方案的收敛性方面产生了显着差异。我不清楚为什么会这样，我也不知道 Stan 是如何实现的 positive_ordered，但它确实有效。
在这种分层方法中，我们通过考虑 h1 ~ normal(mu_h1, sigma_h1)，在超参数 mu_h1 ~ normal(max(y), 10) 上使用先验，在 sigma_h1 ~ cauchy(0, 10)（它是半柯西因为 sigma 被声明为 real<lower=0>）。
老实说，我不确定某些参数的解释（和可解释性）。 h_1 和 h_theta 的估计值非常相似，并且在某种程度上相互抵消。我想这会在拟合模型时产生一些收敛问题，但正如您进一步看到的那样，Rhat 值似乎没问题。尽管如此，由于我对模型、数据及其上下文了解不够，我对这些参数的可解释性仍然持怀疑态度。从统计建模的角度来看，通过将一些其他参数转换为组级参数来扩展模型很简单；但是我想困难会来自于不可区分性和缺乏可解释性。
抛开所有这些要点，这应该为您提供一个如何实现分层模型的实际示例。
斯坦模型
model_code <- "
data {
int N; // Number of observations
int J; // Number of boys
int<lower=1,upper=J> boy[N]; // Boy of observation
real y[N]; // Height
real t[N]; // Time
}
parameters {
real<lower=0> h1[J];
real<lower=0> h_theta;
real<lower=0> theta;
positive_ordered[2] s;
real<lower=0> sigma;
// Hyperparameters
real<lower=0> mu_h1;
real<lower=0> sigma_h1;
}
transformed parameters {
vector[N] mu;
for (i in 1:N)
mu[i] = h1[boy[i]] - 2 * (h1[boy[i]] - h_theta) / (exp(s[1] * (t[i] - theta)) + (exp(s[2] * (t[i] - theta))));
}
model {
h1 ~ normal(mu_h1, sigma_h1); // Partially pool h1 parameters across boys
mu_h1 ~ normal(max(y), 10); // Prior on h1 hyperparameter mu
sigma_h1 ~ cauchy(0, 10); // Half-Cauchy prior on h1 hyperparameter sigma
h_theta ~ normal(max(y), 2); // Prior on h_theta
theta ~ normal(max(t), 2); // Prior on theta
s ~ cauchy(0, 1); // Half-Cauchy priors on s[1] and s[2]
y ~ normal(mu, sigma);
}
"
总而言之：我们通过将成人身高参数建模为 h1 ~ normal(mu_h1, sigma_h1)，其中超参数 mu_h1 和 sigma_h1 表征所有男孩成人身高值的正态分布。我们在超参数上选择弱信息先验，并在所有剩余参数上选择弱信息先验，类似于第一个完全池化示例。
拟合模型
# Fit model
fit <- stan(
model_code = model_code,
data = list(
N = nrow(df),
J = length(unique(df$boy)),
boy = df$boy,
y = df$height,
t = df$age),
iter = 4000)
提取摘要
我们提取所有参数的参数估计；请注意，我们现在有与组（即男孩）一样多的 h1 参数。
# Get summary
summary(fit, pars = c("h1", "h_theta", "theta", "s", "sigma"))$summary
# mean se_mean sd 2.5% 25% 50%
#h1[1] 142.9406153 0.1046670943 2.41580757 138.4272280 141.2858391 142.909765
#h1[2] 143.7054020 0.1070466445 2.46570025 139.1301456 142.0233342 143.652657
#h1[3] 144.0352331 0.1086953809 2.50145442 139.3982034 142.3131167 143.971473
#h1[4] 143.8589955 0.1075753575 2.48015745 139.2689731 142.1666685 143.830347
#h1[5] 144.7359976 0.1109871908 2.55284812 140.0529359 142.9917503 144.660586
#h1[6] 143.9844938 0.1082691127 2.49497990 139.3378948 142.2919990 143.926931
#h1[7] 144.3857221 0.1092604239 2.51645359 139.7349112 142.6665955 144.314645
#h1[8] 143.7469630 0.1070594855 2.46860328 139.1748700 142.0660983 143.697302
#h1[9] 143.6841113 0.1072208284 2.47391295 139.0885987 141.9839040 143.644357
#h1[10] 142.9518072 0.1041206784 2.40729732 138.4289207 141.3114204 142.918407
#h1[11] 143.5352502 0.1064173663 2.45712021 138.9607665 141.8547610 143.483157
#h1[12] 143.0941582 0.1050061258 2.42894673 138.5579378 141.4295430 143.055576
#h1[13] 143.6194965 0.1068494690 2.46574352 138.9426195 141.9412820 143.577920
#h1[14] 143.4477182 0.1060254849 2.44776536 138.9142081 141.7708660 143.392231
#h1[15] 143.1415683 0.1049131998 2.42575487 138.6246642 141.5014391 143.102219
#h1[16] 143.5686919 0.1063594201 2.45328456 139.0064573 141.8962853 143.510276
#h1[17] 144.0170715 0.1080567189 2.49269747 139.4162885 142.3138300 143.965127
#h1[18] 143.4740997 0.1064867748 2.45545200 138.8768051 141.7989566 143.426211
#h_theta 134.3394366 0.0718785944 1.72084291 130.9919889 133.2348411 134.367152
#theta 8.2214374 0.0132434918 0.45236221 7.4609612 7.9127800 8.164685
#s[1] 0.1772044 0.0004923951 0.01165119 0.1547003 0.1705841 0.177522
#s[2] 1.6933846 0.0322953612 1.18334358 0.6516669 1.1630900 1.463148
#sigma 2.2601677 0.0034146522 0.13271459 2.0138514 2.1657260 2.256678
# 75% 97.5% n_eff Rhat
#h1[1] 144.4795105 147.8847202 532.7265 1.008214
#h1[2] 145.2395543 148.8419618 530.5599 1.008187
#h1[3] 145.6064981 149.2080965 529.6183 1.008087
#h1[4] 145.4202919 149.0105666 531.5363 1.008046
#h1[5] 146.3200407 150.0701757 529.0592 1.008189
#h1[6] 145.5551778 149.1365279 531.0372 1.008109
#h1[7] 145.9594956 149.5996605 530.4593 1.008271
#h1[8] 145.3032680 148.8695637 531.6824 1.008226
#h1[9] 145.2401743 148.7674840 532.3662 1.008023
#h1[10] 144.4811712 147.9218834 534.5465 1.007937
#h1[11] 145.1153635 148.5968945 533.1235 1.007988
#h1[12] 144.6479561 148.0546831 535.0652 1.008115
#h1[13] 145.1660639 148.6562172 532.5386 1.008138
#h1[14] 144.9975197 148.5273804 532.9900 1.008067
#h1[15] 144.6733010 148.1130207 534.6057 1.008128
#h1[16] 145.1163764 148.6027096 532.0396 1.008036
#h1[17] 145.5578107 149.2014363 532.1519 1.008052
#h1[18] 145.0249329 148.4886949 531.7060 1.008055
#h_theta 135.4870338 137.6753254 573.1698 1.006818
#theta 8.4812339 9.2516700 1166.7226 1.002306
#s[1] 0.1841457 0.1988365 559.9036 1.005333
#s[2] 1.8673249 4.1143099 1342.5839 1.001562
#sigma 2.3470429 2.5374239 1510.5824 1.001219
可视化成人身高估计值
最后，我们绘制了所有男孩的成人身高估计值 h_1，包括 50% 和 97% 的置信区间。
# Plot h1 values
summary(fit, pars = c("h1"))$summary %>%
as.data.frame() %>%
rownames_to_column("Variable") %>%
mutate(
Variable = gsub("(h1\[|\])", "", Variable),
Variable = df$key[match(Variable, df$boy)]) %>%
ggplot(aes(x = `50%`, y = Variable)) +
geom_point(size = 3) +
geom_segment(aes(x = `2.5%`, xend = `97.5%`, yend = Variable), size = 1) +
geom_segment(aes(x = `25%`, xend = `75%`, yend = Variable), size = 2) +
labs(x = "Median (plus/minus 95% and 50% CIs)", y = "h1")
以下模型是 Preece 和 Baines (1978, Annals of Human Biology) 的模型 1，用于描述人类的成长。
这个模型我的Stan代码如下：
```{stan output.var="test"}
data {
int<lower=1> n;
ordered[n] t; // age
ordered[n] y; // height of human
}
parameters {
positive_ordered[2] h;
real<lower=0, upper=t[n-1]>theta;
positive_ordered[2] s;
real<lower=0>sigma;
}
model {
h[1] ~ uniform(0, y[n]);
h[2] ~ normal(180, 20);
sigma ~ student_t(2, 0, 1);
s[1] ~ normal(5, 5);
s[2] ~ normal(5, 5);
theta ~ normal(10, 5);
y ~ normal(h[2] - (2*(h[2] - h[1]) * inv(exp(s[1]*(t - theta)) + exp(s[2]*(t - theta)))), sigma);
}
```
我现在想创建此模型的分层版本，其参数建模在男孩之间相同（例如测量变异性 sigma）或分层版本（例如 h_1，成人身高）。
至于参数 sigma 和 theta，我想保持这些先验在男孩中与 sigma ~ student_t(2, 0, 1); 和 theta ~ normal(10, 5); 相同。
不幸的是，我几乎没有实施层次建模的经验，而且我一直在努力尝试做贝叶斯教科书中简单的二项式层次模型之外的任何层次示例（参见[=的第 12 章，第 358-359 页） 44=]统计反思 作者：Richard McElreath）。但是，我确实了解分层建模背后的理论，正如安德鲁·盖尔曼 (Andrew Gelman) 的 贝叶斯数据分析 第 5 章和 统计反思 第 12 章所写作者：理查德·麦克埃尔瑞斯。
我很想知道如何在 Stan 中实现这种类型的层次模型。理想情况下，我正在寻找代码旁边的解释，以便我将来可以学习如何独立实现这些类型的层次模型示例。
我的数据样本如下：
age boy01 boy02 boy03 boy04 boy05 boy06 boy07 boy08 boy09 boy10 boy11 boy12 boy13 boy14 boy15 boy16 boy17 boy18
1 1 81.3 76.2 76.8 74.1 74.2 76.8 72.4 73.8 75.4 78.8 76.9 81.6 78 76.4 76.4 76.2 75 79.7
2 1.25 84.2 80.4 79.8 78.4 76.3 79.1 76 78.7 81 83.3 79.9 83.7 81.8 79.4 81.2 79.2 78.4 81.3
3 1.5 86.4 83.2 82.6 82.6 78.3 81.1 79.4 83 84.9 87 84.1 86.3 85 83.4 86 82.3 82 83.3
4 1.75 88.9 85.4 84.7 85.4 80.3 84.4 82 85.8 87.9 89.6 88.5 88.8 86.4 87.6 89.2 85.4 84 86.5
5 2 91.4 87.6 86.7 88.1 82.2 87.4 84.2 88.4 90 91.4 90.6 92.2 87.1 91.4 92.2 88.4 85.9 88.9
6 3 101. 97 94.2 98.6 89.4 94 93.2 97.3 97.3 100. 96.6 99.3 96.2 101. 101. 101 95.6 99.4
7 4 110. 105. 100. 104. 96.9 102. 102. 107. 103. 111 105. 106. 104 106. 110. 107. 102. 104.
8 5 116. 112. 107. 111 104. 109. 109 113. 108. 118. 112 113. 111 113. 117. 115. 109. 112.
9 6 122. 119. 112. 116. 111. 116. 117. 119. 114. 126. 119. 120. 117. 120. 122. 121. 118. 119
10 7 130 125 119. 123. 116. 122. 123. 126. 120. 131. 125. 127. 124. 129. 130. 128 125. 128
我承认小数点不够精确。数据采用 tibble table 的形式，它似乎不响应 R 的常用命令以获得更高的精度。为了保持一致性，最好忽略第 5 行之后的所有行，因为第 1 - 5 行显示原始数据中存在的完整精度。
完整数据中，年龄为
> Children$age
[1] 1.00 1.25 1.50 1.75 2.00 3.00 4.00 5.00 6.00 7.00 8.00 8.50 9.00 9.50 10.00 10.50 11.00 11.50 12.00 12.50
[21] 13.00 13.50 14.00 14.50 15.00 15.50 16.00 16.50 17.00 17.50 18.00
并且有 39 个男孩，以与上述示例相同的宽数据格式列出。
免责声明：首先让我们使用 Stan 拟合一个（非分层的）非线性增长模型。
我们读入了示例数据。
library(tidyverse); df <- read.table(text = " age boy01 boy02 boy03 boy04 boy05 boy06 boy07 boy08 boy09 boy10 boy11 boy12 boy13 boy14 boy15 boy16 boy17 boy18 1 1 81.3 76.2 76.8 74.1 74.2 76.8 72.4 73.8 75.4 78.8 76.9 81.6 78 76.4 76.4 76.2 75 79.7 2 1.25 84.2 80.4 79.8 78.4 76.3 79.1 76 78.7 81 83.3 79.9 83.7 81.8 79.4 81.2 79.2 78.4 81.3 3 1.5 86.4 83.2 82.6 82.6 78.3 81.1 79.4 83 84.9 87 84.1 86.3 85 83.4 86 82.3 82 83.3 4 1.75 88.9 85.4 84.7 85.4 80.3 84.4 82 85.8 87.9 89.6 88.5 88.8 86.4 87.6 89.2 85.4 84 86.5 5 2 91.4 87.6 86.7 88.1 82.2 87.4 84.2 88.4 90 91.4 90.6 92.2 87.1 91.4 92.2 88.4 85.9 88.9 6 3 101. 97 94.2 98.6 89.4 94 93.2 97.3 97.3 100. 96.6 99.3 96.2 101. 101. 101 95.6 99.4 7 4 110. 105. 100. 104. 96.9 102. 102. 107. 103. 111 105. 106. 104 106. 110. 107. 102. 104. 8 5 116. 112. 107. 111 104. 109. 109 113. 108. 118. 112 113. 111 113. 117. 115. 109. 112. 9 6 122. 119. 112. 116. 111. 116. 117. 119. 114. 126. 119. 120. 117. 120. 122. 121. 118. 119 10 7 130 125 119. 123. 116. 122. 123. 126. 120. 131. 125. 127. 124. 129. 130. 128 125. 128", header = T, row.names = 1); df <- df %>% gather(boy, height, -age);我们定义了 Stan 模型。
model <- " data { int N; // Number of observations real y[N]; // Height real t[N]; // Time } parameters { real<lower=0,upper=1> s[2]; real h_theta; real theta; real sigma; } transformed parameters { vector[N] mu; real h1; h1 = max(y); for (i in 1:N) mu[i] = h1 - 2 * (h1 - h_theta) / (exp(s[1] * (t[i] - theta)) + (exp(s[2] * (t[i] - theta)))); } model { // Priors s ~ cauchy(0, 2.5); y ~ normal(mu, sigma); } "这里我们考虑
s[1]和s[2]. 的弱信息（正则化）先验
根据数据拟合模型。
library(rstan); options(mc.cores = parallel::detectCores()); rstan_options(auto_write = TRUE); model <- stan( model_code = model, data = list( N = nrow(df), y = df$height, t = df$age), iter = 4000);6个模型参数汇总：
summary(model, pars = c("h1", "h_theta", "theta", "s", "sigma"))$summary # mean se_mean sd 2.5% 25% 50% #h1 131.0000000 0.000000000 0.0000000 131.0000000 131.0000000 131.0000000 #h_theta 121.6874553 0.118527828 2.7554944 115.4316738 121.1654809 122.2134014 #theta 6.5895553 0.019738319 0.5143429 5.4232740 6.4053479 6.6469534 #s[1] 0.7170836 0.214402086 0.3124318 0.1748077 0.3843143 0.8765256 #s[2] 0.3691174 0.212062373 0.3035039 0.1519308 0.1930381 0.2066811 #sigma 3.1524819 0.003510676 0.1739904 2.8400096 3.0331962 3.1411533 # 75% 97.5% n_eff Rhat #h1 131.0000000 131.0000000 8000.000000 NaN #h_theta 123.0556379 124.3928800 540.453594 1.002660 #theta 6.8790801 7.3376348 679.024115 1.002296 #s[1] 0.9516115 0.9955989 2.123501 3.866466 #s[2] 0.2954409 0.9852540 2.048336 6.072550 #sigma 3.2635849 3.5204101 2456.231113 1.001078那么这是什么意思呢？从
s[1]和s[2]的Rhat值可以看出这两个参数的收敛性存在问题。这是因为s[1]和s[2]是不可区分的：不能同时估计它们。s[1]和s[2]上更强的正则化先验可能会将s参数之一驱动为零。
我不确定我是否理解 s[1] 和 s[2] 的要点。从统计建模的角度来看，您无法在我们正在考虑的简单非线性增长模型中获得两个参数的估计值。
更新
正如这里所承诺的那样，这是一个更新。这变得很长 post，我试图通过添加额外的解释使事情尽可能清楚。
初步评论
使用
positive_ordered作为s的数据类型在解决方案的收敛性方面产生了显着差异。我不清楚为什么会这样，我也不知道 Stan 是如何实现的positive_ordered，但它确实有效。在这种分层方法中，我们通过考虑
h1 ~ normal(mu_h1, sigma_h1)，在超参数mu_h1 ~ normal(max(y), 10)上使用先验，在sigma_h1 ~ cauchy(0, 10)（它是半柯西因为sigma被声明为real<lower=0>）。老实说，我不确定某些参数的解释（和可解释性）。
h_1和h_theta的估计值非常相似，并且在某种程度上相互抵消。我想这会在拟合模型时产生一些收敛问题，但正如您进一步看到的那样，Rhat值似乎没问题。尽管如此，由于我对模型、数据及其上下文了解不够，我对这些参数的可解释性仍然持怀疑态度。从统计建模的角度来看，通过将一些其他参数转换为组级参数来扩展模型很简单；但是我想困难会来自于不可区分性和缺乏可解释性。
抛开所有这些要点，这应该为您提供一个如何实现分层模型的实际示例。
斯坦模型
model_code <- "
data {
int N; // Number of observations
int J; // Number of boys
int<lower=1,upper=J> boy[N]; // Boy of observation
real y[N]; // Height
real t[N]; // Time
}
parameters {
real<lower=0> h1[J];
real<lower=0> h_theta;
real<lower=0> theta;
positive_ordered[2] s;
real<lower=0> sigma;
// Hyperparameters
real<lower=0> mu_h1;
real<lower=0> sigma_h1;
}
transformed parameters {
vector[N] mu;
for (i in 1:N)
mu[i] = h1[boy[i]] - 2 * (h1[boy[i]] - h_theta) / (exp(s[1] * (t[i] - theta)) + (exp(s[2] * (t[i] - theta))));
}
model {
h1 ~ normal(mu_h1, sigma_h1); // Partially pool h1 parameters across boys
mu_h1 ~ normal(max(y), 10); // Prior on h1 hyperparameter mu
sigma_h1 ~ cauchy(0, 10); // Half-Cauchy prior on h1 hyperparameter sigma
h_theta ~ normal(max(y), 2); // Prior on h_theta
theta ~ normal(max(t), 2); // Prior on theta
s ~ cauchy(0, 1); // Half-Cauchy priors on s[1] and s[2]
y ~ normal(mu, sigma);
}
"
总而言之：我们通过将成人身高参数建模为 h1 ~ normal(mu_h1, sigma_h1)，其中超参数 mu_h1 和 sigma_h1 表征所有男孩成人身高值的正态分布。我们在超参数上选择弱信息先验，并在所有剩余参数上选择弱信息先验，类似于第一个完全池化示例。
拟合模型
# Fit model
fit <- stan(
model_code = model_code,
data = list(
N = nrow(df),
J = length(unique(df$boy)),
boy = df$boy,
y = df$height,
t = df$age),
iter = 4000)
提取摘要
我们提取所有参数的参数估计；请注意，我们现在有与组（即男孩）一样多的 h1 参数。
# Get summary
summary(fit, pars = c("h1", "h_theta", "theta", "s", "sigma"))$summary
# mean se_mean sd 2.5% 25% 50%
#h1[1] 142.9406153 0.1046670943 2.41580757 138.4272280 141.2858391 142.909765
#h1[2] 143.7054020 0.1070466445 2.46570025 139.1301456 142.0233342 143.652657
#h1[3] 144.0352331 0.1086953809 2.50145442 139.3982034 142.3131167 143.971473
#h1[4] 143.8589955 0.1075753575 2.48015745 139.2689731 142.1666685 143.830347
#h1[5] 144.7359976 0.1109871908 2.55284812 140.0529359 142.9917503 144.660586
#h1[6] 143.9844938 0.1082691127 2.49497990 139.3378948 142.2919990 143.926931
#h1[7] 144.3857221 0.1092604239 2.51645359 139.7349112 142.6665955 144.314645
#h1[8] 143.7469630 0.1070594855 2.46860328 139.1748700 142.0660983 143.697302
#h1[9] 143.6841113 0.1072208284 2.47391295 139.0885987 141.9839040 143.644357
#h1[10] 142.9518072 0.1041206784 2.40729732 138.4289207 141.3114204 142.918407
#h1[11] 143.5352502 0.1064173663 2.45712021 138.9607665 141.8547610 143.483157
#h1[12] 143.0941582 0.1050061258 2.42894673 138.5579378 141.4295430 143.055576
#h1[13] 143.6194965 0.1068494690 2.46574352 138.9426195 141.9412820 143.577920
#h1[14] 143.4477182 0.1060254849 2.44776536 138.9142081 141.7708660 143.392231
#h1[15] 143.1415683 0.1049131998 2.42575487 138.6246642 141.5014391 143.102219
#h1[16] 143.5686919 0.1063594201 2.45328456 139.0064573 141.8962853 143.510276
#h1[17] 144.0170715 0.1080567189 2.49269747 139.4162885 142.3138300 143.965127
#h1[18] 143.4740997 0.1064867748 2.45545200 138.8768051 141.7989566 143.426211
#h_theta 134.3394366 0.0718785944 1.72084291 130.9919889 133.2348411 134.367152
#theta 8.2214374 0.0132434918 0.45236221 7.4609612 7.9127800 8.164685
#s[1] 0.1772044 0.0004923951 0.01165119 0.1547003 0.1705841 0.177522
#s[2] 1.6933846 0.0322953612 1.18334358 0.6516669 1.1630900 1.463148
#sigma 2.2601677 0.0034146522 0.13271459 2.0138514 2.1657260 2.256678
# 75% 97.5% n_eff Rhat
#h1[1] 144.4795105 147.8847202 532.7265 1.008214
#h1[2] 145.2395543 148.8419618 530.5599 1.008187
#h1[3] 145.6064981 149.2080965 529.6183 1.008087
#h1[4] 145.4202919 149.0105666 531.5363 1.008046
#h1[5] 146.3200407 150.0701757 529.0592 1.008189
#h1[6] 145.5551778 149.1365279 531.0372 1.008109
#h1[7] 145.9594956 149.5996605 530.4593 1.008271
#h1[8] 145.3032680 148.8695637 531.6824 1.008226
#h1[9] 145.2401743 148.7674840 532.3662 1.008023
#h1[10] 144.4811712 147.9218834 534.5465 1.007937
#h1[11] 145.1153635 148.5968945 533.1235 1.007988
#h1[12] 144.6479561 148.0546831 535.0652 1.008115
#h1[13] 145.1660639 148.6562172 532.5386 1.008138
#h1[14] 144.9975197 148.5273804 532.9900 1.008067
#h1[15] 144.6733010 148.1130207 534.6057 1.008128
#h1[16] 145.1163764 148.6027096 532.0396 1.008036
#h1[17] 145.5578107 149.2014363 532.1519 1.008052
#h1[18] 145.0249329 148.4886949 531.7060 1.008055
#h_theta 135.4870338 137.6753254 573.1698 1.006818
#theta 8.4812339 9.2516700 1166.7226 1.002306
#s[1] 0.1841457 0.1988365 559.9036 1.005333
#s[2] 1.8673249 4.1143099 1342.5839 1.001562
#sigma 2.3470429 2.5374239 1510.5824 1.001219
可视化成人身高估计值
最后，我们绘制了所有男孩的成人身高估计值 h_1，包括 50% 和 97% 的置信区间。
# Plot h1 values
summary(fit, pars = c("h1"))$summary %>%
as.data.frame() %>%
rownames_to_column("Variable") %>%
mutate(
Variable = gsub("(h1\[|\])", "", Variable),
Variable = df$key[match(Variable, df$boy)]) %>%
ggplot(aes(x = `50%`, y = Variable)) +
geom_point(size = 3) +
geom_segment(aes(x = `2.5%`, xend = `97.5%`, yend = Variable), size = 1) +
geom_segment(aes(x = `25%`, xend = `75%`, yend = Variable), size = 2) +
labs(x = "Median (plus/minus 95% and 50% CIs)", y = "h1")