WinBUGS Weibull 网络元分析

Posted 2023-02-19

【中文标题】WinBUGS Weibull 网络元分析

【英文标题】：WinBUGS Weibull Network Meta-Analysis

【发布时间】：2015-11-10 12:58:45

【问题描述】：

我目前正在对多项临床试验的生存数据进行荟萃分析。

为此，我使用相同方法从已发布的分析中获取代码。但是，当使用已发布分析中的数据运行此代码时，我无法复制他们的结果。事实上，结果无法收敛到任何合理的估计。

代码本身（不包括数据）应该是正确的，因为它直接来自作者。我认为问题必须与初始值或 采样如何运行的参数，但在玩了很多之后 初始值、老化时间、变薄等...我没有得到有意义的结果。

对于如何运行它（初始值等）以使其正常运行的任何建议，我将不胜感激。或者，如果代码中存在问题，或者数据的设置方式与代码不匹配，那么了解这些会很有用。

附带说明，我正在使用 R2WinBUGs 进行分析，尽管我有 单独使用 WinBUG 也会遇到同样的问题。

该方法的一些额外背景：

其工作方式是通过估计形状和规模的差异 重新参数化的 Weibull 分布的参数 使用随机效应在多项研究中进行治疗。

Weibull 分布被重新参数化，使得 危险率是 a+b*log(t)，其中 a 是尺度参数，b 是 a 形状参数。由此，您可以计算可能性 给定数量的故障中给定数量的功能 一段时间内的患者。

不幸的是，这篇文章是公开的，但如果你可以在这里访问 是链接： http://onlinelibrary.wiley.com/doi/10.1002/jrsm.25/abstract;jsessionid=2BA8F0D9BEF9A33F84975618D33F8DD9.f03t03?userIsAuthenticated=false&deniedAccessCustomisedMessage=

输入模型的变量的快速总结：

NT：包括的单独治疗的数量。

N：主数据集中的行数。 NS：研究数量

s：研究该行数据对应的（编号为1:6）

r：此治疗/研究间隔内失败的患者人数

n：在此间隔开始时有风险的患者数量 治疗/研究

t：这行数据对应的处理（编号1:3）

b：指示哪个治疗是基线，其他治疗与之比较（每行设置为 1）。

bs：作为本研究对照组的治疗

bt：治疗是本研究的研究分支

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

【参考方案1】：

给你一些建议，希望对你有帮助：

我认为https://stats.stackexchange.com/ 可以更好地回答您的问题；虽然如果你想把问题移到那里，你应该重写它而不是发布代码转储。

WinBUGS已经好几年了，距离上次更新已经8年了！我认为你应该忘记它，因为有几个更好的选择。

您可以在Jags 或Stan 中尝试几乎相同的代码，其中两者都可以通过rJags 和RStan 在R 中使用。 Stan 特别重要，因为它使用了 HCMC，它在 WinBUGS 没有的许多情况下收敛。

我认为您应该阅读 John K. Kruschke 的简单书：Doing Bayesian Data Analysis 2e，以便能够自己理解和进行贝叶斯数据分析。

【讨论】：

您好，谢谢您的评论。我实际上已经用 rstan 重写了模型，它与已发布的结果相匹配，差异很小。同意 WinBUGS 已过时，希望更多的学者能做出改变。

【参考方案2】：

最终，我无法让模型在 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 之间的差异，对精度/方差的混淆可能会导致混淆输入数据的内容。对于 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)