はじめに

こんにちは、事業会社で働いているデータサイエンティストです。

ビジネスでもアカデミアでも、皆さんは大量の欠損を含む時系列データを扱ったことがありますか？比較的わかりやすい例として、集計が一時的に停止した重要な指標が挙げられます。実は、因果推論の観点では、処置を受けた時系列において「処置を受けなかった場合の時系列」も欠損データの一種とみなされます。

詳しくはこちらの記事をぜひご確認ください：

今回の記事では、このような欠損データに対して、ベイズ時系列モデルを用いた補完手法を紹介します。もちろん、データサイエンスは魔法ではありません。十分な情報がなければ補完は困難ですが、相関のありそうな他の時系列データを活用することで、より精度の高い補完が可能になります。本記事で紹介する手法は、因果推論におけるDID・Causal Impact・合成統制法といった推定法の発想にも通じるものです。

モデル説明

本記事で説明するモデルは、決して複雑なものではありません。シンプルに説明すると、欠損の有無に関わらず、手元にあるデータから複数の時系列の背後に隠されたトレンドが存在し、そのトレンドによってデータが生成されたと考えるのが本モデルの根幹です。

まず、トレンドの生成から説明しましょう。

T期のトレンドは、T-1期のトレンドを平均とする正規分布から生成されるとします：

$$
trend_{t} \sim Normal(trend_{t - 1}, \sigma_{trend})
$$

分散などの事前分布は省略します。実装のコードはご自身で確認してください。

次に、観測されたデータ $series$ は次のように生成されるとします：

$$
series_{i,t} \sim Normal(\sum_{l = 1}^{ラグ}trend_{t - l} * \beta_{i, l}, \sigma_{i})
$$

ここで、どのデータがどのように欠損しているかを特別に考慮する必要はありません。

モデルの推定が終わったら、欠損している$series$を上の$series_{i,t} \sim Normal(\sum_{l = 1}^{ラグ}trend_{t - l} * \beta_{i, l}, \sigma_{i})$にしたがってサンプリングすれば復元できます。

では、早速実際のデータ分析に入りましょう！

データ

今回利用するデータは、Google トレンドにおける「アルバイト」と「求人」の検索ボリュームです。2020年2月から2025年2月までのデータを対象としますが、2023年1月1日以前の「求人」のデータを意図的に欠損させます。これにより、モデルが「アルバイト」の全データと一部しか残存していない「求人」のデータから、どの程度正確に「求人」データを補完できるかを検証できます。

では、まずダウンロードしたデータを整理して、マスターテーブルも作りましょう：

# GoogleトレンドのデフォルトだとCSVに入らないデータが入るのでskipに2を指定する
df <- readr::read_csv("アルバイト.csv", skip = 2) |>
  `colnames<-`(c("date", "volumn")) |>
  dplyr::mutate(
    search_word = "アルバイト"
  ) |>
  dplyr::bind_rows(
    readr::read_csv("求人.csv", skip = 2) |>
      `colnames<-`(c("date", "volumn")) |>
      dplyr::mutate(
        search_word = "求人"
      )
  ) |>
  dplyr::mutate(
    volumn = volumn/100,
    train_status = dplyr::case_when(
      search_word == "アルバイト" ~ 1,
      search_word == "求人" & date >= as.Date("2023-01-01") ~ 1,
      # 2023-01-01より前の「求人」のデータを学習データから外す
      TRUE ~ 0
    )
  )


search_word_master <- df |>
  dplyr::select(search_word) |>
  dplyr::distinct() |>
  dplyr::arrange(search_word) |>
  dplyr::mutate(
    search_word_id = dplyr::row_number()
  )


date_master <- df |>
  dplyr::select(date) |>
  dplyr::distinct() |>
  dplyr::arrange(date) |>
  dplyr::mutate(
    date_id = dplyr::row_number()
  )

df_with_id <- df |>
  dplyr::left_join(
    date_master, by = "date"
  ) |>
  dplyr::left_join(
    search_word_master, by = "search_word"
  )

モデル実装コード

モデルのコードはこちらになります。

Googleトレンドのデータに合わせた変数名になっていますが、適宜変えていただければと思います！

impute_stan

