# ------------------------------------------------------------
# Summer ozone constrained distributed lag model
#
# Three lag variables enter the model (specified by the matrix
# M below).  These are:
#
# avg0to6    - the mean of lag0, lag1, ..., lag6 ozone
# avg3minus7 - mean of lag0, lag1 and lag2 minus avg0to6
# avg0minus3 - lag0 minus avg3minus7
#
# Since columns 2 and 3 of M sum to zero, and column 1 sums
# to 1, the estimate associated with avg0to6 is the
# distributed lag coefficient of interest.
#
# See [1] for more details.
#
# [1] Michelle L. Bell, Aidan McDermott,
#     Scott L. Zeger, Jonathan M. Samet, Francesca Dominici.
#     Ozone and mortality in 95 U.S. urban communities 
#     from 1987 to 2000. Accepted for publication in the 
#     Journal of the American Medical Association. 
# ------------------------------------------------------------

options (object.size=Inf, memory=Inf)

# change the next two lines to suit your local environment
attach("/users/project6/aphena/APHENA/.Data",pos=2)
source("support.s")

OzoneDL <- function(cityname)
{
  
# Prepare the data taking only variables of interest:
  vars <- c("death","agecat","date","dow","o3tmean",
            "tmpd","dptp","rmtmpd","rmdptp","time",
            "markdeath")

  data <- get(paste(cityname,"8700",sep=""))[,vars]

# Summer model: Apr to Oct -- need month
  month <- floor(data[,"date"]/100) %% 100

# Df are modified by the amount of missing data:
# Set extreme mortality values to NA
# Set winter values to NA
# -------------------------------------------------

  death                          <- data[,"death"]
  death[data[,"markdeath"] == 1] <- NA
  death[month < 4 | month > 10 ] <- NA
  data[,"death"]                 <- death

# Names used in DL variable labels (taken from column names of M)
# --------------------------------------------------------------- 
  names <- c("avg0to6","avg3minus7","avg0minus3")

  M <- matrix(c( 1/7, 1/3 - 1/7,  1 - 1/3,
                 1/7, 1/3 - 1/7,    - 1/3,
                 1/7, 1/3 - 1/7,    - 1/3, 
                 1/7,     - 1/7,        0,
                 1/7,     - 1/7,        0,
                 1/7,     - 1/7,        0,
                 1/7,     - 1/7,        0),
       ncol=3, nrow=7, byrow=T, dimnames=list(NULL,names))

  return(DLmodel(dataframe=data,cause="death",pollutant="o3tmean",M=M))
}


# Vector of cities we wish to analyze
O3List <- c( "la",   "ny",   "chic", "dlft", "hous", "sand",
             "staa", "phoe", "det",  "miam", "phil", "seat",
             "sanj", "clev", "sanb", "pitt", "oakl", "atla",
             "sana", "rive", "denv", "sacr", "stlo", "buff",
             "clmo", "cinc", "stpe", "kan",  "hono", "tamp",
             "memp", "indi", "nwk",  "balt", "salt", "roch",
             "wor",  "orla", "jckv", "fres", "loui", "bost",
             "birm", "dc",   "okla", "prov", "elpa", "taco",
             "aust", "dayt", "jers", "bake", "akr",  "char",
             "nash", "tuls", "gdrp", "no",   "stoc", "albu",
             "syra", "tole", "ral",  "wich", "colo", "batr",
             "mode", "madi", "spok", "ltrk", "knox", "shr",
             "desm", "ftwa", "corp", "jcks", "hunt", "lex",
             "arlv", "kans", "bidd", "cdrp", "clmg", "covt",
             "john", "lafy", "lasv", "linc", "lkch", "milw",
             "mobi", "musk", "oma",  "port", "tucs")

# Collect estimate and se in a nx2 matrix
O3DLResults <- matrix(NA,ncol=2,nrow=length(O3List))

for ( j in 1:length(O3List) ) {
  cat(" ----- city",j,"---",O3List[j],"-----\n")
  modelOut <- OzoneDL(O3List[j])
  print(modelOut$coef)
  
  if ( ! is.null (modelOut$coef) ) {
    O3DLResults[j,] <- modelOut$coef["avg0to6",1:2]
  }
}