多変量解析_05_回帰分析のおさらい
2.9K Views
July 25, 23
スライド概要
神戸大学経営学研究科で2022年度より担当している「統計的方法論特殊研究（多変量解析）」の講義資料「05_回帰分析のおさらい」です。
2025/2/5: 全体的に内容に手を入れました。RmarkdownからQuartoに変更しました。
神戸大学経営学研究科准教授 分寺杏介（ぶんじ・きょうすけ）です。 主に心理学的な測定・教育測定に関する研究を行っています。 講義資料や学会発表のスライドを公開していきます。 ※スライドに誤りを見つけた方は，炎上させずにこっそりお伝えいただけると幸いです。
関連スライド
各ページのテキスト
5.1 回帰分析とは 1 Chapter 5 回帰分析のおさらい 各種分析手法（因子分析・構造方程式モデリング・項目反応理論）の紹介に向けた準備として， 各手法のベースにある回帰分析の理論的な部分をおさらいしつつ，いくつかの重要な概念を説明 しています。 本資料は，神戸大学経営学研究科で 2022 年度より担当している「統計的方法論特殊研究（多 変量解析） 」の講義資料です。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 4 3 2 y 1 0 -1 -2 -3 -4 -3 -2 -1 0 x 1 2 3 図 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 つの点に関してはズ
5.1 回帰分析とは 3 4 3 2 y 1 0 -1 -2 -3 -4 -3 -2 -1 0 x 1 2 3 図 5.2: 3 つの点を通る直線を求める 4 3 2 y 1 0 -1 -2 -3 -4 -3 -2 -1 0 x 1 2 3 図 5.3: 3 つの点を通る直線なんてない レがゼロですが，残り 1 個の点に対しては大きくズレています。この赤い線のように，平均的な ズレを最小化する直線を求める分析法が回帰分析なのです。 点の数が何個であろうと最もベーシックな回帰分析では，直線（一次関数）を求めます。傾き を 𝑏1 ，切片を 𝑏0 とおけば 𝑦 = 𝑏 0 + 𝑏1 𝑥 (5.1) という直線について， 𝑏0 , 𝑏1 の値がいくつのときに最もそれっぽい直線になるかを求めるわけ です。 5.1.1 最小二乗法 実際の分析においては，最小二乗法によって回帰直線を求めることが多いです。これは図 5.4 の青い矢印のように，縦方向でのズレを求めて，これの二乗和が最小になる (𝑏1 , 𝑏0 ) の組み合わ せを求める方法です。 「それっぽい直線」という観点で言えば，図 5.5 の青い矢印のように，赤い直線との直線距離 に関して最小化しても良さそうな気がしますが，最小二乗法では縦方向での距離を最小化してい ます。なぜなのでしょうか。
5.1 回帰分析とは 4 1 y 0 -1 -2 -3 -1 0 x 1 図 5.4: 最小二乗法 1 y 0 -1 -2 -3 -1 0 x 1 図 5.5: 直線距離で考えたらダメなのかい それは，回帰分析が 𝑥 で 𝑦 を説明する分析手法だからです。𝑦 そのものは無いけれど他の変数 がある時に 𝑦 を近似するのが回帰直線なのです。例えば，𝑥 を「体重」，𝑦 を「身長」だとしま しょう。あなたのもとに，あるゲストがやってきます。どういうわけかその人に関して「体重」 の情報だけが入っているとします。あなたはゲストのためにちょうどよい高さの椅子を用意する 必要があります。しかし，ちょうどよい高さの椅子をあてがうためには，ゲストの身長が分かっ ていないといけません*1 。さて，あなたはゲストの身長をどのように見積もるでしょうか。 普通ならば，図 5.6 のように，周りの人の体型などから「 （標準体型で）体重が 70kg くらいの 人ならば，身長はだいたい 170cm くらい」 「 （標準体型で）体重が 60kg くらいの人ならば，身長 はだいたい 162cm くらい」という感じで予測を立てたうえで，ゲストの体重の情報を当てはめ て「ゲストの体重が 65kg ということは，標準体型ならば身長はだいたい 165.5m くらいかなぁ」 といった予測を行うでしょう。回帰分析で行っていることは概ねこんな感じのことですが，ここ で問題になる「予測のズレ」は，身長（つまり 𝑦）に関してのみです。つまり，𝑥 による 𝑦 の近似 （予測）のズレが最小になるように回帰直線を引きたいがために，最小二乗法では縦方向でのズ レのみを考えているわけですね。 *1 もちろん実際には股下やら体型やらによって最適な椅子の高さは決まるわけですが，手元には「体重」の情報し か無いので，できる範囲で予測しましょう，ということです。
5.1 回帰分析とは 5 身長 180 160 140 40 50 60 体重 70 80 図 5.6: 回帰直線による予測 求めたい回帰直線は 𝑦 = 𝑏0 + 𝑏1 𝑥 ですが，一つ一つの点について考えると，ほとんどの場合 等号は成立しません。というのも回帰直線が示すのはあくまでも予測値であって，実際のデータ では予測とのズレがあるためです。そこで回帰分析モデルでは，実際の一つ一つのデータの値 𝑦𝑝 は，回帰直線による予測値 𝑦 ̂𝑝 と誤差 𝜀𝑝 の和であると考えます。ということで，回帰分析モデル の全体は 𝑦𝑝 = 𝑦̂𝑝 + 𝜀𝑝 (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.2 重回帰分析 6 5.1.2 （単）回帰分析の線形代数的表現 因子分析や SEM のところでは，それぞれの分析手法の基本的な考え方を説明するために線形 代数的表現を多用しています。ということで，ここで回帰分析の線形代数による表現を説明して おきます。書き方とその意味に慣れておいてください。 まず，回帰分析モデルを，図 5.1 のところで示した連立方程式的な書き方で書き直すと，次の ようになります。 (𝑥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 {⎡ 𝑦 ⎤ ⎡ 𝑥 ⎤ ⎡ 𝜀 ⎤ 2 {⎢ 2 ⎥ ⎢ ⎥ ⎢ 2 ⎥ ⋮ ⎥ = 𝑏 0 + 𝑏1 ⎢ ⋮ ⎥ + ⎢ ⋮ ⎥ ⎢ ⎨ ⎢𝑥𝑃 −1 ⎥ ⎢𝜀𝑃 −1 ⎥ {⎢𝑦𝑃 −1 ⎥ { 𝑦 ⎩⎣ 𝑃 ⎦ ⎣ 𝑥𝑃 ⎦ ⎣ 𝜀𝑃 ⎦ (5.5) だいぶ見通しが良くなってきました。あとはベクトルを太字で表すことで，以下のように割と 綺麗に書き直すことが出来ます。 2 つの変数に関する長さ 𝑃 のベクトル 𝐱 = [𝑥1 𝑥2 ⋯ ⊤ 𝑥𝑃 ] , 𝐲 = [𝑦1 が与えられたとき，以下の式における誤差ベクトル 𝜺 = [𝜀1 𝜀2 ⋯ 𝑦2 ⋯ ⊤ 𝑦𝑃 ] ⊤ 𝜀𝑃 ] の二乗和 𝜺⊤ 𝜺 が最小になる係数 𝑏0 , 𝑏1 の値はそれぞれいくつか。 𝐲 = 𝑏 0 + 𝑏1 𝐱 + 𝜺 ということで，これでデータが何個あっても一つの式で表すことができるようになりました。 便利ですね。 5.2 重回帰分析 𝑦 を予測するという意味では，もちろん変数は多い方が良いでしょう。 「身長」を予測するうえ で，体重だけでなく「性別」や「体脂肪率」といったことがわかれば，予測はより正確になるは
5.2 重回帰分析 7 ずです。回帰分析でも，説明変数を複数個用いることが出来ます。これを一般的には重回帰分析 と呼びます。説明変数が 𝐾 個あるならば，対応する回帰係数を 𝐾 + 1（切片の分）個用意して， 以下のような直線を考えるわけです*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 直線がイメージしづらいかもしれませんが，説明変数が 1 個の場合の回帰直線が 2 次元平面に引かれていたよ うに，説明変数が 2 個の場合は 3 次元空間の中に一本の直線が引かれます。同様に，説明変数が 𝐾 個ならば 𝐾 + 1 次元空間に一本の直線を引くわけです。
5.2 重回帰分析 8 これ以降，説明変数には 2 つの添字が付きます。2 つ目の添字 𝑘(𝑘 = 1, 2, ⋯ , 𝐾) は，何番目の 説明変数であるかを表すものです。例えば 𝑥𝑝3 （または 𝑥𝑝,3 ）は「𝑝 番目の人の 3 番目の説明変 数の値」（ 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 で重回帰分析 読み込み 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 という回 帰式を表します。 • 2 つの説明変数を *でつなぐと，その変数の主効果と交互作用項を同時に指定でき *4 ちなみに機械学習では one-hot encoding という名で使われるケースもあります。推定できなくなってしまうの は，説明変数の効果が線形和で表される回帰モデルならではの問題なのです。 12
5.2 重回帰分析 13 ます。例えば 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
20
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
出力は上から，
*7 内部的には，summary() で呼び出せる複数の関数が存在しており，与えるオブジェクトのクラスによって異な
る関数が実行されています。lm() の出力は lm クラスのオブジェクトであり，これに対しては summary.lm()
が自動的に呼び出されているのです。なので lm オブジェクトに対する summary() 関数のヘルプを見たい場合
には ?summary.lm としてあげる必要があります。
5.2 重回帰分析 17 Residuals 𝑦 − 𝑦 ̂ の要約統計量です。回帰分析の前提条件として，誤差 𝜀 は正規分布 𝑁 (0, 𝜎𝜀2 ) に従う必要があるため，特に中央値が 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 億個あっても，大半はほぼ無視できるレベル *8 回帰分析モデルでは，これも仮定のひとつとしておかれます。この仮定が守られていない場合，検定結果にバイ アスがかかったりするわけです。
5.2 重回帰分析 18 の影響しか与えていないでしょうし，できることならなるべく少ない変数の数で十分な予測モデ ルを立てることができれば，みんな嬉しいはずです（オッカムの剃刀）。そこで，この 1 億個の 変数の中から，最も説明力が高い 2 つだけを選んできたとしましょう。今，1 億個の説明変数が 「説明力の高い順」に並んでいたとすると，上の重回帰モデルは次のように書き換えることがで きます。 𝑦 = 𝑦̂ + 𝜀 (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 億個モ
5.2 重回帰分析 19 30 y 20 10 20 22 24 yhat 26 28 図 5.7: 予測値と真値の散布図 デルで言えば，最も意味のない変数 𝑥100000000 を追加したとしても， 𝑦 = 𝑦̂ + 𝜀 𝑦 ̂ = 𝑥1 𝑏1 + 𝑥2 𝑏2 + 𝑥100000000 𝑏100000000 + 𝑏0 𝜀 = 𝑥3 𝑏3 + ⋯ + 𝑥99999998 𝑏99999998 + 𝑥99999999 𝑏99999999 となるわけで，わずかでも 𝜎𝜀2 は小さくなるでしょう。 どれだけ少ない変数で十分な説明ができるかという視点では，重相関係数がほぼ向上しない無 駄な変数は取り除きたいところです。そこで使われることがあるのが自由度調整済み決定係数で す。これは，決定係数 𝑅2 に対して説明変数の数によるペナルティを与えたものです。つまり決 定係数が同じになるならば，説明変数の数が少ないモデルのほうが良い値を示します。 2 𝑅𝑎𝑑𝑗 =1− 𝜎𝜀2 𝑛 − 1 𝜎𝑦2 𝑛 − 𝑘 − 1 2 一部の人はこの 𝑅𝑎𝑑𝑗 を，説明変数の数が異なるモデル同士の比較（ある変数を入れるか入れな 2 いかの判断など）に用いているようです。ただ，式を見てわかるように 𝑅𝑎𝑑𝑗 にかかるペナルティ はサンプルサイズ 𝑛 の影響も受けています。つまり，サンプルサイズが大きくなるほど説明変数 の追加に対するペナルティが弱くなるため，ほぼ無意味な変数を加えることによる 𝜎𝜀2 の減少の ほうがペナルティよりも大きくなる可能性が高まります。そして結果的に変数の数が多いモデル のほうが選ばれやすくなる，という問題点があったりします。 …ということで長々と決定係数・重相関係数についてお話をしましたが，基本的には自由度調 2 整済み決定係数 𝑅𝑎𝑑𝑗 はあまり使えるものではない，という認識でも良いような気がします*9 。 5.2.5 標準偏回帰係数 偏回帰係数は，それぞれ「その説明変数が 1 大きくなったら 𝑦 ̂ はどの程度大きくなるか」を表 しているため，そのままでは比較することができません。age と gender の係数を比べて gender *9 良いモデルを選びたい場合には，決定係数の比較よりも交差検証 (cross validation) などを用いるのが良いとさ れているようです。
5.2 重回帰分析 20 のほうが影響が大きい，と判断はできないのです。その理由の一つは，変数のスケールが異なる ためです。age が 1 増えると言うのは，年齢が一つ変わるだけでさほど大きな変化ではありませ ん。対して gender が 1 増えると言うのは，性別が男性から女性に変わることを意味するため， これはかなり大きな変化です。変数の変化の程度を調整するための方法として標準化というもの がありました。標準化したあとの変数は，すべて平均 0，分散 1 になるため，どんな変数でも「1 増える」ということの意味が「1 標準偏差増える」に揃えられます。 標準偏回帰係数とは，説明変数と結果変数のすべてを標準化した上で算出された偏回帰係数の ことです。結果変数も標準化することによって，係数は「説明変数が 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 ということで，これを使って標準偏回帰係数を出してみましょう。 0.1651446
5.2 重回帰分析
21
標準化してから回帰分析
1
2
3
4
5
6
7
8
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# 使う変数だけ標準化
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()
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 一般化線形モデル 22 5.3 一般化線形モデル ここまで説明してきた回帰分析モデルは， • 説明変数の単純な線形結合で結果変数を予測し • 誤差は全データで共通の正規分布 𝑁 (0, 𝜎𝜀2 ) に従う という仮定のもとでのモデルでした。ただ，世の中にはそんなに単純ではないケースが結構あ ります。例えば結果変数が「病気になるかどうか」の場合，取りうる値は 0（病気にならない）か 1（病気になる）の 2 つのみです。ですがこれまでの回帰モデルでは説明変数が増える（減る）ご とに 𝑦 ̂ は際限なく大きく（小さく）なってしまうため，直線を当てはめるのはちょっと嫌な感じ です。 この関係を表したものが 図 5.8 です。X 軸には説明変数として「1 日あたりタバコの平均本 数」を，Y 軸には結果変数として「肺がんになったかどうか」をプロットしています。普通に考 えると，タバコの本数が多いほど肺がんになりやすいので，回帰直線を引くと右上がりになりま す。ですが直線を引いてしまうと，20 本を超えてもなお値は大きくなってしまいます。実際に 結果変数の取りうる値が 0 か 1 であるため，ここでは Y 軸を「肺がんになる確率」として回帰 直線を用いた予測をしようとすると，例えば 1 日にタバコを 30 本吸う人は肺がんになる確率が 160% くらいになってしまいます。これは困りました… 肺がん罹患有無 1.5 1.0 0.5 0.0 0 10 20 1日あたりタバコの平均本数 30 図 5.8: 二値変数に直線を当てはめると… また「誤差は正規分布 𝑁 (0, 𝜎𝜀2 ) に従う」という仮定に関しても問題がありそうです。図 5.8 で は，X 軸が 10 本くらいのところでは罹患している人としていない人がばらばらといるため，そ れなりの分散がありそうな一方で，0 本や 20 本のあたりでは全員が同じ値になっています。つ まり X 軸の値によって誤差分布の分散 𝜎𝜀2 が変動しているように思われるため，やはり直線を当 てはめるのはおかしそうです。 ということで，回帰分析モデルを拡張する必要があります。予測値の変動を説明変数の線形結 合で表した (5.6) 式の右辺はわかりやすいのでこのまま使うとして，直線ではなく確率を表す， つまり取りうる値が 0 から 1 に収まるように変換することを考えます。いくつか選択肢はある
5.3 一般化線形モデル 23 のですが，ここではロジスティック変換をしてみましょう。ロジスティック関数とは， 𝑔(𝑥) = 1 1 + exp(𝑥) (5.25) という関数です。図 5.9 は，𝑥 の値が −10 から 10 まで動いたときの 𝑔(𝑥) の値を表したもの です。このように，ロジスティック関数は 𝑥 の値が何であろうと確かに 𝑔(𝑥) は 0 から 1 までに 収まるようです。 1.00 予測確率 0.75 0.50 0.25 0.00 -10 -5 0 線形結合の和 5 10 図 5.9: ロジスティック変換 ということで，ロジスティック関数の 𝑥 を説明変数の線形和 𝑏0 +𝑥1 𝑏1 +𝑥2 𝑏2 +𝑥3 𝑏3 +⋯+𝑥𝐾 𝑏𝐾 に置き換えてあげれば， 𝑦= 1 1 + exp [−(𝑏0 + 𝑥1 𝑏1 + 𝑥2 𝑏2 + 𝑥3 𝑏3 + ⋯ + 𝑥𝐾 𝑏𝐾 )] (5.26) となり，この右辺 𝑦 は 0 から 1 までの値しか取らなくなります。また，式を変形させると log 𝑦 = 𝑏 0 + 𝑥 1 𝑏1 + 𝑥 2 𝑏2 + 𝑥 3 𝑏3 + ⋯ + 𝑥 𝐾 𝑏𝐾 1−𝑦 (5.27) という形になり，もとのきれいな右辺に戻すことも可能です。 またこの変換では（誤差）分散の不均一性も表現しています。というのは図 5.9 を見るとわか るように，確率が 0.5 付近では説明変数の線形和（X 座標）が少し変わるだけで予測確率が大き く変動するのに対して，確率が 0 または 1 付近の極端なエリアでは，ほとんど変動がなくなって います。 ということで，このように回帰分析モデルを変形させることで正規分布以外の分布にも対応可 能にしたものが一般化線形モデル (generalized linear model; GLM) です*10 。一般化線形 モデルでは，モデル式が 𝑔(𝑦) = 𝑏0 + 𝑥1 𝑏1 + 𝑥2 𝑏2 + 𝑥3 𝑏3 + ⋯ + 𝑥𝐾 𝑏𝐾 (5.28) という形になっています。つまり， 「𝑦 をなんか変換したもの」である 𝑔(𝑦) を 𝑥 の線形和 𝐗𝐛 で 予測しようという試みなわけです。この「𝑦 をなんか変換したもの」が，ロジスティック回帰の *10 これに対して正規分布の普通の重回帰モデルを一般線形モデル (general linear model) と呼んだりします。
5.3 一般化線形モデル 24 場合には 𝑔(𝑦) = log 𝑦 1−𝑦 というロジット変換*11 の式になっており，またふつうの重回帰モデルも 𝑔(𝑦) = 𝑦 という恒等関数によって変換されたものと見れば，これが一般化線形モデルの特殊なケースだと いうことがわかるでしょう。変換の関数 𝑔(𝑦) はリンク関数と呼びます。リンク関数は「二項分 布のときにはロジット変換」というように基本的には 𝑦 の分布に対してお決まりのものがありま す。一般化線形モデルでは変形の形ごとに異なる名称がついており，ロジスティック変換を用い た場合をロジスティック回帰と呼びます。もう少ししっかりとロジスティック回帰モデルの全体 を書くと， 𝑦𝑝 = 1 1 + exp [−(𝑏0 + 𝑥𝑝1 𝑏1 + 𝑥𝑝2 𝑏2 + 𝑥𝑝3 𝑏3 + ⋯ + 𝑥𝑝𝐾 𝑏𝐾 )] + 𝜀𝑝 , 𝜀𝑝 ∼ 𝑁 (0, 𝑦𝑝 (1 − 𝑦𝑝 )) (5.29) となります。つまり誤差分散は確率が 0.5 付近では最も大きく，0 や 1 付近では小さくなること を織り込んだ上でモデルのパラメタ（回帰係数）を推定しているのです。 INFO プロビット回帰 リンク関数を標準正規分布累積確率関数の逆関数 𝑔(𝑦) = Φ−1 (𝑦) と置いた GLM は，プロ ビット回帰と呼びます。ロジスティック回帰とプロビット回帰はどちらも「結果変数が 0 か 1 の二値」というケースに対して適用されるモデルです。そして 2 つの回帰で用いられ るリンク関数の逆関数はよく似た形（左右対称の釣鐘型）をしています。そんな 2 つの回 帰モデルについて，特に経済学など？ では理論的な背景なども考慮して厳密に使い分けら れる節があるようなのですが，この講義で紹介する手法（特に項目反応理論）では，基本 的にその区別には関心がありません。言ってしまえば「どっちでも良い」という感じがあ り，基本的には計算しやすいという理由からロジスティック回帰のほうを使用することが 多いです。 ということで，先程のタバコの例でもロジスティック回帰による回帰直線を当てはまると， 図 5.10 のようになります。 一般化線形モデルでは他にもポアソン回帰やガンマ回帰を始めとしていくつかありますが，こ こでは全ては紹介しません。以降の回では，ロジスティック回帰だけ理解しておけばついていけ るはずです。 *11 ロジスティック関数とロジット関数は互いに逆関数の関係になっています。そのため，説明変数側の変換として 見た場合にはロジスティック変換と呼ばれる一方で，結果変数側の変換（リンク関数）として見る場合にはロ ジット変換と呼ばれるわけです（たぶん） 。実際に回帰分析としての手法名では「ロジスティック回帰」とも「ロ ジット回帰」とも呼ばれることがあるはずです。
5.3 一般化線形モデル
25
肺がん罹患有無
1.00
0.75
0.50
0.25
0.00
0
10
20
1日あたりタバコの平均本数
30
図 5.10: ロジスティック回帰直線をあてはめる
5.3.1 R で一般化線形モデル
R で一般化線形モデルを行う場合は，lm() の代わりに glm() という関数を使います。基本的
には使い方は同じですが，GLM なので「誤差分散の分布の形」および「リンク関数」を指定し
てあげる必要があります。指定の仕方にちょっとクセがありますが，そういうものとして理解し
てください。
GLM（ロジスティック回帰）
1
2
3
4
5
6
7
1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 結果変数は二値変数でないといけないので，「Q1_Aが平均以上か」を表す変数を作成
# datには追加していない
Q1_A_binom <- dat$Q1_A > mean(dat$Q1_A)
# 実は引数dataの中に無い変数でも分析に利用できる
# ↑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
5.3 一般化線形モデル
15
16
17
18
19
20
21
26
(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
result_glm <- glm(Q1_A ~ age + gender, data = dat, family = gaussian(link =
2
"identity"))
summary(result_glm)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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
on 2431
on 2429
degrees of freedom
degrees of freedom
5.4 変数間の関係のグラフィカルモデル 19 20 21 AIC: 14097 Number of Fisher Scoring iterations: 2 当然ですが，lm() のときと同じ結果が得られます。 5.4 変数間の関係のグラフィカルモデル 最後に，今後の授業で登場するグラフィカルモデル（あるいは有向非巡回グラフ; Directed acyclic graph [DAG]）の導入をしておきます。といっても難しい事はありません。単に説明変数 と被説明変数の関係性を矢印で表すことで，モデルの見通しを良くするために使うツールです。 回帰分析の例では，gender および age によって Q1_A の値を予測してみました。これをグラ フィカルモデルで表すと図 5.11 のようになります。 図 5.11: 重回帰モデル このように，グラフィカルモデルでは観測可能な変数を長方形で表し，矢印によって回帰の向 きを表します。また，両方向の矢印は変数間の相関関係を表しています。重回帰分析では，複数 の説明変数が一つの結果変数を説明しようとするため，グラフィカルモデルはこのように複数の 変数から一つの変数に向かって矢印がひかれる形になります。 グラフィカルモデルを自由に拡張して，例えば 図 5.12 のように変数間の関係を広げていった ものが，パス解析と呼ばれるやつです。このモデルを分析する場合，矢印一個一個を回帰分析し ても良いかもしれませんが，全体として変数間のモデルができているのですから，全部まとめ てデータとの適合度などを評価してあげたいところです。ということでパス解析も SEM の一種 です。 SEM（第 7 章）では，このような変数間のグラフィカルモデルをもとに分析用のコードを作 成していくことになります。といっても難しいことはないはずです。データを取る時点で変数間 の理論的な関係性の仮説は考えているはずですから… 27
5.4 変数間の関係のグラフィカルモデル 図 5.12: パス解析 28
参考文献 参考文献 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 29