WinBUGS在统计分析中的应用(第二部分)

第一节 WinBUGS数据分析案例

在这一节中,我将拿一个经典的研究数据,利用WinBUGS给出简单的分析。首先介绍一下这个数据:Seeds

seed O. aegyptiaco 75 	seed O. aegyptiaco 73	Bean 	Cucumber 	Bean 	Cucumber
r 	n 	r/n 	r 	n 	r/n 	r 	n 	r/n 	r 	n 	r/n
10 	39 	0.26 	5 	6 	0.83 	8 	16 	0.5 	3 	12 	0.25
23 	62 	0.37 	53 	74 	0.72 	10 	30 	0.33 	22 	41 	0.54
23 	81 	0.28 	55 	72 	0.76 	8 	28 	0.29 	15 	30 	0.5
26 	51 	0.51 	32 	51 	0.63 	23 	45 	0.51 	32 	51 	0.63
17 	39 	0.44 	46 	79 	0.58 	0 	4 	0 	3 	7 	0.43
10 	13 	0.77

这个数据来自Crowder (1978)。之后Breslow and Clayton (1993) 作为例子,也分析过这个数据。数据反映的是某一品种的豆类种子和某一品种的黄瓜种子,分别放在21个培养皿(plates)中分别培养,在根提取物aegyptiaco 75和aegyptiaco 73的作用下的出芽率的差异。其中r是出芽的个数,n是种子的个数,而r/n是出芽率。我们用random effect logistic regression模型来进行分析(注意,在Bayesian分析中,通常是将covariates看做是服从某一个分布的随机变量,这和一般意义上的GLM, GLME, LME中对于covariates解释是不同的):

$r_i \sim Binomial(p_i, ~n_i)$

$logit(p_i)=\alpha_0 + \alpha_1 x_{1i} + \alpha_2 x_{2i} + \alpha_{12} x_{1i} x_{2i} + b_i$

$b_i \sim Normal(0, \tau)$

其中$x_{1i}$是种子的类型,$x_{2i}$是根提取物的类型,$\alpha_{12} x_{1i} x_{2i}$是交互项, $\alpha_0,~\alpha_1,~\alpha_2,~\alpha_{12},~\tau$是给定的独立的 “noninformative” 先验参数。在Bayesian分析中,通常我们会定义一个DAG图(即Directed Acyclic Graph有向无圈图) 。我们可以在WinBUGS中通过设计DAG图来定义模型。不过这一节中我们还是用WinBUGS中的BUGS语言来定义模型,如何在WinBUGS中通过设计DAG图来定义模型我将在下一节中详细介绍,但是必须要说明的是BUGS语言比DAG图灵活,不过直观性不如后者。

模型

model
{
 for( i in 1 : N ) {
  r[i] ~ dbin(p[i],n[i])
  b[i] ~ dnorm(0.0,tau)
  logit(p[i]) <- alpha0 + alpha1 * x1[i] + alpha2 * x2[i] +
   alpha12 * x1[i] * x2[i] + b[i]
 }
 alpha0 ~ dnorm(0.0,1.0E-6)
 alpha1 ~ dnorm(0.0,1.0E-6)
 alpha2 ~ dnorm(0.0,1.0E-6)
 alpha12 ~ dnorm(0.0,1.0E-6)
 tau ~ dgamma(0.001,0.001)
 sigma <- 1 / sqrt(tau)
}
WinBUGS doodle模型
WinBUGS doodle模型

数据

list(r = c(10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10,   8, 10,   8, 23, 0,  3, 22, 15, 32, 3),
   n = c(39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41, 30, 51, 7),
   x1 = c(0,   0,  0,   0,   0, 0,   0,   0,  0,   0,   0,  1,   1,   1,   1, 1,   1,  1,   1,   1, 1),
   x2 = c(0,   0,  0,   0,   0, 1,   1,   1,  1,   1,   1,  0,   0,   0,   0, 0,   1,  1,   1,   1, 1),
   N = 21)

初始值

list(alpha0 = 0, alpha1 = 0, alpha2 = 0, alpha12 = 0, tau = 1)

经过10000次迭代,我们得到参数的估计如下:

   node    mean     sd MC error   2.50%  median  97.50% start
 alpha0 -0.5546 0.1941 0.007696 -0.9353 -0.5577 -0.1597  1001
 alpha1 0.08497 0.3127  0.01283 -0.5814 0.09742  0.6679  1001
