#first the legal disclaimer # Program name: SpatioTempProgramMurphy.txt # Copyright Yale Pepper Center (2008) author: Terrence E. Murphy # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # You should have received a copy of the GNU General Public License # along with this program. If not, see . # The following Winbugs code was used to analyze the fall-related utilizations of #hospitals and emergency departments in the itreatment and usual care arms of the #Connecticut Collaboration for Fall Prevention. The results of this analysis were # published in the New England Journal of Medicine on July 17, 2008, vol 359: 252-261. # This Bayesian hierarchical model operationalizes an enhancement of the model # suggested by Waller et al in the Journal of the American Statistical Association (1997) # vol 92: 607-617 that was described in Murphy, Tinetti and Allore (2008) in # Contemporary Clinical Trials volume 29: 343-350. This model nests the 8 specific strata, # called ZAS for zcta-age-sex, defined by age group (4 levels) and sex (2 levels) within each zip code tabulation # area (ZCTA) in the two study arms. The treatment arm consists of 58 ZCTAs in the # greater Hartford region and the usual care arm of 53 ZCTAs along the southeastern # coast. Because the units of analysis are spatial in nature (ZCTAs), we used aggregate # statistical information.The stratification served to remove bias from the # two risk factors most strongly associated with probability of fall-related injury and also # raised the sample size. #The enhancement of Waller et al used here is the use of a random error term for #non-spatial heterogeneity for each stratum, in concert with the spatial residual #term calculated for each ZCTA that was the same for each stratum of each #individual ZCTA. # The following code and corresponding data are an incomplete sample of the # CCFP data and will not yield the same results as the first citation listed above. # They are provided to serve as simple examples of code and data structure # for gerontological researchers. model { ################################################################# ## we define model for ZCTAs 1-20 which are a sample of the usual care area ZCTAs ## There are 160 ZAS strata because 20 x 8 = 160 ## T is the number of six month periods of time during the study, which was ## ten because the intervention and evaluation periods consisted of ## three and two years each respectively, with data grouped in ## six month intervals to represent the shortest length of time when it ## could be reasonably expected that a change in community-level ## fall-related utilization might be detected for(per in 1 : T) { for(zas in 1 : 160) { # Poisson likelihood for observed counts f[zas,per] ~ dpois(mu[zas,per]) log(mu[zas,per]) <- log(CMS[zas,per]*0.03869) + alpha[1] ## the constant 0.03869 represents the global, six-month rate of fall-related utilization across ## ZCTAs in both study arms over the five years corresponding to the intervention and evaluation periods + alpha[2]*(BLFR[zas]-0.06938) ## BLFR is each ZCTAs baseline fall rate, measured in the two years prior to intervention ## the constant 0.06938 represents the global, twelve-month rate of fall-related utilization across ## ZCTAs in both study arms in the two year pre-intervention period. The twelve-month rate was ## used simply because it was more readily available from the pre-intervention analysis, the ## six-month rate would be about half this and could be used but this would have implied work ## and time to re-structure the dataset. + alpha[3]*(age[zas]) ## age is ordinal with values from 1 - 4 + alpha[4]*(sex[zas]-1) ## sex in dataset is 1(male) and 2( female), but subtracting 1 from each ## we can calculate risk ratio of female referent to male + alpha[5]*(treatment[zas] ) ## treatment is 0 for usual care and 1 for treatment, again to enable ## risk ratio calculation referent to usual care + alpha[6]*(Eval[per]) ## Eval is a variable that tells whether each of the ten six-month periods ## under consideration is from the intervention period (0) or the evaluation period (1) + alpha[7]*(treatment[zas] )*(Eval[per]) ## the interaction of treatment and PostSat is needed to evaluate whether the rate of fall-related ## utilization is lower in the treatment arm ZCTAs during the time of evaluationi r ##referent to the usual care ZCTAs during the intervention period. + alpha[8]*(PRlt15k[zas]) ## PRlt15k is the proportion of folks 65+ in each ZCTA with annual income < 15k from census + alpha[9]*(PRgt75k[zas]) ## PRgt75k is the proportion of folks 65+ in each ZCTA with annual income > 75k from census + alpha[10]*(PRinst[zas]) ## PRinst is the proportion of folks 65+ in each ZCTA that are institutionalized + alpha[11]*(PRphDis[zas]) ## PRinst is the proportion of folks 65+ in each ZCTA that are non-institutionalized with physical disability + alpha[12]*(PRnonWh[zas]) ## PRinst is the proportion of folks 65+ in each ZCTA that are non-white + theta[zas,per] ## theta is the random effect of undefined heterogeneity for each stratum (zas) + phiC[zas,per] ## phiC is the spatial residual for each ZCTA in the usual care arm (C indicates control) ## Next we calculate the rates of each zas so later we can combine them and estimate ## rates for individual ZCTAs enroute to calculating rates for each study arm Rate[zas,per]<-exp( alpha[1] + alpha[2]*(BLFR[zas]-0.06938) + alpha[3]*(age[zas]) + alpha[4]*(sex[zas]-1) + alpha[5]*(treatment[zas] ) + alpha[6]*(Eval[per]) + alpha[7]*(treatment[zas] )*(Eval[per]) + alpha[8]*(PRlt15k[zas]) + alpha[9]*(PRgt75k[zas]) + alpha[10]*(PRinst[zas]) + alpha[11]*(PRphDis[zas]) + alpha[12]*(PRnonWh[zas]) + theta[zas,per] + phiC[zas,per] ) fhat [zas,per]<- exp( log(CMS[zas,per]*0.03869) + Rate[zas,per]) ## the fhat allows for an analysis of predicted falls per zas and fResid of the corresponding residuals ################################### fResid[zas,per] <- f[zas,per] - fhat[zas,per] } } ################################################################# ########################################################################## ## next we define model ZCTAs 21-40 which are intervention area ZCTAs ## the reason we needed to separately create the nested structure to model ## the usual care region separately from treatment area is that prior analyses ## showed that the model would not converge is a single spatial distribution ## attempted to capture all the spatial variability in both arms. Therefore ## we created a separate CAR (conditionally auto-regressive) distribution ## for each study arm, which converged very quickly. This makes intuitive ## sense, that there may be different unmeasured spatial sources of variability ## in the two arms which are very different geographically. for(per in 1 : T) { for(zas in 161 : 320) { # Poisson likelihood for observed counts f[zas,per] ~ dpois(mu[zas,per]) log(mu[zas,per]) <- log(CMS[zas,per]*0.03869) + alpha[1] + alpha[2]*(BLFR[zas]-0.06938) + alpha[3]*(age[zas]) + alpha[4]*(sex[zas]-1) + alpha[5]*(treatment[zas] ) + alpha[6]*(Eval[per]) + alpha[7]*(treatment[zas] )*(Eval[per]) + alpha[8]*(PRlt15k[zas]) + alpha[9]*(PRgt75k[zas]) + alpha[10]*(PRinst[zas]) + alpha[11]*(PRphDis[zas]) + alpha[12]*(PRnonWh[zas]) + theta[zas,per] + phiT[zas,per] useless[zas,per] <- treatment[zas] + sex[zas] # the code above defined the analytical model, the next block merely calculates the rates based on the # estimated model coefficients for each ZAS which are then added together first within ZCTA and later # those ZCTA within each study arm to compare values from posterior distribution of model terms Rate[zas,per]<-exp( alpha[1] + alpha[2]*(BLFR[zas]-0.06938) + alpha[3]*(age[zas]) + alpha[4]*(sex[zas]-1) + alpha[5]*(treatment[zas] ) + alpha[6]*(Eval[per]) + alpha[7]*(treatment[zas] )*(Eval[per]) + alpha[8]*(PRlt15k[zas]) + alpha[9]*(PRgt75k[zas]) + alpha[10]*(PRinst[zas]) + alpha[11]*(PRphDis[zas]) + alpha[12]*(PRnonWh[zas]) + theta[zas,per] + phiT[zas,per] ) fhat [zas,per]<- exp( log(CMS[zas,per]*0.03869) + Rate[zas,per] ) ############################ fResid[zas,per] <- f[zas,per] - fhat[zas,per] } } ################################################################# ################################################# # CAR prior distribution for spatially correlated heterogeneity # This actually estimates the joint prior of the 53 different control ZCTAs for each of # of the 10 phiC variables, i.e., phiC1 ... phiC10, which correpond # to the 10 six-month periods of time under consideration # the parameters for these multivariate CAR distributions are: # nbrC is a vector of values identifying the specific ZCTAs adjacent # to each ZCTA where C is for the control arm # wtsC is a vector of values with 1 for each adjacent ZCTA # and 0 for non-adjacent ZCTAs, this defines the adjacency matrix # nnbrC is a vector of he number of adjacent neighboring ZCTAs # tau.phi is a vector of values of the random effect known # as the spatial residual at each of the 10 periods of data collection phiC1[1:53]~car.normal(nbrC[], wtsC[],nnbrC[],tau.phi[1]) phiC2[1:53]~car.normal(nbrC[], wtsC[],nnbrC[],tau.phi[2]) phiC3[1:53]~car.normal(nbrC[], wtsC[],nnbrC[],tau.phi[3]) phiC4[1:53]~car.normal(nbrC[], wtsC[],nnbrC[],tau.phi[4]) phiC5[1:53]~car.normal(nbrC[], wtsC[],nnbrC[],tau.phi[5]) phiC6[1:53]~car.normal(nbrC[], wtsC[],nnbrC[],tau.phi[6]) phiC7[1:53]~car.normal(nbrC[], wtsC[],nnbrC[],tau.phi[7]) phiC8[1:53]~car.normal(nbrC[], wtsC[],nnbrC[],tau.phi[8]) phiC9[1:53]~car.normal(nbrC[], wtsC[],nnbrC[],tau.phi[9]) phiC10[1:53]~car.normal(nbrC[], wtsC[],nnbrC[],tau.phi[10]) ################################################### # CAR prior distribution for spatial correlated heterogeneity # This actually estimates the joint prior of the 58 different treatment ZCTAs for each of # of the 10 phiT variables, i.e., phiT1 ... phiT10, correponding # to the 10 six-month periods of time under consideration # note that we define different spatial distributions for the # T (treatment) and control arms of the study phiT1[1:58]~car.normal(nbrT[], wtsT[],nnbrT[],tau.phi[1]) phiT2[1:58]~car.normal(nbrT[], wtsT[],nnbrT[],tau.phi[2]) phiT3[1:58]~car.normal(nbrT[], wtsT[],nnbrT[],tau.phi[3]) phiT4[1:58]~car.normal(nbrT[], wtsT[],nnbrT[],tau.phi[4]) phiT5[1:58]~car.normal(nbrT[], wtsT[],nnbrT[],tau.phi[5]) phiT6[1:58]~car.normal(nbrT[], wtsT[],nnbrT[],tau.phi[6]) phiT7[1:58]~car.normal(nbrT[], wtsT[],nnbrT[],tau.phi[7]) phiT8[1:58]~car.normal(nbrT[], wtsT[],nnbrT[],tau.phi[8]) phiT9[1:58]~car.normal(nbrT[], wtsT[],nnbrT[],tau.phi[9]) phiT10[1:58]~car.normal(nbrT[], wtsT[],nnbrT[],tau.phi[10]) ############################################ ######### ZCTA Rates in Usual Care Arm ! # Calculate rates for 20 of the control ZCTAs by adding up the rates of each one's eight age-sex strata (ZAS) # in each of the 10 six-month periods under consideration # Each of the 20 control ZCTAs is assigned to its own set of 8 age-sex strata and the zas are respectively numbered 1 - 160 for (zcta in 1:20) { RateZCTA[zcta,1]<- (Rate[(zcta-1)*8 + 1,1]+Rate[(zcta-1)*8 +2 ,1]+Rate[(zcta-1)*8 +3 ,1]+Rate[(zcta-1)*8 + 4,1] +Rate[(zcta-1)*8 + 5,1]+Rate[(zcta-1)*8 + 6,1]+Rate[(zcta-1)*8 + 7,1]+Rate[(zcta-1)*8 + 8,1]) RateZCTA[zcta,2]<- (Rate[(zcta-1)*8 + 1,2]+Rate[(zcta-1)*8 +2 ,2]+Rate[(zcta-1)*8 +3 ,2]+Rate[(zcta-1)*8 + 4,2] +Rate[(zcta-1)*8 + 5,2]+Rate[(zcta-1)*8 + 6,2]+Rate[(zcta-1)*8 + 7,2]+Rate[(zcta-1)*8 + 8,2]) RateZCTA[zcta,3]<- (Rate[(zcta-1)*8 + 1,3]+Rate[(zcta-1)*8 +2 ,3]+Rate[(zcta-1)*8 +3 ,3]+Rate[(zcta-1)*8 + 4,3] +Rate[(zcta-1)*8 + 5,3]+Rate[(zcta-1)*8 + 6,3]+Rate[(zcta-1)*8 + 7,3]+Rate[(zcta-1)*8 + 8,3]) RateZCTA[zcta,4]<- (Rate[(zcta-1)*8 + 1,4]+Rate[(zcta-1)*8 +2 ,4]+Rate[(zcta-1)*8 +3 ,4]+Rate[(zcta-1)*8 + 4,4] +Rate[(zcta-1)*8 + 5,4]+Rate[(zcta-1)*8 + 6,4]+Rate[(zcta-1)*8 + 7,4]+Rate[(zcta-1)*8 + 8,4]) RateZCTA[zcta,5]<- (Rate[(zcta-1)*8 + 1,5]+Rate[(zcta-1)*8 +2 ,5]+Rate[(zcta-1)*8 +3 ,5]+Rate[(zcta-1)*8 + 4,5] +Rate[(zcta-1)*8 + 5,5]+Rate[(zcta-1)*8 + 6,5]+Rate[(zcta-1)*8 + 7,5]+Rate[(zcta-1)*8 + 8,5]) RateZCTA[zcta,6]<- (Rate[(zcta-1)*8 + 1,6]+Rate[(zcta-1)*8 +2 ,6]+Rate[(zcta-1)*8 +3 ,6]+Rate[(zcta-1)*8 + 4,6] +Rate[(zcta-1)*8 + 5,6]+Rate[(zcta-1)*8 + 6,6]+Rate[(zcta-1)*8 + 7,6]+Rate[(zcta-1)*8 + 8,6]) RateZCTA[zcta,7]<- (Rate[(zcta-1)*8 + 1,7]+Rate[(zcta-1)*8 +2 ,7]+Rate[(zcta-1)*8 +3 ,7]+Rate[(zcta-1)*8 + 4,7] +Rate[(zcta-1)*8 + 5,7]+Rate[(zcta-1)*8 + 6,7]+Rate[(zcta-1)*8 + 7,7]+Rate[(zcta-1)*8 + 8,7]) RateZCTA[zcta,8]<- (Rate[(zcta-1)*8 + 1,8]+Rate[(zcta-1)*8 +2 ,8]+Rate[(zcta-1)*8 +3 ,8]+Rate[(zcta-1)*8 + 4,8] +Rate[(zcta-1)*8 + 5,8]+Rate[(zcta-1)*8 + 6,8]+Rate[(zcta-1)*8 + 7,8]+Rate[(zcta-1)*8 + 8,8]) RateZCTA[zcta,9]<- (Rate[(zcta-1)*8 + 1,9]+Rate[(zcta-1)*8 +2 ,9]+Rate[(zcta-1)*8 +3 ,9]+Rate[(zcta-1)*8 + 4,9] +Rate[(zcta-1)*8 + 5,9]+Rate[(zcta-1)*8 + 6,9]+Rate[(zcta-1)*8 + 7,9]+Rate[(zcta-1)*8 + 8,9]) RateZCTA[zcta,10]<- (Rate[(zcta-1)*8 + 1,10]+Rate[(zcta-1)*8 +2 ,10]+Rate[(zcta-1)*8 +3 ,10]+Rate[(zcta-1)*8 + 4,10] +Rate[(zcta-1)*8 + 5,10]+Rate[(zcta-1)*8 + 6,10]+Rate[(zcta-1)*8 + 7,10]+Rate[(zcta-1)*8 + 8,10]) } ########################################################### ############# ZCTA Rates in Intervention Arm ! # Calculate rates for 20 of the treatment ZCTAs by adding up the rates of each one's eight age-sex strata (ZAS) # in each of the 10 six-month periods under consideration # The zas are respectively numbered 161 - 320 based on groups of 8 each for ZCTAs 21 - 40 for (zcta in 21:40) { RateZCTA[zcta,1]<- (Rate[(zcta-1)*8 +1,1]+Rate[(zcta-1)*8 +2 ,1]+Rate[(zcta-1)*8 +3 ,1]+Rate[(zcta-1)*8 +4,1] +Rate[(zcta-1)*8 +5,1]+Rate[(zcta-1)*8 +6,1]+Rate[(zcta-1)*8 +7,1]+Rate[(zcta-1)*8 +8,1]) RateZCTA[zcta,2]<- (Rate[(zcta-1)*8 +1,2]+Rate[(zcta-1)*8 +2 ,2]+Rate[(zcta-1)*8 +3 ,2]+Rate[(zcta-1)*8 +4,2] +Rate[(zcta-1)*8 +5,2]+Rate[(zcta-1)*8 +6,2]+Rate[(zcta-1)*8 +7,2]+Rate[(zcta-1)*8 +8,2]) RateZCTA[zcta,3]<- (Rate[(zcta-1)*8 +1,3]+Rate[(zcta-1)*8 +2 ,3]+Rate[(zcta-1)*8 +3 ,3]+Rate[(zcta-1)*8 +4,3] +Rate[(zcta-1)*8 +5,3]+Rate[(zcta-1)*8 +6,3]+Rate[(zcta-1)*8 +7,3]+Rate[(zcta-1)*8 +8,3]) RateZCTA[zcta,4]<- (Rate[(zcta-1)*8 +1,4]+Rate[(zcta-1)*8 +2 ,4]+Rate[(zcta-1)*8 +3 ,4]+Rate[(zcta-1)*8 +4,4] +Rate[(zcta-1)*8 +5,4]+Rate[(zcta-1)*8 +6,4]+Rate[(zcta-1)*8 +7,4]+Rate[(zcta-1)*8 +8,4]) RateZCTA[zcta,5]<- (Rate[(zcta-1)*8 +1,5]+Rate[(zcta-1)*8 +2 ,5]+Rate[(zcta-1)*8 +3 ,5]+Rate[(zcta-1)*8 +4,5] +Rate[(zcta-1)*8 +5,5]+Rate[(zcta-1)*8 +6,5]+Rate[(zcta-1)*8 +7,5]+Rate[(zcta-1)*8 +8,5]) RateZCTA[zcta,6]<- (Rate[(zcta-1)*8 +1,6]+Rate[(zcta-1)*8 +2 ,6]+Rate[(zcta-1)*8 +3 ,6]+Rate[(zcta-1)*8 +4,6] +Rate[(zcta-1)*8 +5,6]+Rate[(zcta-1)*8 +6,6]+Rate[(zcta-1)*8 +7,6]+Rate[(zcta-1)*8 +8,6]) RateZCTA[zcta,7]<- (Rate[(zcta-1)*8 +1,7]+Rate[(zcta-1)*8 +2 ,7]+Rate[(zcta-1)*8 +3 ,7]+Rate[(zcta-1)*8 +4,7] +Rate[(zcta-1)*8 +5,7]+Rate[(zcta-1)*8 +6,7]+Rate[(zcta-1)*8 +7,7]+Rate[(zcta-1)*8 +8,7]) RateZCTA[zcta,8]<- (Rate[(zcta-1)*8 +1,8]+Rate[(zcta-1)*8 +2 ,8]+Rate[(zcta-1)*8 +3 ,8]+Rate[(zcta-1)*8 +4,8] +Rate[(zcta-1)*8 +5,8]+Rate[(zcta-1)*8 +6,8]+Rate[(zcta-1)*8 +7,8]+Rate[(zcta-1)*8 +8,8]) RateZCTA[zcta,9]<- (Rate[(zcta-1)*8 +1,9]+Rate[(zcta-1)*8 +2 ,9]+Rate[(zcta-1)*8 +3 ,9]+Rate[(zcta-1)*8 +4,9] +Rate[(zcta-1)*8 +5,9]+Rate[(zcta-1)*8 +6,9]+Rate[(zcta-1)*8 +7,9]+Rate[(zcta-1)*8 +8,9]) RateZCTA[zcta,10]<- (Rate[(zcta-1)*8 +1,10]+Rate[(zcta-1)*8 +2 ,10]+Rate[(zcta-1)*8 +3 ,10]+Rate[(zcta-1)*8 +4,10]+Rate[(zcta-1)*8 +5,10]+Rate[(zcta-1)*8 +6,10]+Rate[(zcta-1)*8 +7,10]+Rate[(zcta-1)*8 +8,10]) } ########################### Rates for Usual Care Arm################################ for (per in 1:T) { RZ1to8[per ] <- RateZCTA[1,per]+RateZCTA[2,per]+RateZCTA[3,per]+RateZCTA[4,per]+RateZCTA[5,per]+RateZCTA[6,per]+RateZCTA[7,per]+RateZCTA[8,per] RZ9to16[per] <- RateZCTA[9,per]+RateZCTA[10,per]+RateZCTA[11,per]+RateZCTA[12,per]+ RateZCTA[13,per]+RateZCTA[14,per]+RateZCTA[15,per]+RateZCTA[16,per] RZ17to20[per] <- RateZCTA[17,per]+RateZCTA[18,per]+RateZCTA[19,per]+RateZCTA[20,per] UCrate[per] <- (RZ1to8[per] + RZ9to16[per] + RZ17to20[per]) / 20 } UCrateInterv <- (UCrate[1]+UCrate[2]+UCrate[3]+UCrate[4]+UCrate[5]+UCrate[6])/6 UCrateEval <- (UCrate[7]+UCrate[8]+UCrate[9]+UCrate[10])/4 ############################ Rates for Intervention Arm ################################## for (per in 1:T) { RZ21to28[per ] <- RateZCTA[21,per]+RateZCTA[22,per]+RateZCTA[23,per]+RateZCTA[24,per]+RateZCTA[25,per]+RateZCTA[26,per]+RateZCTA[27,per]+RateZCTA[28,per] RZ29to36[per] <- RateZCTA[29,per]+RateZCTA[30,per]+RateZCTA[31,per]+RateZCTA[32,per]+ RateZCTA[33,per]+RateZCTA[34,per]+RateZCTA[35,per]+RateZCTA[36,per] RZ37to40[per] <- RateZCTA[37,per]+RateZCTA[38,per]+RateZCTA[39,per]+RateZCTA[40,per] IVrate[per] <- (RZ21to28[per ]+ RZ29to36[per] + RZ37to40[per] )/20 } IVrateInterv <- (IVrate[1]+IVrate[2]+IVrate[3]+IVrate[4]+IVrate[5]+IVrate[6])/6 IVrateEval <- (IVrate[7]+IVrate[8]+IVrate[9]+IVrate[10])/4 ########################### Credibility Interval for Rate Ratio in Evaluation Period ######################################## RateRatioEval <- (UCrateEval)/(IVrateEval) # Assign zctas for zas combinations 1-160 for tracking of the # spatial residual terms in the control ZCTAs in the periods 1 - 10 for (zcta in 1:20) { for (sa in 1:8) { phiC[((zcta-1)*8 + sa),1]<-phiC1[zcta] phiC[((zcta-1)*8 + sa),2]<-phiC2[zcta] phiC[((zcta-1)*8 + sa),3]<-phiC3[zcta] phiC[((zcta-1)*8 + sa),4]<-phiC4[zcta] phiC[((zcta-1)*8 + sa),5]<-phiC5[zcta] phiC[((zcta-1)*8 + sa),6]<-phiC6[zcta] phiC[((zcta-1)*8 + sa),7]<-phiC7[zcta] phiC[((zcta-1)*8 + sa),8]<-phiC8[zcta] phiC[((zcta-1)*8 + sa),9]<-phiC9[zcta] phiC[((zcta-1)*8 + sa),10]<-phiC10[zcta] } } ########################################################### # Assign zctas for zas combinations 161-320 for tracking of the # spatial residual terms in the treatment ZCTAs in the periods 1 - 10 # we add 160 so that the resultant ZAS go from 161 - 320 for treatment for (zcta in 1:20) { for (sa in 1:8) { phiT[((zcta-1)*8 + sa + 160),1]<-phiT1[zcta] phiT[((zcta-1)*8 + sa + 160),2]<-phiT2[zcta] phiT[((zcta-1)*8 + sa + 160),3]<-phiT3[zcta] phiT[((zcta-1)*8 + sa + 160),4]<-phiT4[zcta] phiT[((zcta-1)*8 + sa + 160),5]<-phiT5[zcta] phiT[((zcta-1)*8 + sa + 160),6]<-phiT6[zcta] phiT[((zcta-1)*8 + sa + 160),7]<-phiT7[zcta] phiT[((zcta-1)*8 + sa + 160),8]<-phiT8[zcta] phiT[((zcta-1)*8 + sa + 160),9]<-phiT9[zcta] phiT[((zcta-1)*8 + sa + 160),10]<-phiT10[zcta] } } ########################################################### # Prior distributions for the uncorrelated heterogeneity # Because we calculate this term in each period, we # set up 10 vaguely informative priors for the ten unique periods for (zasT in 1:320) { theta1[zasT]~ dnorm(0, tau.theta[1]) theta2[zasT]~ dnorm(0, tau.theta[2]) theta3[zasT]~ dnorm(0, tau.theta[3]) theta4[zasT]~ dnorm(0, tau.theta[4]) theta5[zasT]~ dnorm(0, tau.theta[5]) theta6[zasT]~ dnorm(0, tau.theta[6]) theta7[zasT]~ dnorm(0, tau.theta[7]) theta8[zasT]~ dnorm(0, tau.theta[8]) theta9[zasT]~ dnorm(0, tau.theta[9]) theta10[zasT]~ dnorm(0, tau.theta[10]) } for (zasT in 1:320) { theta[zasT,1]<-theta1[zasT] theta[zasT,2]<-theta2[zasT] theta[zasT,3]<-theta3[zasT] theta[zasT,4]<-theta4[zasT] theta[zasT,5]<-theta5[zasT] theta[zasT,6]<-theta6[zasT] theta[zasT,7]<-theta7[zasT] theta[zasT,8]<-theta8[zasT] theta[zasT,9]<-theta9[zasT] theta[zasT,10]<-theta10[zasT] } ################################################# # Because our analysis compared the rates in the evaluation period to those in the # intervention period, we number the six periods from intervention as 0 and the # four from evaluation as 1 for (per in 1:6) { Eval[per] <-0 } for (per in 7:T) { Eval[per] <-1 } ################################################################## # Weights for(k in 1:280) { wtsT[k]<-1 } for(k in 1:206) { wtsC[k]<-1 } #################################################################### # Improper prior distribution for the mean relative risk in the study region for (i in 1:8) { alpha[i] ~ dflat() } # Because the coefficients for the proportions describing the ZCTAs from census # data are very noisy, we use a vaguely informative prior regarding their variances # to limit the range of possible values, this facilitates convergence for (i in 9:12) { alpha[i] ~ dnorm(0,39) } #################################################################### #Hyperprior distributions on inverse variance parameter of random effects for (per in 1:10) { tau.theta[per] ~ dgamma(1,1) tau.phi[per] ~ dgamma(1,1) } } ###################################################################### ####################################################################