Multiple Regression

1何のための手法か

「アクセントが変化した数」を、集落の何が説明するか。

小豆島の9集落について、中井[1998]の記述と2012年の調査とを突き合わせ、2拍名詞のアクセント型がいくつ変わったかを数えた値がある。高拍 H・低拍 L・中拍 M・拍内上昇 R・拍内下降 F をそれぞれ 11・00・0.5 0.5・01・10 という数値に置き換え、1998年側と2012年側で対応するセルの差の絶対値を取り、地点ごとに和をとったものである（原典 pp.49–50、表4-6〜表4-8）。この和が変化数である。

問いは「この変化数の集落ごとの違いは、集落の何と結びついているか」である。原典が説明変数として並べたのは、人口（平成22年・平成12年の国勢調査小地域集計）、世帯数（同）、人口増加率（人口H22÷人口H12）、3種類の中心性（情報中心性・媒介中心性・次数中心性）、標高（基盤地図情報の5mメッシュ数値標高モデル）である（原典 pp.47–48、p.52）。中心性は、バスルートと定期旅客航路のデータから集落間のつながりを 0/1 で表した隣接行列（表4-5）を作り、R の sna パッケージの degree・betweenness・infocent で算出している（原典 pp.50–52）。

ここで重回帰分析が必要になるのは、これらの説明変数が互いに強く相関しているからである。人口が多い集落は世帯数も多く、道路が集中している集落は媒介中心性も次数中心性も高い。単相関を9本並べても、どれがどれの言い換えなのかが分からない。重回帰は「他の変数を固定したときの効き」という形で答えを出そうとするが、その「固定」が効かないほど変数同士が似ているとき、係数は不安定になる。第3節以降はこの問題への対処の話である。

2モデル式と最小二乗

残差二乗和を最小にする係数は、正規方程式の解として一意に書ける。

観測が \(n\) 個、説明変数が \(p\) 個あるとする。\(i\) 番目の観測について

計画行列 \(X\)（第1列が全部 1、以下が各説明変数）と係数ベクトル \(\beta\) を使って \(y = X\beta + \varepsilon\) と書く。最小二乗法は、残差二乗和

3説明変数が互いに相関しているという問題

多重共線性は「推定できない」問題ではなく、「推定値の分散が膨らむ」問題である。

原典の図4-1（表4-9 を標準化したデータの相関係数行列・散布図行列、p.54）では、人口H22・世帯数H22・人口H12・世帯数H12 の4つが互いに 0.998〜0.999 という値を示している。中心性どうしの相関も高く、原典 p.53 は「媒介中心性、次数中心性 同士の相関についても高い」、2拍名詞（＋助詞）の側については p.62 で「情報中心性、次数中心性 同士の相関も高い」と述べている。これだけ似た変数を同時にモデルに入れると何が起きるかを、係数の分散の形で見る。

\(j\) 番目の説明変数 \(x_j\) を、残りのすべての説明変数（定数項を含む）で回帰し、その残差を \(\tilde{x}_j\) と書く。Frisch–Waugh–Lovell の定理により、重回帰の係数 \(\hat{\beta}_j\) は、\(y\) を \(\tilde{x}_j\) だけに回帰した単回帰係数と一致する。

\(x_j\) を他の説明変数で回帰したときの決定係数を \(R_j^2\) とすると、定義から \(\mathrm{RSS}_j = \mathrm{TSS}_j\,(1-R_j^2)\)、ただし \(\mathrm{TSS}_j = \sum_i (x_{ij}-\bar{x}_j)^2\) である。これを代入すると、

4All possible subset selection

\(2^p\) 通りを全部当てはめて、情報量規準の小さい順に並べる。

説明変数が \(p\) 個あるとき、それぞれを「使う／使わない」で分けた組み合わせは \(2^p\) 通りある。空集合を除けば \(2^p - 1\) 通りである。All possible subset selection は、この全部について lm を実行し、決定係数・自由度調整済み決定係数・AIC・BIC などを計算して並べ替える、という手続きである。原典が使ったのは矢野環による R スクリプト papca.R 中の関数 All.possible.subset.selection で、全文が原典 Appendix. 4（A-11〜A-12）に収録されている。

