GLM

自然界の変数は連続値ばかりではない。

個体数は整数しかとらないし、生存/死亡は水準2のカテゴリ変数である。

このような変数でも目的変数に取れるのが、GLM（一般化線形モデル）である。

例題1

実験

ある昆虫の卵の孵化率が温度によってどう変化するか調べるため、卵を10℃、20℃、30℃の3つの温度条件にそれぞれ30個ずつ置き、1週間後にそれぞれの卵が孵化していたかどうか記録した。（ある卵が孵化するかどうかは、他の卵に左右されないとする）

その結果、以下のようなデータが得られた。

データファイル（csv）

> hatching.data

   Hatching Temp.factor Temp.numeric

1         0        Cold           10

2         0        Cold           10

3         0        Cold           10

4         1        Cold           10

5         0        Cold           10

6         0        Cold           10

7         1        Cold           10

8         0        Cold           10

9         0        Cold           10

10        0        Cold           10

11        0        Cold           10

12        0        Cold           10

13        0        Cold           10

14        0        Cold           10

15        0        Cold           10

16        0        Cold           10

17        0        Cold           10

18        1        Cold           10

19        0        Cold           10

20        0        Cold           10

21        1        Cold           10

22        0        Cold           10

23        0        Cold           10

24        0        Cold           10

25        0        Cold           10

26        0        Cold           10

27        0        Cold           10

28        0        Cold           10

29        0        Cold           10

30        0        Cold           10

31        1         Mid           20

32        1         Mid           20

33        1         Mid           20

34        1         Mid           20

35        1         Mid           20

36        1         Mid           20

37        1         Mid           20

38        1         Mid           20

39        1         Mid           20

40        1         Mid           20

41        1         Mid           20

42        1         Mid           20

43        1         Mid           20

44        1         Mid           20

45        1         Mid           20

46        1         Mid           20

47        1         Mid           20

48        1         Mid           20

49        1         Mid           20

50        1         Mid           20

51        1         Mid           20

52        0         Mid           20

53        1         Mid           20

54        1         Mid           20

55        1         Mid           20

56        1         Mid           20

57        1         Mid           20

58        1         Mid           20

59        1         Mid           20

60        1         Mid           20

61        0         Hot           30

62        1         Hot           30

63        1         Hot           30

64        1         Hot           30

65        0         Hot           30

66        1         Hot           30

67        1         Hot           30

68        0         Hot           30

69        1         Hot           30

70        0         Hot           30

71        1         Hot           30

72        0         Hot           30

73        1         Hot           30

74        1         Hot           30

75        1         Hot           30

76        0         Hot           30

77        0         Hot           30

78        1         Hot           30

79        0         Hot           30

80        0         Hot           30

81        1         Hot           30

82        0         Hot           30

83        1         Hot           30

84        1         Hot           30

85        0         Hot           30

86        1         Hot           30

87        0         Hot           30

88        1         Hot           30

89        1         Hot           30

90        1         Hot           30

> 

ここで、Hatchingのゼロイチデータは、0が孵化しなかったこと、1が孵化したことを示している。

また、Temp.factorとTemp.numericはそれぞれ温度処理をカテゴリ変数と連続変数として入力したものである。

実際、Temp.factorは要因ベクトル、Temp.numericは数値ベクトルである。

> is.factor(hatching.data$Temp.factor)

[1] TRUE

> is.numeric(hatching.data$Temp.numeric)

[1] TRUE

> 

両者は同じものを違う形で表現しているだけなのだが、数値か要因かというのは実践的にはとても重大な違いである。

本稿ではそれも合わせて説明する。

さて、このようなデータから、温度が孵化率に影響を与えているか検証するためには、どのような解析を行ったらよいだろうか？

むしろ、元のゼロイチデータをそのままプロットした方がわかりやすいだろう。

par(cex=1.5, las=1, lwd=2, family="sans")

plot(Hatching ~ jitter(Temp.numeric), hatching.data)

points(c(10, 20, 30), 

       with(hatching.data, tapply(Hatching, Temp.factor, function(x) sum(x)/length(x))), #孵化率計算

       pch=16, cex=3)

ここで、jitter関数は変数をランダムに”揺らす”関数で、同じ値が何個も重なっているときの図示に有効である。

また、points関数は点を上書きする関数で、ここでは孵化数/卵数のようにそれぞれの温度グループの孵化率を計算して上書きしている。

GLMの実践（説明変数が連続変数の場合）

Rでは、GLMはglm関数によって作成できる。

作り方はだいたい線形モデルと同じだが、引数familyで誤差分布とリンク関数を指定する点が異なる。

誤差分布とリンク関数は目的変数の性質によって選ぶ。

• 二値データ → 誤差分布：binomial, リンク関数：logit

