説明変数を組み込んだ統計モデルについて説明
ポアソン回帰(Poisson regression)$\cdots$個体ごとに異なる説明変数(個体の属性)によって平均種子数が変化する.
$\rightarrow$第2章ではどの個体の種子数$y_i$も平均$\lambda$のポアソン分布にしたがうと仮定
一般化線形モデル(Generalized Linear Model)$\cdots$$\uparrow$と似たような構造の統計モデルを総称してこう呼ぶ.
植物個体$\cdots$$i$
種子数$\cdots$$y_i$
体サイズ(body size)$\cdots$$x_i$ # 植物の大きさの大小を表す正の実数
処理$C$$\cdots$全個体のうち50個体($i\in{1, 2, \cdots, 50}$)は処理しない
処理$T$$\cdots$残り50個体($i\in{51, 52, \cdots, 100}$)は施肥処理(肥料を加える)する
d <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv")
d
y | x | f |
---|---|---|
6 | 8.31 | C |
6 | 9.44 | C |
6 | 9.50 | C |
12 | 9.07 | C |
10 | 10.16 | C |
4 | 8.32 | C |
9 | 10.61 | C |
9 | 10.06 | C |
9 | 9.93 | C |
11 | 10.43 | C |
6 | 10.36 | C |
10 | 10.15 | C |
6 | 10.92 | C |
10 | 8.85 | C |
11 | 9.42 | C |
8 | 11.11 | C |
3 | 8.02 | C |
8 | 11.93 | C |
5 | 8.55 | C |
5 | 7.19 | C |
4 | 9.83 | C |
11 | 10.79 | C |
5 | 8.89 | C |
10 | 10.09 | C |
6 | 11.63 | C |
6 | 10.21 | C |
7 | 9.45 | C |
9 | 10.44 | C |
3 | 9.44 | C |
10 | 10.48 | C |
... | ... | ... |
10 | 10.54 | T |
8 | 11.30 | T |
8 | 12.40 | T |
7 | 10.18 | T |
5 | 9.53 | T |
6 | 10.24 | T |
8 | 11.76 | T |
9 | 9.52 | T |
9 | 10.40 | T |
6 | 9.96 | T |
7 | 10.30 | T |
10 | 11.54 | T |
6 | 9.42 | T |
11 | 11.28 | T |
11 | 9.73 | T |
11 | 10.78 | T |
5 | 10.21 | T |
6 | 10.51 | T |
4 | 10.73 | T |
5 | 8.85 | T |
6 | 11.20 | T |
5 | 9.86 | T |
8 | 11.54 | T |
5 | 10.03 | T |
9 | 11.88 | T |
8 | 9.15 | T |
6 | 8.52 | T |
8 | 10.24 | T |
7 | 10.86 | T |
9 | 9.97 | T |
d$x
d$y
d$f
class(d)
class(d$y)
class(d$x)
class(d$f)
summary(d)
y x f Min. : 2.00 Min. : 7.190 C:50 1st Qu.: 6.00 1st Qu.: 9.428 T:50 Median : 8.00 Median :10.155 Mean : 7.83 Mean :10.089 3rd Qu.:10.00 3rd Qu.:10.685 Max. :15.00 Max. :12.400
plot(d$x, d$y, pch = c(21, 19)[d$f])
legend("topleft", legend = c("C", "T"), pch = c(21,19))
箱ひげ図(box-whisker plot)
plot(d$f, d$y)
まず最初に個体$i$の体サイズ$x_i$だけに依存する統計モデルについて考えてみる.
ある個体$i$において種子数が$y_i$である確率$p(y_i \mid \lambda_i)$はポアソン分布に従う.
$$p(y_i \mid \lambda_i) = \frac{\lambda^{y_i}_i exp(-\lambda_i)}{y_i!}$$個体ごとに異なる平均$\lambda_i$を説明変数$x_i$の関数として定義しなければならない.
ある個体$i$の平均種子数$\lambda_i$を
$\lambda_i = exp(\beta_1 + \beta_2x_i)$
とする. この時,
curve(exp(-2 + -0.8 * x), from = -4 , to = 5, lwd = 4, lty = 2)
curve(exp(-1 + 0.4 * x), add = TRUE, lwd = 4)
abline(v = 0:0)
text(-2,3,expression(paste(group("{",list(beta[1],beta[2]) ,"}") == group("{",list(-2, 0.8),"}"))))
text(2,2,expression(group("{",list(beta[1],beta[2]) ,"}") == group("{",list(-1, 0.4),"}")))
上の式を変形して,
$log\lambda_i = \beta_1 + \beta_2x_i$
この時,
データ$\bf Y$のもとでの対数尤度
$$log L(\beta_1, \beta_2) = \sum_{i}log \frac{\lambda^{y_i}_i exp(-\lambda_i)}{y_i!}$$
複数のパラメータ$\{\beta_1, \beta_2\}$を同時に扱うので, 最尤推定量の導出は難しい.
→実際のポアソン回帰は大抵数値的な試行錯誤で最尤推定値を探し出す
→解析的に導出出来なくても問題ない
→Rだと簡単にGLMのあてはめが出来る
fit <- glm(y ~ x, data = d, family = poisson)
fit
Call: glm(formula = y ~ x, family = poisson, data = d) Coefficients: (Intercept) x 1.29172 0.07566 Degrees of Freedom: 99 Total (i.e. Null); 98 Residual Null Deviance: 89.51 Residual Deviance: 84.99 AIC: 474.8
summary(fit)
Call: glm(formula = y ~ x, family = poisson, data = d) Deviance Residuals: Min 1Q Median 3Q Max -2.3679 -0.7348 -0.1775 0.6987 2.3760 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 1.29172 0.36369 3.552 0.000383 *** x 0.07566 0.03560 2.125 0.033580 * --- 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.993 on 98 degrees of freedom AIC: 474.77 Number of Fisher Scoring iterations: 4
f.b1 <- function(x){(1/(sqrt(2*pi)))*exp(-(x-3.552)^2/2)}
integrate(f.b1,-Inf,0)$value
integrate(f.b1,-Inf,0)$value * 2
#グラフ
plot(f.b1,from = -2, to = 1)
xvals <- seq(-2,0,length = 10)
yvals <- f.b1(xvals)
polygon(c(xvals,rev(xvals)),c(rep(0,10),rev(yvals)),col="grey")
f.b2 <- function(x){(1/(sqrt(2*pi)))*exp(-(x-2.125)^2/2)}
integrate(f.b2,-Inf,0)$value
integrate(f.b2,-Inf,0)$value * 2
# グラフ
plot(f.b2, from = -2, to = 1)
xvals <- seq(-2,0,length = 10)
yvals <- f.b2(xvals)
polygon(c(xvals,rev(xvals)),c(rep(0,10),rev(yvals)),col="grey")
logLik(fit)
'log Lik.' -235.3863 (df=2)
$\lambda = exp(1.29 + 0.0757x)$
を使う
plot(d$x, d$y, pch = c(21, 19)[d$f])
xx <- seq(min(d$x), max(d$x), length = 100)
#lines(xx, exp(1.29 + 0.0757 * xx), lwd = 2)
# 同じ事をpredict()関数を使ってもできる
yy <- predict(fit, newdata = data.frame(x = xx), type = "response")
lines(xx, yy, lwd = 2, lty=2)
ここで, 無視していた施肥効果$f_i$を説明変数として組み込んだモデルを検討する.
植物の体サイズ$x_i$の効果を無視して施肥効果$f_i$だけが影響するモデルの平均値を
$\lambda_i = exp(\beta_1 + \beta_3d_i)$
として, 説明変数を$f_i$ではなく$d_i$として, 以下の値を取る.
\begin{eqnarray}
d_i=\left\{ \begin{array}{ll}
0 & (f_i = Cの場合) \\
1 & (f_i = Tの場合) \\
\end{array} \right.
\end{eqnarray}
fit.f <- glm(y ~ f, data = d, family = poisson)
fit.f
Call: glm(formula = y ~ f, family = poisson, data = d) Coefficients: (Intercept) fT 2.05156 0.01277 Degrees of Freedom: 99 Total (i.e. Null); 98 Residual Null Deviance: 89.51 Residual Deviance: 89.48 AIC: 479.3
logLik(fit.f)
'log Lik.' -237.6273 (df=2)
説明変数が$x_i$だけのときよりもあてはまりが悪い.
重回帰(multiple regression)$\cdots$複数の説明変数を持つ統計モデルによるあてはめ
fit.all <- glm(y ~ x + f, data = d, family = poisson)
fit.all
Call: glm(formula = y ~ x + f, family = poisson, data = d) Coefficients: (Intercept) x fT 1.26311 0.08007 -0.03200 Degrees of Freedom: 99 Total (i.e. Null); 97 Residual Null Deviance: 89.51 Residual Deviance: 84.81 AIC: 476.6
logLik(fit.all)
'log Lik.' -235.2937 (df=3)
説明変数が$x_i$だけの時よりはあてはまりが良い. →比較は第4章で検討