Geek Answers Handbook

WinBUGS Weibull Network Meta-Analysis

I am currently working on a meta-analysis of survival data in several clinical trials.

To do this, I have code from a published analysis using the same methodology. However, when I run this code using data from a published analysis, I cannot reproduce their results. In fact, the results do not converge to any reasonable estimate.

The code itself (not including data) must be correct, since it comes directly from the authors. I assume that the problem should consist of initial values ​​or parameters of how the sampling is performed, but after playing with a lot of initial values, burn duration, thinning, etc. I did not get significant results.

I would appreciate any suggestions on how to run this (initial values, etc.) to make it work properly. Otherwise, if there is a problem in the code or if the data is configured in a way that does not match the code, it would be useful to know.

As a side note, I am doing an analysis using R2WinBUG, although I have got the same problems using only WinBUG.

A bit of extra background for the method:

How it works, evaluating the difference in the shape and scale of the parameters of the Weibull reparameterized distribution between treatments in several studies using random effects.

The Wybull distribution is reparametrized so that the log hazard level is + b * log (t), where a is the scale parameter and b is the shape parameter. From this, you can calculate the probability function of a given number of failures from a given number of patients at intervals.

Unfortunately, the article is publicly available, but if you can access it here is the link: http://onlinelibrary.wiley.com/doi/10.1002/jrsm.25/abstract;jsessionid=2BA8F0D9BEF9A33F84975618D33F8DD9.f03t03?userIsAuthenticated=false&deniedessccess

A brief summary of the variables introduced into the model:

NT: number of individual treatments included.

N: number of rows in the main data set. NS: Number of studies

s: Examine the correspondence of the data line (this is numbered 1: 6)

r: number of patients under the interval for this treatment / study

n: number of patients at risk at the start of the interval for this treatment / study

t: processing corresponding to this data line (numbered 1: 3)

b: indicates which treatment is the base with which others are compared (1 for each row).

bs: Treatment being the control group of this study

bt: The treatment, which is the research unit of this study

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 (eg 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
+8
bayesian winbugs survival-analysis winbugs14 r2winbugs
source share
2 answers

In the end, I was unable to get the model to work correctly in WinBUGS. However, I was able to transfer the model to STAN and get a very close match. The following is 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]); } }

Please note that due to differences between STAN / BUGS, confusion regarding accuracy / variance can be confusing for data entry. For R and prec2, I downloaded:

 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)
+2
source share

Some suggestions for you that I hope will be helpful:

  • I think your question could be better answered https://stats.stackexchange.com/ ; Although if you want to move the question there, you should rewrite it, and not post a dump of code.

  • WinBUGS is already several years old, and 8 years have passed since its last update! I think you should forget this because there are some better alternatives.

  • You can try almost the same code in Jags or Stan , in which both can be used in R via rJags and RStan . Stan is especially important because he uses HCMC, which converges in many situations that are not in WinBUGS.

  • I think you should read the easy book: Performing Bayesian 2e Data Analysis by John K. Kruschke to be able to understand and do Bayesian data analysis yourself.

+1
source share

More articles:


All Articles