小杉先生の「分散分析アルゴリズム」を

金明哲先生のコード.
教員データ分析入門のコード（anovakun).
http://riseki.php.xdomain.jp/index.php?ANOVA君/ANOVA君の使い方.

で実行。手計算は小杉先生のテキストが詳しい.

1 被験者間要因

1.1 １要因の分散分析（p2）

1.1.1 Jin

a1<- c(28,29,29,25) 
a2<- c(30,29,28,30) 
a3<- c(35,32,38,32) 
bunsan1<- data.frame(A=factor(c(rep("a1",4),rep("a2",4),rep("a3",4))),y=c(a1,a2,a3)) 
bunsan1

A	y
a1	28
a1	29
a1	29
a1	25
a2	30
a2	29
a2	28
a2	30
a3	35
a3	32
a3	38
a3	32

boxplot(y~A,data=bunsan1,col="lightblue")

aov(y~A,data=bunsan1)

## Call:
##    aov(formula = y ~ A, data = bunsan1)
## 
## Terms:
##                        A Residuals
## Sum of Squares  92.66667  38.25000
## Deg. of Freedom        2         9
## 
## Residual standard error: 2.061553
## Estimated effects may be unbalanced

summary(aov(y~A,data=bunsan1))

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## A            2  92.67   46.33    10.9 0.00394 **
## Residuals    9  38.25    4.25                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

当サイト「分散分析と確率モデル」を参考に回帰分析を適用する。
分散分析は回帰分析の特別な場合である。

fit = lm( formula = y ~ A, data = bunsan1 )
f <- fit$ fitted.values;f#予測値

##     1     2     3     4     5     6     7     8     9    10    11    12 
## 27.75 27.75 27.75 27.75 29.25 29.25 29.25 29.25 34.25 34.25 34.25 34.25

a <- fit$ coefficients;a#回帰係数

## (Intercept)         Aa2         Aa3 
##       27.75        1.50        6.50

e <- fit$ residuals;e#残差

##     1     2     3     4     5     6     7     8     9    10    11    12 
##  0.25  1.25  1.25 -2.75  0.75 -0.25 -1.25  0.75  0.75 -2.25  3.75 -2.25

se <- sum(e^2);se#残差平方和

## [1] 38.25

#se/9残差分散
se/9

## [1] 4.25

sqrt(se/9)#残差の標準誤差

## [1] 2.061553

summary(fit)

## 
## Call:
## lm(formula = y ~ A, data = bunsan1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.750 -1.500  0.500  0.875  3.750 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   27.750      1.031  26.921 6.52e-10 ***
## Aa2            1.500      1.458   1.029  0.33034    
## Aa3            6.500      1.458   4.459  0.00158 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.062 on 9 degrees of freedom
## Multiple R-squared:  0.7078, Adjusted R-squared:  0.6429 
## F-statistic:  10.9 on 2 and 9 DF,  p-value: 0.003939

決定係数は、eta^2 と同じ値になっている。

1.1.2 教育データ分析入門.

教育データ分析はanovakunによる。（多重比較を含む）.
当サイト「2_0_anova君」参照.
etaは効果量（イータ）
eta^2は主効果や交互作用の偏差の2乗和をデータ全体の偏差2乗和で割った値。
つまり，その主効果あるいは交互作用の影響が，データ全体のばらつきに比べてどの程度大きいかを示す。.

