glm 関数と回帰分析

R で一般化線形モデルによる回帰分析は lm 関数または glm 関数で行う。このページでは回帰分析を一般化線形モデルとして扱うので、glm 関数のみ取り扱う(lm 関数による回帰分析)。

単回帰

確率変数 X と Y があるとき、Y を X で説明しようとする手法。このとき、Y を従属変数、X を独立変数と呼ぶ。単回帰は一般に Y = aX + b あるいや Y = aX のモデルを立てて、それらの係数 a と b をデータから推定しようとするもの。a と b のことをモデルのパラメーターという。

例えば、2 つの確率変数 X と Y の観測値が以下のように得られたとする。この観測値を利用してモデルのパラメーターを計算する。

x <- c(22.1, 32.9, 33.2, 12.4, 19.4, 26.5, 41.2, 54.2)
y <- c(12.9, 19.1, 13.5,  9.1,  9.8, 11.1, 25.0, 30.7)

m <- glm(y ~ x, family = gaussian)

glm により分析された結果から得られるモデルは summary 関数を利用して概要を確認できる。

summary(m)
## Call:
## glm(formula = y ~ x)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.5565  -1.2076   0.9756   1.5260   2.6741  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.50777    2.63909  -0.192  0.85377    
## x            0.55917    0.08077   6.923  0.00045 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
## 
## (Dispersion parameter for gaussian family taken to be 7.998795)
## 
##     Null deviance: 431.340  on 7  degrees of freedom
## Residual deviance:  47.993  on 6  degrees of freedom
## AIC: 43.036
## 
## Number of Fisher Scoring iterations: 2

summary によって出力される各項の意味は次のようになっている。

Deviance Residuals残差平方和の四分位数
Coefficients予測したパラメーターの値と統計量など
Estimateモデルパラメーターの予測値
Std. Error標準偏差
t valueそのパラメーターが 0 かどうかの検定結果である統計量
P(>|t|)t value に対応する p value
Null devianceNull Deviance = 2(LL(飽和モデル) - LL(帰無モデル)) によって計算される。飽和モデルはデータの個数分だけのパラメーターを持つと想定したモデルを指し、帰無モデルはすべてのデータがただ一つのパラメーターを持つと想定したモデルを指す。
Residual devianceResidual Deviance = 2(LL(飽和モデル) - LL(提案モデル)) によって計算される。提案モデルは解析者が提案したモデルを指す。上の例では 2 つのパラメーター(a と b)を想定しているモデルを指す。
AIC赤池情報量規準
plot(x, y, xlim = c(0, max(x)), ylim = c(0, max(y)))
abline(a = m$coefficients[1], b = m$coefficients[2])
glm関数を利用した回帰分析

重回帰

次に、一つの従属変数 Y を複数の独立変数 P, Q, R で回帰する例を示す。乱数を利用してサンプルデータを生成する。独立変数 R を生成したが、これは直接 Y に影響を与えていないことに注意。

p <- rnorm(10, mean = 5, sd = 2)
q <- rnorm(10, mean = 10, sd = 2)
r <- rnorm(10, mean = 0, sd = 1)

y <- 2 * p + 3 *q + runif(10)

dat <- data.frame(y = y, p = p, q = q, r = r)

ここで重回帰を行うが、重回帰の場合は複数のモデルが想定できる。まず、Y がどの独立変数にも影響されないという、もっとも簡単なモデルを考えることができる。これを帰無モデルとする。次に、Y がすべての独立変数に影響されるというモデルも考えられ、これを飽和モデルとする。最後に、Y が 3 つの独立変数のうち、1 個だけあるいは 2 個だけに影響されるというモデルが考えられ、これを提唱モデルとする。

full.model <- glm(y ~ ., data = dat)
null.model <- glm(y ~ 1, data = dat)
proposed.model.1 <- glm(y ~ p + q, data = dat)
proposed.model.2 <- glm(y ~ p + r, data = dat)

ここで 4 つのモデルを想定しましたが、どのモデルがよいのかを選択するとき、AIC がよく利用される。AIC が小さいモデルがよいモデルとされる。AIC は以下のように計算できる。ただし、l はそのモデルの対数尤度を表し、p はそのモデルに含まれるパラメーター数を表す。

\[\text{AIC} = -2 l + 2 p\]
-2 * as.numeric(logLik(full.model)) + 2 * 5
## [1] -0.6187658

AIC(full.model)
## [1] -0.6187658
AIC(null.model)
## [1] 73.26417
AIC(proposed.model.1)
## [1] -2.558648
AIC(proposed.model.2)
## [1] 70.79138

各モデルの AIC を計算してみると proposed.model.1 の AIC がもっとも小さいことがわかった。つまり、想定している 4 つのモデルのうち「y = ap + bq + c」の回帰モデルがもっともよいと言える。

summary(proposed.model.1)
## Call:
## glm(formula = y ~ p + q, data = dat)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.253468  -0.051425   0.009472   0.134386   0.162452  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.08471    0.30846   0.275    0.792    
## p            1.86876    0.05256  35.554 3.62e-09 ***
## q            3.12618    0.02959 105.646 1.79e-12 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## (Dispersion parameter for gaussian family taken to be 0.02909855)
## 
##     Null deviance: 596.52717  on 9  degrees of freedom
## Residual deviance:   0.20369  on 7  degrees of freedom
## AIC: -2.5586
## 
## Number of Fisher Scoring iterations: 2