時系列データの解析（厚労省公開の医療費データ）

この本に沿って時系列データの解析方法をまとめました。

サンプルデータを使っても面白くないので、厚労省が公開している医療費のデータを使いました。

厚労省の医療費データベース

例によってこのデータはエクセルで公開されていて、そのまま解析できる状態じゃありません。

今回は入院の総医療費だけを扱ったので、その部分だけ加工してcsvにしました。

一応、加工したデータはダウンロードのページに置いてます。

それでは、解析していきます。

まずはデータ読み込みと加工。

Iryouhi <- read.csv("医療費.csv", as.is = T)
Nyuin <- ts(Iryouhi[, 2], frequency = 12, start = c(1984, 4))

read.csv関数でcsvファイルを読み込み、ts関数で時系列データの形に変換します。

時系列プロットにしてイメージを掴みます。

plot(Nyuin, yaxt = "n", ylab = "入院医療費（1ヵ月毎）", xlab = "年月", main = "入院医療費の経時変化")
axis(2, c(5*10^11, 7.5*10^11, 10^12), c("5000億", "7500億", "1兆"))

2010年では一ヵ月に約1兆円かかっています。

単純計算で年間約12兆、入院は全体の3分の1だと考えると全体は36兆。

大体計算は合いそうですね。

ギザギザなっているのは周期的な季節変動。

全体的に医療費が上がっているのは分かるのですが、季節変動を取り除くとどうなるかやってみます。

NyuinTS <- decompose(Nyuin)
plot(NyuinTS)

decompose()関数で勝手に季節変動を取り除いてくれます。

上から元のデータ、季節変動を除いたトレンド、季節変動、ランダムな成分となっています。

トレンドだけを抽出するとこのようになります。

plot(NyuinTS$trend)

移動平均を行っても、同じような結果になります。

NyuinMA12 <- filter(Nyuin, rep(1, 12))
plot(NyuinMA12)

このトレンドは一見するとただ単に毎年上がっているように見えるのですが、2000年までの医療費の上がり方とそれ以降の上がり方が違うように見えます。その事をさらに確認するべく、自己相関プロットを描いてみます。

acf(Nyuin)

右肩上がりで上昇しているトレンドを持つ時系列データでは、自己相関プロットはこのようにLag（間隔）が大きくなるにつれて下がっていくような図になります。

入院医療費が全体的に上昇するトレンドを持っているので、このような図が得られたのです。

これを、入院医療費の期間を区切って描いてみます。

par(mfrow=c(3, 2))
acf(window(Nyuin, start = c(1985, 4), end = c(1990, 3)), main = "1985-1990")
acf(window(Nyuin, start = c(1990, 4), end = c(1995, 3)), main = "1990-1995")
acf(window(Nyuin, start = c(1995, 4), end = c(2000, 3)), main = "1995-2000")
acf(window(Nyuin, start = c(2000, 4), end = c(2005, 3)), main = "2000-2005")
acf(window(Nyuin, start = c(2005, 4), end = c(2010, 3)), main = "2005-2010")
par(mfrow=c(1, 1))

時系列を5年ごとに区切って自己相関プロットを描いています。

結果を見ると明らかに1985年から2000年までの3つのプロットと2000年から2010年までの2つのプロットの傾向が違うことが分かります。

2000年から2005年の間は自己相関が低く、この期間の時系列に上昇傾向のトレンドが弱いことが分かります。この時期を思い返してみると、小泉元首相の構造改革の時期とほぼ重なっているように思えます。

構造改革で小泉さんは「医療費を抑制する」ことを打ち出していましたが、こんなところでそれがデータに現れているのではないでしょうか？

また、2005年から2010年の自己相関もあまり高くありません。

1985年から2000年までの医療費上昇率より、ここ5年間の医療費上昇率の方が低いことが伺えます。

というより、過去の上昇率が異常に高かったのでは？

さらに議論をするためには、多方面からの分析が必要でしょう。

さて、時系列分析に戻ります。

時系列データにモデルを当てはめ、将来の予測を行うこともできます。

NyuinSARIMA <- arima(Nyuin, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 1), period = 12))
nseason         <- 120
NyuinPredSARIMA <- predict(NyuinSARIMA, n.ahead = nseason)
ts.plot(Nyuin, NyuinPredSARIMA$pred, col = c(1, "blue"), main = "SARIMA予測", 
	ylab = "入院医療費（1ヵ月毎）", xlab = "年月")
