### ADAPTING Farrington TO JAPAN DATA ### ### Mortality data # death_week: 疫学週 # week_ending_date: 疫学週の最終日 # count: 観測死亡者数 library(tidyr) library(ggplot2) library(dplyr) library(gcookbook) library(cowplot) library(scales) library(surveillance) library(lubridate) data_shibo <- read.csv("your_directory/data.csv") dl =list() for(i in unique(data_shibo$prefectureEN)){ df.pref <- data_shibo[data_shibo$prefectureEN==i,] ### 重みづけ用 #ind.2020 <- grep(2020, df.pref$week_ending_date) B <- new("sts", observed = df.pref$count, state = rep(0, nrow(df.pref)), start = c(2010, 1), freq = 52) #B <- new("sts", observed = tmp, state = rep(0, nrow(df.pref)), start = c(2010, 1), freq = 52) ## 分析パラメータの設定 ## b:戻る年月 ## w:見る幅 control1 <- list(range = grep("2017|2018|2019|2020", df.pref$week_ending_date), b = 5, w = 3, weightsThreshold = 2.58, pThresholdTrend = 1, noPeriods = 10, populationOffset=FALSE, fitFun="algo.farrington.fitGLM.flexible", trend = TRUE, pastWeeksNotIncluded = 26, alpha = 0.05, thresholdMethod = "muan") ## 分析 fit <- farringtonFlexible(B, control = control1) df <- tidy.sts(fit) %>% mutate(alarm_date = ifelse(alarm, date, NA) %>% as_date) df$EB = attr(fit,'control')$expected y_shift <- 0.01*max(df$observed) p <- ggplot(df, aes(date, observed)) + 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 = upperbound), color = "#1bc066", size = 1.5) + geom_line(aes(y = ceiling(EB)), color = "#FF6F00", size = 1.5) + geom_point(aes(x = alarm_date, y = observed + y_shift), 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$observed + 2*max(df$observed)*0.01)*1.1) , expand = FALSE) + labs(x = "日付", y = "死亡者数\n") + ggtitle( unique(df.pref$prefectureJP[df.pref$prefectureEN==i]) ) dl[[i]]=df quartz(file = paste('Farrington_',unique(df.pref$prefectureJP[df.pref$prefectureEN==i]),'.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$epoch>157,]$observed-x[x$epoch>157,]$EB<0,0,x[x$epoch>157,]$observed-x[x$epoch>157,]$EB)) })) cbind(c(1:47),sapply(dl,function(x){ sum(ifelse(x[x$epoch>157,]$observed-x[x$epoch>157,]$upperbound<0,0,x[x$epoch>157,]$observed-x[x$epoch>157,]$upperbound)) }))