今回は、参考の時系列分析をコロナ禍の死亡数の経年変化に適用する目的で実施してみた。
残念ながら、データが2018年までしか公表されていないが、分析してみて見えそうかどうかを見極めておこうと思う。
過日、データ公表時点で実施しようと思う。
今回得られた予測（赤）と実数（青）のグラフは以下のものとなります。
※結果的に1万人程度の変動なら季節変動を超えて見えそうです
【参考】
①[R] 実データを時系列解析して結果を考察してみる
②データセット一覧(5-4死亡月別にみた年次別死亡数及び死亡率（人口千対）)@統計で見る日本
③b2010_sgo1j.xls@過去の生産・出荷・在庫・在庫率指数（接続指数）@鉱工業指数 集計結果又は推計結果（データダウンロード）
###やったこと
・時系列解析をやってみる
・死亡数に適用
###・時系列解析をやってみる
コードは参考①の方のものをまるまる使わせていただきました。
一応、今回利用したものはおまけに置いておきます。
酒類の総合原指数月次付加価値額生産(IIP)をデータ変更して清涼飲料と肉製品もやってみました。
一つ理解できたことは、年間の大きな変化はそれぞれ相関があり、また年をまたがっての大きな増加や減少はそれぞれ特徴的に存在する。
ということと、その場合予測関数はそれぞれ異なるパラメータが選ばれるということです。
今回は、ここの説明をせずに死亡数のみに注目して記載しようと思います。
###・死亡数に適用
####データは周期性
と仮定して以下でオブジェクト化
IIP <- ts(data,start=c(2008),frequency=12)
この処理をすると、データにjanみたいなデータが付与され、
[1] "head(IIP)"
Jan Feb Mar Apr May Jun
2008 110384 102413 102778 94055 91597 84258
と出力されます。
####要素分解
decompose()関数が変動を以下の3つの変動に分解してくれる
【参考】
・decompose@RDocumentation
- トレンド+周期変動
- トレンド(年平均の変動)
- 周期変動
- 季節変動
- 不規則変動
$Y_t=T_t + S_t + e_t$
または、
$Y_t=T_tS_te_t$
出力は元のデータを加えた4つのグラフを出力します。このグラフから以下のことが読み取れます。
- 日本は冬場の死亡数が多い
- 2011年は2万人程度の大きな不規則変動を表している。2012年も5千人以上の変動がある
- 2010年から2011年にかけて大きく1万人程度増加しているのが二つ目のトレンドとして見える
つまり、COVID-19の場合も、ほんとに存在すればこの程度には見えると思います。
ADF検定
詳細は参考を見てください。どちらも難しいけど、。。
この検定により、以下の結果が得られます。
【参考】
・Rで計量時系列分析：単位根過程、見せかけの回帰、共和分、ベクトル誤差修正モデル
・第７回： 時系列分析@データ解析(pdf)
Augmented Dickey-Fuller Test
data: IIP
Dickey-Fuller = -8.4775, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary
この検定は、参考①によれば、「ADF検定とは、帰無仮説「過程が単位根である」が棄却されれば定常過程であるとする検定です。」
ということで、この系はp-value = 0.01で、「過程が単位根であるという帰無仮説が棄却され、定常過程である」ということが言えそうです。
自己相関
自己相関の図は以下が得られました。強い相関が長期間あることが分かります。
偏自己相関
一方、偏自己相関は、以下のとおり得られました。こちらの相関性は自己相関よりは短期間ですが、少なくとも前年のデータとの相関は95%信頼区間を超えて相関が強いことが分かります。
つまり周期性が強いことが分かります。
####モデルの自動適合
実際、プログラムを動かせば一応フィッティング出来てパラメータが得られ、予測曲線が描画できる。
しかし、以下の参考のWikipediaによれば、このパラメータの持つ意味は大きく、得られた曲線はいろいろな意味を持った曲線であることが分かる。
今回Fittingに利用するクライテリアは、参考①と同じく以下のとおりとしました。
AIC = −2 log(L) + 2(p + q + P + Q + k)
【参考】
・auto.arima; Fit Best ARIMA Model To Univariate Time Series@RDocumentation
・Autoregressive integrated moving average@wikipedia
・Automatic Time Series Forecasting:the forecast Package for R(pdf)
結果は以下のとおりとなりました。
Series: IIP
ARIMA(3,0,0)(1,1,1)[12]
Coefficients:
ar1 ar2 ar3 sar1 sma1
0.4601 0.3585 0.1453 -0.0736 -0.8764
s.e. 0.1115 0.1051 0.0996 0.1174 0.2171
sigma^2 estimated as 11226869: log likelihood=-1036.46
AIC=2084.91 AICc=2085.74 BIC=2101.01
そして、「そこで、残差に自己相関が無い適切なモデルか判断します。」
ということで以下の残差のプロット、残差の自己相関、Liung-Box検定の結果は以下のとおりです。参考①の酒類の生産量に比較すると2011年の残差が大きいこと、また残差の自己相関の信頼区間越えが発生しているので適切でない可能性もあります。また、p-valueは小さめです。しかし一応95%信頼区間以上の値となっており、全体としては妥当なモデルと言えると思います。
※残差の自己相関が一か所信頼区間以上の値になっているのが気になりますが、。。2011年の異常値との関係を別途究明しようと思います
⇒これについてはおまけに記載しました
そして、2018年の1月から12月の死亡数の予測が以下のように求められました。
Jan Feb Mar Apr May Jun Jul Aug
2018 133082.8 117476.9 121778.5 111193.0 108672.5 100391.7 103671.2 105367.6
Sep Oct Nov Dec
2018 102307.3 110265.9 113573.0 123402.7
この予測結果と実績を比較した図が以下のとおりとなりました。
ほぼ、参考①と同程度の予測ができています。
###まとめ
・時系列解析を死亡数の経年変化に適用してみた
・要素分解により、トレンド、季節変動、そして不規則変動に分離でき、不規則変動から異常検知が出来そうである
・未来予測については一応計算でき、実績とほぼ似たような予測値がえられた
・予測関数の精度をさらに検証する必要がある。特に死亡数は2011年辺りに異常値があり、これとの関係を明らかにしたい。
⇒おまけに記載したように2011年のデータ異常の所為という結論が得られました。
###おまけ
####2011年のデータ異常問題
邪道ですが2011年のデータを2012年のデータで置き換えると、「残差のプロット、残差の自己相関、Liung-Box検定の結果」は以下のようになりました。
一応、残差のプロットも小さくなり、残差の自己相関も信頼区間以下になり、検定のp-valueも青い線をクリアできました。
データを入れ替えた結果、deconpose()も以下のように変化しました。
すなわち、不規則変動が小さくなり、毎年の変動が大きく見えるようになりました。
一方、自己相関と偏自己相関は以下のとおり強い相関と周期性を示しています。
そして、以下のようにフィッティング関数もARIMA(4,0,0)(0,1,1)[12]になり、予測値も変更されました。
Series: IIP
ARIMA(4,0,0)(0,1,1)[12]
Coefficients:
ar1 ar2 ar3 ar4 sma1
0.6306 0.1294 0.0569 0.1644 -0.8182
s.e. 0.1004 0.1190 0.1176 0.0980 0.1373
sigma^2 estimated as 5568777: log likelihood=-996.19
AIC=2004.39 AICc=2005.22 BIC=2020.48
その結果、予測値は以下のようになりました。
より予測値は実測値に近づいたように見えます。
ということでやはり、2011年の異常値が悪さをしていたという結論が得られます。
####コードについて
以下は、参考①のものをまるまる引用させていただいています。
なお、画像ファイル格納などは今回追加しております。
入力のcsvファイルは、単に分析する単一列を見出し付きでコピペして、張り付けて作成しました。
data <- read.csv("./industrial/sibou.csv",encoding = "utf-8")
IIP <- ts(data,start=c(2008),frequency=12)
print("dim(IIP)")
print(dim(IIP))
print("head(IIP)")
print(head(IIP))
# 要素分解
pngFileName <- paste("./industrial/IIP_compose_sibou.png",sep="")
png(file=pngFileName,res=400, w=1000, h=1500,pointsize=5)
plot(decompose(IIP))
dev.off()
# ADF検定
library(tseries)
print(adf.test(IIP))
# 自己相関
pngFileName <- paste("./industrial/IIP_acf_sibou.png",sep="")
png(file=pngFileName,res=400, w=1000, h=600,pointsize=4)
acf(IIP,lag.max=120,ylim=c(-.2,1.0)) #96
dev.off()
# 偏自己相関
pngFileName <- paste("./industrial/IIP_pacf_sibou.png",sep="")
png(file=pngFileName,res=400, w=1000, h=600,pointsize=4)
pacf(IIP,lag.max=120,ylim=c(-.2,1.0))
dev.off()
# モデルの自動適合
# auto.arima関数からSARIMA(3,0,0)(1,1,1)[12]が推定される
library(forecast)
Fit <- auto.arima(IIP,ic="aic",trace=F,stepwise=F,approximation=F,allowmean=F,allowdrift=F);Fit
print(Fit)
# モデルの推定
# 推定されたSARIMA(3,0,0)(1,1,1)[12]をarima関数で作成する
Model <- arima(IIP,order=c(3,0,0),seasonal=list(order=c(1,1,1),period = 12));Model
# 残差検定の診断
pngFileName <- paste("./industrial/IIP_tsdiag_sibou.png",sep="")
png(file=pngFileName,res=400, pointsize=3)
tsdiag(Model)
dev.off()
# 予測値の取得
Model.pred <- predict(Model,12)
value <- (Model.pred$pred);value
print(value)
data_p <- read.csv("./industrial/sibou_p.csv",encoding = "utf-8")
IIP_p <- ts(data_p,start=c(2008),frequency=12)
data_n <- read.csv("./industrial/sibou_n.csv",encoding = "utf-8")
IIP_n <- ts(data_n,start=c(2008),frequency=12)
pngFileName <- paste("./industrial/IIP_compose_pred_sibou.png",sep="")
png(file=pngFileName,res=400, w=1000, h=600,pointsize=4)
plot(IIP_p,col="red",ylim=c(80000,140000)) #,ylim=c(800000,1400000)
par(new=T)
plot(IIP_n,col="blue",ylim=c(80000,140000))
dev.off()