多変量解析_05_回帰分析のおさらい

5.1 回帰分析とは Chapter 5 回帰分析のおさらい 本資料は，神戸大学経営学研究科で 2022 年度より担当している「統計的方法論特殊研究（多 変量解析）」の講義資料です。CC BY-NC 4.0 ライセンスの下に提供されています。 作成者連絡先 神戸大学大学院経営学研究科 分寺杏介（ぶんじ・きょうすけ） mail: [email protected] website: https://www2.kobe-u.ac.jp/~bunji/ 統計を学ぶ人の大多数が最初に出会う「多変量解析」は，きっと（重）回帰分析でしょう。実 際に，後半で登場する因子分析・構造方程式モデリング・項目反応理論はいずれも回帰分析 のモデルと関係のあるモデルです。一応本講義では回帰分析については既習であることを前 提としているので，基本的な回帰分析のおさらいをはさみ，後半の内容に必要な部分だけで も前提知識を補いたいと思います。なお，この Chapter では，基本的に今後紹介する手法の 理解に必要・重要と思われる内容に絞って紹介しています。回帰分析は相当幅広く用いられ ているため，本気を出せば回帰分析だけで 1 コマ 15 回の講義を余裕で組めると思います。回 帰分析を深く理解したい方は，恐縮ながら経済学研究科や国際協力研究科などで開講されて いる講義も探してみてください。 Contents 5.1 5.2 5.3 5.4 回帰分析とは . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 重回帰分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 一般化線形モデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 変数間の関係のグラフィカルモデル . . . . . . . . . . . . . . . . . . . . . 1 6 20 26 5.1 回帰分析とは 中学数学では，図5.1のように「2 つの点を通る直線を求める」問題を解いたことがあると思 います。回帰分析はこの延長にある分析手法と言っても良いでしょう。 この問題を（今後の問題設定に合わせて）数式を用いて表現すると， 1

• https://creativecommons.org/licenses/by-nc/4.0/
• https://www2.kobe-u.ac.jp/~bunji/

5.1 回帰分析とは 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 つの点に 関してはズレがゼロですが，残り 1 個の点に対しては大きくズレています。この赤い線のよ 2

5.1 回帰分析とは 4 3 2 y 1 0 -1 -2 -3 -4 -3 -2 -1 0 x 1 2 3 2 3 図5.2 3 つの点を通る直線を求める 4 3 2 y 1 0 -1 -2 -3 -4 -3 -2 -1 0 x 1 図5.3 3 つの点を通る直線なんてない うに，平均的なズレを最小化する直線を求める分析法が回帰分析なのです。 点の数が何個であろうと最もベーシックな回帰分析では，直線（一次関数）を求めます。傾 きを 𝑏1 ，切片を 𝑏0 とおけば 𝑦 = 𝑏 0 + 𝑏1 𝑥 (5.1) という直線について，𝑏0 , 𝑏1 の値がいくつのときに最もそれっぽい直線になるかを求めるわ けです。 5.1.1 最小二乗法 実際の分析においては，最小二乗法によって回帰直線を求めることが多いです。これは 図5.4の青い矢印のように，縦方向でのズレを求めて，これの二乗和が最小になる (𝑏1 , 𝑏0 ) の 組み合わせを求める方法です。 「それっぽい直線」という観点で言えば，図5.5の青い矢印のように，赤い直線との直線距離 に関して最小化しても良さそうな気がしますが，最小二乗法では縦方向での距離を最小化し ています。なぜなのでしょうか。 それは，回帰分析が 𝑥 で 𝑦 を説明する分析手法だからです。𝑦 そのものは無いけれど他の変 3

5.1 回帰分析とは 1 y 0 -1 -2 -3 -1 0 x 図5.4 1 最小二乗法 1 y 0 -1 -2 -3 -1 0 x 1 図5.5 直線距離で考えたらダメなのかい 数がある時に𝑦 を近似するのが回帰直線なのです。例えば，𝑥 を「体重」，𝑦 を「身長」だと しましょう。あなたのもとに，あるゲストがやってきます。どういうわけかその人に関して 「体重」の情報だけが入っているとします。あなたはゲストのためにちょうどよい高さの椅子 を用意する必要があります。しかし，ちょうどよい高さの椅子をあてがうためには，ゲスト の身長が分かっていないといけません*1 。さて，あなたはゲストの身長をどのように見積も るでしょうか。 普通ならば，図5.6のように，周りの人の体型などから「 （標準体型で）体重が 70kg くらいの 人ならば，身長はだいたい 170cm くらい」「（標準体型で）体重が 60kg くらいの人ならば， 身長はだいたい 162cm くらい」という感じで予測を立てたうえで，ゲストの体重の情報を当 てはめて「ゲストの体重が 65kg ということは，標準体型ならば身長はだいたい 165.5m く らいかなぁ」といった予測を行うでしょう。回帰分析で行っていることは概ねこんな感じの ことですが，ここで問題になる「予測のズレ」は，身長（つまり 𝑦 ）に関してのみです。つま り，𝑥 による 𝑦 の近似（予測）のズレが最小になるように回帰直線を引きたいがために，最 小二乗法では縦方向でのズレのみを考えているわけですね。 *1 4 もちろん実際には股下やら体型やらによって最適な椅子の高さは決まるわけですが，ここでは簡略化して「最 適な椅子の高さは身長のみで決まる」という状況を考えてみます。

