利用STAN做贝叶斯回归分析：正态回归

公式中是诸如斜率之类的待估参数集，是我们拥有的数据集。是参数的先验分布，也就是我们对于可取值的已有知识。是数据的似然函数（译者注：换一个角度就是样本的联合分布）， 就是对应的后验分布。这个等式包含了参数的先验信息（即“我”认为的回归系数的正负等等），又利用观测数据的似然函数将先验信息更新以得到后验分布。而所谓的后验分布，就是在给定观测数据和先验分布的情况下，相应参数的条件分布。在频率学派的框架下，似然函数是被研究的重点。但大多数时候我们并不关心给定参数时数据的出现概率，而是专注于诸如回归方程斜率为正的概率，参数大于的概率之类的问题。利用贝叶斯回归并计算其后验概率能直接回答这类问题。

/*
*简单正态回归模型
*/
data {
    　int N; // 样本数量
      　int N2; // new_X 矩阵大小
      　int K; // 模型矩阵的列数
      　real y[N]; // 响应变量y
      　matrix[N, K] X; // 模型矩阵X
     　 matrix[N2, K] new_X; // 预测值矩阵
}
parameters {
      vector[K] beta; // 回归模型中的参数
      real sigma; // 标准差
}
transformed parameters {
      vector[N] linpred;
      linpred <- X*beta;
}
model {  
      beta[1] ~ cauchy(0, 10); // 截距项的先验 （Gelman， 2008）

      for(i in 2:K)
       beta[i] ~ cauchy(0, 2.5); // 斜率项的先验 （Gelman， 2008）

      y ~ normal(linpred, sigma);
}
generated quantities {
      vector[N2] y_pred;
      y_pred <- new_X*beta; // 模型预测的y的值
}

# 加载相应的包
library(rstan)
library(coda)
set.seed(20151204)

# 生成解释变量
dat <- data.frame(x1 = runif(100, -2, 2), x2 = runif(100, -2, 2))

# 设定模型
X <- model.matrix( ~ x1*x2, dat)

# 回归斜率项
betas <- runif(4, -1, 1)

# 仿真数据的标准差
sigma <- 1

# 生成y的仿真数据
y_norm <- rnorm(100, X%*%betas, sigma)

# 生成协变量的数据矩阵
new_X <- model.matrix( ~ x1*x2, expand.grid(x1 = seq(min(dat$x1), max(dat$x1), length = 20), x2 = c(min(dat$x2), mean(dat$x2), max(dat$x2))))

# 设定文件保存路径
setwd("/home/lionel/Desktop/Blog/STAN/")

# 回归模型
m_norm <- stan(file = "normal_regression.stan", data = list(N = 100, N2 = 60, K = 4, y = y_norm, X = X, new_X = new_X), pars = c("beta", "sigma", "y_pred"))

# 画出参数的后验分布图
post_beta <- As.mcmc.list(m_norm, pars = "beta")
plot(post_beta)

# 计算斜率项大于0的后验概率
apply(extract(m_norm, pars="beta")$beta, 2, function(x) length(which(x>0))/4000)

[1] 0 0 1 0

# 绘制参数的相关图
pairs(m_norm, pars = "beta")

# 绘制不同beta值的可信区间
plot(m_norm, pars = c("beta", "sigma"))

new_x <- data.frame(x1 = new_X[, 2], x2 = rep(c("Min", "Mean", "Max"), each = 20))
new_y <- extract(m_norm, pars = "y_pred")
pred <- apply(new_y与1相关指南：
NCCN临床实践指南：急性髓性白血病（2017.V2）
NCCN临床实践指南：乳腺癌（2017.V2）
2017 国际专家组建议：生活方式干预预防和管理代谢综合征
2017 ESPGHAN/NASPGHAN指南：儿童和青少年幽门螺杆菌感染的管理
2016 EAU指南：男性不育症
更多信息请点击：有关1更多指南
, 2, quantile, probs = c(0.025, 0.5, 0.975)) # 参数中位数和95%可信带

# 画图
plot(dat$x1, y_norm, pch = 16)
lines(new_x$x1[1:20], pred[2, 1:20], col = "red", lwd = 3)
lines(new_x$x1[1:20], pred[2, 21:40], col = "orange", lwd = 3)
lines(new_x$x1[1:20], pred[2, 41:60], col = "blue", lwd = 3)
lines(new_x$x1[1:20], pred[1, 1:20], col = "red", lwd = 1, lty = 2)
lines(new_x$x1[1:20], pred[1, 21:40], col = "orange", lwd = 1, lty = 2)
lines(new_x$x1[1:20], pred[1, 41:60], col = "blue", lwd = 1, lty = 2)
lines(new_x$x1[1:20], pred[3, 1:20], col = "red", lwd = 1, lty = 2)
lines(new_x$x1[1:20], pred[3, 21:40], col = "orange", lwd = 1, lty = 2)
lines(new_x$x1[1:20], pred[3, 41:60], col = "blue", lwd = 1, lty = 2)
legend("topright",legend=c("Min","Mean","Max"), lty = 1, col = c("red", "orange", "blue"), bty = "n", title = "Effect of x2 value on\nthe regression")

The bible in Bayesian analysis

Applied textbook focused for ecologist

The STAN reference guide

A book to rethink stats