# Model for all data model { for (i in 1:49){ illit[(i-1)*3+1] ~ dbin(pw[i],total[(i-1)*3+1]) illit[(i-1)*3+2] ~ dbin(pf[i],total[(i-1)*3+2]) illit[(i-1)*3+3] ~ dbin(pb[i],total[(i-1)*3+3]) logit(pw[i]) <-beta0w + gammaw*x[i] + b[i,1] logit(pf[i]) <-beta0f + gammaf*x[i] + b[i,2] logit(pb[i]) <-beta0b + gammab*x[i] + b[i,3] b[i,1:3] ~ dmnorm(mu[1:3],W[1:3,1:3]) } #bstar[1:3] ~ dmnorm(mu[1:3],W[1:3,1:3]) W[1:3, 1:3] ~ dwish(R[1:3, 1:3], 3) D[1:3, 1:3] <- inverse(W[1:3, 1:3]) for (i in 1 : 3) {sigma[i] <- sqrt(D[i, i]) } fi[1] <- exp(beta0w) fi[2] <- exp(beta0f) fi[3] <- exp(beta0b) fi[4] <- beta0w fi[5] <- beta0f fi[6] <- beta0b fi[7] <- beta0f-beta0w fi[8] <- beta0b-beta0w fi[9] <- exp(gammaw) fi[10] <- exp(gammaf) fi[11] <- exp(gammab) fi[12] <- gammaw fi[13] <- gammaf fi[14] <- gammab fi[15] <- exp(-1.96*sigma[1]) fi[16] <- exp(1.96*sigma[1]) fi[17] <- exp(-1.96*sigma[2]) fi[18] <- exp(1.96*sigma[2]) fi[19] <- exp(-1.96*sigma[3]) fi[20] <- exp(1.96*sigma[3]) fi[21] <- exp(beta0w-1.96*sigma[1]) fi[22] <- exp(beta0w+1.96*sigma[1]) fi[23] <- exp(beta0f-1.96*sigma[2]) fi[24] <- exp(beta0f+1.96*sigma[2]) fi[25] <- exp(beta0b-1.96*sigma[3]) fi[26] <- exp(beta0b+1.96*sigma[3]) fi[27] <- exp(beta0w+gammaw-1.96*sigma[1]) fi[28] <- exp(beta0w+gammaw+1.96*sigma[1]) fi[29] <- exp(beta0f+gammaf-1.96*sigma[2]) fi[30] <- exp(beta0f+gammaf+1.96*sigma[2]) fi[31] <- exp(beta0b+gammab-1.96*sigma[3]) fi[32] <- exp(beta0b+gammab+1.96*sigma[3]) beta0w ~ dflat() beta0f ~ dflat() beta0b ~ dflat() gammaw ~ dflat() gammaf ~ dflat() gammab ~ dflat() } INIT 1 list(beta0w=0,beta0f=0,beta0b=0,gammaw=0,gammaf=0,gammab=0,W = structure(.Data = c(1, 0, 0,0, 1, 0,0, 0, 1), .Dim = c(3, 3)),b=structure(.Data = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 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, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), .Dim = c(49, 3))) INIT 2 list(beta0w=-5.32,beta0f=-2.58,beta0b=-3.21,gammaw=1.52,gammaf=0.03,gammab=1.22,W = structure(.Data = c(1.62, 0, 0,0, 3.62, 0,0, 0, 3.64), .Dim = c(3, 3)),b=structure(.Data = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), .Dim = c(49, 3))) list(mu=c(0,0,0),R = structure(.Data = c(0.033, 0,0,0,0.033, 0,0, 0, 0.033), .Dim = c(3, 3)), illit=c(39252,11183,16532,1763,4649,11,1639,2422,12,3762,6924,450,7001,4113,3228,1896,2392,3496,13202,11589,25073,533,1411,4591,65482,1738,95148,42476,9788,10173,93372,450,139105,36246,297,156065,46898,554,163237,14478,3159,65167,101695,1267,28553,87406,754,57251,60959,1335,188673,20070,882,177605,35890,666,60102,71903,6677,139393,27796,1479,12560,46878,7136,90225,932,3085,52,1151,1198,25,381,811,47,5807,7331,403,18733,530,140,960,551,366,952,1547,30,130,909,7,3031,7103,174,2043,3743,49,9840,45600,2148), total=c(538563,98741,896,338000,81458,692,250846,41736,449,2413000,1042692,42648,390286,168699,7840,875250,377644,23612,6930800,3160602,344160,2314400,830946,170804,6086167,1220403,354952,4626715,639060,253328,2501111,134020,93417,4714000,1210933,279000,2958000,821403,140033,2052167,384103,8886,1811000,389744,8000,2044250,164778,14389,2616800,149107,187864,440750,105659,324,409750,65459,545,940500,115400,11538,1400200,69712,54712,158000,16727,26485,1015539,94992,219939,266500,29396,111976,1364208,23486,495563,1148000,50979,90027,1667357,8654,675267,710706,5211,580167,1421152,13850,1305896,762000,58500,749046,1784123,21845,185409,1618630,13000,384235,1269979,15706,720126,743333,7000,765539,1025429,10091,373304,984973,34776,598253,1635059,26411,135054,3348429,97753,673321,310667,71744,1130,287750,29950,595,127000,19310,1119,725875,85244,10333,243286,7794,2333,192000,15306,9150,317333,42972,937,65000,12120,467,1010333,244931,6000,681000,103972,1960,3280000,800000,69290 ),x=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,0,1,1,0,0,0,0,0) # NonJC states only model { for (i in 1:27){ illit[(i-1)*3+1] ~ dbin(pw[i],total[(i-1)*3+1]) illit[(i-1)*3+2] ~ dbin(pf[i],total[(i-1)*3+2]) illit[(i-1)*3+3] ~ dbin(pb[i],total[(i-1)*3+3]) logit(pw[i]) <-beta0w + b[i,1] logit(pf[i]) <-beta0f + b[i,2] logit(pb[i]) <-beta0b + b[i,3] b[i,1:3] ~ dmnorm(mu[1:3],W[1:3,1:3]) } #bstar[1:3] ~ dmnorm(mu[1:3],W[1:3,1:3]) W[1:3, 1:3] ~ dwish(R[1:3, 1:3], 3) D[1:3, 1:3] <- inverse(W[1:3, 1:3]) for (i in 1 : 3) {sigma[i] <- sqrt(D[i, i]) } fi[1] <- exp(beta0w) fi[2] <- exp(beta0f) fi[3] <- exp(beta0b) fi[4] <- beta0w fi[5] <- beta0f fi[6] <- beta0b fi[7] <- beta0f-beta0w fi[8] <- beta0b-beta0w fi[15] <- exp(-1.96*sigma[1]) fi[16] <- exp(1.96*sigma[1]) fi[17] <- exp(-1.96*sigma[2]) fi[18] <- exp(1.96*sigma[2]) fi[19] <- exp(-1.96*sigma[3]) fi[20] <- exp(1.96*sigma[3]) fi[21] <- exp(beta0w-1.96*sigma[1]) fi[22] <- exp(beta0w+1.96*sigma[1]) fi[23] <- exp(beta0f-1.96*sigma[2]) fi[24] <- exp(beta0f+1.96*sigma[2]) fi[25] <- exp(beta0b-1.96*sigma[3]) fi[26] <- exp(beta0b+1.96*sigma[3]) beta0w ~ dflat() beta0f ~ dflat() beta0b ~ dflat() } INIT 1 list(beta0w=0,beta0f=0,beta0b=0,gammaw=0,gammaf=0,gammab=0,W = structure(.Data = c(1, 0, 0,0, 1, 0,0, 0, 1), .Dim = c(3, 3)),b=structure(.Data = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 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, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), .Dim = c(49, 3))) INIT 2 list(beta0w=-5.32,beta0f=-2.58,beta0b=-3.21,W = structure(.Data = c(1.62, 0, 0,0, 3.62, 0,0, 0, 3.64), .Dim = c(3, 3)),b=structure(.Data = c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), .Dim = c(27, 3))) DATA FOR NON JC list(mu=c(0,0,0),R = structure(.Data = c(0.033, 0,0,0,0.033, 0,0, 0, 0.033), .Dim = c(3, 3)), illit=c(8617,8393,43,2366,7820,27,3261,3005,22,9652,111568,2303,2732,24124,635,3501,55136,1157,34654,341345,8604,11572,107192,8711,36517,187942,14908,32387,74131,16213,22510,13536,5605,28284,108984,10044,14790,55034,4201,12313,29960,391,7244,16759,160,8177,5932,777,1763,4649,11,1639,2422,12,3762,6924,450,932,3085,52,1151,1198,25,5807,7331,403,952,1547,30,130,909,7,3031,7103,174,2043,3743,49,9840,45600,2148), total=c(538563,98741,896,338000,81458,692,250846,41736,449,2413000,1042692,42648,390286,168699,7840,875250,377644,23612,6930800,3160602,344160,2314400,830946,170804,6086167,1220403,354952,4626715,639060,253328,2501111,134020,93417,4714000,1210933,279000,2958000,821403,140033,2052167,384103,8886,1811000,389744,8000,2044250,164778,14389,440750,105659,324,409750,65459,545,940500,115400,11538,310667,71744,1130,287750,29950,595,725875,85244,10333,317333,42972,937,65000,12120,467,1010333,244931,6000,681000,103972,1960,3280000,800000,69290)) # JC states only model { for (i in 1:22){ illit[(i-1)*3+1] ~ dbin(pw[i],total[(i-1)*3+1]) illit[(i-1)*3+2] ~ dbin(pf[i],total[(i-1)*3+2]) illit[(i-1)*3+3] ~ dbin(pb[i],total[(i-1)*3+3]) logit(pw[i]) <-beta0w + b[i,1] logit(pf[i]) <-beta0f + b[i,2] logit(pb[i]) <-beta0b + b[i,3] b[i,1:3] ~ dmnorm(mu[1:3],W[1:3,1:3]) } #bstar[1:3] ~ dmnorm(mu[1:3],W[1:3,1:3]) W[1:3, 1:3] ~ dwish(R[1:3, 1:3], 3) D[1:3, 1:3] <- inverse(W[1:3, 1:3]) for (i in 1 : 3) {sigma[i] <- sqrt(D[i, i]) } fi[1] <- exp(beta0w) fi[2] <- exp(beta0f) fi[3] <- exp(beta0b) fi[4] <- beta0w fi[5] <- beta0f fi[6] <- beta0b fi[7] <- beta0f-beta0w fi[8] <- beta0b-beta0w fi[15] <- exp(-1.96*sigma[1]) fi[16] <- exp(1.96*sigma[1]) fi[17] <- exp(-1.96*sigma[2]) fi[18] <- exp(1.96*sigma[2]) fi[19] <- exp(-1.96*sigma[3]) fi[20] <- exp(1.96*sigma[3]) fi[21] <- exp(beta0w-1.96*sigma[1]) fi[22] <- exp(beta0w+1.96*sigma[1]) fi[23] <- exp(beta0f-1.96*sigma[2]) fi[24] <- exp(beta0f+1.96*sigma[2]) fi[25] <- exp(beta0b-1.96*sigma[3]) fi[26] <- exp(beta0b+1.96*sigma[3]) beta0w ~ dflat() beta0f ~ dflat() beta0b ~ dflat() } INIT 1 list(beta0w=0,beta0f=0,beta0b=0,gammaw=0,gammaf=0,gammab=0,W = structure(.Data = c(1, 0, 0,0, 1, 0,0, 0, 1), .Dim = c(3, 3)),b=structure(.Data = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 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, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), .Dim = c(49, 3))) INIT 2 list(beta0w=-5.32,beta0f=-2.58,beta0b=-3.21,W = structure(.Data = c(1.62, 0, 0,0, 3.62, 0,0, 0, 3.64), .Dim = c(3, 3)),b=structure(.Data = c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ), .Dim = c(22, 3))) # JC DATA list(mu=c(0,0,0),R = structure(.Data = c(0.033, 0,0,0,0.033, 0,0, 0, 0.033), .Dim = c(3, 3)), illit=c(39252,11183,16532,7001,4113,3228,1896,2392,3496,13202,11589,25073,533,1411,4591,65482,1738,95148,42476,9788,10173,93372,450,139105,36246,297,156065,46898,554,163237,14478,3159,65167,101695,1267,28553,87406,754,57251,60959,1335,188673,20070,882,177605,35890,666,60102,71903,6677,139393,27796,1479,12560,46878,7136,90225,381,811,47,18733,530,140,960,551,366), total=c(2616800,149107,187864,1400200,69712,54712,158000,16727,26485,1015539,94992,219939,266500,29396,111976,1364208,23486,495563,1148000,50979,90027,1667357,8654,675267,710706,5211,580167,1421152,13850,1305896,762000,58500,749046,1784123,21845,185409,1618630,13000,384235,1269979,15706,720126,743333,7000,765539,1025429,10091,373304,984973,34776,598253,1635059,26411,135054,3348429,97753,673321,127000,19310,1119,243286,7794,2333,192000,15306,9150))