1 Friedman検定（p184）

データを順位化し、カイ二乗統計量に帰着するのがFriedman検定である。
Wilcoxonは順位和を正規分布に近似するものであったが、Friedmanはどのような数理によるのであろうか。
大学院農学生命科学研究科教授 岩田洋佳 様 のRpub「ノンパラメトリック検定」に丁寧な説明がある。
乱塊法のノンパラメトリック版というのが適切な説明であろう。
Engineering Skills 様の　Friedman検定にも丁寧な説明がある。

乱塊法は個人差を分離する被験者内計画であり、内容は２元配置分散分析であるので、乱数を用いて２元配置のデータを生成しシミュレーションしてみよう。

2 シミュレーション.

2.1 有意なデータ.

library(tidyr)
data <- crossing(x = 1:4, y = 1:4);#data
str(data)
## tibble [16 × 2] (S3: tbl_df/tbl/data.frame)
##  $ x: int [1:16] 1 1 1 1 2 2 2 2 3 3 ...
##  $ y: int [1:16] 1 2 3 4 1 2 3 4 1 2 ...
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 <- 1*x + 2*y+ rnorm(16, sd = 0.5)#誤差分散 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)
dat <- data.frame(x=x,y=y,z=z)
str(dat)
## 'data.frame':    16 obs. of  3 variables:
##  $ x: int  1 1 1 1 2 2 2 2 3 3 ...
##  $ y: int  1 2 3 4 1 2 3 4 1 2 ...
##  $ z: num  2.9 4.86 6.77 9.18 4.06 ...
dat1 <- spread(dat,key=y,value=z)
dat1
x 1 2 3 4
1 2.899245 4.856536 6.766449 9.184857
2 4.061834 6.547841 7.560887 10.030950
3 5.672720 6.727266 9.924311 10.634314
4 5.674074 8.120242 10.238387 12.966805

2.2 被験者内計画（乱塊法）

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:33:49 2024.
## 
##  
## << DESCRIPTIVE STATISTICS >>
## 
## ---------------------------
##   A   n     Mean    S.D. 
## ---------------------------
##  a1   4   4.5770  1.3521 
##  a2   4   6.5630  1.3372 
##  a3   4   8.6225  1.7202 
##  a4   4  10.7042  1.6213 
## ---------------------------
## 
## 
## << SPHERICITY INDICES >>
## 
## == Mendoza's Multisample Sphericity Test and Epsilons ==
## 
## -------------------------------------------------------------------------
##  Effect  Lambda  approx.Chi  df      p         LB     GG     HF     CM 
## -------------------------------------------------------------------------
##       A  0.0161      4.7412   5 0.5114 ns  0.3333 0.5704 1.2532 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  24.9277   3  8.3092                              
## ------------------------------------------------------------
##       A  83.5793   3 27.8598  93.3955   0.0000 *** 0.7517 
##   s x A   2.6847   9  0.2983                              
## ------------------------------------------------------------
##   Total 111.1916  15  7.4128                              
##                 +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   4.5770  1.3521 
##  a2   4   6.5630  1.3372 
##  a3   4   8.6225  1.7202 
##  a4   4  10.7042  1.6213 
## ---------------------------
## 
## ----------------------------------------------------------
##   Pair     Diff  t-value  df       p   adj.p            
## ----------------------------------------------------------
##  a1-a3  -4.0455  17.4899   3  0.0004  0.0024  a1 < a3 * 
##  a2-a4  -4.1413  14.2012   3  0.0008  0.0038  a2 < a4 * 
##  a1-a4  -6.1273  12.7596   3  0.0010  0.0042  a1 < a4 * 
##  a1-a2  -1.9860   5.9650   3  0.0094  0.0283  a1 < a2 * 
##  a2-a3  -2.0595   4.5910   3  0.0194  0.0388  a2 < a3 * 
##  a3-a4  -2.0817   4.5036   3  0.0204  0.0388  a3 < a4 * 
## ----------------------------------------------------------
## 
## 
## output is over --------------------///

要因yによる変動は有意となっている。

2.2.1 要因yで有意でないデータ.