alpha12 -0.8229 0.4321  0.01785  -1.697 -0.8218 0.01641  1001
 alpha2   1.356 0.2743  0.01236  0.8257   1.347   1.909  1001
  sigma  0.2731 0.1437 0.007956 0.04133  0.2654  0.5862  1001

第二节 结合SAS做比较分析

下面我们用同样的数据在SAS中进行分析,看看结果相差多少:

首先生成数据:

data seeds;
 input plate r n seed extract;
 datalines;
1   10  39  0   0
2   23  62 0   0
3   23  81 0   0
4   26  51 0   0
5   17  39 0   0
6    5   6  0   1
7   53  74 0   1
8   55  72 0   1
9   32  51 0   1
10  46  79 0   1
11  10  13 0   1
12   8  16 1   0
13  10  30 1   0
14   8  28 1   0
15  23  45 1   0
16   0   4  1   0
17   3  12 1   1
18  22  41 1   1
19  15  30 1   1
20  32  51 1   1
21   3   7  1   1
;
run;
data seeds;
 set seeds;
 interact=seed*extract;
run;
proc print; run;

然后调用GENMOD过程,拟合经典的logistic regression回归模型

proc genmod data=seeds;
/* Basic logistic regression WHITHOUT random plate effect */
title ‘Classical logistic regression’;
 model r/n= seed extract interact / dist=B;
run;

得到

Analysis Of Parameter Estimates
Parameter	DF	Estimate	Standard Error	Wald 95% Confidence Limits	Chi-Square	Pr > ChiSq
Intercept	1	-0.5582		0.126		-0.8052	-0.3112			19.62	     	<.0001
seed		1	0.1459		0.2232		-0.2915	0.5833			0.43		0.5132
extract		1	1.3182		0.1775		0.9704	1.666			55.17	      	<.0001
interact	1	-0.7781		0.3064		-1.3787	-0.1775			6.45		0.0111
Scale		0	1		0		1	1

调用NLMIXED过程,拟合经典的logistic regression回归模型

proc nlmixed data=seeds;
/* Cassical logistic regression using NLMIXED */
title 'Classical logistic regression with NLMIXED';
parms b0=-0.55 b1=0.15 b2=1.32 b12=-0.78;
logitp=b0+b1*seed+b2*extract+b12*interact;
p=exp(logitp)/(1+exp(logitp));
model r ~ binomial(n,p) ;
run;

得到:

Parameter Estimates
Parameter	Estimate	Standard Error	DF	t Value	Pr > |t|	Alpha	Lower	Upper	Gradient
b0		-0.5582		0.126		21	-4.43	0.0002		0.05	-0.8202	-0.2961	-0.00000229
b1		0.1459		0.2232		21	0.65	0.5203		0.05	-0.3182	0.61	-8.82E-07
b2		1.3182		0.1775		21	7.43	<.0001		0.05	0.9491	1.6872	-0.00000161
b12		-0.7781		0.3064		21	-2.54	0.0191		0.05	-1.4154	-0.1408	-6.61E-07

调用NLMIXED过程,拟合经典带随机截距的logistic regression回归模型

proc nlmixed data=seeds;
/* Logistic regression + RANDOM plate INTERCEPT */
title 'Logistic regression with random plate intercept (NLMIXED)';
parms b0=-0.55 b1=0.15 b2=1.32 b12=-0.78;
logitp=b0+b1*seed+b2*extract+b12*interact+alpha;
p=exp(logitp)/(1+exp(logitp));
random alpha ~ normal(0,varu) subject=plate out=seedmixed;
model r ~ binomial(n,p) ;
run;

得到:

Parameter Estimates
Parameter	Estimate	Standard Error	DF	t Value	Pr > |t|	Alpha	Lower	Upper	Gradient
b0		-0.5484		0.1666		20	-3.29	0.0036		0.05	-0.896	-0.2009	-0.00087
b1		0.097		0.278		20	0.35	0.7308		0.05	-0.483	0.677	0.00022
b2		1.337		0.2369		20	5.64	<.0001		0.05	0.8428	1.8313	-0.00018
b12		-0.8104		0.3852		20	-2.1	0.0482		0.05	-1.6139	-0.007	0.00037
1.07		0.295 varu	0.05581		0.05	20	9		0.05	-0.0527	0.1643	0.001515

当然了winBUGS的强大之处并不在于此,而是在处理诸如GLME(有些文献称GLMM),空间数据模型等计算复杂的模型,之后还会继续讨论。

