Sys.setlocale(locale = "C")
estimates <- read.csv("http://www.biostat.jhsph.edu/MCAPS/estimates-subset.csv")

## Reformat FIPS codes
estimates <- transform(estimates, fips = formatC(fips, flag = "0", width =5))

library(tlnise)

## Get national risk estimates for heart failure
HF <- subset(estimates, outcome == "heart failure")

initTLNise()
pooled <- with(HF, tlnise(beta, var, prnt = FALSE, seed = 1234))

print(pooled$gamma)

## Get national risk estimates for COPD
COPD <- subset(estimates, outcome == "COPD")

initTLNise()
pooled <- with(COPD, tlnise(beta, var, prnt = FALSE, seed = 4321))

print(pooled$gamma)

## Compare repiratory infection risk estimates in East and West
east <- subset(estimates, outcome == "respiratory infection"
               & regionEW == "East")
west <- subset(estimates, outcome == "respiratory infection"
               & regionEW == "West")

initTLNise()
pooled.east <- with(east, tlnise(beta, var, prnt = FALSE, seed = 1))

initTLNise()
pooled.west <- with(west, tlnise(beta, var, prnt = FALSE, seed = 2))

rbind(pooled.east$gamma, pooled.west$gamma)