図1 投与量と死亡率
<< 目次へ
最終更新日 ： 2016.9.17
2値データの解析
2値データとは性別（男性か女性か）、生死（生か死か）、回答（YesかNoか）、有無（有るか無いか）といった2つのカテゴリカルな値を持つデータのことである。だからこれは、2つの水準をもったカテゴリカル型の変数であるともいえる。しかし、2値データはカテゴリカル型の変数であるとはいえど、少し特殊な場合であり、本来、カテゴリカルなデータとは名義尺度であるが、2値データは間隔尺度としての尺度水準をもつ。そのため2値データは、例えば、平均値や相関係数を求めることもできる。
さて、今までは従属変数が連続型である場合（分散分析モデルや重回帰モデルがあった）とカテゴリカル型である、いわゆる計数データ（対数線形モデルがあった）である場合の解析方法を述べてきた。ここでは従属変数が2値データである場合の解析方法、すなわちロジスティック回帰分析の方法について解説することにする。
ロジスティック回帰分析の基本として、1つの従属変数（もちろんこれは2値データである）と1つの独立変数（ここでは連続型のデータを扱う）に対して、ロジスティック曲線を当てはめる方法を紹介する。これは殺虫剤の効果を調べるための致死量（LD50）を算出するために用いられる典型的は方法である。
例えば、マウスにある薬物を投与して投与量ごとにどれだけのマウスが死亡したかという仮想データが表1のようにあったとしよう。ここではどのような薬物を投与したかということと、投与量の単位は考えないことにしよう。
| 投与量 | 死亡数 | 生存数 | マウスの数 | 死亡率 |
| 5.42 | 3 | 53 | 56 | 0.05 |
| 5.61 | 10 | 47 | 57 | 0.18 |
| 5.78 | 15 | 44 | 59 | 0.25 |
| 5.95 | 25 | 28 | 53 | 0.47 |
| 6.12 | 49 | 11 | 60 | 0.82 |
| 6.28 | 50 | 6 | 56 | 0.89 |
| 6.43 | 58 | 1 | 59 | 0.98 |
| 6.58 | 57 | 0 | 57 | 1.00 |
表1のデータより、横軸（x軸）に投与量を、そして縦軸（y軸）に死亡率をとってプロットすると図1のようになる。いうまでもないが、死亡率は0から1までの範囲の値をとり、0というのは1匹も死亡しなかったことを意味し、1というのは全てのマウスが死亡したことを意味している。また、死亡率が0.5というのは「その投与量で半分のマウスが死亡した」ということから、LD50と呼ばれている。このLD50というのが、いわゆる致死量のことである。似たような用語でED50というものもあるが、これは薬の効果量のことであり、例えば、風邪薬において、どのくらいの容量だと熱を下げる効果が得られるかというようなものである。
図1 投与量と死亡率
図1にロジスティック曲線を当てはめることを考えよう。ロジスティック曲線とは式1のような関数が描く曲線のことである。つまり、αとβの2つのパラメータを推定することになる。パラメータを推定する方法はいくつかあるが、その方法をここで解説することはしない。ここではRのglm()を用いて2つのパラメータを推定する方法を紹介しよう。
・・・ 式1
Rのglm()の使い方を説明する前に、少しだけパラメータ推定のためのヒントを提示しておこう。
ここで 
とすると、
結局、観測値yiを log((1 / y) - 1) と変換し、それについて線形回帰分析を行えば良いことになる。ただし、分子をγとして、このγから推定するとなるともう少し話は難しくなる。なぜ分子を1と置く場合とγとして新たなパラメータとして推定する場合とがあるのかというと、現実には上限値が1でない場合もあるからである。
今回の例のように、理論的にy軸の値の上限値が1.0と決まっている場合意外は分子をγとおいて推定する必要がある。例えば、携帯電話の契約件数のように何万という単位でy軸の値（契約件数）が変化する場合などが当てはまる。
それでは実際にRのglm()を用いて解析してみよう。
# 投与量Dose, 死亡率DeathRateという変数を用意する
> Dose
[1] 5.42 5.61 5.78 5.95 6.12 6.28 6.43 6.58
> DeathRate
[1] 0.05 0.18 0.25 0.47 0.82 0.89 0.98 1.00
# 投与量と死亡率のプロット
# 引数のcexはデータポイントの大きさを指定するためのもの
> plot(Dose, DeathRate, cex=1.5)
# 続いて生存数Surv, 死亡数Deathという変数を用意する
> Surv
[1] 53 47 44 28 11 6 1 0
> Death
[1] 3 10 15 25 49 50 58 57
# glm()を使って解析するためには次のような行列データを用意する
# 1列目に死亡数, 2列目に生存数を指定する
# 関数cbind()は2つのベクトルを1つの行列にまとめる関数である
> Y <- cbind(Death, Surv)
> Y
Death Surv
[1,] 3 53
[2,] 10 47
[3,] 15 44
[4,] 25 28
[5,] 49 11
[6,] 50 6
[7,] 58 1
[8,] 57 0
# glmを用いてモデル解析を行う
# 引数familyにbinomialを指定することを忘れずに
> result <- glm(Y ~ Dose, family = binomial)
> summary(result)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -37.8887 3.3053 -11.46 <2e-16 ***
Dose 6.4031 0.5561 11.52 <2e-16 ***
---
# predict()の引数にtype="response"と指定すると,
# ロジスティック回帰を行った予測値が返される
> (x <- predict(result, type="response"))
1 2 3 4 5 6
0.03977859 0.12268535 0.29344456 0.55226264 0.78555662 0.91075150
7 8
0.96385189 0.98584987
# 予測値のデータポイントを描く
> points(Dose, x, pch=3, col="red", cex=1.5)
# ロジスティック曲線を描く
> x.jiku <- seq(5.4, 6.6, 0.01) # 5.4から6.6まで0.01きざみでデータを生成
> y.jiku <- 1 / (1 + exp(37.889 - 6.403*x.jiku)) # 曲線を描くy軸の値
> lines(x.jiku, y.jiku, col="blue", lwd=2) # lines()で曲線を描く
# このような関数を用意して予測値を求めることもできる
> f <- function(x) 1 / (1 + exp(37.889 - 6.403*x))
> f(Dose)
[1] 0.03974404 0.12258558 0.29324824 0.55202354 0.78539034 0.91066974
[7] 0.96381623 0.98583535
図2 ロジスティック曲線を当てはめた結果
なお、ロジスティック曲線において、S字カーブのちょうど真中の座標は(x, y) = (1/2, α/β)となることを覚えておくべきである。これによれば、簡単にLD50なりED50を求めることができるからである。今回の例だと、致死量であるLD50は(1/2, 37.889/6.403) = (0.5, 5.917382)と計算できるから、LD50 = 5.917といえる。
■ Rのglm()とシンプレックス法
先の例ではRのglm()を利用してロジスティック曲線を描くための関数のパラメータを求めたが、新婦レックス法によるパラメータ推定を利用することもできる。両者の使い分けの注意点を説明するために、ここでは文献【1】のp103-表4.1のデータを両方のやり方で解析してみることにする。
| いびきの程度 | 心臓病あり | 心臓病なし | ありの割合 |
| まったくない（0） | 24 | 1355 | 0.017 |
| ときどき（2） | 35 | 603 | 0.055 |
| ほぼ毎晩（4） | 21 | 192 | 0.099 |
| 毎晩（5） | 30 | 224 | 0.118 |
まずはx軸にいびきの程度を、y軸に心臓病ありの割合をとってプロットしてみる（図3）。
図3 いびきと心臓病の割合
グラフを見てみると、そもそもこれはロジスティック曲線を当てはめるより直線を当てはめたほうがよさそうである（実際、このデータポイントによりフィットするのは線形回帰モデルだろう）。ただ、ここでは練習のためにロジスティック曲線を当てはめてみることにする。
まずは先ほど紹介したRのglm()を用いた解析例を示す。これによって得られたロジスティック曲線を当てはめたものが図4である。
> Y
[,1] [,2]
[1,] 24 1355
[2,] 35 603
[3,] 21 192
[4,] 30 224
> Snorer
[1] 0 2 4 5
> result <- glm(Y ~ Snorer, binomial)
> summary(result)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.86625 0.16621 -23.261 < 2e-16 ***
Snorer 0.39734 0.05001 7.945 1.94e-15 ***
---
> (y <- predict(result, type="response"))
> points(Snorer, y, pch=3, col="red", cex=1.5)
> x.jiku <- seq(0, 5, 0.1)
> y.jiku <- 1 / (1 + exp(3.866 - 0.397*x.jiku))
> lines(x.jiku, y.jiku, col="green", lwd=2)
図4 Rのglm()で求めたロジスティック曲線
続いてシンプレックス法によって求めたロジスティック曲線を当てはめてみる。シンプレックス法によって得た曲線を当てはめて描いたものが図5である。
図5 シンプレックス法で求めたロジスティック曲線
図4と図5を見比べると明らかにシンプレックス法で求めたロジスティック曲線のほうがよくフィットしていることが確認できる。この場合はシンプレックス法で求めた方がよいのだろうか。
ここでx軸とy軸の範囲を拡張してみる。x軸の値を0から25までとし、y軸の値を0.0から0.1までとして心臓病ありの割合をプロットする。そして、Rのglm()によって得られた曲線とシンプレックス法によって得られた曲線を2本描いてみる（図5）。すると、たしかに現在グラフウィンドウ上にプロットされているデータポイントへの当てはまりはsimplex（青いライン）のほうがよさそうである。
ところが比率というのは、理論的には必ず0から1の範囲の値をとる（0は1つも起きない、1はすべて起きるといったように）。したがって、生存率を求める場合などにおいては、シンプレックス法ではなくてRのglm()を用いるべきである。理由は図5を見れば明らかだろうが、シンプレックス法では分子を1ではなくγとしておき、パラメータガンマを任意に設定するため、今回のように1に到達しないような曲線が得られたり、あるいは逆に1を越えるような曲線が得られたりするからである。これはパラメータを推定する際にγを1であるとしてもダメである。
まとめると、理論的に0から1の範囲で値をとるようなデータの場合はRのglm()を用いるべきである。そうではない場合にはシンプレックス法によってパラメータを推定することが望ましいといえる。
図6 glmとsimplexの違い
■ ロジスティック回帰分析モデルにおける独立変数
ロジスティック回帰分析モデルを用いるにあたって、独立変数にはどのような種類の変数を用いてよいのか疑問に思うかもしれない。先に結論をいうと、どのような変数を用いてもよい。例えば、今回の最初の例では1つの連続型変数を用いたが、後の例では1つの順序データの変数を用いた（これは少し特殊な例である）。
このほかにも連続型の変数を2つ用いてもよいし、連続型とカテゴリカル型が混在していてもよい。そのような場合、すなわちロジスティックモデルにおいて複数の独立変数を用いるようなモデルを特に多重ロジスティック回帰とよぶ。
例えば、2値データであるYに対してX1, X2, ..., Xnという独立変数があったとして、
Y = X1 + X2 + X1*X2 （X1とX2はカテゴリカル型）
というモデル、あるいは、
Y = X1 + X2 + X1*X2 （X1とX2は連続型）
というモデル、さらに、
Y = X1 + X2 + X1*X2 （X1はカテゴリカル型でX2は連続型）
といったモデルでも扱うことができる。従属変数が連続型であったとき（線形回帰モデル）は上から順に分散分析モデル、重回帰モデル、共分散分析モデルという区別があった。ところが、従属変数が2値データである場合（ロジスティックモデル）ではどれもまとめて多重ロジスティック回帰モデルといって、特に区別しない。
参考文献
【1】 Alan Agresti著, 渡邉ほか訳「カテゴリカルデータ解析入門」サイエンティスト社（2005）
【2】 Annette J. Dobson著, 田中ほか訳「一般化線形モデル入門」共立出版（2008）
【3】 桶井・桶井「図解でわかる回帰分析」日本実業出版社（2004）
