こちらでインタラクティブなグラフを操作できます
https://yuta-vet.shinyapps.io/shiny_app_deploy/?_ga=2.106368643.1882585030.1624518165-33520689.1624518165
初めに
割と避けてきたトピックですが、オリンピックの開催が決定したことだし、時系列データとして予測してみたいと思います。Rを使ってアルゴリズムは、Auto Arima、Neural Network、TBATSを使います。
まずは、東京都のサイトからデータを取ってきます。日にちごとの感染者数が載っていませんので加工するひと手間があります。（意図的？）
https://catalog.data.metro.tokyo.lg.jp/dataset/t000010d0000000068/resource/c2d997db-1450-43fa-8037-ebb11ec28d4c
日本語が入っている状態でcsvで読み込むとエラーがでることが多いので、xlsxにしてから読み込ませます。列名として公表_年月日を選択する場合は ` ` で囲む必要があります。Pの隣にあるやつです。シングルクォーテーションではありません。
日ごとにグループにして、その日のデータ数をカウントします。
library(tidyverse)
library(forecast)
df <- readxl::read_xlsx("130001_tokyo_covid19_patients.xlsx")
# df <- df$公表_年月日 はエラーがでる
df <- as.data.frame(df$`公表_年月日`) # 1列だとdataframeじゃなくなる
colnames(df) <- "date"
df <- df %>% group_by(date) %>% summarise(number=n())この状態だと、陽性判明者数が０の日のデータは現れてきません。今後、予測に使うアルゴリズムたちは欠損値を嫌います。欠損値は、interpolateした方がいいです。
まず、すべての日にちが入ったデータフレームを作って、そこに今あるデータをレフトジョインさせます。
date_df <- as.data.frame(data.frame(matrix(rep(NA, 1), nrow=1))[numeric(0), ])
colnames(date_df) <- c("date")
date_series <- seq(from=as.Date("2020/1/24", format="%Y/%m/%d"),
by =1,
to = as.Date("2021/6/23", format="%Y/%m/%d"))
for(i in 1:517){
date_df[i,1] <- as.character(date_series[i])
}
date_df$date <- as.Date(date_df$date,format="%Y-%m-%d")
df2 <- left_join(date_df,df,by="date")
nrow(date_df) # 517
nrow(df2) # 517
nrow(df) # 491
sum(is.na(df2$number)) # 26これでdf2にはすべての日にちが入った（２０２０年１月２４日～２０２１年６月２３日）データが完成しました。ただ、まだNAがいるのでinterpolateします。interpolateする１番簡単な方法は、dataframeからtsの形にデータを変えて、tsclean( ts データ)をすれば一発なのですが、それではわかりにくいので自分でやってみます。その日の数値がNAだった場合は、前後の日にちのデータを平均した値を入れます。もし、後ろの日にちもNAだった場合には前の日のデータを入れます。
for( i in 2:(nrow(df2)-1) ){
if ( is.na(df2$number[i]) ) {
if ( is.na(df2$number[i+1]) ){
df2$number[i] <- df2$number[i-1]
}else{
df2$number[i] <- (df2$number[i+1]+df2$number[i-1])/2
}
}
}
sum(is.na(df2$number)) # 0さっそくグラフを見てみます。
df2$date <- as.Date(df2$date, format="%Y-%m-%d")
ggplot(df2, aes(x=date,y=number))+
geom_line()+
theme_minimal()+
scale_x_date(date_breaks ="21 day")+
theme(axis.text.x = element_text(angle = 90))
６月２３日においては、第３波が終わったところにいます。
さて、AutoArimaを使って予測してみたいと思います。まずは、tsのデータ型にします。c(4,6) ：4は第4週で、６は、日、月、火、水、木、金の金を指します。
tsdata <- ts(df2$number, start = c(4,6), frequency =7) autoarimaに学習させます。forecast パッケージの関数です。forecastで６０日後まで予想させます。autoplot、autolayerもforecastパッケージです。
arima_fit <- auto.arima(tsdata)
autoplot(tsdata)+
autolayer(forecast(arima_fit,h=60),color="lightgreen")
autoplotは簡単なのですが、tsにしてしまうと日付の情報がなくなってしまいます。そのためデータフレームに戻してggplotでグラフにします。prediction interval 80%を示しています。
df_arima <- as.data.frame(forecast(arima_fit,h=60))
df_arima$date <- seq(from=as.Date("2021/6/24"), by=1,length.out=60)
colnames(df_arima) <- c("PointForecast_arima","Lo80_arima","Hi80_arima","Lo95_arima","Hi95_arima","date")
ggplot()+
geom_line(data=df2, aes(x=date,y=number))+
geom_line(data=df_arima,
aes(x=date,y=PointForecast_arima),color="red")+
geom_ribbon(data=df_arima,
aes(x=date,y=PointForecast_arima,ymin = Lo80_arima, ymax = Hi80_arima), alpha = 0.2)+
scale_x_date(date_breaks ="21 day")+
theme(axis.text.x = element_text(angle = 90))+
ggtitle("ARIMA FORECAST")+
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5))
予測のところを拡大したいです。df2を日付でfilterします。
ggplot()+
geom_line(data=df2 %>% filter(date > "2021/3/1" & date <= "2021/6/23"),
aes(x=date,y=number))+
geom_line(data=df_arima,
aes(x=date,y=PointForecast_arima),color="red")+
geom_ribbon(data=df_arima,
aes(x=date,y=PointForecast_arima,ymin = Lo80_arima, ymax = Hi80_arima), alpha = 0.2)+
scale_x_date(date_breaks ="7 day")+
theme(axis.text.x = element_text(angle = 90))+
ggtitle("ARIMA FORECAST")+
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5))
増えてしまう予想ですね。ワクチン等を加味するとまた違うかもしれません。xregに変数として入れてautoarimaすればできます。オリンピック加味すると、全然違うでしょう。