rstanarmを使ってロジスティクス回帰(事前情報とサンプルサイズ別の結果)
library(tidyverse)
library(caret)
library(GGally)
library(ggplot2)
library(corrplot)
library(bayesplot)
theme_set(bayesplot::theme_default(base_family = "sans"))
library(rstanarm)
options(mc.cores = 1)
library(loo)
library(projpred)
SEED=14124869
#df1 少サンプル
nnn<-40
v1_vec<-round(rnorm(nnn,60,10))
v2_vec<-rbinom(nnn,1,5/10)
v3_vec<-rbinom(nnn,1,5/10)
v4_vec<-rbinom(nnn,1,5/10)
pp<-0.1+v2_vec*0.15
outcome_vec<-apply(matrix(pp),1,function(x) rbinom(1,1,x))
df1<-data.frame(v1_vec,v2_vec,v3_vec,v4_vec,outcome_vec)
#df2 大サンプル
nnn<-1600
v1_vec<-round(rnorm(nnn,60,10))
v2_vec<-rbinom(nnn,1,5/10)
v3_vec<-rbinom(nnn,1,5/10)
v4_vec<-rbinom(nnn,1,5/10)
pp<-0.1+v2_vec*0.15
outcome_vec<-apply(matrix(pp),1,function(x) rbinom(1,1,x))
df2<-data.frame(v1_vec,v2_vec,v3_vec,v4_vec,outcome_vec)
#少サンプル 、無情報事前分布
df<-df1
names(df)<-c("v1","v2","v3","v4","outcome")
df[1]<-scale(df[1])
n=dim(df)[1]
p=dim(df)[2]
corrplot(cor(df[,c(1:4)]))
df$outcome<-factor(df$outcome)
x<-model.matrix(outcome ~. -1, data = df)
y<-df$outcome
(reg_formula <- formula(paste("outcome ~", paste(names(df)[1:(p-1)], collapse = " + "))))
#無情報
t_prior <- normal(location = 0, scale = 1)
post1 <- stan_glm(reg_formula, data = df,
family = binomial(link = "logit"),
prior = t_prior, prior_intercept = t_prior, QR=TRUE,
seed = SEED, refresh = 0)
pplot <- plot(post1, "areas", prob = 0.95, prob_outer = 1)
pplot + geom_vline(xintercept = 0)
round(coef(post1), 2)
round(posterior_interval(post1, prob = 0.9), 2)
#有情報
t_prior <- normal(location = c(0,4,4,0), scale = 1)
t_prior_intercept <- normal(location = 0, scale = 1)
post1 <- stan_glm(reg_formula, data = df,
family = binomial(link = "logit"),
prior = t_prior, prior_intercept = t_prior_intercept, QR=TRUE,
seed = SEED, refresh = 0)
pplot <- plot(post1, "areas", prob = 0.95, prob_outer = 1)
pplot + geom_vline(xintercept = 0)
round(coef(post1), 2)
round(posterior_interval(post1, prob = 0.9), 2)
logistic_model <- glm(outcome ~v1 + v2 +v3 +v4, data = df, family = "binomial")
summary(logistic_model)
#大サンプル 、無情報事前分布
df<-df2
names(df)<-c("v1","v2","v3","v4","outcome")
df[1]<-scale(df[1])
n=dim(df)[1]
p=dim(df)[2]
corrplot(cor(df[,c(1:4)]))
df$outcome<-factor(df$outcome)
x<-model.matrix(outcome ~. -1, data = df)
y<-df$outcome
(reg_formula <- formula(paste("outcome ~", paste(names(df)[1:(p-1)], collapse = " + "))))
#無情報
t_prior <- normal(location = 0, scale = 1)
post1 <- stan_glm(reg_formula, data = df,
family = binomial(link = "logit"),
prior = t_prior, prior_intercept = t_prior, QR=TRUE,
seed = SEED, refresh = 0)
pplot <- plot(post1, "areas", prob = 0.95, prob_outer = 1)
pplot + geom_vline(xintercept = 0)
round(coef(post1), 2)
round(posterior_interval(post1, prob = 0.9), 2)
#有情報
t_prior <- normal(location = c(0,4,4,0), scale = 1)
t_prior_intercept <- normal(location = 0, scale = 1)
post1 <- stan_glm(reg_formula, data = df,
family = binomial(link = "logit"),
prior = t_prior, prior_intercept = t_prior_intercept, QR=TRUE,
seed = SEED, refresh = 0)
pplot <- plot(post1, "areas", prob = 0.95, prob_outer = 1)
pplot + geom_vline(xintercept = 0)
round(coef(post1), 2)
round(posterior_interval(post1, prob = 0.9), 2)
logistic_model <- glm(outcome ~v1 + v2 +v3 +v4, data = df, family = "binomial")
summary(logistic_model)
小サンプル、無情報
![](https://assets.st-note.com/img/1713443043346-NMMoR2eIVH.png?width=1200)
![](https://assets.st-note.com/img/1713443156969-DhuOJMNZQF.png?width=1200)
小サンプル 事前情報あり
![](https://assets.st-note.com/img/1713443080226-8m0GBpx7yP.png?width=1200)
![](https://assets.st-note.com/img/1713443178792-CrJ96FQxvb.png?width=1200)
![](https://assets.st-note.com/img/1713443243274-OHDGIjh1Pc.png?width=1200)
大サンプル 無情報
![](https://assets.st-note.com/img/1713443269020-Rinwkaa0as.png?width=1200)
![](https://assets.st-note.com/img/1713443279620-RX9N8T2aET.png?width=1200)
大サンプル 有情報
![](https://assets.st-note.com/img/1713443321808-aFOAji4HFf.png?width=1200)
![](https://assets.st-note.com/img/1713443331760-RY35ZOk55I.png?width=1200)
![](https://assets.st-note.com/img/1713443349388-rqfFdfHOtI.png?width=1200)
この記事が気に入ったらサポートをしてみませんか?