source("data/anovakun_485.txt")
#(dat[,-1])#被験者間計画のデータ形式　As
anovakun(bunsan1,"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 Sun Jan 28 22:21:48 2024.
## 
##  
## << DESCRIPTIVE STATISTICS >>
## 
## ---------------------------
##   A   n     Mean    S.D. 
## ---------------------------
##  a1   4  27.7500  1.8930 
##  a2   4  29.2500  0.9574 
##  a3   4  34.2500  2.8723 
## ---------------------------
## 
## 
## << ANOVA TABLE >>
## 
## ------------------------------------------------------------
##  Source       SS  df      MS  F-ratio  p-value      eta^2 
## ------------------------------------------------------------
##       A  92.6667   2 46.3333  10.9020   0.0039 **  0.7078 
##   Error  38.2500   9  4.2500                              
## ------------------------------------------------------------
##   Total 130.9167  11 11.9015                              
##                 +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   4  27.7500  1.8930 
##  a2   4  29.2500  0.9574 
##  a3   4  34.2500  2.8723 
## ---------------------------
## 
## ----------------------------------------------------------
##   Pair     Diff  t-value  df       p   adj.p            
## ----------------------------------------------------------
##  a1-a3  -6.5000   4.4590   9  0.0016  0.0047  a1 < a3 * 
##  a2-a3  -5.0000   3.4300   9  0.0075  0.0150  a2 < a3 * 
##  a1-a2  -1.5000   1.0290   9  0.3303  0.3303  a1 = a2   
## ----------------------------------------------------------
## 
## 
## output is over --------------------///

1.2 二要因の分散分析(p5).

1.2.1 Jin

a1<- c(28,29,29,25) 
a2<- c(30,29,28,30) 
a3<- c(35,32,38,32) 
bunsan3<-data.frame(A=factor(c(rep("a1",6),rep("a2",6))),B=factor(rep(c(rep("b1",3),rep("b2",3)),2)),y= c(a1,a2,a3))
#head(bunsan3)
bunsan3

A	B	y
a1	b1	28
a1	b1	29
a1	b1	29
a1	b2	25
a1	b2	30
a1	b2	29
a2	b1	28
a2	b1	30
a2	b1	35
a2	b2	32
a2	b2	38
a2	b2	32

#交互作用なし
summary(aov(y~A+B,data=bunsan3))

##             Df Sum Sq Mean Sq F value Pr(>F)  
## A            1  52.08   52.08   6.271 0.0336 *
## B            1   4.08    4.08   0.492 0.5009  
## Residuals    9  74.75    8.31                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#交互作用あり
aov(y~A*B,data=bunsan3)

## Call:
##    aov(formula = y ~ A * B, data = bunsan3)
## 
## Terms:
##                        A        B      A:B Residuals
## Sum of Squares  52.08333  4.08333 10.08333  64.66667
## Deg. of Freedom        1        1        1         8
## 
## Residual standard error: 2.84312
## Estimated effects may be unbalanced

summary(aov(y~A*B,data=bunsan3))

##             Df Sum Sq Mean Sq F value Pr(>F)  
## A            1  52.08   52.08   6.443 0.0348 *
## B            1   4.08    4.08   0.505 0.4974  
## A:B          1  10.08   10.08   1.247 0.2965  
## Residuals    8  64.67    8.08                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fm1 <- y~B*A #学年*学部
#二元配置分散分析
source("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/socialStatisticsBasic.R", encoding="UTF-8")
output.twoway.anova <- twoway.factorial.anova(fm1,bunsan3,interaction=T)

interaction.plot(x.factor=bunsan3$B,trace.factor=bunsan3$A,response=bunsan3$y, type="b",  pch=c(1,2))

1.2.2 教育データ分析入門.

1.2.2.1 ファイルの読み込み

#dat.2 <- read.csv("data/data_ch5-2.csv", header = TRUE)
head(bunsan3)

A	B	y
a1	b1	28
a1	b1	29
a1	b1	29
a1	b2	25
a1	b2	30
a1	b2	29

table(bunsan3$A)

## 
## a1 a2 
##  6  6

1.2.2.2 二元配置分散分析

anovakun(bunsan3, "ABs", 2, 2, auto = TRUE, holm = TRUE, eta = TRUE)

## 
## [ ABs-Type Design ]
## 
## This output was generated by anovakun 4.8.5 under R version 4.2.2.
## It was executed on Sun Jan 28 22:21:49 2024.
## 
##  
## << DESCRIPTIVE STATISTICS >>
## 
## -------------------------------
##   A   B   n     Mean    S.D. 
## -------------------------------
##  a1  b1   3  28.6667  0.5774 
##  a1  b2   3  28.0000  2.6458 
##  a2  b1   3  31.0000  3.6056 
##  a2  b2   3  34.0000  3.4641 
## -------------------------------
## 
## 
## << ANOVA TABLE >>
## 
## ------------------------------------------------------------
##  Source       SS  df      MS  F-ratio  p-value      eta^2 
## ------------------------------------------------------------
##       A  52.0833   1 52.0833   6.4433   0.0348 *   0.3978 
##       B   4.0833   1  4.0833   0.5052   0.4974 ns  0.0312 
##   A x B  10.0833   1 10.0833   1.2474   0.2965 ns  0.0770 
##   Error  64.6667   8  8.0833                              
## ------------------------------------------------------------
##   Total 130.9167  11 11.9015                              
##                 +p < .10, *p < .05, **p < .01, ***p < .001
## 
## 
## output is over --------------------///

群間変動＝Aの効果＋Bの効果＋交互作用.
群内変動＝Error.

繰り返しがなければ群内変動はない。この場合、交互作用を誤差として扱う。

1.2.3 青木

青木繁伸先生の例題で実行してみる。

1.2.3.1 各水準の繰返し数が 1 である場合.

http://aoki2.si.gunma-u.ac.jp/lecture/TwoWayANOVA/TwoWay1.html

1.2.3.2 Jin

a1<- c(9,17,12,16) 
a2<- c(1,21,16,11) 
a3<- c(7,19,6,9) 
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))
#head(bunsan3)
bunsan3