text(1990, c(5*10^11, 7.5*10^11, 10^12, 1.25*10^12), c("5000億", "7500億", "1兆", "1.25兆"))

これは季節ARIMA（SARIMA）モデルで将来の予測をしたものです。

ARIMAモデルはAR（自己相関）, 差分の回数, MA（移動平均）の3つをモデル化したもので、それぞれの次数によって予測結果は変わります。ですので、今回の予測結果が全てではありませんが、どちらにせよこのまま医療費が上がっていくことは予想されるでしょう。

最後に、その他の時系列解析のコードを載せておきます。

#------ホルト・ウィンタース法（汎用的な時系列モデル）
NyuinHW1 <- HoltWinters(Nyuin)
NyuinHW2 <- HoltWinters(Nyuin, seasonal = "mult")
NyuinHW3 <- HoltWinters(Nyuin, gamma = 0, beta = 0)
plot(NyuinHW1)
plot(NyuinHW2)
plot(NyuinHW3)

nseason      <- 1200
NyuinPred1    <- predict(NyuinHW1, n.ahead = nseason)
NyuinPred2    <- predict(NyuinHW2, n.ahead = nseason)
NyuinPred3    <- predict(NyuinHW3, n.ahead = nseason)
ts.plot(Nyuin, NyuinPred1, col = c(1, "blue"), main = "HW1予測")
ts.plot(Nyuin, NyuinPred2, col = c(1, "blue"), main = "HW2予測")
ts.plot(Nyuin, NyuinPred3, col = c(1, "blue"), main = "HW3予測")


#------ARIMA、SARIMA（→信頼区間は？）

#---order:AR（自己相関）, 差分の回数, MA（移動平均）
#---予測のプロットを見ると、このARIMAはちょっとおかしい

NyuinARIMA  <- arima(Nyuin, order = c(1, 1, 1))
NyuinSARIMA <- arima(Nyuin, order = c(1, 1, 1), seasonal = list(order = c(1, 1, 1), period = 12))

plot(Nyuin - NyuinARIMA$residuals, main = "ARIMA")
plot(Nyuin - NyuinSARIMA$residuals, main = "SARIMA")

#---残差の検定（Ljung-Box検定）→有意でなければホワイト・ノイズ
plot(NyuinARIMA$residuals, main = "残差（ARIMA）")
Box.test(NyuinARIMA$residuals)
Box.test(NyuinSARIMA$residuals)

#---lagを変化させて、「s次まで自己相関が無い」ことをチェックする
tsdiag(NyuinARIMA)
tsdiag(NyuinSARIMA)


nseason         <- 120
NyuinPredARIMA  <- predict(NyuinARIMA,  n.ahead = nseason)
NyuinPredSARIMA <- predict(NyuinSARIMA, n.ahead = nseason)
ts.plot(Nyuin, NyuinPredARIMA$pred, col = c(1, "blue"), main = "ARIMA予測")
ts.plot(Nyuin, NyuinPredSARIMA$pred, col = c(1, "green"), main = "SARIMA予測", 
	ylab = "入院医療費（1ヵ月毎）", xlab = "年月")
text(1990, c(5*10^11, 7.5*10^11, 10^12, 1.25*10^12), c("5000億", "7500億", "1兆", "1.25兆"))



#------単位根検定（修正ディッキー・フラー検定、ADF検定）
#---Call: lm()とあるように、線形回帰モデルを利用している
library(urca)
summary(ur.df(Nyuin))
summary(ur.df(Nyuin, type = "trend"))
summary(ur.df(Nyuin, type = "trend", lag = 12))
summary(ur.df(Nyuin, type = "drift", lag = 12))




#------スペクトル
#---ペリオドグラム（高速フーリエ変換、spans=で修正Daniell平滑化を行う）
spectrum(Nyuin, method = "pgram")


#---ARモデル
#---order=でARの次数を指定できる、指定しなければAICが最小の次数を選ぶ
spectrum(Nyuin, method = "ar")


#---フィルター後のスペクトル
NyuinFilter <- filter(Nyuin, rep(1, 12))
NyuinFilter <- NyuinFilter[!is.na(NyuinFilter)]
spectrum(NyuinFilter, method = "ar")


#---※スルツキー効果：差分をとることで見せかけの周期性が見られる現象