##############################################################################
## County-Specific Estimates

## Estimate county-specific RR and variances
## An .Rdata file is saved that contains two lists (county RR and county vars
## Indexed by a matrix "info" that gives the location for each outcome/lag/model

library (splines)
library (tsModel)
library (tlnise)

#Set working directory
setwd ("~/CoarsePM/For Roger - paper")
#Load model functions
source("~/CoarsePM/Codes/Utilities.R")

### Read in county fips 
### 108 counties with collocated coarse PM > 210 days
### 204 counties for PM2.5 estimates ONLY
fips.108 <- read.csv ("fips.108.csv", colClasses=c("character"))[,1]
fips.204 <- read.csv("fips.204.csv", colClasses=c("character"))[,1]

#Outcomes: Aggregated cardiovascular and respiratory emergency admissions
outcomes <- c("CVD", "Resp")

#Df for time (8 PER YEAR)
d <- 8*7

#Stratify By Age?
age75p <- FALSE
age65to74 <- FALSE

#Set up lists to record results
n <- (108*2 + (204-108)*3)*length(outcomes)*3
betas <- vector("list", n)
vars <- vector("list", n)
info <- matrix (NA, n,4) #To extract estimates from lists
i <- 0
err <- NULL

# Loop over county

for (siteName in fips.204){
        tmp.data <- .readRDS (paste("~/CoarsePM/Data/MCAPS1999-2005/", siteName, ".rds", sep = ""))
        names (tmp.data)[123] <- "pmCoarse"
        
        admitCVD <-  tmp.data$admitp3 + tmp.data$admitp4 + tmp.data$admitp5 + tmp.data$admitp6 + tmp.data$admitp7 
        admitResp <- tmp.data$admitp8 + tmp.data$admitp9
           
        #Extract variables from RDS file       
        admissions <- data.frame (tmp.data$date, tmp.data$agecat, tmp.data$denom, tmp.data$tmpd, tmp.data$rmtmpd,
                                  tmp.data$dow, tmp.data$dptp, tmp.data$rmdptp,
                                  tmp.data$pmCoarse, tmp.data$pm10c, tmp.data$pm25c,
                                  admitCVD, admitResp, tmp.data$admitp11, tmp.data$pm25tmean2)
        names (admissions) <- c("date", "agecat", "denom", "tmpd", "rmtmpd", "dow", "dptp", "rmdptp",
                                "pmCoarse", "pm10c", "pm25c", "admitCVD", "admitResp", "admitp11", "pm25tmean2")

        
        if (siteName %in% fips.108){ #What models to fie#
            models <- c("model.pm25", "model.coarse", "model.pm25.coarse")} else {
            models <- c("model.pm25") }
        
        for (model in models){  #Loop over models#
          
          p2 <- NULL
          if (model == "model.pm25") { p1 <- "pm25tmean2"}  #County  PM2.5
          if (model == "model.coarse") {p1 <- "pmCoarse"} #Collocated PM2.5
          if (model == "model.pm25.coarse") {p1 <- "pmCoarse"; p2 <- "pm25c"} #Bi-pollutant model
          
          polls <- structure(c(paste( "Lag(", p1, ", 0, agecat)", sep = ""), 
                               paste( "Lag(", p1, ", 1, agecat)", sep = ""), 
                               paste( "Lag(", p1, ", 2, agecat)", sep = "")
                               ), names = c(paste("Lag", 0:2))
                             )
                    
          polls2 <- NULL
          if (!is.null (p2)){
            polls2 <- structure(c(
                                  paste( "Lag(", p2, ", 0, agecat)", sep = ""  ), 
                                  paste( "Lag(", p2, ", 1, agecat)", sep = ""  ), 
                                  paste( "Lag(", p2, ", 2, agecat)", sep = ""  )
                                  ), names = c(paste("Lag", 0:2))
                                )
          }
                    
        for(outcome in outcomes) {  #Cycle over outcome#
          
          for(poll in polls) {   #Cycle of lag
            cat("Emergency", model, outcome, poll, siteName, d/7, "\n")
            which.lag <- which (polls == poll)
            poll2 <- polls2[which.lag]
            
            rval <- tryCatch({f <- fitSingleSite(admissions, outcome = outcome,
                                                 pollutant = poll, df.Time = d,
                                                 subset = T,poll2 = poll2 )
                              junk <- postProcess(f, "pm")
                              list (estimates = c(junk$beta), var = as.matrix(junk$var), info=c(model, siteName, outcome, which.lag -1) )
                              
                            },  condition = function(cond) {
                            cat("\t", trunc(d / 4), as.character(cond))
                            cond
                          })
                if (is.null(rval$message)){
                  i <- i + 1
                  betas[[i]] <- rval$estimates
                  vars[[i]] <- rval$var
                  info[i,] <- rval$info
                } else{
                  err <- c(err, paste(model, siteName, outcome, poll, as.character(d/7),":", rval$message, sep = " ") )
                }
    
          }  ## End cycle over lags
        }  ## End cycle over outcomes
        } ## End cycle over models
      } ## End cycle over counties


info <- data.frame (info)
names (info) <- c("model", "fips", "outcome", "lag")

betas <- betas[complete.cases (info)]
vars <- vars[complete.cases (info)]
info <- info[complete.cases (info),]

#Assign geographic locations and county names to info
#Excluding Anchorage AK and Honolulu HI
info$regionEW <- NA
info$CountyName <- NA

Countyinfo <- read.csv ("~/CoarsePM/Data/CountyInfo.csv", colClasses = c("character", "numeric",rep("character", 4) ))

for (siteName in fips.204 ) {
    info$regionEW [ info$fips == siteName ] <- as.character(Countyinfo$regionEW[Countyinfo$fips == siteName][1])
    info$CountyName [ info$fips == siteName ] <- as.character(Countyinfo$CountyName[Countyinfo$fips == siteName][1])
    }

save(info, vars, betas, file = "Results.Emergency.Rdata" )