### 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)) }))