【ベーシックDay２】気象データアナリスト講座【データミックス】

本記事では、(株)データミックスの「気象データアナリスト養成講座」を受けた感想や、講義メモをまとめます。

今回は、2023年3月25日に受けた、ベーシックステップ２日目の内容です。

ロジスティック回帰

講義９日目の内容は「ロジスティック回帰」でした。

講義８日目では「線形回帰」を扱いましたが、比べると以下のような違いがあります。

例えば「試験の合格/不合格」や「広告をクリックする/しない」を目的変数とするのがロジスティック回帰です。

ロジスティック回帰の曲線を関数であらわすために、線形回帰の式を変換していきます。

具体的には、「－∞から∞」ではなく「0から1」の範囲におさまるような関数を考えていきます。

ロジスティック回帰の式は上記の通りです。

線形回帰の式をもとにして、「0から1」の範囲におさまるような関数になります。

「0から1」の範囲におさまっているかどうかは、以下の式を考えると確認できます。

交差エントロピー誤差と最尤推定

線形回帰の「最小二乗法」や「残差の二乗和」と同じように、ロジスティック回帰にも損失関数や計算方法があります。

交差エントロピー誤差関数は、たとえば「合格の人には高い推定値を出す/ 不合格の人には低い推定値を出す」を定式化したものです。

交差エントロピー誤差を最小化するようにモデルのパラメータを求める方法が最尤推定です。

オッズ比

「ある事象の起こりやすさ」を表す指標として、オッズがあります。

2つのオッズの比をオッズ比といいます。

例えば、「喫煙の有無」と「肺がんのリスク」に関するデータがあったとき、オッズとオッズ比は以下のようになります。

上記のデータでは、リスク比は1.25、オッズ比は1.5になりました。

「(a) << (b)」かつ「(c) << (d)」のとき、オッズ比とリスク比はほぼ等しくなります。

オッズ比は「e^β」であらわすことができます。

ロジスティック回帰のモデル式を変形して、オッズ比との関係を考えてみます。

ロジスティック回帰の式を変形していくと、オッズ比として解釈可能になります。

Ｒ言語でロジスティック回帰

ここからは、Rを使ってロジスティック回帰分析をやってみます。

実際の講義では、PimaIndiansDiabetes（ピマインディアン女性の糖尿病データ）のデータを使って進められました。

今回は、Rのデータセットにある「birthwt」のデータを使って分析します。

データのセットアップ

はじめに「birthwt」のデータを読み込み、どんなデータになっているのか見ていきます。

「birthwt」は「低出生体重児に関連する危険因子」のデータです。

library(MASS)
data(birthwt) ### データセットを読み込む
head(birthwt) ### 見出しを確認
View(birthwt) ### 表が出力される
birthwt %>% glimpse ### 横方向にデータの一覧を見れる

「View(birthwt) 」の結果は以下の通りです。

今回、目的変数(y)とするのは「low：出生時体重が2.5kg未満の指標（2.5kg未満は1）」のデータです。

それぞれの説明変数（age～bwt）がどのくらい説明できるのか、ロジスティクス回帰モデルを作ってみます。

EDA（要約統計量）～前処理

まずは「birthwt」の要約統計量などをチェックして、データの前処理をおこないます。

table(birthwt$low)
birthwt %>% summary

「low」には、0のデータが130個、1のデータが59個あることがわかりました。

「summary」の結果は以下の通りです。

ロジスティクス回帰では、欠損値（Naデータ）があるとエラーになります。

そのため欠損値があるときは前処理が必要です。

以下のコードで欠損値の処理ができます。

df_na_omit <- na.omit(birthwt)
df_na_omit %>% summary

今回のデータを見てみると、「ptl：早産経験の有無（経験有りが1）」に「2」と「3」のデータがあります。

ptlは「0または1」のデータのはずなので、「2」と「3」のデータは削除します。

## ptlが1より大きいデータを除外する
birthwt_sub <- birthwt %>% filter(ptl <= 1)

ロジスティクス回帰

ロジスティクス回帰モデルを作ります。

