Age-Sex Standardized Poisson-Based Interrupted Time Series Analysis

## Housekeeping
rm(list=ls())


## Install packages
if(!require("pacman")) install.packages("pacman")
pacman::p_load(WordR, MASS, plyr, foreign, data.table,lubridate, rJava, WordR, ReporteRs, fakeR, car,
  flextable, WordR, rJava, ReporteRs, ggplot2, dplyr, devtools, formula.tools, stringr, rlang, lmtest)
install_github("tlcaputi/tlcPack"); library(tlcPack)


## Set working directory where the data is
setwd("C:/Users/tcapu/Google Drive/UCC/manuscript/SmokingBanLungCancer/analyze")


y <- readRDS("./input/data_new.rds")

y <- y %>%
      rename(
        lungMortality = lungMortalityCases,
        lungIncidence = lungIncidenceCases,
        brainMortality = brainMortalityCases,
        brainIncidence = brainIncidenceCases,
        smoke_rate = total_smoking_rate,
        smoke_frequency = total_smoking_frequency,
        age_sex_smoke_rate = age_sex_smoking_rate,
        age_sex_smoke_frequency = age_sex_smoking_frequency,
        smoke_price = real_price
      )


y <- y %>%
      mutate(lungMortality = ifelse(lungMortality=='..', 0, lungMortality)) %>%
      mutate(
        year = as.numeric(year),
        population = as.numeric(population),
        lungMortality = as.numeric(lungMortality),
        lungIncidence = as.numeric(lungIncidence),
        brainMortality = as.numeric(brainMortality),
        brainIncidence = as.numeric(brainIncidence),
        int2004 = ifelse(year<2004, "preban", "postban"),
        year_from_baseline = year - min(year)
      ) %>% filter(sex!="Both sexes", year <=2015)


y <- y  %>% mutate(
  lungIncidence = ifelse(is.na(lungIncidence), 0, lungIncidence),
  lungMortality = ifelse(is.na(lungMortality), 0, lungMortality),
  brainIncidence = ifelse(is.na(brainIncidence), 0, brainIncidence),
  brainMortality = ifelse(is.na(brainMortality), 0, brainMortality)
)


y <- y %>% group_by(year, sex) %>% mutate(totalPopulation = sum(population)) %>%
            ungroup()

y <- y %>% mutate(
          lungIncidenceRate = lungIncidence/population,
          lungMortalityRate = lungMortality/population,
          brainIncidenceRate = brainIncidence/population,
          brainMortalityRate = brainMortality/population,
          weighting = population/totalPopulation
        ) %>% mutate(
          lungIncidenceSumUp = lungIncidenceRate * weighting,
          lungMortalitySumUp = lungMortalityRate * weighting,
          brainIncidenceSumUp = brainIncidenceRate * weighting,
          brainMortalitySumUp = brainMortalityRate * weighting
        )

y <- y %>% group_by(year) %>% mutate(
  lungIncidence_ASR = sum(lungIncidenceSumUp),
  lungMortality_ASR = sum(lungMortalitySumUp),
  brainIncidence_ASR = sum(brainIncidenceSumUp),
  brainMortality_ASR = sum(brainMortalitySumUp)
) %>% ungroup()


## Now we select only the variables we need for analysis
y <- y %>%
    select(
      sex, # Male or Female
      year, # Year
      year_from_baseline, # Year since 1994
      ageGroup, # 5-year age groups from 40 to 85 and then 85+
      population, # number of Irish people in the sex X age group
      pm2p5, # Yearly particulate matter 2.5 concentration
      pm2p5_lag3, # 3-year lagged particulate matter 2.5 concentration
      smoke_rate, # Annual Irish smoking prevalence
      smoke_frequency, # Annual avg. number of cigs per day
      lungMortality, # Cases of lung mortality within age X sex X year
    ) %>% arrange(ageGroup, sex, year) %>% ungroup()

y %>% head()

basefit_formula <- as.formula("lungMortality ~ ageGroup + sex + offset(log(population)) + smoke_rate + smoke_frequency")

basefit <- glm(basefit_formula, data = y, family = poisson)

summary(basefit)

durbinWatsonTest(basefit, max.lag=3)

lag_names <- c()
  y$age_sex_group <- paste0(y$ageGroup, y$sex)

  for(k in unique(y$age_sex_group)){
    for (j in 1:3){
      nam <- paste0("lungMortality_lag", j)
      y[y$age_sex_group==k,nam] <- lag(as.numeric(unlist(y[y$age_sex_group==k, "lungMortality"])), j)
      lag_names <- u(lag_names, nam)
    }
  }

mod_w_lag_formula <- as.formula(paste(c(basefit_formula, lag_names), collapse="+"))
mod_w_lag <- lm(mod_w_lag_formula, data = na.omit(y))

step1 <- step(mod_w_lag, na.rm=T)
lags <- grep("lag", names(step1$coefficients), value=T)


if(length(lags)>0){
    formula_w_lags <- as.formula(paste(c(as.character(basefit_formula), lags), collapse="+"))
} else {
    formula_w_lags <- basefit_formula
}

