######################################################################
## Pool Estimates

## 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"

## A subdirectory "Emergency.pooled.results" will be created where
## outcome-lag-region specific posterior estimates will be stored
## as individual .Rdata files

#Load county estimates
load ("Results.Emergency.Rdata")

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

info$lag <- as.numeric (as.character(info$lag))
info$outcome <- as.character (info$outcome)
fips.include <- as.character(unique (info$fips))

### estimation of second-stage model parameters
regions <- c("All", "East", "West")

dir.create ("Emergency.pooled.results")
setwd ("Emergency.pooled.results")

find.max <- function (x){ #Find maximum variance
    y <- max(diag(x))}


library (tlnise)
for (region in regions){
  for (model in models){
    for(outcome_  in outcomes){   
        for(lag_  in lags){

          if (region != "All"){
            use <- info$model == model & info$outcome == outcome_ & info$lag == lag_ & info$regionEW == region
            use[is.na(use)] <- FALSE} else{
            use <- info$model == model & info$outcome == outcome_ & info$lag == lag_}
            
            param.beta <- betas[use]
            param.var <- vars[use]
            param.beta <-  matrix (unlist (param.beta), nrow = sum(use), ncol = length(param.beta[[1]]), byrow = T )
            param.var <- array(unlist(param.var), c(dim(param.var[[1]]), sum(use)) )
                   
            param.beta <- na.omit(param.beta)                 ## rm.na for tlnise
            param.var <- na.omit(param.var)
            
          
            param.beta[1:sum(use),] <- param.beta[order(apply (param.var, 3, find.max)),]       ## sort by variance
            param.var[,,1:sum(use)] <- param.var[,,order(apply (param.var, 3, find.max))]
            id <- apply (param.var, 3, find.max)

            param.beta <- as.matrix(param.beta[id < 2,])
            param.var <- param.var[,,id < 2]

            removed <- NULL
            estimate <- NULL
            fips.rm <- NULL
            beta.rm <- NULL
            var.rm <- NULL
            
            estimate <- tlnise(param.beta,param.var,
                               rep(1,nrow(param.beta)),
                               seed=1234, prnt = FALSE)
            estimate <- list( gamma=estimate$gamma, theta=estimate$theta,
                             SDtheta=estimate$SDtheta,A=estimate$A,
                             rtA=estimate$rtA, Dgamma=estimate$Dgamma,
                             Vtheta=estimate$Vtheta, B0=estimate$B0,
                             lr=estimate$lr, fips.used= info$fips[use], 
                             beta=param.beta, var=param.var,
                             removed=removed)  

            ## save estimates in external file
            filename <- paste(model,outcome_ ,lag_, region, "RData",sep=".")
            save(estimate,file=filename)

        } ## end for lag_
    } ## end for outcome_
  } ## end of models
} ## end of region