試しにすべての変数を使ってモデルを作ってみました。

model_logistic <- glm(low ~., data = birthwt_sub, family = binomial)

このロジスティクス回帰モデルを実行すると、以下のエラーになりました。

Warning messages:
1:  glm.fit: アルゴリズムは収束しませんでした  
2:  glm.fit: 数値的に 0 か 1 である確率が生じました  

どうやら「bwt：出生体重（単位はグラム）」が含まれていると、エラーになるようです。

なので「bwt」を除いて再度ロジスティクス回帰モデルを作っていきます。

birthwt_result <- birthwt_sub[,c(1,2,3,4,5,6,7,8,9)]
model_logistic <- glm(low~., data = birthwt_result, family = binomial)
model_logistic %>% summary

結果は以下の通りです。

線形回帰モデルと似たような出力結果となりました。

ptl, ht, smokeのp値は、有意水準0.05以下で有意とされます。

予測結果の確率

ロジスティクス回帰モデルは0から1の間の確率で出力されるので、その値を具体的に見てみます。

##予測結果の確率(0-1の間)
pred_proba <- predict(model_logistic, type = "response")
pred_proba %>% head

各データに対して、「出生時体重が2.5kg未満：1」となる確率が返ってきました。

「行番号85のデータ：22％」「行番号86のデータ：12％」...であるとわかります。

この結果をグラフ化して見てみます。

## 予測結果(確率)をdata frameに追加
df_pred <- birthwt_sub %>% 
  mutate(pred_proba = predict(model_logistic, type = "response") )

## low変数を因子型に変換
df_pred$low <- factor(df_pred$low, levels = c(0, 1), labels = c("2.5kg以上", "2.5kg未満"))

## 予測結果のヒストグラム (予測結果の確率 w. 実際の色)
df_pred %>% 
  ggplot(aes(x=pred_proba, colour=low))+
  geom_histogram(aes(y=..density.., fill=low), alpha=.3, position="identity") + geom_density(aes(colour = low), size =1)+
  geom_vline(xintercept=0.5, linetype="dashed")

グラフを見ると、赤と青が重なっているところが多くあります。

本当は「左に赤いグラフ、右に青いグラフ」となるのが理想ですが、今回の結果ではハッキリ分かれませんでした。

特に青いグラフが左側にも多くあります。

「出生時体重が2.5kg未満」であるにもかかわらず、予測の確率が50%を切るデータが多いです。

今後、予測モデルの改善や、他に使えるデータがないか探す必要があると思います。

相関と因果

講義の後半は「相関と因果」でした。

回帰モデルなど、あらゆるモデルは「バイアスがない理想の環境」を想定してつくられます。

しかし現実では、何らかのバイアスが入り込んでしまいます。

なるべくバイアスを減らして、理想環境に近づけることが必要です。

バイアスの種類

相関や因果に発生するバイアスは、大きく４種類あります。

それぞれのバイアスについて見ていきます。

本当は関連性がないのに、偶然にも相関関係が⾒られることがあります。

たとえば「気温が年々上がっている」と「海賊が年々減っている」というのは、グラフにすると相関関係があるように見えてしまいます。

しかし「気温が上がっているから海賊が減っている」わけではありません。

これが偶然のバイアスです。

因果の⽅向性が反対になっていることがあります。

たとえば、芸能人がある商品を紹介しているとき、本当は「商品が良いから紹介している」はずなのに、「芸能人か紹介しているから商品が良く見える」ことがあります。

これが反対のバイアスです。

事象Ａと事象Ｂに因果関係があるように見えて、じつは「別の要因X」が存在していることがあります。

このとき「別の要因X」を交絡因子といいます。

たとえば「アイスの売上」と「プールの事故者数」は、グラフにすると相関があるように見えます。

しかし実際には「気温」という交絡因子があり、「気温」と「アイスの売上」、「気温」と「プールの事故者数」で、それぞれ相関関係があります。

交絡因子によるバイアスがあるときは、交絡因子を固定することで、本来の関係性が見えてきます。

サンプルデータの抽出を誤ってしまうことで、因果関係があるように見えてしまうことがあります。