model_w_lags <- glm(formula_w_lags, data = y, family = poisson)

summary(model_w_lags)

## Perform a GLM and return only the interaction term
glm_int <- function(formula, data, family){
  mod <- glm(formula, data, family=family)
  variables <- all.vars(mod$terms)
  mod_coeff <- data.frame(summary(mod)$coefficients)
  int_variable <- grep("year_from_baseline:int", rownames(mod_coeff), value=T)
  mod_coeff$rownames <- rownames(mod_coeff)
  z_value <- mod_coeff[grep("year_from_baseline:int", rownames(mod_coeff)),"z.value"]
  mat <- as.data.frame(matrix(NA, nrow=1, ncol=2))
  names(mat) <- c("int_year", "z_val")
  mat$int_year <- as.numeric(gsub("year_from_baseline:int|postban", "", int_variable))
  mat$z_val <- as.numeric(z_value)
  return(mat)
}

## Take a list of glm_int objects and return the most significant negative value.
most_significant_value <- function(list_of_models){
  data <- Reduce(function(x,y) {rbind.data.frame(x, y)}, list_of_models)
  which_is_min <- which(data$z_val==min(data$z_val, na.rm=T))
  min_year <- data[which_is_min, "int_year"]
  if(min(data$z_val, na.rm=T) > - 1.96) warning("No significant results!")
  print(paste0(min_year, ": z-value is", min(data$z_val, na.rm=T)))
  return(as.numeric(min_year))
}


for(i in 2003:2010){
  interruption_name <- paste0("int", i)
  values <- ifelse(y$year < i, "preban", "postban")
  y[,interruption_name] <- factor(values, levels=c("preban", "postban"))
}



# Test models with interruption points between 0 and 6 years after the ban
  int2004_model <- glm_int(as.formula(paste0(formula_w_lags, "+year_from_baseline*int2004")), y, family=poisson)
  int2005_model <- glm_int(as.formula(paste0(formula_w_lags, "+year_from_baseline*int2005")), y, family=poisson)
  int2006_model <- glm_int(as.formula(paste0(formula_w_lags, "+year_from_baseline*int2006")), y, family=poisson)
  int2007_model <- glm_int(as.formula(paste0(formula_w_lags, "+year_from_baseline*int2007")), y, family=poisson)
  int2008_model <- glm_int(as.formula(paste0(formula_w_lags, "+year_from_baseline*int2008")), y, family=poisson)
  int2009_model <- glm_int(as.formula(paste0(formula_w_lags, "+year_from_baseline*int2009")), y, family=poisson)
  int2010_model <- glm_int(as.formula(paste0(formula_w_lags, "+year_from_baseline*int2010")), y, family=poisson)

  models_to_test <- list(
    int2004_model,
    int2005_model,
    int2006_model,
    int2007_model,
    int2008_model,
    int2009_model,
    int2010_model
  )


int <- most_significant_value(models_to_test)
interruption <- paste0("int", int)

print(interruption)

final_model_formula <- as.formula(paste0(formula_w_lags, "+ year_from_baseline *", interruption))
final_model <- glm(final_model_formula, data = y, family = poisson)
print(as.character(final_model_formula))

print(summary(final_model))

## Sex Group

sex_final_formula <- as.formula(
    paste0("lungMortality ~ ageGroup + offset(log(population)) + smoke_rate + smoke_frequency + ",
          paste(lags, collapse="+"),
          " + year_from_baseline * ", interruption)
)


for (i in unique(y$sex)){
    print(i)
    temp_mod <- glm(sex_final_formula, data = y, family = poisson)
    print(summary(temp_mod))
}


## Age Groups

age_final_formula <- as.formula(
    paste0("lungMortality ~ sex + offset(log(population)) + smoke_rate + smoke_frequency + ",
          paste(lags, collapse="+"),
          " + year_from_baseline * ", interruption)
)


# for (i in unique(y$ageGroup)){
#     print(i)
#     temp_mod <- glm(age_final_formula, data = y, family = poisson)
#     print(summary(temp_mod))
# }

counterfactual <- y
counterfactual[,interruption] <- "preban"

  y <- y %>% mutate(
    crude_rate_orig = (lungMortality)/population,
    with_smoking_ban = predict(final_model, type="response", y),
    without_smoking_ban = predict(final_model, type="response", counterfactual),
    with_smoking_ban_rate = predict(final_model, type="response", y)/population,
    without_smoking_ban_rate = predict(final_model, type="response", counterfactual)/population
  )




  ## aggregate dataset from age X sex X year observations to year observations to create a figure
  fig_dataset <- y %>% group_by(year) %>%
                transmute(
                  crude_rate = sum(lungMortality)/sum(population), #weighted.mean(original_rate, population, na.rm=T),
                  crude_rate_wtd = weighted.mean(crude_rate_orig, population, na.rm=T),
                  with_smoking_ban_rate = weighted.mean(with_smoking_ban_rate, population, na.rm=T),
                  without_smoking_ban_rate = weighted.mean(without_smoking_ban_rate, population, na.rm=T),
                  without_smoking_ban_rate_direct = sum(without_smoking_ban)/sum(population),
                  with_smoking_ban_rate_direct = sum(with_smoking_ban)/sum(population)
                ) %>% ungroup()

  fig_dataset <- fig_dataset %>%
                      select(year, crude_rate, crude_rate_wtd, with_smoking_ban_rate, without_smoking_ban_rate,
                             with_smoking_ban_rate_direct, without_smoking_ban_rate_direct) %>%
                      unique()
  post_interruption <- fig_dataset %>% filter(year >= int - 1)


