glm関数の使い方

R で一般化線形モデルを利用してデータを分析する場合は、glm 関数を利用する。

説明変数が一つの時

説明変数が一つの時のモデルの作成方法。例えば、「温度(x)とイネの伸長速度(y)」、「A 群のマウスの体重(x)と B 群のマウスの体重(y)」などのようなデータを glm で一般化線形モデルにあてはめる。

以下、正規分布に従う乱数を発生させて、サンプルデータとして使う。

x <- rnorm(20, 20, 2)     # 平均10、標準偏差2の乱数を20個つくる
y <- rnorm(20, 10, 2)     # 平均20、標準偏差2の乱数を20個つくる

# yを目的変数とし、xを説明変数としてモデルを作成
r <- glm(y ~ x)

# 解析結果を確認
summary(r)
## Call:
## glm(formula = y ~ x)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2242  -1.2713  -0.3496   0.6797   4.7749  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  18.9721     4.2603   4.453 0.000307 ***
## x            -0.4572     0.2105  -2.172 0.043463 *  
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## (Dispersion parameter for gaussian family taken to be 3.111351)
## 
##     Null deviance: 70.682  on 19  degrees of freedom
## Residual deviance: 56.004  on 18  degrees of freedom
## AIC: 83.351
## 
## Number of Fisher Scoring iterations: 2

「Coefficients:」の「t value」と「Pr(>|t|)」は、左の項目の係数がゼロに等しいかどうかの検定結果を示している。例えば、「(Intercept)」がゼロになる確率を示している。この場合は 0.000307 と非常に可能性が小さいことわかる。また、x の Pr 値は 0.043463 であるから、危険率 5% で考えた時、x の係数がゼロにならないということになる。つまり、y = β0 + β1x のモデルにおいて、x と y には相関がある。

また、β0 と β1 の値は、「Coefficients:」の「Estimate Std.」からわかる。この例では、切片が 18.9721、x の係数が -0.4572 であることがわかる。実際に散布図と回帰線を書いてみると、次のようになることが確認できる。

plot(x, y)                                            # 散布図
abline(a = r$coefficients[1], b = r$coefficients[2])  # 回帰線
glm関数による回帰分析

説明変数が2つの場合

説明変数が一つの場合、y = β0 + β1x だけを考えればよいが、説明変数が二つある場合は y = β0 + β1x1 + β2x2 のモデルにおいて、次のような場合が考えられる。

  1. β1 = 0, β2 = 0
  2. β1 ≠ 0, β2 = 0
  3. β1 = 0, β2 ≠ 0
  4. β1 ≠ 0, β2 =≠0

このほかに、y = β0 + β1x1 × x2 などのモデルが考えれる。

# サンプルデータを作成
# y が x1 + x2 で説明されるような数値にする。
x1 <- rnorm(100, 10, 2)
x2 <- rnorm(100, 20, 2)
y <- x1 + x2 + rnorm(100)

y と x1 + x2 をプロットしてみると以下のようになる。

GLMモデル選択とAICの説明用

次に、glm 関数を利用して様々なモデルを構築する。説明変数が増えると、仮定できるモデルの種類も増えていきます。複数のモデルの中から、最適なモデルを選択するときに AIC と呼ばれる指標がよく使われる。一般的に、AIC が低いほどモデルが良いとされている。例えば、サンプルデータは y = x1 + x2 となるように作ってあるので、以下のモデルのうち、model4 がもっともフィットすると考えられる。実際に AIC の値を見ると、確かに model4 の AIC が最も低いことがわかる。

model1 <- glm(y ~ 1)
model2 <- glm(y ~ x1)
model3 <- glm(y ~ x2)
model4 <- glm(y ~ x1 + x2)
model5 <- glm(y ~ x1 * x2)

# 各モデルの AIC を計算
AIC(model1)
## [1] 539.9367
AIC(model2)
## [1] 474.5637
AIC(model3)
## [1] 421.4786
AIC(model4)
## [1] 270.5629
AIC(model5)
## [1] 270.9836