### ADAPTING EuroMOMO TO JAPAN DATA ### 
### Mortality data

# death_week: 疫学週
# week_ending_date: 疫学週の最終日
# count: 観測死亡者数


library(ggplot2)
library(lubridate)

# Confirmed data
dat <- read.csv("your_directory/data.csv")

# Seasonality, 2 each for 1 year and half a year
library(dplyr)
dat.M <- dat %>% mutate(sin52 = sin((2*pi/(365.25/7)) * death_week),
                          cos52 = cos((2*pi/(365.25/7)) * death_week),
                          sin26 = sin((4*pi/(365.25/7)) * death_week),
                          cos26 = cos((4*pi/(365.25/7)) * death_week))
dat.M$year <- as.numeric(substr(dat.M$week_ending_date,6,9))
dat.M <- subset(dat.M, year %in% 2012:2020)

# Prefecture as list
dl <- list(0)
for (i in seq(47)){
  dpref <- subset(dat.M, prefecture==i)
  dl[i] <- list(dpref)
}

for (i in seq(47)){
  # Baseline 
  di <- dl[[i]]
  # Fit trend and seasonality, a priori, unlike FluMOMO
  m <- glm(count ~ death_week + sin52 + cos52 + sin26 + cos26, 
           family = quasipoisson, data=di[di$death_week<526,])
  
  di$EB <- exp(predict.glm(m, newdata=di, se.fit=TRUE)$fit)
  di$VlogB <- predict.glm(m, newdata=di, se.fit=TRUE)$se.fit^2
  di$Vcount <- with(di, EB * max(1, sum(residuals(m, type = "deviance")^2)/df.residual(m))) 

  # Baseline estimation variances - not on log scale
  di$VB <- with(di, (exp(VlogB)-1)*exp(2*log(EB) + VlogB))
  
  
  # Baseline 2/3 residual confidence intervals (FluMOMO)
  di$RVB <- with(di, ((2/3)*(EB^(2/3-1))^2)*Vcount + ((2/3)*(EB^(2/3-1))^2)*VB)
  di[with(di, is.na(RVB) | is.infinite(RVB)), "RVB"] <- 0
  di$EB_95L <- with(di, 
                    pmax(0,sign((sign(EB)*abs(EB)^(2/3))
                                -1.96*sqrt(RVB))*abs((sign(EB)*abs(EB)^(2/3))
                                                     -1.96*sqrt(RVB))^(3/2)))
  di$EB_95U <- with(di, 
                    sign((sign(EB)*abs(EB)^(2/3))
                         +1.96*sqrt(RVB))*abs((sign(EB)*abs(EB)^(2/3))
                                              +1.96*sqrt(RVB))^(3/2))
  
  # Excess by one-sided 5% confidence interval (vs. Farrington)
  di$EB_95one <- with(di, 
                      sign((sign(EB)*abs(EB)^(2/3))
                           +1.645*sqrt(RVB))*abs((sign(EB)*abs(EB)^(2/3))
                                                 +1.645*sqrt(RVB))^(3/2))
  
 
  # to prepare plots
  di$alarm = ifelse(di$count>di$EB_95one,1,0)
  di$week_ending_date = as.Date(di$week_ending_date,"%d%b%Y")
  di = di[di$year >2016,]

  df <- di %>%
  	mutate(alarm_date = as_date(ifelse(alarm, week_ending_date, NA)))

  p <- ggplot(df, aes(week_ending_date, count)) + 
  	theme_classic(base_family = "HiraKakuProN-W3") +
  	theme(text = element_text(size = 18), 
  	      panel.grid.major.y = element_line(colour = "gray90"),
  	      axis.text.x = element_text(angle = 40, vjust = 1, hjust = 1)) +
  geom_col(width = 5, fill = '#2292f0') +
  geom_line(aes(y = EB_95one), color = "#1bc066", size = 1.5) +
  geom_line(aes(y = ceiling(EB)), color = "#FF6F00", size = 1.5) +
  geom_point(aes(x = alarm_date, y = count + max(count)*0.01), color = '#c01b75', shape = 3, stroke = 1.5) +
  	scale_x_date(date_labels = '%Y/%m', date_breaks = '3 month', expand = c(0, 0)) +
  	coord_cartesian(ylim = c(0, max(df$count + 2*max(df$count)*0.01)*1.1), expand = FALSE) +
  	labs(x = '日付', y = '死亡者数\n')+
  	ggtitle(di$prefectureJP[1])

  dl[[i]] = di
  quartz(file = paste('EuroMOMO_',di$prefectureJP[1],'.pdf',sep=""), type = "pdf", family = "HiraKakuProN-W3", width = 30, height = 20)
  print(p)
  dev.off()
}







# to create summary table
cbind(c(1:47),sapply(dl,function(x){
  sum(ifelse(x[x$death_week>521,]$count-x[x$death_week>521,]$EB<0,0,x[x$death_week>521,]$count-x[x$death_week>521,]$EB))
}))

cbind(c(1:47),sapply(dl,function(x){
  sum(ifelse(x[x$death_week>521,]$count-x[x$death_week>521,]$EB_95one<0,0,x[x$death_week>521,]$count-x[x$death_week>521,]$EB_95one))
}))