5.1 回帰分析とは 身長 180 160 140 40 50 60 体重 70 80 図5.6 回帰直線による予測 求めたい回帰直線は 𝑦 = 𝑏0 + 𝑏1 𝑥 ですが，一つ一つの点について考えると，ほとんどの場合 等号は成立しません。というのも回帰直線が示すのはあくまでも予測値であって，実際のデ ータでは予測とのズレがあるためです。そこで回帰分析モデルでは，実際の一つ一つのデー タの値 𝑦𝑝 は，回帰直線による予測値 𝑦𝑝̂ と誤差 𝜀𝑝 の和であると考えます。ということで，回 帰分析モデルの全体は 𝑦𝑝 = 𝑦𝑝̂ + 𝜀𝑝 (5.2) 𝑦𝑝̂ = 𝑏0 + 𝑏1 𝑥𝑝 (5.3) 𝜀𝑝 ∼ 𝑁 (0, 𝜎𝜀 ) (5.4) となり，𝑃 をデータの総数とすると，最小二乗法では 𝑃 𝑃 ∑(𝑦𝑝 − 𝑦𝑝̂ )2 = ∑ [𝑦𝑝 − (𝑏0 + 𝑏1 𝑥𝑝 )] 𝑝=1 2 (5.5) 𝑝=1 を最小化するような (𝑏0 , 𝑏1 ) の組み合わせを求めている，というわけです。 最尤法 最小二乗法と同じくらい用いられる推定法が最尤法です。(5.2)-(5.4) 式の回帰分析モ デルの式を変形させると， 𝑦𝑝 ∼ 𝑁 (𝑏0 + 𝑏1 𝑥𝑝 , 𝜎𝜀 ) という形にもなります。したがって，データから最も尤度が大きくなる (𝑏0 , 𝑏1 , 𝜎𝜀 ) の 組を見つけてあげることで，回帰直線を最尤推定することも可能なわけです。 最小二乗法とどちらが良いかという話ではなく，最適化の目的関数が異なるというこ となので，回帰直線を求める方法は最小二乗法だけじゃないということは頭の片隅に 入れておいても良い知識だと思います。 5.1.2 （単）回帰分析の線形代数的表現 因子分析や SEM のところでは，それぞれの分析手法の基本的な考え方を説明するために線 形代数的表現を多用しています。ということで，ここで回帰分析の線形代数による表現を説 明しておきます。書き方とその意味に慣れておいてください。 5

