【R】ロジスティック回帰を行う

# はじめに

勉強がてら、これまでに触れたことの無かったロジスティック回帰を試してみたので、その際の記録を記事としてまとめてみました。
また、以前に[散布図の書き方](https://qiita.com/nkojima/items/c008adec58c7834edfc1)を記事にしたので、これを活用して作図までしてみました。

# 解析に使うデータ

* 以下のデータは、ある植物の生育地点における落葉層の厚さ（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|

# データの読み込み

* ここでは外部ファイルから読み込まず、ベタ書きでデータフレームを作成しています。
  * 外部ファイルから読み込んでデータフレームに入れる場合は、[こちらの記事](https://qiita.com/nkojima/items/131fdce3a22ccc1760b4)をご覧ください。

```r:データの読み込み
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`と指定しています。
* さらにsummary関数で結果を表示します。

```r:ロジスティック回帰の実行
ans <- glm(flowering~a0_layer, data=df, family=binomial(link="logit"))
summary(ans)
```

* `(Intercept)`行は切片を表しています。
* `Estimate`列は偏回帰係数、`Std. Error`は偏回帰係数の標準誤差、`Pr(>|z|)`は偏回帰係数に対する有意性検定のp値を表しています。
  * 以下の数式だと、切片が`b0`、偏回帰係数は`b1`や`b2`となります。

```math
p = \frac{1}{1+\exp\{-(b_{0}+b_{1}x_{1}+b_{2}x_{2}+...\}}
```

* ここではp値が0.2となっているため、「説明変数（落葉層の厚さ）の目的変数（開花の有無）に対する影響は統計的に有意ではない」となります。



```text: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）の値に応じて、点の形状を「白抜きの丸」と「塗りつぶしの丸」に変えています。

```r:グラフの作図
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)
```


![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/312166/48e1b717-b7fc-d2c0-f633-fa8edf7f5fdd.png)


# 参考資料

* 青木繁伸 (2009) Rによる統計解析. オーム社.
* [【R】データの読み込み](https://qiita.com/nkojima/items/131fdce3a22ccc1760b4)
* [【R】散布図を描く方法](https://qiita.com/nkojima/items/c008adec58c7834edfc1)
* [ロジスティック回帰](https://oku.edu.mie-u.ac.jp/~okumura/stat/logistic.html)
* [plot の要点 - R システムの plot を用いた種々の散布図と折れ線グラフ](https://www.kkaneko.jp/pro/r/rplot.html)
* [軸の目盛りを変える - Rでグラフに数値を表示したり、軸の目盛間隔、軸ラベルを変更する方法](https://toukeier.hatenablog.com/entry/%3Fp%3D826#%E8%BB%B8%E3%81%AE%E7%9B%AE%E7%9B%9B%E3%82%8A%E3%82%92%E5%A4%89%E3%81%88%E3%82%8B)
* [[R]plot関数のpch一覧](https://zenn.dev/rinc5/articles/c3d6d64af4c46f)
* [統計解析ソフトRを用いたデータ解析-回帰分析偏](https://waidai-csc.jp/updata/2023/08/Rseminar-03.pdf)
* [一変数関数のグラフを描く - Rのグラフィックス](https://www.f.waseda.jp/sakas/R/Rgraphics17.html#fgraph)
* [Qiitaでの様々な数式の書き方 ～ TeX記法を使ったサンプルコード付き](https://qiita.com/ttabata/items/92a4a72b042e6a2dbbad)