はじめに
AIを使うと、データ分析のPythonコードを、爆速で勝手に書いてくれるようになりました。
CSVを読み込み、前処理をして、回帰分析を実行し、結果を表にまとめる。こうした一連のコードが、かなり短時間でそれらしく出てきます。
ただ、コードが動いたことと、その分析が自分の意図に合っていることは別です。
今回も、AIが生成した回帰分析コードを読んでいたら、「これ、勝手にこんなことまでやっていいのか？」と思う箇所がありました。特に、あとで出てくる HC3 ロバスト標準誤差です。「ロバスト標準誤差を使ってください」なんて一言も指示した覚えがありません。
便利ではあります。ただ、こういうところが怖くもあります。
この記事は、「全部分かっている人が教える記事」ではありません。AIが生成した回帰分析コードを、自分なりに読み直しながら、「このコードは何をしようとしているのか」「自分がやりたい分析に合っているのか」を確認してみる記録です。というか、愚痴です。
今回は、YouTube動画データを使った初期的な重回帰分析コードを題材にします。テキストマイニング的なタイトル特徴量の部分は扱わず、重回帰分析の基本的な処理だけを見ます。
このYouTube分析シリーズの記事はこちらです。今後の記事もここに追加していきます。
- 第1回: YouTube動画はどんな要素が再生パフォーマンスに関連しているのか？
- 第2回: YouTube動画は「同じチャンネル内で当たる」のか？ロジスティック回帰で分析する
- 第3回: この記事
- 第4回: YouTube動画は「どのチャンネルか」でどこまで変わるのか？固定効果で回帰分析を見直す
- 第5回: 固定効果の次に、階層ベイズで「チャンネル差」をもう少し丁寧に見る
- 第6回: 広告クリエイティブを個人で集めようとして分かったこと
もともとは、第1回の記事でYouTube動画の再生パフォーマンスを分析したときに生成したPythonコードです。
この記事では、その分析結果そのものではなく、「その分析コードが何をしていたのか」を読み直します。
全体像
今回のコードは、おおまかにこういう流れです。
生データCSV
YouTube動画メタデータ
↓
pandasで読み込み
DataFrame化
↓
分析用の列を作る
日時・動画尺・統制変数
↓
モデル式を組む
A0からA3へ段階的に追加
↓
statsmodelsに渡す
formula + DataFrame
↓
OLSを推定
.fit()
↓
HC3ロバスト標準誤差
標準誤差・p値を補正
↓
係数表・適合表を出力
CSV / Markdown
この図で見ると、今回の分析は「CSVを入れたら回帰結果が出る」という一本道ではありません。
途中で、目的変数をどう作るか、説明変数をどう作るか、どの変数を統制するか、標準誤差をどう扱うか、という判断が入っています。たぶん、ここが一番大事です。
前半の pandas や numpy の処理は、比較的読みやすいです。どの列からどの列を作っているかを追えば、だいたい雰囲気はつかめます。
問題は後半です。statsmodels に入れた瞬間、Pythonの表操作から、統計モデルの世界に入ります。
コードが何をしようとしているかを見る
今回、最初に思ったのは、コードを完全に理解できるかどうかよりも、そのコードが何をしようとしているかを見る方が大事なのでは、ということです。
たとえば、「動画の特徴と1日あたり再生数の関連を見たい」と思っていました。だとすると、目的変数が1日あたり再生数ベースになっていることは、その意図と合っています。
一方で、もし目的変数が累積再生数そのものになっていたら、公開からの日数の違いにかなり引っ張られるはずです。その場合、コードは動いていても、やりたいこととは少しズレているかもしれません。
この確認は、必ずしも高度な数式を知らないとできないものではありません。
むしろ、まずはこういう問いを立てるくらいで、かなり大事なところに近づけるのではと思いました。
- この目的変数で、本当に見たいものを見ているのか
- この説明変数は、何を統制しようとしているのか
- この変数変換は、どんな仮説をコードに落としているのか
- この結果から、言ってよいことと言いすぎなことは何か
「コードを読めるか」というより、「このコードがやろうとしていることを説明できるか」に近いです。
動画尺をどうモデルに入れるか
たとえば、動画尺の処理はこんな感じです。
df["duration_min"] = df["duration_seconds"] / 60
df["duration_min_c"] = df["duration_min"] - df["duration_min"].mean()
df["duration_min_c2"] = df["duration_min_c"] ** 2
ここでは、動画尺をそのまま入れるのではなく、平均中心化した duration_min_c と、その二乗項 duration_min_c2 を作っています。
これは単なる前処理というより、「動画尺と再生パフォーマンスの関係は直線だけではないかもしれない」という仮説をコードに入れている部分です。
長ければ長いほどよい、という単純な関係ではなく、短すぎても長すぎても違うかもしれない。そのような曲線的な関係を、二乗項で表そうとしている。
このあたりは、コードとしては数行ですが、分析設計としてはかなり意味があります。
.fit()、え、これだけっすか
回帰分析の中心は、次の1行です。
result = smf.ols(formula=formula, data=df).fit()
正直、最初に見ると「え、これだけっすか」となります。
回帰分析というと、もっと行列を作ったり、係数を計算したり、いろいろ書くのかと思いきや、実際のコードは .fit() で終わっています。
もちろん、短いこと自体は悪くありません。ライブラリを使うというのは、そういうことです。
ただ、ここを「なんか回帰してくれるやつ」として済ませると、かなりブラックボックスになります。
たとえば、モデル式はこう書かれています。
log_views_per_day ~ log_subscriber_count
これは、左側が目的変数、右側が説明変数です。
ざっくり言えば、次のようなモデルです。
log_views_per_day = 切片 + 係数 * log_subscriber_count + 誤差
statsmodels の formula API は、この文字列を読み取り、内部で目的変数ベクトルと説明変数行列に変換してくれます。切片も自動で入ります。
ここまでは、なんとか追えます。
ただし、ここから先の「標準誤差がどう計算されるか」は、もう少し意識して見ないと流してしまいそうです。
モデルを段階的に積み上げる
このコードでは、1つのモデルだけでなく、A0からA3まで複数のモデルを推定しています。
formulas = {
"A0_subscriber_only": "log_views_per_day ~ log_subscriber_count",
"A1_duration": (
"log_views_per_day ~ log_subscriber_count + "
"duration_min_c + duration_min_c2"
),
"A2_posting_controls": (
"log_views_per_day ~ duration_min_c + duration_min_c2 + "
"log_subscriber_count + log_channel_video_count_in_sample + "
"days_since_publish + "
"C(published_hour_jst) + C(published_weekday_jst)"
),
"A3_genre": (
"log_views_per_day ~ duration_min_c + duration_min_c2 + "
"log_subscriber_count + log_channel_video_count_in_sample + "
"days_since_publish + "
"C(published_hour_jst) + C(published_weekday_jst) + C(genre)"
),
}
これはかなり大事そうです。というか、ここを見ないと、何を比較しているのか分からなくなります。
いきなり全部の変数を入れるのではなく、説明変数を段階的に増やしています。
| モデル | 見ていること |
|---|---|
| A0 | 登録者数だけでどの程度説明できるか |
| A1 | 動画尺を足すと説明力が増えるか |
| A2 | 公開時刻・曜日・公開後日数を足すとどうなるか |
| A3 | ジャンル差を統制するとどうなるか |
モデルを複雑にしたときに、何が変わったのかを追いやすくするための設計です。
これも、コードとしては辞書にformulaを書いているだけですが、分析としては「問いを段階的に深める」構造になっています。
ロバスト標準誤差って何やねん
一番引っかかったのがここです。
result = result.get_robustcov_results(cov_type="HC3")
ols(...).fit() までは、重回帰分析の基本としてまだ追いやすいです。しかし、そのあとに get_robustcov_results(cov_type="HC3") が出てくると、急に分からなくなります。
最初の感想は、かなり素朴に「ロバスト標準誤差って何やねん」でした。
しかも、これを自分で指定した記憶がない。
AIに「重回帰分析して」と頼んだら、勝手に HC3 が入っていた。これは、いい感じに気を利かせてくれたとも言えますが、別の見方をすれば、分析上かなり大事な設定を勝手に足しているとも言えます。
ここで、「なるほど、AIは賢いな」で済ませるのは少し怖いです。いや、怖いというより、普通に待ってくれという感じです。
調べた範囲では、普通のOLS標準誤差とロバスト標準誤差の違いは、ざっくりこうです。
| 項目 | 普通のOLS標準誤差 | ロバスト標準誤差 |
|---|---|---|
| 係数の推定値 | 同じ | 同じ |
| 標準誤差 | 誤差分散が一定という前提で計算 | 分散不均一があっても崩れにくいように計算 |
| t値・p値 | 普通の標準誤差に基づく | ロバスト標準誤差に基づく |
| 信頼区間 | 普通の標準誤差に基づく | ロバスト標準誤差に基づく |
図にすると、イメージはこんな感じです。
これはYouTubeの実データではありません。あくまで、誤差分散が一定の場合と、分散不均一がある場合の違いを見やすくするために、Pythonで適当に作ったサンプルデータです。
左は、誤差のばらつきがだいたい一定のケースです。普通のOLS標準誤差とHC3ロバスト標準誤差はかなり近い値になります。
右は、xが大きくなるほど点の散らばりも大きくなるケースです。こういう分散不均一があると、普通のOLS標準誤差とロバスト標準誤差の差が出やすくなります。
今回作った例では、傾きの推定値そのものは変わりません。一方で、傾きの標準誤差は次のように変わりました。
| ケース | 傾きの推定値 | 普通のOLS標準誤差 | HC3ロバスト標準誤差 |
|---|---|---|---|
| 誤差分散がほぼ一定 | 0.758 | 0.031 | 0.031 |
| 分散不均一あり | 0.733 | 0.046 | 0.055 |
ここで扱っているのは、厳密には「標準偏差」ではなく「標準誤差」です。標準偏差はデータそのもののばらつき、標準誤差は推定された係数の不確実性を表すものです。
ここで大事そうなのは、「ロバストと書いてあるから安心」としないことです。というか、そんなことできるわけないじゃないですか。
調べた範囲では、HC3 は係数を変える処理ではなさそうです。変わるのは標準誤差やp値の見積もりで、チャンネル内相関までは直接扱っていない。ひとまず、このくらいまで分かれば、結果の読み方は少し変わりそうです。
では、HC3 は数理的にどう計算されているのか。
ここまで来ると、正直、いったんブラックボックスにさせてください、という気持ちになります。というか、明らかにやばい方向に進んでいます。
ここから根掘り葉掘り進めると、残差、レバレッジ、分散共分散行列の補正、みたいな話がどんどん出てくるはずです。さすがに今回はそこまでは追いません。
ただ、今回のYouTubeデータについては、ロバスト標準誤差を使う意図自体は合っていそうだと思いました。
今回のデータに即して考えると、理由はこうです。
普通のOLS標準誤差は、ざっくり言えば「説明変数の値が違っても、誤差のばらつき方はだいたい同じ」という前提に乗っています。
でも、YouTubeのデータでそれを仮定するのは、少し無理がありそうです。
登録者数が小さいチャンネルでは、再生数のレンジも比較的狭いかもしれません。一方で、登録者数が大きいチャンネルでは、同じチャンネル内でも動画ごとの当たり外れがかなり大きくなりそうです。さらに、ジャンルによっても再生数の分布や外れ値の出方は違うはずです。
つまり、モデルが予測しきれなかったズレの大きさが、どの動画でも同じくらい、とは考えにくい。
そう考えると、「誤差分散が一定」という前提で標準誤差を計算するよりも、分散不均一があるかもしれない前提で、標準誤差を少し慎重に見積もる方が自然です。
もちろん、ロバスト標準誤差を使えばすべて解決、という話ではありません。チャンネル内の相関や、そもそものモデル指定の問題までは別に考える必要があります。
ただ、少なくとも今回のような、チャンネル規模もジャンルも混ざったYouTube動画データに対して、「普通のOLS標準誤差だけをそのまま信じるより、ロバスト標準誤差で確認する」という判断は、かなり妥当に見えます。
なので、ここでは「内部計算式までは追わない。ただし、今回のデータのばらつき方を考えると、標準誤差を少し慎重に見るための設定としては納得できる」くらいで止めておきます。
まとめ
やってみて、やっぱり本当にいい時代になったなと思いました。
理論さえある程度知っていれば、AIを駆動させることで、ものすごい勢いで分析コードを書けます。前処理をして、モデルを組んで、表を出して、図まで作る。そこまで一気に進むと、ある程度は「やり遂げた」という実感も得られます。
これを一から実装のノウハウとして学んでいたら、一つのことをやるだけでも何十時間とかかっていたはずです。pandas の書き方、statsmodels の使い方、結果表の整形、グラフの作り方。どれも、それ単体で普通に重い。
その意味では、AIが実装部分をかなり肩代わりしてくれるのは、かなり幸せなことだと思います。
一方で、確認するためのコストはむしろ増えています。
自分が理解できない高度なものをAIが勝手に実装してくると、その分だけ、あとから確認しなければならない理論が増えます。今回で言えば、HC3 ロバスト標準誤差がまさにそれでした。
頼んでいないのに、それっぽく高度な設定が入っている。便利ではある。でも、それを読む側は、「これは何をしているのか」「自分の分析に必要なのか」「どこまで信じてよいのか」を確認しなければならない。
つまり、実装コストは下がったけれど、確認コストは残る。場合によっては増える。
それでも、実装の苦労とどちらが楽なのかと言われれば、今のやり方の方が幸せだと感じます。
実際に作ってみることで、理論の学習も進みます。コードが先に動いてしまうからこそ、「これは何をしているのか」と後から調べる入口ができます。
もちろん、どこまで妥協なく確認するかは、自分との戦いです。趣味ならどこまでも掘れますが、仕事ならどこかで線を引かないといけない。その線をどこに引くかも含めて、たぶん分析の面白さなのだと思います。
生成されたPythonコードを、完璧ではなくても、自分の言葉で説明できるところまで読む。そして、「このコードは、自分がやりたかったことに合っているか」と問い直す。
そのくらいの姿勢でAIと付き合いながら、これからも勉強を続けていこうかなと思いました。