##############################################
##Regression Calibration Function

##Input takes the county's FIPS
##Return a dataframe to be analyzed similarly as in
##"CountyEstimate.R"

library (splines)
Reg.Cal <- function (siteName){

    load (paste ("~/CoarsePM/Data/Monitor.Specific/",siteName, ".monitor.Rdata", sep = ""))
    
    admissions.junk <- .readRDS (paste("~/CoarsePM/Data/MCAPS1999-2005/", siteName, ".rds", sep = ""))
    names (admissions.junk)[123] <- "pmCoarse"

    admissions <- data.frame (admissions.junk$date, admissions.junk$tmpd, admissions.junk$rmtmpd,admissions.junk$dow, 
                                admissions.junk$dptp, admissions.junk$rmdptp, admissions.junk$agecat)
    names (admissions) <- c("date", "tmpd", "rmtmpd", "dow", "dptp", "rmdptp", "agecat")
   
    admissions <- subset (admissions, agecat == "75p")

    dow1 <- ifelse (admissions$dow == "Monday", 1, 0)
    dow2 <- ifelse (admissions$dow == "Tuesday",1, 0)
    dow3 <- ifelse (admissions$dow == "Wednesday", 1, 0)
    dow4 <- ifelse (admissions$dow == "Thursday", 1, 0)
    dow5 <- ifelse (admissions$dow == "Friday", 1, 0)
    dow6 <- ifelse (admissions$dow == "Saturday", 1, 0)
    
    Z <- cbind (dow1, dow2, dow3, dow4, dow5, dow6
           , ns (admissions$date, 3*7)
            , ns(admissions$tmpd, 6), ns(admissions$dptp, 3)
            ,ns(admissions$rmtmpd, 6),  ns(admissions$rmdptp,3))

    #Number of collocated monitors
    num.site <-  ncol (monitor.data[[3]])-1
    
    coll.id <- names (monitor.data[[3]])
    W.pm25 <- monitor.data[[2]][, coll.id]
    W.coarse <- monitor.data[[3]]
    W.pm10 <- monitor.data[[1]][, coll.id]
    
    #Remove Missing Values
    id.use <- apply (!is.na(W.coarse[,-1]), 1, sum, na.rm = T)!=0 & complete.cases (Z)
    W.pm25 <- W.pm25[id.use,]
    W.coarse <- W.coarse[id.use,]
    W.pm10 <- W.pm10[id.use,]
    Z <- Z[id.use,] 

    #Remove redundent columns in Z
    id <- !is.na(lm (rep(1, nrow (Z))~Z-1)$coef)
    Z <- Z[, id]
                  
    #To compute SigmaUU
    bar.pm25 <- apply (W.pm25[,-1], 1, mean, na.rm = T)
    bar.coarse <- apply (W.coarse[,-1], 1, mean, na.rm = T)
    bar.pm10 <- apply (W.pm10[,-1], 1, mean, na.rm = T)
    resid.pm25 <- as.matrix(W.pm25[,-1] - bar.pm25)
    resid.coarse <- as.matrix(W.coarse[,-1] - bar.coarse)
    resid.pm10 <- as.matrix (W.pm10[,-1] - bar.pm10)
    
    ki.coarse <- apply (!is.na(W.coarse[,-1]), 1, sum, na.rm = T)
    ki.pm25 <- apply (!is.na(W.pm25[,-1]), 1, sum, na.rm = T)
    ki.pm10 <- apply (!is.na(W.pm10[,-1]), 1, sum, na.rm = T)
    
    K.coarse <- sum (ki.coarse -1)
    K.pm25 <- sum (ki.pm25 -1)
    K.pm10 <- sum (ki.pm10 -1)
    
    u.pm25 <-  sum(resid.pm25^2, na.rm = T)/K.pm25
    u.coarse <- sum(resid.coarse^2, na.rm = T)/K.coarse
    u.pm10 <- sum (resid.pm10^2, na.rm= T)/K.pm10
    
    v.pm25 <- sum (ki.pm25) - sum (ki.pm25^2)/sum(ki.pm25)
    v.coarse <- sum (ki.coarse) - sum (ki.coarse^2)/sum(ki.coarse)
    v.pm10 <- sum (ki.pm10) - sum (ki.pm10^2)/sum (ki.pm10)
    
    n <- nrow (W.pm25)
    
    Mu.x.pm25 <- sum (ki.pm25*bar.pm25, na.rm = T)/sum (ki.pm25)
    Mu.x.coarse <- sum (ki.coarse*bar.coarse, na.rm = T)/sum(ki.coarse) 
    Mu.x.pm10 <- sum (ki.pm10*bar.pm10, na.rm = T)/sum (ki.pm10)
    Mu.w.pm25 <- Mu.x.pm25
    Mu.w.coarse <- Mu.x.coarse
    Mu.w.pm10 <- Mu.x.pm10
    
    #To compute SigmaXX
    xx.pm25 <- ( sum( (bar.pm25 - Mu.w.pm25)^2*ki.pm25, na.rm = T ) - (n-1)*u.pm25 )/v.pm25
    xx.coarse <- ( sum( (bar.coarse - Mu.w.coarse)^2*ki.coarse, na.rm = T ) - (n-1)*u.coarse )/v.coarse
    xx.pm10 <- ( sum( (bar.pm10 - Mu.w.pm10)^2*ki.pm10, na.rm = T ) - (n-1)*u.pm10 )/v.pm10
   
    #To compute SigmaXZ
    bar.Z <- apply (Z, 2, mean, na.rm = T)
    resid.Z <- Z - matrix( bar.Z, ncol = ncol (Z), nrow = nrow (Z), byrow = T)
    zz <- (t(resid.Z) %*% resid.Z)/(n-1)

    
    xz.pm25 <- t(ki.pm25*(bar.pm25 - Mu.w.pm25)) %*% resid.Z / v.pm25
    xz.pm10 <- t(ki.pm10*(bar.pm10 - Mu.w.pm10)) %*% resid.Z / v.pm10
    xz.coarse <- t(ki.coarse*(bar.coarse-Mu.w.coarse)) %*% resid.Z / v.coarse
    
    zx.pm25 <- t(xz.pm25)
    zx.pm10 <- t(xz.pm10)
    zx.coarse <- t(xz.coarse)
     
    #Calibration
    expect.pm25 <- NULL
    expect.coarse <- NULL
    expect.pm10 <- NULL
    for (i in 1:length (bar.coarse)){

      #Calibrate PM2.5
      project.matrix <- cbind (xx.pm25, xz.pm25) %*%
        solve ( rbind ( cbind (xx.pm25 + u.pm25/ki.pm25[i], xz.pm25), cbind (zx.pm25, zz)), tol =1e-25)
      
      observed <- c( bar.pm25[i]-Mu.w.pm25, Z[i,]-bar.Z)
      expect.pm25 <- cbind (expect.pm25, Mu.w.pm25 + project.matrix %*% as.matrix (observed))

      #Calibrate PM10
      project.matrix <- cbind (xx.pm10, xz.pm10) %*%
        solve ( rbind ( cbind (xx.pm10 + u.pm10/ki.pm10[i], xz.pm10), cbind (zx.pm10, zz)), tol =1e-25)
      
      observed <- c( bar.pm10[i]-Mu.w.pm10, Z[i,]-bar.Z)
      expect.pm10 <- cbind (expect.pm10, Mu.w.pm10 + project.matrix %*% as.matrix (observed))
      
      #Calibrate Coarse
      project.matrix <- cbind (xx.coarse, xz.coarse) %*%
        solve ( rbind ( cbind (xx.coarse + u.coarse/ki.coarse[i], xz.coarse), cbind (zx.coarse, zz)), tol =1e-25)
      
      observed <- c( bar.coarse[i]-Mu.w.coarse, Z[i,]-bar.Z)
      expect.coarse <- cbind (expect.coarse, Mu.w.coarse + project.matrix %*% as.matrix (observed))
    }
    
    expect.pm25 <- t(expect.pm25)
    expect.coarse <- t(expect.coarse)
    expect.pm10 <- t(expect.pm10)
    
    check.dat <- data.frame (monitor.data[[1]]$date[id.use],bar.pm10, bar.pm25, bar.coarse, expect.pm10, expect.pm25, expect.coarse)
    names (check.dat) <- c("date", "pm10.Ave", "pm25.Ave", "pmCoarse.Ave", "pm10.RC","pm25.RC", "pmCoarse.RC")

    admitCVD <- admissions.junk$admitp6 + admissions.junk$admitp7 + admissions.junk$admitp3 + admissions.junk$admitp4 + admissions.junk$admitp5
    admitResp <- admissions.junk$admitp8 + admissions.junk$admitp9

    admissions <- data.frame (admissions.junk$date, admissions.junk$agecat, admissions.junk$denom, admissions.junk$tmpd, admissions.junk$rmtmpd,admissions.junk$dow, admissions.junk$dptp, admissions.junk$rmdptp, admissions.junk$pmCoarse, admissions.junk$pm10c, admissions.junk$pm25c,  admitCVD, admitResp, admissions.junk$admitp11, admissions.junk$pm25tmean2)

    names (admissions) <- c("date", "agecat", "denom", "tmpd", "rmtmpd", "dow", "dptp", "rmdptp", "pmCoarse", "pm10c", "pm25c", "admitCVD", "admitResp", "admitp11", "pm25tmean2")

    final.data <- merge (admissions, check.dat, "date", "date", all.x = TRUE)

    return (final.data)
}


##############################
#
#CountyEstimate

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 
fips.108 <- read.csv ("fips.108.csv", colClasses=c("character"))[,1]

#Outcomes: Aggregated cardiovascular and respiratory emergency admissions
outcomes <- c("CVD", "Resp")
models <- c("model.pm25.coarse")

#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.108){
  
  junk <- load (paste ("~/CoarsePM/Data/Monitor.Specific/",siteName, ".monitor.Rdata", sep = ""))
  if ( ncol (monitor.data[[3]]) > 2) {
  
    admissions <- Reg.Cal (siteName)
        
    for (model in models){  #Loop over models#
      
      p2 <- NULL
      if (model == "model.pm25.coarse") {p1 <- "pmCoarse.RC"; p2 <- "pm25.RC"} #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 of If         
      }  ## 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 as.character(unique (info$fips)) ) {
    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.RC.Rdata" )