【译文】利用STAN做贝叶斯回归分析：Part 2 非正态回归
作者 Lionel Hertzogn
前一篇文章已经介绍了怎样在R中调用STAN对正态数据进行贝叶斯回归。本文则将利用三个样例来演示怎样在R中利用STAN拟合非正态模型。
三个样例各自是negative binomial回归（过离散的泊松数据）。gamma回归（右偏的连续数据）和beta-binomial回归（过离散的二项数据）。
相关的STAN代码及一些说明会贴在本文末尾。
负二项回归
泊松分布经常使用于计数数据建模，它如果了数据的方差等于均值。当方差比均值大的时候，我们称数据存在过离散并改用负二项分布来建模。
如今如果我们手头有服从负二项分布的响应变量y以及k个解释变量X。改用方程形式表示即：
yi∼NB(μi,ϕ)
E(yi)=μi
Var(yi)=μi+μ2i/ϕ
log(μi)=β0+β1∗X1i+…+βk∗Xki
负二项分布共同拥有两个參数：μ，是必须为正数的期望，因此我们能够利用一个对数链接函数把线性模型（即解释变量乘上回归系数）映射到均值μ（见第四个方程）。
ϕ是过离散度參数。值越小意味着离散度越大，离泊松分布越远，当ϕ变得越来越大，负二项分布看起来会越来越像泊松分布。
让我们仿真一些数据并用STAN来建模
# 载入包
library(arm) # 利用里面的invlogit函数
library(emdbook) # 利用里面的rbetabinom function函数
library(rstanarm) # 利用launch_shinystan函数
# 生成服从负二项分布的仿真数据
# 解释变量
N <- 100 # 样本量
dat <- data.frame(x1 = runif(N, -2, 2), x2 = runif(N, -2, 2))
# 模型
X <- model.matrix( ~ x1*x2, dat)
K <- dim(X)[2] # 回归系数维度
# 回归的斜率
betas <- runif(K, -1, 1)
# 设定仿真数据的过离散度
ph i<- 5
# 生成响应变量
y_nb < -rnbinom(100, size = phi, mu = exp(X%*%betas))
# 拟合模型
m_nb<-stan(file = "neg_bin.stan", data = list(N=N, K=K, X=X, y=y_nb), pars = c("beta","phi","y_rep"))
# 利用shinystan诊断模型
launch_shinystan(m_nb)Shinystan界面
上述代码的最后一行命令会在你的浏览器中打开一个窗体，上面有一些选项供预计、诊断和探索模型。
一些选项超出了我有限的知识范围（如对数后验vs样本阶梯大小），所以我通常关注回归系数的后验分布（Diagnose -> NUTS (plots) -> By model parameter），出现的柱形图应该比較接近正态。我也会关注后验预測检验（Diagnose -> PPcheck -> Distribution of observed data vs replications）。y_rep的分布应当与观測数据同样。
模型看起来不错，如今我们能够利用模型的回归系数画出预測回归线以及对应的置信区间。
# 计算后验预測值和对应的可信区间
post <- as.array(m_nb)
# 我们来看看x2在x3的三个取值上的变化性
new_X <- model.matrix( ~ x1*x2, expand.grid(x2=seq(-2, 2, length=10), x1=c(min(dat$x1), mean(dat$x1), max(dat$x1))))
# 计算每一个样本的预測值
pred <- apply(post[, , 1:4], c(1, 2), FUN = function(x) new_X%*%x)
# 每条链会存在不同的矩阵，并将信息重组存于一个矩阵内
dim(pred) <- c(30, 4000)
# 提取预測中位数和95%可信区间
pred_int <- apply(pred, 1, quantile, probs=c(0.025, 0.5, 0.975))
# 绘图
plot(dat$x2, y_nb, pch=16)
lines(new_X[1:10,3], exp(pred_int[2,1:10]), col="orange", lwd=5)
lines(new_X[1:10,3], exp(pred_int[1,1:10]), col="orange", lwd=3, lty=2)
lines(new_X[1:10,3], exp(pred_int[3,1:10]), col="orange", lwd=3, lty=2)
lines(new_X[1:10,3], exp(pred_int[2,11:20]), col="red", lwd=5)
lines(new_X[1:10,3], exp(pred_int[1,11:20]), col="red", lwd=3, lty=2)
lines(new_X[1:10,3], exp(pred_int[3,11:20]), col="red", lwd=3, lty=2)
lines(new_X[1:10,3], exp(pred_int[2,21:30]), col="blue", lwd=5)
lines(new_X[1:10,3], exp(pred_int[1,21:30]), col="blue", lwd=3, lty=2)
lines(new_X[1:10,3], exp(pred_int[3,21:30]), col="blue", lwd=3, lty=2)
legend("topleft", legend=c("Min", "Mean", "Max"), ncol=3, col = c("orange", "red", "blue"), lwd = 3, bty="n", title="Value of x1")输出结果例如以下：
一如既往，我们应当关注这类模型链接空间和响应空间的差异。
上述模型在链接空间里做出了预測，如果我们想把结果绘制于响应空间，我们要利用链接函数的反函数（本例中是指数函数）来将预測值变换回去。
Gamma分布
有时我们收集到的数据呈现出右偏的特点，比方体型、植物生物量等数据。
这一类数据能够利用对数正态分布（将响应变量做对数变换）或者gamma分布建模。此处我将利用gamma分布拟合数据，如果我们有服从gamma分布的对应变量y和一些解释变量X，方程形式例如以下：
yi∼Gamma(α,β)
E(yi)=α/β
Var(yi)=α/β2
Mmmmmmm（译者注：语气词，作者卖萌）与负二项分布不同。我们不能直接利用链接函数把线性模型映射到gamma分布的某一个參数上（比方μ）。所以我们须要利用一些微积分的手段来改变模型的參数形式：
E(yi)=α/β=μ
Var(yi)=α/β2=ϕ
整理可得：
α=μ2/ϕ
β=μ/ϕ
上式中μ是分布的期望,ϕ是离散度參数,形式上和先前的负二项分布一致.由于α和β必须是正数,我们能够再次在线性模型上套一个对数链接函数:
log(μi)=β0+β1∗X1i+…+βk∗Xki
让我们生成一些数据来拟合这个模型：
# 生成服从gamma分布的数据
mus <- exp(X%*%betas)
y_g <- rgamma(100, shape = mus**2/phi, rate = mus/phi)
# 模型
m_g <- stan(file = "gamma.stan", data = list(N=N, K=K, X=X, y=y_g), pars = c("betas","phi","y_rep"))
# 模型检验
launch_shinystan(m_g)我们再一次确认了模型是正确的,页面上的全部结果看上去都很棒.由于我们利用了和前一个样例一样的链接函数，我们仅仅须要把上述代的m_nb改为m_g就能够将结果可视化。
Beta Binomial分布
最后，当我们从一个确定次数的试验中收集到了一些数据（比方一次掷10次硬币。每次正面朝上的比例），响应变量一般会服从二项分布。如今我们能够得到大量的这类数据，而且就像泊松数据可能会过离散。二项分布数据也会遇到相似情况。此时，我们就须要使用Beta-Binomial模型：
yi∼BetaBinomial(N,α,β)
E(yi)=N∗α/(α+β)
把方差表达式先放在一边(点击这里看看这公式多丑)，常数N(试验次数)也能够无论。和gamma分布的样例一样，我们能够将方程整理为别的形式，新方程的期望为μ。离散參数为ϕ:
α=μ∗ϕ
β=(1−μ)∗ϕ
这次μ代表概率因而必须落在0。1之间。我们能够利用logit链接函数将线性买模型映射到μ
logit(μ)=β0+β1∗X1i+…+βk∗Xki
再来见证仿真的威力吧：
# 生成Beta-Binomial数据
W <- rep(20, 100) # 试验次数
y_bb <- rbetabinom(100, prob = invlogit(X%*%betas), size = W, theta = phi)
# 模型
m_bb <- stan(file = "beta_bin.stan", data = list(N=N, W=W, K=K, X=X, y=y_bb), pars=c("betas", "phi", "y_rep"))
# 模型检验
launch_shinystan(m_bb)一切看起来如此美妙。让我们把结果画出来：
# 计算后验预測值和置信区间
post<-as.array(m_bb) # 后验抽样
# 后验预測值
pred<-apply(post[,,1:4],c(1,2),FUN = function(x) new_X%*%x)
# 每条链都存在重组信息的不同矩阵里
dim(pred)<-c(30,4000)
# 得到后验越策值的中位数和95%置信区间
pred_int<-apply(pred,1,quantile,probs=c(0.025,0.5,0.975))
# 绘图
plot(dat$x2, y_bb, pch = 16)
lines(new_X[1:10,3], 20*invlogit(pred_int[2,1:10]), col="orange", lwd=5)
lines(new_X[1:10,3], 20*invlogit(pred_int[1,1:10]), col="orange", lwd=3, lty=2)
lines(new_X[1:10,3], 20*invlogit(pred_int[3,1:10]), col="orange", lwd=3, lty=2)
lines(new_X[1:10,3], 20*invlogit(pred_int[2,11:20]), col="red", lwd=5)
lines(new_X[1:10,3], 20*invlogit(pred_int[1,11:20]), col="red", lwd=3, lty=2)
lines(new_X[1:10,3], 20*invlogit(pred_int[3,11:20]), col="red", lwd=3, lty=2)
lines(new_X[1:10,3], 20*invlogit(pred_int[2,21:30]), col="blue", lwd=5)
lines(new_X[1:10,3], 20*invlogit(pred_int[1,21:30]), col="blue", lwd=3, lty=2)
lines(new_X[1:10,3], 20*invlogit(pred_int[3,21:30]), col="blue", lwd=3, lty=2)
legend("topleft", legend = c("Min", "Mean", "Max"), ncol=3, col = c("orange", "red", "blue"), lwd = 3, bty="n", title="Value of x1")输出例如以下：
注意exp函数要改为invlogit函数，由于这里用的链接函数不同。
零星想法
通过本文我们学习了怎样拟合非正态模型。
STAN是一个很灵活的工具，能够拟合许不同分布和多类參数形式（參看參考手冊），其可能性仅仅受到你的如果的限制（也许还有你的数学水平。。。）。在这里我声明一下，rsranarm包能够让你在不写出模型形式的情况下拟合STAN模型。仅仅须要写一些诸如glm函数的典型的R代码（请參考此文）。那么，我们还须要学习STAN吗？这取决于你的目的，如果你的模型比較经常使用。那么利用rstanarm包能够节约大量时间，让你把精力集中于科研。而且rstanarm包内已有的函数执行速度很快。但如果有一天你认为你会在模型上有所突破，须要拟合一个个性化的模型，那么STAN会是一种很灵活和高效的建模语言。
模型代码
/*
*简单负二项回归样例
*使用负二项回归的另外一种參数形式
*细节參看Stan參考手冊 section 40.1-3
*/
data {
int N; //the number of observations
int K; //the number of columns in the model matrix
int y[N]; //the response
matrix[N,K] X; //the model matrix
}
parameters {
vector[K] beta; //the regression parameters
real phi; //the overdispersion parameters
}
transformed parameters {
vector[N] mu;//the linear predictor
mu <- exp(X*beta); //using the log link
}
model {
beta[1] ~ cauchy(0,10); //prior for the intercept following Gelman 2008
for(i in 2:K)
beta[i] ~ cauchy(0,2.5);//prior for the slopes following Gelman 2008
y ~ neg_binomial_2(mu,phi);
}
generated quantities {
vector[N] y_rep;
for(n in 1:N){
y_rep[n] <- neg_binomial_2_rng(mu[n],phi); //posterior draws to get posterior predictive checks
}
}Gamma
/*
*简单gamma分布样例
*注意我使用了对数链接函数，大多数时候它比典型反向链接更有意义
*/
data {
int N; //the number of observations
int K; //the number of columns in the model matrix
real y[N]; //the response
matrix[N,K] X; //the model matrix
}
parameters {
vector[K] betas; //the regression parameters
real phi; //the variance parameter
}
transformed parameters {
vector[N] mu; //the expected values (linear predictor)
vector[N] alpha; //shape parameter for the gamma distribution
vector[N] beta; //rate parameter for the gamma distribution
mu <- exp(X*betas); //using the log link
alpha <- mu .* mu / phi;
beta <- mu / phi;
}
model {
betas[1] ~ cauchy(0,10); //prior for the intercept following Gelman 2008
for(i in 2:K)
betas[i] ~ cauchy(0,2.5);//prior for the slopes following Gelman 2008
y ~ gamma(alpha,beta);
}
generated quantities {
vector[N] y_rep;
for(n in 1:N){
y_rep[n] <- gamma_rng(alpha[n],beta[n]); //posterior draws to get posterior predictive checks
}
}Beta-Binomial
/*
*简单beta-binomial样例
*/
data {
int N; //the number of observations
int K; //the number of columns in the model matrix
int y[N]; //the response
matrix[N,K] X; //the model matrix
int W[N]; //the number of trials per observations, ie a vector of 1 for a 0/1 dataset
}
parameters {
vector[K] betas; //the regression parameters
real phi; //the overdispersion parameter
}
transformed parameters {
vector[N] mu; //the linear predictor
vector[N] alpha; //the first shape parameter for the beta distribution
vector[N] beta; //the second shape parameter for the beta distribution
for(n in 1:N)
mu[n] <- inv_logit(X[n,]*betas); //using logit link
alpha <- mu * phi;
beta <- (1-mu) * phi;
}
model {
betas[1] ~ cauchy(0,10); //prior for the intercept following Gelman 2008
for(i in 2:K)
betas[i] ~ cauchy(0,2.5);//prior for the slopes following Gelman 2008
y ~ beta_binomial(W,alpha,beta);
}
generated quantities {
vector[N] y_rep;
for(n in 1:N){
y_rep[n] <- beta_binomial_rng(W[n],alpha[n],beta[n]); //posterior draws to get posterior predictive checks
}
}