####################################################################################
## County-Specific Estimates: 2nd Stage Analysis

## First load the .Rdata obtained from "CountyEstimate.R" 
## It contains the county-specific RRs, "betas", and variances, "vars", in two lists
## indexed by the matrix "info"

load ("Results.Emergency.Rdata")
urbanicity <- read.csv ("CountyInfo.csv", colClasses = c("character","numeric",rep("character",4),"numeric"))[, c("fips", "Purban")]

info <- merge (info, urbanicity, "fips", "fips") 

models <- c("model.pm25", "model.coarse", "model.pm25.coarse")
lags <- c(0, 1, 2)
outcomes <- c("CVD", "Resp")

result.final.T <- NULL #Record results
library (tlnise)

for (model in models){
  for (outcome in outcomes){
    for (lag in lags){
      print (paste (model, outcome,lag))

      #Single Pollutant Model 
      if (model == "model.pm25" | model == "model.coarse"){
        use <- info$model == model & info$lag == lag & info$outcome == outcome
        param.beta = matrix(unlist (betas[use]), nrow=sum(use), ncol=1, byrow =T) 
        param.var = array(unlist (vars[use]), c(1, 1, sum(use)))
        param.urb <- info[use, ]$Purban*100
        
        fit <- tlnise (param.beta, param.var, cbind(rep(1,sum (use)), param.urb), prnt = F)
        
        CI.1 <- round( (fit$gamma[2,1]-1.96*fit$gamma[2,2])*1000, 3)
        CI.2 <- round( (fit$gamma[2,1]+1.96*fit$gamma[2,2])*1000, 3)
        Est <- round(  (fit$gamma[2,1])*1000,3)
      
        result.final.T <- rbind (result.final.T, c(model, outcome, lag, Est, CI.1, CI.2))
      }

      #Bi-pollutant model
      if (model == "model.pm25.coarse"){
        use <- info$model == "model.pm25.coarse"
        use <- info$model == model & info$lag == lag & info$outcome == outcome
        param.beta = matrix(unlist (betas[use]), nrow=sum(use), ncol=2, byrow =T) 
        param.var = array(unlist (vars[use]), c(2, 2, sum(use)))
        param.urb <- info[use, ]$Purban*100
       
        fit <- tlnise (param.beta, param.var, cbind(rep(1,sum (use)), param.urb), prnt = F)
        CI.1 <- round( (fit$gamma[3:4,1]-1.96*fit$gamma[3:4,2])*1000, 3)
        CI.2 <- round( (fit$gamma[3:4,1]+1.96*fit$gamma[3:4,2])*1000, 3)
        Est <- round(  (fit$gamma[3:4,1])*1000,3)

        result.final.T <- rbind (result.final.T, c("pm25.joint", outcome, lag, Est[1], CI.1[1], CI.2[1]))
        result.final.T <- rbind (result.final.T, c("coarse.joint", outcome, lag, Est[2], CI.1[2], CI.2[2]))
      }
    } #End for Lag   
  } #End for outcome
} #End for model

result.final.T <- data.frame (result.final.T)
names (result.final.T) <- c("Model", "Outcome", "Lag", "Slope", "Lower", "Upper")

write.csv (result.final.T, file = "2level.ubanicity.csv", row.names = F)