多変量解析_05_回帰分析のおさらい
4.2K Views
July 25, 23
スライド概要
神戸大学経営学研究科で2022年度より担当している「統計的方法論特殊研究（多変量解析）」の講義資料「05_回帰分析のおさらい」です。
2026/05/03: ロジスティック回帰の説明を少し修正，GLMにおける誤差分散の表現についての説明を修正，ロジスティック回帰の回帰係数の解釈の説明を追加，回帰は因果ではない話を追加，タイポを修正しました。
2025/02/05: 全体的に内容に手を入れました。RmarkdownからQuartoに変更しました。
神戸大学経営学研究科准教授 分寺杏介（ぶんじ・きょうすけ）です。 主に心理学的な測定・教育測定に関する研究を行っています。 講義資料や学会発表のスライドを公開していきます。 ※スライドに誤りを見つけた方は，炎上させずにこっそりお伝えいただけると幸いです。
関連スライド
各ページのテキスト
5.1 回帰分析とは 1 Chapter 5 回帰分析のおさらい 各種分析手法（因子分析・構造方程式モデリング・項目反応理論）の紹介に向けた準備として，各手 法のベースにある回帰分析の理論的な部分をおさらいしつつ，いくつかの重要な概念を説明して います。 本資料は，神戸大学経営学研究科で 2022 年度より担当している「統計的方法論特殊研究（多変 量解析）」および京都大学教育学部で 2025 年度より担当している「心理・教育測定演習」の講義資 料です。CC BY-NC 4.0 ライセンスの下に提供されています。 作成者連絡先 神戸大学大学院経営学研究科 分寺杏介 （ぶんじ・きょうすけ） mail: [email protected] website: https://www2.kobe-u.ac.jp/~bunji/ 統計を学ぶ人の大多数が最初に出会う「多変量解析」は，きっと（重）回帰分析でしょう。実際に， 後半で登場する因子分析・構造方程式モデリング・項目反応理論はいずれも回帰分析のモデルと 関係のあるモデルです。一応本講義では回帰分析については既習であることを前提としているの で，基本的な回帰分析のおさらいをはさみ，後半の内容に必要な部分だけでも前提知識を補いたい と思います。なお，この Chapter では，基本的に今後紹介する手法の理解に必要・重要と思われる 内容に絞って紹介しています。回帰分析は相当幅広く用いられているため，本気を出せば回帰分 析だけで 1 コマ 15 回の講義を余裕で組めると思います。回帰分析を深く理解したい方は，恐縮な がら経済学研究科や国際協力研究科などで開講されている講義も探してみてください。 5.1 回帰分析とは 中学数学では，図 5.1 のように「2 つの点を通る直線を求める」問題を解いたことがあると思い ます。 回帰分析はこの延長にある分析手法と言っても良いでしょう。
5.1 回帰分析とは 2 図 5.1: 2 つの点を通る直線を求める この問題を （今後の問題設定に合わせて） 数式を用いて表現すると， 一次関数 𝑦 = 𝑏0 + 𝑏1 𝑥 が (𝑥𝑝 , 𝑦𝑝 ) = (−1, −3), (1, 1) の 2 点を通るとする。 このとき， 係数 𝑏0 , 𝑏1 の値はそれぞれいくつか。 と言った感じでしょうか。なお，これ以降は一つ一つのデータ（行方向）を表すために，回答者 (person) の頭文字をとって添字 𝑝 をつけることにします。つまり，(𝑥𝑝 , 𝑦𝑝 ) = (−1, −3), (1, 1) は 丁寧に書くと (𝑥1 , 𝑦1 ) = (−1, −3), (𝑥2 , 𝑦2 ) = (1, 1) となり，それぞれ 1 人目，2 人目の変数 𝑥, 𝑦 の値を表しているとします。 さて， 先程の書き方をもう少し変えてみると， (𝑥𝑝 , 𝑦𝑝 ) = (−1, −3), (1, 1) の 2 点に関する以下の連立方程式を満たす係数 𝑏0 , 𝑏1 の値はそ れぞれいくつか。 { −3 = 𝑏0 + 𝑏1 × −1 1 = 𝑏 0 + 𝑏1 × 1 と表すこともできるでしょう。そして 図 5.1 の場合，正解は (𝑏0 , 𝑏1 ) = (−1, 2) となり，一次関 数で表せば 𝑦 = −1 + 2𝑥 となります。 ここまでは 2 つの点のみについて考えてきましたが，点の数が増えたらどうでしょうか。図 5.2 のように， 「3 つの点を通る直線を求める」問題があったとしても，すべての点が一直線上に並ん でいるなら一つの答えを出すことができます。 このケースでも， 正解は 𝑦 = −1 + 2𝑥 ですね。 では，図 5.3 のように 3 つの点が一直線上に並んでいない場合は？当然ながらすべての点を通 る直線は存在しません。回帰分析の問題設定はこういった状況です。すべての点を通らないならば それなりに最もそれっぽい直線を何らかの方法で決定して引いてあげる必要があります。 どんな線が最もそれっぽい直線かですが，図 5.3 の赤い線と青い線を比較すると，ほとんどの 人は赤い線のほうがそれっぽいと感じるのではないでしょうか。青い線の方は，3 つのうち 2 つ の点を通過しているにも関わらず，なぜか赤い線のほうがそれっぽいですね。これは，赤い線のほ うが平均的にはズレが小さいためです。 確かに青い線は， 通過している 2 つの点に関してはズレが ゼロですが，残り 1 個の点に対しては大きくズレています。この赤い線のように，平均的なズレを
5.1 回帰分析とは 3 図 5.2: 3 つの点を通る直線を求める 図 5.3: 3 つの点を通る直線なんてない 最小化する直線を求める分析法が回帰分析なのです。 点の数が何個であろうと最もベーシックな回帰分析では，直線（一次関数）を求めます。傾きを 𝑏1 ，切片を 𝑏0 とおけば 𝑦 = 𝑏 0 + 𝑏1 𝑥 (5.1) という直線について，𝑏0 , 𝑏1 の値がいくつのときに最もそれっぽい直線になるかを求めるわけで す。 5.1.1 最小二乗法 実際の分析においては，最小二乗法によって回帰直線を求めることが多いです。これは図 5.4 の青い矢印のように，縦方向でのズレを求めて，これの二乗和が最小になる (𝑏1 , 𝑏0 ) の組み合わ せを求める方法です。 「それっぽい直線」という観点で言えば，図 5.5 の青い矢印のように，赤い直線との直線距離に 関して最小化しても良さそうな気がしますが，最小二乗法では縦方向での距離を最小化していま す。 なぜなのでしょうか。 それは，回帰分析が 𝑥 で 𝑦 を説明する分析手法だからです。𝑦 そのものは無いけれど他の変数 がある時に 𝑦 を近似するのが回帰直線なのです。例えば，𝑥 を「体重」 ，𝑦 を「身長」だとしましょ う。あなたのもとに，あるゲストがやってきます。どういうわけかその人に関して「体重」の情報だ
5.1 回帰分析とは 4 図 5.4: 最小二乗法 図 5.5: 直線距離で考えたらダメなのかい けが入っているとします。あなたはゲストのためにちょうどよい高さの椅子を用意する必要があ ります。しかし，ちょうどよい高さの椅子をあてがうためには，ゲストの身長が分かっていないと いけません*1 。 さて， あなたはゲストの身長をどのように見積もるでしょうか。 普通ならば，図 5.6 のように，周りの人の体型などから「 （標準体型で）体重が 70kg くらいの人 ならば，身長はだいたい 170cm くらい」 「 （標準体型で）体重が 60kg くらいの人ならば，身長はだ いたい 162cm くらい」 という感じで予測を立てたうえで，ゲストの体重の情報を当てはめて 「ゲス トの体重が 65kg ということは，標準体型ならば身長はだいたい 165.5m くらいかなぁ」といった 予測を行うでしょう。回帰分析で行っていることは概ねこんな感じのことですが，ここで問題にな る「予測のズレ」は，身長（つまり 𝑦 ）に関してのみです。つまり，𝑥 による 𝑦 の近似（予測）のズレが 最小になるように回帰直線を引きたいがために，最小二乗法では縦方向でのズレのみを考えてい るわけですね。 求めたい回帰直線は 𝑦 = 𝑏0 + 𝑏1 𝑥 ですが，一つ一つの点について考えると，ほとんどの場合等 号は成立しません。というのも回帰直線が示すのはあくまでも予測値であって，実際のデータでは 予測とのズレがあるためです。 そこで回帰分析モデルでは， 実際の一つ一つのデータの値 𝑦𝑝 は， 回 *1 もちろん実際には股下やら体型やらによって最適な椅子の高さは決まるわけですが， 手元には「体重」の情報しか無 いので，できる範囲で予測しましょう，ということです。
5.1 回帰分析とは 5 図 5.6: 回帰直線による予測 帰直線による予測値 𝑦𝑝̂ と誤差 𝜀𝑝 の和であると考えます。 ということで， 回帰分析モデルの全体は 𝑦𝑝 = 𝑦𝑝̂ + 𝜀𝑝 (5.2) 𝑦𝑝̂ = 𝑏0 + 𝑏1 𝑥𝑝 𝜀𝑝 ∼ 𝑁 (0, 𝜎𝜀2 ) となり，𝑃 をデータの総数とすると， 最小二乗法では 𝑃 𝑃 ∑(𝑦𝑝 − 𝑦𝑝̂ )2 = ∑ [𝑦𝑝 − (𝑏0 + 𝑏1 𝑥𝑝 )] 𝑝=1 2 (5.3) 𝑝=1 を最小化するような (𝑏0 , 𝑏1 ) の組み合わせを求めている，というわけです。 INFO 最尤法 最小二乗法と同じくらい用いられる推定法が最尤法です。(5.2) 式の回帰分析モデルの式を 変形させると， 𝑦𝑝 ∼ 𝑁 (𝑏0 + 𝑏1 𝑥𝑝 , 𝜎𝜀2 ) (5.4) という形にもなります。 したがって， データから最も尤度が大きくなる (𝑏0 , 𝑏1 , 𝜎𝜀 ) の組を見 つけてあげることで， 回帰直線を最尤推定することも可能なわけです。 最小二乗法とどちらが良いかという話ではなく，最適化の目的関数が異なるということな ので，回帰直線を求める方法は最小二乗法だけじゃないということは頭の片隅に入れてお いても良い知識だと思います。 5.1.2 （単）回帰分析の線形代数的表現 因子分析や SEM のところでは，それぞれの分析手法の基本的な考え方を説明するために線形 代数的表現を多用しています。ということで，ここで回帰分析の線形代数による表現を説明してお きます。 書き方とその意味に慣れておいてください。 まず，回帰分析モデルを，図 5.1 のところで示した連立方程式的な書き方で書き直すと，次のよ うになります。
5.2 重回帰分析 6 (𝑥1 , 𝑦1 ), (𝑥2 , 𝑦2 ), ⋯ , (𝑥𝑃 −1 , 𝑦𝑃 −1 ), (𝑥𝑃 , 𝑦𝑃 ) の 𝑃 個の点に関する以下の連立方程式につい て， 誤差の和 ∑𝑝=1 𝜀2𝑝 が最小になる係数 𝑏0 , 𝑏1 の値はそれぞれいくつか。 𝑃 ⎧ 𝑦1 = 𝑏0 + 𝑏1 𝑥1 + 𝜀1 { { 𝑦2 = 𝑏0 + 𝑏1 𝑥2 + 𝜀2 { ⋮ ⎨ {𝑦𝑃 −1 = 𝑏0 + 𝑏1 𝑥𝑃 −1 + 𝜀𝑃 −1 { { ⎩ 𝑦𝑃 = 𝑏0 + 𝑏1 𝑥𝑃 + 𝜀𝑃 毎回全部式を書いていると面倒なので，ここで線形代数の力を借ります。連立方程式のうち，す べての式に共通しているのが 𝑏0 , 𝑏1 の部分です。一旦ここはそのままにしておいて，その他の（𝑝 ごとに異なる） 部分をベクトルでまとめて書き直してみます。 1 𝑥1 𝜀1 ⎧ 𝑦1 {⎡ 𝑦 ⎤ ⎡1⎤ ⎡ 𝑥 ⎤ ⎡ 𝜀 ⎤ 2 {⎢ 2 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ 2 ⎥ ⋮ ⎥ = 𝑏 0 ⎢ ⋮ ⎥ + 𝑏1 ⎢ ⋮ ⎥ + ⎢ ⋮ ⎥ ⎢ ⎨ ⎢1⎥ ⎢𝑥𝑃 −1 ⎥ ⎢𝜀𝑃 −1 ⎥ {⎢𝑦𝑃 −1 ⎥ { ⎩ ⎣ 𝑦𝑃 ⎦ ⎣1⎦ ⎣ 𝑥𝑃 ⎦ ⎣ 𝜀𝑃 ⎦ (5.5) なお切片 𝑏0 には，ベクトル同士の足し算を成立させるために 𝑃 人分の 1 が並んだベクトル（以後 太字の 𝟏 と表記する） をかけています。 だいぶ見通しが良くなってきました。あとはベクトルを太字で表すことで，以下のように割と綺 麗に書き直すことが出来ます。 2 つの変数に関する長さ 𝑃 のベクトル 𝐱 = [𝑥1 𝑥2 ⋯ ⊤ 𝑥𝑃 ] , 𝐲 = [𝑦1 が与えられたとき，以下の式における誤差ベクトル 𝜺 = [𝜀1 𝜀2 ⋯ 𝑦2 ⋯ ⊤ 𝑦𝑃 ] ⊤ 𝜀𝑃 ] の二乗和 𝜺⊤ 𝜺 が最小になる係数 𝑏0 , 𝑏1 の値はそれぞれいくつか。 𝐲 = 𝑏 0 𝟏 + 𝑏1 𝐱 + 𝜺 ということで，これでデータが何個あっても一つの式で表すことができるようになりました。便 利ですね。 5.2 重回帰分析 𝑦 を予測するという意味では，もちろん変数は多い方が良いでしょう。 「身長」を予測するうえ で，体重だけでなく「性別」や「体脂肪率」といったことがわかれば，予測はより正確になるはずで す。回帰分析でも，説明変数を複数個用いることが出来ます。これを一般的には重回帰分析と呼び ます。説明変数が 𝐾 個あるならば，対応する回帰係数を 𝐾 + 1（切片の分）個用意して，以下のよ
5.2 重回帰分析 7 うな直線を考えるわけです*2 。 𝑦 = 𝑏 0 + 𝑏 1 𝑥1 + 𝑏 2 𝑥2 + ⋯ + 𝑏 𝑘 𝑥𝑘 + ⋯ + 𝑏 𝐾 𝑥𝐾 (5.6) ここでのポイントは，重回帰分析では変数の効果が線形和（重み付け和）になっているという点 です。回帰係数は「その説明変数の値が 1 大きくなったら 𝑦 の予測値がどれだけ大きくなるか」を 示していますが， 最もベーシックな重回帰分析モデルの場合 • 説明変数の値がいくつであってもその効果は一定 • 他の説明変数の値は完全に無関係 であることが仮定されています。もし高次の効果（e.g., 薬は適量ならば飲むほど健康に良いが， 飲み過ぎるとかえって良くない＝二次の効果）や交互作用を考慮したい場合には，直接それを表す 項を追加する必要があります。 例えば 𝑦 = 𝑏0 + 𝑥1 𝑏1 + 𝑥2 𝑏2 + 𝑥21 𝑏3 + 𝑥1 𝑥2 𝑏4 (5.7) といった感じです。 また，重回帰分析で得られる回帰係数は，厳密には偏回帰係数と呼ばれます。これは「他の説明変 数の値が固定されているときに，その説明変数の値が 1 大きくなったら 𝑦 の予測値がどれだけ大 きくなるか」を表しています。もう少し正確に言うと，その説明変数以外で 𝑦 を予測した際の誤差 に対して，その説明変数のみで単回帰分析をした際の回帰係数が偏回帰係数です。例えば 2 変数 𝑥1 , 𝑥2 の重回帰における 𝑥1 の偏回帰係数 𝑏1 は， 𝑦 = 𝑏01 + 𝑏2 𝑥2 + 𝜀1 (Step1) 𝜀1 = 𝑏02 + 𝑏1 𝑥1 + 𝜀2 (Step2) (5.8) という二段階で推定することで，単回帰の組み合わせに分解することができるのです（Step 1 の 左辺は 𝑦 ̂ ではなく 𝑦 である点に注意） 。 これにより，偏回帰係数は他の変数の影響を完全に取り除 いた （コントロールした） 上でのある変数の影響の大きさを表している， と解釈できるわけです。 5.2.1 重回帰分析の線形代数的表現 説明変数の数が増えても，同様に線形代数による表現を利用します。(5.6) 式に基づく予測を， まず 𝑝 番目のデータ一つについて考えてみます。説明変数が何個になっても，観測値 𝑦𝑝 を予測値 𝑦𝑝̂ と誤差 𝜀𝑝 に分ける点 (5.2 式) は変わりません。 𝑦𝑝 = 𝑦𝑝̂ + 𝜀𝑝 = 𝑏0 + 𝑏1 𝑥𝑝1 + 𝑏2 𝑥𝑝2 + ⋯ + 𝑏𝐾 𝑥𝑝𝐾 + 𝜀𝑝 (5.9) これ以降， 説明変数には 2 つの添字が付きます。2 つ目の添字 𝑘(𝑘 = 1, 2, ⋯ , 𝐾) は， 何番目の説明 変数であるかを表すものです。例えば 𝑥𝑝3 （または 𝑥𝑝,3 ）は「𝑝 番目の人の 3 番目の説明変数の値」 *2 直線がイメージしづらいかもしれませんが， 説明変数が 1 個の場合の回帰直線が 2 次元平面に引かれていたよう に，説明変数が 2 個の場合は 3 次元空間の中に一本の直線が引かれます。 同様に， 説明変数が 𝐾 個ならば 𝐾 + 1 次 元空間に一本の直線を引くわけです。
5.2 重回帰分析 8 （dat[p,3]） を指しています。 添字の順序は，R でデータフレームの行と列の要素を指定するとき の順番（apply() 関数の引数 MARGIN の指定） だと覚えると良いかもしれません （個人の感想） 。 まずは (5.9) 式の 𝑏1 𝑥𝑝1 + 𝑏2 𝑥𝑝2 + ⋯ + 𝑏𝐾 𝑥𝑝𝐾 （切片と誤差項以外）の部分について，ベク トルの積を使って簡略化していきます。これは 𝑝 番目のデータにおける説明変数ベクトル 𝐱𝑝 = [𝑥𝑝1 , 𝑥𝑝2 , ⋯ , 𝑥𝑝𝐾 ] と，回帰係数ベクトル 𝐛 = [𝑏1 , 𝑏2 , ⋯ , 𝑏𝐾 ]⊤ を使うと， 𝑏1 ⎡𝑏 ⎤ 𝐱𝑝 𝐛 = [𝑥𝑝1 𝑥𝑝2 ⋯ 𝑥𝑝𝐾 ] ⎢ 2 ⎥ ⎢ ⋮ ⎥ ⎣ 𝑏𝐾 ⎦ = 𝑥𝑝1 𝑏1 + 𝑥𝑝2 𝑏2 + ⋯ + 𝑥𝑝𝐾 𝑏𝐾 (5.10) と，かなりシンプルに表すことが出来ました。ここでもう少し工夫してみます。せっかくなの で，切片についても同じようにまとめてしまいましょう。 (5.9) 式の 𝑏0 の部分は， 𝑏0 × 1 と見 ることも可能です。 𝑏1 には 𝑥𝑝1 がかかっていたのと同じように考えるならば， 𝑏0 の部分は 𝑥𝑝0 = 1 ということです。このように考えて，説明変数ベクトルと回帰係数ベクトルをそれぞれ 𝐱𝑝 = [1, 𝑥𝑝1 , 𝑥𝑝2 , ⋯ , 𝑥𝑝𝐾 ]，𝐛 = [𝑏0 , 𝑏1 , 𝑏2 , ⋯ , 𝑏𝐾 ]⊤ と書き換えてやれば， 𝑏0 ⎡𝑏 ⎤ ⎢ 1⎥ 𝐱𝑝 𝐛 = [1 𝑥𝑝1 𝑥𝑝2 ⋯ 𝑥𝑝𝐾 ] ⎢ 𝑏2 ⎥ ⎢ ⋮ ⎥ ⎣𝑏𝐾 ⎦ = 𝑏0 + 𝑥𝑝1 𝑏1 + 𝑥𝑝2 𝑏2 + ⋯ + 𝑥𝑝𝐾 𝑏𝐾 (5.11) = 𝑦𝑝̂ というように， 切片項まで含めて回帰直線による予測値 𝑦𝑝̂ を簡略化できました。 あとはこれを 𝑃 個のデータ全体に拡張するだけです。 といっても考え方は簡単で，𝑃 × (𝐾 + 1) の説明変数全体 （+ 切片） を 𝐱1 1 ⎡ 𝐱 ⎤ ⎡1 2 𝐗=⎢ ⎥=⎢ ⎢ ⋮ ⎥ ⎢⋮ ⎣𝐱𝑃 ⎦ ⎣1 𝑥11 𝑥21 ⋮ 𝑥𝑃 1 𝑥12 𝑥22 ⋮ 𝑥𝑃 2 ⋯ ⋯ ⋱ ⋯ 𝑥1𝐾 𝑥2𝐾 ⎤ ⎥ ⋮ ⎥ 𝑥𝑃 𝐾 ⎦ (5.12)
5.2 重回帰分析 9 と表せば， これを用いてすべてのデータに関する予測値 𝐲̂ を 𝑦1̂ ⎡𝑦 ̂ ⎤ 𝐲̂ = ⎢ 2 ⎥ = 𝐗𝐛 ⎢⋮⎥ ⎣𝑦𝑝̂ ⎦ 𝑏 𝑥1𝐾 ⎡ 0 ⎤ 𝑏 𝑥2𝐾 ⎤ ⎢ 1 ⎥ ⎥ ⎢ 𝑏2 ⎥ ⋮ ⎥⎢ ⎥ ⋮ 𝑥𝑃 𝐾 ⎦ ⎣𝑏𝐾 ⎦ 𝑏0 + 𝑥11 𝑏1 + 𝑥12 𝑏2 + ⋯ + 𝑥1𝐾 𝑏𝐾 ⎡ 𝑏 +𝑥 𝑏 +𝑥 𝑏 +⋯+𝑥 𝑏 ⎤ 21 1 22 2 2𝐾 𝐾 ⎥ =⎢ 0 ⎢ ⎥ ⋮ ⎣𝑏0 + 𝑥𝑃 1 𝑏1 + 𝑥𝑃 2 𝑏2 + ⋯ + 𝑥𝑃 𝐾 𝑏𝐾 ⎦ 1 ⎡1 =⎢ ⎢⋮ ⎣1 𝑥11 𝑥21 ⋮ 𝑥𝑃 1 𝑥12 𝑥22 ⋮ 𝑥𝑃 2 ⋯ ⋯ ⋱ ⋯ (5.13) とまとめることが可能です。あとは誤差ベクトルを付け加えたら，以下のようにスッキリとした表 現の出来上がりです。 長さ 𝑃 の結果変数ベクトル 𝐲 と 𝑃 × (𝐾 + 1) の説明変数 （＋切片） 行列 𝐗 が与えられたと き，以下の式における誤差ベクトル 𝜺 の二乗和 𝜺⊤ 𝜺 が最小になる回帰係数ベクトル 𝐛 の値 はいくつか。 𝐲 = 𝐗𝐛 + 𝜺 因子分析や SEM のところでは， このような感じで 「変数の行列」×「係数のベクトル （または行 列） 」の積を用いて手法のコンセプトを紹介していきます。ということで，今のところはここで紹介 した記法に慣れておいてください。 LIGHTBULB 正規方程式 上記のように重回帰分析を線形代数で表現すると，最小二乗法の解も解析的に求めるこ とが出来ます。いま，求めたいのは 𝜺⊤ 𝜺 が最小になるような 𝐛 の値でした。そこで， 𝜺 を (𝐲 − 𝐗𝐛) に戻して考えると，𝜺⊤ 𝜺 = (𝐲 − 𝐗𝐛)⊤ (𝐲 − 𝐗𝐛) が最小になるような 𝐛 を求め たら良いわけです。 最小値になる 𝐛 を求めたいわけなので，ちょっと大変かもしれませんが，行列を展開して 𝐛 で偏微分したものイコール 0 という式をとれば*3 ， 𝐗⊤ 𝐗𝐛 = 𝐗⊤ 𝐲 というかたちになります。この方程式のことを正規方程式と呼びます。あとは 𝐗⊤ 𝐗 が正則 であれば， 両辺に逆行列をかけることで， この式を 𝐛 = の形に変形させ， 𝐛 = (𝐗⊤ 𝐗)−1 𝐗⊤ 𝐲 という解が得られます。
5.2 重回帰分析 10 ちなみに「𝐗⊤ 𝐗 が正則でない」ときには逆行列が求められないため，解を得ることが出来ま せん。このような状態になるのは，行列がランク落ちしている場合ですが，ランク落ちが発 生するのは 𝐗 の特定の変数が線形従属になっている場合，つまり他の変数の線形和で表せ てしまう場合です。回帰分析ではよく「多重共線性に気をつけましょう」という話がありま すが，この理由として，完全な多重共線性がある場合には 𝐗⊤ 𝐗 が正則でなくなってしま うために， 解が求められなくなってしまうから， という説明も可能なのです。 5.2.2 （おまけ）ダミー変数 回帰係数は「その説明変数の値が 1 大きくなったら 𝑦 の予測値がどれだけ大きくなるか」を示 しています。では，説明変数が名義変数の場合はどうしたら良いでしょうか。まずはカテゴリ数が 2 の場合（e.g., 生物学的な性：男女）を考えてみます。通常，このようなデータでは一方のカテゴリ に 1 を，もう一方のカテゴリに 0 をあてます。今回は「男性は 0，女性は 1」としておきましょう。こ の 「性別 (𝑥𝑝2 )」 と 「体重 (𝑥𝑝1 )」 から 「身長 (𝑦𝑝 )」 を予測する重回帰モデルを立てると (5.14) 𝑦𝑝̂ = 𝑥𝑝1 𝑏1 + 𝑥𝑝2 𝑏2 + 𝑏0 となります。このモデルを男性と女性それぞれに分けてみてみると，男性は 𝑥2 = 0，女性は 𝑥2 = 1 なので + 𝑏0 (男性) 𝑦𝑝̂ = 𝑥𝑝1 𝑏1 + 𝑏2 + 𝑏0 (女性) 𝑦𝑝̂ = 𝑥𝑝1 𝑏1 (5.15) となります。したがって，男性モデルでは切片が 𝑏0 となっており，女性モデルでは切片が 𝑏2 + 𝑏0 となっているわけです。回帰係数の解釈通り，𝑏2 は 𝑥2 の値が 1 大きくなった際に 𝑦 の予測値が どれだけ大きくなるかを意味するわけですが，より正確に言えば体重 𝑥1𝑝 の値が同じ男性と女性 がいた場合， 身長は平均的にどれだけ異なるかを表しているわけですね。 カテゴリ数が 3 の場合はこうはいきません。先程の例について，性別の代わりに血液型を 𝑥2 と してみます。血液型は A, B, O, AB をそれぞれ 1, 2, 3, 4 とコーディングしたものだとすると，重 回帰モデルは ⎧𝑥𝑝1 𝑏1 + 𝑏2 + 𝑏0 { {𝑥𝑝1 𝑏1 + 2𝑏2 + 𝑏0 𝑦𝑖̂ = ⎨𝑥𝑝1 𝑏1 + 3𝑏2 + 𝑏0 { { ⎩𝑥𝑝1 𝑏1 + 4𝑏2 + 𝑏0 (𝐴) (𝐵) (𝑂) (5.16) (𝐴𝐵) となります。 つまり A 型-B 型の差と B 型-O 型の差と O 型-AB 型の差がすべて同じと， 明らか におかしなことになっています。このように，名義変数についてそのまま対応可能なカテゴリ数は 2 までなので，カテゴリ数が 3 以上の場合には，そのまま数値化して投入するのではなく，二値の ダミー変数化することで対応します。 血液型の場合，4 カテゴリなので 3 つのダミー変数 (𝑥2𝐴 , 𝑥2𝐵 , 𝑥2𝑂 ) および対応する回帰係数 *3 このあたりの展開は 「正規方程式」で検索したらすぐ見つかります
5.2 重回帰分析 11 (𝑏2𝐴 , 𝑏2𝐵 , 𝑏2𝑂 ) を用意します。これらのダミー変数はそれぞれ • 𝑥2𝐴 は A 型の場合 1，それ以外の場合 0 になる • 𝑥2𝐵 は B 型の場合 1，それ以外の場合 0 になる • 𝑥2𝑂 は O 型の場合 1，それ以外の場合 0 になる と設定し，4 つの血液型に対してダミー変数はそれぞれ以下のような値を取ります。 血液型 𝑥2𝐴 𝑥2𝐵 𝑥2𝑂 A型 1 0 0 B型 0 1 0 O型 0 0 1 AB 型 0 0 0 こうすることで重回帰モデルは ⎧𝑥𝑝1 𝑏1 + 𝑏2𝐴 +𝑏0 { {𝑥𝑝1 𝑏1 + 𝑏2𝐵 +𝑏0 𝑦̂ = ⎨𝑥𝑝1 𝑏1 + 𝑏2𝑂 +𝑏0 { { +𝑏0 ⎩𝑥𝑝1 𝑏1 (𝐴) (𝐵) (𝑂) (5.17) (𝐴𝐵) と書くことができ，各回帰係数は AB 型との差として解釈可能になります。ここでのポイントは， ダミー変数はカテゴリ数 −1 個だという点です。なぜカテゴリ数ぶんのダミー変数があると困る のか考えてみましょう。例えば上の 3 つのダミー変数に加えて「AB 型なら 1」となる 𝑥2𝐴𝐵 およ び回帰係数 𝑏2𝐴𝐵 があるとします。最小二乗法で推定した結果，(𝑏2𝐴 , 𝑏2𝐵 , 𝑏2𝑂 , 𝑏2𝐴𝐵 ) がそれぞれ (1, 2, 3, 4)，また切片 𝑏0 が 0 だったとしましょう。これを重回帰モデルで表すと ⎧𝑥𝑝1 𝑏1 + 1+0 { {𝑥𝑝1 𝑏1 + 2+0 𝑦̂ = ⎨ {𝑥𝑝1 𝑏1 + 3+0 { ⎩𝑥𝑝1 𝑏1 + 4+0 (𝐴) (𝐵) (𝑂) (5.18) (𝐴𝐵) となります。ですが，この式には別の解が考えられます。例えば (𝑏2𝐴 , 𝑏2𝐵 , 𝑏2𝑂 , 𝑏2𝐴𝐵 ) がそれぞれ (−4, −3, 2, 1)，切片 𝑏0 が 5 だとしても，重回帰式は ⎧𝑥1𝑝 𝑏1 − 4+5 { {𝑥1𝑝 𝑏1 − 3+5 𝑦̂ = ⎨𝑥1𝑝 𝑏1 − 2+5 { { ⎩𝑥1𝑝 𝑏1 − 1+5 (𝐴) (𝐵) (𝑂) (5.19) (𝐴𝐵) となり，結果的に (5.18) 式と全く同じになってしまいました。このように，回帰分析ではダミー変 数がカテゴリ数の数だけあると解が一意に定まらなくなってしまうため，カテゴリ数-1 個だけ用
5.2 重回帰分析 意する必要があるのです*4 。 ダミー変数がカテゴリ数より少ないとなると，必然的にダミー変数が与えられない（全てのダミ ー変数が 0 であることによって表される） カテゴリが発生します。 このカテゴリのことを基準カテ ゴリや参照カテゴリと呼んだりしますが，この基準カテゴリの決め方も重要になります。というの も回帰係数は基準カテゴリとの差であるためです。回帰分析では，回帰係数の検定（𝐻0 ∶ 𝑏 = 0）が 行われるため， 「あるカテゴリが基準カテゴリと平均的に差があるか」は検証できる一方で，基準カ テゴリ以外のカテゴリ同士での差の検定はそのままでは出来ません。実験的な環境であれば対照 群を基準カテゴリにおくのが良いでしょう。 5.2.3 R で重回帰分析 ↓本 Chapter で使用するファイルのダウンロードはこちらから chapter04.rds 読み込み 1 dat <- readRDS("chapter04.rds") ここらで少し，R での重回帰分析のやり方を紹介しておきましょう。R では，lm() という関 数があります （linear model の略） 。 今回は試しに 「年齢 age」 と 「性別 gender」 の 2 つを説明変数と して， 第 1 因子の 5 項目の合計点 （Q1_A） を結果変数とします。 lm() を始めとする R のいくつかの関数では，説明変数と結果変数の関係を表すために，以前少 しだけ登場した特殊な記法 （formula） を使用します。 lm() で回帰分析 1 lm(Q1_A ~ age + gender, data = dat) 第一引数には回帰モデルの式 (formula) を与えます。上の例では Q1_A ~ age + gender が第 一引数です。見方としては重回帰モデル式 𝑦 = 𝑥1 𝑏1 + 𝑥2 𝑏2 と同じですが，両辺をチルダ (~) でつ なぐのが特徴です。とりあえずこの授業の中ではチルダで結んだら回帰ということを覚えておき ましょう。 実は SEM で使用する lavaan パッケージの記法でもチルダは回帰を表します。 INFO formula で使える記号 基本的には重回帰分析のように + で説明変数を足していけば良いのですが，他にもいくつ かの記号があります。 • 2 つの説明変数を: でつなぐと，その変数間の交互作用項を指定できます。例えば x1 と x2 があるときに，y ~ x1:x2 とすると，これは 𝑦 = 𝑥1 𝑥2 𝑏12 + 𝑏0 という回帰式 *4 ちなみに機械学習では one-hot encoding という名で使われるケースもあります。 推定できなくなってしまうの は，説明変数の効果が線形和で表される回帰モデルならではの問題なのです。 12
5.2 重回帰分析 13 を表します。 • 2 つの説明変数を*でつなぐと，その変数の主効果と交互作用項を同時に指定できま す。 例えば y ~ x1*x2 とすると， これは 𝑦 = 𝑥1 𝑏1 + 𝑥2 𝑏2 + 𝑥1 𝑥2 𝑏12 + 𝑏0 という回帰 式を表します。 • 0 を足すと，切片項 𝑏0 が 0 に固定されます。原点を通る回帰をしたい場合に使えるか もしれません。例えば y ~ x1+x2+0 とすると，これは 𝑦 = 𝑥1 𝑏1 + 𝑥2 𝑏2 という回帰 式を表します。 • 変数の和や積などを説明変数として使いたい場合は，I() をかませることができま す。例えば y ~ I(x1+x2) とすると，これは 𝑦 = (𝑥1 + 𝑥2 )𝑏12 + 𝑏0 という回帰式を 表します。ただこれを使うくらいなら，事前にデータフレームに新しい変数を作成し ておいたほうがスマートな気がします。 それでは出力を見てみましょう。 1 2 3 4 5 6 7 Call: lm(formula = Q1_A ~ age + gender, data = dat) Coefficients: (Intercept) 18.05453 age 0.07083 gender 1.89117 下に出ている Coefficients が（回帰）係数を指しています。(Intercept) が切片を，それ以外 のものは各説明変数に対する回帰係数を表しています。 したがって， 重回帰モデルの式は 𝑦𝑝̂ = 18.05453 + 0.07083𝑥𝑝,age + 1.89117𝑥𝑝,gender となります。 変数 gender は二値変数のため， 男性 （gender==1） と女性 (gender==2) で分けると， 𝑦𝑝̂ = { 18.05453 + 0.07083𝑥𝑝,age + 1.89117 (男性) 18.05453 + 0.07083𝑥𝑝,age + 3.78234 (女性) ということですね。二値カテゴリ変数のコーディングに関して，通常は 0 と 1 にするのですが，こ れは 0 のカテゴリを基準カテゴリと見なすためです。ただ回帰係数という観点で言えば，0 と 1 にしないといけないということはなく，差が 1 ならば何でも良い*5 ことになります。もし dat で の gender の値が全員 1 ずつ小さい値で（つまり男性は 0，女性は 1 と）コーディングされていた としても，gender に対する回帰係数は変わりません。 少し切片が変わって 𝑦𝑝̂ = (18.05453 + 1.89117) + 0.07083𝑥𝑝,age + 1.89117𝑥𝑝,gender = 19.9457 + 0.07083𝑥𝑝,age + 1.89117𝑥𝑝,gender *5 回帰係数の検定という観点では差が 1 である必要すらも無いのですが， 回帰係数の解釈を考えるとやはり 0/1 のコ ーディングが一番しっくりきます。
5.2 重回帰分析 14 となります。 男女ごとに回帰式を見ると 𝑦𝑝̂ = { 19.9457 + 0.07083𝑥𝑝,age (男性) 19.9457 + 0.07083𝑥𝑝,age + 1.89117 (女性) となるわけです。 INFO 切片の役割 切片は「すべての説明変数の値が 0 のときの予測値」を表しています。したがって dat で言 えば「age が 0 ＝ 0 歳で，gender が 0 の人」の予測値になりますが当然そんな人はいませ ん。このように，多くの場合では切片の値それ自体は使えない（検定する意味もない）もので す。 もし切片に役割を持たせたい場合には，すべての量的変数を中心化すると良いです。中心 化，つまり平均値を引くことで各変数の値が 0 というのが「その変数が平均値の人」を表す ようになります。つまり全ての量的変数を中心化した状態で行った回帰分析の切片は「すべ ての変数が平均値の人」 の予測値になるわけです。 なお，gender のようなカテゴリカル変数は中心化してもしなくても良いです。例えば男女 比が 1:1 のときに中心化を行うと 「男性は-0.5， 女性は 0.5」 という感じになります。 こうなる とむしろ 0 が意味を持たなくなってしまうので，あえて中心化しないことで，切片項は「量 的変数はすべて平均値の男性での平均値」を意味するようになります。一方で，カテゴリカ ル変数も全て中心化した場合，切片項は純粋に全体平均値を表すようになります。というこ とでどちらにせよ解釈可能な意味を持つので，用途に応じて中心化するかを決めるのが良 いでしょう。 とりあえず回帰分析ができたので，次は回帰係数の検定をしてみます。といっても lm() には既 に回帰係数の検定が用意されているので難しいことはありません。 回帰係数の検定 1 2 3 # まずlm()の結果をオブジェクトに代入 result_lm <- lm(Q1_A ~ age + gender, data = dat) summary(result_lm) R では，関数の出力を含むあらゆるものをオブジェクト（変数）に代入することが出来ます。これ により，例えば「因子分析を行い，結果から因子負荷行列を取り出し，その行ごとの和をとる」みた いなことが全てコードとして書けるようになるわけです。同様に，特定の関数の出力を受け取って 別の処理を行うような関数もたくさん存在しています。ここでは関数の出力自体も変数として扱 えるという感覚を理解しましょう。 何かしらの分析を行う関数の出力オブジェクトは，たいてい list 型で与えられます。list 型 はほぼ何でも，いくつ入れても良い型です。例えば list 型である result_lm（lm() の出力）の中 身を見てみると， 以下のような状態になっています。
5.2 重回帰分析 15 list の中身 1 2 # 長いので省略しています str(result_lm) 1 2 3 4 5 6 List of 12 $ coefficients : Named num [1:3] 18.0545 0.0708 1.8912 ..- attr(*, "names")= chr [1:3] "(Intercept)" "age" "gender" $ residuals : Named num [1:2432] -1.079 -2.112 -4.041 -0.041 -1.15 ... ..- attr(*, "names")= chr [1:2432] "1" "2" "3" "4" ... $ effects : Named num [1:2432] -1.15e+03 4.00e+01 -4.38e+01 8.35e-02 -1.19 ... ..- attr(*, "names")= chr [1:2432] "(Intercept)" "age" "gender" "" ... $ rank : int 3 $ fitted.values: Named num [1:2432] 21.1 23.1 23 23 21.1 ... ..- attr(*, "names")= chr [1:2432] "1" "2" "3" "4" ... $ assign : int [1:3] 0 1 2 $ qr :List of 5 ..$ qr : num [1:2432, 1:3] -49.3153 0.0203 0.0203 0.0203 0.0203 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : chr [1:2432] "1" "2" "3" "4" ... .. .. ..$ : chr [1:3] "(Intercept)" "age" "gender" .. ..- attr(*, "assign")= int [1:3] 0 1 2 ..$ qraux: num [1:3] 1.02 1.02 1.02 ..$ pivot: int [1:3] 1 2 3 ..$ tol : num 1e-07 ..$ rank : int 3 ..- attr(*, "class")= chr "qr" [list output truncated] - attr(*, "class")= chr "lm" 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 一番上に List of 12 とあることから，このオブジェクトはいろいろなものが 12 個くっつい ている，ということがわかります。上の出力では，$記号で区切られた一つ一つが list 内の各要 素です。つまり一つ目が coefficients で，これは長さ 3 の名前付き実数型 (Named num) のベク トルとなっています。その下の residuals や effect は同じく名前付き実数型ベクトルですが長 さは 2432 のようです。また，少し下にある qr は List of 5 となっており，字下げされた 5 つの 要素 (qr, qraux, pivot, tol, rank) が存在しています。これは，list の中に list が入れ子 になっている，という状態を表しています。このように，list 型オブジェクトは型の違いも気に せずにいろいろなものをくっつけて一つのオブジェクトとして扱うことが出来ます。 list 型オブジェクトの中身は，それぞれの要素名の前に$がついていることから分かるように$ 記号を使って取り出すことが出来ます*6 。例えば以下のコードでは，result_lm の中から回帰係 *6 データフレームの列も$で取り出せるということは， データフレームは全ての要素が「同じ長さのベクトル」のリス トである，という見方もできるわけです。リストとして扱うことは多くは無いですが，知っておくとどこかで役立つ かもしれません。
5.2 重回帰分析
16
数の部分 (coefficients) だけをベクトルとして取り出すことができます。この行為に意味があ
るかはともかく，大抵のオブジェクトでは分析結果から中身を取り出せるという感覚だけでも理
解しておいてください。コレ，分析の手順の自動化率を高める際などにかなり重要になってきま
す。
結果の一部だけを取り出す
1
result_lm$coefficients
1
2
(Intercept)
18.05453074
age
0.07082883
gender
1.89116882
今回は lm() 関数の出力を result_lm というオブジェクトに代入しました。そしてこれに対し
て summary() 関数を実行しました。Chapter 2 では summary() 関数のことを「要約統計量を出
してくれる関数」と紹介しましたが，実際には summary() 関数はもっと広く「与えられたオブジェ
クトの型に合わせた何かしらのサマリーを出してくれる」
という関数なのです*7 。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
Call:
lm(formula = Q1_A ~ age + gender, data = dat)
Residuals:
Min
1Q
-19.1742 -2.6150
Median
0.6091
3Q
3.1820
Max
8.8502
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.054531
0.394846 45.726
<2e-16 ***
age
0.070829
0.008116
8.727
<2e-16 ***
gender
1.891169
0.189308
9.990
<2e-16 ***
--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.386 on 2429 degrees of freedom
Multiple R-squared: 0.07011,
Adjusted R-squared: 0.06934
F-statistic: 91.56 on 2 and 2429 DF, p-value: < 2.2e-16
出力は上から，
Residuals 𝑦 − 𝑦 ̂ の要約統計量です。回帰分析の前提条件として，誤差 𝜀 は正規分布 𝑁 (0, 𝜎𝜀2 ) に
*7 内部的には，summary() で呼び出せる複数の関数が存在しており，
与えるオブジェクトのクラスによって異なる
関数が実行されています。lm() の出力は lm クラスのオブジェクトであり，これに対しては summary.lm() が自
動的に呼び出されているのです。なので lm オブジェクトに対する summary() 関数のヘルプを見たい場合には
?summary.lm としてあげる必要があります。
5.2 重回帰分析 17 従う必要があるため， 特に中央値が 0 から大きく離れていたら注意が必要かもしれません。 Coefficients 回帰係数です。左から「係数の推定値」 「標準誤差」 「𝑡 統計量」 「𝑝 値」が表示されて います。 今回の場合， いずれも 𝑝 < .05 なので， 統計的に有意ということですね。 Multiple R-squared 重相関係数です。この後説明します。 Adjusted R-squared 自由度調整済み重相関係数です。この後説明します。 F-statistic 「切片以外の回帰係数が全て 0 である」という帰無仮説に対する 𝐹 統計量です。と いうことはこの検定が有意ではなかった場合， モデルとしての意味が無いことになります。 5.2.4 決定係数・重相関係数 summary(lm()) で出てくる Multiple R-squared は決定係数と呼ばれます。回帰分析モデル では， 結果変数 𝑦 は予測値 𝑦 ̂ と誤差 𝜀 によって (5.20) 𝑦 = 𝑦̂ + 𝜀 と表されます。 当然ながら， 予測値と誤差の間は無関係であるはずなので*8 ，𝑦 の分散は (5.21) 𝜎𝑦2 = 𝜎𝑦2̂ + 𝜎𝜀2 というように，予測値と測定誤差の分散の単純な和として表すことができます。信頼性係数のとき と同じ感じですね。 このとき， 決定係数 𝑅2 は 𝑅2 = 𝜎𝑦2̂ =1− 2 𝜎𝑦 𝜎𝜀2 𝜎𝑦2 (5.22) と定義されます。つまり，決定係数 𝑅2 は結果変数の分散のうち，説明変数で説明可能な分散の割 合ということです。 ここで， 𝑦 に影響（個人差）を与えるすべての説明変数が分かっている状況を考えてみましょ う。例えば Q1_A の値を決定づけるものは，幼いときの親の関わり方や小学校での友人関係，ある いは昨日の晩御飯や体温，回答の瞬間の気分なんかも影響するかもしれません。ほぼゼロだとして もゼロではない限りは重回帰モデルに含める意味があります。仮にちょうど 1 億個の説明変数を 使えば 𝑦 の個人差が完全に説明可能だったとすると，重回帰モデルは，以下のようになります（左 辺が 𝑦 ̂ ではなく 𝑦 である点に注意） 。 ∗ ∗ 𝑦𝑝 = 𝑏0 + 𝑥𝑝,1 𝑏1∗ + 𝑥𝑝,2 𝑏2∗ + 𝑥𝑝,3 𝑏3∗ + ⋯ + 𝑥𝑝,99999999 𝑏99999999 + 𝑥𝑝,100000000 𝑏100000000 (5.23) ただ， 実際問題として 1 億個の説明変数による回帰モデルが判明したとしても， このままでは変数 の数が多すぎて実社会には応用できません。1 億個あっても，大半はほぼ無視できるレベルの影 響しか与えていないでしょうし，できることならなるべく少ない変数の数で十分な予測モデルを 立てることができれば， みんな嬉しいはずです （オッカムの剃刀） 。 そこで，この 1 億個の変数の中 から，最も説明力が高い 2 つだけを選んできたとしましょう。今，1 億個の説明変数が「説明力の *8 回帰分析モデルでは， これも仮定のひとつとしておかれます。この仮定が守られていない場合，検定結果にバイアス がかかったりするわけです。
5.2 重回帰分析 18 高い順」 に並んでいたとすると， 上の重回帰モデルは次のように書き換えることができます。 𝑦 = 𝑦̂ + 𝜀 (5.24) 𝑦 ̂ = 𝑏 0 + 𝑥 1 𝑏1 + 𝑥 2 𝑏2 𝜀 = 𝑥3 𝑏3 + ⋯ + 𝑥99999998 𝑏99999998 + 𝑥99999999 𝑏99999999 + 𝑥100000000 𝑏100000000 もちろん，(5.24) 式の場合には，観測された 𝑦 を完全には説明できません。しかし，もしも説明 変数 2 つで 𝑦 の個人差のうちの結構な割合が説明できるとしたら，実社会への応用はかなり現実 味を帯びてきます。そもそも人間の行動を完全に予測するなんてのは無理な話なので， 「効率よく 大体予測できたら良い」くらいの気持ちで考えておけばよいわけです。このように，回帰分析では， どれだけ少ない変数で十分な説明ができるかが一つの重要な指標となります。 また，Multiple R-squared はその名の通り重相関係数の二乗という意味合いもあります。重 相関係数とは，𝐲 と 𝐲̂ の相関のことです。実際に相関係数を出してみましょう。なお predict() 関数は，lm() の出力を第一引数に， データを第二引数に与えると， 各行の説明変数の値と lm() で 得られた回帰係数を用いて 𝐲̂ を計算してくれる関数です。 重相関係数の二乗の計算 1 2 3 4 # 予測値を出す yhat <- predict(result_lm, dat) # yhatとyの相関 cor(yhat, dat$Q1_A)^2 1 [1] 0.07010747 確かに summary(lm()) で出てくる Multiple R-squared と一致する値が得られました。重相 関係数が大きいほど， 説明変数による予測 𝑦 ̂ は正確であるということです。 ちなみに重相関係数が 1 になるのは，𝑦 = 𝑦 ̂ のときのみです。今回の例では，𝑅2 = 0.0701 とかなり低い値を示しまし た。実際に散布図を見ても，図 5.7 のように予測はほとんど意味をなしていません。これは，Q1_A の得点の変動を説明するのに年齢 age と性別 gender は少なくとも線形ではほぼ意味をなさない （説明できてない）ということです。このように，重相関係数はモデルの適合度の指標としても用い ることができます。具体的にどの程度あれば十分という基準は特にないのですが，あまりに低い場 合には，回帰係数の検定が有意であったとしても「でも，結局誤差のほうが大きいんだよね？」とい う感じで評価されてしまうかもしれません。 回帰分析では，説明変数の数を増やすと必ず重相関係数は高くなります。仮に完全に無関係な変 数 （乱数など） を加えたとしても， 重相関係数が低下することはありません。 先程の 1 億個モデルで 言えば， 最も意味のない変数 𝑥100000000 を追加したとしても， 𝑦 = 𝑦̂ + 𝜀 𝑦 ̂ = 𝑥1 𝑏1 + 𝑥2 𝑏2 + 𝑥100000000 𝑏100000000 + 𝑏0 𝜀 = 𝑥3 𝑏3 + ⋯ + 𝑥99999998 𝑏99999998 + 𝑥99999999 𝑏99999999 となるわけで， わずかでも 𝜎𝜀2 は小さくなるでしょう。
5.2 重回帰分析 19 図 5.7: 予測値と真値の散布図 どれだけ少ない変数で十分な説明ができるかという視点では，重相関係数がほぼ向上しない無 駄な変数は取り除きたいところです。そこで使われることがあるのが自由度調整済み決定係数で す。これは，決定係数 𝑅2 に対して説明変数の数によるペナルティを与えたものです。つまり決定 係数が同じになるならば， 説明変数の数が少ないモデルのほうが良い値を示します。 2 𝑅𝑎𝑑𝑗 =1− 𝜎𝜀2 𝑛 − 1 𝜎𝑦2 𝑛 − 𝑘 − 1 2 一部の人はこの 𝑅𝑎𝑑𝑗 を，説明変数の数が異なるモデル同士の比較（ある変数を入れるか入れない 2 かの判断など）に用いているようです。ただ，式を見てわかるように 𝑅𝑎𝑑𝑗 にかかるペナルティは サンプルサイズ 𝑛 の影響も受けています。つまり，サンプルサイズが大きくなるほど説明変数の 追加に対するペナルティが弱くなるため，ほぼ無意味な変数を加えることによる 𝜎𝜀2 の減少のほ うがペナルティよりも大きくなる可能性が高まります。そして結果的に変数の数が多いモデルの ほうが選ばれやすくなる， という問題点があったりします。 …ということで長々と決定係数・重相関係数についてお話をしましたが，基本的には自由度調 2 整済み決定係数 𝑅𝑎𝑑𝑗 はあまり使えるものではない， という認識でも良いような気がします*9 。 5.2.5 標準偏回帰係数 偏回帰係数は，それぞれ「その説明変数が 1 大きくなったら 𝑦 ̂ はどの程度大きくなるか」を表し ているため， そのままでは比較することができません。age と gender の係数を比べて gender の ほうが影響が大きい，と判断はできないのです。その理由の一つは，変数のスケールが異なるため です。age が 1 増えると言うのは，年齢が一つ変わるだけでさほど大きな変化ではありません。対 して gender が 1 増えると言うのは，性別が男性から女性に変わることを意味するため，これはか なり大きな変化です。変数の変化の程度を調整するための方法として標準化というものがありま した。標準化したあとの変数は，すべて平均 0，分散 1 になるため，どんな変数でも「1 増える」とい うことの意味が 「1 標準偏差増える」 に揃えられます。 *9 良いモデルを選びたい場合には， 決定係数の比較よりも交差検証 (cross validation) などを用いるのが良いとされ ているようです。
5.2 重回帰分析 20 標準偏回帰係数とは，説明変数と結果変数のすべてを標準化した上で算出された偏回帰係数の ことです。 結果変数も標準化することによって， 係数は 「説明変数が 1 標準偏差大きくなったら，𝑦 が標準偏差いくつ分大きくなるか」を表すことになり，異なる変数の間でも大小の比較が可能そう な気がしてきます。 R には標準化を行う関数として scale() というものがあります。 [scale] 変数の平均・標準偏差を変える 1 2 x <- 1:10 scale(x) # 標準化 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 [,1] [1,] -1.4863011 [2,] -1.1560120 [3,] -0.8257228 [4,] -0.4954337 [5,] -0.1651446 [6,] 0.1651446 [7,] 0.4954337 [8,] 0.8257228 [9,] 1.1560120 [10,] 1.4863011 attr(,"scaled:center") [1] 5.5 attr(,"scaled:scale") [1] 3.02765 1 2 # 平均値を引いてから標準偏差で割るのと同じ結果 (x - mean(x)) / sd(x) 1 2 [1] -1.4863011 -1.1560120 -0.8257228 -0.4954337 -0.1651446 [7] 0.4954337 0.8257228 1.1560120 1.4863011 ということで， これを使って標準偏回帰係数を出してみましょう。 標準化してから回帰分析 1 2 3 4 5 6 7 # 使う変数だけ標準化 dat_std <- scale(dat[, c("Q1_A", "age", "gender")]) # matrix型になってしまうので明示的にデータフレームに変更 dat_std <- data.frame(dat_std) # 標準化したデータで回帰 result_lm <- lm(Q1_A ~ age + gender, data = dat_std) # あとはsummary() 0.1651446
5.3 一般化線形モデル
8
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
21
summary(result_lm)
Call:
lm(formula = Q1_A ~ age + gender, data = dat_std)
Residuals:
Min
1Q
-4.2178 -0.5752
Median
0.1340
3Q
0.7000
Max
1.9468
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.364e-16 1.956e-02
0.000
1
age
1.709e-01 1.958e-02
8.727
<2e-16 ***
gender
1.956e-01 1.958e-02
9.990
<2e-16 ***
--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9647 on 2429 degrees of freedom
Multiple R-squared: 0.07011,
Adjusted R-squared: 0.06934
F-statistic: 91.56 on 2 and 2429 DF, p-value: < 2.2e-16
age の係数が 1.709e-01，つまり 1.709 × 10−1 = 0.1709，そして gender の係数が 1.956e-01，
つまり 1.956 × 10−1 = 0.1956 となりました。この標準偏回帰係数を比較することで，やはり
gender のほうが Q1_A に与える影響は大きそうだ，ただあまり差はなさそうだ，ということがわ
かりました。
標準偏回帰係数による変数の比較は，常に良い方法ではないとされています。標準化という手続
きは，単純に平均値を引いてから標準偏差で割るだけです。したがって標準化前後では変数の分布
は変わりません。age と gender を見ても，gender は二値変数のため，分布の形は age とは完全
に別物でしょう。そのような状況下で，2 つの変数の「1 標準偏差」の意味は同じなのでしょうか。
King（1986）では，比較する変数の間に共通の単位があるような場合には標準化は有効だが，そう
でなければ比べようがないとしています。安易に標準偏回帰係数を使うのは危険かもしれない，
ということです。
5.3 一般化線形モデル
ここまで説明してきた回帰分析モデルは，
• 説明変数の単純な線形結合で結果変数を予測し
• 誤差は全データで共通の正規分布 𝑁 (0, 𝜎𝜀2 ) に従う
という仮定のもとでのモデルでした。ただ，世の中にはそんなに単純ではないケースが結構あり
ます。
例えば結果変数が
「病気になるかどうか」
の場合，
取りうる値は 0（病気にならない）
か 1（病
5.3 一般化線形モデル 22 気になる）の 2 つのみです。ですがこれまでの回帰モデルでは説明変数が増える（減る）ごとに 𝑦 ̂ は際限なく大きく （小さく） なってしまうため， 直線を当てはめるのはちょっと嫌な感じです。 この関係を表したものが 図 5.8 です。X 軸には説明変数として 「1 日あたりタバコの平均本数」 を，Y 軸には結果変数として 「肺がんになったかどうか」 をプロットしています。 普通に考えると， タバコの本数が多いほど肺がんになりやすいので，回帰直線を引くと右上がりになります。ですが 直線を引いてしまうと，20 本を超えてもなお値は大きくなってしまいます。実際に結果変数の取 りうる値が 0 か 1 であるため，ここでは Y 軸を「肺がんになる確率」として回帰直線を用いた予測 をしようとすると，例えば 1 日にタバコを 30 本吸う人は肺がんになる確率が 160% くらいにな ってしまいます。 これは困りました… 図 5.8: 二値変数に直線を当てはめると… また「誤差は正規分布 𝑁 (0, 𝜎𝜀2 ) に従う」という仮定に関しても問題がありそうです。図 5.8 で は，X 軸が 10 本くらいのところでは罹患している人としていない人がばらばらといるため， それ なりの分散がありそうな一方で，0 本や 20 本のあたりでは全員が同じ値になっています。つまり X 軸の値によって誤差分布の分散 𝜎𝜀2 が変動しているように思われるため，やはり直線を当ては めるのはおかしそうです。 ということで，回帰分析モデルを拡張する必要があります。予測値の変動を説明変数の線形結合 で表した (5.6) 式の右辺はわかりやすいのでこのまま使うとして，直線ではなく確率を表す，つま り取りうる値が 0 から 1 に収まるように変換することを考えます。いくつか選択肢はあるのです が， ここではロジスティック変換をしてみましょう。 ロジスティック関数とは， 𝑔(𝑥) = 1 1 + exp(−𝑥) (5.25) という関数です。図 5.9 は，𝑥 の値が −10 から 10 まで動いたときの 𝑔(𝑥) の値を表したもの です。このように，ロジスティック関数は 𝑥 の値が何であろうと確かに 𝑔(𝑥) は 0 から 1 までに収 まるようです。 ということで， ロジスティック関数の 𝑥 を説明変数の線形和 𝑏0 +𝑥1 𝑏1 +𝑥2 𝑏2 +𝑥3 𝑏3 +⋯+𝑥𝐾 𝑏𝐾 に置き換えてあげれば， 𝑦̂ = 1 1 + exp [−(𝑏0 + 𝑥1 𝑏1 + 𝑥2 𝑏2 + 𝑥3 𝑏3 + ⋯ + 𝑥𝐾 𝑏𝐾 )] (5.26)
5.3 一般化線形モデル 23 図 5.9: ロジスティック変換 となり， この 𝑦 ̂ は 0 から 1 までの値しか取らなくなります。 また， 式を変形させると log 𝑦̂ = 𝑏 0 + 𝑥 1 𝑏1 + 𝑥 2 𝑏2 + 𝑥 3 𝑏3 + ⋯ + 𝑥 𝐾 𝑏𝐾 1 − 𝑦̂ (5.27) という形になり， もとのきれいな右辺に戻すことも可能です。 ということで，このように回帰分析モデルを変形させることで正規分布以外の分布にも対応可 能にしたものが一般化線形モデル (generalized linear model; GLM) です*10 。一般化線形 モデルでは， モデル式が 𝑔(𝑦)̂ = 𝑏0 + 𝑥1 𝑏1 + 𝑥2 𝑏2 + 𝑥3 𝑏3 + ⋯ + 𝑥𝐾 𝑏𝐾 (5.28) という形になっています。つまり， 「予測値 𝑦 ̂ をなんか変換したもの」である 𝑔(𝑦)̂ を 𝑥 の線形和 𝐗𝐛 で表現しようという試みなわけです。この「𝑦 ̂ をなんか変換したもの」が，ロジスティック回帰 の場合には 𝑔(𝑦)̂ = log 𝑦̂ 1 − 𝑦̂ (5.29) というロジット変換*11 の式になっており，またふつうの重回帰モデルも 𝑔(𝑦)̂ = 𝑦 ̂ という恒等関数によって変換されたものと見れば，これが一般化線形モデルの特殊なケースだと いうことがわかるでしょう。変換の関数 𝑔(𝑦)̂ はリンク関数と呼びます。リンク関数は「二項分布の ときにはロジット変換」 というように基本的には 𝑦 ̂ の分布に対してお決まりのものがあります。 一 般化線形モデルでは変形の形ごとに異なる名称がついており，ロジスティック変換を用いた場合 をロジスティック回帰と呼びます。もう少ししっかりとロジスティック回帰モデルの全体を書く *10 これに対して正規分布の普通の重回帰モデルを一般線形モデル (general linear model) と呼んだりします。 *11 ロジスティック関数とロジット関数は互いに逆関数の関係になっています。 そのため，説明変数側の変換として見 た場合にはロジスティック変換と呼ばれる一方で，結果変数側の変換（リンク関数）として見る場合にはロジット変 換と呼ばれるわけです（たぶん）。実際に回帰分析としての手法名では「ロジスティック回帰」とも「ロジット回帰」と も呼ばれることがあるはずです。
5.3 一般化線形モデル と， 𝑦𝑝 = 24 1 + 𝜀𝑝 1 + exp [−(𝑏0 + 𝑥𝑝1 𝑏1 + 𝑥𝑝2 𝑏2 + 𝑥𝑝3 𝑏3 + ⋯ + 𝑥𝑝𝐾 𝑏𝐾 )] (5.30) = 𝑦𝑝̂ + 𝜀𝑝 となります。 なお 𝑦𝑝 は 0 か 1 の値しか取り得ないため， 誤差 𝜀𝑝 は 1 ⎧ {− 1+exp[−(𝑏0 +𝑥𝑝1 𝑏1 +𝑥𝑝2 𝑏2 +𝑥𝑝3 𝑏3 +⋯+𝑥𝑝𝐾 𝑏𝐾 )] 𝜀𝑝 = ⎨ 1 − 1+exp[−(𝑏 +𝑥 𝑏 +𝑥 1𝑏 +𝑥 𝑏 +⋯+𝑥 𝑏 )] { ⎩ 0 𝑝1 1 𝑝2 2 𝑝3 3 𝑝𝐾 𝐾 (𝑦𝑝 = 0の場合) (𝑦𝑝 = 1の場合) (5.31) という値をとる 2 値変数と見ることができます。この 𝜀𝑝 に入る具体的な値を考えてみましょう。 まず 𝑦 ̂ = 0.99 の場合には， 𝜀𝑝 = { −0.99 0.01 (𝑦𝑝 = 0の場合) (𝑦𝑝 = 1の場合) (5.32) となります。 そして， 「𝑦𝑝 = 1の場合」 に該当する確率は 𝑦 ̂ = 0.99 のため， この場合 • |𝜀𝑝 | が 0.99 になる確率は 0.01（𝑦𝑝 = 1の場合） • |𝜀𝑝 | が 0.01 になる確率は 0.99（𝑦𝑝 = 0の場合） となります。𝑦𝑝 = 1 の場合には誤差が 0.99 と相当大きくなるのですが， そもそもこのような誤 差が観測される確率が 0.01 しかないので，基本的には誤差は 0.01 となります。結果として，|𝜀𝑝 | の期待値 （または 𝜀𝑝 の分散） は小さめの値となります。 一方で 𝑦 ̂ = 0.5 の場合には， 𝜀𝑝 = { −0.5 0.5 (𝑦𝑝 = 0の場合) (𝑦𝑝 = 1の場合) (5.33) となります。 そして， 「𝑦𝑝 = 1の場合」 に該当する確率は 𝑦 ̂ = 0.5 のため， この場合 • |𝜀𝑝 | が 0.5 になる確率は 0.5（𝑦𝑝 = 1の場合） • |𝜀𝑝 | が 0.5 になる確率は 0.5（𝑦𝑝 = 0の場合） となり，|𝜀𝑝 | の期待値 （または 𝜀𝑝 の分散） は先程よりも大きな値となるわけです。 このように 𝜀𝑝 の分散は，予測確率 𝑦 ̂ が 0.5 付近では最も大きく，0 や 1 付近では小さくなるように変動します。 このような特性があるため GLM では，予測確率に応じて誤差分散が変動することを織り込んだ 上でモデルのパラメタ （回帰係数） を推定できるのです。 INFO プロビット回帰 リンク関数を標準正規分布累積確率関数の逆関数 𝑔(𝑦) = Φ−1 (𝑦) と置いた GLM は，プロ ビット回帰と呼びます。ロジスティック回帰とプロビット回帰はどちらも「結果変数が 0 か 1 の二値」というケースに対して適用されるモデルです。そして 2 つの回帰で用いられるリ ンク関数の逆関数はよく似た形（左右対称の釣鐘型）をしています。そんな 2 つの回帰モデ ルについて，特に経済学など？では理論的な背景なども考慮して厳密に使い分けられる節
5.3 一般化線形モデル 25 があるようなのですが，この講義で紹介する手法（特に項目反応理論）では，基本的にその区 別には関心がありません。言ってしまえば「どっちでも良い」という感じがあり，基本的には 計算しやすいという理由からロジスティック回帰のほうを使用することが多いです。 ということで，先程のタバコの例でもロジスティック回帰による回帰直線を当てはまると， 図 5.10 のようになります。 図 5.10: ロジスティック回帰直線をあてはめる 一般化線形モデルでは他にもポアソン回帰やガンマ回帰を始めとしていくつかありますが，こ こでは全ては紹介しません。以降の回では，ロジスティック回帰だけ理解しておけばついていける はずです。 5.3.1 回帰係数の解釈 回帰係数の解釈は，基本的には普通の重回帰分析と同じです。すなわち，ある説明変数の回帰係 数が正の値であれば， その説明変数が大きくなるほど予測値 𝑦 ̂ も大きくなることを意味します。 逆 に回帰係数が負の値であれば，その説明変数が大きくなるほど予測値 𝑦 ̂ は小さくなることを意味 します。ただし GLM の場合，予測値 𝑦 ̂ は説明変数の線形結合をリンク関数で変換したものであ るため，回帰係数の解釈もリンク関数を考慮する必要があります。ロジスティック回帰の場合，リ 𝑦̂ ンク関数は (5.29) 式に示されていたロジット変換でした。この式の右辺の 1− オッズ 𝑦̂ の部分は， （odds）と呼ばれるものです。オッズは「ある事象が起こる確率 𝑦 」 ̂ を「その事象が起こらない確率 1 − 𝑦」 ̂ で割ったもので，例えば 𝑦 ̂ = 0.8 のときにはオッズは 4 になります。この値は，例えばある 事象について成功／失敗のような二値の結果のいずれかが実現すると考えると，成功確率が 80% のときには， 「成功」 という事象が 「失敗」 という事象の 4 倍の確率で起こる， ということを意味して います。 (5.29) 式の両辺を exp で変換すると，ロジスティック回帰の回帰係数は，説明変数が 1 単位増 えると予測確率 𝑦 ̂ のオッズが exp(𝑏) 倍になるという解釈になると分かります。 ある説明変数 𝑥 の 回帰係数 𝑏1 が 0.3 だった場合，説明変数が 1 単位増えると予測確率のオッズが exp(0.3) ≈ 1.35 倍になる， ということです。 このように， ロジスティック回帰の回帰係数は， 説明変数が 1 単位増え ると予測確率のオッズがどれくらい変化するかを表すものとして解釈されます。このような解釈
5.3 一般化線形モデル 26 は， ロジスティック回帰の回帰係数をオッズ比と呼ぶこともある理由の一つです。 INFO そもそもなぜオッズで解釈するのか 直感的には，回帰係数を「説明変数が 1 単位増えると予測確率 𝑦 ̂ がどれくらい変化するか」 と解釈したいところですが，なぜオッズの変動の大きさ（オッズ比）で解釈する必要がある のでしょうか。それは，確率の変化のインパクトがベースラインによって異なるためです。 例えば，被説明変数を 「病気の罹患の有無」だとした場合に， 「𝑦 ̂ が 0.01 から 0.06 に変化した 場合」と「𝑦 ̂ が 0.5 から 0.55 に変化した場合」を比べてみましょう。予測確率の変動幅はどち らも 0.05 ですが，明らかに前者のほうが変化のインパクトが大きくなっているような気が します。なにしろ確率が 6 倍になっているのですから。そしてオッズ比は，この「変化のイン パクト」のようなものを，確率そのものよりは良く反映できているのです。実際に上の 2 つ のケースでオッズの変化を見てみると， 0.06 1−0.06 • 𝑦 ̂ が 0.01 から 0.06 に変化した場合のオッズの変化（オッズ比）は 0.01 ≈ 6.32 倍 1−0.01 0.55 1−0.55 • 𝑦 ̂ が 0.5 から 0.55 に変化した場合のオッズの変化（オッズ比）は 0.5 ≈ 1.22 倍 1−0.5 となります。 図 5.11 は， 予測確率 𝑦 ̂ と対数オッズ （5.28 式の左辺） の関係を表したものです。 この図から， オッズの変化のインパクトは 𝑦 ̂ が 0.5 を中心として対称的であることや，𝑦 ̂ が 0 や 1 に近 づくほどオッズの変化のインパクトが大きくなっていることがわかります。 図 5.11: 𝑦 ̂ と対数オッズの関係 改めて，ロジスティック回帰における回帰係数を「オッズの変動幅（オッズ比）」として解釈 することを踏まえて「病気の罹患の有無」の予測を例に考えてみましょう。例えば「ある薬の 投薬量」という説明変数の回帰係数が −0.3 だった場合，この薬の投薬量が 1 単位増えると 予測確率のオッズが exp(−0.3) ≈ 0.741 倍になると解釈されます。ただしこれが実際の確 率をどれくらい変化させるかは， ベースラインの予測確率 𝑦 ̂ によって異なります。 •（他の説明変数をもとに予測された）ベースラインの予測確率 𝑦 ̂ が 0.01 だった場合
5.3 一般化線形モデル 27 には， オッズが 0.741 倍になると予測確率は 0.0074 くらいになります。 • ベースラインの予測確率 𝑦 ̂ が 0.5 だった場合には，オッズが 0.741 倍になると予測 確率は 0.426 くらいになります。 • またベースラインの予測確率 𝑦 ̂ が 0.9 だった場合には，オッズが 0.741 倍になると 予測確率は 0.870 くらいになります。 このように，オッズ比は同じであっても，実際の確率の変化の大きさはベースラインの予測 確率 𝑦 ̂ によって大きく異なる点に注意して解釈してください。 また 図 5.11 からわかる重要な点として， ロジスティック回帰は暗に • 確率の変化のインパクトがベースラインによって異なり • そのインパクトは 𝑦 ̂ = 0.5（オッズ 1）を中心として対称（0.5 → 0.3 と 0.5 → 0.7 は 同じ） であり • 確率 0 や 1 に近づくほどインパクトが大きくなる という前提を置いていることになります。これは，通常の回帰分析で「説明変数の影響が，𝑦 ̂ の値によらず常に一定」という前提を置いているのと同じようなものだと言えます。したが って，例えば「S 字カーブではあるが左右対称ではない（0 → 0.5 の変化の仕方と 0.5 → 1 の 変化の仕方が異なる）」ような場合には，もしかしたらロジスティック回帰はあまり適切な モデルではないのかもしれません。具体例としては，上の「投薬量」の例でベースラインの予 測確率が非常に高いケースでは，本来は薬がきちんと効くことで罹患確率が顕著に低下す ることが予想されるとしても， ロジスティック回帰の前提のもとでは 0.9 → 0.87 と， あまり 変わらないという予測になってしまいます。もちろんこれは「もともとの罹患確率が 0.9 と ほぼ確実に罹患すると思われる人なので，投薬してももう手遅れ」ということを表している 可能性もありケースバイケースではあるのですが，もしかしたらロジスティック回帰の前 提があまり適切ではないケースもあるかもしれない，ということは頭の片隅に置いておく と良いかもしれません。 5.3.2 R で一般化線形モデル R で一般化線形モデルを行う場合は，lm() の代わりに glm() という関数を使います。基本的 には使い方は同じですが，GLM なので「誤差分散の分布の形」および「リンク関数」を指定してあ げる必要があります。指定の仕方にちょっとクセがありますが，そういうものとして理解してくだ さい。 GLM（ロジスティック回帰） 1 2 3 4 # 結果変数は二値変数でないといけないので，「Q1_Aが平均以上か」を表す変数を作成 # datには追加していない Q1_A_binom <- dat$Q1_A > mean(dat$Q1_A) # 実は引数dataの中に無い変数でも分析に利用できる
5.3 一般化線形モデル
5
6
7
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
28
# ↑lm()でも可能です
result_glm <- glm(Q1_A_binom ~ age + gender, data = dat, family =
binomial(link = "logit"))
summary(result_glm)
Call:
glm(formula = Q1_A_binom ~ age + gender, family = binomial(link = "logit"),
data = dat)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.759336
0.189153 -9.301 < 2e-16 ***
age
0.025611
0.003922
6.530 6.60e-11 ***
gender
0.706794
0.088404
7.995 1.29e-15 ***
--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 3358.4
Residual deviance: 3246.0
AIC: 3252
on 2431
on 2429
degrees of freedom
degrees of freedom
Number of Fisher Scoring iterations: 4
まず引数 family には誤差分散の分布の形を指定します。ロジスティック回帰の場合は誤差分
散は 𝜀 ∼ 𝑁 (0, 𝑦(1 − 𝑦)) ということで二項分布
（ベルヌーイ分布）
が指定されます。
また，
リンク関
数は誤差分散の分布の関数の引数として指定する必要があります。ロジスティック回帰の場合は
ロジット変換をしているため，結果的に引数 family は binomial(link="logit") となるわけで
す。もしここでプロビット回帰をしたい場合には，binomial(link="probit") としたら良いわ
けですね。
また，当然 glm() でも正規分布の普通の重回帰分析を行うことも出来ます。この場合，family
には正規分布を示す gaussian を，リンク関数には恒等関数を表す identity を指定してあげて
ください。指定しなくてもデフォルトでそうなっていますが，ここでは明示的に指定しておきま
す。
GLM（正規分布もできるよ）
1
2
result_glm <- glm(Q1_A ~ age + gender, data = dat, family = gaussian(link =
"identity"))
summary(result_glm)
5.4 変数間の関係のグラフィカルモデル
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
29
Call:
glm(formula = Q1_A ~ age + gender, family = gaussian(link = "identity"),
data = dat)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 18.054531
0.394846 45.726
<2e-16 ***
age
0.070829
0.008116
8.727
<2e-16 ***
gender
1.891169
0.189308
9.990
<2e-16 ***
--Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 19.23359)
Null deviance: 50241
Residual deviance: 46718
AIC: 14097
on 2431
on 2429
degrees of freedom
degrees of freedom
Number of Fisher Scoring iterations: 2
当然ですが，lm() のときと同じ結果が得られます。
5.4 変数間の関係のグラフィカルモデル
最後に，今後の授業で登場するグラフィカルモデル（あるいは有向非巡回グラフ; Directed
acyclic graph [DAG]）の導入をしておきます。といっても難しい事はありません。単に説明変数と
被説明変数の関係性を矢印で表すことで，
モデルの見通しを良くするために使うツールです。
回帰分析の例では，gender および age によって Q1_A の値を予測してみました。これをグラフ
ィカルモデルで表すと図 5.12 のようになります。
図 5.12: 重回帰モデル
このように，グラフィカルモデルでは観測可能な変数を長方形で表し，矢印によって回帰の向
きを表します。また，両方向の矢印は変数間の相関関係を表しています。重回帰分析では，複数の
説明変数が一つの結果変数を説明しようとするため，グラフィカルモデルはこのように複数の変
数から一つの変数に向かって矢印がひかれる形になります。
なお，本資料の範囲においては，グラフィカルモデルの矢印はあくまでも回帰の関係を表すもの
5.4 変数間の関係のグラフィカルモデル であって，因果関係を表すものではないことに注意してください。例えば図 5.12 では，age から Q1_A（協調性）に向かって矢印が引かれていますが，これは「年齢が高いほど Q1_A の値が大きく なる傾向がある」ということを表しているだけであって， 「年を取ることで Q1_A の値が大きくな る」 という因果関係を表しているわけではない， ということです。 そして，基本的には回帰分析の結果から因果関係を推測することはできないということも覚え ておいてください。もしも age から Q1_A の回帰係数が有意であったとしても，それは「年齢が高 いほど Q1_A の値が大きくなる傾向がある」ということを表しているだけであって， 「 （一人の人間 が）年齢を重ねることで Q1_A が高まっていく」という因果関係を表しているわけではないので す。ですが，実際に世の中では回帰分析の結果をあたかも因果関係のように解釈してしまっている ケースは結構あるように思います。直感的に考えると， 「age から Q1_A の回帰係数が有意になる」 ようなデータにおいては，たぶん「Q1_A から age に回帰した場合も同じように有意になる」こと が多いでしょう。age と Q1_A の間に相関関係があるならば， どちらを説明変数にしても有意な回 帰係数が出る可能性が高いのです。 …と考えると，回帰分析の結果を因果関係として解釈することがいかに間違っているかがわか るのではないでしょうか。回帰分析はあくまでも「変数間の関係性を表すモデル」であって， 「変数 間の因果関係を表すモデル」 ではないのです。 図 5.13: パス解析 グラフィカルモデルを自由に拡張して，例えば 図 5.13 のように変数間の関係を広げていった ものが，パス解析と呼ばれるやつです。このモデルを分析する場合，矢印一個一個を回帰分析して も良いかもしれませんが，全体として変数間のモデルができているのですから，全部まとめてデー タとの適合度などを評価してあげたいところです。 ということでパス解析も SEM の一種です。 SEM（第 7 章）では，このような変数間のグラフィカルモデルをもとに分析用のコードを作成 していくことになります。といっても難しいことはないはずです。データを取る時点で変数間の理 論的な関係性の仮説は考えているはずですから… 30
参考文献 参考文献 King, G. (1986). How not to lie with statistics: avoiding common mistakes in quantitative political science. American Journal of Political Science, 30(3), 666. https://doi.org/10.2 307/2111095 31