{"id":723,"date":"2015-07-03T08:11:07","date_gmt":"2015-07-03T08:11:07","guid":{"rendered":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/?p=723"},"modified":"2025-02-26T13:21:37","modified_gmt":"2025-02-26T13:21:37","slug":"stan-vs-openbugs-controlled-from-stata","status":"publish","type":"post","link":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2015\/07\/03\/stan-vs-openbugs-controlled-from-stata\/","title":{"rendered":"Stan vs OpenBUGS (controlled from Stata)"},"content":{"rendered":"

A rather long posting this week for which I apologise. I would have split this\u00a0posting over two weeks but I am about to go on my summer holiday and I wanted\u00a0to avoid a long gap between the two halves.<\/p>\n

I have spent\u00a0the last few\u00a0weeks discussing how Stan<\/strong> can be called from Stata<\/strong> and now the time has come to test the process with a realistic problem. To create a\u00a0comparison, I thought that I would take a problem that is challenging for OpenBUGS and see whether Stan offers any advantages. In the spirit of this blog, the question ought to be,
\n\u00a0\u00a0\u00a0 is it as easy to run Stan from Stata, as it is to run OpenBUGS from Stata?<\/strong>
\nbut the answer to that question is clearly, yes, so inevitably we will get drawn into the issue,
\nis Stan\u00a0as good as\u00a0OpenBUGS?<\/strong><\/p>\n

For my test problem I have chosen a network meta-analysis as described in,<\/p>\n

White I.R., Barrett J.K., Jackson D, Higgins J.P.T.
\nConsistency and inconsistency in network meta-analysis: model estimation using multivariate meta-regression<\/strong>
\nResearch Synthesis Methods<\/em> 2012;3:111-125<\/p>\n

In that paper the authors\u00a0present the data for a worked example and give the WinBUGS code for fitting their model; so much of my work is done for me.\u00a0What is more,\u00a0their analysis was performed with three chains, each with a burn-in of 30,000 and a run length of 150,000 so the model is clearly challenging enough to make a good test.<\/p>\n

I\u2019ve blogged about network meta-analysis before here<\/a> and here<\/a>. I have to admit to having some doubts about its usefulness\u00a0but I was impressed that\u00a0the paper by White et al. addresses one of\u00a0my main concerns, how\u00a0do we detect\u00a0when\u00a0a network meta-analysis is inconsistent?<\/p>\n

The example used in the paper\u00a0is a meta-analysis of\u00a028 trials of drug treatments given to people who had had a heart attack. The outcome was the number of deaths in the first month after treatment. 8 different drugs were compared in the trials but the designs varied; for instance,\u00a0some compared streptokinase with alteplase, some compared streptokinase with urokinase, yet others compared alteplase with urokinase and so on.<\/p>\n

Well, if streptokinase is better than alteplase by one unit on some outcome scale and alteplase is better than urokinase by one unit on the same scale, then streptokinase\u00a0ought to\u00a0be better than urokinase by two units. If this is so, then the data are said to be consistent. It is possible however that the streptokinase-urokinase trials were different in some systematic way, perhaps they were conducted more recently or they all come from the USA.\u00a0If they do differ systematically then\u00a0the directly measured effect\u00a0may not be consistent with the indirectly implied effect. White et al. refer to the set of treatments compared in a study as the \u2018design\u2019 and handle inconsistency by treating it as an interaction between the treatment effect and the design.<\/p>\n

There is much more to the model that they propose but I do not think that it is necessary for us to understand all of the details as our primary interest is in seeing if Stan can handle it.<\/p>\n

Let\u2019s first fit the model with design by treatment interaction using OpenBUGS and the code provided in White et al. I have removed a few lines from White et al.’s program that calculate derived quantities that are not essential for the estimation. The data, taken from Table 1 of the paper, are in a file called thrombolytic.dta<\/strong> that is available from my webpage along with all of the\u00a0Stata code .<\/p>\n

Here is the full Stata program.<\/p>\n

use “thrombolytic.dta”, clear<\/span>
\n * Number of treatments in each study<\/span>
\n sort study<\/span>
\n by study: egen nt = count(trt>0)<\/span>
\n * offsets – row numbers for start of each study or design<\/span>
\n sort design study trt<\/span>
\n gen os = _n*(study != study[_n-1])<\/span>
\n gen os_design = _n*(design != design[_n-1])<\/span>
\n by design: egen b = min(trt)<\/span>
\n replace os = 59 in 58<\/span>
\n * Write data to file<\/span>
\n wbslist (vect study design trt b r n, format(%5.0f) name(study d t b r n) ) \/\/\/<\/span>
\n\u00a0\u00a0 (vect os if os != 0 , format(%5.0f) name(offset) ) \/\/\/<\/span>
\n\u00a0\u00a0 (vect os_design nt if os_design != 0 , format(%5.0f) name(offset.design num.ests) ) \/\/\/<\/span>
\n\u00a0\u00a0 (zero=c(7{0})) (A=58, S=28, T=8, D=13) using data.txt, replace<\/span>
\n * Prepare three sets of starting values for 3 chains<\/span>
\n matrix RE = J(28,8,0)<\/span>
\n forvalues i=1\/28 {<\/span>
\n\u00a0\u00a0 matrix RE[`i’,1] = .<\/span>
\n }<\/span>
\n wbslist (matrix RE , format(%5.0f) )\u00a0<\/span>\u00a0( mu=c(28{-1}) ) <\/span>( tau=0.2) \/\/\/<\/span>
\n\u00a0\u00a0\u00a0 using init1.txt, replace<\/span>
\n wbslist (matrix RE , format(%5.0f) ) <\/span>( mu=c(28{-2}) ) <\/span>( tau=1) \/\/\/<\/span>
\n\u00a0\u00a0\u00a0 using init2.txt, replace<\/span>
\n wbslist (matrix RE , format(%5.0f) ) <\/span>( mu=c(28{-1.5}) ) <\/span>( tau=0.5) \/\/\/<\/span>
\n\u00a0\u00a0\u00a0 using init3.txt, replace<\/span>
\n * Write the model<\/span>
\n wbsmodel thrombo_bayes.do model.txt<\/span>
\n \/*<\/span>
\n model {<\/span>
\n\u00a0\u00a0 for(i in 1:S) {<\/span>
\n\u00a0\u00a0\u00a0\u00a0 eff.study[i, b[offset[i]], b[offset[i]]] <- 0<\/span>
\n\u00a0\u00a0\u00a0\u00a0 for(k in (offset[i] + 1):(offset[i + 1]-1)) {<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 eff.study[i,t[k],b[k]] <- eff.des[d[k],t[k]] + RE[i,t[k]] – RE[i,b[k]]<\/span>
\n\u00a0\u00a0\u00a0\u00a0 }<\/span>
\n\u00a0\u00a0 }<\/span>
\n # Random effects for heterogeneity<\/span>
\n for(i in 1:S) {<\/span>
\n\u00a0\u00a0 RE[i,1] <- 0<\/span>
\n RE[i,2:T] ~ dmnorm(zero[], Prec[,])<\/span>
\n }<\/span>
\n # Prec is the inverse of the structured heterogeneity matrix<\/span>
\n for(i in 1:(T-1)) {<\/span>
\n\u00a0\u00a0 for(j in 1:(T-1)){<\/span>
\n\u00a0\u00a0\u00a0\u00a0 Prec[i,j] <- 2*(equals(i,j)-1\/T)\/(tau*tau)<\/span>
\n\u00a0\u00a0 }<\/span>
\n }<\/span>
\n for(i in 1:A) {<\/span>
\n\u00a0\u00a0 logit(p[i]) <- mu[study[i]] + eff.study[study[i],t[i],b[i]]<\/span>
\n\u00a0\u00a0 r[i] ~ dbin(p[i],n[i])<\/span>
\n }<\/span>
\n # Priors<\/span>
\n for(i in 1:S) {<\/span>
\n mu[i] ~ dnorm(0,0.01)<\/span>
\n }<\/span>
\n tau ~ dunif(0,2)<\/span>
\n for(i in 1:D) {<\/span>
\n\u00a0\u00a0 for(k in (offset.design[i] + 1):(offset.design[i] + num.ests[i])) {<\/span>
\n\u00a0\u00a0\u00a0\u00a0 eff.des[i,t[k]] ~ dnorm(0,0.01)<\/span>
\n\u00a0\u00a0 }<\/span>
\n }<\/span>
\n }<\/span>
\n *\/<\/span>
\n * Prepare the script<\/span>
\n wbsscript using script.txt , replace openbugs \/\/\/<\/span>
\n\u00a0\u00a0 model(model.txt) data(data.txt) init(init1.txt+init2.txt+init3.txt) \/\/\/<\/span>
\n\u00a0\u00a0 set(tau eff.des RE mu) burn(30000) update(150000) thin(20) \/\/\/<\/span>
\n\u00a0\u00a0 log(openlog.txt) coda(thromb)<\/span>
\n * Run the analysis<\/span>
\n wbsrun using script.txt , openbugs<\/span>
\n * check the log file<\/span>
\n type openlog.txt<\/span>
\n * read results and store<\/span>
\n wbscoda using thromb , openbugs clear chains(1 2 3)<\/span>
\n save thrombmcmc123.dta, replace<\/span>
\n * calculate the interaction terms<\/span>
\n gen w1 = eff_des_3_3 – eff_des_2_3 \/\/ wAC3<\/span>
\n gen w2 = eff_des_4_4 – eff_des_1_4 \/\/ wAD4<\/span>
\n gen w3 = eff_des_9_6 – eff_des_5_6 + eff_des_1_2 \/\/ wAF9<\/span>
\n gen w4 = eff_des_10_7 – eff_des_6_7 + eff_des_1_2 \/\/ wAG10<\/span>
\n gen w5 = eff_des_12_7 – eff_des_6_7 + eff_des_2_3 \/\/ wAG12<\/span>
\n gen w6 = eff_des_7_8 – eff_des_2_8 \/\/ wAH7<\/span>
\n gen w7 = eff_des_11_8 – eff_des_2_8 + eff_des_1_2 \/\/ wAH11<\/span>
\n gen w8 = eff_des_13_8 – eff_des_2_8 + eff_des_2_3 \/\/ wAH13<\/span>
\n gen dAE = eff_des_8_5 + eff_des_1_2<\/span>
\n * various convergence checks<\/span>
\n mcmctrace w1 w3 w5 w7 , gopt( xline(7500) xline(15000) yscale(range(-4 6)) \/\/\/<\/span>
\n\u00a0\u00a0 ylabel(-4(2)6) xscale(range(0 22500)) xlabel(0(5000)20000) )<\/span>
\n mcmcstats_v2 w* , newey(200)<\/span>
\n gen w7_1 = w7 if chain == 1<\/span>
\n gen w7_2 = w7 if chain == 2<\/span>
\n gen w7_3 = w7 if chain == 3<\/span>
\n mcmcac w7_1 w7_2 w7_3 , gopt(yscale(range(0 0.5)) <\/span>ylabel(0(0.2)0.4) ) acopt(lag(50))<\/span><\/p>\n

Here we have 3 chains each with a burnin of 30,000 and a run length of 150,000 which makes 540,000 iterations. On my laptop, OpenBUGS took about 1s for 1000 updates so that the entire run took about 9 minutes. Here are the results for the interaction terms based on the combined chains. The results are very similar to those given in the paper and the only interaction that is clearly non-zero is w7<\/p>\n

----------------------------------------------------------------------------\r\nParameter\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\r\n----------------------------------------------------------------------------\r\nw1\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 22500\u00a0\u00a0 -0.165\u00a0\u00a0 0.366\u00a0\u00a0 0.0111\u00a0\u00a0 -0.103 ( -1.038,\u00a0\u00a0 0.452 )\r\nw2\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 22500\u00a0\u00a0  0.479\u00a0\u00a0 0.808\u00a0\u00a0 0.0085\u00a0\u00a0  0.470 ( -1.101,\u00a0\u00a0 2.092 )\r\nw3\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 22500\u00a0\u00a0 -0.193\u00a0\u00a0 0.504\u00a0\u00a0 0.0103\u00a0\u00a0 -0.158 ( -1.300,\u00a0\u00a0 0.775 )\r\nw4\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 22500\u00a0\u00a0  0.387\u00a0\u00a0 0.791\u00a0\u00a0 0.0091\u00a0\u00a0  0.374 ( -1.136,\u00a0\u00a0 1.992 )\r\nw5\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 22500\u00a0\u00a0  0.034\u00a0\u00a0 0.762\u00a0\u00a0 0.0107\u00a0\u00a0  0.034 ( -1.490,\u00a0\u00a0 1.545 )\r\nw6\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 22500\u00a0\u00a0 -0.049\u00a0\u00a0 0.468\u00a0\u00a0 0.0095\u00a0\u00a0 -0.050 ( -0.991,\u00a0\u00a0 0.896 )\r\nw7\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 22500\u00a0\u00a0  1.263\u00a0\u00a0 0.634\u00a0\u00a0 0.0138\u00a0\u00a0  1.247 (\u00a0\u00a00.049,\u00a0\u00a0 2.561 )\r\nw8\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 22500\u00a0\u00a0 -0.321\u00a0\u00a0 0.505\u00a0\u00a0 0.0064\u00a0\u00a0 -0.319 ( -1.330,\u00a0\u00a0 0.674 )\r\n----------------------------------------------------------------------------<\/pre>\n
To inspect the mixing I looked at some trace plots. As usual they are not\u00a0particularly informative. The vertical lines\u00a0separate the three chains.<\/p>\n
\"nw_trace1\"<\/a><\/p>\n
So next I checked the Brooks-Gelman-Rubin plots. Here is the plot for w7.<\/p>\n
\"nw_bgr1\"<\/a><\/p>\n
This plot compares the within chain variability to the variability of the pooled chains, R being the ratio of the two. Ideally R should be 1.0. Perhaps the first 1,000 thinned iterations have not quite settled, which would argue for a slightly longer burnin or better starting values.<\/p>\n
To check this I ran\u00a0a version of the Geweke test as implemented in mcmcgeweke<\/strong>. This compares w7 as calculated from the first 10% of each chain with the value from the final 50%. Not only is there\u00a0 evidence of drift, but the three estimates from the second halves of the chains are still subject to quite a lot of sampling error.<\/p>\n\n\n\n\n\n\n\n
<\/td>\nFirst 10%<\/td>\n<\/td>\nLast 50%<\/td>\n<\/td>\n<\/td>\n<\/tr>\n
<\/td>\nMean<\/td>\nSt Error<\/td>\nMean<\/td>\nSt Error<\/td>\np-value<\/td>\n<\/tr>\n
Chain 1<\/td>\n1.342<\/td>\n0.051<\/td>\n1.249<\/td>\n0.014<\/td>\n0.080<\/td>\n<\/tr>\n
Chain 2<\/td>\n1.149<\/td>\n0.028<\/td>\n1.267<\/td>\n0.014<\/td>\n0.0001<\/td>\n<\/tr>\n
Chain 3<\/td>\n1.159<\/td>\n0.037<\/td>\n1.322<\/td>\n0.022<\/td>\n0.0002<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n
The reason that such a long chain is required is the strong autocorrelation. For me the autocorrelation plot is key. Remember that the stored values have already been thinned by 20, yet when we look at the autocorrelations for the three chains\u00a0using mcmcac<\/strong> there is still noticeable correlation up to a lag of about 30 (that is, about 600 of the original iterations).<\/p>\n
\"nw_ac1\"<\/a>
\nThere are convergence issues that remain with the OpenBUGS analysis but it is time to try the same model in Stan<\/strong>. I\u2019ll start by presenting a basic Stata program that performs the analysis and then I\u2019ll discuss the changes in the code that are needed when we switch from OpenBUGS to Stan.<\/p>\n
Here is the Stata program<\/p>\n
use “thrombolytic.dta”, clear<\/span>
\n * Number of treatments in each study<\/span>
\n sort study<\/span>
\n by study: egen nt = count(trt>0)<\/span>
\n * offsets – row numbers for start of each study or design<\/span>
\n sort design study trt<\/span>
\n gen os = _n*(study != study[_n-1])<\/span>
\n gen os_design = _n*(design != design[_n-1])<\/span>
\n by design: egen b = min(trt)<\/span>
\n replace os = 59 in 58<\/span>
\n *——————————————————–<\/span>
\n * Write the data file: jags & Stan use the same format<\/span>
\n *——————————————————–<\/span>
\n wbslist_v2 (vect study design trt b r n, format(%5.0f) name(study d t b r n) ) \/\/\/<\/span>
\n (vect os if os != 0 , format(%5.0f) name(offset) ) \/\/\/<\/span>
\n (vect os_design nt if os_design != 0 , format(%5.0f) name(offset_design num_ests) ) \/\/\/<\/span>
\n (zero <- c(7{0.0})) (A <-58 _n S <- 28 _n T <- 8 _n D <- 13) \/\/\/<\/span>
\n using thromData.txt, replace jags<\/span>
\n *—————————————<\/span>
\n * Write the model file<\/span>
\n *—————————————<\/span>
\n wbsmodel thrombo_stan.do thromModel.stan<\/span>
\n \/*<\/span>
\n data {<\/span>
\n\u00a0\u00a0 int A;<\/span>\u00a0\u00a0\u00a0int S;<\/span>\u00a0\u00a0 int T;<\/span>\u00a0\u00a0 int D;<\/span>
\n\u00a0\u00a0 int study[A];<\/span>\u00a0\u00a0 int d[A];\u00a0\u00a0\u00a0 <\/span>int t[A];<\/span>\u00a0\u00a0 int b[A];<\/span>\u00a0\u00a0 int r[A];<\/span>\u00a0\u00a0 int n[A];<\/span>
\n\u00a0\u00a0 int offset[29];<\/span>\u00a0\u00a0 int offset_design[13];<\/span>\u00a0\u00a0 int num_ests[13];<\/span>
\n\u00a0\u00a0 vector[7] zero;<\/span>
\n }<\/span>
\n parameters {<\/span>
\n\u00a0\u00a0 real<lower=0.0,upper=2.0> tau;<\/span>
\n\u00a0\u00a0 real eff_des[A,A];<\/span>
\n\u00a0\u00a0 vector[7] RE[S];<\/span>
\n\u00a0\u00a0 real mu[S];<\/span>
\n }<\/span>
\n transformed parameters {<\/span>
\n\u00a0\u00a0 matrix[7,7] Prec;<\/span>
\n\u00a0\u00a0 real eff_study[S,A,A];<\/span>
\n # Prec is the inverse of the structured heterogeneity matrix<\/span>
\n\u00a0\u00a0 for(i in 1:7) {<\/span>
\n\u00a0\u00a0\u00a0\u00a0 Prec[i,i] <- 1.75\/(tau*tau);<\/span>
\n\u00a0\u00a0\u00a0\u00a0 for(j in 1:7){<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 if( i != j ) Prec[i,j] <- -0.25\/(tau*tau);<\/span>
\n\u00a0\u00a0\u00a0\u00a0 }<\/span>
\n\u00a0\u00a0 }<\/span>
\n\u00a0\u00a0 for(i in 1:S) {<\/span>
\n\u00a0\u00a0\u00a0\u00a0 eff_study[i, b[offset[i]], b[offset[i]]] <- 0;<\/span>
\n\u00a0\u00a0\u00a0\u00a0 for(k in (offset[i] + 1):(offset[i + 1]-1)) {<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 if( t[k] == 1 ) {<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 if( b[k] == 1 ) eff_study[i,t[k],b[k]] <- eff_des[d[k],t[k]];<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 else eff_study[i,t[k],b[k]] <- eff_des[d[k],t[k]] – RE[i][b[k]-1];<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 }<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 else {<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 if( b[k] == 1 ) eff_study[i,t[k],b[k]] <- eff_des[d[k],t[k]] + RE[i][t[k]-1];<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 else eff_study[i,t[k],b[k]] <- eff_des[d[k],t[k]] + RE[i][t[k]-1] – RE[i][b[k]-1];<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0 }<\/span>
\n\u00a0\u00a0 \u00a0}<\/span>
\n\u00a0 }<\/span>
\n }<\/span>
\n model {<\/span>
\n # Random effects for heterogeneity<\/span>
\n\u00a0\u00a0 for(i in 1:S) {\u00a0 <\/span>
\n\u00a0\u00a0\u00a0\u00a0 RE[i] ~ multi_normal_prec(zero, Prec);<\/span>
\n\u00a0\u00a0 }<\/span>
\n\u00a0\u00a0 for(i in 1:A) {<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0 r[i] ~ binomial_logit(n[i],mu[study[i]] + eff_study[study[i],t[i],b[i]]);<\/span>
\n\u00a0\u00a0\u00a0}<\/span>
\n # Priors<\/span>
\n\u00a0\u00a0 for(i in 1:S) {<\/span>
\n\u00a0\u00a0\u00a0\u00a0 mu[i] ~ normal(0.0,10.0);<\/span>
\n\u00a0\u00a0 }<\/span>
\n\u00a0\u00a0 tau ~ uniform(0.0,2.0);<\/span>
\n\u00a0\u00a0 for(i in 1:D) {<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0 for(k in (offset_design[i] + 1):(offset_design[i] + num_ests[i])) {<\/span>
\n\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 eff_des[i,t[k]] ~ normal(0.0,10.0);<\/span>
\n\u00a0\u00a0 \u00a0\u00a0 }<\/span>
\n\u00a0\u00a0 }<\/span>
\n }<\/span>
\n *\/<\/span>
\n *—————————————<\/span>
\n * Prepare a Stan batch file<\/span>
\n *—————————————<\/span>
\n wbsscript_v3 using throm.bat , replace \/\/\/<\/span>
\n model(thromModel.stan) data(thromData.txt) \/\/\/<\/span>
\n updates(3000) burnin(500) coda(outthromprec) stan<\/span>
\n *—————————————<\/span>
\n * Run the batch file<\/span>
\n *—————————————<\/span>
\n wbsrun_v3 using throm.bat , stan<\/span>
\n type log.txt<\/span>
\n *—————————————<\/span>
\n * read the results into Stata<\/span>
\n *—————————————<\/span>
\n wbscoda_v3 using outthromprec.csv , stan clear<\/span><\/p>\n
The Stata code for preparing the data file is almost unchanged. Stan uses the same data file format as JAGS, so we merely set the jags option to wbslist_v2<\/strong>. Note though that when we do not use a keyword, wbslist_v2 writes our text\u00a0directly to the file and in the JAGS\/Stan format each item must be placed on a newline, so _n replaces the comma of OpenBUGS as our separator.<\/p>\n
We do not need initial values as Stan is designed to set them internally.<\/p>\n
After specifying the model we create a script file with wbsscript_v3<\/strong> much as\u00a0with OpenBUGS\u00a0except that we set the stan option. Remember that the file that is created is actually a batch file and not a script file because Stan works by compiling the code in the model file and turning it into an executable program that\u00a0is run via the Windows operating system.<\/p>\n
The results are written to a comma delimited file that could be read by Excel. However, the file contains a header of summary information about the run so we need to use wbscoda_v3<\/strong> to extract the simulations and to read them into Stata.<\/p>\n
This leaves us with the Stan\u00a0model file. Of course there are many ways of writing any piece of computer code and my experience of Stan is rather limited so my code may not be the most efficient, but at least it works.<\/p>\n
The first block defines the data and sets up space for the information that will be read from the data file. Remember that Stan does not see the data until after the program has been compiled so it needs to\u00a0be told\u00a0the sizes of the variables. The next block defines the parameters, i.e. everything for which you give a prior. Stan will choose random starting values for these parameters. Next come the transformed parameters. These are\u00a0quantities that can be calculated from the basic parameters. For instance, the precision matrix, Prec, can be calculated from tau. Stan will\u00a0obtain the starting values for the transformed parameters by deriving them from the starting values for the basic parameters.<\/p>\n
When the elements of Prec are calculated from tau, the OpenBUGS code used the equals() function, while in Stan I code the test explicitly using \u201c== \u201c<\/p>\n
You will notice some other differences in the code. In OpenBUGS we had a matrix of parameters RE in which the first column was fixed at zero and the other columns were estimated.\u00a0Use of\u00a0fixed and estimated values in the same structure is not\u00a0allowed in Stan so\u00a0I formed a parameter matrix RE with one less column and treated references to column 1 as a special case in the subsequent code. In fact I have not defined RE to be a matrix at all, but rather I have created a set of vectors corresponding to the different columns. This is somewhat more natural since each vector(column) is drawn separately\u00a0from a multivariate normal distribution.<\/p>\n
Finally, notice that the distributions have different names in Stan and there are many variations on each distribution. So a multivariate normal, multi_normal(), is used if you know the mean and covariance matrix and multi_normal_prec() is used if you know the mean and precision matrix. binomial() is used when we know n and p, but binomial_logit() is used when we know n and logit(p).<\/p>\n
The Stan syntax takes a little getting used to if you are a long-standing BUGS user, but in truth it is closer to a standard programming language and much more flexible than WinBUGS or OpenBUGS.<\/p>\n
Anyway the code runs, but don\u2019t\u00a0start it if you are in a hurry. It is incredibly slow! 3000 updates and a 500 burnin took about 1.5 hours.<\/p>\n
Although\u00a0Stan is very slow, the estimates that it produces are of very high quality. The 3,000 updates have a\u00a0similar sem to that given by 450,000 OpenBUGS updates.<\/p>\n
—————————————————————————-
\nParameter\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 n\u00a0\u00a0\u00a0\u00a0 mean\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 sd\u00a0\u00a0 \u00a0sem\u00a0\u00a0 median\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 95% CrI
\n—————————————————————————-
\nw7\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 3000\u00a0\u00a0 1.249\u00a0\u00a0 0.646\u00a0\u00a0 0.0141\u00a0\u00a0 1.239 ( -0.048,\u00a0\u00a0 2.515 )
\n—————————————————————————-<\/p>\n
What is more there is no indication of early drift<\/p>\n
\"nw_trace2\"<\/a><\/p>\n
The visual impression of the trace plot is confirmed by the Geweke analysis, which shows a\u00a0mean (se) of 1.269 (0.0245) in the first 10% and 1.234(0.0172) in the final 50% with a p-value 0.242.<\/p>\n
Most importantly,\u00a0the autocorrelation is negligible, so we have the equivalent of a random sample.<\/p>\n
\"nw_ac2\"<\/p>\n
 <\/p>\n
I hit one other problem. Stan wrote all of the parameters to the output file and\u00a0there are so many that neither Excel nor Stata could cope. So to obtain the estimate of w7 I had to open part of the output\u00a0file in Excel and then manually copy and paste the relevant columns into Stata.<\/p>\n
This particular model creates estimates for every treatment comparison in every study. Hence the huge output file. I suspect that it would be possible to reparameterize the model with far fewer active parameters, which might greatly improve the performance of Stan. One of the difficulties that Stan had in this example\u00a0is that the code that I wrote is based on a structure and parameterization that was intended for OpenBUGS.<\/p>\n
For me, there is no clear winner even though Stan was\u00a0much slower. I have been using WinBUGS for years and I know more or less how to write efficient WinBUGS code.\u00a0Perhaps I need to put the same investment into Stan before I\u00a0will be\u00a0in a position to judge how good it is.<\/p>\n
If you want to download any of the programs referred to in this posting they are available from my webpage at https:\/\/www2.le.ac.uk\/Members\/trj<\/a><\/p>\n","protected":false},"excerpt":{"rendered":"
A rather long posting this week for which I apologise. I would have split this\u00a0posting over two weeks but I am about to go on my summer holiday and I wanted\u00a0to avoid a long gap between the two halves. I have spent\u00a0the last few\u00a0weeks discussing how Stan can be called from Stata and now the […]<\/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":[69,24,48,56],"class_list":["post-723","post","type-post","status-publish","format-standard","hentry","category-uncategorized","tag-autocorrelation","tag-network-meta-analysis","tag-openbugs","tag-stan"],"_links":{"self":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/723","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=723"}],"version-history":[{"count":20,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/723\/revisions"}],"predecessor-version":[{"id":751,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/723\/revisions\/751"}],"wp:attachment":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/media?parent=723"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/categories?post=723"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/tags?post=723"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}