data {
  int time_type;
  int search_word_type;
  int lag_type;
  
  int N;
  array[N] int time;
  array[N] int search_word;
  array[N] real volumn;
}
parameters {
  real<lower=0> sigma_begin_trend;
  real<lower=0> sigma_trend;
  vector[time_type] trend;
  
  matrix<lower=0,upper=1>[search_word_type, lag_type] sigma_beta;
  matrix[search_word_type, lag_type] beta;
  vector<lower=0>[search_word_type] sigma_volumn;
}
model {
  sigma_begin_trend ~ gamma(0.01, 0.01);
  sigma_trend ~ gamma(0.01, 0.01);
  
  trend[1] ~ normal(0, sigma_begin_trend);
  for (i in 2:time_type){
    trend[i] ~ normal(trend[i - 1], sigma_trend);
  }
  
  to_vector(sigma_beta) ~ gamma(0.01, 0.01);
  to_vector(beta) ~ normal(0, to_vector(sigma_beta));
  
  for (i in 1:N){
    real m = 0.0;
    for (l in 1:lag_type){
      if (time[i] - l < 1) {
        m = m;
      } else {
        m += beta[search_word[i], l] * trend[time[i] - l];
      }
    }
    volumn[i] ~ normal(m, sigma_volumn[search_word]);
  }
}
generated quantities {
  array[search_word_type] vector[time_type] predict;
  for (i in 1:search_word_type){
    for (j in 1:time_type){
      real m = 0.0;
      for (l in 1:lag_type){
        if (j - l < 1){
          m = m;
        } else {
          m += beta[i, l] * trend[j - l];
        }
      }
      predict[i, j] = normal_rng(m, sigma_volumn[i]);
    }
  }
}

モデル推定

ここで、データをlistに整形して、Stanのモデルをコンパイルします：

data_list <- list(
  time_type = nrow(date_master),
  search_word_type = nrow(search_word_master),
  lag_type = 3,
  
  N = nrow(df_with_id[which(df_with_id$train_status == 1),]),
  time = df_with_id$date_id[which(df_with_id$train_status == 1)],
  search_word = df_with_id$search_word_id[which(df_with_id$train_status == 1)],
  volumn = df_with_id$volumn[which(df_with_id$train_status == 1)]
)

m_impute_init <- cmdstanr::cmdstan_model("impute.stan")

では、ここで実際に推定します：

> m_impute_estimate <- m_impute_init$variational(
     data = data_list
 )
------------------------------------------------------------ 
EXPERIMENTAL ALGORITHM: 
  This procedure has not been thoroughly tested and may be unstable 
  or buggy. The interface is subject to change. 
------------------------------------------------------------ 
Gradient evaluation took 0.001993 seconds 
1000 transitions using 10 leapfrog steps per transition would take 19.93 seconds. 
Adjust your expectations accordingly! 
Begin eta adaptation. 
Iteration:   1 / 250 [  0%]  (Adaptation) 
Iteration:  50 / 250 [ 20%]  (Adaptation) 
Iteration: 100 / 250 [ 40%]  (Adaptation) 
Iteration: 150 / 250 [ 60%]  (Adaptation) 
Iteration: 200 / 250 [ 80%]  (Adaptation) 
Success! Found best value [eta = 1] earlier than expected. 
Begin stochastic gradient ascent. 
  iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes  
   100      -362161.995             1.000            1.000 
   200        67430.270             3.685            6.371 
   300       -19030.075             3.971            4.543 
   400        66801.591             3.300            4.543 
   500       168158.952             2.760            1.285 
   600       149391.016             2.321            1.285 
   700       -24000.042             3.022            1.285 
   800       186125.578             2.785            1.285 
   900       204797.144             2.486            1.129 
  1000       172978.329             2.256            1.129 
  1100       165409.917             2.160            1.129   MAY BE DIVERGING... INSPECT ELBO 
  1200       219351.507             1.548            0.603   MAY BE DIVERGING... INSPECT ELBO 
  1300       219398.493             1.093            0.246   MAY BE DIVERGING... INSPECT ELBO 
  1400       223578.956             0.967            0.184   MAY BE DIVERGING... INSPECT ELBO 
  1500       222163.667             0.907            0.126   MAY BE DIVERGING... INSPECT ELBO 
  1600       153057.020             0.940            0.184   MAY BE DIVERGING... INSPECT ELBO 
  1700       217591.157             0.247            0.184 
  1800       208635.495             0.138            0.091 
  1900       186339.196             0.141            0.120 
  2000       222801.604             0.139            0.120 
  2100       220690.562             0.136            0.120 
  2200       214941.011             0.114            0.043 
  2300       222558.094             0.117            0.043 
  2400       205890.330             0.123            0.081 
  2500       226311.134             0.132            0.090 
  2600       211470.682             0.093            0.081 
  2700       224441.388             0.070            0.070 
  2800       223101.494             0.066            0.070 
  2900       178675.866             0.079            0.070 
  3000       184041.057             0.065            0.058 
  3100       216623.372             0.079            0.070 
  3200       200619.321             0.085            0.080 
  3300       219890.168             0.090            0.081 
  3400       223310.615             0.084            0.080 
  3500       213703.819             0.079            0.070 
  3600       218541.446             0.074            0.058 
  3700       210543.925             0.072            0.045 
  3800       214769.663             0.074            0.045 
  3900       213825.695             0.049            0.038 
  4000       222794.546             0.050            0.040 
  4100       214793.941             0.039            0.038 
  4200       221871.108             0.034            0.037 
  4300       215968.264             0.028            0.032 
  4400       216378.540             0.027            0.032 
  4500       220822.748             0.024            0.027 
  4600       219218.552             0.023            0.027 
  4700       219732.754             0.019            0.020 
  4800       218793.745             0.018            0.020 
  4900       217565.969             0.018            0.020 
  5000       218753.417             0.014            0.007   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  12.0 seconds.