参考文献:

[1] Crowder M J (1978) Beta-binomial Anova for proportions. Applied Statistics. 27, 34-37.

[2] Breslow N E and Clayton D G (1993) Approximate inference in generalized linear mixed models. Journal of the American Statistical Association. 88, 9-25.

最后再送出一本书,供大家研究参考

Disease Mapping with WINBUGS and MLWin(djvu格式)

WinBUGS在统计分析中的应用 第二部分完

WinBUGS在统计分析中的应用(第二部分)》有18个想法

  1. 谢谢楼主的讲座和书!
    to priss111:谢老大和楼主都很牛–这么说就对啦!

  2. 很感谢楼主的文章,让我受益很多!
    我是winbugs的新手,我在做时遇到了一个问题,向您请教一下。模型是这样的: y=e+J*P, 期中 e ~ dnorm(u0,tau0), J~dnorm(u1,tau1), P~dpois(lanta), y的数据是已知的,我想要求出参数u0,u1,tau0,tau1,and lanta. 假设先验分布: u0,u1~dnorm, tau0,tau1~dgamma, lanta~ dbeta, and y~ dnorm,这些分布是有文献支持的. 但是我发现这样做出来的结果不好:首先我先假设参数的值,产生y的数据集,然后用y来求出这些参数,对比一下,发现差距比较大。我想是否我代码编错了?

    model
    {
    for(i in 1:N) {
    yy[i] ~ dnorm (y[i],tau)
    cq[i] ~ dnorm (u1,tau1)
    bq[i] ~ dnorm(u0,tau0)
    dq[i] ~ dpois(lanta)
    y[i] <- bq[i]+cq[i]*dq[i]
    }
    u1 ~dnorm (0.01, 0.0002)
    u0 ~dnorm (1, 0.0002)
    tau0~dgamma (4, 2)
    tau1 ~dgamma (4, 2)
    lanta~dbeta (2,1)
    tau<- 1/(1/tau0+1/tau1*pow(1-exp(-lanta),2))
    }
    1. 问题不是出在你的model上,而是出在你的实验设计上。你先假设参数的值,产生y的数据集,然后用y来求出这些参数,理论上来说,只有y无限多,且产生这些y的参数遍历其所有空间时,才能返回来求出这些参数的真值。举一个简单的例子:在normal distribution P(y|theta, tau)中,我已知theta=0,tau=1,用这样的参数来产生y;假如我产生了5次data, 非常非常巧,我5次data的值都是0,因此你的数据集y={0, 0, 0, 0, 0}.那么你能用这样的数据集来估计出真正的theta和tau么?答案是否定的,你只能估计出正确的theta,而tau却非常的不一样。所以说,你产生的data要能正确地反映出你model,你才能正确的估计它的参数。这是一个parameter identifiability 的问题。希望这样的解答对你有用。

      1. icwei您好,也很想请教您一个问题:

        您说的”问题不是出在你的model上,而是出在你的实验设计上”
        这句话中的‘出在实验设计上’怎么理解,能不能再通俗的解释一下。

        从您这句话中也受到启发。 因为我也遇到了一个这样的问题,
        当然,不排除还有其他问题的可能性。

      2. 我的意思是原作者有一个model,他为了检验这个model的正确性,便设计了一个实验,这个实验就是先假设所有参数已知,用这些参数产生data,然后再用这些data来估计这些参数,对比求得的参数和原参数的值来证明model是否正确。这样的实验没有大问题,但实施的时候要小心,你需要的检验你所产生的这些data是否能反映model的所有特性。

        实际上,原作者提出的验证model的方法并不正确,即使他能够得到和原参数一样的值也并不能说明这个model对于实际的data就是正确的。他犯了一个基本的错误,既他用frequentist的方法来验证bayesian的正确性。在frenquentist的概念中,model是存在真值的,所以你可以提出一个hypothesis去检验估计值是否等于真值;而在bayesian中,我们是用credit interval 来表示所求的参数,既真值并不存在,而是与我们的prior knowledge相关的一个范围。

        实际问题中,我们更多的是只得到一些数据,然后用这些数据去估计我们的参数。正确的判断参数是否能够反映系统正确性的方法是看posterior predictive distribution, 既P(y_pred|y). 如果你的posterior predictive distribution能够正确的模拟你的data,那说明这个model所估计的参数是真正能够反映你系统的特性的;反之不能。

      3. 不好意思,我学统计学的时候用英文学的,所以可能用中文表达起来有些难懂。举个不完全恰当但简单的例子吧,像我上面说的,我有dataset={0, 0, 0, 0, 0},他们是从一个有噪音的系统产生的,噪音符合normal distribution,但均值和方差均不为0。那么我用一个normal distribution的model去fit这些data,得到的结果为mean=0, variance=0.这个model可以完美的fit我当前的data.可是当我用这个model求posterior predictive distribution的时候,我所得到的是P(y_pred|y)=N(0, 0), 其所相对应得predictive dataset为y_pred={0, 0, 0, 0, 0, …..}. 可是当我们让这个系统继续产生更多的data时,我们会发现我们model产生的值并不能很好的反映系统的特性(尽管估计的参数能够很好的fit我们原始的dataset)。

        那么,我们的疑问是:我们的model错了么???但在这个特殊的例子里,我们可以说,model没有错,只是他估计的参数错了。那么参数不正确是由什么引起的呢?那就是我们的实验设计。这个实验设计为只得到5个数据,而这5个数据不足以反映系统的全部特性,因此才估计出错误的参数。

        因此,我们说,参数是否能代表系统的特性,不但model要能反映真实系统的特性,数据也要能够反映真实系统的特性,这涉及到实验设计的问题。通常情况下,只要考虑posterior predictive distribution是否能反映系统产生的未来数据的特性就可以了。

  3. 谢谢icwen这么认真的回复,从中又受到一些启发。 哦,我基本上明白: 对于bayesian中的参数估计,需要用后验分布来检验所估计的参数是否能真实的反应实际情形;而其他类型参数估计的模拟实验,则直接用所估计到的参数与模拟时设定的参数进行比较,越吻合说明该估计方法对所估计的参数即有效也稳健。

    不知道这样理解对不对。 另外,还想请教您一个问题:不只您是不是学数学的背景?

    1. 你理解的基本正确,但注意不要混淆posterior distribution(后验分布)和posterior predictive distribution(后验估计分布?不知道这个中文名是否正确)的概念。他们的区别是:假设我们的参数为a,数据为y,那么后验分布是指P(a|y),而后验估计分布是指P(y_pred|y),它等于P(y_pred|y)=integral(P(y_pred|a)*P(a|y))da,既你对未来数据的预测是基于当前model的likelihood,以后验分布作为piror,进行加权而得到的。他是对你估计参数的所有值得一个综合性考量。

      我不是纯数学背景出身,我大本学的是自控和信号处理,从博士开始专攻统计信号处理,博士毕业后开始转入统计学在生物医学方面的数学建模,现在主要兴趣在Bayesian statistics, mcmc, hierarchical graphical model, differential equation 以及 stochastic system。有兴趣的话可以和我联系,我们可以共同探讨。我的email是 chen.wei@mrc-bsu.cam.ac.uk

      1. icwei,您好!有个问题请教您:我有个问题想请教您:对于无信息先验分布(non – informative prior)是怎么定义的?我看了一些材料,多处都是按照以下方式定义的。但是对于无信息先验分布,方差不是应该足够大吗(In case of non-informative priors, the variances should have been specified as very large or instead uniform distributions covering the entire range of plausible values should have been used. )?能否帮我解释下?Thanks.
        model
        {
        for (i in 1:N)
        {
        pd[i]<-pow(d[i],c)
        pred[i]<-1.3+a*(1-exp(-b*pd[i]))

        h[i]~dnorm(pred[i],prec)
        }

        # Priors
        a ~ dnorm(0, 1.0E-6)
        b ~ dnorm(0, 1.0E-6)
        c ~ dnorm(0, 1.0E-6)
        prec~dgamma(0.001,0.001)

        }

  4. 感谢楼主的文章和各位同学的讨论,在下受益匪浅。
    在我最近的一篇文章中,我打算小试一下贝叶斯模型,在经过很多次尝试后,都已失败告终,感觉很深奥,所以想和各位探讨一二。
    我的问题是找出各林分因子对树木死亡率的影响,即stand basal area(SBA),the ratio of focal species basal area to stand basal area (FSBA), stand age (SA)对 tree mortality rate 影响。假设我总共有PN个sample plots,每个plots,我有N个树是起初的活的株数,有S棵树是在一段时间(time)后的活树数量(不包括后来长大的),一般情况,每个plot的树木死亡率通过(ln(N)-ln(S))/time来计算获得,可是我的有些plot(一般以上的plots),树木没有死亡,或者是全部死亡,即N=S,or S=0,所以我想参照Condit R.(2006)和 Kraft J. N.(2010)的hieracichal bayies算死亡率(mp)的方法,把以上两种情况的plots都用在实际分析中。
    他们的idea是这样的,对[i]th species (plots, in my cases),树木的存活率(sp[i]),N[i],S[i]满足binomial distribution,即 S[i]~dbin(sp[i],N[i]),and, sp[i]=exp(-mp[i]*time),mp为死亡率,这样就实现了mp=(ln(N)-ln(S))/time这个运算过程,然后他们又假设among species,死亡率是满足lognormal distribution,并且得到数据上的验证,即有strong right skews,所以 mp[i]~dlnorm(logmu,logsd),那么,我现在想和大家探讨的第一个问题是,我这样写出他们的winbugs model是不是正确的?
    model
    { for(i in 1:speciesnumber)
    { S[i]~dbin(sp[i],N[i])
    sp[i]<-exp(-mp[i]*time[i])
    mp[i]~dlnorm(logmu,tau)
    }
    logmu~dnorm(1,0.00000001)
    tau~dgamma(1,1)
    }
    接下来,就到我的case中了,我检查的一下我的among plots mortality rate distribution,我发现也stong right skews,所以我也想用了上面的方法得mp,然后我通过mp=a+b*SBA+c*FSBA+d*SA,得出各系数a,b,c,d。我的model是这样写的:
    model
    { for (i in 1:PN)
    { S[i]~dbin(sp[i],N[i])
    sp[i]<-exp(-mp[i]*time)
    mp[i]~dlnorm(logmu,tau)
    mp[i]<-a+b*SBA+c*FSBA+d*SA
    }
    logmu~dnorm(1,0.0000001)
    tau~dgamma(1,1)
    a~dnorm(1,0.0000001)
    b~dnorm(1.0.0000001)
    c~dnorm(1.0.0000001)
    d~dnorm(1.0.0000001)
    }
    可是我这样写,却被提示“multiple definisitons of node mp[1]”,我想请问一下,我的model到底出了什么错误?
    一万个感谢!

    Condit, R., Ashton, P., Bunyavejchewin, S., Dattaraja, H.S., Davies, S., Esufali, S., Ewango, C., Foster, R., Gunatilleke, I.A.U.N., Gunatilleke, C.V.S., Hall, P., Harms, K.E., Hart, T., Hernandez, C., Hubbell, S., Itoh, A., Kiratiprayoon, S., LaFrankie, J., de Lao, S.L., Makana, J.R., Noor, M.N.S., Kassim, A.R., Russo, S., Sukumar, R., Samper, C., Suresh, H.S., Tan, S., Thomas, S., Valencia, R., Vallejo, M., Villa, G., and Zillio, T. 2006. The importance of demographic niches to tree diversity. Science 313(5783): 98-101.
    Kraft, N.J.B., Metz, M.R., Condit, R.S., and Chave, J. 2010. The relationship between wood density and mortality in a global tropical forest data set. New Phytologist 188(4): 1124-1136. doi: 10.1111/j.1469-8137.2010.03444.x

  5. 我想要用MCMC方法估计参数,看到您这篇文章以后做了个程序,可是原理不是很清楚,不知道对不对,cigam总是不出结果。
    分布函数是F(x)=(1-F(x) )G(x-k)+F(k)=1-〖[1+ε((x-k)/σ)]〗^(-1⁄ε)
    ε服从Pareto(a,c)

    σ服从gamma(a,b)
    不知道这样写程序是否正确呢?

    model
    {
    for(i in 1:N)
    {
    p[i] <-tao*pow(1+cigam*tao*(x[i]-k),-1/cigam-1)
    }
    tao~ dgamma(0.01,0.01)
    cigam~ dpar(0.1,0.1)
    }
    data
    list(N=53,k=308,x=c(309,310,311,313,313,316,319,320,328,330,332,332,338,341,341,349,351,354,356,365,368,370,371,375,384,384,385,389,389,391,396,408,417,418,418,435,447,448,460,463,465,473,500,507,521,522,554,618,702,754,842,858,1152))
    initial values
    list(tao=1,cigam=1)

    谢谢,麻烦你啦嘻嘻

  6. 带入您的代码,发现如下提示:
    this chain contains uninitialized variables
    请更正…谢谢!

发表评论

邮箱地址不会被公开。 必填项已用*标注