rstanarmを使ってロジスティクス回帰(事前情報とサンプルサイズ別の結果)パラメータ1つ

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<-rbinom(nnn,1,5/10)
	pp<-0.1+v1_vec*0.15
	outcome_vec<-apply(matrix(pp),1,function(x) rbinom(1,1,x))
	df1<-data.frame(v1_vec,outcome_vec)
 #df2  大サンプル
	nnn<-1600
	v1_vec<-rbinom(nnn,1,5/10)
	pp<-0.1+v1_vec*0.15
	outcome_vec<-apply(matrix(pp),1,function(x) rbinom(1,1,x))
	df2<-data.frame(v1_vec,outcome_vec)


 #少サンプル 、無情報事前分布
df<-df1
names(df)<-c("v1","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 = 4, 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 , data = df, family = "binomial")
summary(logistic_model)




 #大サンプル 、無情報事前分布
df<-df2
names(df)<-c("v1","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 = 4, 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, data = df, family = "binomial")
summary(logistic_model)

小サンプル 無情報





小サンプル 情報あり






大サンプル 無情報




大サンプル 情報あり





この記事が気に入ったらサポートをしてみませんか?