たとえば、お店の売上を増やしたくてクーポンメールを送った場合を考えます。

上記の場合、「クーポンメールのおかげで売上が増えた」と考えてしまいます。

しかし実際には、クーポンメールを送った客は「ヘビーユーザー層」かもしれません。

すると「クーポンメールのおかげで売上が増えた」とは言えなくなります。

このように、⺟集団の代表的なデータを選択できていないのが「選択バイアス」です。

バイアスへの対処

バイアスをなるべく減らすためには、「ドメイン知識（その業界に特化した知見）」と「データ」の両方から判断していく必要があります。

ドメイン知識として仮説を立てるには、ひとつひとつの変数とそれぞれの関係性を描きだすのが有効です。

相関・因果の検証

相関・因果の検証をおこなうためには、A/Bテストなどの実験をおこなう方法と、手元にあるデータから因果関係を探る方法（疑似実験）があります。

それぞれの検証では、妥当性のチェックが必要です。

実験

因果関係を解釈する方法として、実験があります。

実験するには、コストがかかったり、過去に起きたことだから実験できない状況のときもありますが、可能なら実験したほうが確実なデータが得られます。

たとえば「広告を見せた人と見せなかった人で、売上はどのくらい変わるのか？」という検証を考えましょう。

ポイントとなるのは、各グループをランダムに抽出することです。

年齢や職業、趣味、家族構成、運動量など、結果に影響を与えそうな要素に偏りがないように抽出します。

そうすると、純粋に広告の効果を測ることができます。

ただし、できるだけ選択バイアスがないようにすると、サンプルが小さくなってしまう可能性があるので、その点は注意が必要です。

疑似実験

データを集めたあとに、擬似的に実験効果を把握するのが「疑似実験」です。

疑似実験には、おもに３つの手法があります。

交絡因子があるのは、「アイスの売上」と「プールの事故者数」の関係のように疑似相関がある場合です。

気温を冬と夏で分けて考えるのと同じように、交絡因子ごとにデータを分割することで各データを比較します。

交絡因子を変数としてモデルに入れることで、疑似実験する方法もあります。

たとえば「広告を見せるグループと見せないグループ」のモデリングであれば、時間帯ごと（0～8時、8時～16時、16時～24時など）に分けて、時間帯を変数としてモデルに組み込みます。

「介入された人達」と「介入されなかった人達」の特徴を使って傾向スコアを予測し、「傾向スコアが同じ人達は似た人たち」とみなして、差分を見ていく方法です。

作業としては、「傾向スコアを予測するモデルを作る」ことになります。

予測する対象は「介入する確率」です。

傾向スコアを予測するモデルは、ロジスティック回帰で作ることができます。

傾向スコアは、必要な情報をすべて取れているとき（全交絡因⼦が観測済みの場合）のみ使える手法です。

DAGを描いてみる

DAG(Directed Acyclic Graph)とは、因果関係を可視化してバイアスの特定や処置に活用するツールです。

⽮印と形を使って、因果関係をあらわしていきます。

交絡因子があるときは、交絡因子を加味してモデルに入れます。

因果関係（A→Y）と外因Xの関係性

ここまで交絡因子について見てきましたが、外因Xには他の因子もあります。

因果関係（A→Y）に対して、外因Xが影響しているかどうかを考えます。

このとき、バックドア基準という考え方があります。

現状が「Open/Close」のどちらになっているかを確認し、Closeの状態にしていきます。

取得できたデータが同じでも、「何を知りたいか」によって、中間因子/交絡因子/共通因子は変わります。

AとYの因果関係を考えるとき、Aは複数ある因子Xのうち”注目した因子”に過ぎないので、人間の関心に合わせて「どの因⼦同⼠が独⽴なのか」を考慮する必要があります。

まとめ：ロジスティック回帰、相関と因果

今回は、(株)データミックス「気象データアナリスト養成講座」の、ベーシックステップ２日目の内容をまとめました。

今回の講義では、ロジスティック回帰の具体的な方法を学びました。

来週の講義もがんばりたいと思います。