WinBUGS Weibull网络元分析 [英] WinBUGS Weibull Network Meta-Analysis
问题描述
我目前正在进行多项临床试验中生存数据的荟萃分析.
I am currently working on a meta-analysis of survival data across several clinical trials.
为此，我使用相同的方法从发布的分析中获取代码.但是，当使用已发布的分析中的数据运行此代码时，我无法复制其结果.实际上，结果无法收敛到任何合理的估计.
To do this, I have code from a published analysis using the same methodology. However, when running this code using the data from the published analysis, I am unable to replicate their results. In fact, the results fail to converge to any kind of reasonable estimate.
代码本身(不包括数据)应正确无误，因为它直接来自作者.我认为问题必须与初始值或 采样运行方式的参数，但在多次播放之后 初始值，老化时间，变薄等...我没有得到有意义的结果.
The code itself (not including the data) should be correct as it comes directly from the authors. I assume the problem has to do w/ initial values or parameters of how the sampling is run, but after playing w/ many initial values, the length of burn in, thinning, etc... I haven't gotten meaningful results.
任何人都希望获得有关如何运行此命令(初始值等)以使其正常运行的建议.或者，如果代码中有问题，或者如果以与代码不匹配的方式设置数据，那么了解这一点将很有用.
I would appreciate anyone suggestions about how to run this (initial values, etc...) to get it to run properly. Alternately, if there are problems in the code or if the data is set up in a way that doesn't match the code, that would be useful to know.
作为旁注，尽管我有R2WinBUG，但仍在进行分析 单独使用WinBUG会遇到同样的问题.
As a side note, I am doing the analyses using R2WinBUGs, though I have gotten the same kind of problems using WinBUGs alone.
方法的一些额外背景:
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.
重新对Weibull分布进行参数化，以使对数的对数 危险率是a + b * log(t)，其中a是比例参数，b是a 形状参数.由此，您可以计算出可能性 给定数量的故障中的给定数量的故障的功能 间隔的病人.
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.
不幸的是，该文章是公开的，但是如果您可以在此处访问 是链接: http://onlinelibrary.wiley./10.1002/jrsm.25/abstract;jsessionid=2BA8F0D9BEF9A33F84975618D33F8DD9.f03t03?userIsAuthenticated=false&deniedAccessCustomisedMessage=
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=
输入模型的变量的快速摘要:
A quick summary of the variables entered into the model:
NT:包含的单独治疗次数.
NT: Number of separate treatments included.
N:主数据集中的行数. NS:研究数量
N: Number of rows in the main dataset. NS: Number of studies
s:研究数据行所对应的(编号为1:6)
s: Study that the line of data corresponds to (this is numbered 1:6)
r:此治疗/研究间隔失败的患者人数
r: number of patients failing in interval for this treatment/study
n:间隔开始时处于风险中的患者数量 治疗/研究
n: number of patients at risk at start of interval for this treatment/study
t:此数据行对应的处理方式(编号为1:3)
t: Treatment that this line of data corresponds to (numbered 1:3)
b:指示哪种治疗是与其他治疗进行比较的基准(每行设置为1).
b: Indicates which treatment is the baseline one to which others are compared (set to 1 for each line).
bs:治疗是这项研究的控制臂
bs: Treatment that is the control arm of this study
bt:治疗是这项研究的研究重点
bt: Treatment that is the research arm of this study
WinBUGS代码(包括数据):
WinBUGS code (including data):
#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
推荐答案
最终，我无法使模型在WinBUGS中正常运行.但是，我能够将模型移植到STAN并获得非常接近的匹配.请参见下面的STAN代码:
Ultimately, I was unable to get the model to run properly in WinBUGS. However, I was able to port the model over to STAN and get a very close match. See below for the STAN code:
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之间的差异，对精度/方差的混淆可能会使输入数据混乱.对于R和prec2，我加载了:
Note that due to differences between STAN/BUGS, confusion over precision/variance can make it confusing what to enter for data. For R and prec2, I loaded:
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)
这篇关于WinBUGS Weibull网络元分析的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