12秒で終わりました！

推定結果を保存します：

m_impute_summary <- m_impute_estimate$summary()

モデル推定結果

ここでは、まず可視化用のデータフレイムを作成します：

predict_result_df <- m_impute_summary |> 
  dplyr::filter(stringr::str_detect(variable, "^predict\\[")) |>
  dplyr::mutate(
    id = variable |>
      purrr::map(
        \(x){
          as.integer(stringr::str_split(x, "\\[|\\]|,")[[1]][2:3])
        }
      )
  ) |>
  tidyr::unnest_wider(id, names_sep = "_") |>
  dplyr::filter(id_1 == 2) |>
  dplyr::bind_cols(
    date = date_master$date,
    volumn = df_with_id |>
      dplyr::filter(search_word == "求人") |>
      dplyr::pull(volumn)
    ) |>
  dplyr::mutate(
    train_status = dplyr::case_when(
      date >= as.Date("2023-01-01") ~ 1,
      TRUE ~ 0
    ) |>
      as.factor()
  )

では、実際に可視化しましょう！

predict_result_df |>
  ggplot2::ggplot() +
  ggplot2::geom_point(ggplot2::aes(x = date, y = volumn, color = train_status)) + 
  ggplot2::geom_line(ggplot2::aes(x = date, y = mean)) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = date, ymin = q5, ymax = q95), fill = "blue", alpha = 0.3)

train_statusが1のデータは学習に使用した実データ、train_statusが0のデータはモデルに与えていない欠損データを示します。線は予測の平均値を表し、青い区間は予測の信用区間を示しています。

完璧ではないものの、「アルバイト」のデータと一部の「求人」データを活用することで、「求人」の欠損データをかなり高い精度で復元できたのではないかと思います。

次に、X軸に実データ、Y軸にモデルの予測値をプロットします：

predict_result_df |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = volumn, y = mean, color = train_status)) + 
  ggplot2::geom_abline(intercept = 0, slope = 1, color = "red") +
  ggplot2::xlim(c(0, 1)) + 
  ggplot2::ylim(c(0, 1))

この図を見るとさらに分かりやすいですが、実データと予測値はほぼ 45 度線付近に分布しており、数値データの予測モデルとして理想的な結果になっています。

結論

いかがでしたか？

本記事で紹介したモデルを活用すれば、集計組織の都合で生じた大量の欠損データはもちろん、因果推論において必然的に欠損する「処置群に処置がなかった場合」のデータも復元可能です。これにより、高い精度と透明性を持った因果推論が実現でき、意思決定の強力な武器となります。

最後に、私たちと一緒に、データサイエンスの力で社会を改善したい方はこちらをご確認ください：