A	B	y
a1	b1	9
a1	b2	17
a1	b3	12
a1	b4	16
a2	b1	1
a2	b2	21
a2	b3	16
a2	b4	11
a3	b1	7
a3	b2	19
a3	b3	6
a3	b4	9

#交互作用なし
summary(aov(y~A+B,data=bunsan3))

##             Df Sum Sq Mean Sq F value Pr(>F)  
## A            2  21.50   10.75   0.659 0.5510  
## B            3 268.67   89.56   5.492 0.0372 *
## Residuals    6  97.83   16.31                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1.2.3.3 anovakun.

anovakun(bunsan3, "ABs", 3, 4, auto = TRUE, holm = TRUE, eta = TRUE)

## 
## [ ABs-Type Design ]
## 
## This output was generated by anovakun 4.8.5 under R version 4.2.2.
## It was executed on Sun Jan 28 22:21:49 2024.
## 
##  
## << DESCRIPTIVE STATISTICS >>
## 
## -----------------------------
##   A   B   n     Mean  S.D. 
## -----------------------------
##  a1  b1   1   9.0000    NA 
##  a1  b2   1  17.0000    NA 
##  a1  b3   1  12.0000    NA 
##  a1  b4   1  16.0000    NA 
##  a2  b1   1   1.0000    NA 
##  a2  b2   1  21.0000    NA 
##  a2  b3   1  16.0000    NA 
##  a2  b4   1  11.0000    NA 
##  a3  b1   1   7.0000    NA 
##  a3  b2   1  19.0000    NA 
##  a3  b3   1   6.0000    NA 
##  a3  b4   1   9.0000    NA 
## -----------------------------
## 
## 
## << ANOVA TABLE >>
## 
## ------------------------------------------------------------
##  Source       SS  df      MS  F-ratio  p-value      eta^2 
## ------------------------------------------------------------
##       A  21.5000   2 10.7500                       0.0554 
##       B 268.6667   3 89.5556                       0.6924 
##   A x B  97.8333   6 16.3056                       0.2521 
##   Error   0.0000   0     NaN                              
## ------------------------------------------------------------
##   Total 388.0000  11 35.2727                              
##                 +p < .10, *p < .05, **p < .01, ***p < .001
## 
## 
## output is over --------------------///

anovakunは交互作用を計算してしまうので、aovの方が良い。

1.2.4 情報処理実習(京都大学).

https://kyoto-edu.sakura.ne.jp/ryukoku/?&year=2019&id=thu04&course=statistics&content=R04_twoway-anova&showsAns=true

https://cogpsy.educ.kyoto-u.ac.jp/personal/Kusumi/datasem15/ishiguro1.pdf

井関龍太 http://riseki.php.xdomain.jp/index.php/

data <- read.csv("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/self-study2.csv",na.strings=99999, fileEncoding = "utf-8")
fm1 <- time~grade*faculty #学年*学部
fm2 <- time~faculty*grade #学部*学年

#二元配置分散分析
source("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/socialStatisticsBasic.R", encoding="UTF-8")
output.twoway.anova <- twoway.factorial.anova(fm1,data,interaction=T)

1.2.5 まとめ.

要因の効果を除いた変動が誤差であり、考えるモデルによって誤差は異なる。有意な結果が得られなければ、モデルを再検討し分析し直すことになる。

2 被験者内計画