並べ替えの基準は AIC である。対数尤度 \(\log L\) とパラメータ数 \(k\) について

5実際のデータと、変数を絞っていく過程

9集落 × 2行、説明変数9個。ここから2変数まで絞られる。

原典 p.52 の表4-9（2拍名詞のみ）をそのまま示す。18行あるが、集落は9つである。原典の脚注13 は「同一集落内での調査地点数が2地点もしくは音声の揺れにより同一地点に2パターンの音声が認められる場合、重回帰分析をする際に該当の地点の重みが大きくなる。これを防ぐため、調査地点が1地点および音声のパターンが1つである調査地点のデータ列を2つにした」と説明している。説明変数の値は、同じ集落の a と b のあいだで完全に一致する。

原典 表4-9「重回帰分析のデータ（2拍名詞のみ）」（p.52）。原典 脚注12 により、[変化数] は表4-8 の「和」に相当し、[人口増加率] は 人口H22÷人口H12 である。

原典が踏んだ手順は次のとおりである（原典 pp.53–58、Appendix. 5）。

1. 標準化　表4-9 を scale で標準化する（data1 → 標準化後）。
2. 相関の確認　相関係数行列と散布図行列（図4-1）を描き、人口H22・世帯数H22・人口H12・世帯数H12 の相関がいずれも 0.9 以上であることを確認する。
3. 1回目の総当たり　説明変数9個について All.possible.subset.selection を実行する（表4-10）。人口系4変数のうち、AIC の低い組み合わせに現れやすい世帯数H22 のみを残し、人口H22・人口H12・世帯数H12 を消去して data1.2 を作る。
4. 2回目の総当たり　data1.2（説明変数6個）について再度実行する（表4-11）。AIC 最小は 世帯数H22 + 情報中心性 + 媒介中心性 + 次数中心性 + 標高。
5. VIF と感度分析　このモデルについて vif と sensitivity を実行する。
6. 3回目の総当たり　VIF が 10 を超えた 世帯数H22・媒介中心性 を消去して data1.3 を作り、再度実行する（表4-14）。AIC 最小は 情報中心性 + 次数中心性。

原典 表4-11「Data1.2 についてのモデル選択（2拍名詞のみ）」（p.55）の上位5件のみを抜粋。原典は20件を掲載している。AIC 最小の行と3番目の行では AIC 差が約 6.1 ある一方、cAIC では 36.912 対 34.900 と順位が逆転している点は、原典の表からそのまま読み取れる。

6結果

残ったのは、情報中心性（負）と次数中心性（正）の2つである。

第2段階（data1.2、2拍名詞のみ）。AIC で選ばれたモデルの回帰式は、原典 p.54 に次の形で書かれている。

推定値・標準誤差・\(t\) 値・\(p\) 値は原典 表4-12「data1.2 の重回帰分析結果」（p.56）。VIF は原典 p.55 の本文および Appendix. 5（A-17）の vif(RegModel.1) 出力。原典は「VIF が 10 になると多重共線性が起きているとする指標があり、この指標に則れば、次数中心性、情報中心性、標高が適切な変数である」と述べている。

第3段階（data1.3、2拍名詞のみ）。VIF が 10 を超えた 世帯数H22 と 媒介中心性 を落として総当たりをやり直すと、AIC 最小は 情報中心性 + 次数中心性 の2変数モデルになる。

原典 表4-15「data1.3 についての重回帰分析の結果」（p.57）、VIF は原典 p.58 および Appendix. 5（A-19）。決定係数 0.767、自由度修正済み決定係数 0.736、\(F\) 値 24.7（\(p\) 値 \(1.79\times10^{-5}\)、自由度 2 と 15）、残差の標準誤差 0.514（原典 p.58）。
なお原典 p.58 の本文中の回帰式は「変化数 = 0.000 + 1.901 次数中心性 − 1.595 情報中心性」となっており、表4-15 および Appendix. 5 の summary 出力（1.90 / −1.59、\(t\) = 6.99 / −5.85）と末尾の桁が一致しない。ここでは表4-15 の値を載せた。

