😺
ベイジアンlog-binomial回帰の実装
ゴール
- Bayesian log-binomialモデルの実装を行う
やること
- log-binomialモデルのメトロポリス・ヘイスティング法による実装
-
glm()によるRRのとの比較 - log-binomilモデルの回帰標準化によるリスク差の推定
データセットの作成
最近公開している記事で用いたデータセットを引き続き使用する。
N <- 500 # sample size
v1 <- rnorm(N, mean = 60, 10)
v2 <- rbinom(N, size = 1, prob = 0.5)
v3 <- rbinom(N, size = 1, prob = 0.5)
v4 <- rbinom(N, size = 1, prob = 0.5)
X <- do.call(cbind, list(1,v1,v2,v3,v4))
colnames(X) <- c('const','v1', 'v2', 'v3', 'v4')
P <- ncol(X)
true_beta <- c(-2.2, 0, 0.8, 0, 0) # true beta in the model
Y <- unlist(lapply(invlogit(X%*%true_beta), rbinom, n = 1, size = 1))
mean(Y)
X[,2] <- scale(X[,2])
mydata <- data.frame(X[,-1],Y)
p1_true <- invlogit(t(c(1,0,1,0,0))%*%c(-2.2, 0, 0.8, 0, 0))
p0_true <- invlogit(t(c(1,0,0,0,0))%*%c(-2.2, 0, 0.8, 0, 0))
> p1_true-p0_true
[,1]
[1,] 0.09806562
> p1_true/p0_true
[,1]
[1,] 1.983109
log-binomialモデルの関数作成
log-binomialモデルの関数を作成する。log-binomialモデルではリンク関数が対数関数であるため、二項分布のパラメータである
事前分布は多変量正規分布にした。
## logbinomial model
logbinom_mvnorm <- function(beta, prior_mean, prior_sd){
eta <- exp(X%*%beta)
# log prior
lnp <- LaplacesDemon::dmvn(beta, prior_mean, prior_sd, log=TRUE)
# log-binomial model
# assign very low likelihood to p over 1
if(any(eta > 1)) {
lnf <- (-1e6) # very small value
} else {
lnf <- dbinom(Y, size = 1, prob = eta, log=TRUE)
}
return(lnp + sum(lnf))
}
MCMCの実行
チェーン4本でメトロポリス・ヘイスティング法によるMCMCを行う。MH()関数をどのように作成したかなどについてはベイジアンロジスティック回帰からリスク差・リスク比を出すスクラッチ実装 - メトロポリス・ヘイスティング法に記載しているので割愛。
# Run MCMC -----------------------------
# MCMC setting
draw <- 20000
burn <- 5000
num_sample <- draw-burn
C <- 4
# step size of parameters
# to be manually adjusted to have roughly the acceptance rate
# of 20%
step_beta <- 0.2
# prior distributions
prior_mean <- rep(0,P) # mvnorm
prior_sd <- 5 * diag(P) # mvnorm
# run MCMC using parallelizing
cl<- makeCluster(C)
registerDoSNOW(cl)
out_list <- list(NULL)
out_list <-
foreach(i=1:C) %dopar% {
MH(draw=draw,
X = X,
Y = mydata$Y,
burnin = burn,
num_params=P,
step_beta=step_beta,
loglike=logbinom_mvnorm,
prior_mean=prior_mean,
prior_sd=prior_sd,
degree_freedom=NULL
)
}
stopCluster(cl)
# acceptance ratio
for(i in 1:C){
print(out_list[[i]]$acceptance)
}
# Convert results list to an array
out_array <- array(dim = c(num_sample, P, C),
dimnames = list(num_iter = 1:num_sample,
param = paste0('beta', seq(1:P)-1),
chain = 1:C))
for (i in 1:C) {
out_array[,,i] <- out_list[[i]]$result
}
# Transpose the array so that posterior::as_draws_array() works
# The draws_array format assumes the dimensions in the order of (iterations, chains, variables)
out_transposed <- aperm(out_array, c(1, 3, 2))
MCMCの結果
収束診断なども色々と割愛して結果だけ見ていく。
> mcmc_draws <- posterior::as_draws_array(out_transposed)
> posterior::summarise_draws(mcmc_draws)
# A tibble: 5 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 beta0 -2.38 -2.37 0.244 0.244 -2.77 -1.98 1.00 791. 1510.
2 beta1 -0.0727 -0.0731 0.111 0.113 -0.253 0.110 1.00 3995. 5033.
3 beta2 0.576 0.578 0.230 0.228 0.196 0.949 1.00 968. 1481.
4 beta3 0.0228 0.0247 0.219 0.215 -0.348 0.378 1.00 1431. 2227.
5 beta4 0.0273 0.0318 0.222 0.224 -0.341 0.391 1.00 1349. 2377.
exp()した結果の確認。
真のリスク比が約2なので、うまく推定はできていそう。
# Check RR -------
beta2 <- out_transposed[,,3]
colnames(beta2) <- paste0("chain",seq(1:4))
> bind_rows(
+ apply(beta2,2,
+ function(x) {
+ y <- exp(x)
+ list(median = median(y),
+ lwr = quantile(y, prob = 0.025),
+ upr = quantile(y, prob = 0.975))
+ }
+ )
+ )
# A tibble: 4 × 3
median lwr upr
<dbl> <dbl> <dbl>
1 1.78 1.13 2.87
2 1.74 1.10 2.77
3 1.78 1.15 2.72
4 1.83 1.16 2.88
> median(exp(unlist(beta2)))
[1] 1.781635
> quantile(exp(unlist(beta2)), prob = 0.025)
2.5%
1.137019
> quantile(exp(unlist(beta2)), prob = 0.975)
97.5%
2.814798
glm()との比較
ベイズによる方法と同じような推定値を得た。
# GLMでlog-binomial回帰やって確認
# 切片がサンプルデータに入っていたのでそこだけ除いたdfを作る
data = data.frame(X[,-1],Y)
mylogbin <- glm(Y~., data, family = binomial(link = 'log'))
> exp(summary(mylogbin)$coefficients[3,1])
[1] 1.806095
> exp(confint(mylogbin))
Waiting for profiling to be done...
2.5 % 97.5 %
(Intercept) 0.05552346 0.1488225
v1 0.74585671 1.1591524
v2 1.14746200 2.9245436
v3 0.67066731 1.6137626
v4 0.67439436 1.6190663
リスク差の推定
ベイジアンロジスティック回帰からリスク差・リスク比を出すスクラッチ実装 - メトロポリス・ヘイスティング法においてはロジスティック回帰の標準化によってリスク差とリスク比を計算した。
よくよく考えると、ロジスティック回帰でなくとも回帰標準化によって欲しい推定値をえられるのでは？と思い、log-binomialからリスク差を推定してみた。
真のリスク差が0.1なので、悪くはない結果と言えそう。
# RD using regression standardization ------------
X_0 <- X
X_0[,'v2'] <- 0
X_1 <- X
X_1[,'v2'] <- 1
res_list <- list(NULL)
for(i in 1:C) {
res_list[[i]] <- out_list[[i]]$result
}
result <-
reduce(res_list,
.f = ~ matrix(unlist(.x), ncol = P)
)
p0 <- c(exp(X_0%*%t(result)))
p1 <- c(exp(X_1%*%t(result)))
> summary(p1-p0)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.05879 0.05286 0.07379 0.07567 0.09628 0.58959
> quantile(p1-p0, prob = c(0.025, 0.975))
2.5% 97.5%
0.01423722 0.14865164
> bayestestR::hdi(p1-p0)
95% HDI: [0.01, 0.15]
感想
log-binomial回帰を実施したい場合、ベイズの方が
Discussion