被験者内計画は対応のあるデータで構成される。このため、個人差の変動を誤差から分離することが可能になる。それにより検出力の高い検定が可能になる。

2.1 被験者内の１要因分散分析（p9）

被験者と時期を要因とする２元配置分散分析として計算され、繰り返しがないので交互作用が誤差になる。（p11参照）

id <- c("a","b","c","d","e","f")
a1<- c(25,23,22,26,28,32) 
a2<- c(25,30,25,30,23,34) 
a3<- c(32,29,32,32,35,37) 
bunsan1<- data.frame(ID=factor(id),前期=a1,中期=a2,後期=a3)
bunsan1

ID	前期	中期	後期
a	25	25	32
b	23	30	29
c	22	25	32
d	26	30	32
e	28	23	35
f	32	34	37

anovakunを用いて被験者内の分散分析を行う。

sは個人差.
Aは要因による効果.
s*Aは誤差（交互作用）. 事後分析として多重比較が行われている。

anovakun(bunsan1[-1], "sA",3, auto = TRUE, 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 Sun Jan 28 22:21:50 2024.
## 
##  
## << DESCRIPTIVE STATISTICS >>
## 
## ---------------------------
##   A   n     Mean    S.D. 
## ---------------------------
##  a1   6  26.0000  3.6332 
##  a2   6  27.8333  4.1673 
##  a3   6  32.8333  2.7869 
## ---------------------------
## 
## 
## << SPHERICITY INDICES >>
## 
## == Mendoza's Multisample Sphericity Test and Epsilons ==
## 
## -------------------------------------------------------------------------
##  Effect  Lambda  approx.Chi  df      p         LB     GG     HF     CM 
## -------------------------------------------------------------------------
##       A  0.0727      4.1943   2 0.1228 ns  0.5000 0.6062 0.6963 0.5000 
## -------------------------------------------------------------------------
##                               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 
## ------------------------------------------------------------
##       s 123.7778   5 24.7556                              
## ------------------------------------------------------------
##       A 150.1111   2 75.0556  11.0556   0.0029 **  0.4392 
##   s x A  67.8889  10  6.7889                              
## ------------------------------------------------------------
##   Total 341.7778  17 20.1046                              
##                 +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   6  26.0000  3.6332 
##  a2   6  27.8333  4.1673 
##  a3   6  32.8333  2.7869 
## ---------------------------
## 
## ----------------------------------------------------------
##   Pair     Diff  t-value  df       p   adj.p            
## ----------------------------------------------------------
##  a1-a3  -6.8333   9.7179   5  0.0002  0.0006  a1 < a3 * 
##  a2-a3  -5.0000   2.6600   5  0.0449  0.0898  a2 = a3   
##  a1-a2  -1.8333   1.1033   5  0.3201  0.3201  a1 = a2   
## ----------------------------------------------------------
## 
## 
## output is over --------------------///

bunsan2<- data.frame(A=factor(c(rep("前期",6),rep("中期",6),rep("後期",6))),y=c(a1,a2,a3)) 
bunsan2

A	y
前期	25
前期	23
前期	22
前期	26
前期	28
前期	32
中期	25
中期	30
中期	25
中期	30
中期	23
中期	34
後期	32
後期	29
後期	32
後期	32
後期	35
後期	37

aov(y~A,data=bunsan2)

## Call:
##    aov(formula = y ~ A, data = bunsan2)
## 
## Terms:
##                        A Residuals
## Sum of Squares  150.1111  191.6667
## Deg. of Freedom        2        15
## 
## Residual standard error: 3.574602
## Estimated effects may be unbalanced

summary(aov(y~A,data=bunsan2))

##             Df Sum Sq Mean Sq F value Pr(>F)  
## A            2  150.1   75.06   5.874 0.0131 *
## Residuals   15  191.7   12.78                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

時期による差があるか検定する場合、個人の誤差を考慮しなければ、誤差分散は 12.78 と大きくなるので P値は 0.0131 となる。
被験者内計画では個人の変動を分離するので、誤差分散は 6.7889 と小さくなる。
このため、時期による差のP値は 0.0029 となり検出力が上がり、異常を検出しやすくなる。
これは乱傀法に他ならない。

2.2 被験者内の２要因分散分析（p12）

当サイト「2_0_anova君」参照.

id <- c("a","b","c","d","e","f")
a1b1<- c(22,27,25,30,31,32) 
a1b2<- c(27,32,25,28,33,34) 
a2b1<- c(30,25,29,26,29,31) 
a2b2<- c(31,27,32,29,32,32) 
bunsan3<- data.frame(ID=factor(id),a1b1,a1b2,a2b1,a2b2)
bunsan3

ID	a1b1	a1b2	a2b1	a2b2
a	22	27	30	31
b	27	32	25	27
c	25	25	29	32
d	30	28	26	29
e	31	33	29	32
f	32	34	31	32

群間変動（セル間変動）

Aは教員の違いによる効果　2.0417
Bは教育効果要因（前後）による効果　26.0417
A*Bは交互作用による効果　0.0417
合計　2.0417+26.0417+0.0417=28.1251

群内変動（セル内変動）
　誤差変動 190.49（自由度 20）は大きすぎて検出力が落ちる。
　誤差変動をさらに分解する。

sは個人差 85.8750
s*Aは教員要因による誤差.
　教師要因による変動は 171.13 であり、教師による効果と個人差を引くと、誤差変動 171.13-85.875-2.0417=83.21 である。
s*Bは教育効果要因による誤差.
　教育効果要因による変動は 118.13 であり、教育効果による効果と個人差を引くと、誤差変動 118.13-26.04-85.88=6.21 である。
s* A*Bは交互作用による誤差.
　15.21.
合計の郡内変動=個人差+要因Aの誤差+要因Bの誤差+交互作用の誤差.
　85.875+83.21+6.21+15.21=190.49.

被験者内計画では、誤差の中で個人変動を分離して算出することで、要因による検出力を高めている。考え方は乱塊法と同じであるが、計算の手間が多い。

小杉先生のテキスト p16-p17 に丁寧にまとめられている。

anovakunでは、事後分析として多重比較が行われている。

anovakun(bunsan3[-1], "sAB",2,2, auto = TRUE, holm = TRUE, eta = TRUE)

## 
## [ sAB-Type Design ]
## 
## This output was generated by anovakun 4.8.5 under R version 4.2.2.
## It was executed on Sun Jan 28 22:21:50 2024.
## 
##  
## << DESCRIPTIVE STATISTICS >>
## 
## -------------------------------
##   A   B   n     Mean    S.D. 
## -------------------------------
##  a1  b1   6  27.8333  3.8687 
##  a1  b2   6  29.8333  3.6560 
##  a2  b1   6  28.3333  2.3381 
##  a2  b2   6  30.5000  2.0736 
## -------------------------------
## 
## 
## << SPHERICITY INDICES >>
## 
## == Mendoza's Multisample Sphericity Test and Epsilons ==
## 
## -------------------------------------------------------------------------
##  Effect  Lambda  approx.Chi  df      p         LB     GG     HF     CM 
## -------------------------------------------------------------------------
##  Global  0.0005     11.1986   5 0.0547 +   0.3333 0.4968 0.6595 0.4038 
##       A  1.0000     -0.0000   0            1.0000 1.0000 1.0000 1.0000 
##       B  1.0000     -0.0000   0            1.0000 1.0000 1.0000 1.0000 
##   A x B  1.0000     -0.0000   0            1.0000 1.0000 1.0000 1.0000 
## -------------------------------------------------------------------------
##                               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 
## --------------------------------------------------------------
##         s  85.8750   5 17.1750                              
## --------------------------------------------------------------
##         A   2.0417   1  2.0417   0.1227   0.7404 ns  0.0093 
##     s x A  83.2083   5 16.6417                              
## --------------------------------------------------------------
##         B  26.0417   1 26.0417  20.9732   0.0059 **  0.1191 
##     s x B   6.2083   5  1.2417                              
## --------------------------------------------------------------
##     A x B   0.0417   1  0.0417   0.0137   0.9114 ns  0.0002 
## s x A x B  15.2083   5  3.0417                              
## --------------------------------------------------------------
##     Total 218.6250  23  9.5054                              
##                   +p < .10, *p < .05, **p < .01, ***p < .001
## 
## 
## output is over --------------------///

3 データの構造式

分散分析をする際には、データの構造式を明確にする必要がある。交互作用が有意でない場合は構造式を変更し分析をし直すことになり、簡単にはいかない。
これまで出てきた分析について構造式を整理しておく。

3.1 被験者間計画.

１要因.
全変動＝（要因）＋誤差＝（群間変動）＋群内変動.
２要因.
全変動＝（要因A＋要因B＋交互作用）＋誤差＝（群間変動）＋群内変動.

3.2 被験内計画.

１要因.
全変動＝（要因）＋個人差＋誤差＝（群間変動）＋群内変動.
２要因.
全変動＝（要因A＋要因B＋交互作用）＋個人差＋要因A誤差＋要因B誤差＋交互作用誤差
　　　＝（群間変動）＋群内変動.

3.3 混合計画.

A被験者間・B被験者内.
全変動＝（要因A＋個人差）＋要因B＋交互作用＋被験者内誤差
　　　＝（群間変動）＋群内変動.

4 混合計画

被験者内計画では、個人差を誤差から分離することで検出力の高い検定を行うことができた。
被験者間計画では、個人差は間要因の差を検出するための基準として用いられる。
混合計画は間要因と内要因を含むので、

個人差を誤差として間要因の差を検定する.
（誤差分散はプールした分散を用いる）
内要因の誤差を用いて、交互作用と内要因を検定する。

このように、間要因検定の基準とする誤差と、内要因検定の基準とする誤差が異なるのが特徴である。

id <- c("a","b","c","d","e","f","g","h","i","j","k","l")
a1<- c(28,29,29,25,30,29) 
b1<- c(28,30,35,32,38,32) 
a2<- c(33,27,30,31,30,28) 
b2<- c(32,32,35,41,36,41) 
bunsan4<- data.frame(ID=factor(id),クラス=factor(c(rep("A",6),rep("B",6))),事前=c(a1,b1),事後=c(a2,b2))
bunsan4

ID	クラス	事前	事後
a	A	28	33
b	A	29	27
c	A	29	30
d	A	25	31
e	A	30	30
f	A	29	28
g	B	28	32
h	B	30	32
i	B	35	35
j	B	32	41
k	B	38	36
l	B	32	41

anovakun(bunsan4[-1], "AsB",2,2, 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 Sun Jan 28 22:21:50 2024.
## 
##  
## << DESCRIPTIVE STATISTICS >>
## 
## -------------------------------
##   A   B   n     Mean    S.D. 
## -------------------------------
##  a1  b1   6  28.3333  1.7512 
##  a1  b2   6  29.8333  2.1370 
##  a2  b1   6  32.5000  3.5637 
##  a2  b2   6  36.1667  4.0702 
## -------------------------------
## 
## 
## << SPHERICITY INDICES >>
## 
## == Mendoza's Multisample Sphericity Test and Epsilons ==
## 
## -------------------------------------------------------------------------
##  Effect  Lambda  approx.Chi  df      p         LB     GG     HF     CM 
## -------------------------------------------------------------------------
##       B  0.7546      0.5067   1 0.4766 ns  1.0000 1.0000 1.0000 1.0000 
## -------------------------------------------------------------------------
##                               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 165.3750   1 165.3750  15.7375   0.0027 **  0.4166 
##     s x A 105.0833  10  10.5083                              
## ---------------------------------------------------------------
##         B  40.0417   1  40.0417   5.0420   0.0486 *   0.1009 
##     A x B   7.0417   1   7.0417   0.8867   0.3686 ns  0.0177 
## s x A x B  79.4167  10   7.9417                              
## ---------------------------------------------------------------
##     Total 396.9583  23  17.2591                              
##                    +p < .10, *p < .05, **p < .01, ***p < .001
## 
## 
## output is over --------------------///

4.1 被験者内変動.

被験者内変動は、
　　\(\dfrac{(-5)^2+2^2+\cdots+(-9)^2}{2}=\dfrac{253}{2}=126.5\).
　　前後の要因による変動は、40.04.
　　交互作用による変動は、7.04.
　　被験者内変動からこれらを減じたものが被験者ない誤差になる。
　　即ち、　126.5-40.04-7.04=79.42.
複雑なので論証に自信はないが、
　　全変動＝間変動＋（内変動）.
　　　　　＝間要因＋個人差＋（内要因＋交互作用＋内誤差）.
このように変動を分解していると思う。

4.2 全変動の分解

　　個人ー全平均＝(個人-個平均)＋(個平均-全平均)＝(内変動)＋(間変動)＝126+270.
　　　　　　　　
　　(内変動)＝(個人ー前後別の平均)＋(前後別の平均ー全平均)＝(個人ー前後別の平均)＋40
　　　　　＝(個人ー前後クラス別平均)＋(前後クラス別平均ー前後別の平均)＋40
　　　　　＝79.416+7.041+40.041＝126
　　　　　（内変動の分解は自信がない）

変動を分解することが分散分析の本質であり、モデルによって分解の仕方が異なるということだと思う。

小杉先生がp21で説明されている通り次のことが重要である。

間要因のの変動は個人差を間要因誤差としてと比較すること.
内要因の変動と交互作用は被験者内誤差と比較すること.

4.3 参考.

a1<- c(28,29,29,25,30,29) 
b1<- c(28,30,35,32,38,32) 
a2<- c(33,27,30,31,30,28) 
b2<- c(32,32,35,41,36,41) 
bunsan5<- data.frame(ID=factor(id),A=factor(c(rep("A1",6),rep("B1",6),rep("A2",6),rep("B2",6))),
 B=factor(c(rep("前",12),rep("後",12))),y=c(a1,b1,a2,b2))
bunsan5

ID	A	B	y
a	A1	前	28
b	A1	前	29
c	A1	前	29
d	A1	前	25
e	A1	前	30
f	A1	前	29
g	B1	前	28
h	B1	前	30
i	B1	前	35
j	B1	前	32
k	B1	前	38
l	B1	前	32
a	A2	後	33
b	A2	後	27
c	A2	後	30
d	A2	後	31
e	A2	後	30
f	A2	後	28
g	B2	後	32
h	B2	後	32
i	B2	後	35
j	B2	後	41
k	B2	後	36
l	B2	後	41

aov(y~B+A,data=bunsan5)

## Call:
##    aov(formula = y ~ B + A, data = bunsan5)
## 
## Terms:
##                         B         A Residuals
## Sum of Squares   40.04167 172.41667 184.50000
## Deg. of Freedom         1         2        20
## 
## Residual standard error: 3.037269
## 1 out of 5 effects not estimable
## Estimated effects may be unbalanced

5 交互作用が有意でないとき.

5.1 データ.

a1<- c(76,78,77,79,77,78,74,75) 
a2<- c(73,72,75,75,75,73,72,74) 

bunsan6<-data.frame(A=factor(c(rep("a1",8),rep("a2",8))),B=factor(rep(c(rep("b1",2),rep("b2",2),rep("b3",2),rep("b4",2)),2)),y= c(a1,a2))
#head(bunsan6)
bunsan6

A	B	y
a1	b1	76
a1	b1	78
a1	b2	77
a1	b2	79
a1	b3	77
a1	b3	78
a1	b4	74
a1	b4	75
a2	b1	73
a2	b1	72
a2	b2	75
a2	b2	75
a2	b3	75
a2	b3	73
a2	b4	72
a2	b4	74

summary(aov(y~A*B,data=bunsan6))

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## A            1  39.06   39.06  32.895 0.000436 ***
## B            3  17.19    5.73   4.825 0.033392 *  
## A:B          3   4.69    1.56   1.316 0.334956    
## Residuals    8   9.50    1.19                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#anovakun(bunsan6, "ABs", 2, 4, auto = TRUE, holm = TRUE, eta = TRUE)

5.2 プーリング.

summary(aov(y~A+B,data=bunsan6))

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## A            1  39.06   39.06  30.286 0.000185 ***
## B            3  17.19    5.73   4.442 0.028160 *  
## Residuals   11  14.19    1.29                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A	B	y
a1	b1	28
a1	b1	29
a1	b1	29
a1	b2	25
a1	b2	30
a1	b2	29
a2	b1	28
a2	b1	30
a2	b1	35
a2	b2	32
a2	b2	38
a2	b2	32

A	B	y
a1	b1	9
a1	b2	17
a1	b3	12
a1	b4	16
a2	b1	1
a2	b2	21
a2	b3	16
a2	b4	11
a3	b1	7
a3	b2	19
a3	b3	6
a3	b4	9

ID	クラス	事前	事後
a	A	28	33
b	A	29	27
c	A	29	30
d	A	25	31
e	A	30	30
f	A	29	28
g	B	28	32
h	B	30	32
i	B	35	35
j	B	32	41
k	B	38	36
l	B	32	41

ID	A	B	y
a	A1	前	28
b	A1	前	29
c	A1	前	29
d	A1	前	25
e	A1	前	30
f	A1	前	29
g	B1	前	28
h	B1	前	30
i	B1	前	35
j	B1	前	32
k	B1	前	38
l	B1	前	32
a	A2	後	33
b	A2	後	27
c	A2	後	30
d	A2	後	31
e	A2	後	30
f	A2	後	28
g	B2	後	32
h	B2	後	32
i	B2	後	35
j	B2	後	41
k	B2	後	36
l	B2	後	41

A	B	y
a1	b1	76
a1	b1	78
a1	b2	77
a1	b2	79
a1	b3	77
a1	b3	78
a1	b4	74
a1	b4	75
a2	b1	73
a2	b1	72
a2	b2	75
a2	b2	75
a2	b3	75
a2	b3	73
a2	b4	72
a2	b4	74

A	B	y
a1	b1	28
a1	b1	29
a1	b1	29
a1	b2	25
a1	b2	30
a1	b2	29
a2	b1	28
a2	b1	30
a2	b1	35
a2	b2	32
a2	b2	38
a2	b2	32

A	B	y
a1	b1	9
a1	b2	17
a1	b3	12
a1	b4	16
a2	b1	1
a2	b2	21
a2	b3	16
a2	b4	11
a3	b1	7
a3	b2	19
a3	b3	6
a3	b4	9

ID	クラス	事前	事後
a	A	28	33
b	A	29	27
c	A	29	30
d	A	25	31
e	A	30	30
f	A	29	28
g	B	28	32
h	B	30	32
i	B	35	35
j	B	32	41
k	B	38	36
l	B	32	41

ID	A	B	y
a	A1	前	28
b	A1	前	29
c	A1	前	29
d	A1	前	25
e	A1	前	30
f	A1	前	29
g	B1	前	28
h	B1	前	30
i	B1	前	35
j	B1	前	32
k	B1	前	38
l	B1	前	32
a	A2	後	33
b	A2	後	27
c	A2	後	30
d	A2	後	31
e	A2	後	30
f	A2	後	28
g	B2	後	32
h	B2	後	32
i	B2	後	35
j	B2	後	41
k	B2	後	36
l	B2	後	41

A	B	y
a1	b1	76
a1	b1	78
a1	b2	77
a1	b2	79
a1	b3	77
a1	b3	78
a1	b4	74
a1	b4	75
a2	b1	73
a2	b1	72
a2	b2	75
a2	b2	75
a2	b3	75
a2	b3	73
a2	b4	72
a2	b4	74

A	B	y
a1	b1	28
a1	b1	29
a1	b1	29
a1	b2	25
a1	b2	30
a1	b2	29
a2	b1	28
a2	b1	30
a2	b1	35
a2	b2	32
a2	b2	38
a2	b2	32

A	B	y
a1	b1	9
a1	b2	17
a1	b3	12
a1	b4	16
a2	b1	1
a2	b2	21
a2	b3	16
a2	b4	11
a3	b1	7
a3	b2	19
a3	b3	6
a3	b4	9

ID	クラス	事前	事後
a	A	28	33
b	A	29	27
c	A	29	30
d	A	25	31
e	A	30	30
f	A	29	28
g	B	28	32
h	B	30	32
i	B	35	35
j	B	32	41
k	B	38	36
l	B	32	41

ID	A	B	y
a	A1	前	28
b	A1	前	29
c	A1	前	29
d	A1	前	25
e	A1	前	30
f	A1	前	29
g	B1	前	28
h	B1	前	30
i	B1	前	35
j	B1	前	32
k	B1	前	38
l	B1	前	32
a	A2	後	33
b	A2	後	27
c	A2	後	30
d	A2	後	31
e	A2	後	30
f	A2	後	28
g	B2	後	32
h	B2	後	32
i	B2	後	35
j	B2	後	41
k	B2	後	36
l	B2	後	41

A	B	y
a1	b1	76
a1	b1	78
a1	b2	77
a1	b2	79
a1	b3	77
a1	b3	78
a1	b4	74
a1	b4	75
a2	b1	73
a2	b1	72
a2	b2	75
a2	b2	75
a2	b3	75
a2	b3	73
a2	b4	72
a2	b4	74