Stan の私的メモ。Stan のインストールはこちら。
今回は一般化線形モデルを glm() と Stan によるベイジアンモデリングで対比させる形で書きました。
一般線形モデル
一般線形モデル(General linear model), 又は単に線形モデル(linear model)は線形予測子が正規分布, リンク関数が恒等リンク関数の場合の回帰モデルです。
以下の架空のデータを使います。
d <- data.frame(
x = c(1,2,3,4,5,6,7,8,9,10),
y = c(9,8,10,9,11,12,11,12,14,15)
)
lm()の場合,
d.lm <- lm(y ~ x, data = d)glm()の場合,
d.glm <- glm(
y ~ x,
family = gaussian(link = "identity"),
data = d
)
lm(), glm()の結果は以下
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.33333 0.60168 12.188 1.9e-06 ***
x 0.68485 0.09697 7.063 0.000106 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for gaussian family taken to be 0.7757576)
Null deviance: 44.9000 on 9 degrees of freedom
Residual deviance: 6.2061 on 8 degrees of freedom
AIC: 29.608
stanによるベイジアンモデリングの場合,
data {
int N;
vector[N] x;
vector[N] y;
}
parameters {
real alpha;
real e;
real s;
}
model {
for(i in 1:N)
y[i] ~ normal(x[i] * alpha + e, s);
}
d.list <- list(
N = 10,
x = d$x,
y = d$y
)
d.fit <- stan(
file = 'model/normal.stan',
data = d.list,
iter = 1000,
chains = 4
)
stanの結果は以下
Inference for Stan model: stan_code.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha 0.68 0.01 0.12 0.46 0.61 0.68 0.75 0.93 422 1.01
e 7.33 0.04 0.73 5.87 6.89 7.33 7.77 8.76 423 1.00
s 1.04 0.02 0.32 0.63 0.82 0.97 1.19 1.79 419 1.01
lp__ -2.55 0.07 1.43 -6.44 -3.11 -2.18 -1.53 -0.97 377 1.01
一般化線形モデル
一般化線形モデル (Generalized linear model, GLM)は, 正規分布以外の分布を線形予測子に指定できるように拡張したモデルです。
ロジスティック回帰
ロジスティック回帰で使う架空のデータ。y は何かしらの生存数のようなイメージ。
d <- data.frame(
N = 10,
x = c(1,2,3,4,5,6,7,8,9,10),
y = c(0,0,1,2,2,4,6,6,8,9)
)
線形予測子にbinomial, リンク関数にlogitを指定したロジスティック回帰。目的変数には, 0/1 の値を持つ変数, または cbind() で度数を表す変数を指定する。
d.glm <- glm(
cbind(y, N-y) ~ x,
family = binomial(link = "logit"),
data = d
)
glm()の結果は以下
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.4802 0.8666 -5.170 2.34e-07 ***
x 0.6586 0.1273 5.173 2.31e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 49.4047 on 9 degrees of freedom
Residual deviance: 2.2893 on 8 degrees of freedom
AIC: 25.565
Number of Fisher Scoring iterations: 4
stanの場合,
data {
int N;
int x[N];
int y[N];
}
parameters {
real beta0;
real beta1;
}
model {
for(i in 1:N)
y[i] ~ binomial(N, inv_logit(beta0 + beta1 * x[i]));
}
d.list <- list(
N = 10,
x = d$x,
y = d$y
)
d.fit <- stan(
file = 'model/binomial.stan',
data = d.list,
iter = 1000,
chains = 4
)
stanの結果は以下。
Inference for Stan model: binomial.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
beta0 -4.63 0.06 0.90 -6.59 -5.19 -4.56 -3.99 -3.06 270 1.01
beta1 0.68 0.01 0.13 0.45 0.59 0.67 0.76 0.97 257 1.01
lp__ -43.87 0.06 1.03 -46.65 -44.29 -43.55 -43.14 -42.88 327 1.02
ポアソン回帰
ポアソン回帰で使う架空のカウントデータ。
d <- data.frame(
x = c(1,2,3,4,5,6,7,8,9,10),
y = c(44,32,26,20,14,9,7,4,2,1)
)
線形予測子にpoisson, リンク関数にlogを指定したポアソン回帰,
d.glm <- glm(
y ~ x,
family = poisson(link = "log"),
data = d
)
glm()の結果は以下
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.19000 0.13638 30.723 <2e-16 ***
x -0.33717 0.03535 -9.537 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 118.9260 on 9 degrees of freedom
Residual deviance: 2.6949 on 8 degrees of freedom
AIC: 47.803
stanの場合,
data {
int N;
int x[N];
int y[N];
}
parameters {
real beta0;
real beta1;
}
model {
for(i in 1:N)
y[i] ~ poisson_log(beta0 + beta1 * x[i]);
}
d.list <- list(
N = 10,
x = d$x,
y = d$y
)
d.fit <- stan(
file = 'model/poisson.stan',
data = d.list,
iter = 1000,
chains = 4
)
stanの結果は以下
Inference for Stan model: stan_code.
4 chains, each with iter=1000; warmup=500; thin=1;
post-warmup draws per chain=500, total post-warmup draws=2000.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
beta0 4.18 0.01 0.14 3.90 4.09 4.19 4.28 4.46 514 1.00
beta1 -0.34 0.00 0.04 -0.41 -0.36 -0.34 -0.31 -0.27 477 1.00
lp__ 337.93 0.05 1.10 335.03 337.57 338.24 338.66 338.94 430 1.01
一般化線形混合モデル
一般化線形混合モデル (Generalized linear mixed model, GLMM)は, 一般化線形モデルの拡張モデルです。
個体差のあるロジスティック回帰
みどり本にでてくる植物の種子データを使います。
d <- as.data.frame(read.csv(file = "data/data.csv"))
# > head(d)
# N y x id
# 1 8 0 2 1
# 2 8 1 2 2
# 3 8 2 2 3
# 4 8 4 2 4
# 5 8 1 2 5
# 6 8 0 2 6
データの概略です。
- N : 各植木鉢に植えられている植物がもつ種子の総数
- y : 生存種子の数
- x : 各植木鉢に植えられている植物ごとの葉の数
- id: 個体ID
{glmmML} の glmmML() を使います。
d.glmm <- glmmML(
cbind(y, N-y) ~ x,
family = binomial,
cluster = id,
data = d
)
glmmML()の結果が以下。
coef se(coef) z Pr(>|z|)
(Intercept) -4.190 0.8777 -4.774 1.81e-06
x 1.005 0.2075 4.843 1.28e-06
Scale parameter in mixing distribution: 2.408 gaussian
Std. Error: 0.2202
LR p-value for H_0: sigma = 0: 2.136e-55
Residual deviance: 269.4 on 97 degrees of freedom AIC: 275.4
stanの場合,
data {
int N;
int n;
int x[n];
int y[n];
}
parameters {
real beta0;
real beta1;
vector[n] r;
real s_r;
}
model {
for(i in 1:n) {
y[i] ~ binomial(N, inv_logit(beta0 + beta1 * x[i] + r[i]));
r[i] ~ normal(0, s_r);
}
}
d.list <- list(
N = 8,
n = 100,
x = d$x,
y = d$y
)
d.fit <- stan(
file = 'model/glmm.stan',
data = d.list,
iter = 1000,
chains = 4
)
stanの結果が以下。Scale parameterの一致性が少し微妙。
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
beta0 -4.20 0.06 0.98 -6.23 -4.83 -4.15 -3.55 -2.36 308 1.00
beta1 1.01 0.01 0.23 0.58 0.85 1.01 1.16 1.50 315 1.01
r[1] -2.06 0.04 1.76 -5.93 -3.15 -1.88 -0.82 0.92 2000 1.00
r[2] -0.12 0.03 1.19 -2.58 -0.85 -0.05 0.72 2.00 1200 1.00
r[3] 0.81 0.03 1.01 -1.31 0.17 0.88 1.48 2.69 1053 1.00
r[4] 2.01 0.03 0.89 0.30 1.42 1.99 2.60 3.79 758 1.00
r[5] -0.16 0.04 1.25 -3.15 -0.82 -0.05 0.69 2.03 1108 1.00
r[6] -2.05 0.04 1.74 -5.69 -3.20 -1.92 -0.82 0.92 2000 1.00
r[7] -2.06 0.04 1.78 -6.15 -3.10 -1.91 -0.77 0.91 2000 1.00
r[8] 3.81 0.04 1.07 1.88 3.08 3.73 4.44 6.08 746 1.00
r[9] -0.13 0.03 1.17 -2.66 -0.83 -0.08 0.65 1.85 1339 1.00
r[10] 3.08 0.03 0.96 1.30 2.43 3.05 3.71 5.01 889 1.00
r[11] -2.08 0.04 1.81 -6.23 -3.10 -1.94 -0.83 0.77 2000 1.00
r[12] -2.09 0.04 1.80 -6.18 -3.22 -1.90 -0.81 0.80 2000 1.00
r[13] -2.07 0.04 1.82 -5.97 -3.17 -1.84 -0.78 0.88 2000 1.00
r[14] -2.11 0.04 1.83 -6.04 -3.24 -1.92 -0.83 0.98 2000 1.00
r[15] -2.09 0.04 1.84 -6.25 -3.18 -1.91 -0.80 0.99 2000 1.00
r[16] -0.15 0.03 1.18 -2.62 -0.90 -0.06 0.68 1.96 2000 1.00
r[17] 1.46 0.03 0.94 -0.43 0.83 1.47 2.08 3.30 779 1.00
r[18] 2.51 0.03 0.89 0.83 1.90 2.51 3.10 4.30 811 1.00
r[19] 1.44 0.03 0.94 -0.52 0.85 1.47 2.08 3.21 784 1.00
r[20] 3.80 0.04 1.10 1.86 3.05 3.71 4.49 6.20 976 1.00
...
s_r 2.62 0.02 0.33 2.02 2.39 2.60 2.82 3.32 408 1.01
lp__ -441.60 0.50 9.59 -461.04 -447.78 -441.19 -434.75 -424.00 370 1.01
CodeはGitHubにあります。
[1] R で個体差のあるロジスティック回帰1 - glmmML
[2] 分布から見た線形モデル・GLM・GLMM
[3] Stanのマニュアルの8章～12章の私的メモ
[4] Bayesian regression with STAN Part 2: Beyond normality