Rで行う基本的な統計検定とその選択フローチャート5
今回は生物学で用いられている基本的な統計検定をRで行う際の手順と実際のコードを特集します。自身のデータセットの特徴と実験目的から適切な統計検定法を選択し、それをRで実行できるようになるための指針となれば幸いです。
第一回：平均値の差の検定（2群）
第二回：一元配置の分散分析
第三回：一元配置の多重比較
第四回：二元配置の分散分析とその後の多重比較
第五回：比率の差の検定
第六回：サンプルサイズの設計法
第七回：回帰分析
第五回 回帰分析
Rで重回帰をやるためのコードを紹介する。
今回用いるデータセットはクマノミの産卵に関する仮想データで産卵間隔（interval）、平均水温（water-temp）、平均摂餌量(feed)、平均塩濃度の差分（Δsalinity）、平均pH(pH)、メスの体長(female-SL)、オスの体長(male-SL)、産卵基質板の大きさ(plate-area)、初産からの経過期間(mpFS)、前回の産卵間隔（lastSPint）からなる。以下はデータの一部。
このデータから産卵間隔を目的変数にして重回帰分析を行って次の産卵日を予測したい。
さっそくRで実行してみる。
dat <- data.frame(read.csv("~/Desktop/interval.csv"))
# 各変数の相関関係を確認し、可視化する
# 絶対値0.7以上でどちらかの変数を削除することを検討する
# corrplotで可視化：https://tarutaru-blogs.com/r-multiple-regression-analysis/
library(corrplot)
corrplot(cor(dat), method = "color", addCoef.col = TRUE)
各変数同士の相関関係を表すヒートマップを作成した。色が濃いほど相関関係が強いことを表す。
メスとオスの体長、メスの体長とmpFS(month post First Spawning)、オスの体長とmpFSが強く相関している。年月が経てば大きくなっていくので当然だろう。これら3つはinterval daysにもあまり相関がないようなので、重回帰解析からは省略（3つのうち1つを残す）してもよいだろう。
psychパッケージを使って変数の関係を可視化することもできる
参照：https://navaclass.com/multiple-regression/
# psychで可視化
install.packages("psych")
library(psych)
psych::pairs.panels(dat)
ヒートマップは同じ情報が2回出てきていたが、psychでは各変数のヒストグラム、変数同士の相関係数と散布図も表示してくれる。
とにかくオスメスの体長データをデータフレームから削除して（mpFSを残して）、重回帰を実行してみる。
dat <- dat[, colnames(dat) != "male.SL"]
dat <- dat[, colnames(dat) != "female.SL"]
result <- lm(interval ~ ., data = dat)
summary(result)
結果
Call:
lm(formula = interval ~ ., data = dat)
Residuals:
Min 1Q Median 3Q Max
-18.9058 -2.8537 0.5875 2.4376 20.2136
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.322e+02 1.150e+02 3.758 0.000836 ***
water.temp 1.626e+00 3.016e+00 0.539 0.594164
feed -4.821e-04 8.606e-02 -0.006 0.995572
Δsalinity 3.892e-01 4.955e+00 0.079 0.937961
plate.area -1.712e-03 2.141e-02 -0.080 0.936867
DO -5.207e-01 9.867e-01 -0.528 0.602020
pH -5.752e+01 1.053e+01 -5.462 8.83e-06 ***
mpFS -1.752e-01 2.283e-01 -0.768 0.449441
lastSPint -7.195e-02 9.440e-02 -0.762 0.452571
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6.56 on 27 degrees of freedom
Multiple R-squared: 0.7844, Adjusted R-squared: 0.7206
F-statistic: 12.28 on 8 and 27 DF, p-value: 3.183e-07
pHがintervalに有意に寄与していそうとわかる。ここでは右端のp値をもとに判定しているが、これは偏回帰係数の信頼区間に0が含まれているかどうかを検定している（＝偏回帰係数が0でないかどうか）。ちなみにEstimateは偏回帰係数と呼ばれ、回帰式の係数であるが、変数同士の偏回帰係数の大小を比べても寄与率の違いは分からない（変数自体の大小に左右されるので）。どの変数がどれくらい寄与しているのかを調べるには、標準回帰係数を比較する必要がある（後述）。
分散拡大要因を調べておく。vifは10を超えないように、イメージとしてはVIFが3で「ちょっとまずい」、5で「まあまあまずい」、10で「かなりまずい」でいいらしい
参照：https://best-biostatistics.com/correlation_regression/multi-co.html
library(car)
vif(result)
結果
water.temp feed Δsalinity plate.area DO pH mpFS lastSPint
1.100880 1.803353 2.135544 2.130872 1.379255 3.279187 2.589269 1.118274
vifは大丈夫そう。
全ての変数を使った重回帰モデルが構築できたところで、次はどの変数を使ってどの変数を使わないかを決めていく。ステップワイズ法という1つずつ変数を減らしてモデルを評価する方法で行う。
# 変数を選択して最適なモデルを作成する
result.step<-step(result)
summary(result.step)
結果
Start: AIC=143.08
interval ~ water.temp + feed + Δsalinity + plate.area + DO +
pH + mpFS + lastSPint
Df Sum of Sq RSS AIC
- feed 1 0.00 1162.0 141.08
- Δsalinity 1 0.27 1162.2 141.08
- plate.area 1 0.28 1162.2 141.09
- DO 1 11.98 1173.9 141.44
- water.temp 1 12.51 1174.5 141.46
- lastSPint 1 25.00 1187.0 141.84
- mpFS 1 25.35 1187.3 141.85
<none> 1162.0 143.08
- pH 1 1283.77 2445.7 167.87
Step: AIC=141.08
interval ~ water.temp + Δsalinity + plate.area + DO + pH + mpFS +
lastSPint
Df Sum of Sq RSS AIC
- Δsalinity 1 0.28 1162.2 139.09
- plate.area 1 0.29 1162.2 139.09
- DO 1 12.39 1174.3 139.46
- water.temp 1 12.54 1174.5 139.46
- lastSPint 1 25.12 1187.1 139.85
- mpFS 1 31.06 1193.0 140.03
<none> 1162.0 141.08
- pH 1 1899.26 3061.2 173.95
Step: AIC=139.08
interval ~ water.temp + plate.area + DO + pH + mpFS + lastSPint
Df Sum of Sq RSS AIC
- plate.area 1 0.1 1162.4 137.09
- DO 1 12.2 1174.5 137.46
- water.temp 1 12.3 1174.6 137.47
- lastSPint 1 25.5 1187.7 137.87
- mpFS 1 34.6 1196.8 138.14
<none> 1162.2 139.09
- pH 1 3342.9 4505.1 185.86
Step: AIC=137.09
interval ~ water.temp + DO + pH + mpFS + lastSPint
Df Sum of Sq RSS AIC
- DO 1 12.2 1174.6 135.47
- water.temp 1 12.6 1174.9 135.48
- lastSPint 1 25.4 1187.7 135.87
- mpFS 1 50.8 1213.1 136.63
<none> 1162.4 137.09
- pH 1 3473.3 4635.7 184.89
Step: AIC=135.47
interval ~ water.temp + pH + mpFS + lastSPint
Df Sum of Sq RSS AIC
- water.temp 1 7.2 1181.7 133.68
- lastSPint 1 26.2 1200.8 134.26
<none> 1174.6 135.47
- mpFS 1 74.3 1248.9 135.67
- pH 1 4174.9 5349.5 188.04
Step: AIC=133.68
interval ~ pH + mpFS + lastSPint
Df Sum of Sq RSS AIC
- lastSPint 1 26.1 1207.8 132.47
<none> 1181.7 133.68
- mpFS 1 75.9 1257.7 133.93
- pH 1 4167.9 5349.6 186.05
Step: AIC=132.47
interval ~ pH + mpFS
Df Sum of Sq RSS AIC
<none> 1207.8 132.47
- mpFS 1 70.9 1278.7 132.52
- pH 1 4170.6 5378.4 184.24
> summary(result.step)
Call:
lm(formula = interval ~ pH + mpFS, data = dat)
Residuals:
Min 1Q Median 3Q Max
-18.791 -1.978 0.399 1.901 21.239
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 468.9783 41.7925 11.222 8.37e-13 ***
pH -57.4220 5.3793 -10.675 3.08e-12 ***
mpFS -0.1826 0.1312 -1.391 0.173
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6.05 on 33 degrees of freedom
Multiple R-squared: 0.7759, Adjusted R-squared: 0.7623
F-statistic: 57.14 on 2 and 33 DF, p-value: 1.912e-11
何回も変数の増減を行ってモデルを評価している。モデルの評価基準はAIC(Akaike's Information Criterion）を使って、AICが小さいほど良いモデルと判定する。
最終的に採択されたモデルを> summary(result.step)で表示している。
この場合は、intervalを最もよく説明するモデルはpHとmpFSを使っている。
偏回帰係数は両方負であるので、
pHが低いほど産卵間隔は長くなり、初産から期間が経っているほど産卵間隔は短くなることがわかる。
とにかく以下のような回帰式が得られた。
interval = 468.9783 - 57.422 x pH - 0.1826 x mpFS
前回の産卵日から今日までの10日間の平均pHは7.8、初産から28ヶ月目なので、
今回の産卵間隔は15.9739daysなので16日間隔、つまり6日後です！！
これで産卵日が予測しやすくなって、実験が捗りますねぇ（そううまくはいかない）
では、pHと初産からの期間のどちらがより産卵間隔に寄与するのか、次は標準回帰係数を調べていく。
偏回帰係数は独立変数がどのような値をとるかに依存するので、そのままでは数値の大小を比較して影響度を論じることはできない。そこで標準回帰係数を計算する。標準回帰係数bは偏回帰係数B、独立変数の標準偏差Si、従属変数の標準偏差Syを用いて、b＝B・Si/Syで表される。しかしいちいち計算するのは面倒なので先に全ての独立変数を標準化してから回帰分析を行う。
# 標準回帰係数を求めて各変数の重要度を定量化する
zdat <- data.frame(scale(dat))
zresult <- lm(interval ~ ., data = zdat)
zresult.step<-step(zresult)
summary(zresult.step)
結果
Call:
lm(formula = interval ~ pH + mpFS, data = zdat)
Residuals:
Min 1Q Median 3Q Max
-1.51417 -0.15939 0.03215 0.15319 1.71148
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.456e-15 8.125e-02 0.000 1.000
pH -8.822e-01 8.265e-02 -10.675 3.08e-12 ***
mpFS -1.150e-01 8.265e-02 -1.391 0.173
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4875 on 33 degrees of freedom
Multiple R-squared: 0.7759, Adjusted R-squared: 0.7623
F-statistic: 57.14 on 2 and 33 DF, p-value: 1.912e-11
データフレームを標準化し、先ほど同様に重回帰、ステップワイズ法で変数選択を行った。
summaryで表示されているEstimateが標準回帰係数になる。これをみるとpHのほうが標準回帰係数の絶対値が大きいので、intervalにより寄与していると考えられる。
以下では標準回帰係数の大小を可視化する。
#重回帰分析の各説明変数の係数及び切片を抽出し、可視化する
install.packages("coefplot",dependency=TRUE)
library("coefplot")
coef(summary(zresult))
coefplot(zresult)
これらからクマノミの産卵日の間隔はpHが低いと長くなり、初産から期間が経てば短くなるという傾向が得られた（もちろん仮想データに基づいているので実際どうかは知りません）。
参照：
http://mizumot.com/handbook/?page_id=484
https://www.gixo.jp/blog/3186/
https://www1.doshisha.ac.jp/~mjin/R/Chap_15/15.html
https://navaclass.com/multiple-regression/
https://tarutaru-blogs.com/r-multiple-regression-analysis/
https://istat.co.jp/ta_commentary/multiple_02
https://www.gixo.jp/blog/3186/