２元配置分散分析
1 2元配置分散分析
1.1 青木
青木繁伸先生の例題で実行してみる。
source("data/anovakun_485.txt")1.1.1 各水準の繰返し数が 1 である場合.
http://aoki2.si.gunma-u.ac.jp/lecture/TwoWayANOVA/TwoWay1.html
1.2 Jin
a1<- c(9,6,6,11)
a2<- c(4,1,4,7)
a3<- c(5,2,2,3)
bunsan3<-data.frame(A=factor(c(rep("a1",4),rep("a2",4),rep("a3",4))),B=factor(rep(c(rep("b1",1),rep("b2",1),rep("b3",1),rep("b4",1)),3)),y= c(a1,a2,a3))
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 1.0.1
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.3.0 ✔ stringr 1.5.0
## ✔ readr 2.1.4 ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()#wide型への変換
bunsan4 <-spread(bunsan3,key=B,value=y)
bunsan3| A | B | y |
|---|---|---|
| a1 | b1 | 9 |
| a1 | b2 | 6 |
| a1 | b3 | 6 |
| a1 | b4 | 11 |
| a2 | b1 | 4 |
| a2 | b2 | 1 |
| a2 | b3 | 4 |
| a2 | b4 | 7 |
| a3 | b1 | 5 |
| a3 | b2 | 2 |
| a3 | b3 | 2 |
| a3 | b4 | 3 |
bunsan4| A | b1 | b2 | b3 | b4 |
|---|---|---|---|---|
| a1 | 9 | 6 | 6 | 11 |
| a2 | 4 | 1 | 4 | 7 |
| a3 | 5 | 2 | 2 | 3 |
#交互作用なし
summary(aov(y~A+B,data=bunsan3))## Df Sum Sq Mean Sq F value Pr(>F)
## A 2 56 28 14 0.0055 **
## B 3 30 10 5 0.0452 *
## Residuals 6 12 2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1interaction.plot(x.factor=bunsan3$B,trace.factor=bunsan3$A,response=bunsan3$y, type="b", pch=c(1,2))
市原p175の図を参照すると、interaction.plotが適度に交差している場合、行間・列間に差があるとは言えない。
1.3 参考 anovakun.
1.3.1 被験者間計画.
要因Bによる効果を検出したいとする。
要因Aは個人差であるが、個人差を残差に含めて分散分析を行うと、
anovakun(bunsan3[-1], "As", 4, auto = TRUE, holm = TRUE, eta = TRUE)##
## [ As-Type Design ]
##
## This output was generated by anovakun 4.8.5 under R version 4.2.2.
## It was executed on Sat Jan 27 17:34:03 2024.
##
##
## << DESCRIPTIVE STATISTICS >>
##
## --------------------------
## A n Mean S.D.
## --------------------------
## a1 3 6.0000 2.6458
## a2 3 3.0000 2.6458
## a3 3 4.0000 2.0000
## a4 3 7.0000 4.0000
## --------------------------
##
##
## << ANOVA TABLE >>
##
## -----------------------------------------------------------
## Source SS df MS F-ratio p-value eta^2
## -----------------------------------------------------------
## A 30.0000 3 10.0000 1.1765 0.3777 ns 0.3061
## Error 68.0000 8 8.5000
## -----------------------------------------------------------
## Total 98.0000 11 8.9091
## +p < .10, *p < .05, **p < .01, ***p < .001
##
##
## output is over --------------------///残差分散が大きいため検出力が低い。
1.3.2 被験者内計画.
被験者内計画では、残差から個人差を除いて検定するので検出力が上がる。いわゆる乱塊法である。
anovakun(bunsan4[-1], "sA", 4, eta = TRUE)##
## [ sA-Type Design ]
##
## This output was generated by anovakun 4.8.5 under R version 4.2.2.
## It was executed on Sat Jan 27 17:34:04 2024.
##
##
## << DESCRIPTIVE STATISTICS >>
##
## --------------------------
## A n Mean S.D.
## --------------------------
## a1 3 6.0000 2.6458
## a2 3 3.0000 2.6458
## a3 3 4.0000 2.0000
## a4 3 7.0000 4.0000
## --------------------------
##
##
## << SPHERICITY INDICES >>
##
## == Mendoza's Multisample Sphericity Test and Epsilons ==
##
## -------------------------------------------------------------------------
## Effect Lambda approx.Chi df p LB GG HF CM
## -------------------------------------------------------------------------
## A 0.0000 53.5173 5 0.0000 *** 0.3333 0.5000 1.6667 0.3333
## -------------------------------------------------------------------------
## LB = lower.bound, GG = Greenhouse-Geisser
## HF = Huynh-Feldt-Lecoutre, CM = Chi-Muller
##
##
## << ANOVA TABLE >>
##
## -----------------------------------------------------------
## Source SS df MS F-ratio p-value eta^2
## -----------------------------------------------------------
## s 56.0000 2 28.0000
## -----------------------------------------------------------
## A 30.0000 3 10.0000 5.0000 0.0452 * 0.3061
## s x A 12.0000 6 2.0000
## -----------------------------------------------------------
## Total 98.0000 11 8.9091
## +p < .10, *p < .05, **p < .01, ***p < .001
##
##
## << POST ANALYSES >>
##
## < MULTIPLE COMPARISON for "A" >
##
## == Shaffer's Modified Sequentially Rejective Bonferroni Procedure ==
## == The factor < A > is analysed as dependent means. ==
## == Alpha level is 0.05. ==
##
## --------------------------
## A n Mean S.D.
## --------------------------
## a1 3 6.0000 2.6458
## a2 3 3.0000 2.6458
## a3 3 4.0000 2.0000
## a4 3 7.0000 4.0000
## --------------------------
##
## ----------------------------------------------------------
## Pair Diff t-value df p adj.p
## ----------------------------------------------------------
## a1-a2 3.0000 Inf 2 0.0000 0.0000 a1 > a2 *
## a2-a4 -4.0000 2.6186 2 0.1201 0.3604 a2 = a4
## a3-a4 -3.0000 2.5981 2 0.1217 0.3651 a3 = a4
## a1-a3 2.0000 2.0000 2 0.1835 0.5505 a1 = a3
## a2-a3 -1.0000 1.0000 2 0.4226 0.8453 a2 = a3
## a1-a4 -1.0000 0.6547 2 0.5799 0.8453 a1 = a4
## ----------------------------------------------------------
##
##
## output is over --------------------///個人による変動が除かれるため、要因Bの効果は有意となっている。
2 Friedman検定（p184）
データを順位化し、カイ二乗統計量に帰着するのがFriedman検定である。
2.1 乱数シミュレーション1.
2.2 乱数シミュレーション2.
2.2.1 データ.
library(tidyr)
data <- crossing(x = 1:4, y = 1:4);#data
data| x | y |
|---|---|
| 1 | 1 |
| 1 | 2 |
| 1 | 3 |
| 1 | 4 |
| 2 | 1 |
| 2 | 2 |
| 2 | 3 |
| 2 | 4 |
| 3 | 1 |
| 3 | 2 |
| 3 | 3 |
| 3 | 4 |
| 4 | 1 |
| 4 | 2 |
| 4 | 3 |
| 4 | 4 |
library(tidyr)
library(gapminder)
options(rgl.printRglwidget = TRUE)
#WebGLに変換してRStudio, VSCode, ウェブブラウザなど出力
library(rgl)
x <- data$x
y <- data$y
#z <- 0.2*x - 0.3*y+ rnorm(16, sd = 0.1)#誤差分散 0.1^2
z <- 0*x - 0.3*y+ rnorm(16, sd = 0.1)#誤差分散 0.1^2
#z <- 0*x - 0*y+ rnorm(16, sd = 0.1)#誤差分散 0.1^2#plot3d(x = x, y = y, z = z,type = "s", col = "blue", size = 1,zlim=c(-0.6,0.6))
dat <- data.frame(x=x,y=y,z=z)
dat| x | y | z |
|---|---|---|
| 1 | 1 | -0.2595007 |
| 1 | 2 | -0.6179110 |
| 1 | 3 | -0.7260247 |
| 1 | 4 | -1.2874619 |
| 2 | 1 | -0.5933815 |
| 2 | 2 | -0.4906335 |
| 2 | 3 | -0.8298534 |
| 2 | 4 | -1.3966036 |
| 3 | 1 | -0.3057194 |
| 3 | 2 | -0.5999264 |
| 3 | 3 | -0.6869141 |
| 3 | 4 | -1.2929743 |
| 4 | 1 | -0.4896250 |
| 4 | 2 | -0.6732608 |
| 4 | 3 | -0.6914906 |
| 4 | 4 | -1.2350756 |
dat1 <- spread(dat,key=y,value=z)source("data/anovakun_485.txt")
#(dat1[,-1])#被験内間計画のデータ形式 sA
anovakun((dat1[,-1]),"sA",4,holm=TRUE,eta=TRUE)##
## [ sA-Type Design ]
##
## This output was generated by anovakun 4.8.5 under R version 4.2.2.
## It was executed on Sat Jan 27 17:34:04 2024.
##
##
## << DESCRIPTIVE STATISTICS >>
##
## ---------------------------
## A n Mean S.D.
## ---------------------------
## a1 4 -0.4121 0.1565
## a2 4 -0.5954 0.0765
## a3 4 -0.7336 0.0665
## a4 4 -1.3030 0.0676
## ---------------------------
##
##
## << SPHERICITY INDICES >>
##
## == Mendoza's Multisample Sphericity Test and Epsilons ==
##
## -------------------------------------------------------------------------
## Effect Lambda approx.Chi df p LB GG HF CM
## -------------------------------------------------------------------------
## A 0.0085 5.4701 5 0.4278 ns 0.3333 0.5665 1.2299 0.3333
## -------------------------------------------------------------------------
## LB = lower.bound, GG = Greenhouse-Geisser
## HF = Huynh-Feldt-Lecoutre, CM = Chi-Muller
##
##
## << ANOVA TABLE >>
##
## ---------------------------------------------------------
## Source SS df MS F-ratio p-value eta^2
## ---------------------------------------------------------
## s 0.0304 3 0.0101
## ---------------------------------------------------------
## A 1.7749 3 0.5916 60.7643 0.0000 *** 0.9376
## s x A 0.0876 9 0.0097
## ---------------------------------------------------------
## Total 1.8929 15 0.1262
## +p < .10, *p < .05, **p < .01, ***p < .001
##
##
## << POST ANALYSES >>
##
## < MULTIPLE COMPARISON for "A" >
##
## == Holm's Sequentially Rejective Bonferroni Procedure ==
## == The factor < A > is analysed as dependent means. ==
## == Alpha level is 0.05. ==
##
## ---------------------------
## A n Mean S.D.
## ---------------------------
## a1 4 -0.4121 0.1565
## a2 4 -0.5954 0.0765
## a3 4 -0.7336 0.0665
## a4 4 -1.3030 0.0676
## ---------------------------
##
## ---------------------------------------------------------
## Pair Diff t-value df p adj.p
## ---------------------------------------------------------
## a3-a4 0.5695 43.2451 3 0.0000 0.0002 a3 > a4 *
## a1-a4 0.8910 12.9375 3 0.0010 0.0050 a1 > a4 *
## a2-a4 0.7076 9.8233 3 0.0022 0.0090 a2 > a4 *
## a1-a3 0.3215 5.1852 3 0.0139 0.0418 a1 > a3 *
## a2-a3 0.1381 1.9813 3 0.1419 0.2837 a2 = a3
## a1-a2 0.1834 1.7982 3 0.1700 0.2837 a1 = a2
## ---------------------------------------------------------
##
##
## output is over --------------------///