はじめに
勉強がてら、これまでに触れたことの無かったロジスティック回帰を試してみたので、その際の記録を記事としてまとめてみました。
また、以前に散布図の書き方を記事にしたので、これを活用して作図までしてみました。
解析に使うデータ
- 以下のデータは、ある植物の生育地点における落葉層の厚さ（a0_layer）と開花の有無（flowering）を示したものです。
| 株1 | 株2 | 株3 | 株4 | 株5 | 株6 | 株7 | |
|---|---|---|---|---|---|---|---|
| 落葉層の厚さ (cm) | 12 | 10 | 14 | 13 | 11 | 20 | 17 |
| 開花の有無 | 1 | 1 | 1 | 0 | 1 | 0 | 0 |
データの読み込み
- ここでは外部ファイルから読み込まず、ベタ書きでデータフレームを作成しています。
- 外部ファイルから読み込んでデータフレームに入れる場合は、こちらの記事をご覧ください。
データの読み込み
df <- data.frame(
a0_layer = c(12, 10, 14, 13, 11, 20, 17),
flowering = c(1, 1, 1, 0, 1, 0, 0))
ロジスティック回帰の実行
-
glm関数を使って解析をする際には、第1引数で目的変数~説明変数1+説明変数2...の形でデータを指定します。- ここでは説明変数は1つだけなので、
目的変数~説明変数の形でflowering~a0_layerと指定しています。
- ここでは説明変数は1つだけなので、
- さらにsummary関数で結果を表示します。
ロジスティック回帰の実行
ans <- glm(flowering~a0_layer, data=df, family=binomial(link="logit"))
summary(ans)
-
(Intercept)行は切片を表しています。 -
Estimate列は偏回帰係数、Std. Errorは偏回帰係数の標準誤差、Pr(>|z|)は偏回帰係数に対する有意性検定のp値を表しています。- 以下の数式だと、切片が
b0、偏回帰係数はb1やb2となります。
- 以下の数式だと、切片が
p = \frac{1}{1+\exp\{-(b_{0}+b_{1}x_{1}+b_{2}x_{2}+...\}}
- ここではp値が0.2となっているため、「説明変数（落葉層の厚さ）の目的変数（開花の有無）に対する影響は統計的に有意ではない」となります。
summary関数の結果
Call:
glm(formula = flowering ~ a0_layer, family = binomial(link = "logit"),
data = df)
Deviance Residuals:
1 2 3 4 5 6 7
0.54846 0.21522 1.22754 -1.55323 0.34595 -0.07311 -0.30866
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 13.4308 10.5179 1.277 0.202
a0_layer -0.9677 0.7910 -1.223 0.221
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 9.5607 on 6 degrees of freedom
Residual deviance: 4.4868 on 5 degrees of freedom
AIC: 8.4868
Number of Fisher Scoring iterations: 6
グラフの作図
-
plot関数を使って散布図を作成した後、ロジスティック曲線をcurve関数で描画しています。-
curve関数は1変数関数の描画に使える関数ですが、ここでは説明変数が1つなのでcurve関数を使っています。
-
-
ylimでY軸の最大値と最小値を指定して、yaxpでY軸の目盛りの最初と最後の値に加えて区間数を指定しています。 -
pchで開花の有無（flowering）の値に応じて、点の形状を「白抜きの丸」と「塗りつぶしの丸」に変えています。
グラフの作図
plot(
df$a0_layer, df$flowering,
xlab="Thickness of A0 layer (cm)", ylab="Flowering",
ylim=c(0,1), yaxp=c(0,1,1),
pch=ifelse(df$flowering==1, 16, 1),
cex=1.5
)
curve(1 / (1 + exp(-(13.4308 + (-0.9677*x)))), add=TRUE)
