5章_分散分析・多重比較
1 繰り返しのない１元配置分散分析
dat <- read.csv("data/data_ch5-1.csv", header = TRUE)
head(dat)| Student | class | score |
|---|---|---|
| S001 | A | 76 |
| S002 | A | 54 |
| S003 | A | 62 |
| S004 | A | 46 |
| S005 | A | 53 |
| S006 | A | 64 |
str(dat)## 'data.frame': 87 obs. of 3 variables:
## $ Student: chr "S001" "S002" "S003" "S004" ...
## $ class : chr "A" "A" "A" "A" ...
## $ score : int 76 54 62 46 53 64 42 96 87 92 ...table(dat$class)##
## A B C
## 29 29 292 等分散の検定
leveneTest(dat$score,dat$class,center=mean)## Warning in leveneTest.default(dat$score, dat$class, center = mean): dat$class
## coerced to factor.| Df | F value | Pr(>F) | |
|---|---|---|---|
| group | 2 | 0.8338617 | 0.4379319 |
| 84 | NA | NA |
ルビーン検定とは2群以上の分散の均質性を検定する。等分散は棄却されない。
3 箱ひげ図と蜂群図
par(family = "HiraKakuProN-W3") #日本語フォントの指定
library(beeswarm)
boxplot(dat$score~dat$class,ylim=c(0,100),main="３クラスの比較")
beeswarm(dat$score~dat$class,ylim=c(0,100),main="３クラスの比較",pch=16,add=TRUE)
4 anovakun by井関龍太
#http://riseki.php.xdomain.jp/index.php?ANOVA君
source("data/anovakun_485.txt")
#(dat[,-1])#被験者間計画のデータ形式 As
anovakun(dat[,-1],"As",3,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 Wed Jun 28 13:37:23 2023.
##
##
## << DESCRIPTIVE STATISTICS >>
##
## ----------------------------
## A n Mean S.D.
## ----------------------------
## a1 29 60.7931 21.4715
## a2 29 64.0690 20.3977
## a3 29 74.1034 16.7040
## ----------------------------
##
##
## << ANOVA TABLE >>
##
## ----------------------------------------------------------------
## Source SS df MS F-ratio p-value eta^2
## ----------------------------------------------------------------
## A 2789.6782 2 1394.8391 3.6195 0.0311 * 0.0793
## Error 32371.3103 84 385.3727
## ----------------------------------------------------------------
## Total 35160.9885 86 408.8487
## +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 independent means. ==
## == Alpha level is 0.05. ==
##
## ----------------------------
## A n Mean S.D.
## ----------------------------
## a1 29 60.7931 21.4715
## a2 29 64.0690 20.3977
## a3 29 74.1034 16.7040
## ----------------------------
##
## -----------------------------------------------------------
## Pair Diff t-value df p adj.p
## -----------------------------------------------------------
## a1-a3 -13.3103 2.5819 84 0.0116 0.0347 a1 < a3 *
## a2-a3 -10.0345 1.9464 84 0.0549 0.1099 a2 = a3
## a1-a2 -3.2759 0.6354 84 0.5269 0.5269 a1 = a2
## -----------------------------------------------------------
##
##
## output is over --------------------///eta^2 は相関比（correlation ratio）. 順序がないカテゴリカルデータと連続データの「相関」が相関比である。 統計ER.
#水準3、holmの方法で多重比較、効果量ηを表示
#anovakun(データ, “要因計画の型”, 各要因の水準数，…)
#被験者内計画のデータ形式 sA
5 繰り返しのある２元配置分散分析
5.1 ファイルの読み込み
dat.2 <- read.csv("data/data_ch5-2.csv", header = TRUE)
head(dat.2)| student | class | pre | post | delayed |
|---|---|---|---|---|
| S001 | A | 31 | 48 | 30 |
| S002 | A | 39 | 51 | 44 |
| S003 | A | 56 | 67 | 58 |
| S004 | A | 47 | 44 | 50 |
| S005 | A | 29 | 33 | 47 |
| S006 | A | 37 | 41 | 43 |
table(dat.2$class)##
## A B
## 30 305.2 二元配置分散分析
anovakun(dat.2[, -1], "AsB", 2, 3, auto = TRUE, holm = TRUE, eta = TRUE)##
## [ AsB-Type Design ]
##
## This output was generated by anovakun 4.8.5 under R version 4.2.2.
## It was executed on Wed Jun 28 13:37:24 2023.
##
##
## << DESCRIPTIVE STATISTICS >>
##
## --------------------------------
## A B n Mean S.D.
## --------------------------------
## a1 b1 30 37.7333 9.4465
## a1 b2 30 49.5667 11.0068
## a1 b3 30 45.6333 10.3340
## a2 b1 30 38.3667 10.5813
## a2 b2 30 40.1333 12.0279
## a2 b3 30 40.5000 11.3068
## --------------------------------
##
##
## << SPHERICITY INDICES >>
##
## == Mendoza's Multisample Sphericity Test and Epsilons ==
##
## -------------------------------------------------------------------------
## Effect Lambda approx.Chi df p LB GG HF CM
## -------------------------------------------------------------------------
## B 0.0007 14.1756 5 0.0145 * 0.5000 0.9994 1.0350 1.0314
## -------------------------------------------------------------------------
## LB = lower.bound, GG = Greenhouse-Geisser
## HF = Huynh-Feldt-Lecoutre, CM = Chi-Muller
##
##
## << ANOVA TABLE >>
##
## == Adjusted by Greenhouse-Geisser's Epsilon for Suggested Violation ==
##
## --------------------------------------------------------------------
## Source SS df MS F-ratio p-value eta^2
## --------------------------------------------------------------------
## A 970.6889 1 970.6889 3.3048 0.0742 + 0.0412
## s x A 17035.9556 58 293.7234
## --------------------------------------------------------------------
## B 1491.7444 2 746.3374 26.1221 0.0000 *** 0.0633
## A x B 765.4111 2 382.9442 13.4032 0.0000 *** 0.0325
## s x A x B 3312.1778 115.93 28.5711
## --------------------------------------------------------------------
## Total 23575.9778 179 131.7094
## +p < .10, *p < .05, **p < .01, ***p < .001
##
##
## << POST ANALYSES >>
##
## < MULTIPLE COMPARISON for "B" >
##
## == Holm's Sequentially Rejective Bonferroni Procedure ==
## == The factor < B > is analysed as dependent means. ==
## == Alpha level is 0.05. ==
##
## ----------------------------
## B n Mean S.D.
## ----------------------------
## b1 60 38.0500 9.9497
## b2 60 44.8500 12.3807
## b3 60 43.0667 11.0467
## ----------------------------
##
## ----------------------------------------------------------
## Pair Diff t-value df p adj.p
## ----------------------------------------------------------
## b1-b2 -6.8000 7.0112 58 0.0000 0.0000 b1 < b2 *
## b1-b3 -5.0167 5.1767 58 0.0000 0.0000 b1 < b3 *
## b2-b3 1.7833 1.8056 58 0.0762 0.0762 b2 = b3
## ----------------------------------------------------------## Warning in sprintf("%s", paste0("< SIMPLE EFFECTS for \"", part.info1, "\"
## INTERACTION >"), : one argument not used by format '%s'## < SIMPLE EFFECTS for "A x B" INTERACTION >## Warning in is.na(charmatch("Error", parttab$source.col)) && !is.na(partepsi):
## 'length(x) = 20 > 1' in coercion to 'logical(1)'## --------------------------------------------------------------------------
## Effect Lambda approx.Chi df p LB GG HF CM
## --------------------------------------------------------------------------
## B at a1 0.0145 8.1778 2 0.0168 * 0.5000 0.7979 0.8370 0.8255
## B at a2 0.0475 5.8830 2 0.0528 + 0.5000 0.8407 0.8866 0.8744
## --------------------------------------------------------------------------
## LB = lower.bound, GG = Greenhouse-Geisser
## HF = Huynh-Feldt-Lecoutre, CM = Chi-Muller
##
## ---------------------------------------------------------------------
## Source SS df MS F-ratio p-value eta^2
## ---------------------------------------------------------------------
## A at b1 6.0167 1 6.0167 0.0598 0.8077 ns 0.0010
## Er at b1 5834.8333 58 100.6006
## ---------------------------------------------------------------------
## A at b2 1334.8167 1 1334.8167 10.0429 0.0024 ** 0.1476
## Er at b2 7708.8333 58 132.9109
## ---------------------------------------------------------------------
## A at b3 395.2667 1 395.2667 3.3692 0.0716 + 0.0549
## Er at b3 6804.4667 58 117.3184
## ---------------------------------------------------------------------
## B at a1 2179.0889 1.6 1365.5040 39.1960 0.0000 *** 0.1915
## s x B at a1 1612.2444 46.28 34.8378
## ---------------------------------------------------------------------
## B at a2 78.0667 1.68 46.4302 1.3318 0.2706 ns 0.0070
## s x B at a2 1699.9333 48.76 34.8633
## ---------------------------------------------------------------------
## +p < .10, *p < .05, **p < .01, ***p < .001
##
##
## < MULTIPLE COMPARISON for "B at a1" >
##
## == Holm's Sequentially Rejective Bonferroni Procedure ==
## == The factor < B at a1 > is analysed as dependent means. ==
## == Alpha level is 0.05. ==
##
## -----------------------------------------------------------
## Pair Diff t-value df p adj.p
## -----------------------------------------------------------
## b1-b2 -11.8333 8.6553 29 0.0000 0.0000 b1 < b2 *
## b1-b3 -7.9000 7.7556 29 0.0000 0.0000 b1 < b3 *
## b2-b3 3.9333 2.4150 29 0.0223 0.0223 b2 > b3 *
## -----------------------------------------------------------
##
## output is over --------------------///5.3 スタック形式に変更
x <- stack(dat.2[, 3 : 5])
attach(x)
head(x)| values | ind |
|---|---|
| 31 | pre |
| 39 | pre |
| 56 | pre |
| 47 | pre |
| 29 | pre |
| 37 | pre |
# データフレームの作成
y <- data.frame(dat.2$class, x)
head(y)| dat.2.class | values | ind |
|---|---|---|
| A | 31 | pre |
| A | 39 | pre |
| A | 56 | pre |
| A | 47 | pre |
| A | 29 | pre |
| A | 37 | pre |
# 因子の型に変更
y$dat.2.class <- factor(y$dat.2.class)
# 水準の順序を指定
y$ind <- factor(y$ind, levels = c("pre", "post", "delayed"))
# データフレームの列名を変更
names(y) <- c("class", "score", "test")
# 交互作用確認プロット
interaction.plot(y$test, y$class, y$score, type = "b", pch = c(1, 2), xlab = "Test", ylab = "Score", trace.label = "Class")
attach(y)
res <- aov(score ~ class+test)
summary(res)## Df Sum Sq Mean Sq F value Pr(>F)
## class 1 971 970.7 8.092 0.00497 **
## test 2 1492 745.9 6.218 0.00246 **
## Residuals 176 21114 120.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ff <- summary(res)[[1]]["F value"][1,]#F値の取り出し
1-1/ff## [1] 0.8764142alpha(dat.2[, 3:5])## Number of categories should be increased in order to count frequencies.##
## Reliability analysis
## Call: alpha(x = dat.2[, 3:5])
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.89 0.89 0.85 0.73 8.1 0.025 42 10 0.72
##
## 95% confidence boundaries
## lower alpha upper
## Feldt 0.83 0.89 0.93
## Duhachek 0.84 0.89 0.93
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## pre 0.87 0.88 0.78 0.78 7.0 0.032 NA 0.78
## post 0.83 0.83 0.72 0.72 5.0 0.043 NA 0.72
## delayed 0.81 0.82 0.69 0.69 4.5 0.047 NA 0.69
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## pre 60 0.87 0.89 0.79 0.75 38 9.9
## post 60 0.92 0.91 0.85 0.80 45 12.4
## delayed 60 0.92 0.92 0.86 0.81 43 11.0``` #被験者間は２水準、被験者内は３水準 # スタック形式に変更 #https://yaginogogo.hatenablog.jp/entry/2016/04/22/011327