5.2 重回帰分析 まず，回帰分析モデルを，図5.1のところで示した連立方程式的な書き方で書き直すと，次の ようになります。 (𝑥1 , 𝑦1 ), (𝑥2 , 𝑦2 ), ⋯ , (𝑥𝑃 −1 , 𝑦𝑃 −1 ), (𝑥𝑃 , 𝑦𝑃 ) の 𝑃 個の点に関する以下の連立方程式につ 𝑃 いて，誤差の和 ∑𝑝=1 𝜀2𝑝 が最小になる係数 𝑏0 , 𝑏1 の値はそれぞれいくつか。 ⎧ 𝑦1 { { 𝑦2 { ⎨ {𝑦𝑃 −1 { { ⎩ 𝑦𝑃 = 𝑏0 + 𝑏1 𝑥1 + 𝜀1 = 𝑏0 + 𝑏1 𝑥2 + 𝜀2 ⋮ = 𝑏0 + 𝑏1 𝑥𝑃 −1 + 𝜀𝑃 −1 = 𝑏0 + 𝑏1 𝑥𝑃 + 𝜀𝑃 毎回全部式を書いていると面倒なので，ここで線形代数の力を借ります。連立方程式のうち， すべての式に共通しているのが 𝑏0 , 𝑏1 の部分です。一旦ここはそのままにしておいて，その 他の（𝑝 ごとに異なる）部分をベクトルでまとて書き直してみます。 𝜀1 𝑥1 ⎧ 𝑦1 {⎡ 𝑦 ⎤ ⎡ 𝑥 ⎤ ⎡ 𝜀 ⎤ {⎢ 2 ⎥ ⎢ 2 ⎥ ⎢ 2 ⎥ ⋮ ⎥ = 𝑏 0 + 𝑏1 ⎢ ⋮ ⎥ + ⎢ ⋮ ⎥ ⎨⎢ ⎢𝑥𝑃 −1 ⎥ ⎢𝜀𝑃 −1 ⎥ {⎢𝑦𝑃 −1 ⎥ { ⎩⎣ 𝑦𝑃 ⎦ ⎣ 𝑥𝑃 ⎦ ⎣ 𝜀𝑃 ⎦ だいぶ見通しが良くなってきました。あとはベクトルを太字で表すことで，以下のように割 と綺麗に書き直すことが出来ます。 2 つの変数に関する長さ 𝑃 のベクトル 𝐱 [𝑦1 [𝜀1 𝑦2 ⋯ 𝑦𝑃 ] ⊤ = [ 𝑥1 ⊤ 𝑥2 ⋯ 𝑥𝑃 ] , 𝐲 = が 与 え ら れ た と き， 以 下 の 式 に お け る 誤 差 ベ ク ト ル 𝛆 = ⊤ 𝜀2 ⋯ 𝜀𝑃 ] の二乗和 𝛆⊤ 𝛆 が最小になる係数 𝑏0 , 𝑏1 の値はそれぞれいくつか。 𝐲 = 𝑏 0 + 𝑏1 𝐱 + 𝛆 ということで，これでデータが何個あっても一つの式で表すことができるようになりました。 便利ですね。 5.2 重回帰分析 𝑦 を予測するという意味では，もちろん変数は多い方が良いでしょう。「身長」を予測するう えで，体重だけでなく「性別」や「体脂肪率」といったことがわかれば，予測はより正確に なるはずです。回帰分析でも，説明変数を複数個用いることが出来ます。これを一般的には 重回帰分析と呼びます。説明変数が 𝐾 個あるならば，対応する回帰係数を 𝐾 + 1（切片の 分）個用意して，以下のような直線を考えるわけです*2 。 𝑦 = 𝑏0 + 𝑏1 𝑥1 + 𝑏2 𝑥2 + ⋯ + 𝑏𝑘 𝑥𝑘 + ⋯ + 𝑏𝐾 𝑥𝐾 *2 6 (5.6) 直線がイメージしづらいかもしれませんが，説明変数が 1 個の場合の回帰直線が 2 次元平面に引かれていた ように，説明変数が 2 個の場合は 3 次元空間の中に一本の直線が引かれます。同様に，説明変数が 𝐾 個なら ば 𝐾 + 1 次元空間に一本の直線を引くわけです。

5.2 重回帰分析 ここでのポイントは，重回帰分析では変数の効果が線形和（重み付け和）になっているとい う点です。回帰係数は「その説明変数の値が 1 大きくなったら 𝑦 の予測値がどれだけ大きく なるか」を示していますが，最もベーシックな重回帰分析モデルの場合 • 説明変数の値がいくつであってもその効果は一定 • 他の説明変数の値は完全に無関係 であることが仮定されています。もし高次の効果（e.g., 薬は適量ならば飲むほど健康に良い が，飲み過ぎるとかえって良くない＝二次の効果）や交互作用を考慮したい場合には，直接 それを表す項を追加する必要があります。例えば 𝑦 = 𝑏0 + 𝑥1 𝑏1 + 𝑥2 𝑏2 + 𝑥21 𝑏3 + 𝑥1 𝑥2 𝑏4 といった感じです。 また，重回帰分析で得られる回帰係数は，厳密には偏回帰係数と呼ばれます。これは「他の 説明変数の値が固定されているときに，その説明変数の値が 1 大きくなったら 𝑦 の予測値が どれだけ大きくなるか」を表しています。もう少し正確に言うと，その説明変数以外で 𝑦 を 予測した際の誤差に対して，その説明変数のみで単回帰分析をした際の回帰係数が偏回帰係 数です。例えば 2 変数 𝑥1 , 𝑥2 の重回帰における 𝑥1 の偏回帰係数 𝑏1 は， 𝑦 = 𝑏01 + 𝑏2 𝑥2 + 𝜀1 (Step1) 𝜀1 = 𝑏02 + 𝑏1 𝑥1 + 𝜀2 (Step2) という二段階で推定することで，単回帰の組み合わせに分解することができるのです（Step 1 の左辺は 𝑦 ̂ ではなく 𝑦 である点に注意）。これにより，偏回帰係数は他の変数の影響を完全 に無視した上でのある変数の影響の大きさを見積もることができるわけですね。 5.2.1 （おまけ）ダミー変数 回帰係数は「その説明変数の値が 1 大きくなったら 𝑦 の予測値がどれだけ大きくなるか」を 示しています。では，説明変数が名義変数の場合はどうしたら良いでしょうか。まずはカテ ゴリ数が 2 の場合（e.g., 生物学的な性：男女）を考えてみます。通常，このようなデータで は一方のカテゴリに 1 を，もう一方のカテゴリに 0 をあてます。今回は「男性は 0，女性は 1」としておきましょう。この「性別 (𝑥𝑝2 )」と「体重 (𝑥𝑝1 )」から「身長 (𝑦𝑝 )」を予測する重 回帰モデルを立てると 𝑦𝑝̂ = 𝑥𝑝1 𝑏1 + 𝑥𝑝2 𝑏2 + 𝑏0 (5.7) となります。このモデルを男性と女性それぞれに分けてみてみると，男性は 𝑥2 = 0，女性は 𝑥2 = 1 なので + 𝑏0 (男性) 𝑦𝑝̂ = 𝑥𝑝1 𝑏1 + 𝑏2 + 𝑏0 (女性) 𝑦𝑝̂ = 𝑥𝑝1 𝑏1 (5.8) となります。したがって，男性モデルでは切片が 𝑏0 となっており，女性モデルでは切片が 𝑏2 + 𝑏0 となっているわけです。回帰係数の解釈通り，𝑏2 は 𝑥2 の値が 1 大きくなった際に 𝑦 の予測値がどれだけ大きくなるかを意味するわけですが，より正確に言えば体重𝑥1𝑝 の値が 同じ男性と女性がいた場合，身長は平均的にどれだけ異なるかを表しているわけですね。 7

5.2 重回帰分析 カテゴリ数が 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.9) (𝐴𝐵) となります。つまりA型-B型の差とB型-O型の差とO型-AB型の差がすべて同じということ になってしまい，明らかにおかしなことになっています。このように，名義変数ではカテゴ リ数が 2 までならば対応可能ということで，カテゴリ数が多い場合にはそのまま数値化して 投入するのではなく，二値のダミー変数化することで対応します。 血液型の場合，4 カテゴリなので 3 つのダミー変数 (𝑥2𝐴 , 𝑥2𝐵 , 𝑥2𝑂 ) および対応する回帰係数 (𝑏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.10) (𝐴𝐵) と書くことができ，各回帰係数はAB型との差として解釈可能になります。ここでのポイント は，ダミー変数はカテゴリ数 −1 個だという点です。なぜカテゴリ数ぶんのダミー変数がある と困るのか考えてみましょう。例えば上の 3 つのダミー変数に加えて「AB 型なら 1」となる 𝑥24 および回帰係数 𝑏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 (𝐴𝐵) 8 (5.11)

5.2 重回帰分析 となります。ですが，この式には別の解が考えられます。例えば (𝑏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.12) となり，結果的に先程のモデルと全く同じになってしまいました。このように，回帰分析で はダミー変数がカテゴリ数の数だけあると解が一意に定まらなくなってしまうため，カテゴ リ数-1 個だけ用意する必要があるのです*3 。 ダミー変数がカテゴリ数より少ないとなると，必然的にダミー変数が与えられない（全ての ダミー変数が 0 であることによって表される）カテゴリが発生します。このカテゴリのこと を基準カテゴリや参照カテゴリと呼んだりしますが，この基準カテゴリの決め方も重要にな ります。というのも回帰係数は基準カテゴリとの差であるためです。回帰分析では，回帰係 数の検定（𝐻0 ∶ 𝑏 = 0）が行われるため， 「あるカテゴリが基準カテゴリと平均的に差がある か」は検証できる一方で，基準カテゴリ以外のカテゴリ同士での差の検定はそのままでは出 来ません。実験的な環境であれば対照群を基準カテゴリにおくのが良いでしょう。 5.2.2 重回帰分析の線形代数的表現 説明変数の数が増えても，同様に線形代数による表現を利用します。(5.6) 式に基づく予測 を，まず 𝑝 番目のデータ一つについて考えてみます。説明変数が何個になっても，観測値 𝑦𝑝 を予測値 𝑦𝑝̂ と誤差 𝜀𝑝 に分ける点 (5.2式) は変わりません。 𝑦𝑝 = 𝑦𝑝̂ + 𝜀𝑝 = 𝑏0 + 𝑏1 𝑥𝑝1 + 𝑏2 𝑥𝑝2 + ⋯ + 𝑏𝐾 𝑥𝑝𝐾 + 𝜀𝑝 (5.13) これ以降，説明変数には 2 つの添字が付きます。2 つ目の添字 𝑘(𝑘 = 1, 2, ⋯ , 𝐾) は，何番目 の説明変数であるかを表すものです。例えば 𝑥𝑝3 （または 𝑥𝑝,3 ）は「𝑝 番目の人の 3 番目の 説明変数の値」（dat[p,3]）を指しています。添字の順序は，Rでデータフレームの行と列 の要素を指定するときの順番（apply() 関数の引数 MARGIN の指定）だと覚えると良いかも しれません（個人差アリ）。 まずは (5.13) 式の 𝑏1 𝑥𝑝1 + 𝑏2 𝑥𝑝2 + ⋯ + 𝑏𝐾 𝑥𝑝𝐾 （切片と誤差項以外）の部分について，ベ クトルの積を使って簡略化していきます。これは 𝑝 番目のデータにおける説明変数ベクトル 𝐱𝑝 = [𝑥𝑝1 , 𝑥𝑝2 , ⋯ , 𝑥𝑝𝐾 ] と，回帰係数ベクトル 𝐛 = [𝑏1 , 𝑏2 , ⋯ , 𝑏𝐾 ]⊤ を使うと， 𝑏1 ⎡𝑏 ⎤ 𝐱𝑝 𝐛 = [𝑥𝑝1 𝑥𝑝2 ⋯ 𝑥𝑝𝐾 ] ⎢ 2 ⎥ ⎢ ⋮ ⎥ ⎣𝑏𝐾 ⎦ = 𝑥𝑝1 𝑏1 + 𝑥𝑝2 𝑏2 + ⋯ + 𝑥𝑝𝐾 𝑏𝐾 *3 9 (5.14) ちなみに機械学習では one-hot encoding という名で使われるケースもあります。推定できなくなってしま うのは，説明変数の効果が線形和で表される回帰モデルならではの問題なのです。

5.2 重回帰分析 と，かなりシンプルに表すことが出来ました。ここでもう少し工夫してみます。せっかくな ので，切片についても同じようにまとめてしまいましょう。(5.13) 式の 𝑏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.15) = 𝑦𝑝̂ というように，切片項まで含めて回帰直線による予測値 𝑦𝑝̂ を簡略化できました。 あとはこれを 𝑃 個のデータ全体に拡張するだけです。といっても考え方は簡単で， 𝑃 ×(𝐾 +1) の説明変数全体（+ 切片）を 𝐱1 1 𝑥11 ⎡ 𝐱 ⎤ ⎡1 𝑥 2 21 𝐗=⎢ ⎥=⎢ ⎢ ⋮ ⎥ ⎢⋮ ⋮ 𝐱 1 𝑥 ⎣ 𝑃⎦ ⎣ 𝑃1 𝑥12 𝑥22 ⋮ 𝑥𝑃 2 ⋯ 𝑥1𝐾 ⋯ 𝑥2𝐾 ⎤ ⎥ ⋱ ⋮ ⎥ ⋯ 𝑥𝑃 𝐾 ⎦ と表せば，これを用いてすべてのデータに関する予測値 𝐲̂ を 𝑦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 𝑥11 ⎡1 𝑥 21 =⎢ ⎢⋮ ⋮ ⎣1 𝑥𝑃 1 𝑥12 𝑥22 ⋮ 𝑥𝑃 2 とまとめることが可能です。あとは誤差ベクトルを付け加えたら，以下のようにスッキリと した表現の出来上がりです。 2 つの変数に関する長さ 𝑃 の結果変数ベクトル 𝐲 と 𝑃 × (𝐾 + 1) の説明変数行列 𝐗 が 与えられたとき，以下の式における誤差ベクトル 𝛆 の二乗和 𝛆⊤ 𝛆 が最小になる回帰係数 ベクトル 𝐛 の値はいくつか。 𝐲 = 𝐗𝐛 + 𝛆 因子分析や SEM のところでは，このような感じで「変数の行列」×「係数のベクトル（また は行列） 」の積を用いて手法のコンセプトを紹介していきます。ということで，今のところは ここで紹介した記法に慣れておいてください。 10

5.2 重回帰分析 正規方程式 上記のように重回帰分析を線形代数で表現すると，最小二乗法の解も解析的に求める ことが出来ます。いま，求めたいのは 𝛆⊤ 𝛆 が最小になるような 𝐛 の値でした。そこ で，𝛆 を (𝐲 − 𝐗𝐛) に戻して考えると，𝛆⊤ 𝛆 = (𝐲 − 𝐗𝐛)⊤ (𝐲 − 𝐗𝐛) が最小になるよ うな 𝐛 を求めたら良いわけです。 最小値になる 𝐛 を求めたいわけなので，ちょっと大変かもしれませんが，行列を展開 して 𝐛 で偏微分したものイコール 0 という式をとればa ， 𝐗⊤ 𝐗𝐛 = 𝐗⊤ 𝐲 というかたちになります。この方程式のことを正規方程式と呼びます。あとは 𝐗⊤ 𝐗 が正則であれば，両辺に逆行列をかけることで，この式を 𝐛 = の形に変形させ， 𝐛 = (𝐗⊤ 𝐗)−1 𝐗⊤ 𝐲 という解が得られます。 ちなみに「𝐗⊤ 𝐗 が正則でない」ときには逆行列が求められないため，解を得ること が出来ません。このような状態になるのは，行列がランク落ちしている場合ですが， ランク落ちが発生するのは 𝐗 の特定の変数が線形従属になっている場合，つまり他 の変数の線形和で表せてしまう場合です。回帰分析ではよく「多重共線性に気をつけ ましょう」という話がありますが，この理由として，完全な多重共線性がある場合に は 𝐗⊤ 𝐗 が正則でなくなってしまうために，解が求められなくなってしまうから，と いう説明も可能なのです。 a このあたりの展開は「正規方程式」で検索したらすぐ見つかります 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 と同じですが，両辺をチルダ (~) でつなぐのが特徴です。とりあえずこの授業の中ではチルダで結んだら回帰ということ 11

5.2 重回帰分析 を覚えておきましょう。実は SEM で使用する lavaan パッケージの記法でもチルダは回帰 を表します。 formula で使える記号 基本的には + で説明変数を入れていけば良いのですが，他にもいくつかの記号があり ます。 • 2 つの説明変数を: でつなぐと，その変数間の交互作用項を指定できます。例 えば x1 と x2 があるときに，y ~ x1:x2 とすると，これは 𝑦 = 𝑥1 𝑥2 𝑏12 + 𝑏0 という回帰式を表します。 • 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 Call: 3 lm(formula = Q1_A ~ age + gender, data = dat) 4 5 Coefficients: 6 (Intercept) age gender 7 18.05453 0.07083 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 ならば何でも良い*4 ことになりま *4 12 回帰係数の検定という観点では差が 1 である必要すらも無いのですが，回帰係数の解釈を考えるとやはり 0/1 のコーディングが一番しっくりきます。

5.2 重回帰分析 す。もし dat での gender の値が全員 1 ずつ小さい値で（つまり男性は 0，女性は 1 と）コ ーディングされていたとしても，gender に対する回帰係数は変わりません。少し切片が変 わって 𝑦𝑝̂ = (18.05453 + 1.89117) + 0.07083𝑥𝑝,age + 1.89117𝑥𝑝,gender = 19.9457 + 0.07083𝑥𝑝,age + 1.89117𝑥𝑝,gender となります。男女ごとに回帰式を見ると 𝑦𝑃̂ = { 19.9457 + 0.07083𝑥𝑝,age (男性) 19.9457 + 0.07083𝑥𝑝,age + 1.89117 (女性) となるわけです。 切片の役割 切片は「すべての説明変数の値が 0 のときの予測値」を表しています。したがって dat で言えば「age が 0 ＝ 0 歳で，gender が 0 の人」の予測値になりますが当然そ んな人はいません。このように，多くの場合では切片の値それ自体は使えない（検定 する意味もない）ものです。 もし切片に役割を持たせたい場合には，すべての量的変数を中心化すると良いです。 中心化，つまり平均値を引くことで各変数の値が 0 というのが「その変数が平均値の 人」を表すようになります。つまり全ての量的変数を中心化した状態で行った回帰分 析の切片は「すべての変数が平均値の人」の予測値になるわけです。 なお，gender のようなカテゴリカル変数は中心化してもしなくても良いです。例え ば男女比が 1:1 のときに中心化を行うと「男性は-0.5，女性は 0.5」という感じになり ます。こうなるとむしろ 0 が意味を持たなくなってしまうので，あえて中心化しない ことで，切片項は「量的変数はすべて平均値の男性での平均値」を意味するようにな ります。一方で，カテゴリカル変数も全て中心化した場合，切片項は純粋に全体平均 値を表すようになります。ということでどちらにせよ解釈可能な意味を持つので，用 途に応じて中心化するかを決めるのが良いでしょう。 とりあえず回帰分析ができたので，次は回帰係数の検定をしてみます。といっても lm() に は既に回帰係数の検定が用意されているので難しいことはありません。 回帰係数の検定 1 # まずlm()の結果をオブジェクトに代入 2 result_lm <- lm(Q1_A ~ age + gender, data = dat) 3 summary(result_lm) R では，関数の出力を含むあらゆるものをオブジェクト（変数）に代入することが出来ます。 これにより，例えば「因子分析を行い，結果から因子負荷行列を取り出し，その行ごとの和 をとる」みたいなことが全てコードとして書けるようになるわけです。同様に，特定の関数 の出力を受け取って別の処理を行うような関数もたくさん存在しています。ここでは関数の 出力自体も変数として扱えるという感覚を理解しましょう。 13

5.2 重回帰分析 何かしらの分析を行う関数の出力オブジェクトは，たいてい list 型で与えられます。list 型はほぼ何でも，いくつ入れても良い型です。例えば list 型である result_lm（lm() の 出力）の中身を見てみると，以下のような状態になっています。 list の中身 1 # 長いので省略しています 2 str(result_lm) 1 List of 12 2 3 4 $ 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 ... 5 6 ..- 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 ... 7 ..- attr(*, "names")= chr [1:2432] "(Intercept)" "age" "gender" "" ... 8 $ rank 9 $ fitted.values: Named num [1:2432] 21.1 23.1 23 23 21.1 ... 10 : int 3 ..- attr(*, "names")= chr [1:2432] "1" "2" "3" "4" ... 11 $ assign : int [1:3] 0 1 2 12 $ qr :List of 5 13 ..$ qr 14 .. ..- attr(*, "dimnames")=List of 2 15 .. .. ..$ : chr [1:2432] "1" "2" "3" "4" ... 16 .. .. ..$ : chr [1:3] "(Intercept)" "age" "gender" 17 .. ..- attr(*, "assign")= int [1:3] 0 1 2 18 ..$ qraux: num [1:3] 1.02 1.02 1.02 19 ..$ pivot: int [1:3] 1 2 3 20 ..$ tol 21 ..$ rank : int 3 22 ..- attr(*, "class")= chr "qr" 23 24 : num [1:2432, 1:3] -49.3153 0.0203 0.0203 0.0203 0.0203 ... : num 1e-07 [list output truncated] - attr(*, "class")= chr "lm" 一番上に List of 12 とあることから，このオブジェクトはいろいろなものが 12 個くっつ いている，ということがわかります。一つ目が coefficients で，これは長さ 3 の名前付き 実数型 (Named num) のベクトルだということです。その下の residuals や effect は同じ く名前付き実数型ベクトルですが長さは 2314 のようです。また，少し下にある qr は List of 5 となっています。つまり list の中に list が入れ子になっている，ということなので す。このように，list 型オブジェクトは型の違いも気にせずにいろいろなものをくっつけて 一つのオブジェクトとして扱うことが出来ます。 list 型オブジェクトの中身は，それぞれの要素名の前に $ がついていることから分かるよ 14

5.2 重回帰分析
うに $ 記号を使って取り出すことが出来ます*5 。例えば以下のコードでは，result_lm の中
から回帰係数の部分 (coefficients) だけをベクトルとして取り出すことができます。この
行為に意味があるかはともかく，大抵のオブジェクトでは分析結果から中身を取り出せると
いう感覚だけでも理解しておいてください。コレ，分析の手順の自動化率を高める際などに
かなり重要になってきます。
結果の一部だけを取り出す
1

result_lm$coefficients

1

(Intercept)

age

gender

2

18.05453074

0.07082883

1.89116882

今回は lm() 関数の出力を result_lm というオブジェクトに代入しました。そしてこれに対
して summary() 関数を実行しました。Chapter 2 では summary() 関数のことを「要約統計
量を出してくれる関数」と紹介しましたが，実際には summary() 関数はもっと広く「与え
られたオブジェクトの型に合わせた何かしらのサマリーを出してくれる」という関数なので
す*6 。
1
2

Call:

3

lm(formula = Q1_A ~ age + gender, data = dat)

4
5

Residuals:

6

Min

1Q

Median

3Q

Max

7

-19.1742

-2.6150

0.6091

3.1820

8.8502

8
9

Coefficients:

10

Estimate Std. Error t value Pr(>|t|)

11

(Intercept) 18.054531

0.394846

45.726

<2e-16 ***

12

age

0.070829

0.008116

8.727

<2e-16 ***

13

gender

1.891169

0.189308

9.990

<2e-16 ***

14

---

15

Signif. codes:

0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

16
17

Residual standard error: 4.386 on 2429 degrees of freedom

18

Multiple R-squared:

19

F-statistic: 91.56 on 2 and 2429 DF,

0.07011,

Adjusted R-squared:

0.06934

p-value: < 2.2e-16

出力は上から
*5

データフレームの列も $ で取り出せるということは，データフレームは全ての要素が「同じ長さのベクトル」
のリストである，という見方もできるわけです。リストとして扱う必然性は無いですが。
*6 内部的には，summary() で呼び出せる複数の関数が存在しており，与えるオブジェクトのクラスによって異な
る関数が実行されています。lm() の出力は lm クラスのオブジェクトであり，これに対しては summary.lm()
が自動的に呼び出されているのです。なので lm オブジェクトに対する summary() 関数のヘルプを見たい場
合には?summary.lm としてあげる必要があります。

15

5.2 重回帰分析 Residuals 𝑦 − 𝑦 ̂ の要約統計量です。回帰分析の前提条件として，誤差 𝜀 は正規分布 𝑁 (0, 𝜎𝜀2 ) に従う必要があるため，特に中央値が 0 から大きく離れていたら注意が必要 かもしれません。 Coefficients 回帰係数です。左から「係数の推定値」 「標準誤差」 「𝑡 統計量」 「𝑝 値」が表 示されています。今回の場合，いずれも 𝑝 < .05 なので，統計的に有意ということで すね。 Multiple R-squared 重相関係数です。この後説明します。 Adjusted R-squared 自由度調整済み重相関係数です。この後説明します。 F-statistic 「切片以外の回帰係数が全て 0 である」という帰無仮説に対する 𝐹 統計量で す。ということはこの検定が有意ではなかった場合，モデルとしての意味が無いこと になります。 5.2.4 決定係数・重相関係数 summary(lm()) で出てくる Multiple R-squared は決定係数と呼ばれます。回帰分析モデ ルでは，結果変数 𝑦 は予測値 𝑦 ̂ と誤差 𝜀 によって 𝑦 = 𝑦̂ + 𝜀 (5.16) と表されます。当然ながら，予測値と誤差の間は無関係であるはずなので*7 ，𝑦 の分散は 𝜎𝑦2 = 𝜎𝑦2̂ + 𝜎𝜀2 (5.17) というように，予測値と測定誤差の分散の単純な和として表すことができます。信頼性係数 のときと同じ感じですね。このとき，決定係数 𝑅2 は 𝑅2 = 𝜎𝑦2̂ 𝜎𝑦2 =1− 𝜎𝜀2 𝜎𝑦2 (5.18) と定義されます。つまり，決定係数 𝑅2 は結果変数の分散のうち，説明変数で説明可能な分 散の割合ということです。 ここで，𝑦 に影響を与えるすべての説明変数が分かっている状況を考えてみましょう。例え ば Q1_A の値を決定づけるものは，幼いときの親の関わり方や小学校での友人関係，あるい は昨日の晩御飯や体温，回答の瞬間の気分なんかも影響するかもしれません。ほぼゼロだと してもゼロではない限りは重回帰モデルに含める意味があります。仮に説明変数がちょうど 1 億個で 𝑦 を完全に説明可能だったとすると，重回帰モデルは，以下のようになります（左 辺が 𝑦 ̂ ではなく 𝑦 である点に注意） 。 𝑦𝑝 = 𝑏0 + 𝑥𝑝,1 𝑏1 + 𝑥𝑝,2 𝑏2 + 𝑥𝑝,3 𝑏3 + ⋯ + 𝑥𝑝,99999999 𝑏99999999 + 𝑥𝑝,100000000 𝑏100000000 ただ，実際問題として 1 億個の説明変数による回帰モデルが判明したとしても，このままで は変数の数が多すぎて実社会には応用できません。1 億個あっても，大半はほぼ無視できる *7 16 回帰分析モデルでは，これも仮定のひとつとしておかれます。この仮定が守られていない場合，検定結果に バイアスがかかったりするわけです。

5.2 重回帰分析 レベルの影響しか与えていないでしょうし，できることならなるべく少ない変数の数で十分 な予測モデルを立てることができれば，みんな嬉しいはずです（オッカムの剃刀）。そこで， この 1 億個の変数の中から，最も説明力が高い 2 つだけを選んできたとしましょう。今，1 億個の説明変数が「説明力の高い順」に並んでいたとすると，上の重回帰モデルは次のよう に書き換えることができます。 𝑦 = 𝑦̂ + 𝜀 𝑦 ̂ = 𝑏 0 + 𝑥 1 𝑏1 + 𝑥 2 𝑏2 𝜀 = 𝑥3 𝑏3 + ⋯ + 𝑥99999998 𝑏99999998 + 𝑥99999999 𝑏99999999 + 𝑥100000000 𝑏100000000 回帰分析では，どれだけ少ない変数で十分な説明ができるかが問われているわけです。 また，Multiple R-squared はその名の通り重相関係数の二乗という意味合いもあります。 重相関係数とは，𝐲 と 𝐲̂ の相関のことです。実際に相関係数を出してみましょう。なお predict() 関数は，lm() の出力を第一引数に，データを第二引数に与えると，各行の説明 変数の値と lm() で得られた回帰係数を用いて 𝐲̂ を計算してくれる関数です。 重相関係数の二乗の計算 1 # 予測値を出す 2 yhat <- predict(result_lm, dat) 3 # yhatとyの相関 4 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 Warning: The `size` argument of `element_line()` is deprecated as of 2 3.4.0. 3 i Please use the `linewidth` argument instead. 4 This warning is displayed once every 8 hours. 5 Call `lifecycle::last_lifecycle_warnings()` to see where this warning 6 was generated. ggplot2 回帰分析では，説明変数の数を増やすと必ず重相関係数は高くなります。仮に完全に無関係 な変数（乱数など）を加えたとしても，重相関係数が低下することはありません。先程の 1 17

5.2 重回帰分析 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 調整済み決定係数 𝑅𝑎𝑑𝑗 はあまり使えるものではない，という認識でも良いような気がしま す*8 。 5.2.5 標準偏回帰係数 偏回帰係数は，それぞれ「その説明変数が 1 大きくなったら 𝑦 ̂ はどの程度大きくなるか」を 表しているため，そのままでは比較することができません。age と gender の係数を比べて gender のほうが影響が大きい，と判断はできないのです。その理由の一つは，変数のスケ *8 18 良いモデルを選びたい場合には，決定係数の比較よりも交差検証 (cross validation) などを用いるのが良い とされているようです。

5.2 重回帰分析 ールが異なるためです。age が 1 増えると言うのは，年齢が一つ変わるだけでさほど大きな 変化ではありません。対して gender が 1 増えると言うのは，性別が男性から女性に変わる ことを意味するため，これはかなり大きな変化です。変数の変化の程度を調整するための方 法として標準化というものがありました。標準化したあとの変数は，すべて平均 0，分散 1 になるため，どんな変数でも「1 増える」ということの意味が「1 標準偏差増える」というこ とになります。 標準偏回帰係数とは，説明変数と結果変数のすべてを標準化した上で算出された偏回帰係数 のことです。結果変数も標準化することによって，係数は「説明変数が 1 標準偏差大きくな ったら，𝑦 が標準偏差いくつ分大きくなるか」を表すことになり，異なる変数の間でも大小 の比較が可能そうな気がしてきます。 R には標準化を行う関数として scale() というものがあるので，これを使って標準偏回帰 係数を出してみましょう。 scale() 関数 1 x <- 1:10 2 scale(x) # 標準化 1 [,1] 2 [1,] -1.4863011 3 [2,] -1.1560120 4 [3,] -0.8257228 5 [4,] -0.4954337 6 [5,] -0.1651446 7 [6,] 0.1651446 8 [7,] 0.4954337 9 [8,] 0.8257228 10 [9,] 1.1560120 11 [10,] 1.4863011 12 attr(,"scaled:center") 13 [1] 5.5 14 attr(,"scaled:scale") 15 [1] 3.02765 1 (x - mean(x)) / sd(x) # 平均値を引いてから標準偏差で割る 1 [1] -1.4863011 -1.1560120 -0.8257228 -0.4954337 -0.1651446 2 [7] 0.4954337 0.8257228 1.1560120 1.4863011 標準化してから回帰分析 19 1 # 使う変数だけ標準化 2 dat_std <- scale(dat[, c("Q1_A", "age", "gender")]) 0.1651446

5.3 一般化線形モデル

3

# matrix型になってしまうので明示的にデータフレームに変更

4

dat_std <- data.frame(dat_std)

5

# 標準化したデータで回帰

6

result_lm <- lm(Q1_A ~ age + gender, data = dat_std)

7

# あとはsummary()

8

summary(result_lm)

1
2

Call:

3

lm(formula = Q1_A ~ age + gender, data = dat_std)

4
5

Residuals:

6
7

1Q

Median

3Q

Max

-4.2178 -0.5752

Min

0.1340

0.7000

1.9468

8
9

Coefficients:

10

Estimate Std. Error t value Pr(>|t|)

11

(Intercept) 5.364e-16

1.956e-02

0.000

1

12

age

1.709e-01

1.958e-02

8.727

<2e-16 ***

13

gender

1.956e-01

1.958e-02

9.990

<2e-16 ***

14

---

15

Signif. codes:

0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

16
17

Residual standard error: 0.9647 on 2429 degrees of freedom

18

Multiple R-squared:

19

F-statistic: 91.56 on 2 and 2429 DF,

0.07011,

Adjusted R-squared:

0.06934

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 一般化線形モデル
ここまで説明してきた回帰分析モデルは，

• 説明変数の単純な線形結合で結果変数を予測し

20

5.3 一般化線形モデル • 誤差は全データで共通の正規分布𝑁 (0, 𝜎𝜀 ) に従う という仮定のもとでのモデルでした。ただ，世の中にはそんなに単純ではないケースが結構 あります。例えば結果変数が「病気になるかどうか」の場合，取りうる値は 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, 𝜎𝜀 ) に従う」という仮定に関しても問題がありそうです。図5.8で は，X 軸が 10 本くらいのところでは罹患している人としていない人がばらばらといるため， それなりの分散がありそうな一方で，0 本や 20 本のあたりでは全員が同じ値になっていま す。つまり X 軸の値によって誤差分布の分散 𝜎𝜀 が変動しているように思われるため，やは り直線を当てはめるのはおかしそうです。 ということで，回帰分析モデルを拡張する必要があります。予測値の変動を説明変数の線形 結合で表した (5.6) 式 𝑦 = 𝑏0 + 𝑥1 𝑏1 + 𝑥2 𝑏2 + 𝑥3 𝑏3 + ⋯ + 𝑥𝐾 𝑏𝐾 の右辺はわかりやすいのでこのまま使うとして，直線ではなく確率を表す，つまり取りうる 値が 0 から 1 に収まるように変換することを考えます。いくつか選択肢はあるのですが，こ こではロジスティック変換をしてみましょう。ロジスティック関数とは， 𝑔(𝑥) = 21 1 1 + exp(𝑥) (5.19)

5.3 一般化線形モデル という関数です。図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.20) となり，この右辺 𝑦 は 0 から 1 までの値しか取らなくなります。また，式を変形させると log 𝑦 = 𝑏 0 + 𝑥 1 𝑏1 + 𝑥 2 𝑏2 + 𝑥 3 𝑏3 + ⋯ + 𝑥 𝐾 𝑏𝐾 1−𝑦 (5.21) という形になり，もとのきれいな右辺に戻すことも可能です。 またこの変換では（誤差）分散の不均一性も表現しています。というのは図5.9を見るとわか るように，確率が 0.5 付近では説明変数の線形和（X 座標）が少し変わるだけで予測確率が 大きく変動するのに対して，確率が 0 または 1 付近の極端なエリアでは，ほとんど変動がな くなっています。 ということで，このように回帰分析モデルを変形させることで正規分布以外の分布にも対応 可能にしたものが一般化線形モデル (generalized linear model; GLM) です*9 。一般化線 形モデルでは，モデル式が 𝑔(𝑦) = 𝑏0 + 𝑥1 𝑏1 + 𝑥2 𝑏2 + 𝑥3 𝑏3 + ⋯ + 𝑥𝐾 𝑏𝐾 (5.22) という形になっています。つまり，「𝑦 をなんか変換したもの」𝑔(𝑦) を 𝑥 の線形和で予測し ようという試みなわけです。この「𝑦 をなんか変換したもの」が，ロジスティック回帰の場 合には 𝑔(𝑦) = log *9 22 𝑦 1−𝑦 (5.23) これに対して正規分布の普通の重回帰モデルを一般線形モデル (general linear model) と呼んだりします。

5.3 一般化線形モデル というロジット変換*10 の式になっており，またふつうの重回帰モデルも 𝑔(𝑦) = 𝑦 (5.24) という恒等関数によって変換されたものと見れば，これが一般化線形モデルの特殊なケース だということがわかるでしょう。変換の関数 𝑔(𝑦) はリンク関数と呼びます。リンク関数は 「二項分布のときにはロジット変換」というように基本的には 𝑦 の分布に対してお決まりのも のがあります。一般化線形モデルでは変形の形ごとに異なる名称がついており，ロジスティ ック変換を用いた場合をロジスティック回帰と呼びます。もう少ししっかりとロジスティッ ク回帰モデルの全体を書くと， 𝑦𝑝 = 1 +𝜀 , 1 + exp [−(𝑏0 + 𝑥𝑝1 𝑏1 + 𝑥𝑝2 𝑏2 + 𝑥𝑝3 𝑏3 + ⋯ + 𝑥𝑝𝐾 𝑏𝐾 )] 𝑝 𝜀𝑝 ∼ 𝑁 (0, 𝑦𝑝 (1−𝑦𝑝 )) となります。つまり誤差分散は確率が 0.5 付近では最も大きく，0 や 1 付近では小さくなる ことを織り込んだ上でモデルのパラメタ（回帰係数）を推定しているのです。 プロビット回帰 リンク関数を標準正規分布累積確率関数の逆関数 𝑔(𝑦) = Φ−1 (𝑦) と置いた GLM は， プロビット回帰と呼びます。ロジスティック回帰とプロビット回帰はどちらも「結果 変数が 0 か 1 の二値」というケースに対して適用されるモデルです。そして 2 つの回 帰で用いられるリンク関数の逆関数はよく似た形（左右対称の釣鐘型）をしています。 そんな 2 つの回帰モデルについて，特に経済学など？では理論的な背景なども考慮し て厳密に使い分けられる節があるようなのですが，この講義で紹介する手法（特に項 目反応理論）では，基本的にその区別には関心がありません。言ってしまえば「どっ ちでも良い」という感じがあり，基本的には計算しやすいという理由からロジスティ ック回帰のほうを使用することが多いです。 ということで，先程のタバコの例でもロジスティック回帰による回帰直線を当てはまると， 図5.10のようになります。 一般化線形モデルでは他にもポアソン回帰やガンマ回帰を始めとしていくつかありますが， ここでは全ては紹介しません。以降の回では，ロジスティック回帰だけ理解しておけばつい ていけるはずです。 5.3.1 R で一般化線形モデル R で一般化線形モデルを行う場合は，lm() の代わりに glm() という関数を使います。基本 的には使い方は同じですが，GLM なので「誤差分散の分布の形」および「リンク関数」を指 定してあげる必要があります。指定の仕方にちょっとクセがありますが，そういうものとし て理解してください。 *10 23 ロジスティック関数とロジット関数は互いに逆関数の関係になっています。そのため，説明変数側の変換と して見た場合にはロジスティック変換と呼ばれる一方で，結果変数側の変換（リンク関数）として見る場合 にはロジット変換と呼ばれるわけです（たぶん） 。実際に回帰分析としての手法名では「ロジスティック回帰」 とも「ロジット回帰」とも呼ばれることがあるはずです。

5.3 一般化線形モデル
1.00

肺がん罹患有無

0.75

0.50

0.25

0.00
0

図5.10

10
20
1日あたりタバコの平均本数

30

ロジスティック回帰直線をあてはまる

GLM（ロジスティック回帰）
1

# 結果変数は二値変数でないといけないので，「Q1_Aが平均以上か」を表す変数を作成

2

# datには追加していない

3

Q1_A_binom <- dat$Q1_A > mean(dat$Q1_A)

4

# 実は引数dataの中に無い変数でも分析に利用できる

5

result_glm <- glm(Q1_A_binom ~ age + gender, data = dat, family =
binomial(link = "logit"))

6

summary(result_glm)

1
2

Call:

3

glm(formula = Q1_A_binom ~ age + gender, family = binomial(link =
"logit"),

4

data = dat)

5
6

Deviance Residuals:

7

Min

1Q

Median

3Q

Max

8

-2.0010

-1.2168

0.8228

1.0760

1.6097

9
10

Coefficients:

11

Estimate Std. Error z value Pr(>|z|)

12

(Intercept) -1.759336

0.189153

13

age

0.025611

0.003922

6.530 6.60e-11 ***

14

gender

0.706794

0.088404

7.995 1.29e-15 ***

15

---

16

Signif. codes:

-9.301

< 2e-16 ***

0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

17
18

(Dispersion parameter for binomial family taken to be 1)

19

24

20

Null deviance: 3358.4

on 2431

degrees of freedom

21

Residual deviance: 3246.0

on 2429

degrees of freedom

22

AIC: 3252

5.3 一般化線形モデル

23
24

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 =

2

summary(result_glm)

gaussian(link = "identity"))

1
2

Call:

3

glm(formula = Q1_A ~ age + gender, family = gaussian(link = "identity"),

4

data = dat)

5
6

Deviance Residuals:

7

Min

1Q

Median

3Q

Max

8

-19.1742

-2.6150

0.6091

3.1820

8.8502

9
10

Coefficients:

11

Estimate Std. Error t value Pr(>|t|)

12

(Intercept) 18.054531

0.394846

45.726

<2e-16 ***

13

age

0.070829

0.008116

8.727

<2e-16 ***

14

gender

1.891169

0.189308

9.990

<2e-16 ***

15

---

16

Signif. codes:

0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

17
18

(Dispersion parameter for gaussian family taken to be 19.23359)

19
20

Null deviance: 50241

on 2431

degrees of freedom

21

Residual deviance: 46718

on 2429

degrees of freedom

22

AIC: 14097

23
24

25

Number of Fisher Scoring iterations: 2

5.4 変数間の関係のグラフィカルモデル 当然ですが，lm() のときと同じ結果が得られます。 5.4 変数間の関係のグラフィカルモデル 最後に，今後の授業で登場するグラフィカルモデル（あるいは有向非巡回グラフ; Directed acyclic graph [DAG]）の導入をしておきます。といっても難しい事はありません。回帰分析 の例では，gender および age によって Q1_A の値を予測してみました。これをグラフィカ ルモデルで表すと図5.11のようになります。 gender Q1_A age 図5.11 重回帰モデル このように，グラフィカルモデルでは観測可能な変数を長方形で表し，矢印によって回帰の 向きを表します。また，両方向の矢印は変数間の相関関係を表しています。重回帰分析では， 複数の説明変数が一つの結果変数を説明しようとするため，グラフィカルモデルはこのよう に複数の変数から一つの変数に向かって矢印がひかれる形になります。 グラフィカルモデルを自由に拡張して，例えば図5.12のように変数間の関係を広げていった ものが，パス解析と呼ばれるやつです。矢印一個一個を回帰分析しても良いかもしれません が，全体として変数間のモデルができているのですから，全部まとめてデータとの適合度な どを評価してあげたいところです。ということでパス解析も SEM の一種です。 SEM では，このような変数間のグラフィカルモデルをもとに分析用のコードを作成してい くことになります。といっても難しいことはないはずです。データを取る時点で変数間の理 論的な関係性の仮説は考えているはずですから… gender *** age Q1_A *** 図5.12 26 パス解析

参考文献 参考文献 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.2307/2111095 27

• https://doi.org/10.2307/2111095