RStan:時系列データの欠損値推定

# 結果のプロット
ggplot(results_df, aes(x = time)) +
  # 1. 観測値のデータ（比較用）
  geom_line(aes(y = y_true), color = "darkgray", linewidth = 1) +

  # 2. 推定された潜在状態(mu)の95%信用区間
  geom_ribbon(aes(ymin = mu_lower, ymax = mu_upper), fill = "skyblue", alpha = 0.3) +

  # 3. 推定された潜在状態(mu)の平均値
  geom_line(aes(y = mu_mean), color = "blue", linewidth = 1) +

  # 4. 観測されたデータ点
  geom_point(aes(y = y_obs), color = "black", size = 2.5, na.rm = TRUE) +

  # 5. 欠損観測値の推定値（平均値）
  geom_point(data = . %>% filter(is_missing), aes(y = y_imputed_mean), color = "red", size = 3) +

  # 6. 欠損観測値の推定値の95%信用区間
  geom_errorbar(data = . %>% filter(is_missing), aes(ymin = y_imputed_lower, ymax = y_imputed_upper), color = "red", width = 0.5) +

  # ラベルとテーマ
  labs(
    title = "時系列データの欠損値推定結果 (RStan)",
    subtitle = "赤点が推定された欠損値、青いリボンが潜在状態の95%信用区間",
    x = "時間 (Time)",
    y = "値 (Value)"
  ) +
  theme_minimal()

# 欠損していた箇所のデータのみを抽出
missing_comparison_df <- results_df %>%
  filter(is_missing == TRUE) %>% # is_missingがTRUEの行（欠損点）のみをフィルタリング
  select(
    time,
    y_true, # 実際の観測値（欠損前の値）
    y_imputed_mean, # 推定された欠損観測値（事後平均値）
    y_imputed_lower, # 推定値の95%信用区間の下限
    y_imputed_upper # 推定値の95%信用区間の上限
  )

# 散布図を作成して比較
ggplot(missing_comparison_df, aes(x = y_true, y = y_imputed_mean)) +
  # 推定値の95%信用区間をエラーバーで表示
  geom_errorbar(aes(ymin = y_imputed_lower, ymax = y_imputed_upper), color = "gray", width = 0) +
  # 推定値の平均を点でプロット
  geom_point(color = "red", size = 3, alpha = 0.8) +
  # 理想的な直線 (y = x) を破線で追加
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "blue") +
  # 軸の範囲を調整して見やすくする
  # coord_fixed(
  #   xlim = range(c(missing_comparison_df$y_true, missing_comparison_df$y_imputed_lower)),
  #   ylim = range(c(missing_comparison_df$y_true, missing_comparison_df$y_imputed_upper))
  # ) +
  labs(
    title = "観測値と推定値の比較（欠損点のみ）",
    x = "実際にもっていた値 (y_true)",
    y = "Stanによる推定値 (y_imputed)",
    caption = "エラーバーは推定値の95%信用区間を示します"
  ) +
  theme_minimal()

# 欠損点のみのデータで、時系列順にプロット
ggplot(missing_comparison_df, aes(x = factor(time))) +
  # 推定値の95%信用区間
  geom_errorbar(aes(ymin = y_imputed_lower, ymax = y_imputed_upper), width = 0.4, color = "skyblue") +
  # 推定値の平均値
  geom_point(aes(y = y_imputed_mean), shape = 4, size = 4, color = "blue") +
  # 実際の値（真値）
  geom_point(aes(y = y_true), shape = 16, size = 3, color = "red") +

  # 凡例のためのダミーaes
  geom_point(aes(y = y_imputed_mean, color = "Stanによる推定値"), alpha = 0) +
  geom_point(aes(y = y_true, color = "実際の観測値"), alpha = 0) +
  scale_color_manual(
    name = "凡例",
    values = c("Stanによる推定値" = "blue", "実際の観測値" = "red"),
    guide = guide_legend(
      override.aes = list(
        shape = c(4, 16),
        size = c(4, 3),
        alpha = c(1, 1),
        color = c("blue", "red")
      )
    )
  ) +
  labs(
    title = "欠損点における観測値と推定分布の比較",
    x = "時間 (欠損していた時点)",
    y = "値"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")