放送大学

心理統計法('17) その9

放送大学「心理統計法(‘17)」第11章

授業ではrstanパッケージを使ってますが、MCMCpackパッケージを使ってみます。初学者ですので間違いが多々あると思います。

（Rスクリプトはここ）
早稲田大学文学部文学研究科 豊田研究室

（参考）
ベイジアンMCMCによる統計モデル
Using the ggmcmc package
Bayesian Inference With Stan ~番外編~
Exercise 2: Bayesian A/B testing using MCMC Metropolis-Hastings
政治学方法論 I

1
2
3
library(MCMCpack)
library(knitr)

餡の選好問題

1つの2項分布に関する推測

1
2
3
4
5
6
7
8
9
10
11
12
x <- 305                               #正反応数 
n <- 500                               #データ数 
#
set.seed(1234)
posterior <- MCbinomialbeta(x,n,mc=100000)
#summary(posterior)
df<-unlist(posterior[,1])
ryou<-data.frame(c(round(mean(df),3) , round(sd(df),3) ,
	round(quantile(df,c(0.025,0.05,0.5,0.95,0.975)),3)))
colnames(ryou)<-"p"
rownames(ryou)<-c("EAP","post.sd","0.025","0.05","0.5","0.95","0.975")
kable(t(ryou))
EAP post.sd 0.025 0.05 0.5 0.95 0.975
p 0.61 0.022 0.567 0.574 0.61 0.645 0.652

予測分布

1
2
3
4
5
6
7
xaste<-rbinom(length(df),n, df)
#
ryou<-data.frame(c(round(mean(xaste),0) , round(sd(xaste),1) ,
	round(quantile(xaste,c(0.025,0.05,0.5,0.95,0.975)),0)))
colnames(ryou)<-"x*"
rownames(ryou)<-c("EAP","post.sd","0.025","0.05","0.5","0.95","0.975")
kable(t(ryou))
EAP post.sd 0.025 0.05 0.5 0.95 0.975
x* 305 15.4 274 279 305 330 335

オッズ

1
2
3
4
5
6
7
odds<-df/(1-df)
#
ryou<-data.frame(c(round(mean(odds),2) , round(sd(odds),2) ,
	round(quantile(odds,c(0.025,0.05,0.5,0.95,0.975)),2)))
colnames(ryou)<-"odds"
rownames(ryou)<-c("EAP","post.sd","0.025","0.05","0.5","0.95","0.975")
kable(t(ryou))
EAP post.sd 0.025 0.05 0.5 0.95 0.975
odds 1.57 0.14 1.31 1.35 1.56 1.82 1.87
つぶ餡が好きな人は0.65より大きい
1
2
p<-table(df > 0.65)
p[2]/length(df)

TRUE
0.03086

もう一度同じ調査をするとつぶ餡が好きな人は330人より多い
1
2
p<-table(xaste > 330)
p[2]/length(xaste)

TRUE
0.04706

つぶ餡好きはこし餡好きの1.6倍より多い
1
2
p<-table(odds > 1.6)
p[2]/length(odds)

TRUE
0.39434

2つの2項分布に関する推測

ブランド認知問題1

1
2
3
4
5
6
7
8
9
10
11
x <- c(85,31)                          #正反応数 
n <- c(123,121)                        #データ数 
#
post1 <- MCbinomialbeta(x[1],n[1],mc=100000)
post2 <- MCbinomialbeta(x[2],n[2],mc=100000)
#
df<-data.frame(unlist(post1[,1]),unlist(post2[,1]) )
colnames(df)<-c("p1","p2")
ryou<-rbind(EAP=apply(df, 2, mean),post.sd=apply(df, 2, sd),
apply(df, 2, quantile, probs=c(0.025,0.05, 0.5,0.95, 0.975)))
kable(t(round(ryou,3)))
EAP post.sd 0.025 0.05 0.5 0.95 0.975
p1 0.688 0.041 0.604 0.619 0.689 0.754 0.766
p2 0.260 0.040 0.186 0.197 0.259 0.327 0.341

比率の差 , 比率の比 , オッズ・オッズ比

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#### 比率の差
#
hisa<-df[,1] -  df[,2]
#
#### 比率の比
#
hihi<-df[,1] /  df[,2]
#
#### オッズ・オッズ比
#
oddss1<-df[,1]/(1 - df[,1])
#
oddss2<-df[,2]/(1 - df[,2])
#
oddshi<- oddss1 / oddss2 # (df[,1]/(1 - df[,1])) /  (df[,2]/(1 - df[,2]))
#
seisei<-cbind(hisa,hihi,oddss1,oddss2,oddshi)
#
ryou<-rbind(EAP=apply(seisei, 2, mean),post.sd=apply(seisei, 2, sd),
apply(seisei, 2, quantile, probs=c(0.025,0.05, 0.5,0.95, 0.975)))
colnames(ryou)<-c("比率の差","比率の比","oddss1","oddss2","オッズ比")
kable(t(round(ryou,3)))
EAP post.sd 0.025 0.05 0.5 0.95 0.975
比率の差 0.428 0.057 0.313 0.332 0.429 0.520 0.537
比率の比 2.710 0.460 1.962 2.053 2.658 3.544 3.760
oddss1 2.265 0.447 1.527 1.624 2.219 3.070 3.270
oddss2 0.355 0.074 0.229 0.246 0.349 0.487 0.517
オッズ比 6.654 1.949 3.688 4.029 6.370 10.259 11.260
女性の認知率 ー 男性の認知率 ＞ ０．３
1
2
ninchisa<-table(df[,1]-df[,2] > 0.3)
ninchisa[2]/nrow(df)

TRUE
0.98541

女性の認知率 / 男性の認知率 ＞ 3
1
2
ninchihi<-table(df[,1] / df[,2] > 3)
ninchihi[2]/nrow(df)

TRUE
0.23597

オッズ比　＞　８
1
table(oddshi > 8)[2] / length(oddshi)

TRUE
0.21336

g個の2項分布に関する推測

お年玉問題

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
x <- c(42,31,29,20)                     #正反応数 
n <- c(51,49,50,48)                     #データ数 
#  g <- length(x)                          #群の数
#
#
post1 <- MCbinomialbeta(x[1],n[1],mc=100000)
post2 <- MCbinomialbeta(x[2],n[2],mc=100000)
post3 <- MCbinomialbeta(x[3],n[3],mc=100000)
post4 <- MCbinomialbeta(x[4],n[4],mc=100000)
#
df<-data.frame(p1=unlist(post1[,1]),p2=unlist(post2[,1]) ,
	p3=unlist(post3[,1]),p4=unlist(post4[,1]))
colnames(df)<-c("p1","p2","p3","p4")
ryou<-rbind(EAP=apply(df, 2, mean),post.sd=apply(df, 2, sd),
apply(df, 2, quantile, probs=c(0.025,0.05, 0.5,0.95, 0.975)))
kable(t(round(ryou,3)))
EAP post.sd 0.025 0.05 0.5 0.95 0.975
p1 0.811 0.053 0.697 0.718 0.815 0.892 0.904
p2 0.627 0.067 0.493 0.515 0.629 0.735 0.753
p3 0.576 0.068 0.442 0.463 0.577 0.686 0.706
p4 0.420 0.069 0.289 0.308 0.419 0.536 0.558

pjとpj’との比較

1
2
3
4
5
6
7
8
9
hikaku<-matrix(rep(NA,16),ncol=4)
rownames(hikaku)<-colnames(hikaku)<-c("高校生(p1)","大学生[前](p2)","大学生[後](p3)","社会人(p4)")
for (i in 1:4){
	for (j in 1:4){
		hikaku[i,j]<-round(sum((df[,i]-df[,j]>0)=="TRUE")/nrow(df),3)
	}
}
#
kable(hikaku)
高校生(p1) 大学生前 大学生後 社会人(p4)
高校生(p1) 0.000 0.984 0.996 1.000
大学生前 0.016 0.000 0.701 0.983
大学生後 0.004 0.299 0.000 0.946
社会人(p4) 0.000 0.017 0.054 0.000

連言命題が正しい確率

u_p4<p3 u_p3<p2 u_p2<p1

1
sum((df[,3]-df[,4]>0)=="TRUE" & (df[,2]-df[,3]>0)=="TRUE" & (df[,1]-df[,2]>0)=="TRUE")/nrow(df)

0.63403

u_p4<p3 u_p4<p2 u_p2<p1 * u_p3<p1

1
sum((df[,3]-df[,4]>0)=="TRUE" & (df[,2]-df[,4]>0)=="TRUE" & (df[,1]-df[,2]>0)=="TRUE" & (df[,1]-df[,3]>0)=="TRUE")/nrow(df)

0.91629

u_p2<p1 u_p3<p1 u_p4<p1

1
sum((df[,1]-df[,2]>0)=="TRUE" & (df[,1]-df[,3]>0)=="TRUE" & (df[,1]-df[,4]>0)=="TRUE")/nrow(df)

0.98091

u_p4<p1 u_p4<p2 u_p4<p3

1
sum((df[,1]-df[,4]>0)=="TRUE" & (df[,2]-df[,4]>0)=="TRUE" & (df[,3]-df[,4]>0)=="TRUE")/nrow(df)

0.93452