たまに、２つの相関係数が有意に異なるのかや、１つの重回帰モデル中の２つの回帰係数が有意に異なるかを示せると、主張が通りやすいという場面がある。
まぁ、あまり必要になることはないのだが、相関係数の差の検定や回帰係数の差の検定について、日本語でググるとRでパッと使える方法がまとまっているわけではなさそうだったので、ここに書いておこう。
相関係数の差の検定
相関係数の差の検定は、Rの{psych}パッケージのr.test()関数で簡単にできる。
下記に関数の解説があるが、
https://www.rdocumentation.org/packages/psych/versions/2.0.12/topics/r.testr.test function - RDocumentation
Test4のtable1が崩れていたり、「Case A: where Masculinity at time 1 (M1) correlates with Verbal Ability .5 (r12), femininity at time 1 (F1) correlates with Verbal ability r13 =.4」とあるところはM1とF1が逆になっちゃってたり、出力の説明が不足していたりするので要注意。そこから参照されている論文Steiger(1980)をみると、いくつかの疑問は解決する。
1) この関数は1つの相関係数が有意かどうか検定したい場合にも使え、たとえば変数1・2ともにサンプルサイズが100で相関係数が0.66だとすると、
library(psych) r.test(n=100, r12=0.66)
というふうにする。多くの場合両側で検定すると思うが、片側で検定したい場合は
twotailed = F
という引数を付ければよい。
2) ここから2つの相関係数の差の話に移る。まずは、2つの独立した（サンプルを共有していない）相関係数の差を検定する場合。サンプルサイズは異なっても良い。
変数1と2はサンプルサイズが100で相関係数が0.3、変数3と4はサンプルサイズが120で相関係数が0.25だった場合、
r.test(n=100, r12=0.30, r34=0.25, n2=120)
サンプルサイズが等しいのであれば、n2は付けなくてもよい（デフォルトでnと同じになる）。
3) サンプルを共有している変数が3つ（1・2・3）あって、1と2、1と3というふうに1つの変数を共有する形で相関係数同士の差を検定したい場合。
r.test(n, r12, r23, r13)
というふうにする。この場合、r12とr13を比べた検定結果が返される。このとき、r23も引数として与える点に注意。（使い方の説明が不親切だが、冒頭の解説リンクのtest4のCaseAと同じことで、そちらについては論文Steiger(1980)と照合すると分かる。）
4) サンプルは共有してるけど変数は共有していないような、2組の相関係数を比べる場合。（冒頭の解説リンクでいうとtest4のCaseBに相当）
r.test(n, r12, r34, r23, r13, r14, r24)
冒頭の解説リンクは使い方の説明が不親切なのだが、論文Steiger(1980)と照合すると分かる。差を検定したい関心のある相関係数がr12とr34で、この2つの相関係数の差の検定結果を返している。
他に参考になるページ
http://www.snap-tck.com/room04/c01/stat/stat05/stat0501.html統計学入門−第5章
母相関係数の差の検定の概要と結果 :: 【公式】株式会社アイスタット｜統計分析研究所
回帰係数の差の検定
重回帰分析の回帰係数の差の検定については、たとえば複数の回帰係数の信頼区間を出して、大きい方の下限値と小さい方の上限値が重なるかどうかという基準で検定すると、（たとえば5%水準の検定をしたいときに95%信頼区間を用いると）厳しすぎる検定になるらしい。*1
で、使えるケースが多少限られるものの、心理統計の分野で分散分析の延長で共分散分析を習うときに出てくる「平行性の検定」の考え方で、「交互作用が有意になるかどうか」という観点で検定すると簡単にできる。*2
d <- airquality # 練習用データセットの読み込み d <- d[complete.cases(d),] # 欠損値削除 summary(lm(Ozone~Wind+Temp, data = d))
このようにすると以下のような結果が得られる。
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -67.3220 23.6210 -2.850 0.00524 ** Wind -3.2948 0.6711 -4.909 0.0000032617283 *** Temp 1.8276 0.2506 7.294 0.0000000000529 ***
この-3.29と1.83が有意に異なるかを知りたいという話で、まあ標準誤差とかをみても余裕で異なりそうではあるが、以下のように（従属変数であるOzoneは残して）ロング型にデータを変更し、交互作用の検定をすればいい。これは要するに、2つの変数を1つの変数にまとめて、変数名を新たな変数（ダミー変数）として設け、いわば「値」と「変数名」の交互作用をみるような感じ。
d2 <- d %>% select(Ozone, Wind, Temp) %>% # 変数を限定 pivot_longer(-Ozone, names_to = 'Variable', values_to = 'Value') # ロング型に変換 summary(lm(Ozone~Variable*Value, data = d2))
以下のような結果が得られる。
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -147.6461 19.7613 -7.471 0.0000000000019 *** VariableWind 246.6873 21.0072 11.743 < 0.0000000000000002 *** Value 2.4391 0.2522 9.673 < 0.0000000000000002 *** VariableWind:Value -8.1679 0.7210 -11.329 < 0.0000000000000002 ***
「VariableWind:Value」のところが交互作用（Tempが基準=0になって要するにWindダミーになってるという意味）で、有意になってるので、さっきの回帰係数の差は有意ってことになります。
後で知ったのですが、回帰係数の差の検定については、{car}パッケージのlinearHypothesis関数を使って、簡単に調べることができる。これは、2つの説明変数の回帰係数が等しいという帰無仮説を検定するもの。
あまりよく考えてないんですが、従属変数を複数もつ多変量回帰（Multivariate Regression。carパッケージのManovaで検定したりする）の場合、「すべてのモデルにおいて両変数の回帰係数に差がない」という帰無仮説を検定してるような気がする。モデルの数が増えると、帰無仮説がほとんど棄却されるようになります。上のように交互作用項を入れた上でManovaする方法の場合、交互作用はそれほど有意になりやすくはないです。
# 必要なパッケージの読み込み library(car) # airqualityデータセットを使用 data("airquality") # NA値の除去 airquality_complete <- na.omit(airquality) cor(airquality_complete) # 線形回帰モデルの作成 model <- lm(Ozone ~ Solar.R + Wind + Temp, data = airquality_complete) # Solar.R と Wind の回帰係数が等しいという仮説を検定 linearHypothesis(model, "Solar.R = Wind") # Multivariate Regression（多変量回帰：従属変数が複数ある）でもできる model_multi <- lm(cbind(Ozone, Wind) ~ Solar.R + Temp + Month + Day, data = airquality_complete) # Manovaで説明変数の総合的な有意性を確認 Manova(model_multi) # Solar.R と Tempの回帰係数が等しいという仮説を検定（棄却できる） linearHypothesis(model_multi, "Solar.R = Temp") # Month と Dayの回帰係数が等しいという仮説を検定（棄却できない） linearHypothesis(model_multi, "Month = Day")