## Set parameters for Jupyter Notebook
par(mar=c(5.1,5.1,4.1,2.1))
options(repr.plot.width=8.5, repr.plot.height=5.5)

  plot(fig_dataset$year, fig_dataset$crude_rate_wtd, xlab="Year", ylab="Lung Cancer Mortality\n(per 1,000 people aged 40+)",
   ylim=c(min(c(fig_dataset$without_smoking_ban_rate,fig_dataset$with_smoking_ban_rate,fig_dataset$with_smoking_ban_rate_direct),na.rm=T),
          max(c(fig_dataset$without_smoking_ban_rate,fig_dataset$with_smoking_ban_rate,fig_dataset$with_smoking_ban_rate_direct),na.rm=T))
  )

  lines(fig_dataset$year, fig_dataset$with_smoking_ban_rate, type="l",col="blue")
  lines(post_interruption$year, post_interruption$without_smoking_ban_rate, type="l",lty="dashed",col="red")
  abline(v=int,col="black")
  abline(v=2004,lty="dashed",col="darkgreen")
  legend("bottomleft",legend=c("Predicted Rate Without Smoking Ban","Standardized-Adjusted Rate",
                            paste0(int," Modeled Interruption"),"2004 Smoking Ban", "Observed Crude Rate"),
         col=c("red","blue","black","darkgreen","black"),lty=c(2,1,1,2, NA),pch=c(NA,NA,NA,NA, 1))

  ## Calculate fitted and observed values and estimate differences using bootstrapping
  fitted_values <- as.vector(filter(y, !is.na(without_smoking_ban),
                    !is.na(with_smoking_ban)) %>% ungroup() %>% select(without_smoking_ban)  ) %>%
                    unlist() %>% as.numeric()
  observed_values <- as.vector(filter(y, !is.na(without_smoking_ban),
                    !is.na(with_smoking_ban)) %>% ungroup() %>% select(with_smoking_ban) ) %>%
                    unlist() %>% as.numeric()

## Set seed to replicate
set.seed(1234)

  B <- replicate(1000, expr={
    ix <- sample(1:length(observed_values), length(observed_values), replace=T)
    (sum(fitted_values[ix])-sum(observed_values[ix])) # /sum(comb_preds$mean[ix])
  })


### Create Summary


    ## Rounding function
    r2 <- function(x) round(x,2)

    ## Creates a string that looks like a confidence interval
    ci95 <- function(a, b, c) paste0(r2(a), " (95%CI ", r2(b), " to ", r2(c), ")")



  post_interruption_cases <- r2(sum(y[y$year >=int, "lungMortality"],na.rm=T))
  total_difference <- r2(sum(y$without_smoking_ban - y$with_smoking_ban, na.rm=T))
  total_difference_lo95 <- r2(quantile(B,0.025, na.rm=T)) #sum(d$wos_raw - d$wsm_raw, na.rm=T))
  total_difference_hi95 <- r2(quantile(B,0.975, na.rm=T)) # sum(d$wos_raw - d$wsm_raw, na.rm=T))


  percent_cases_avoided <- ci95(
    total_difference/sum(y[y$year >=int, ]$without_smoking_ban, na.rm=T) * 100,
    total_difference_lo95/sum(y[y$year >=int, ]$without_smoking_ban, na.rm=T) * 100,
    total_difference_hi95/sum(y[y$year >=int, ]$without_smoking_ban, na.rm=T) * 100
  )

  total_cases_avoided <- ci95(
    total_difference,
    total_difference_lo95,
    total_difference_hi95
  )

  yearly_cases_avoided <- ci95(
    total_difference/(max(y$year - int + 1, na.rm=T)),
    total_difference_lo95/(max(y$year - int + 1, na.rm=T)),
    total_difference_hi95/(max(y$year - int + 1, na.rm=T))
  )


  numSaved <- paste0("The total number of cases avoided is ",total_cases_avoided,".")
  numSavedYr <- paste0("The number of cases avoided per year is ",yearly_cases_avoided,".")
  intYear <- paste0("The interruption year is ",int)
  postIntDeaths <- paste0("The total number of post-interruption cases is ",post_interruption_cases,".")
  pctPostIntDeaths <- paste0("The ratio of cases avoided to post-interruption cases (*100) is ",percent_cases_avoided,".")


## Print summary
    print(intYear)
    print(numSaved)
    print(numSavedYr)
    print(postIntDeaths)
    print(pctPostIntDeaths)