| クラスサイズ | 授業者 | 児童数 | 机間指導を受けた児童数 |
|---|---|---|---|
| 小 | A | 12 | 11 |
| 大 | A | 30 | 15 |
| 小 | B | 10 | 9 |
| 大 | B | 33 | 19 |
data <- matrix(c(1, 2, 1, 2,
1, 1, 2, 2,
12, 30, 10, 33,
11, 15, 9, 19),
nrow = 4, ncol = 4)
data <- data.frame(data)
colnames(data) <- c("cs", "tch", "n.pupil", "n.kikan")
data## cs tch n.pupil n.kikan
## 1 1 1 12 11
## 2 2 1 30 15
## 3 1 2 10 9
## 4 2 2 33 194時間（10人程度2時間・30人程度2時間，授業者A 2時間・授業者B 2時間）における机間指導を受けた児童数\(x_i (i=1,2,3,4)\)が，各授業における児童数\(n_i\)，母比率\(p_i\)の2項分布にしたがうと仮定したとき，データ\(\boldsymbol{x} = (x_1, x_2, x_3, x_4)\)，母数\(\boldsymbol{\theta} = (p_1, p_2, p_3, p_4)\)の尤度は，4個の2項分布の積
\[ f(\boldsymbol{x}|\boldsymbol{\theta}) = f(x_1, x_2, x_3, x_4|p_1, p_2, p_3, p_4) = f(x_1|p_1) \times f(x_2|p_2) \times f(x_3|p_3) \times f(x_4|p_4) \]
となる。事前分布\(f(p_1), f(p_2), f(p_3),f(p_4)\)は確率の定義域に対する区間\([0,1]\)の一様分布を仮定し，同時事前分布を
\[ f(\boldsymbol{\theta}) = f(p_1) \times f(p_2) \times f(p_3) \times f(p_4) \]
と仮定し，事後分布を
\[ f(\boldsymbol{\theta}|\boldsymbol{x}) \varpropto f(\boldsymbol{x}|\boldsymbol{\theta}) f(\boldsymbol{\theta}) \]
と導き，MCMCで母数の事後分布を近似させることとした。
また，10人程度クラス（データの1, 3行目）の方が30人程度クラス（2, 4行目）の方が机間指導を受ける児童の割合が高いという仮説が真のときに1，偽のときに0をとる生成量
\[ u^{(t)}_{p1>p2} \times u^{(t)}_{p1>p4} \times u^{(t)}_{p3>p2} \times u^{(t)}_{p3>p4} \]
と，授業者A（データの1, 2行目）の方が授業者B（3, 4行目）の方が机間指導を受ける児童の割合が高いという仮説が真のときに1，偽のときに0をとる生成量
\[ u^{(t)}_{p1>p3} \times u^{(t)}_{p1>p4} \times u^{(t)}_{p2>p3} \times u^{(t)}_{p2>p4} \]
を定義した。
# g X 2 のクロス表の推測
# 母比率の推測と比較
## 「はじめての統計データ分析」第6章 6.4.3節 p. 149 と
## http://www.asakura.co.jp/G_27_2.php?id=232の
## scrA/part6.Rの l.31-38
## scrA/Bi03.R，stan/Bi03.stan
## scrA/myfunc/E1Ind.Rの l.57-71 を参考にした。
################################################################
# データ
################################################################
data <- matrix(c(1, 2, 1, 2,
1, 1, 2, 2,
12, 30, 10, 33,
11, 15, 9, 19),
nrow = 4, ncol = 4)
data <- data.frame(data)
colnames(data) <- c("cs", "tch", "n.pupil", "n.kikan")
x <- data$n.kikan #机間指導児童数
n <- data$n.pupil #児童数
g <- length(x)
dat <- list(g = g, x = x, n = n)
################################################################
# stanのコード
################################################################
stancode <- "
//対応のないg群の比率に関する統計的推測
data {
int<lower=0> g; //群の数
int<lower=0> x[g]; //正反応数
int<lower=0> n[g]; //データ数
}
parameters {
real<lower=0,upper=1> p[g]; //母比率
}
transformed parameters {
}
model {
for (i in 1:g){
x[i] ~ binomial(n[i],p[i]);}
}
generated quantities{
int<lower=0> xaste[g];
real log_lik;
real<lower=0,upper=1> U2[g,g]; //2水準比較
for (i in 1:g){ U2[i,i] =0;}
for (i in 1:(g-1)){
for (j in (i+1):g){
U2[i,j] = p[i]-p[j]>0 ? 1 : 0;
U2[j,i] = !(U2[i,j]); } }
log_lik = 0.0;
for (i in 1:g){
xaste[i] = binomial_rng(n[i],p[i]); //予測分布
log_lik = log_lik + binomial_lpmf(x[i]|n[i],p[i]);} //対数確率
}
"
################################################################
# 推定するパラメータ
################################################################
parameters <-c ("p", #"xaste",
"U2", "log_lik"
)
################################################################
# stanのコードを読み込んで推定する自作関数
################################################################
Bi03 <- function(x,
n,
prob=c(0.025, 0.05, 0.5, 0.95, 0.975),
see = 777,
cha = 4,
war = 50000,
ite = 100000,
fi = NA
)
{
library(rstan)
scr <- stancode # Stanファイル名
par <- parameters # sampling変数
fit <- stan(model_code = scr,
data = dat,
pars = parameters,
seed = see,
chains = cha,
warmup = war,
iter = ite,
fit = fi
)
ext <- extract(fit, par);
out <- list(fit = fit,
par = par,
prob=prob,
p=ext$p,
# xaste=ext$xaste,
U2=ext$U2,
log_lik=ext$log_lik
);
class(out) <- "Bi03"
return(invisible(out))
}
################################################################
# 結果の出力に関する自作関数
################################################################
print.Bi03 <- function(x, digits=3)
{
log_lik <- x$log_lik;
p <- x$p;
# xaste <- x$xaste;
U2 <- x$U2;
waic <- (-2) * (log(mean(exp(log_lik)))) + 2*(var(log_lik))
print(x$fit, pars = x$par, digits_summary = digits, probs = x$prob)
print(paste("waic=",round(waic,digits),sep=""))
out <- list(waic=waic)
return(invisible(out))
}
################################################################
# 推定の実施と結果の表示
################################################################
Bi03out <- Bi03(x, n, fi = NA) #g個の2項分布に関する推測
print(Bi03out,2)
################################################################
# 連言命題が正しい確率を求める関数
################################################################
printIJ<-function(x, digits=3, IJ)
{
pro <- 1
for (i in 1:nrow(IJ)){
pro <- pro*x[,IJ[i,1],IJ[i,2]]
}
print(round(mean(pro),digits))
}
# 2クラスあるクラスサイズの大小で机間指導を受ける児童の割合が違う
## 小規模の方が割合が高い
printIJ(Bi03out$U2,3,IJ=rbind(c(1,2),c(3,4),c(1,4),c(3,2)))
# 2人の指導者が実施した2クラスの授業で机間指導を受ける児童の割合が違う
## 1番目の指導者の方が割合が高い
printIJ(Bi03out$U2,3,IJ=rbind(c(1,3),c(1,4),c(2,4),c(2,3)))# 推定結果
print(Bi03out,2)## Inference for Stan model: f1051ddcd32148ff12c0c193c30a7630.
## 4 chains, each with iter=1e+05; warmup=50000; thin=1;
## post-warmup draws per chain=50000, total post-warmup draws=2e+05.
##
## mean se_mean sd 2.5% 5% 50% 95% 97.5% n_eff Rhat
## p[1] 0.86 0 0.09 0.64 0.68 0.87 0.97 0.98 234385 1
## p[2] 0.50 0 0.09 0.33 0.36 0.50 0.64 0.67 214249 1
## p[3] 0.83 0 0.10 0.59 0.64 0.85 0.97 0.98 237582 1
## p[4] 0.57 0 0.08 0.41 0.43 0.57 0.70 0.73 215090 1
## U2[1,1] 0.00 NaN 0.00 0.00 0.00 0.00 0.00 0.00 NaN NaN
## U2[1,2] 0.99 0 0.08 1.00 1.00 1.00 1.00 1.00 177457 1
## U2[1,3] 0.56 0 0.50 0.00 0.00 1.00 1.00 1.00 205978 1
## U2[1,4] 0.98 0 0.13 1.00 1.00 1.00 1.00 1.00 170895 1
## U2[2,1] 0.01 0 0.08 0.00 0.00 0.00 0.00 0.00 177457 1
## U2[2,2] 0.00 NaN 0.00 0.00 0.00 0.00 0.00 0.00 NaN NaN
## U2[2,3] 0.01 0 0.12 0.00 0.00 0.00 0.00 0.00 174895 1
## U2[2,4] 0.28 0 0.45 0.00 0.00 0.00 1.00 1.00 179093 1
## U2[3,1] 0.44 0 0.50 0.00 0.00 0.00 1.00 1.00 205978 1
## U2[3,2] 0.99 0 0.12 1.00 1.00 1.00 1.00 1.00 174895 1
## U2[3,3] 0.00 NaN 0.00 0.00 0.00 0.00 0.00 0.00 NaN NaN
## U2[3,4] 0.96 0 0.18 0.00 1.00 1.00 1.00 1.00 171941 1
## U2[4,1] 0.02 0 0.13 0.00 0.00 0.00 0.00 0.00 170895 1
## U2[4,2] 0.72 0 0.45 0.00 0.00 1.00 1.00 1.00 179093 1
## U2[4,3] 0.04 0 0.18 0.00 0.00 0.00 0.00 1.00 171941 1
## U2[4,4] 0.00 NaN 0.00 0.00 0.00 0.00 0.00 0.00 NaN NaN
## log_lik -7.75 0 1.37 -11.21 -10.41 -7.44 -6.15 -6.05 128526 1
##
## Samples were drawn using NUTS(diag_e) at Mon Oct 7 16:03:11 2019.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
## [1] "waic=18.09"# 2クラスあるクラスサイズの大小で机間指導を受ける児童の割合が違う
## 小規模の方が割合が高い
printIJ(Bi03out$U2,3,IJ=rbind(c(1,2),c(3,4),c(1,4),c(3,2)))## [1] 0.94# 2人の指導者が実施した2クラスの授業で机間指導を受ける児童の割合が違う
## 1番目の指導者の方が割合が高い
printIJ(Bi03out$U2,3,IJ=rbind(c(1,3),c(1,4),c(2,4),c(2,3)))## [1] 0.009