WinBUGS 威布尔网络元分析
WinBUGS Weibull Network Meta-Analysis
我目前正在对多个临床试验的生存数据进行荟萃分析。
为此,我使用相同的方法从已发布的分析中获取代码。但是,当 运行 使用已发布分析中的数据使用此代码时,我无法复制他们的结果。事实上,结果无法收敛到任何一种合理的估计。
代码本身(不包括数据)应该是正确的,因为它直接来自作者。我认为问题必须与初始值或
采样的参数是 运行,但在玩了很多之后
初始值、老化时间、变薄等...我还没有得到有意义的结果。
如果有人能就如何 运行 这个(初始值等...)正确地 运行 提出建议,我将不胜感激。或者,如果代码中存在问题或者数据的设置方式与代码不匹配,了解这些信息将很有用。
附带说明一下,我正在使用 R2WinBUGs 进行分析,尽管我有
单独使用 WinBUGs 遇到了同样的问题。
A bit of extra background on the method:
The way this works is by estimating the difference in shape and scale
parameters of a reparameterized Weibull distribution between
treatments across multiple studies using random effects.
The Weibull distribution is reparameterized such that log of the
hazard rate is a+b*log(t) where a is a scale parameter and b is a
shape parameter. From this, you can calculate the the likelihood
function of a given number of failures out of a given number of
patients over an interval.
Unfortunately, the article is public, but if you can get access here
is the link:
http://onlinelibrary.wiley.com/doi/10.1002/jrsm.25/abstract;jsessionid=2BA8F0D9BEF9A33F84975618D33F8DD9.f03t03?userIsAuthenticated=false&deniedAccessCustomisedMessage=
输入模型的变量的快速摘要:
NT: Number of separate treatments included.
N: Number of rows in the main dataset.
NS: Number of studies
s: Study that the line of data corresponds to (this is numbered 1:6)
r: number of patients failing in interval for this treatment/study
n: number of patients at risk at start of interval for this
treatment/study
t: Treatment that this line of data corresponds to (numbered 1:3)
b: Indicates which treatment is the baseline one to which others are compared (set to 1 for each line).
bs: Treatment that is the control arm of this study
bt: Treatment that is the research arm of this study
WinBUGS代码(含数据):
#Winbugs code for random effects networks meta-analysis model
Model
{
for (i in 1:N)
{ # N=number of data points in dataset
#likelihood
r[i]~ dbin(p[i],n[i])
p[i]<-1-exp(-h[i]*dt[i]) # hazard h over interval [t,t+dt] # expressed as deaths per unit person-time (e.g. months)
#random effects model
log(h[i])<-nu[i]+log(time[i])*theta[i]
nu[i]<-mu[s[i],1]+delta[s[i],1]*(1-equals(t[i],b[i]))
theta[i]<-mu[s[i],2]+ delta[s[i],2]*(1-equals(t[i],b[i]))
}
for(k in 1 :NS)
{ # NS=number of studies in dataset
delta[k,1:2]~dmnorm(md[k,1:2],omega[1:2,1:2])
md[k,1]<-d[ts[k],1]-d[bs[k],1]
md[k,2]<-d[ts[k],2]-d[bs[k],2]
}
# priors
d[1,1]<-0
d[1,2]<-0
for(j in 2 :NT)
{ # NT=number of treatments
d[j,1:2] ~ dmnorm(mean[1:2],prec2[,])
}
for(k in 1 :NS)
{
mu[k,1:2] ~ dmnorm(mean[1:2],prec2[,])
}
omega[1:2, 1:2] ~ dwish(R[1:2,1:2],2)
}
# Winbugs data set
list(N=242, NS=6, NT=3,
mean=c(0,0),
prec2 = structure(.Data = c(
0.0001,0,
0,0.0001), .Dim = c(2,2)),
R = structure(.Data = c(
0.01,0,
0,0.01), .Dim = c(2,2))
)
s[] r[] n[] t[] b[] time[] dt[]
1 15 152 3 1 3 3
1 11 140 3 1 6 3
1 8 129 3 1 9 3
1 9 121 3 1 12 3
1 9 112 3 1 15 3
1 3 83 3 1 18 3
1 4 80 3 1 21 3
1 5 76 3 1 24 3
1 2 71 3 1 27 3
1 2 41 3 1 30 3
1 1 39 3 1 33 3
1 3 38 3 1 36 3
1 2 35 3 1 39 3
1 1 33 3 1 42 3
1 3 32 3 1 45 3
1 3 29 3 1 48 3
1 2 26 3 1 51 3
1 1 24 3 1 54 3
1 1 23 3 1 57 3
1 1 22 3 1 60 3
1 10 149 1 1 3 3
1 11 140 1 1 6 3
1 9 128 1 1 9 3
1 5 119 1 1 12 3
1 6 114 1 1 15 3
1 3 72 1 1 18 3
1 5 70 1 1 21 3
1 4 65 1 1 24 3
1 7 61 1 1 27 3
1 2 34 1 1 30 3
1 2 32 1 1 33 3
1 3 30 1 1 36 3
1 2 27 1 1 39 3
1 2 25 1 1 42 3
1 1 23 1 1 45 3
1 2 22 1 1 48 3
1 1 19 1 1 51 3
1 2 19 1 1 54 3
1 1 17 1 1 57 3
1 0 16 1 1 60 3
2 4 125 2 1 3 3
2 4 121 2 1 6 3
2 2 117 2 1 9 3
2 5 114 2 1 12 3
2 2 109 2 1 15 3
2 3 107 2 1 18 3
2 2 104 2 1 21 3
2 4 94 2 1 24 3
2 4 90 2 1 27 3
2 3 81 2 1 30 3
2 4 78 2 1 33 3
2 3 61 2 1 36 3
2 5 58 2 1 39 3
2 1 48 2 1 42 3
2 2 47 2 1 45 3
2 3 41 2 1 48 3
2 0 38 2 1 51 3
2 3 29 2 1 54 3
2 3 26 2 1 57 3
2 2 18 2 1 60 3
2 0 16 2 1 63 3
2 1 10 2 1 66 3
2 0 9 2 1 69 3
2 0 3 2 1 72 3
2 0 3 2 1 75 3
2 0 3 2 1 78 3
2 15 196 1 1 3 3
2 9 179 1 1 6 3
2 10 170 1 1 9 3
2 9 162 1 1 12 3
2 9 153 1 1 15 3
2 5 141 1 1 18 3
2 5 136 1 1 21 3
2 10 121 1 1 24 3
2 5 111 1 1 27 3
2 7 92 1 1 30 3
2 7 85 1 1 33 3
2 4 71 1 1 36 3
2 6 67 1 1 39 3
2 4 53 1 1 42 3
2 5 49 1 1 45 3
2 6 36 1 1 48 3
2 3 30 1 1 51 3
2 2 26 1 1 54 3
2 2 24 1 1 57 3
2 0 13 1 1 60 3
2 1 13 1 1 63 3
2 1 11 1 1 66 3
2 1 10 1 1 69 3
2 0 6 1 1 72 3
2 0 6 1 1 75 3
2 0 6 1 1 78 3
3 6 113 2 1 3 3
3 4 105 2 1 6 3
3 3 101 2 1 9 3
3 1 97 2 1 12 3
3 9 96 2 1 15 3
3 4 84 2 1 18 3
3 2 80 2 1 21 3
3 4 74 2 1 24 3
3 3 70 2 1 27 3
3 2 59 2 1 30 3
3 0 57 2 1 33 3
3 6 51 2 1 36 3
3 2 45 2 1 39 3
3 1 37 2 1 42 3
3 3 36 2 1 45 3
3 1 27 2 1 48 3
3 1 26 2 1 51 3
3 2 25 2 1 54 3
3 7 116 1 1 3 3
3 6 111 1 1 6 3
3 4 105 1 1 9 3
3 3 99 1 1 12 3
3 9 96 1 1 15 3
3 5 85 1 1 18 3
3 5 80 1 1 21 3
3 3 68 1 1 24 3
3 7 65 1 1 27 3
3 8 48 1 1 30 3
3 4 40 1 1 33 3
3 2 33 1 1 36 3
3 0 31 1 1 39 3
3 1 28 1 1 42 3
3 2 27 1 1 45 3
3 3 20 1 1 48 3
3 1 17 1 1 51 3
3 0 16 1 1 54 3
4 10 167 2 1 3 3
4 5 149 2 1 6 3
4 6 145 2 1 9 3
4 3 138 2 1 12 3
4 4 135 2 1 15 3
4 5 128 2 1 18 3
4 2 122 2 1 21 3
4 2 120 2 1 24 3
4 7 104 2 1 27 3
4 9 98 2 1 30 3
4 5 89 2 1 33 3
4 2 57 2 1 36 3
4 2 55 2 1 39 3
4 4 53 2 1 42 3
4 2 49 2 1 45 3
4 2 26 2 1 48 3
4 1 24 2 1 51 3
4 1 23 2 1 54 3
4 1 11 2 1 57 3
4 0 10 2 1 60 3
4 0 10 2 1 63 3
4 2 164 1 1 3 3
4 5 153 1 1 6 3
4 7 148 1 1 9 3
4 6 141 1 1 12 3
4 12 135 1 1 15 3
4 6 119 1 1 18 3
4 4 113 1 1 21 3
4 3 109 1 1 24 3
4 5 98 1 1 27 3
4 2 94 1 1 30 3
4 2 92 1 1 33 3
4 4 55 1 1 36 3
4 3 50 1 1 39 3
4 1 48 1 1 42 3
4 2 47 1 1 45 3
4 1 22 1 1 48 3
4 1 21 1 1 51 3
4 0 20 1 1 54 3
4 1 7 1 1 57 3
4 0 6 1 1 60 3
4 0 6 1 1 63 3
5 12 152 2 1 3 3
5 7 135 2 1 6 3
5 9 128 2 1 9 3
5 8 120 2 1 12 3
5 7 112 2 1 15 3
5 1 77 2 1 18 3
5 3 76 2 1 21 3
5 2 73 2 1 24 3
5 4 71 2 1 27 3
5 5 45 2 1 30 3
5 3 40 2 1 33 3
5 2 37 2 1 36 3
5 3 35 2 1 39 3
5 3 32 2 1 42 3
5 3 32 2 1 45 3
5 1 32 2 1 48 3
5 9 149 1 1 3 3
5 4 131 1 1 6 3
5 5 127 1 1 9 3
5 8 122 1 1 12 3
5 11 114 1 1 15 3
5 5 76 1 1 18 3
5 5 71 1 1 21 3
5 5 66 1 1 24 3
5 6 61 1 1 27 3
5 3 35 1 1 30 3
5 4 32 1 1 33 3
5 1 28 1 1 36 3
5 1 27 1 1 39 3
5 6 26 1 1 42 3
5 5 26 1 1 45 3
5 0 26 1 1 48 3
6 22 179 2 1 3 3
6 13 151 2 1 6 3
6 3 138 2 1 9 3
6 5 135 2 1 12 3
6 1 130 2 1 15 3
6 5 104 2 1 18 3
6 7 99 2 1 21 3
6 6 92 2 1 24 3
6 6 66 2 1 27 3
6 7 60 2 1 30 3
6 4 53 2 1 33 3
6 0 30 2 1 36 3
6 2 29 2 1 39 3
6 3 27 2 1 42 3
6 1 24 2 1 45 3
6 0 16 2 1 48 3
6 1 15 2 1 51 3
6 0 14 2 1 54 3
6 1 14 2 1 57 3
6 0 14 2 1 60 3
6 13 178 1 1 3 3
6 7 160 1 1 6 3
6 7 153 1 1 9 3
6 10 146 1 1 12 3
6 10 136 1 1 15 3
6 2 97 1 1 18 3
6 5 95 1 1 21 3
6 3 90 1 1 24 3
6 5 57 1 1 27 3
6 2 52 1 1 30 3
6 6 50 1 1 33 3
6 3 37 1 1 36 3
6 1 34 1 1 39 3
6 2 33 1 1 42 3
6 4 31 1 1 45 3
6 0 17 1 1 48 3
6 0 17 1 1 51 3
6 1 17 1 1 54 3
6 0 16 1 1 57 3
6 0 16 1 1 60 3
END
ts[] bs[]
3 1
2 1
2 1
2 1
2 1
2 1
END
给您的一些建议,希望对您有所帮助:
我认为您的问题可以在 https://stats.stackexchange.com/ 中得到更好的回答;尽管如果您想将问题移到那里,您应该重写它而不是发布代码转储。
WinBUGS已经好几年了,距离上次更新已经8年了!我想你应该忘了它,因为有几个更好的选择。
您可以在 Jags or Stan, in which both are usable in R via rJags and RStan 中尝试几乎相同的代码。 Stan 特别重要,因为它使用 HCMC,它在许多情况下收敛,而 WinBUGS 不收敛。
我认为你应该阅读通俗易懂的书:Doing Bayesian Data Analysis 2e 来自 John K. Kruschke 至能够自己理解并进行贝叶斯数据分析
最终,我无法在 WinBUGS 中将模型正确地 运行。然而,我能够将模型移植到 STAN 并获得非常接近的匹配。请参阅下面的 STAN 代码:
data {
int<lower=1> N;
int<lower=1> NS;
int<lower=1> NT;
cov_matrix[2] prec2;
matrix[2,2] R;
vector[2] means;
int<lower=0> bs[NS];
int<lower=0> ts[NS];
int<lower=0> s[N];
int<lower=0> t[N];
int<lower=0> n[N];
int<lower=0> r[N];
real<lower=0> dt[N];
real<lower=0> time[N];
}
parameters {
vector[2] mu[NS];
vector[2] delta[NS];
vector[2] dj[NT-1];
cov_matrix[2] omega;
}
transformed parameters{
real<lower=0,upper=1> p[N];
real<lower=0> h[N];
real nu[N];
real theta[N];
vector[2] md[NS];
vector[2] d[NT];
d[1][1] <- 0;
d[1][2] <- 0;
for(j in 2:NT){
d[j] <- dj[j-1];
}
for(k in 1:NS){
md[k] <- d[ts[k]] - d[bs[k]];
}
for(i in 1:N){
if(t[i] == 1){
nu[i] <- mu[s[i]][1];
theta[i] <- mu[s[i]][2];
}else{
nu[i] <- mu[s[i]][1] + delta[s[i]][1];
theta[i] <- mu[s[i]][2] + delta[s[i]][2];
}
h[i] <- exp(nu[i] + log(time[i]) * theta[i]);
p[i] <- 1 - exp(- h[i] * dt[i]);
}
}
model {
omega ~ inv_wishart(2,R);
for(j in 1:(NT-1)){
dj[j] ~ multi_normal(means,prec2);
}
for(k in 1:NS){
delta[k] ~ multi_normal(md[k],omega);
mu[k] ~ multi_normal(means,prec2);
}
for(i in 1:N){
r[i] ~ binomial(n[i],p[i]);
}
}
generated quantities{
vector[N] log_lik;
for (l in 1:N) {
log_lik[l] <- binomial_log(r[l], n[l], p[l]);
}
}
请注意,由于 STAN/BUGS 之间的差异,对 precision/variance 的混淆会使输入数据的内容变得混乱。对于 R 和 prec2,我加载了:
dataList$R <- matrix(c(0.01,0,0,0.01),nrow=2,ncol=2,byrow=TRUE)
dataList$prec2 <- matrix(c(10^4,0,0,10^4),nrow=2,ncol=2,byrow=TRUE)
我目前正在对多个临床试验的生存数据进行荟萃分析。
为此,我使用相同的方法从已发布的分析中获取代码。但是,当 运行 使用已发布分析中的数据使用此代码时,我无法复制他们的结果。事实上,结果无法收敛到任何一种合理的估计。
代码本身(不包括数据)应该是正确的,因为它直接来自作者。我认为问题必须与初始值或 采样的参数是 运行,但在玩了很多之后 初始值、老化时间、变薄等...我还没有得到有意义的结果。
如果有人能就如何 运行 这个(初始值等...)正确地 运行 提出建议,我将不胜感激。或者,如果代码中存在问题或者数据的设置方式与代码不匹配,了解这些信息将很有用。
附带说明一下,我正在使用 R2WinBUGs 进行分析,尽管我有 单独使用 WinBUGs 遇到了同样的问题。
A bit of extra background on the method:
The way this works is by estimating the difference in shape and scale parameters of a reparameterized Weibull distribution between treatments across multiple studies using random effects.
The Weibull distribution is reparameterized such that log of the hazard rate is a+b*log(t) where a is a scale parameter and b is a shape parameter. From this, you can calculate the the likelihood function of a given number of failures out of a given number of patients over an interval.
Unfortunately, the article is public, but if you can get access here is the link: http://onlinelibrary.wiley.com/doi/10.1002/jrsm.25/abstract;jsessionid=2BA8F0D9BEF9A33F84975618D33F8DD9.f03t03?userIsAuthenticated=false&deniedAccessCustomisedMessage=
输入模型的变量的快速摘要:
NT: Number of separate treatments included.
N: Number of rows in the main dataset. NS: Number of studies
s: Study that the line of data corresponds to (this is numbered 1:6)
r: number of patients failing in interval for this treatment/study
n: number of patients at risk at start of interval for this treatment/study
t: Treatment that this line of data corresponds to (numbered 1:3)
b: Indicates which treatment is the baseline one to which others are compared (set to 1 for each line).
bs: Treatment that is the control arm of this study
bt: Treatment that is the research arm of this study
WinBUGS代码(含数据):
#Winbugs code for random effects networks meta-analysis model
Model
{
for (i in 1:N)
{ # N=number of data points in dataset
#likelihood
r[i]~ dbin(p[i],n[i])
p[i]<-1-exp(-h[i]*dt[i]) # hazard h over interval [t,t+dt] # expressed as deaths per unit person-time (e.g. months)
#random effects model
log(h[i])<-nu[i]+log(time[i])*theta[i]
nu[i]<-mu[s[i],1]+delta[s[i],1]*(1-equals(t[i],b[i]))
theta[i]<-mu[s[i],2]+ delta[s[i],2]*(1-equals(t[i],b[i]))
}
for(k in 1 :NS)
{ # NS=number of studies in dataset
delta[k,1:2]~dmnorm(md[k,1:2],omega[1:2,1:2])
md[k,1]<-d[ts[k],1]-d[bs[k],1]
md[k,2]<-d[ts[k],2]-d[bs[k],2]
}
# priors
d[1,1]<-0
d[1,2]<-0
for(j in 2 :NT)
{ # NT=number of treatments
d[j,1:2] ~ dmnorm(mean[1:2],prec2[,])
}
for(k in 1 :NS)
{
mu[k,1:2] ~ dmnorm(mean[1:2],prec2[,])
}
omega[1:2, 1:2] ~ dwish(R[1:2,1:2],2)
}
# Winbugs data set
list(N=242, NS=6, NT=3,
mean=c(0,0),
prec2 = structure(.Data = c(
0.0001,0,
0,0.0001), .Dim = c(2,2)),
R = structure(.Data = c(
0.01,0,
0,0.01), .Dim = c(2,2))
)
s[] r[] n[] t[] b[] time[] dt[]
1 15 152 3 1 3 3
1 11 140 3 1 6 3
1 8 129 3 1 9 3
1 9 121 3 1 12 3
1 9 112 3 1 15 3
1 3 83 3 1 18 3
1 4 80 3 1 21 3
1 5 76 3 1 24 3
1 2 71 3 1 27 3
1 2 41 3 1 30 3
1 1 39 3 1 33 3
1 3 38 3 1 36 3
1 2 35 3 1 39 3
1 1 33 3 1 42 3
1 3 32 3 1 45 3
1 3 29 3 1 48 3
1 2 26 3 1 51 3
1 1 24 3 1 54 3
1 1 23 3 1 57 3
1 1 22 3 1 60 3
1 10 149 1 1 3 3
1 11 140 1 1 6 3
1 9 128 1 1 9 3
1 5 119 1 1 12 3
1 6 114 1 1 15 3
1 3 72 1 1 18 3
1 5 70 1 1 21 3
1 4 65 1 1 24 3
1 7 61 1 1 27 3
1 2 34 1 1 30 3
1 2 32 1 1 33 3
1 3 30 1 1 36 3
1 2 27 1 1 39 3
1 2 25 1 1 42 3
1 1 23 1 1 45 3
1 2 22 1 1 48 3
1 1 19 1 1 51 3
1 2 19 1 1 54 3
1 1 17 1 1 57 3
1 0 16 1 1 60 3
2 4 125 2 1 3 3
2 4 121 2 1 6 3
2 2 117 2 1 9 3
2 5 114 2 1 12 3
2 2 109 2 1 15 3
2 3 107 2 1 18 3
2 2 104 2 1 21 3
2 4 94 2 1 24 3
2 4 90 2 1 27 3
2 3 81 2 1 30 3
2 4 78 2 1 33 3
2 3 61 2 1 36 3
2 5 58 2 1 39 3
2 1 48 2 1 42 3
2 2 47 2 1 45 3
2 3 41 2 1 48 3
2 0 38 2 1 51 3
2 3 29 2 1 54 3
2 3 26 2 1 57 3
2 2 18 2 1 60 3
2 0 16 2 1 63 3
2 1 10 2 1 66 3
2 0 9 2 1 69 3
2 0 3 2 1 72 3
2 0 3 2 1 75 3
2 0 3 2 1 78 3
2 15 196 1 1 3 3
2 9 179 1 1 6 3
2 10 170 1 1 9 3
2 9 162 1 1 12 3
2 9 153 1 1 15 3
2 5 141 1 1 18 3
2 5 136 1 1 21 3
2 10 121 1 1 24 3
2 5 111 1 1 27 3
2 7 92 1 1 30 3
2 7 85 1 1 33 3
2 4 71 1 1 36 3
2 6 67 1 1 39 3
2 4 53 1 1 42 3
2 5 49 1 1 45 3
2 6 36 1 1 48 3
2 3 30 1 1 51 3
2 2 26 1 1 54 3
2 2 24 1 1 57 3
2 0 13 1 1 60 3
2 1 13 1 1 63 3
2 1 11 1 1 66 3
2 1 10 1 1 69 3
2 0 6 1 1 72 3
2 0 6 1 1 75 3
2 0 6 1 1 78 3
3 6 113 2 1 3 3
3 4 105 2 1 6 3
3 3 101 2 1 9 3
3 1 97 2 1 12 3
3 9 96 2 1 15 3
3 4 84 2 1 18 3
3 2 80 2 1 21 3
3 4 74 2 1 24 3
3 3 70 2 1 27 3
3 2 59 2 1 30 3
3 0 57 2 1 33 3
3 6 51 2 1 36 3
3 2 45 2 1 39 3
3 1 37 2 1 42 3
3 3 36 2 1 45 3
3 1 27 2 1 48 3
3 1 26 2 1 51 3
3 2 25 2 1 54 3
3 7 116 1 1 3 3
3 6 111 1 1 6 3
3 4 105 1 1 9 3
3 3 99 1 1 12 3
3 9 96 1 1 15 3
3 5 85 1 1 18 3
3 5 80 1 1 21 3
3 3 68 1 1 24 3
3 7 65 1 1 27 3
3 8 48 1 1 30 3
3 4 40 1 1 33 3
3 2 33 1 1 36 3
3 0 31 1 1 39 3
3 1 28 1 1 42 3
3 2 27 1 1 45 3
3 3 20 1 1 48 3
3 1 17 1 1 51 3
3 0 16 1 1 54 3
4 10 167 2 1 3 3
4 5 149 2 1 6 3
4 6 145 2 1 9 3
4 3 138 2 1 12 3
4 4 135 2 1 15 3
4 5 128 2 1 18 3
4 2 122 2 1 21 3
4 2 120 2 1 24 3
4 7 104 2 1 27 3
4 9 98 2 1 30 3
4 5 89 2 1 33 3
4 2 57 2 1 36 3
4 2 55 2 1 39 3
4 4 53 2 1 42 3
4 2 49 2 1 45 3
4 2 26 2 1 48 3
4 1 24 2 1 51 3
4 1 23 2 1 54 3
4 1 11 2 1 57 3
4 0 10 2 1 60 3
4 0 10 2 1 63 3
4 2 164 1 1 3 3
4 5 153 1 1 6 3
4 7 148 1 1 9 3
4 6 141 1 1 12 3
4 12 135 1 1 15 3
4 6 119 1 1 18 3
4 4 113 1 1 21 3
4 3 109 1 1 24 3
4 5 98 1 1 27 3
4 2 94 1 1 30 3
4 2 92 1 1 33 3
4 4 55 1 1 36 3
4 3 50 1 1 39 3
4 1 48 1 1 42 3
4 2 47 1 1 45 3
4 1 22 1 1 48 3
4 1 21 1 1 51 3
4 0 20 1 1 54 3
4 1 7 1 1 57 3
4 0 6 1 1 60 3
4 0 6 1 1 63 3
5 12 152 2 1 3 3
5 7 135 2 1 6 3
5 9 128 2 1 9 3
5 8 120 2 1 12 3
5 7 112 2 1 15 3
5 1 77 2 1 18 3
5 3 76 2 1 21 3
5 2 73 2 1 24 3
5 4 71 2 1 27 3
5 5 45 2 1 30 3
5 3 40 2 1 33 3
5 2 37 2 1 36 3
5 3 35 2 1 39 3
5 3 32 2 1 42 3
5 3 32 2 1 45 3
5 1 32 2 1 48 3
5 9 149 1 1 3 3
5 4 131 1 1 6 3
5 5 127 1 1 9 3
5 8 122 1 1 12 3
5 11 114 1 1 15 3
5 5 76 1 1 18 3
5 5 71 1 1 21 3
5 5 66 1 1 24 3
5 6 61 1 1 27 3
5 3 35 1 1 30 3
5 4 32 1 1 33 3
5 1 28 1 1 36 3
5 1 27 1 1 39 3
5 6 26 1 1 42 3
5 5 26 1 1 45 3
5 0 26 1 1 48 3
6 22 179 2 1 3 3
6 13 151 2 1 6 3
6 3 138 2 1 9 3
6 5 135 2 1 12 3
6 1 130 2 1 15 3
6 5 104 2 1 18 3
6 7 99 2 1 21 3
6 6 92 2 1 24 3
6 6 66 2 1 27 3
6 7 60 2 1 30 3
6 4 53 2 1 33 3
6 0 30 2 1 36 3
6 2 29 2 1 39 3
6 3 27 2 1 42 3
6 1 24 2 1 45 3
6 0 16 2 1 48 3
6 1 15 2 1 51 3
6 0 14 2 1 54 3
6 1 14 2 1 57 3
6 0 14 2 1 60 3
6 13 178 1 1 3 3
6 7 160 1 1 6 3
6 7 153 1 1 9 3
6 10 146 1 1 12 3
6 10 136 1 1 15 3
6 2 97 1 1 18 3
6 5 95 1 1 21 3
6 3 90 1 1 24 3
6 5 57 1 1 27 3
6 2 52 1 1 30 3
6 6 50 1 1 33 3
6 3 37 1 1 36 3
6 1 34 1 1 39 3
6 2 33 1 1 42 3
6 4 31 1 1 45 3
6 0 17 1 1 48 3
6 0 17 1 1 51 3
6 1 17 1 1 54 3
6 0 16 1 1 57 3
6 0 16 1 1 60 3
END
ts[] bs[]
3 1
2 1
2 1
2 1
2 1
2 1
END
给您的一些建议,希望对您有所帮助:
我认为您的问题可以在 https://stats.stackexchange.com/ 中得到更好的回答;尽管如果您想将问题移到那里,您应该重写它而不是发布代码转储。
WinBUGS已经好几年了,距离上次更新已经8年了!我想你应该忘了它,因为有几个更好的选择。
您可以在 Jags or Stan, in which both are usable in R via rJags and RStan 中尝试几乎相同的代码。 Stan 特别重要,因为它使用 HCMC,它在许多情况下收敛,而 WinBUGS 不收敛。
我认为你应该阅读通俗易懂的书:Doing Bayesian Data Analysis 2e 来自 John K. Kruschke 至能够自己理解并进行贝叶斯数据分析
最终,我无法在 WinBUGS 中将模型正确地 运行。然而,我能够将模型移植到 STAN 并获得非常接近的匹配。请参阅下面的 STAN 代码:
data {
int<lower=1> N;
int<lower=1> NS;
int<lower=1> NT;
cov_matrix[2] prec2;
matrix[2,2] R;
vector[2] means;
int<lower=0> bs[NS];
int<lower=0> ts[NS];
int<lower=0> s[N];
int<lower=0> t[N];
int<lower=0> n[N];
int<lower=0> r[N];
real<lower=0> dt[N];
real<lower=0> time[N];
}
parameters {
vector[2] mu[NS];
vector[2] delta[NS];
vector[2] dj[NT-1];
cov_matrix[2] omega;
}
transformed parameters{
real<lower=0,upper=1> p[N];
real<lower=0> h[N];
real nu[N];
real theta[N];
vector[2] md[NS];
vector[2] d[NT];
d[1][1] <- 0;
d[1][2] <- 0;
for(j in 2:NT){
d[j] <- dj[j-1];
}
for(k in 1:NS){
md[k] <- d[ts[k]] - d[bs[k]];
}
for(i in 1:N){
if(t[i] == 1){
nu[i] <- mu[s[i]][1];
theta[i] <- mu[s[i]][2];
}else{
nu[i] <- mu[s[i]][1] + delta[s[i]][1];
theta[i] <- mu[s[i]][2] + delta[s[i]][2];
}
h[i] <- exp(nu[i] + log(time[i]) * theta[i]);
p[i] <- 1 - exp(- h[i] * dt[i]);
}
}
model {
omega ~ inv_wishart(2,R);
for(j in 1:(NT-1)){
dj[j] ~ multi_normal(means,prec2);
}
for(k in 1:NS){
delta[k] ~ multi_normal(md[k],omega);
mu[k] ~ multi_normal(means,prec2);
}
for(i in 1:N){
r[i] ~ binomial(n[i],p[i]);
}
}
generated quantities{
vector[N] log_lik;
for (l in 1:N) {
log_lik[l] <- binomial_log(r[l], n[l], p[l]);
}
}
请注意,由于 STAN/BUGS 之间的差异,对 precision/variance 的混淆会使输入数据的内容变得混乱。对于 R 和 prec2,我加载了:
dataList$R <- matrix(c(0.01,0,0,0.01),nrow=2,ncol=2,byrow=TRUE)
dataList$prec2 <- matrix(c(10^4,0,0,10^4),nrow=2,ncol=2,byrow=TRUE)