一般化線形分析
佐藤圭一
2023-09-12
1 ポワソン回帰分析
1.1 例題
ある植物の種子数について調べたデータがあるとする。この種子数(y)に対する（１）体サイズ(x)・（２）施肥処理（肥料のあり・なし）(f)との関係について調べたい。種子数を従属変数として、ポワソン回帰分析を行う。
データの読み込み
# データを読み込む
dat <- read.csv("seed.data.csv", header = T)
dat$f <- factor(dat$f)
# y（種子数）の分布を見てみる
hist(dat$y)
# 種子数（y)と体サイズ（x）の関係
plot(
y = dat$y, x = dat$x,
col = dat$f) # 肥料の有無で色を変える 
ポワソン回帰分析
out <- glm(y ~ x + f,
data = dat,
family = poisson) # 正確に書いた場合にはfaimilly = poisson(link = "log")
summary(out)##
## Call:
## glm(formula = y ~ x + f, family = poisson, data = dat)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3977 -0.7337 -0.2023 0.6795 2.4317
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.26311 0.36963 3.417 0.000633 ***
## x 0.08007 0.03704 2.162 0.030620 *
## fT -0.03200 0.07438 -0.430 0.667035
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 89.507 on 99 degrees of freedom
## Residual deviance: 84.808 on 97 degrees of freedom
## AIC: 476.59
##
## Number of Fisher Scoring iterations: 4モデルの比較
library(texreg)
# モデル比較
out0 <- glm(y ~ 1, data = dat, family = poisson)
out1 <- glm(y ~ x, data = dat, family = poisson)
out2 <- glm(y ~ f, data = dat, family = poisson)
out3 <- glm(y ~ x + f, data = dat, family = poisson)
screenreg(list(out0, out1, out2, out3))##
## ==================================================================
## Model 1 Model 2 Model 3 Model 4
## ------------------------------------------------------------------
## (Intercept) 2.06 *** 1.29 *** 2.05 *** 1.26 ***
## (0.04) (0.36) (0.05) (0.37)
## x 0.08 * 0.08 *
## (0.04) (0.04)
## fT 0.01 -0.03
## (0.07) (0.07)
## ------------------------------------------------------------------
## AIC 477.29 474.77 479.25 476.59
## BIC 479.89 479.98 484.46 484.40
## Log Likelihood -237.64 -235.39 -237.63 -235.29
## Deviance 89.51 84.99 89.48 84.81
## Num. obs. 100 100 100 100
## ==================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.052 ロジスティック回帰分析
2.1 分析例
データの読み込み
# データの読み込み
dat2 <- read.csv("seed.data_logit.csv", header = T)
plot(y = dat2$y, x = dat2$x)
# 半数以上の種が生きていれば１、そうでなければ０というベクトルをつくる
y_bin <- rep(NA, nrow(dat2))
y_bin <- ifelse(dat2$y/dat2$N>=0.5, 1, 0)
dat2$y_bin <- y_bin
# データを見てみる
head(dat2)## N y x f y_bin
## 1 8 1 9.76 C 0
## 2 8 6 10.48 C 1
## 3 8 5 10.83 C 1
## 4 8 6 10.94 C 1
## 5 8 1 9.37 C 0
## 6 8 1 8.81 C 0ロジスティック回帰分析
out2 <- glm(y_bin ~ x + f, data = dat2,
family = binomial(link = "logit"))
# 結果
summary(out2)##
## Call:
## glm(formula = y_bin ~ x + f, family = binomial(link = "logit"),
## data = dat2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.76402 -0.06042 0.02537 0.15247 2.38021
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -52.466 12.803 -4.098 4.17e-05 ***
## x 5.346 1.306 4.095 4.23e-05 ***
## fT 4.944 1.421 3.480 0.000501 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 122.173 on 99 degrees of freedom
## Residual deviance: 34.675 on 97 degrees of freedom
## AIC: 40.675
##
## Number of Fisher Scoring iterations: 8オッズ比を知りたい場合
exp(out2$coefficients) # 指数化する## (Intercept) x fT
## 1.638198e-23 2.098506e+02 1.403971e+022.2 参考：二値データではないデータに対するロジスティック回帰分析
今回のデータのような形式の場合には、上記のように二値に変換しなくとも、元のデータにそのままロジットを当てはめることもできる。(このデータの場合、今回のようなデータに対しては、このような方法の方が情報ロスがなく、本来望ましい。
dat2$y_live <- dat2$y
dat2$y_dead <- dat2$N - dat2$y # 全体の個数から生きている数を引いて、死んでいる数を計算
# 従属変数（y_live, y_dead)のデータ形式を見てみる
head(dat2[ ,c("y_live", "y_dead")])## y_live y_dead
## 1 1 7
## 2 6 2
## 3 5 3
## 4 6 2
## 5 1 7
## 6 1 7# 分析
out3 <- glm(cbind(y_live, y_dead) ~ x + f, data = dat2,
family = binomial(link = "logit"))
summary(out3)##
## Call:
## glm(formula = cbind(y_live, y_dead) ~ x + f, family = binomial(link = "logit"),
## data = dat2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2584 -0.7948 0.1753 0.8757 3.1589
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -19.5361 1.4138 -13.82 <2e-16 ***
## x 1.9524 0.1389 14.06 <2e-16 ***
## fT 2.0215 0.2313 8.74 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 499.23 on 99 degrees of freedom
## Residual deviance: 123.03 on 97 degrees of freedom
## AIC: 272.21
##
## Number of Fisher Scoring iterations: 53 順序ロジスティック回帰分析と多項ロジスティック回帰分析
模擬データの作成
x <- c(8.0, 11.2, 18.7, 23.9, 12.8, 21.4, 8.9, 4.6, 13.6, 6.4)
y <- c("a", "a", "a", "b", "b", "b", "b", "c", "c", "c")
sample.dat <- data.frame(x = x, y = y)3.1 多項ロジスティック回帰分析
library(nnet)
out4 <- multinom(y ~ x, data = sample.dat)## # weights: 9 (4 variable)
## initial value 10.986123
## iter 10 value 8.862490
## final value 8.862107
## convergedsummary(out4)## Call:
## multinom(formula = y ~ x, data = sample.dat)
##
## Coefficients:
## (Intercept) x
## b -1.525252 0.1237146
## c 2.332104 -0.2323193
##
## Std. Errors:
## (Intercept) x
## b 2.194618 0.1424393
## c 2.390913 0.2326663
##
## Residual Deviance: 17.72421
## AIC: 25.724213.2 順序ロジスティック回帰分析
# 先に従属変数を順序尺度で定義する
sample.dat$y.order <- factor(sample.dat$y,
levels = c("a", "b", "c"),
ordered = T)
sample.dat$y.order## [1] a a a b b b b c c c
## Levels: a < b < c# 分析
library(MASS)
out5 <- polr(y.order ~ x, data = sample.dat,
method = "logistic") # リンク関数
summary(out5)## Call:
## polr(formula = y.order ~ x, data = sample.dat, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## x -0.07918 0.09412 -0.8413
##
## Intercepts:
## Value Std. Error t value
## a|b -2.0015 1.5548 -1.2873
## b|c -0.2043 1.4046 -0.1454
##
## Residual Deviance: 21.04751
## AIC: 27.04751