{"id":158,"date":"2014-06-06T08:59:19","date_gmt":"2014-06-06T08:59:19","guid":{"rendered":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/?p=158"},"modified":"2025-02-26T13:21:38","modified_gmt":"2025-02-26T13:21:38","slug":"model-free-genetic-meta-analysis","status":"publish","type":"post","link":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/2014\/06\/06\/model-free-genetic-meta-analysis\/","title":{"rendered":"Model-free Genetic Meta-analysis"},"content":{"rendered":"

Some years ago some colleagues and I wrote a paper about the Bayesian meta-analysis of genetic association studies (Minelli et al Statistics in Medicine<\/em> 2005;24:3845-3861) and periodically since then people have contacted me asking for the WinBUGS code that we used in that paper and those requests have prompted me to reformulate\u00a0the code\u00a0in terms of the functions for calling WinBUGS\u00a0for Stata as described in the book \u2018Bayesian Analysis with Stata<\/em>\u2019.<\/p>\n

Here is the background to the paper although of course you can get more details from the original.<\/p>\n

Many genetic studies relate a binary outcome (disease yes\/no) to a specific genetic variant. Since\u00a0our autosomal\u00a0chromosomes occur\u00a0in pairs, when a particular genetic variant can take one of two forms, say little a and big A, then that individual will we either aa, aA or AA (there is no way to distinguish aA from Aa). In effect then we have a 2×3 table of numbers of people.<\/p>\n\n\n\n\n\n
<\/td>\naa<\/td>\naA<\/td>\nAA<\/td>\n<\/tr>\n
Disease = Yes<\/td>\n<\/td>\n<\/td>\n<\/td>\n<\/tr>\n
Disease = No<\/td>\n<\/td>\n<\/td>\n<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

We model these data in terms of the odds ratio of disease, but there are two odds ratios one for ORaA<\/sub> = aA vs aa and one for ORAA<\/sub> = AA vs aa. Very often researchers assume that ORAA<\/sub> = ORaA<\/sub>*ORaA<\/sub> and for many diseases and variants this works well because it is like assuming that one\u2019s risk increases equally with each extra copy of the A variant.<\/p>\n

However, the assumption that ORAA<\/sub> = ORaA<\/sub>*ORaA<\/sub> does not always apply and in the paper we\u00a0suggested that we could use a meta-analysis to estimate both the\u00a0form of the relationship between ORs\u00a0and their sizes. We expressed the relationship in terms of a parameter<\/p>\n

\u03bb = logORaA<\/sub>\/logORAA<\/sub><\/p>\n

What is more if we have an estimate of \u03bb for a particular disease and variant, even if that estimate is subject to uncertainty, we can then use the information about ORaA<\/sub> to improve our estimate of the ORAA<\/sub>. Since a known or assumed relationship between ORaA<\/sub> and ORAA<\/sub> is called the mode of inheritance or\u00a0sometimes the\u00a0genetic model, we referred to our analysis as \u2018genetic model free\u2019.<\/p>\n

We fitted this model using WinBUGS and investigated the impact of using different prior distributions using published and simulated data.\u00a0Here I will show how we can automate such an investigation.<\/p>\n

Let us start with some real data. Kato et al (Journal of Hypertension<\/em> 1999;17:757-763) reported a meta-analysis of studies of a variant in the angiotensinogen gene and hypertension. Here is a summary of their data.<\/p>\n\n\n\n\n\n\n\n\n\n\n\n
Study<\/strong><\/td>\n<\/td>\n\u00a0Cases<\/strong><\/strong><\/td>\n<\/td>\n\u00a0<\/strong><\/td>\nControls<\/strong><\/td>\n\u00a0<\/strong><\/td>\n<\/tr>\n
<\/td>\nMM<\/strong><\/td>\nTM<\/strong><\/td>\nTT<\/strong><\/td>\nMM<\/strong><\/td>\nTM<\/strong><\/td>\nTT<\/strong><\/td>\n<\/tr>\n
Hata, 1994<\/em><\/td>\n2<\/td>\n20<\/td>\n83<\/td>\n3<\/td>\n34<\/td>\n44<\/td>\n<\/tr>\n
Iwai, 1994<\/em><\/td>\n3<\/td>\n17<\/td>\n62<\/td>\n4<\/td>\n30<\/td>\n49<\/td>\n<\/tr>\n
Nishiuma, 1995<\/em><\/td>\n3<\/td>\n30<\/td>\n31<\/td>\n17<\/td>\n84<\/td>\n48<\/td>\n<\/tr>\n
Morise, 1995<\/em><\/td>\n5<\/td>\n23<\/td>\n52<\/td>\n5<\/td>\n32<\/td>\n63<\/td>\n<\/tr>\n
Sato, 1997<\/em><\/td>\n8<\/td>\n39<\/td>\n133<\/td>\n12<\/td>\n62<\/td>\n119<\/td>\n<\/tr>\n
Kamitani, 1994<\/em><\/td>\n6<\/td>\n34<\/td>\n68<\/td>\n9<\/td>\n48<\/td>\n47<\/td>\n<\/tr>\n
Kato, 1999<\/em><\/td>\n20<\/td>\n214<\/td>\n483<\/td>\n18<\/td>\n134<\/td>\n363<\/td>\n<\/tr>\n<\/tbody>\n<\/table>\n

 <\/p>\n

So here the two variants are labelled as M and T and in the first study 83 people with hypertension had the TT genotype. The key question is whether the T variant is more or less common in cases than in controls, with the genetic model being of secondary interest.<\/p>\n

In the first study the researchers recruited 81 controls and found that they divided in the ratio 3:34:44. This we can model these data by a multinomial distribution with n=81 and probabilities p01<\/sub>,p02<\/sub>,p03<\/sub>. We choose to parameterize \u03b21=1, \u03b22=p02<\/sub>\/p01<\/sub> and \u03b23= p03<\/sub>\/p01<\/sub>, then<\/p>\n

p01=1\/(1+ \u03b22+ \u03b23),\u00a0\u00a0\u00a0 p02= \u03b22\/(1+ \u03b22+ \u03b23)\u00a0\u00a0\u00a0 and\u00a0\u00a0\u00a0\u00a0 p03= \u03b23\/(1+ \u03b22+ \u03b23).<\/p>\n

Similarly the 105 cases were in the ratio 2:20:83 and\u00a0if these categories have probabilities q01<\/sub>,q02<\/sub>,q03<\/sub>, we can parameterize them in terms of \u03b11=1, \u03b12=q02<\/sub>\/q01<\/sub> and \u03b13= q03<\/sub>\/q01<\/sub>. But \u03b12\/\u03b22 is ORTM<\/sub> and \u03b13\/\u03b23 is ORTT<\/sub> so\u00a0we\u00a0define \u03b42=logORTM<\/sub> and \u03b43=logORTT<\/sub> and\u00a0the formulae for the probabilities in cases can be written,<\/p>\n

q01=1\/(1+ \u03b22exp(\u03b42)+ \u03b23 exp(\u03b43))<\/p>\n

q02= \u03b22exp(\u03b42)\/(1+ \u03b22exp(\u03b42)+ \u03b23 exp(\u03b43))<\/p>\n

q03= \u03b23 exp(\u03b43)\/(1+ \u03b22exp(\u03b42)+ \u03b23 exp(\u03b43))<\/p>\n

Now the problem is parameterized in terms of \u03b22, \u03b23, \u03b42, \u03b43 and the genetic model as described by \u03bb = \u03b42\/\u03b43.<\/p>\n

In a meta-analysis we would expect that \u03b22, \u03b23 could well be different for each study because genetic variants are sometimes more common in one population than another. Similarly \u03b43 might vary with population because frequently genes interact with aspects of the environment that differ between studies, however we might be willing to assume that the mode of inheritance, \u03bb, is the same for all studies.<\/p>\n

We treated the \u03b2\u2019s as fixed effects, 2 for each study. \u03b43 as a random effect N(\u03b8,\u03c42<\/sup>) and \u03bb as a constant. A basic WinBUGS program for fitting this model is,<\/p>\n

model {<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 for( i in 1:7 ) {<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 p[i,1] <- 1\/(1+b[i,1]+b[i,2])<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 p[i,2] <- b[i,1]\/(1+b[i,1]+b[i,2])<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 p[i,3] <- b[i,2]\/(1+b[i,1]+b[i,2])<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 q[i,1] <- 1\/(1+b[i,1]*exp(lambda*d[i])+b[i,2]*exp(d[i]))<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 q[i,2] <- b[i,1]*exp(lambda*d[i])\/ <\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (1+b[i,1]*exp(lambda*d[i])+b[i,2]*exp(d[i]))<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 q[i,3] <- b[i,2]*exp(d[i])\/ <\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (1+b[i,1]*exp(lambda*d[i])+b[i,2]*exp(d[i]))<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 d[i] ~ dnorm(theta,tau2)<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 ncont[i,1:3] ~ dmulti(p[i,1:3],tcont[i])<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 ncase[i,1:3] ~ dmulti(q[i,1:3],tcase[i])<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 b[i,1] ~ dnorm(0,0.0001)<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 b[i,2] ~ dnorm(0,0.0001)<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 }<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 lambda ~ dbeta(1,1)<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 theta ~ dnorm(0,0.0001)<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 tau2 ~ dgamma(0.001,0.001)<\/strong><\/span><\/p>\n

}<\/strong><\/span><\/p>\n

WinBUGS cannot fit this model and crashes but OpenBUGS runs it without a problem. Working within Stata we read the data from a file call kato.dta which has the same structure as the table shown above and variable names case1, case2 etc. Here is the entire code for fitting the model and viewing the results.<\/p>\n

use kato.dta, clear<\/strong><\/span><\/p>\n

gen tcase = case1+case2+case3<\/strong><\/span><\/p>\n

gen tcont = control3+control2+control1<\/strong><\/span><\/p>\n

wbslist (struct case1 case2 case3 , name(ncase) format(%4.0f) ) \/\/\/<\/strong><\/span><\/p>\n

\u00a0\u00a0 (struct control1 control2 control3 , name(ncont) format(%4.0f) ) \/\/\/<\/strong><\/span><\/p>\n

\u00a0\u00a0 (vect tcase tcont) using data.txt, replace<\/strong><\/span><\/p>\n

matrix b=J(7,2,0.5)<\/strong><\/span><\/p>\n

wbslist (lambda=0.5,theta=0,tau2=1,d=c(7{0})) (matrix b , format(%4.1f)) \/\/\/<\/strong><\/span><\/p>\n

\u00a0\u00a0 using init.txt , replace<\/strong><\/span><\/p>\n

tomodel modelfree.do model.txt<\/strong><\/span><\/p>\n

\/*<\/strong><\/span><\/p>\n

model {<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 for( i in 1:7 ) {<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 p[i,1] <- 1\/(1+b[i,1]+b[i,2])<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 p[i,2] <- b[i,1]\/(1+b[i,1]+b[i,2])<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 p[i,3] <- b[i,2]\/(1+b[i,1]+b[i,2])<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 q[i,1] <- 1\/(1+b[i,1]*exp(lambda*d[i])+b[i,2]*exp(d[i]))<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 q[i,2] <- b[i,1]*exp(lambda*d[i])\/ <\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0(1+b[i,1]*exp(lambda*d[i])+b[i,2]*exp(d[i]))<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 q[i,3] <- b[i,2]*exp(d[i])\/ <\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (1+b[i,1]*exp(lambda*d[i])+b[i,2]*exp(d[i]))<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 d[i] ~ dnorm(theta,tau2)<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 ncont[i,1:3] ~ dmulti(p[i,1:3],tcont[i])<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 ncase[i,1:3] ~ dmulti(q[i,1:3],tcase[i])<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 b[i,1] ~ dnorm(0,0.0001)<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 b[i,2] ~ dnorm(0,0.0001)<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 }<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 lambda ~ dbeta(1,1)<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 theta ~ dnorm(0,0.0001)<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 tau2 ~ dgamma(0.001,0.001)<\/strong><\/span><\/p>\n

}<\/strong><\/span><\/p>\n

*\/<\/strong><\/span><\/p>\n

wbsscript using script.txt, replace \/\/\/<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 model(model.txt) data(data.txt) init(init.txt) log(log.txt) \/\/\/<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 burnin(1000) update(10000) set(lambda theta tau2 b) \/\/\/<\/strong><\/span><\/p>\n

\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 coda(mf) openbugs<\/strong><\/span><\/p>\n

wbsrun using script.txt \u00a0\u00a0 , openbugs<\/strong><\/span><\/p>\n

wbscoda using mf , clear openbugs<\/strong><\/span><\/p>\n

gen OR3 = exp(theta)<\/strong><\/span><\/p>\n

gen OR2 = exp(lambda*theta)<\/strong><\/span><\/p>\n

gen sd = 1\/sqrt(tau2)<\/strong><\/span><\/p>\n

mcmcstats *<\/strong><\/span><\/p>\n

and here are the results<\/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
 b_1_1\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 18.437 11.803\u00a0\u00a0 0.7491\u00a0\u00a0 14.840 (\u00a0\u00a0 5.743, 51.549 )<\/strong><\/pre>\n
 b_1_2\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 28.600\u00a0\u00a0 18.929\u00a0\u00a0 1.7735\u00a0\u00a0 22.805 (\u00a0\u00a0 8.778, 83.620 )<\/strong><\/pre>\n
 b_2_1\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 11.348\u00a0\u00a0 5.987\u00a0\u00a0 0.2278\u00a0\u00a0 9.906 (\u00a0\u00a0 4.217, 27.469 )<\/strong><\/pre>\n
 b_2_2\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 20.655\u00a0\u00a0 11.225\u00a0\u00a0 0.6528\u00a0\u00a0 18.350 (\u00a0\u00a0 7.726, 50.830 )<\/strong><\/pre>\n
 b_3_1\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 6.615\u00a0\u00a0 1.868\u00a0\u00a0 0.0263\u00a0\u00a0 6.284 (\u00a0\u00a0 3.882, 11.079 )<\/strong><\/pre>\n
 b_3_2\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 3.879\u00a0\u00a0 1.185\u00a0\u00a0 0.0181\u00a0\u00a0 3.678 (\u00a0\u00a0 2.176,\u00a0\u00a0 6.757 )<\/strong><\/pre>\n
 b_4_1\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 7.459\u00a0\u00a0 3.066\u00a0\u00a0 0.0735\u00a0\u00a0 6.801 (\u00a0\u00a0 3.438, 15.370 )<\/strong><\/pre>\n
 b_4_2\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 14.312\u00a0\u00a0 6.138\u00a0\u00a0 0.1968\u00a0\u00a0 12.965 (\u00a0\u00a0 6.585, 30.110 )<\/strong><\/pre>\n
 b_5_1\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 5.798\u00a0\u00a0 1.577\u00a0\u00a0 0.0250\u00a0\u00a0 5.535 (\u00a0\u00a0 3.404,\u00a0\u00a0 9.478 )<\/strong><\/pre>\n
 b_5_2\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 11.550\u00a0\u00a0 3.170\u00a0\u00a0 0.0560\u00a0\u00a0 11.070 (\u00a0\u00a0 6.803, 19.010 )<\/strong><\/pre>\n
 b_6_1\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 6.534\u00a0\u00a0 2.117\u00a0\u00a0 0.0320\u00a0\u00a0 6.150 (\u00a0\u00a0 3.508, 11.730 )<\/strong><\/pre>\n
 b_6_2\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 6.929\u00a0\u00a0 2.412\u00a0\u00a0 0.0413\u00a0\u00a0 6.516 (\u00a0\u00a0 3.471, 12.870 )<\/strong><\/pre>\n
 b_7_1\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 9.813\u00a0\u00a0 1.793\u00a0\u00a0 0.0383\u00a0\u00a0 9.597 (\u00a0\u00a0 6.883, 13.910 )<\/strong><\/pre>\n
 b_7_2\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 24.927\u00a0\u00a0 5.075\u00a0\u00a0 0.1376\u00a0\u00a0 24.260 ( 16.790, 36.808 )<\/strong><\/pre>\n
 lambda\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 0.179\u00a0\u00a0 0.130\u00a0\u00a0 0.0025\u00a0\u00a0 0.155 (\u00a0\u00a0 0.007,\u00a0\u00a0 0.474 )<\/strong><\/pre>\n
 tau2\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 7.486\u00a0\u00a0 8.726\u00a0\u00a0 0.1713\u00a0\u00a0 5.117 (\u00a0\u00a0 0.887, 29.250 )<\/strong><\/pre>\n
 theta\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 0.541\u00a0\u00a0 0.252\u00a0\u00a0 0.0035\u00a0\u00a0 0.527 (\u00a0\u00a0 0.076,\u00a0\u00a0 1.083 )<\/strong><\/pre>\n
 OR3\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 1.775\u00a0\u00a0 0.502\u00a0\u00a0 0.0069\u00a0\u00a0 1.694 (\u00a0\u00a0 1.079,\u00a0\u00a0 2.954 )<\/strong><\/pre>\n
 OR2\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 1.117\u00a0\u00a0 0.135\u00a0\u00a0 0.0024\u00a0\u00a0 1.075 (\u00a0\u00a0 1.001,\u00a0\u00a0 1.465 )<\/strong><\/pre>\n
 sd\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 10000\u00a0\u00a0 0.488\u00a0\u00a0 0.235\u00a0\u00a0 0.0042\u00a0\u00a0 0.442 (\u00a0\u00a0 0.185,\u00a0\u00a0 1.062 )<\/strong><\/pre>\n
--------------------------------------------------------------------------<\/strong><\/pre>\n
I will show a couple of convergence checks before commenting on these results.<\/p>\n
\"traceModelFree\"<\/a><\/p>\n
\"densityModelFree\"<\/a><\/p>\n
Mixing is fine and there is no obvious problem with the convergence, so let\u2019s return to the findings.<\/p>\n
The genetic model has an average lambda of about 0.18. On this scale 0.5 is the value that most people assume corresponding to ORTT<\/sub> = ORTM<\/sub>*ORTM<\/sub>. This value is not well supported by the analysis and if anything we are closer to lambda=0 when ORTM<\/sub>=1, which is known as the recessive model because we only see increased risk when both variants are T. Here the average effect of TT across the studies is to increase the odds of disease risk by 1.8 compared with MM individuals. The standard deviation is, however, quite large, which suggests that there is a lot of heterogeneity between studies. Perhaps, this variant interacts with something in the environment that differs across the populations or perhaps it is methodological, for instance it might relate to the inclusion criteria of the studies.<\/p>\n
This particular analysis assumed a Beta distribution, B(1,1), for the prior on lambda this is flat between 0 and 1 but importantly excludes the possibility that lambda is outside that range. So this prior rules out the possibility the ORTM<\/sub>> ORTT<\/sub>. In the paper we looked at the impact of several different priors and I am sure that you can see\u00a0other possibilities,\u00a0such as\u00a0assuming a fixed ORTT<\/sub> or for putting structure on the betas.<\/p>\n
In the paper we then simulated lots of datasets based on the pattern in the Kato study and investigated the fit. Suppose that we want to base our simulation on the\u00a0posterior mean values of the betas,<\/p>\n
use kato.dta, clear<\/strong><\/span><\/p>\n
gen tcase = case1+case2+case3<\/strong><\/span><\/p>\n
gen tcont = control3+control2+control1<\/strong><\/span><\/p>\n
input b1 b2 <\/strong><\/span><\/p>\n
18.4 28.6 <\/strong><\/span><\/p>\n
11.3 20.7<\/strong><\/span><\/p>\n
6.6 3.9<\/strong><\/span><\/p>\n
7.5 14.3<\/strong><\/span><\/p>\n
 5.8 11.6 <\/strong><\/span><\/p>\n
 6.5 6.9<\/strong><\/span><\/p>\n
 9.8 24.9<\/strong><\/span><\/p>\n
local lambda = 0.18<\/strong><\/span><\/p>\n
local sd = 0.49<\/strong><\/span><\/p>\n
local theta = 0.54<\/strong><\/span><\/p>\n
we can then simulate a new data set and write that to our data file<\/span><\/p>\n
gen p1 = 1\/(1+b1+b2)<\/strong><\/span><\/p>\n
gen p2 = b1\/(1+b1+b2)<\/strong><\/span><\/p>\n
gen p3 = b2\/(1+b1+b2)<\/strong><\/span><\/p>\n
gen d = rnormal(`theta’,`sd’)<\/strong><\/span><\/p>\n
gen OR2 = exp(d)<\/strong><\/span><\/p>\n
gen OR1 = exp(`lambda’*d)<\/strong><\/span><\/p>\n
gen q1 = 1\/(1+b1*OR1+b2*OR2)<\/strong><\/span><\/p>\n
gen q2 = b1*OR1\/(1+b1*OR1+b2*OR2)<\/strong><\/span><\/p>\n
gen q3 = b2*OR2\/(1+b1*OR1+b2*OR2)<\/strong><\/span><\/p>\n
gen rcase1 = rbinomial(tcase,q1)<\/strong><\/span><\/p>\n
gen rcase2 = rbinomial(tcase-rcase1,q2\/(q2+q3))<\/strong><\/span><\/p>\n
gen rcase3 = tcase-rcase1-rcase2<\/strong><\/span><\/p>\n
gen rcont1 = rbinomial(tcont,p1)<\/strong><\/span><\/p>\n
gen rcont2 = rbinomial(tcont-rcont1,p2\/(p2+p3))<\/strong><\/span><\/p>\n
gen rcont3 = tcont-rcont1-rcont2<\/strong><\/span><\/p>\n
wbslist (struct rcase1 rcase2 rcase3 , name(ncase) format(%4.0f) ) \/\/\/<\/strong><\/span><\/p>\n
\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (struct rcont1 rcont2 rcont3 , name(ncont) format(%4.0f) ) \/\/\/<\/strong><\/span><\/p>\n
\u00a0\u00a0 \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 (vect tcase tcont) using data.txt, replace<\/strong><\/span><\/p>\n
This code creates a new dataset based on the Kato study and write it to the text file data.txt.\u00a0We could\u00a0place\u00a0this code in\u00a0a loop so that we generate lots of datasets and each time we could run OpenBUGS, read the results and store whatever we choose. Then we could run the do file and go and have tea, because it will take a while if we want a lot of simulations (we\u00a0analysed\u00a01,000 random datasets\u00a0for each of\u00a0several different choices of prior distribution).<\/p>\n","protected":false},"excerpt":{"rendered":"
Some years ago some colleagues and I wrote a paper about the Bayesian meta-analysis of genetic association studies (Minelli et al Statistics in Medicine 2005;24:3845-3861) and periodically since then people have contacted me asking for the WinBUGS code that we used in that paper and those requests have prompted me to reformulate\u00a0the code\u00a0in terms of […]<\/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":[],"class_list":["post-158","post","type-post","status-publish","format-standard","hentry","category-uncategorized"],"_links":{"self":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/158","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=158"}],"version-history":[{"count":8,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/158\/revisions"}],"predecessor-version":[{"id":171,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/posts\/158\/revisions\/171"}],"wp:attachment":[{"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/media?parent=158"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/categories?post=158"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/staffblogs.le.ac.uk\/bayeswithstata\/wp-json\/wp\/v2\/tags?post=158"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}