• ゼロ以上の整数データ → 誤差分布：poisson, リンク関数：log

• ゼロより大きい実数データ → 誤差分布：Gamma, リンク関数：log or inverse

誤差分布・リンク関数とはなんぞや？というのはとりあえず置いておく。

温度上昇が孵化率に与える影響のモデルをGLMで作成する方法は、以下。

m.numeric <- glm(Hatching ~ Temp.numeric, hatching.data, family=binomial)

summary(m.numeric)

孵化するか否かは二値データなので、誤差分布はbinomialである。

family=binomialとするだけで、誤差分布をbinomial, リンク関数をlogitに自動的に指定できる。

保存したモデルをsummaryに渡すことで、モデルのサマリーが出力される。

注目すべきはCoefficientsの部分である。

Temp.numericのEstimateが正値（0.103）であり、そのP値は< 0.001であるから、温度と孵化する/しないとの関係の傾きは有意に正、といえる。

もう一つチェックしなければならないのはResidual devianceである。

これは線形モデルで言う残差のようなもので、自由度当たりのresidual devianceが1.5を超えると結果がちょっとアヤシイ…ということになっている（過分散）。

109.30/88 = 1.24であるから、今回は合格ということにしておく。

（ただ、今回のように実データが正真正銘0/1しかとらない二項分布GLM（Bernoulli GLM）では過分散は起こりえないようだ（McCullagh and Nelder 1989））

温度の効果が有意なのか？という問いには、尤度比検定を使うことが多い。

一番簡単なのは、car::Anovaを使う方法である。

スクリプトは以下である。

library(car)

Anova(m.numeric, test="LR")

算出される統計量はカイ二乗値である。

カイ二乗値 = 13.866, 自由度 = 1, P < 0.001なので、温度が孵化率に与える影響は有意である。

以上の解析から、温度の上昇は孵化率に正の影響を与える、と結論できる、のだが…。

結果の図示

GLMの結果の図示はやや面倒である。

なぜなら、GLMの結果は基本的にリンク関数による変換後のスケールで出力されるからである。

実データと合わせてレポートするには、これを実数スケールに戻さないといけない。

つまり、結果にリンク関数の逆関数をかます必要がある。

GLMの結果を図示する工程は以下である。

1. predict関数によって予測値とSEを計算させる。

2. 計算した予測値、予測値±SE（95%信頼区間を求めたければ、予測値±SE×1.96）をリンク関数の逆関数によって戻し変換（back transform）する。

3. lines関数によって実データプロットに上書きする。

predict関数は第一引数に予測に使うモデル、第二に予測したい説明変数のデータフレームを入れる。

また、se.fit=TRUEとすることで、SEが同時に計算される。

binomialの場合、リンク関数はロジットなのだが、これの逆関数は1 / ( 1 + exp( -η ))という形をしている（ここでηは線形予測子）。

たくさん言葉が出て来てワケがわからないかもしれないが、とりあえず、この関数のηの部分に予測値を渡せば、binomialのGLMの結果は確率スケールに戻すことができる。

#1

pred.logit <- predict(m.numeric, data.frame(Temp.numeric=10:30), se.fit=TRUE) #10℃から30℃までの孵化率（ただしロジットスケール）を予測

#2

pred.prob <- 1/(1+exp(-pred.logit$fit)) #予測値を実数スケールに戻す

pred.li <- 1/(1+exp(-(pred.logit$fit-pred.logit$se.fit*1.96))) #信頼区間下限を実数スケールに戻す

pred.ui <- 1/(1+exp(-(pred.logit$fit+pred.logit$se.fit*1.96))) #信頼区間上限を実数スケールに戻す

#3

par(cex=1.5, las=1, lwd=2, family="sans")

plot(Hatching~jitter(Temp.numeric), hatching.data, xlab="Temperature", ylab="Hatching probability")

points(c(10, 20, 30), 

       with(hatching.data, tapply(Hatching, Temp.factor, function(x) sum(x)/length(x))), #孵化率計算

       pch=16, cex=3)

lines(10:30, pred.prob, lty=1) #予測値の上書き

lines(10:30, pred.li, lty=2) #信頼区間下限の上書き

lines(10:30, pred.ui, lty=2) #信頼区間上限の上書き

上のような図を描くことができる。

実線がGLMによって予測される孵化確率、破線で囲まれた区間が予測値の95%信頼区間ということになる。

しかし…

GLMの曲線は実データといまいちフィットしていない。

たしかに孵化率は温度上昇にしたがって上昇傾向ではあるのだが、現実には20℃から30℃では孵化率はむしろ低下している。

このような「上がってから少し下がる」ような傾向をこのモデルはうまく表現しているとは言い難い。