2拍名詞（＋助詞）のケース。同じ9集落・同じ説明変数で、目的変数だけを「助詞を後続させた形も含めた変化数」（原典 表4-18）に替えた分析も行われている。第2段階（data2.2）で AIC 最小となったのは 世帯数H22 + 人口増加率 + 情報中心性 + 媒介中心性 + 次数中心性 + 標高 で、決定係数 0.896、自由度修正済み決定係数 0.839、\(F\) 値 15.77（\(p\) 値 \(8.00\times10^{-5}\)）、残差の標準誤差 0.401（原典 p.63）。

原典 表4-21「data2.2 の重回帰分析結果」（p.64）、VIF は原典 p.63 本文。\(p\) 値が 0.05 を下回らないのは人口増加率のみ。

ここでも VIF が 10 を超えた 世帯数H22 と 媒介中心性 を落として総当たりをやり直すと（data2.3）、AIC 最小はやはり 情報中心性 + 次数中心性 になる。決定係数 0.727、自由度修正済み決定係数 0.69、\(F\) 値 20（\(p\) 値 \(5.94\times10^{-5}\)）、残差の標準誤差 0.556、VIF は両変数とも 4.738（原典 pp.66–67）。

原典 表4-24「data2.3 についての重回帰分析の結果」（p.67）。

7感度分析（Sensitivity Analysis）

説明変数に揺さぶりをかけ、係数がどれだけ動くかを見る。

原典は、VIF とは別の多重共線性の診断として Sensitivity Analysis を併用している。使ったのは R の accuracy パッケージ中の sensitivity 関数で、原典 p.53 に「説明変数に変動を与えて、多重共線性を判定する Sensitivity Analysis を行うパッケージ accuracy 中の関数 sensitivity を用いた」とある。Appendix. 5 の実際の呼び出しは data1.2.sens <- sensitivity(data1.2, lm, 変化数 ~ ., ptb.R = 2000) で、出力の見出しは "Sensitivity of coefficients over 2000 perturbations:" である。原典 p.55 の本文も「2000回の変動を与えて、チェックを行った」と書いている。

原典 表4-13「data1.2 に対する Sensitivity Analysis の結果（1）」（p.56）。「不安定」欄はいずれも空欄で、原典 p.57 は「どの変数も不安定（Unstable）でないという結果が得られた」と述べている。なお、同じ出力が Appendix. 5（A-18）にも digits = 2 で印字されており、2.5%／97.5% の値は末尾1桁が表4-13 と食い違う箇所がある（例：世帯数H22 の 2.5% が表4-13 では −1.100、Appendix. 5 では −1.099）。data2.2・data1.3・data2.3 についても同様の出力が原典 表4-22・表4-16・表4-25 に掲載されており、いずれも「不安定」の判定は出ていない。

8限界

総当たりで選んでから検定する、という手順そのものに問題がある。

9出典

このページの数値はすべて、小野原彩香『数理的アプローチからの言語変化と外言語的要素との関わりに関する研究』（同志社大学大学院文化情報学研究科、2013年度 博士論文）の第4章および巻末 Appendix から引用したものである。本文ページ番号は論文本体の通し番号、Appendix のページ番号は「A-」を冠した巻末資料の番号による。

• 小野原彩香（2013）『数理的アプローチからの言語変化と外言語的要素との関わりに関する研究』同志社大学大学院文化情報学研究科 博士論文。
  　第4章第4節「方法 ―中心性、重回帰分析―」pp.49–52 ／ 第5節「結果」pp.53–70 ／ 第6節「考察」pp.71–73。
  　表4-9（分析データ、p.52）、表4-10・表4-11（モデル選択、p.55）、表4-12・表4-13（係数と感度分析、p.56）、表4-14・表4-15（同、p.57）、表4-18（＋助詞のデータ、p.62）、表4-21（係数、p.64）、表4-22（感度分析、p.65）、表4-24・表4-25（同、p.67）。
  　Appendix. 4「事例2 All.possible.subset.selection のスクリプト」A-11〜A-12 ／ Appendix. 5「事例2 重回帰分析 R スクリプト（2拍名詞）」A-13〜 ／ Appendix. 6「事例2 重回帰分析 R スクリプト（2拍名詞＋助詞）」。