rstanarmを使ってロジスティクス回帰（事前情報とサンプルサイズ別の結果）

2024年4月18日 21:40

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)