source("data/anovakun_485.txt")
w <- 1*x + 0*y+ rnorm(16, sd = 0.1)#誤差分散 0.1^2
dat <- data.frame(x=x,y=y,w=w)
str(dat)
## 'data.frame':    16 obs. of  3 variables:
##  $ x: int  1 1 1 1 2 2 2 2 3 3 ...
##  $ y: int  1 2 3 4 1 2 3 4 1 2 ...
##  $ w: num  1.04 1.05 1.01 1.01 1.99 ...
dat2 <- spread(dat,key=y,value=w)
dat2
x 1 2 3 4
1 1.038865 1.053322 1.005396 1.014327
2 1.990383 2.080404 2.221801 2.154945
3 3.046071 2.966603 3.099124 2.889309
4 4.056317 4.011248 4.171194 3.863273 
anovakun((dat2[,-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:33:49 2024.
## 
##  
## << DESCRIPTIVE STATISTICS >>
## 
## --------------------------
##   A   n    Mean    S.D. 
## --------------------------
##  a1   4  2.5329  1.3052 
##  a2   4  2.5279  1.2606 
##  a3   4  2.6244  1.3418 
##  a4   4  2.4805  1.2021 
## --------------------------
## 
## 
## << SPHERICITY INDICES >>
## 
## == Mendoza's Multisample Sphericity Test and Epsilons ==
## 
## -------------------------------------------------------------------------
##  Effect  Lambda  approx.Chi  df      p         LB     GG     HF     CM 
## -------------------------------------------------------------------------
##       A  0.0236      4.3016   5 0.5662 ns  0.3333 0.6074 1.4966 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 19.5523   3 6.5174                              
## ----------------------------------------------------------
##       A  0.0434   3 0.0145   2.1093   0.1693 ns  0.0022 
##   s x A  0.0617   9 0.0069                              
## ----------------------------------------------------------
##   Total 19.6575  15 1.3105                              
##               +p < .10, *p < .05, **p < .01, ***p < .001
## 
## 
## output is over --------------------///

要因yによる変動は有意でない。

2.2.2 順位によるFriedman検定.

ノンパラメトリックな手法であるFriedmanは順位統計量を用いる。

r<-c(rank(z[1:4]),rank(z[5:8]),rank(z[9:12]),rank(z[13:16]))
#r<-c(rank(w[1:4]),rank(w[5:8]),rank(w[9:12]),rank(w[13:16]))
datr <- mutate(dat,r=r)
str(datr)
## 'data.frame':    16 obs. of  4 variables:
##  $ x: int  1 1 1 1 2 2 2 2 3 3 ...
##  $ y: int  1 2 3 4 1 2 3 4 1 2 ...
##  $ w: num  1.04 1.05 1.01 1.01 1.99 ...
##  $ r: num  1 2 3 4 1 2 3 4 1 2 ...
boxplot(z~y,data=datr)

boxplot(r~y,data=datr)

dat1 <- datr[,-3]
#dat2 <- dat1[,-2]
dat3 <-spread(dat1,key=y,value=r)
dat3
x 1 2 3 4
1 1 2 3 4
2 1 2 3 4
3 1 2 3 4
4 1 2 3 4

2.2.3 Friedman統計量.

岩田先生の資料によれば、
　　　\(S =\displaystyle \dfrac{12}{bk(k+1)}\sum_{j=1}^k \biggl[ R_j - \frac{b(k+1)}{2}\biggr]^2\) 　がFriedman統計量である。
Engineering Skills様によれば、次のように変形することができる。
　　　\(\displaystyle {{\chi}_{0}}^{2} =\frac {12}{nk(k+1)} \sum_{i=1}^{k}{ {R_i}^2 } - 3n(k+1)\).
　　　\(\displaystyle {{\chi}_{0}}^{2} =\frac {12n}{k(k+1)} \left( \sum_{i=1}^{k}{ {r_i}^2 } - \frac{k{(k+1)}^2}{4} \right)\).
　　　 適合度検定なので、岩田先生の式がとりつきやすい。
変形は高校数学で可能であるが、ここでは立ち入らない。
市原先生は２番目の式で計算しているので、第２式で算出しよう。

colSums(dat3[,-1])
##  1  2  3  4 
##  4  8 12 16
R <- sum(colSums(dat3[,-1])^2)
R
## [1] 480
chi <- (12/(4*4*5))*R-3*4*5
chi
## [1] 12

Friedman統計量は12となり、大きな値となっている。

2.2.4 要因yによる効果がない場合.

library(tidyr)
data <- crossing(x = 1:4, y = 1:4)
x <- data$x
y <- data$y
x1 <- rank(runif(4)) ;x2 <-rank(runif(4));x3 <- rank(runif(4)) ;x4 <-rank(runif(4))
z <-c(x1,x2,x3,x4)
dat <- data.frame(x=x,y=y,z=z)
dat1 <-spread(dat,key=y,value=z)[,-1]
dat1
1 2 3 4
3 2 1 4
2 3 4 1
3 2 4 1
4 3 1 2 
apply(dat1, 2, sum)
##  1  2  3  4 
## 12 10 10  8
R <-sum(apply(dat1, 2, sum)^2)
R
## [1] 408
chi <- (12/(4*4*5))*R-3*4*5
chi
## [1] 1.2
ch <-6*var(apply(dat1, 2, mean))*3/2.5
ch
## [1] 1.2

chiとchは同じ統計量になるようであるが、確かめには立ち入らない。

要因による効果がない場合、Friedman統計量は小さくなっている。

2.2.5 シミュレーション

以下に帰無仮説のもとでのシミュレーションを示す。
12が有意な値であることがわかる。
一様乱数runifを用いる。

library(tidyr)
data <- crossing(x = 1:4, y = 1:4);#data
x <- data$x
y <- data$y
s = c();t =c()
for (i in 1:1000) {
  x1 <- rank(runif(4)) ;x2 <-rank(runif(4));x3 <- rank(runif(4)) ;x4 <-rank(runif(4))
  z <-c(x1,x2,x3,x4)
  dat <- data.frame(x=x,y=y,z=z)
  dat1 <-spread(dat,key=y,value=z)[,-1]
  R <-sum(apply(dat1, 2, sum)^2)
  chi <- (12/(4*4*5))*R-3*4*5
  ch <-6*var(apply(dat1, 2, mean))*3/2.5
  s <-c(s,chi)
  t <-c(t,ch)
}
hist(s,prob=T)
curve(dchisq(x,3),add=T)

hist(t,prob=T)
curve(dchisq(x,3),add=T)

mean(s)
## [1] 3.081
var(s)
## [1] 4.765685
mean(t)
## [1] 3.081
var(t)
## [1] 4.765685

https://power-of-awareness.com/entry/2022/09/27/050000

https://oceanone.hatenablog.com/entry/2021/08/11/020933

期待値はカイ二乗統計量のものだが、カイ二乗分布との一致は良いとはいえない。分散は異なっている。
近似するには k>4 である必要がある。 k=5 としてみる。

library(tidyr)
data <- crossing(x = 1:5, y = 1:5);#data
x <- data$x
y <- data$y
s = c()
for (i in 1:1000) {
  x1 <- rank(runif(5)) ;x2 <-rank(runif(5));x3 <- rank(runif(5)) ;x4 <-rank(runif(5));x5 <-rank(runif(5))
  z <-c(x1,x2,x3,x4,x5)
  dat <- data.frame(x=x,y=y,z=z)
  dat1 <-spread(dat,key=y,value=z)[,-1]
  R <-sum(apply(dat1, 2, sum)^2)
  chi <- (12/(5*5*6))*R-3*5*6
  s <-c(s,chi)
}
hist(s,prob=T)
curve(dchisq(x,4),add=T)

mean(s)
## [1] 4.0168
var(s)
## [1] 6.744305

ヒストグラムはカイ二乗分布に近くなり、期待値も理論値に近いが、分散はあてはまりがよくないようである。

3 例6-2（p186）.

3.1 long型データ.

a1<- c(293,240,208,260) 
a2<- c(154,145,130,167) 
a3<- c(180,130,150,160) 
a4<- c(195,160,180,172) 
data<-data.frame(A=factor(c(rep("a1",4),rep("a2",4),rep("a3",4),rep("a4",4))),B=factor(rep(c(rep("b1",1),rep("b2",1),rep("b3",1),rep("b4",1)),4)),y= c(a1,a2,a3,a4))
data
A B y
a1 b1 293
a1 b2 240
a1 b3 208
a1 b4 260
a2 b1 154
a2 b2 145
a2 b3 130
a2 b4 167
a3 b1 180
a3 b2 130
a3 b3 150
a3 b4 160
a4 b1 195
a4 b2 160
a4 b3 180
a4 b4 172

3.2 wide型データ.

dat1 <- spread(data,key=B,value=y)
dat1
A b1 b2 b3 b4
a1 293 240 208 260
a2 154 145 130 167
a3 180 130 150 160
a4 195 160 180 172

3.3 患者による被験者内計画.

これは患者をブロック要因とする乱塊法と同じである。

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:33:55 2024.
## 
##  
## << DESCRIPTIVE STATISTICS >>
## 
## -----------------------------
##   A   n      Mean     S.D. 
## -----------------------------
##  a1   4  205.5000  60.7426 
##  a2   4  168.7500  49.0535 
##  a3   4  167.0000  34.1955 
##  a4   4  189.7500  47.0912 
## -----------------------------
## 
## 
## << SPHERICITY INDICES >>
## 
## == Mendoza's Multisample Sphericity Test and Epsilons ==
## 
## -------------------------------------------------------------------------
##  Effect  Lambda  approx.Chi  df      p         LB     GG     HF     CM 
## -------------------------------------------------------------------------
##       A  0.0754      2.9673   5 0.7454 ns  0.3333 0.6214 1.6018 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 26005.5000   3 8668.5000                              
## ----------------------------------------------------------------
##       A  4042.5000   3 1347.5000   4.9642   0.0266 *   0.1244 
##   s x A  2443.0000   9  271.4444                              
## ----------------------------------------------------------------
##   Total 32491.0000  15 2166.0667                              
##                     +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  205.5000  60.7426 
##  a2   4  168.7500  49.0535 
##  a3   4  167.0000  34.1955 
##  a4   4  189.7500  47.0912 
## -----------------------------
## 
## -----------------------------------------------------------
##   Pair      Diff  t-value  df       p   adj.p            
## -----------------------------------------------------------
##  a2-a4  -21.0000   5.6805   3  0.0108  0.0649  a2 = a4   
##  a1-a2   36.7500   3.6556   3  0.0354  0.1768  a1 = a2   
##  a1-a3   38.5000   2.4362   3  0.0928  0.3713  a1 = a3   
##  a3-a4  -22.7500   1.6930   3  0.1890  0.5671  a3 = a4   
##  a1-a4   15.7500   1.5785   3  0.2126  0.5671  a1 = a4   
##  a2-a3    1.7500   0.1343   3  0.9017  0.9017  a2 = a3   
## -----------------------------------------------------------
## 
## 
## output is over --------------------///

<< SPHERICITY INDICES >> 球面性の検定であるが、有意となっていないので分散分析が可能と判断できる。

3.3.1 aov関数との比較.

str(data)
## 'data.frame':    16 obs. of  3 variables:
##  $ A: Factor w/ 4 levels "a1","a2","a3",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ B: Factor w/ 4 levels "b1","b2","b3",..: 1 2 3 4 1 2 3 4 1 2 ...
##  $ y: num  293 240 208 260 154 145 130 167 180 130 ...
boxplot(y~B,data=data)

summary(aov(y~A,data=data))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## A            3  26006    8669   16.04 0.000169 ***
## Residuals   12   6485     540                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(y~B,data=data))
##             Df Sum Sq Mean Sq F value Pr(>F)
## B            3   4043    1348   0.568  0.646
## Residuals   12  28448    2371
summary(aov(y~A+B,data=data))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## A            3  26006    8669  31.935 3.98e-05 ***
## B            3   4043    1348   4.964   0.0266 *  
## Residuals    9   2443     271                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(aov(y~B+A,data=data))
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## B            3   4043    1348   4.964   0.0266 *  
## A            3  26006    8669  31.935 3.98e-05 ***
## Residuals    9   2443     271                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

乱塊法（被験者内計画）は2元配置分散分析と同じ計算をしている。ブロック間変動を誤差から取り除いて１元配置分散分析を行うことになり検出力が上がる。

3.3.2 Friedman検定.

分散分析は正規性と等分散性を前提としている。
ノンパラメトリックであるFriedman検定をRの関数で実行してみる。

friedman.test(y=as.matrix(dat1[,-1]))
## 
##  Friedman rank sum test
## 
## data:  as.matrix(dat1[, -1])
## Friedman chi-squared = 8.1, df = 3, p-value = 0.04399

P値は大きめになっており、検定が厳し目になっていることがわかる。

4 例6-3（p187）.

4.1 long型データ.

a1<- c(493,582,530,497,711,565) 
a2<- c(929,730,493,524,414,877) 
a3<- c(497,581,382,330,534,1134) 
a4<- c(432,615,451,392,1437,1442) 
a5<- c(729,818,463,440,437,550) 
data<-data.frame(A=factor(c(rep("a1",6),rep("a2",6),rep("a3",6),rep("a4",6),rep("a5",6))),B=factor(rep(c(rep("b1",1),rep("b2",1),rep("b3",1),rep("b4",1),rep("b5",1),rep("b6",1)),5)),y= c(a1,a2,a3,a4,a5))
data
A B y
a1 b1 493
a1 b2 582
a1 b3 530
a1 b4 497
a1 b5 711
a1 b6 565
a2 b1 929
a2 b2 730
a2 b3 493
a2 b4 524
a2 b5 414
a2 b6 877
a3 b1 497
a3 b2 581
a3 b3 382
a3 b4 330
a3 b5 534
a3 b6 1134
a4 b1 432
a4 b2 615
a4 b3 451
a4 b4 392
a4 b5 1437
a4 b6 1442
a5 b1 729
a5 b2 818
a5 b3 463
a5 b4 440
a5 b5 437
a5 b6 550

4.2 wide型データ.

dat1 <- spread(data,key=B,value=y)
dat1
A b1 b2 b3 b4 b5 b6
a1 493 582 530 497 711 565
a2 929 730 493 524 414 877
a3 497 581 382 330 534 1134
a4 432 615 451 392 1437 1442
a5 729 818 463 440 437 550 
anovakun((dat1[-1]),"sA",6,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:33:55 2024.
## 
##  
## << DESCRIPTIVE STATISTICS >>
## 
## ------------------------------
##   A   n      Mean      S.D. 
## ------------------------------
##  a1   5  616.0000  208.4970 
##  a2   5  665.2000  104.9748 
##  a3   5  463.8000   54.9791 
##  a4   5  436.6000   78.5035 
##  a5   5  706.6000  424.7332 
##  a6   5  913.6000  381.7228 
## ------------------------------
## 
## 
## << SPHERICITY INDICES >>
## 
## == Mendoza's Multisample Sphericity Test and Epsilons ==
## 
## -------------------------------------------------------------------------
##  Effect  Lambda  approx.Chi  df      p         LB     GG     HF     CM 
## -------------------------------------------------------------------------
##       A  0.0000     89.7376  14 0.0000 *** 0.2000 0.2998 0.4394 0.2000 
## -------------------------------------------------------------------------
##                               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  232275.1333   4  58068.7833                              
## --------------------------------------------------------------------
##       A  763391.7667   5 152678.3533   2.3013   0.0834 +   0.3287 
##   s x A 1326872.0667  20  66343.6033                              
## --------------------------------------------------------------------
##   Total 2322538.9667  29  80087.5506                              
##                         +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   5  616.0000  208.4970 
##  a2   5  665.2000  104.9748 
##  a3   5  463.8000   54.9791 
##  a4   5  436.6000   78.5035 
##  a5   5  706.6000  424.7332 
##  a6   5  913.6000  381.7228 
## ------------------------------
## 
## ------------------------------------------------------------
##   Pair       Diff  t-value  df       p   adj.p            
## ------------------------------------------------------------
##  a2-a4   228.6000   4.8746   4  0.0082  0.1229  a2 = a4   
##  a2-a3   201.4000   4.0862   4  0.0150  0.2103  a2 = a3   
##  a4-a6  -477.0000   2.4590   4  0.0698  0.9070  a4 = a6   
##  a3-a6  -449.8000   2.4177   4  0.0729  0.9070  a3 = a6   
##  a1-a4   179.4000   2.3543   4  0.0781  0.9070  a1 = a4   
##  a3-a4    27.2000   1.7092   4  0.1626  1.0000  a3 = a4   
##  a1-a3   152.2000   1.7017   4  0.1640  1.0000  a1 = a3   
##  a5-a6  -207.0000   1.4744   4  0.2144  1.0000  a5 = a6   
##  a4-a5  -270.0000   1.3277   4  0.2550  1.0000  a4 = a5   
##  a1-a6  -297.6000   1.3159   4  0.2586  1.0000  a1 = a6   
##  a3-a5  -242.8000   1.2621   4  0.2755  1.0000  a3 = a5   
##  a2-a6  -248.4000   1.2619   4  0.2755  1.0000  a2 = a6   
##  a1-a2   -49.2000   0.7597   4  0.4897  1.0000  a1 = a2   
##  a1-a5   -90.6000   0.3464   4  0.7465  1.0000  a1 = a5   
##  a2-a5   -41.4000   0.1919   4  0.8572  1.0000  a2 = a5   
## ------------------------------------------------------------
## 
## 
## output is over --------------------///

<< SPHERICITY INDICES >> 球面性の検定は有意となっているので、ノンパラメトリックな検定を適用する必要がある。

バートレット検定でも、

bartlett.test(formula=data$y~data$B)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  data$y by data$B
## Bartlett's K-squared = 21.034, df = 5, p-value = 0.0007981

等分散仮説は棄却される。

4.2.1 Friedman検定.

friedman.test(y=as.matrix(dat1[,-1]))
## 
##  Friedman rank sum test
## 
## data:  as.matrix(dat1[, -1])
## Friedman chi-squared = 10.943, df = 5, p-value = 0.05252

参考サイト.
https://zenn.dev/tmizuho/books/3d511e017bfd23/viewer/61b380