ggplot2でforecastインスタンスを描画する

最近時系列データを forecast パッケージを使って処理している。便利！

library(forecast)

d <- AirPassengers
d.arima <- forecast::auto.arima(d)
d.forecast <- forecast(d.arima, level = c(95), h = 50)

plot(d.forecast)

f:id:sinhrks:20141004134746p:plain

この結果を、総称関数 plot と同じように ggplot で描画したいが、そのまま ggplot に渡すとエラーになる。

library(ggplot2)
ggplot(data = d.forecast) + geom_line(aes(x = x, y = fitted))
# エラー: ggplot2 doesn't know how to deal with data of class forecast

ggplot では内部で fortify という関数を使い、入力データを data.frame に変換している。ので forecast 用にも fortify を定義してやればよい。

forecastには forecast オブジェクト用の as.data.frame があるが、これは予測値/信頼区間のみしか data.frame 化しない。こちらに原系列を追加してやる必要がある。

fortify.forecast <- function(forecast.data) {
  require(dplyr)
  forecasted <- as.data.frame(forecast.data)
  forecasted$Time <- as.Date(time(forecast.data$mean))

  fitted <- data.frame(Time = as.Date(time(forecast.data$fitted)),
                       Original = forecast.data$x,
                       Fitted = forecast.data$fitted)

  rownames(fitted) <- NULL
  rownames(forecasted) <- NULL

  dplyr::rbind_list(fitted, forecasted)
}

head(ggplot2::fortify(d.forecast))
#         Time Original   Fitted Point Forecast Lo 95 Hi 95
# 1 1949-01-01      112 111.9353             NA    NA    NA
# 2 1949-02-01      118 117.9664             NA    NA    NA
# 3 1949-03-01      132 131.9662             NA    NA    NA
# 4 1949-04-01      129 128.9774             NA    NA    NA
# 5 1949-05-01      121 120.9892             NA    NA    NA
# 6 1949-06-01      135 134.9782             NA    NA    NA

日時, 原系列, 推定値を Time, Original, Fitted カラムとして追加するようにした。これで ggplot から通常の data.frame と同じように扱える。スペースを含む列名はバッククオートでエスケープする。

ggplot(data = d.forecast) +
  geom_line(mapping = aes_string(x = 'Time', y = 'Original')) +
  geom_line(mapping = aes_string(x = 'Time', y = '`Point Forecast`'), colour='blue') +
  geom_ribbon(mapping = aes_string(x = 'Time', ymin = '`Lo 95`', ymax = '`Hi 95`'), alpha = 0.5)

f:id:sinhrks:20141004135949p:plain