Embedded Files
On this page, I will walk you through the steps to recreate the results in my "Estimating county-level mortality rates using highly censored data from CDC WONDER" paper. This walkthrough assumes you already have R and WinBUGS installed on your computer. For instructions on how to install WinBUGS, click here. For an example of how to use WinBUGS, click here.
Step 1: Follow Steps 1-5(1) on the WONDER page.
Step 2: Open WinBUGS, and open the following four files in "WONDER/Code/BUGS"
BUGS_MCAR_model.txt - The WinBUGS model code.
HD_data.txt - The list of HD data
shp_data.txt - The adjacency structure required for the MCAR model. This should be loaded in as data.
initial_values.txt - A list of initial values for select model parameters.
Note that by default, WinBUGS uses the ".odc" file extension, so you will need to change that to "Text (*.txt)" or "All files (*.*)" to find these files in the folder. Alternatively, you can open these files outside of WinBUGS and simply copy/paste them into new windows within WinBUGS.
Step 3: Check the model, load in both lists of data, compile the code, load in the list of initial values, and then generate initial values for the remaining model parameters. There should be a message in the bottom-left corner that says "model is initialized" -- if so, you're ready to go.
Step 4: We're now ready to fit the model. First, I recommend "set"-ing "Y" and "beta0" -- this will allow you to see that you are indeed fitting a truncated Poisson model (i.e., the Y's will all be less than 10, as shown below) and [visually] assess convergence. I then recommend running this code first for 1,000 iterations, then for another 4,000. This will allow you to gauge how long the code will ultimately take to run, as well as see if there are any serious issues (e.g., horrible convergence for beta0). Once you're convinced that the model is working well, I recommend "set"-ing "lambda" and running another 5,000 iterations.
Step 5: Once you've obtained your 5,000 iterations' worth of samples for lambda, we need to export these samples out of WinBUGS. To do this, we will use the "coda" button the Sample Monitor Tool. This will (eventually) open two windows in WinBUGS -- one file contains the 5,000 posterior samples for each of the 3,109 x 6 lambda parameters (i.e., a *ton* of rows), while the other file simply tells you which rows of the first file correspond to which parameter (e.g., rows 1-5000 are lambda[1,1], rows 5001:10000 are lambda[1,2], etc.). My code assumes you save these two files as "coda.txt" and "coda_index.txt" in the "WONDER" directory.
Step 6: We now return to R and open "WONDER/Code/2_results_bugs.R". Some highlights of this file are:
The following code reads in the "coda.txt" file we exported from WinBUGS and formats it as a 3-dimensional array.
bugs=read.table('coda.txt',sep='\t')
nsims=dim(bugs)[1]/(Ns*Ng)
lami=array(bugs[,2],dim=c(nsims,Ng,Ns))
If you have also explored the output from fitting this model (or the Poisson-gamma model) in R, you should note that the order of these dimensions
(sample x age x county)
is the opposite of the order used previously
(county x age x sample)
The next piece of code creates the maps of the age-adjusted rates. As I might have described elsewhere, I use a "paint-by-numbers" approach to coloring my maps -- if you are familiar with an alternative approach, by all means ignore this and do whatever works for you :)
png('Figures/Maps/aa_bugs.png', #<- where the map is being saved
height=2*2*650,width=2*2*1000,res=600)
par(mar=c(0,0,0,0)+.1,cex=1)
tcuts=quantile(aalci[1,]*100000,1:(ncols-1)/ncols)
tcolb=array(rep(aalci[1,]*100000,each=ncols-1) > tcuts,
dim=c(ncols-1,Ns))
tcol =apply(tcolb,2,sum)+1
plot(us,col=cols[tcol],border='lightgray',lwd=.5)
plot(states,border='black',add=TRUE,lwd=.5)
legend('bottomright',legend=c(paste(
c('Below',round(tcuts[-(ncols-1)],0),'Over'),
c(' ',rep( ' - ',ncols-2),' '),
c(round(tcuts,0),round(tcuts[ncols-1],0)),sep='')),
fill=cols,title='Deaths per 100,000',bty='n',cex=1.5/2,
border='lightgray')
dev.off()
In particular, "tcuts" -- the cut-points for my coloring -- is defined using quantiles of the distribution of age-adjusted rates, and then I construct "tcol" by calculating the number of cut-points each county's rate is greater than (+1, so the lowest rates are shaded in color #1). I then use the "plot(us,...)" code to map the counties with the appropriate colors and overlay a map of the state boundaries.
Page updated
Google Sites
Report abuse