{"id":198,"date":"2014-07-25T09:07:56","date_gmt":"2014-07-25T09:07:56","guid":{"rendered":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/?p=198"},"modified":"2025-02-26T13:21:38","modified_gmt":"2025-02-26T13:21:38","slug":"network-meta-analysis-with-winbugs-and-stata","status":"publish","type":"post","link":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2014\/07\/25\/network-meta-analysis-with-winbugs-and-stata\/","title":{"rendered":"Network Meta-analysis with WinBUGS and Stata"},"content":{"rendered":"

In my previous posting (https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2014\/07\/18\/network-meta-analysis-in-stata\/<\/a>) I introduced an example of Network Meta-analysis taken from a NICE report that estimated the incidence of diabetes amongst people being treated for high blood pressure by one of six classes of drug. In that posting I analysed the data using Stata\u2019s glm command and now I should like to analyse the same data using Stata and WinBUGS.<\/p>\n

Programming the model in WinBUGS requires a little ingenuity because the trials do not include all of the treatments\u00a0with some trials comparing two treatments and others comparing three. I will\u00a0follow the NICE report in using d[1,\u2026,6] to represent the treatment effects and mu[1,\u2026,22] to represent study effects. Treatment 1 (diuretics) is chosen as the baseline, so d[1] is set to zero.<\/p>\n

In the report vague priors are used throughout which negates one of the main benefits of a Bayesian analysis. However meta-analysis presents a particular problem for Bayesians\u00a0because it is so difficult to produce informative priors that are not influenced, directly or indirectly, by the studies that we want to combine.<\/p>\n

Once we have stacked the data as described in the previous posting, we can program the Bayesian analysis\u00a0in Stata by preparing a do file as follows,<\/p>\n

tomodel diabetes_fixed.do model.txt<\/span><\/strong>\r\n\/*<\/span><\/strong>\r\nmodel{ <\/span><\/strong>\r\n\u00a0\u00a0 for(i in 1:48){ <\/span><\/strong>\r\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 r[i] ~ dbin(p[i],n[i]) <\/span><\/strong>\r\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 cloglog(p[i]) <- mu[id[i]]+d[t[i]]<\/span><\/strong>\r\n\u00a0\u00a0 }<\/span><\/strong>\r\n\u00a0\u00a0 d[1]<-0 <\/span><\/strong>\r\n\u00a0\u00a0 for(k in 2:6){ d[k] ~ dnorm(0,0.0001) } <\/span><\/strong>\r\n\u00a0\u00a0 for(k in 1:22){ mu[k] ~ dnorm(0,0.0001) } <\/span><\/strong>\r\n}<\/span><\/strong>\r\n*\/<\/span><\/strong>\r\nwbslist (vect id t r n , for(\"%6.0f\")) using data.txt , replace<\/span><\/strong>\r\nwbslist (mu=c(22{0}),d=c(NA,5{0})) using init.txt , replace\u00a0<\/span><\/strong>\r\nwbsscript using script.txt , replace\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \/\/\/<\/span><\/strong>\r\n\u00a0\u00a0 model(model.txt) data(data.txt) init(init.txt) \/\/\/<\/span><\/strong>\r\n\u00a0\u00a0 log(log.txt) coda(diab) burn(20000) update(20000) thin(2) set(mu d) dic<\/span><\/strong>\r\nwbsrun using script.txt <\/span><\/strong>\r\ntype log.txt<\/span><\/strong>\r\nwbscoda using diab\u00a0 , clear <\/span><\/strong>\r\nsave diabetes_fit.dta, replace<\/span><\/strong><\/pre>\n
mcmcstats d_2 d_3 d_4 d_5 d_6<\/strong><\/span><\/pre>\n
–wbsmodel<\/strong>– is used to write the model to a text file (see my posting https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2014\/03\/21\/keeping-good-records\/ <\/a>\u00a0). \u2013wbslist<\/strong>– is used to write the data and initial values to a file. Notice\u00a0how the first element of d is initialized to be missing (NA in WinBUGS) since it does not have a prior, but rather it has its value set to zero within the program.<\/p>\n
The results calculated by WinBUGS are,<\/p>\n
----------------------------------------------------------------------------<\/strong><\/pre>\n
Parameter\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 n\u00a0\u00a0\u00a0\u00a0 mean\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 sd\u00a0\u00a0\u00a0\u00a0 sem\u00a0\u00a0 median\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 95% CrI<\/strong><\/pre>\n
----------------------------------------------------------------------------<\/strong><\/pre>\n
 d_2\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 -0.247\u00a0\u00a0 0.057\u00a0\u00a0 0.0014\u00a0\u00a0 -0.247 ( -0.360, -0.133 )<\/strong><\/pre>\n
 d_3\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 -0.056\u00a0\u00a0 0.056\u00a0\u00a0 0.0017\u00a0\u00a0 -0.057 ( -0.163,\u00a0\u00a0 0.053 )<\/strong><\/pre>\n
 d_4\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 -0.253\u00a0\u00a0 0.053\u00a0\u00a0 0.0017\u00a0\u00a0 -0.254 ( -0.357, -0.148 )<\/strong><\/pre>\n
 d_5\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 -0.359\u00a0\u00a0 0.054\u00a0\u00a0 0.0013\u00a0\u00a0 -0.359 ( -0.464, -0.252 )<\/strong><\/pre>\n
 d_6\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 -0.453\u00a0\u00a0 0.063\u00a0\u00a0 0.0015\u00a0\u00a0 -0.452 ( -0.575, -0.330 )<\/strong><\/pre>\n
----------------------------------------------------------------------------<\/strong><\/pre>\n
For comparison here are the results from glm as reported in the previous posting<\/p>\n
------------------------------------------------------------------------------<\/strong><\/pre>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 |\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 OIM<\/strong><\/pre>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 r |\u00a0\u00a0\u00a0\u00a0 Coef.\u00a0\u00a0 Std. Err.\u00a0\u00a0\u00a0\u00a0 z\u00a0\u00a0 P>|z|\u00a0\u00a0\u00a0\u00a0 [95% Conf. Interval]<\/strong><\/pre>\n
-------------+----------------------------------------------------------------<\/strong><\/pre>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a02 | -.2469671\u00a0\u00a0 .0562828\u00a0\u00a0 -4.39\u00a0\u00a0 0.000\u00a0\u00a0 -.3572794\u00a0\u00a0 -.1366548<\/strong><\/pre>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0  3 |\u00a0\u00a0 -.056982\u00a0\u00a0 .0557689\u00a0\u00a0 -1.02\u00a0\u00a0 0.307\u00a0\u00a0 -.1662869\u00a0\u00a0\u00a0\u00a0 .052323<\/strong><\/pre>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0  4 | -.2530259\u00a0\u00a0 .0536676\u00a0\u00a0 -4.71\u00a0\u00a0 0.000\u00a0\u00a0 -.3582125\u00a0\u00a0 -.1478392<\/strong><\/pre>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0  5 |\u00a0\u00a0 -.358473 \u00a0.0533383\u00a0\u00a0 -6.72\u00a0\u00a0 0.000\u00a0\u00a0 -.4630142\u00a0\u00a0 -.2539318<\/strong><\/pre>\n
\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0  6 | -.4527086\u00a0\u00a0 .0630479\u00a0\u00a0 -7.18\u00a0\u00a0 0.000\u00a0\u00a0 -.5762803\u00a0\u00a0 -.3291369<\/strong><\/pre>\n
------------------------------------------------------------------------------<\/strong><\/pre>\n
These are effectively equivalent as one would expect.<\/p>\n
We should, of course, check convergence, but there does not seem to be a problem. For instance, the Geweke p-values for the 5 treatment effects are all non-significant (0.18, 0.19, 0.21, 0.10, 0.28 respectively). As a graphical summary, here are the mahalanobis plots for the 5 treatments differences and the 22 study effect. Neither indicates a problem with convergence.<\/p>\n
mcmcmahal d_2 d_3 d_4 d_5 d_6<\/p>\n
\"NM2mahal1\"<\/a><\/p>\n
mcmcmahal mu_1-mu_22<\/p>\n
\"NM2mahal2\"<\/a><\/p>\n
Convergence does not imply that we have a good model, so it makes sense to look at the residuals. When calculating the residuals the chain is thinned from length 20000 to length 1000 to avoid needlessly long calculations when constructing the residual plots.<\/p>\n
I have chosen to calculate the residuals on the scale of r the number of cases of diabetes. In the model file below the vector f[] contains the predicted number of cases and these are read into Stata as f_1, f_2,\u2026<\/p>\n
wbsmodel diabetes_fixed.do modelf.txt<\/strong><\/span>\r\n \/*<\/span><\/strong>\r\n model{<\/span><\/strong>\r\n for(i in 1:48){<\/span><\/strong>\r\n r[i] ~ dbin(p[i],n[i])<\/span><\/strong>\r\n cloglog(p[i]) <- mu[id[i]]+d[t[i]]<\/span><\/strong>\r\n f[i] ~ dbin(p[i],n[i])<\/span><\/strong>\r\n }<\/span><\/strong>\r\n d[1]<-0<\/span><\/strong>\r\n for(k in 2:6){ d[k] ~ dnorm(0,0.0001) }<\/span><\/strong>\r\n for(k in 1:22){ mu[k] ~ dnorm(0,0.0001) }<\/span><\/strong>\r\n }<\/span><\/strong>\r\n *\/<\/span><\/strong>\r\n wbsscript using script.txt , replace \/\/\/<\/span><\/strong>\r\n model(modelf.txt) data(data.txt) init(init.txt) \/\/\/<\/span><\/strong>\r\n log(log.txt) coda(diab) burn(20000) update(20000) \/\/\/<\/span><\/strong>\r\n thin(20) set(f)<\/span><\/strong>\r\n wbsrun using script.txt<\/span><\/strong>\r\n type log.txt<\/span><\/strong>\r\n wbscoda using diab , clear<\/span><\/strong>\r\n save pred.dta, replace<\/span><\/strong><\/pre>\n
mcmccheck , dfile(stack_diabetes.dta) d(r) pfile(pred.dta) \/\/\/<\/span><\/strong>\r\n p(f_) plot(summary) save(diabetes_res.dta,replace)<\/span><\/strong><\/pre>\n
\"NM2summary1\"<\/a><\/p>\n
This plot confirms that we have a reasonable model (perhaps the variability is slightly larger when fewer cases of diabetes are observed, but if so the effect is not strong).<\/p>\n
To plot the coefficients and their credible intervals as we did previously for the confidence intervals after the glm analysis,
\nuse diabetes_fit.dta, clear<\/span><\/strong>
\nmcmcstats d_2 d_3 d_4 d_5 d_6 , save(est.dta,replace)<\/span><\/strong>
\nuse est.dta, clear<\/span><\/strong>
\ngen id = _n<\/span><\/strong>
\ntwoway (scatter mean id ) (rcap ub lb id ) , leg(off) \/\/\/<\/span><\/strong>
\nxlabel(1 “Placebo” 2 “b-blocker” 3 “CCB” 4 “ACE” 5 “ARB”) \/\/\/<\/span><\/strong>
\nyline(0) text(0.05 4 “Diuretic”) xtitle(Treatment)<\/span><\/strong><\/p>\n
\"NM2CrI\"<\/a><\/p>\n
In the NICE report the authors\u00a0go on to fit a random effects model that allows the treatment effects to vary randomly by study. Here is a simplified model that is almost identical to the one used in the report.<\/p>\n
model{ <\/span><\/strong>\r\n\u00a0\u00a0 for(i in 1:8){ <\/span><\/strong>\r\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 r[i] ~ dbin(p[i],n[i]) <\/span><\/strong>\r\n\u00a0\u00a0       e[i] <- 0<\/span><\/strong>\r\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 cloglog(p[i]) <- mu[id[i]]+d[t[i]]<\/span><\/strong>\r\n\u00a0\u00a0 }<\/span><\/strong>\r\n\u00a0\u00a0 for(i in 9:48){ <\/span><\/strong>\r\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 r[i] ~ dbin(p[i],n[i]) <\/span><\/strong>\r\n\u00a0\u00a0       e[i] ~ dnorm(0,tau)<\/span><\/strong>\r\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 cloglog(p[i]) <- mu[id[i]]+d[t[i]]+e[i]<\/span><\/strong>\r\n\u00a0\u00a0 }<\/span><\/strong>\r\n\u00a0\u00a0 d[1]<-0 <\/span><\/strong>\r\n\u00a0\u00a0 for(k in 2:6){ d[k] ~ dnorm(0,0.0001) } <\/span><\/strong>\r\n\u00a0\u00a0 for(k in 1:22){ mu[k] ~ dnorm(0,0.0001) } <\/span><\/strong>\r\n\u00a0\u00a0 tau ~ dgamma(0.001,0.001)<\/span><\/strong>\r\n}<\/span><\/strong><\/pre>\n
A random effect e[] is added that allows the treatment effect d[] to vary randomly from study to study. The data are sorted so that the first 8 lines\u00a0refer to trials in which\u00a0subjects were treated with diuretics, that is to say with the baseline treatment. The DIC and parameter estimates for this model are almost identical to those given in the report for their model, which gives e[] a more compex variance structure.<\/p>\n
The use of a random effects model reduces the DIC from 426 to 412, which\u00a0by most guidelines would be considered as showing a better fit. In the report they certainly use this reduction to justify the random effects model. However if we look at the residuals we get a slightly different impression.<\/p>\n
\"NM2summary2\"<\/a><\/p>\n
The residual are not as dispersed as we would expect under a binomial model and it looks to me like an example of over-fitting.<\/p>\n
The DIC is widely used and almost as widely criticised.\u00a0Its use\u00a0has been questioned\u00a0when we have a model with a large number of random effects and a small sample size. In this example there are 48 observations and the model with random effects has pD=37 suggesting the equivalent of 37 parameters. In this case I\u00a0do not believe the DIC and I would use the earlier fixed effects model.<\/p>\n
The WinBUGS code for the models recommended in the NICE report can be found at http:\/\/www.nicedsu.org.uk\/TSD2%20General%20meta%20analysis%20corrected%2015April2014.pdf<\/span><\/a>\u00a0(the diabetes example\u00a0starts on page 66). Those models can be fitted from within Stata using commands very similar to those given in the posting.<\/p>\n","protected":false},"excerpt":{"rendered":"
In my previous posting (https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2014\/07\/18\/network-meta-analysis-in-stata\/) I introduced an example of Network Meta-analysis taken from a NICE report that estimated the incidence of diabetes amongst people being treated for high blood pressure by one of six classes of drug. In that posting I analysed the data using Stata\u2019s glm command and now I should like to […]<\/p>\n","protected":false},"author":134,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[25,24,4,5],"class_list":["post-198","post","type-post","status-publish","format-standard","hentry","category-uncategorized","tag-dic","tag-network-meta-analysis","tag-stata","tag-winbugs"],"_links":{"self":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/198","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/users\/134"}],"replies":[{"embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/comments?post=198"}],"version-history":[{"count":11,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/198\/revisions"}],"predecessor-version":[{"id":225,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/198\/revisions\/225"}],"wp:attachment":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/media?parent=198"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/categories?post=198"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/tags?post=198"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}