GLMは普通の線形モデルと同じように、説明変数が連続変数の場合、説明変数と目的変数の関係は単調増加or減少を仮定している。

なので、「上がってから少し下がる」ような（非線形な）動きは表現できないのである。

さて、どうすればモデルを改善できるだろうか？

GLMの実践（説明変数がカテゴリ変数の場合）

幸い、今回の実験では、各温度条件で反復が30卵確保されている。

このようなデザインであれば、説明変数をカテゴリ変数に置き換えてモデルを作ることができる。

説明変数がカテゴリであれば、上がってから下がるような傾向を表現することも可能である。

カテゴリの水準ごとにモデルの係数が推定されるからである。

要は、ANOVAにおける線形モデルの作り方と同じである。

今度は、要因ベクトルの方を説明変数にしてみよう。

m.factor <- glm(Hatching ~ Temp.factor, hatching.data, family=binomial)

summary(m.factor)

今度は、Coefficientsのところに切片を除いて2つの係数が推定されている。

これはそれぞれ、ColdとMidの差、ColdとHotの差を表す係数である。

Midの孵化率はHotより高いので、Temp.factorMidの係数5.24の方がTemp.factorHotの係数2.28より高いのは辻褄が合っている。

尤度比検定も行おう。

Anova(m.factor, test="LR")

カイ二乗値 = 50.45, 自由度 = 2, P < 0.001なので、Temp.factorの効果は有意であるといえる。

（カテゴリの水準数が3なので、自由度が3-1 = 2になっていることに注意。説明変数が連続変数の場合と比較してみよう）

よって、温度処理は孵化率に有意な影響を与えている、と結論づけられる。

事後検定

説明変数をカテゴリにした場合、One-way ANOVAのときと全く同じ方法で事後検定も行うことができる。

スクリプトは以下である。

library(lsmeans)

lsm.factor <- lsmeans(m.factor, "Temp.factor")

lsm.factor

pairs(lsm.factor, adjust="tukey")

2行目のスクリプトでは、カテゴリの各水準の平均値（ロジットスケール）を推定している。

元のデータのとり得る値が0と1しかないので、平均値というと奇異に聞こえるかもしれないが、そういうものとして理解してほしい。

これをpairsで囲むことにより、これらの平均値を総当たりでペアごとに検定してくれる（Tukey法で有意水準を調整した）。

今回の場合、すべての組み合わせの間で有意な差があるようである。

結果の図示

これもOne-way ANOVAで説明したのとほとんど同じ方法で描くことが可能である。

つまり、lsmeansのオブジェクトをsummaryに渡してlsmeanとSEを取り出し、実データのスケール（つまり確率スケール）に戻し変換して図示すればよい。

今回はbinomialなので、先ほども説明した関数”1 / ( 1 + exp( -η ))”を使って変換する。

summary.lsm.factor <- summary(lsm.factor)

m.prob <- 1/(1+exp(-summary.lsm.factor$lsmean))

e.l.prob <- 1/(1+exp(-summary.lsm.factor$lsmean+summary.lsm.factor$SE))

e.u.prob <- 1/(1+exp(-summary.lsm.factor$lsmean-summary.lsm.factor$SE))

par(las=1, family="sans", cex=1.5, lwd=1.5)

plot(c(10, 20, 30), m.prob, pch=16, cex=2, ylim=c(0, 1.1), xlim=c(7, 33), 

     xlab="Temperature", ylab="Hatching probability", lab=c(3, 6, 7))

arrows(c(10, 20, 30), e.l.prob, c(10, 20, 30), e.u.prob, angle=90, length=0.15, code=3)

text(x=c(10, 20, 30), y=e.u.prob+0.07, labels=c("a", "b", "c"))

「上がってから少し下がる」関係がGLMによってうまく表現されていることがわかる。

また、エラーバーの長さが上下で異なることにも注意。

ロジットスケールから実数スケールに戻しているので、このようなことが起こる。

この解析結果から、卵の孵化率は温度が高ければ高いほどよくなるわけではなく、20℃でもっとも高くなり、それ以上の温度ではむしろ低下すると結論づけられる。

GLMでもこのような”ANOVA的な”解析を行うことが可能である。

まとめ：

• 0/1データはGLM（binomial）による解析が可能

• Residual devianceとresidual degree of freedomの比をチェックする必要

• 上の比が1.5未満なら、尤度比検定によって説明変数の有意性を検証

• 説明変数が連続値かカテゴリかに注意する必要

例題2

野外調査

モンシロチョウの個体数がキャベツ畑の数に影響されるか調べるため、調査地をランダムに15か所選定し、一定時間内におけるモンシロチョウの目撃数を記録した。また、15か所の調査地の半径1km範囲内に存在するキャベツ畑の数を数え上げた。

その結果、以下のようなデータが得られた。