##### Section 2: WinBUGS model for ovarian cancer data model{ # Model for the outcomes for both tests for(i in 1:51){ ty0CA199[i] <- log(y0CA199[i]) # Transform controls ty0CA125[i] <- log(y0CA125[i]) ty0CA199[i] ~ dnorm(mu0199,tau0199) # Log normal model ty0CA125[i] ~ dnorm(mu0125,tau0125) # for controls } for(i in 1:90){ ty1CA199[i] <- log(y1CA199[i]) # Transform cases ty1CA125[i] <- log(y1CA125[i]) ty1CA199[i] ~ dnorm(mu1199,tau1199) # Log normal model ty1CA125[i] ~ dnorm(mu1125,tau1125) # for cases } # Diffuse priors mu0199~dnorm(0,0.01) tau0199~dgamma(0.01,0.01) mu0125~dnorm(0,0.01) tau0125~dgamma(0.01,0.01) mu1199~dnorm(0,0.01) tau1199~dgamma(0.01,0.01) mu1125~dnorm(0,0.01) tau1125~dgamma(0.01,0.01) AUC199 <- phi( (mu1199 - mu0199)/sqrt(1/tau0199 + 1/tau1199)) AUC125 <- phi( (mu1125 - mu0125)/sqrt(1/tau0125 + 1/tau1125)) for(i in 1:101){ # Get ROC values over the grid x[] = phi^{-1}(t[]) ROC199[i]<- phi(sqrt(tau1199)*(mu1199-mu0199)+sqrt(tau1199/tau0199)*x[i]) ROC125[i]<- phi(sqrt(tau1125)*(mu1125-mu0125)+sqrt(tau1125/tau0125)*x[i]) }# Predictive probabilities of infection over the grid y[] for(i in 1:321){ f1199[i] <- sqrt(tau1199)*exp((-tau1199/2)*pow(y[i]-mu1199,2)) f0199[i] <- sqrt(tau0199)*exp((-tau0199/2)*pow(y[i]-mu0199,2)) f1125[i] <- sqrt(tau1125)*exp((-tau1125/2)*pow(y[i]-mu1125,2)) f0125[i] <- sqrt(tau0125)*exp((-tau0125/2)*pow(y[i]-mu0125,2)) probMale125[i] <- f1125[i]*pm/(f1125[i]*pm + f0125[i]*(1-pm)) probFemale125[i] <- f1125[i]*pf/(f1125[i]*pf + f0125[i]*(1-pf)) probMale199[i] <- f1199[i]*pm/(f1199[i]*pm + f0199[i]*(1-pm)) probFemale199[i] <- f1199[i]*pf/(f1199[i]*pf + f0199[i]*(1-pf)) } } ##### Section 3: WinBUGS model, data, and inits for the regression example with continuous y model{ for(i in 1:n){ z[i] ~ dbern(pi[i]) y[i] ~ dnorm(mmu[i], ttau[i]) mmu[i] <- z[i]*mu[1] + (1-z[i])*mu[2] ttau[i] <- z[i]*tau[1] + (1-z[i])*tau[2] logit(pi[i]) <- beta[1] + beta[2]*x[i] } tildepi[1] ~ dbeta(a[1],b[1]) tildepi[2] ~ dbeta(a[2],b[2]) beta[1] <- logit(tildepi[1]) beta[2] <- logit(tildepi[2]) - logit(tildepi[1]) for(i in 1:2){ mu[i] ~ dnorm(d[i],0.001) tau[i] ~ dgamma(0.001,0.001) } v[1] ~ dnorm(mu[1],tau[1]) #Pred Dens Popn 1 v[2] ~ dnorm(mu[2],tau[2]) #Pred Dens Popn 0 # Get predictive probabilities of \$D\$ for(i in 1:n){ f1[i] <- sqrt(tau[1])*exp((-tau[1]/2)*pow(y[i]-mu[1],2)) f2[i] <- sqrt(tau[2])*exp((-tau[2]/2)*pow(y[i]-mu[2],2)) prob[i] <- f1[i]*pi[i]/(f1[i]*pi[i] + f2[i]*(1-pi[i])) } } list(mu=c(55,45), tau=c(1,1), tildepi=c(.95,.05)) list(z=c(0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)) list(n=300,a=c(19,1),b=c(1,19),d=c(20,10), x = c(0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), y = c(46.9688885299112, 52.3361408968752, 55.0587095330899, 67.063801419522, 68.8339157013614, 35.1997003210818, 43.8434567496886, 62.5678556871719, 51.7401497932518, 40.0210644491513, 55.5803936018858, 77.8540400706143, 53.2577006714074, 60.978096870618, 44.8963836121789, 45.8611825455399, 31.8055612025579, 44.0913994633267, 57.722105348682, 74.6109500714428, 52.7449019892315, 39.1355235642223, 46.0213939033561, 60.391334720546, 55.9977718278595, 44.1538782443230, 39.8951487485588, 57.8750840531612, 52.741948596222, 62.5174297594526, 39.7587406576789, 42.1347008412321, 29.9739391852203, 57.737985500304, 60.9027770410002, 40.8896197976564, 39.3549237757351, 42.9913777609811, 47.4636851909366, 40.19975650674, 46.8049462243584, 24.9840101346073, 53.913147855837, 31.7877258787095, 46.8465351493863, 44.6579914658539, 65.5782255615336, 47.3727177437165, 61.4277824354533, 41.3232437613357, 35.55101258649, 47.2046690388387, 54.7556206451842, 45.41255030363, 42.8147384529643, 22.7165793798118, 46.2008902700476, 55.2319677884315, 53.0094398983507, 48.3747036387213, 66.3509396059549, 43.739428059669, 83.3584841203558, 63.8976098188874, 42.9241764802072, 47.580661251016, 66.5911443421293, 57.8933498142899, 58.7485101903814, 41.1459589104188, 45.0476882605214, 48.17446507216, 41.1737178765139, 42.6121526136231, 48.7625758570532, 34.7531217958322, 47.8557167896222, 58.0750301395481, 32.8325393816697, 47.1903811877905, 40.163918721602, 45.2280622030011, 56.0167740512132, 43.8191866357756, 54.5691711543537, 37.7168944517151, 34.9276809305544, 62.3866423609381, 50.1037075520153, 58.4355197239564, 40.2671573484903, 60.7805867474096, 42.8168780864038, 55.4999551097814, 59.9444187794951, 44.5238375557945, 64.1243050143463, 24.5124935705296, 55.0657551062564, 58.9778393827595, 29.1429226326105, 31.4571459000624, 36.7874450804490, 34.5896554597068, 40.750392806842, 30.0467801951236, 40.7489409160689, 36.3809046657188, 58.2365214172514, 49.4283022810525, 33.9009170390867, 48.8656471648306, 33.020910338002, 41.4326955179173, 51.7548175724649, 39.979249177136, 55.3342775942271, 47.0222040653348, 59.6656425570817, 46.8155625838758, 52.9454435060603, 41.9878708909172, 26.4391187698536, 46.7718284534991, 46.8378047376953, 41.1191896022904, 42.2769904851212, 32.5460505660805, 55.2361044591861, 27.2078694980604, 35.9630078835694, 33.2292633602079, 41.6574609333793, 27.3834733152, 36.9821495866163, 54.1524350368548, 36.2217513370123, 31.8023503054429, 48.5910824045674, 36.3423467217617, 23.7152203775348, 51.896720662421, 49.0922660578311, 61.939437189927, 48.6708629284579, 43.7515213645761, 31.6912919345410, 42.6274611575775, 47.4228070397373, 44.5473638282546, 39.5870019393295, 50.8869434593182, 50.6223583666625, 40.5110515231160, 32.7256715242291, 35.2715751255559, 34.3933905356635, 32.004668642066, 41.4421841859785, 36.4075760228735, 43.3017576155405, 43.2893845266867, 36.035282852267, 42.4105732554219, 43.0190461497302, 32.312692090048, 28.2851310897362, 39.3906583601767, 37.9559659430498, 36.9724028427918, 28.5775671471887, 26.4029930319749, 25.8990848376682, 30.753554487765, 38.3609784614964, 35.7225535772337, 60.3281246520733, 41.9787824083943, 52.4280506016782, 40.3210195605944, 43.4956744280586, 46.2467042503358, 41.5944337046651, 28.0434248253336, 43.257898810863, 50.4003324736337, 37.4314642486052, 39.5386248666131, 29.8912644825252, 22.769206083817, 44.0710870114558, 28.0969782004725, 48.1169564193978, 45.7481824181866, 27.9146110824603, 39.0397759395597, 36.4960026167076, 40.1374871963424, 34.5431268184916, 46.7976471195881, 33.6964177698529, 58.0351050129309, 43.6274635621055, 26.5279249236766, 35.1883181962994, 47.6689047578843, 52.4641974360265, 49.7837014026475, 32.9192405411713, 57.4101669218539, 37.1894216284287, 34.8265181312403, 45.4235939524782, 23.0698515464671, 28.6360956987204, 45.3501290192816, 42.5864893361019, 47.2288184942705, 26.7834412894001, 50.1376097691537, 41.9241131506460, 34.9265683332031, 34.0555733161908, 42.1871614544994, 58.7520051047603, 20.6597945839066, 37.0059402377429, 55.1002441280777, 43.6682558700408, 44.0347669147137, 50.9780002317417, 37.4367349743465, 28.3676745207623, 48.809809867479, 27.7786770201964, 48.9781051215074, 29.4726186109028, 30.9009637589159, 39.6642946948769, 53.2769007089319, 41.0543914126919, 47.4113125630637, 41.3876752720242, 46.5457869044804, 34.5549470388404, 49.5040805196856, 43.0760514235605, 36.1391362732467, 32.6672931431525, 52.4161653503593, 52.2073058573022, 55.0187630400466, 60.1764373110256, 47.7763123062371, 43.5228153909299, 36.5087925071478, 42.2058517214953, 36.5612145931967, 28.8107075612686, 38.4372181214345, 45.9028943324997, 50.3490752334236, 53.8057938466392, 34.661751608125, 23.2465775481229, 50.9503934503450, 49.1443814623675, 39.0080927086418, 44.7646345585957, 21.2309912343723, 45.1648900374699, 38.164479174638, 48.523777652212, 25.7757424033611, 36.4532815990225, 46.0198156332465, 50.95234528742, 37.071586659518, 24.5320355627188, 19.1501760022900, 53.4340318219755, 23.8494224018, 35.0166384167196, 44.5539754278086, 29.5616357081007, 42.7567609118962, 41.7929758492584, 56.5246962019851, 39.5328891883512, 24.4217078546884, 38.5956517992541, 34.5088896020053, 48.6931763462509, 40.7392340133278, 57.8453253696438, 39.5016466623566, 14.5639898571131, 44.3014497066879, 28.1193738972429, 48.283793128681)) ##### Section 4: WinBUGS model, data, and inits for the smoking cessation example model {for (i in 1:N) { y[i] ~ dbern(q[i]) # specification of the LR model with error q[i] <- pi[i]*Se+(1-pi[i])*(1-Sp) logit(pi[i]) <- beta[1] + x2[i]*beta[2] + x3[i]*beta[3] + x4[i]*beta[4] + x5[i]*beta[5] + x6[i]*beta[6] } Se ~ dbeta(99, 1) # the priors are specified for Se and Sp Sp ~ dbeta(14, 2) p[1] ~ dbeta(8, 15) # spec of the prior on 6 probs of smoking p[2] ~ dbeta(10, 15) # cessation for 6 covariate combinations p[3] ~ dbeta(3, 13) p[4] ~ dbeta(8, 10) p[5] ~ dbeta(4, 15) p[6] ~ dbeta(6, 15) # relate the regression coefficients to the # probabilities on which prior was specified # this results in an induced informative prior # on the regression coefficients beta[1] <- xinv[1,1]*logit(p[1]) + xinv[1,2]*logit(p[2]) + xinv[1,3]*logit(p[3]) + xinv[1,4]*logit(p[4]) + xinv[1,5]*logit(p[5]) + xinv[1,6]*logit(p[6]) beta[2] <- xinv[2,1]*logit(p[1]) + xinv[2,2]*logit(p[2]) + xinv[2,3]*logit(p[3]) + xinv[2,4]*logit(p[4]) + xinv[2,5]*logit(p[5]) + xinv[2,6]*logit(p[6]) beta[3] <- xinv[3,1]*logit(p[1]) + xinv[3,2]*logit(p[2]) + xinv[3,3]*logit(p[3]) + xinv[3,4]*logit(p[4]) + xinv[3,5]*logit(p[5]) + xinv[3,6]*logit(p[6]) beta[4] <- xinv[4,1]*logit(p[1]) + xinv[4,2]*logit(p[2]) + xinv[4,3]*logit(p[3]) + xinv[4,4]*logit(p[4]) + xinv[4,5]*logit(p[5]) + xinv[4,6]*logit(p[6]) beta[5] <- xinv[5,1]*logit(p[1]) + xinv[5,2]*logit(p[2]) + xinv[5,3]*logit(p[3]) + xinv[5,4]*logit(p[4]) + xinv[5,5]*logit(p[5]) + xinv[5,6]*logit(p[6]) beta[6] <- xinv[6,1]*logit(p[1]) + xinv[6,2]*logit(p[2]) + xinv[6,3]*logit(p[3]) + xinv[6,4]*logit(p[4]) + xinv[6,5]*logit(p[5]) + xinv[6,6]*logit(p[6]) } # inits list(Se = 0.95, Sp = 0.93, p = c(0.33, 0.39, 0.14, 0.44, 0.18, 0.26)) # data list(y=c(0,1,0,1,0,0,0,0,0,1,0,0,0,0,1,1,0,0,0,0,0,0,1,0,0,1,0,0,0,1,0, 0,0,0,1,0,0,1,0,0,0,1,0,1,1,1,0,1,0,1,0,0,0,1,0,0,0,0,1,0,1,0,1,1,0,0, 0,0,0,0,1,0,1,0,1,0,1,1,1,1,0,0,0,1,0,1,0,0,0,0,0,0,1,0,1,1,0,1,1,0,1, 0,0,0,0,1,0,1,1,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,1,1,0,0,1,0,1,0,0,0,1, 0,0,0,1,0,0,0,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0,1,1,0,0,1,0,0,0,1,0,0,0,0, 0,0,1,0,0,1,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,0,0,1, 1,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,0,0, 1,0,0,0,1,0,1,0,0,0,0,1,1,1,0,0,0,1,1,0,0,0,1,0,1,1,0,0,0,0,1,0,1,1,0, 0,0,1,1,0,1,0,1,1,0,0,1,0,0,1,1,0,0,1,1,0,1,0,0,0,0,0,0,1,1,1,0,1,0,0, 0,0,1,0,1,0,0,0,0,1,0,0,0,1,0,1,1,0,1,0,0,0,0,1,1,0,0,0,0,0,1,0,0,1,0, 0,0,1,0,0,0,0,0,0,0,0,1,0,1,1), x2 = c(0,1,1,1,1,1,1,0,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,1,0,0, 0,1,1,1,0,1,1,1,1,1,1,0,0,1,1,0,1,1,1,1,0,0,1,1,0,0,1,0,1,1,0,1,1,1,1, 1,1,1,1,0,0,1,1,1,1,1,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,0,1,0,1,1,1,1,1, 0,0,0,0,1,0,1,1,1,1,1,1,0,1,1,1,1,0,0,0,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1, 0,0,1,1,1,0,1,1,1,1,0,1,1,1,1,0,0,1,1,0,1,1,1,1,0,0,1,1,0,0,1,1,0,0,1, 1,1,1,0,0,1,1,0,1,0,1,1,0,1,1,1,1,0,1,1,1,0,1,0,1,0,0,1,1,1,1,1,1,0,1, 1,1,1,0,1,0,0,0,1,1,0,1,1,0,1,1,1,0,1,0,1,1,0,1,0,1,0,0,0,1,1,1,1,1,1, 1,1,1,1,0,1,1,1,1,0,1,1,1,0,0,1,1,1,1,1,1,1,0,1,1,0,1,1,0,1,0,0,1,0,1, 0,0,0,0,0,1,0,1,0,1,1,0,1,1,0,1,1,1,0,0,0,1,1,1,1,1,0,1,1,1,0,0,0,0,0, 0,1,0,0,1,1,0,1,0,1,1,1,0,0,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,0,0,1,1,0,1,1,1,0,0,0,0,1), x3 = c(1,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,1,0, 1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,1,0,1,0,0,0,0,0,0,0, 0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0, 1,0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0, 1,1,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,1,0,0,0,0,0,1,0,0,1,1,0,0,1,1,0, 0,0,0,1,1,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,1,1,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0, 1,1,0,1,0,0,1,0,1,0,0,0,0,0,1,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,1,0,1, 1,0,1,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,1,0,0,0,0,0,0,0,0,0,1,0,0), x4 = c(0,1,1,1,0,1,1,1,0,1,1,0,1,0,0,1,0,0,0,0,1,1,1,0,0,0,1,1,1,1,1,0, 0,0,1,0,0,0,1,1,0,1,0,0,0,1,1,0,0,0,1,0,1,0,0,1,1,0,0,1,0,1,0,0,1,1,1, 0,1,1,0,0,0,1,1,1,1,1,0,1,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0,1,1,0,0,1,1,1, 0,0,0,1,0,1,1,1,0,1,0,0,1,1,0,0,0,0,0,1,1,1,0,0,0,1,1,0,0,1,1,0,0,0,1, 1,1,1,1,0,0,0,1,1,0,1,0,0,0,1,1,0,0,1,0,0,1,0,0,0,0,0,0,0,1,1,1,0,1,0, 0,1,1,0,1,0,0,0,1,0,0,1,0,0,0,1,0,0,1,0,1,0,1,1,0,0,1,1,1,1,1,0,0,0,0, 1,0,1,0,0,0,0,0,1,0,0,0,1,0,0,1,1,0,1,0,1,0,0,0,0,1,0,0,0,0,1,0,1,1,1, 1,0,0,1,0,1,1,0,1,0,1,1,1,1,0,1,0,0,1,0,0,1,1,1,1,1,0,1,1,0,1,1,0,0,1, 1,0,0,1,1,0,0,0,0,0,1,1,0,0,0,1,0,1,0,0,1,0,1,0,1,0,1,1,1,1,0,1,1,0,0, 1,1,0,0,1,1,0,1,0,1,1,1,1,0,1,1,0,1,1,0,0,1,1,0,1,1,1,0,1,1,1,0,1,1,0, 1,0,1,1,1,1,1,1,0,0,0,1,0,1), x5 = c(0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,1,1,1,0,0,0,0,1,1,0,0,0,0,0,0, 1,0,0,1,0,1,0,0,1,0,1,0,0,0,0,0,0,1,0,0,0,0,1,0,0,1,1,0,0,0,0,0,0,0,0, 1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,1,1,1,0,0,1,1,0,0,1,1,0,0,0, 1,0,1,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,1,0,0,1,1,1,0, 0,0,0,0,0,0,1,0,0,1,0,1,0,1,0,0,0,0,0,1,1,0,1,0,0,1,1,1,1,0,0,0,0,0,0, 0,0,0,1,0,1,0,1,0,0,0,0,0,0,1,0,1,0,0,1,0,1,0,0,0,1,0,0,0,0,0,1,0,0,0, 0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0, 0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0, 0,1,0,0,0,0,1,1,1,0,0,0,1,1,1,0,1,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1, 0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,1,0,0,1,0,0,0,1,0,0,0,0,0,0,1, 0,1,0,0,0,0,0,0,0,0,0,0,0,0), x6 = c(0,1,1,1,0,1,1,0,1,1,1,0,1,1,0,1,0,1,1,1,1,0,1,0,1,1,1,1,1,1,1,0, 0,0,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,0,1,1,1,1,0,0,1,0,0,0,1,1,0,1,0, 1,0,1,1,0,1,1,0,1,1,1,1,0,1,1,1,1,1,1,0,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1, 1,1,0,0,1,0,1,0,1,1,1,0,1,0,1,0,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1, 0,1,1,1,0,0,1,1,0,1,1,1,0,0,1,1,1,1,1,1,0,1,0,0,0,0,1,0,0,1,1,0,0,0,0, 0,0,1,1,1,1,0,0,1,0,0,1,1,0,1,1,0,0,1,1,1,1,1,0,0,0,1,1,1,0,1,0,1,1,1, 1,0,1,1,1,1,0,1,1,0,1,1,1,0,0,0,1,0,0,1,1,1,1,1,1,1,1,1,1,1,0,1,1,0,1, 0,1,1,1,0,1,1,1,1,1,0,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0,1,1,0,1,0,0,1,0,0, 1,1,1,0,1,1,1,1,1,0,1,0,1,1,1,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0, 1,1,1,1,0,0,1,0,1,1,0,0,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,0,1,1,0,0, 0,1,1,1,1,0,0,1,1,1,1,1,1,1), xinv=structure(.Data=c(-1.0, 1.94E-16, 3.33E-16, 2.57E-16, 1.0, 1.0, -7.77E-16, 2.50E-16, 1.0, -4.86E-17, 3.33E-16, 1.0,-1.0, 1.0, 1.0, -3.75E-16, -3.33E-16, 1.0,1.0, -5.83E-16, -5.55E-16, -9.02E-17, -1.0, 8.88E-16,4.16E-16, -6.11E-16, -3.89E-16, 1.0, -1.0, 5.41E-16,1.0, 1.67E-16, -1.0, 3.33E-16, 5.55E-17, 7.77E-16),.Dim=c(6,6)), N=361) ##### Section 5: WinBUGS model, data, and inits for classical swine fever data model{ x[1:4] ~ dmulti(p[1:4], n) p[1] <- pi*(Sefat1*Sefat2+covDp) + (1-pi)*((1-Spfat1)*(1-Spfat2)+covDn) p[2] <- pi*(Sefat1*(1-Sefat2)-covDp) + (1-pi)*((1-Spfat1)*Spfat2-covDn) p[3] <- pi*((1-Sefat1)*Sefat2-covDp) + (1-pi)*(Spfat1*(1-Spfat2)-covDn) p[4] <- pi*((1-Sefat1)*(1-Sefat2)+covDp) + (1-pi)*(Spfat1*Spfat2+covDn) ls <- (Sefat1-1)*(1-Sefat2) us <- min(Sefat1,Sefat2) - Sefat1*Sefat2 lc <- (Spfat1-1)*(1-Spfat2) uc <- min(Spfat1,Spfat2) - Spfat1*Spfat2 pi ~ dbeta(13.322, 6.281) ### Mode=0.70, 95% sure > 0.50 Sefat1 ~ dbeta(9.628,3.876) ### Mode=0.75, 95% sure > 0.50 Spfat1 ~ dbeta(15.034, 2.559) ### Mode=0.90, 95% sure > 0.70 Sefat2 ~ dbeta(9.628, 3.876) ### Mode=0.75, 95% sure > 0.50 Spfat2 ~ dbeta(15.034, 2.559) ### Mode=0.90, 95% sure > 0.70 covDn ~ dunif(lc, uc) covDp ~ dunif(ls, us) rhoD <- covDp / sqrt(Sefat1*(1-Sefat1)*Sefat2*(1-Sefat2)) rhoDc <- covDn / sqrt(Spfat1*(1-Spfat1)*Spfat2*(1-Spfat2)) } list(n=214, x=c(121,6,16,71)) list(pi=0.7, Sefat1=0.75, Spfat1=0.90, Sefat2=0.75, Spfat2=0.90) ##### Section 5: WinBUGS model, data, and inits for brucellosis data model{ for(i in 1:K){ x[i, 1:3] ~ dmulti(p[i, 1:3], n[i]) p[i,1] <- pi[i]*(Sebapa*Seriv+covDp) + (1-pi[i])*((1-Spbapa)*(1-Spriv)+covDn) p[i,2] <- pi[i]*(Sebapa*(1-Seriv)-covDp) + (1-pi[i])*((1-Spbapa)*Spriv-covDn) p[i,3] <- 1-p[i,1]-p[i,2] pi[i] ~ dbeta(alpha, beta) } Sebapa ~ dbeta(88.3,1.9) ## Mode=0.99, 95% sure > 0.95 Spbapa ~ dbeta(13.3,6.3) ## Mode=0.70, 95% sure > 0.50 Seriv ~ dbeta(130.7,15.4) ## Mode=0.90, 95% sure > 0.85 Spriv ~ dbeta(99.7,6.2) ## Mode=0.95, 95% sure > 0.90 alpha <- mu*psi beta <- psi*(1-mu) mu ~ dbeta(16.9,48.6) ## Mode=0.25; 95% sure < 0.35 psi ~ dgamma(7.23, 1.28) ls <- (Sebapa-1)*(1-Seriv) lc <- (Spbapa-1)*(1-Spriv) us <- min(Sebapa,Seriv)-Sebapa*Seriv uc <- min(Spbapa, Spriv) - Spbapa*Spriv covDn ~ dunif(lc, uc) covDp ~ dunif(ls, us) rhoDp <- covDp / sqrt(Sebapa*(1-Sebapa)*Seriv*(1-Seriv)) rhoDn <- covDn / sqrt(Spbapa*(1-Spbapa)*Spriv*(1-Spriv)) pi21 ~ dbeta(alpha,beta) a[1] <- 1-step(pi21-0.50) a[2] <- 1-step(pi21-0.10) a[3] <- step(pi21-0.90) a[4] <- step(pi21-0.75) } list( Sebapa=0.99, Spbapa=0.70, Seriv=0.90, Spriv=0.95, mu=0.25, psi=7) list(K=20, n=c(67,12,22,71,18,80,147,45,37,118,12,62,11,20,56,60,26,54,17,17), x=structure(.Data=c(4,2,61, 1,1,10, 2,2,18, 2,0,69, 2,0,16, 2,0,78, 80,5,62, 4,5,36, 10,0,27, 10,12,96, 1,0,11, 1,1,60, 3,0,8, 3,7,10, 38,0,18, 1,0,59, 3,0,23, 13,3,38, 7,0,10, 1,1,15),.Dim=c(20,3)))