data## kibo tch rin q1 q2 qt freq
## 1 0 0 0 1 0 0 5
## 2 0 0 1 1 0 0 2
## 3 0 0 0 0 1 0 3
## 4 0 0 1 0 1 0 4
## 5 0 0 0 0 0 1 6
## 6 0 0 1 0 0 1 4
## 7 0 1 0 1 0 0 1
## 8 0 1 1 1 0 0 7
## 9 0 1 0 0 1 0 3
## 10 0 1 1 0 1 0 9
## 11 0 1 0 0 0 1 2
## 12 0 1 1 0 0 1 1
## 13 1 0 0 1 0 0 5
## 14 1 0 1 1 0 0 0
## 15 1 0 0 0 1 0 8
## 16 1 0 1 0 1 0 5
## 17 1 0 0 0 0 1 5
## 18 1 0 1 0 0 1 2
## 19 1 1 0 1 0 0 2
## 20 1 1 1 1 0 0 4
## 21 1 1 0 0 1 0 4
## 22 1 1 1 0 1 0 4
## 23 1 1 0 0 0 1 4
## 24 1 1 1 0 0 1 6| kibo | tch | rin | q1 | q2 | qt | freq |
|---|---|---|---|---|---|---|
| 動きHz | 動きHz | 動きHz | ||||
| 学級規模 | 教師 | 隣接 | 第1三分位未満 | 第3三分位超 | 三分位範囲内 | 机間指導頻度 |
| (小0/大1) | (A0/B1) | (なし0/あり1) | (未満を1/それ以外0) | (超を1/それ以外0) | (内を1/それ以外0) | (実数) |
学級規模の大(30程度)小(10程度)で，机間指導を受ける児童の直前の動きの相対的な大きさによって，机間指導の受けやすさがどう異なるか。
授業を受けた学級の規模が\(i\)，授業者が\(j\)，隣接した座席に直前に机間指導を受けた児童がいたかが\(k\)，机間指導を受ける直前10秒間の身体の揺れの周波数の相対的な大きさが\(l\)であった児童に対する机間指導件数\(y_{ijkl}\)は
\[ y_{ijkl} \sim Poisson(\mu_{ijkl}) \]
\[ log(\mu_{ij}) = \lambda + \lambda^A_i + \lambda^B_j + \lambda^C_k + \lambda^D_l + \lambda^{AB}_{ij} + \lambda^{BC}_{jk} + \lambda^{AD}_{il} \]
\[ \sum_{l=1}^L \lambda^D_l = \sum_{i=1}^I \lambda^{AD}_{il} = \sum_{l=1}^L \lambda^{AD}_{il} = 0 \]
\[ \lambda , \lambda^A_i , \lambda^B_j , \lambda^C_k , \lambda^D_l , \lambda^{AB}_{ij} , \lambda^{BC}_{jk} , \lambda^{AD}_{il} \sim Uniform (-10, 10) \]
data <- list(Y = data$freq,
kibo = data$kibo,
tch = data$tch,
rin = data$rin,
lo = data$q1,
mid = data$qt,
hi = data$q2,
N = nrow(data)
)
stancode2 <- "
data {
int N;
int Y[N];
int kibo[N];
int tch[N];
int rin[N];
int lo[N];
int mid[N];
int hi[N];
}
parameters {
real mu;
real a;
real b;
real c;
real d1;
real d2;
real ab;
real bc;
real ad1;
real ad2;
}
transformed parameters {
real Lmbd[N];
real LogLmbd[N];
real d3;
real ad3;
d3 = -(d1 + d2);
ad3 = -(ad1 + ad2);
for (i in 1:N)
LogLmbd[i] = mu + a * kibo[i]
+ b * tch[i]
+ c * rin[i]
+ d1 * lo[i]
+ d2 * mid[i]
+ d3 * hi[i]
+ ab * kibo[i] * tch[i]
+ bc * tch[i] * rin[i]
+ ad1 * kibo[i] * lo[i]
+ ad2 * kibo[i] * mid[i]
+ ad3 * kibo[i] * hi[i]
;
for (i in 1:N)
Lmbd[i] = exp(LogLmbd[i]);
}
model {
mu ~ uniform (-10, 10);
a ~ uniform (-10, 10);
b ~ uniform (-10, 10);
c ~ uniform (-10, 10);
d1 ~ uniform (-10, 10);
d2 ~ uniform (-10, 10);
ab ~ uniform (-10, 10);
bc ~ uniform (-10, 10);
ad1 ~ uniform (-10, 10);
ad2 ~ uniform (-10, 10);
for (i in 1:N)
Y[i] ~ poisson (Lmbd[i]);
}
generated quantities{
real p_ad1;
real p_ad2;
real p_ad3;
p_ad1 = step(-ad1);
p_ad2 = step(ad2);
p_ad3 = step(ad3);
}
"
parameters <- c("mu",
"a", "b", "c",
"d1", "d2", "d3",
"ab",
"bc",
"ad1", "ad2", "ad3",
"p_ad1", "p_ad2", "p_ad3"
)library(rstan)
fit3 <- stan(model_code = stancode2,
data = data,
pars = parameters,
iter = 100000,
chains = 4,
thin = 1,
warmup = 50000,
seed = 777
)print(fit3, digits_summary = 2)## Inference for Stan model: bb78def2d3907da58364a20ced4adf02.
## 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% 25% 50% 75% 97.5% n_eff Rhat
## mu 1.59 0.00 0.23 1.11 1.44 1.60 1.75 2.03 102193 1
## a 0.02 0.00 0.29 -0.56 -0.18 0.02 0.22 0.60 110919 1
## b -0.71 0.00 0.38 -1.46 -0.96 -0.71 -0.45 0.02 102925 1
## c -0.65 0.00 0.30 -1.26 -0.85 -0.64 -0.44 -0.06 123517 1
## d1 -0.03 0.00 0.21 -0.46 -0.17 -0.03 0.11 0.38 113406 1
## d2 -0.18 0.00 0.22 -0.63 -0.32 -0.17 -0.03 0.24 112842 1
## d3 0.21 0.00 0.20 -0.19 0.08 0.21 0.35 0.60 243999 1
## ab 0.00 0.00 0.41 -0.81 -0.28 0.00 0.28 0.81 116396 1
## bc 1.32 0.00 0.44 0.48 1.03 1.32 1.61 2.19 113145 1
## ad1 -0.34 0.00 0.31 -0.97 -0.55 -0.34 -0.13 0.27 113214 1
## ad2 0.26 0.00 0.30 -0.33 0.05 0.26 0.46 0.86 110156 1
## ad3 0.08 0.00 0.28 -0.47 -0.11 0.08 0.27 0.64 246558 1
## p_ad1 0.86 0.00 0.35 0.00 1.00 1.00 1.00 1.00 131770 1
## p_ad2 0.80 0.00 0.40 0.00 1.00 1.00 1.00 1.00 132339 1
## p_ad3 0.62 0.00 0.49 0.00 0.00 1.00 1.00 1.00 222930 1
## lp__ 39.04 0.01 2.26 33.76 37.75 39.37 40.69 42.45 81315 1
##
## Samples were drawn using NUTS(diag_e) at Sat Apr 27 19:11:27 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).