【译文】利用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 }}

注:原文刊载于datascience+网站

链接:http://datascienceplus.com/bayesian-regression-with-stan-beyond-normality/