1 平均の比較

1.1 比較の原理

2つのグループの間で平均を比較したい場合にはt検定というものを行う．これは厳密には「2つのグループは元々は同じ母集団から抽出されたサンプルである」という仮説（これを帰無仮説と呼ぶ）の元で，両グループの平均値の間に有意といえるほどの差があるかどうかを検定する手法である．

「元々同じ母集団である」という仮説の元では，理想的には両群の平均値は一致する，すなわち差は0になるはずである．実際にはデータにはばらつきは付きものなので，完全に0となることはないだろうが，差の値は0に近い値になる確率は高いだろし，0から遠い値になる確率は低いだろう．

そこで，実際に2つのグループのデータから、それぞれの平均値の差を求めたうえで，「2つのグループは同じ母集団から抽出されたサンプルである」という仮説のもとで，そのような差が出てくる確率はいくらであったかを算出する，というのがt検定の数学的な中身である．

そして，その確率が0.05，すなわち5%未満であったならば「2つは同じ母集団から抽出されたサンプルである」という仮定のもとでは極めて珍しい（発生確率が5%未満でしかない）値が出た，ということになる．ここで解釈の転換をして，「そんな珍しい値が出るなんてそうそうあることではない．そもそも元の仮説が誤っていたのではないか」と考える．つまり，「2つは同じ母集団から抽出されたサンプルである」という仮説を否定するのである．この仮説を否定するということは，「2つのグループは互いに性質の異なる（それゆえ母平均も異なる）母集団から抽出されたサンプルである」という仮説を支持する，ということである．こうして，2つのグループには「母集団のレベルで差がある」と主張するのがt検定である．

1.2 t検定の方法（基本）

Rでt検定を行うにはt.test()関数を用いる．この関数は、平均の信頼区間で用いた関数であるが、本来は名前のとおりt検定を行う関数である．引数は以下の通りである．

第1引数：比較先のベクトルデータ
第2引数：比較元のベクトルデータ
第3引数alternative（省略可）:「両側検定」を行うのか「片側検定」を行うの指定．デフォルトでは両側検定となっていおり，通常は両側検定を行うのが一般的であるため，省略するのが一般的である．特に片側検定にしたい場合には，比較先の平均が比較元の平均より「高い」ということを検証したい場合には “greater”，「低い」ということを検証したい “less”を指定する．

．いずれもベクトルデータで与える必要がある点に注意が必要である．

両側検定と片側検定

両側検定は「データAの平均がデータBの平均に比べて有意に大きい」「データAの平均がデータBの平均に比べて有意に小さい」「両群に差があるとは言えない」という3つのいずれになるかを知りたい場合に用いる．

一方，片側検定では，事前に「データAの方がデータBに比べて大きい（小さい）はずである」ということが理論的に予想されている中で，「有意に大きい（小さい）」もしくは「大きい（小さい）といえない」という2つのいずれになるかを知りたい場合に用いる．

すなわち，片側検定では逆の結論は全く想定していないことになる．数学的には，両側検定の時に算出される確率（P値）はは片側検定のときのP値の2倍となる．すなわち，両側検定にした場合，有意になりにくくなる（5%を上回りやすくなる）．このため，「有意になってほしい」という思いが強いと，片側検定にしたくなるが，片側検定をする場合，それを選択した理由をきちんと述べねばならない．

一般的には両側検定を用いることが多く，t.test()に限らず，Rで行う統計分析で有意検定を行っているものは両側検定がデフォルトで実行される．

ある農作物に与える肥料を新しく試作した．その効果を確かめるために，新しい肥料を与えた畝と通常の肥料を与えた畝を20ずつ用意し，それぞれの収穫量を測定した．その結果，次のようなデータ（リンク）が得られた．このデータから新しい肥料を使用した場合に収穫量が向上するのかどうかを検証せよ．

この例題の場合，t検定を行うコードは以下の通りとなる．

mydata_ttest <- read.csv( # データ読み込み オブジェクト名は自由につけてよい
# データは以下の設定ではpracticeというフォルダに入っているコトになっているが，
# 特にそうしたフォルダに入れていない場合には ./example_ttest.csv とすればよい．
  "./practice/example_ttest.csv", 
  header=T, 
  fileEncoding = "UTF8"
)

# データの内容を確認
head(mydata_ttest)

##   通常の肥料を使用 新しい肥料を使用
## 1            51.85            53.47
## 2            42.18            46.09
## 3            46.82            54.14
## 4            48.16            61.07
## 5            47.02            64.48
## 6            44.47            52.85

# すべてのデータを見てみたければ
# 右上ペインのEvironmentのmydataをクリックするか，
# 以下のコマンドを実行．
# view(mydata)

# それぞれのベクトルデータを取得
新肥料 <- mydata_ttest$新しい肥料を使用
旧肥料 <- mydata_ttest$通常の肥料を使用


t.test(新肥料, 旧肥料)#Alternativeは省略

## 
##  Welch Two Sample t-test
## 
## data:  新肥料 and 旧肥料
## t = 3.9995, df = 36.981, p-value = 0.0002922
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   3.791649 11.578351
## sample estimates:
## mean of x mean of y 
##    53.645    45.960

2行目のp-value = 0.0002922という箇所が結果を示す箇所であり，また最後の

mean of x mean of y 
   53.645    45.960

は，それぞれのデータの平均値である．

このデータではP値は0.05を下回ってり，「二つのグループは同じ母集団から抽出されたサンプルである」という帰無仮説を棄却された．すなわち，両者の平均は有意な違いがあるという結論となり，新しい肥料を用いることで収穫量の向上が望めることが示されたこととなる．

1.3 対応関係のあるデータでのt検定

対応関係とは，例えば2つのデータが「同じ人の1回目の測定と2回目の測定」というように，同じサンプルから条件を変えて取得したデータである場合には，データの中で対応関係（ペア関係）が成立する．このようなデータの場合には，第4引数として, paired = TRUEというオプションを設定する必要がある．

対応関係は，データの並びによってとられる．すなわち，

データ1の1つ目のデータ⇔データ2の1つ目のデータ
データ1の2つ目のデータ⇔データ2の2つ目のデータ
データ1の3つ目のデータ⇔データ2の3つ目のデータ
・・・

という対応関係がとられる．

なお，対応関係のあるデータでは，当然ながら二つのデータの点数は「同じ」であることが前提となる．従って，paired=TRUEの設定をしている時に，一方のデータの中にNAが1つでも含まれていた場合には，t.test()はエラーを返してくる．これを防ぐためには，さらに第5引数na.rm=TRUEのオプション設定を行う．このオプションを設定していた場合，一方のデータにNAが含まれていた場合，NAに対応したもう一方の値も除去される．

15人の協力者に，新しく試作した食パンと従来の食パンの食べ比べをしてもらい，おいしさを10段階で評価してもらったところ，次のようなデータ（リンク）が得られた．このデータをもとに，試作品と従来品とで美味しさの違いがあるのかを検証せよ．

この例題の場合，同じ人に対してデータを取っているため，二つのベクトルの間には対応関係が存在する．このため，pairedの設定が必要になる．は以下の通りとなる．

mydata_ttest_paired <- read.csv(
  "./practice/example_ttest_paired.csv", 
  header=T, 
  fileEncoding="UTF8"
) # データ読み込み
# もしデータの内容を確認したければ，右上ペインのEvironmentのmydataをクリックするか，以下のコマンドを実行．
# view(mydata)

t.test(mydata_ttest_paired$従来品, mydata_ttest_paired$試作品, paired=TRUE, na.rm=T)

## 
##  Paired t-test
## 
## data:  mydata_ttest_paired$従来品 and mydata_ttest_paired$試作品
## t = -0.24203, df = 14, p-value = 0.8123
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -1.314899  1.048232
## sample estimates:
## mean difference 
##      -0.1333333

今回のデータの場合，p値は0.8123となっており，0.05を大きく上回っており，試作品と従来品とでおいしさに有意な差はない，という結論となる．また，pairedの設定をしたので，表題もPaired t-testとなり，最後の結果も両群の平均値ではなく，群間の差のみが算出される．

1.4 ロング形式データでのt検定

ロング形式のデータの場合，filter()を使って，分析に掛けるベクトルデータをそれぞれ取り出せば，基本のt検定の方法で分析ができるが，いちいち個別にベクトルデータを取りださなくても，測定変数と属性変数を~（チルダ）で結んだ「式」をt.test()関数の第1引数に与えることでも分析ができる．なお，属性変数のデータは要因型である必要がある．この場合，第2引数にはそれぞれの変数を含んだデータフレームのオブジェクト名を与える．

第3引数（alternative）は，上記の通りである．

以下は例題1のデータをロング形式に変換したうえで，式形式でのt検定を行う例である．

library(tidyverse) # pivot_longer関数を使うにはtidyverseを読み込む必要がある

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# ロング形式のデータを作成
long_data <- pivot_longer(
  mydata_ttest,
  cols=colnames(mydata_ttest),
  names_to = "肥料",
  values_to= "収穫量" #測定変数に与える列名
)

# 肥料列を要因型に変換
long_data$肥料 <- as.factor(long_data$肥料)

# 内容確認
summary(long_data)

##                肥料        収穫量     
##  新しい肥料を使用:20   Min.   :31.72  
##  通常の肥料を使用:20   1st Qu.:44.65  
##                        Median :51.30  
##                        Mean   :49.80  
##                        3rd Qu.:54.38  
##                        Max.   :64.48

#t検定
t.test(収穫量~肥料, long_data)

## 
##  Welch Two Sample t-test
## 
## data:  収穫量 by 肥料
## t = 3.9995, df = 36.981, p-value = 0.0002922
## alternative hypothesis: true difference in means between group 新しい肥料を使用 and group 通常の肥料を使用 is not equal to 0
## 95 percent confidence interval:
##   3.791649 11.578351
## sample estimates:
## mean in group 新しい肥料を使用 mean in group 通常の肥料を使用 
##                         53.645                         45.960

1.4.1 式, data形式では対応ありデータの分析はできない！！

ただし，ロング形式のデータとなっている場合に，R 4.4.1（2024.10.9時点の最新）では，式, data形式での分析はできない．対応ありデータ(第4引数paired=TRUEとなるデータ）の場合には，以下の例のようにそれぞれのデータを取り出してt.testにかける必要がある．

## ロング形式のデータを作成
long_data_paired <- pivot_longer(
  mydata_ttest_paired,
  cols=colnames(mydata_ttest_paired),
  names_to = "食パン",
  values_to= "評価" #測定変数に与える列名
)

# 確認
long_data_paired$食パン <- as.factor(long_data_paired$食パン)
summary(long_data_paired)

##     食パン        評価       
##  試作品:15   Min.   : 4.000  
##  従来品:15   1st Qu.: 6.000  
##              Median : 7.000  
##              Mean   : 6.867  
##              3rd Qu.: 8.000  
##              Max.   :10.000

# 以下の形式だとエラーになる．
t.test(評価~食パン, long_data_paired, paired=TRUE, na.rm=TRUE)

## Error in t.test.formula(評価 ~ 食パン, long_data_paired, paired = TRUE, : cannot use 'paired' in formula method

# ロングデータからの食パンごとのデータを取り出し
data1 <- filter(long_data_paired, 食パン=="従来品")
data2 <- filter(long_data_paired, 食パン=="試作品")


# t検定
t.test(data1$評価, data2$評価, paired=TRUE)

## 
##  Paired t-test
## 
## data:  data1$評価 and data2$評価
## t = -0.24203, df = 14, p-value = 0.8123
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -1.314899  1.048232
## sample estimates:
## mean difference 
##      -0.1333333

t検定の3つの種類

t検定には大きく以下の3つの種類がある．

対応するデータでのt検定
Studentのt検定
Welchのt検定

ExcelやGoogle Spread SheetでT検定を行う場合も，「検定の種類」という項目があるが，これは上記の3つのどれを行うかを選択するものである．

対応するデータでのt検定は，すでに説明した通り，同じサンプルから取得したデータを比較する場合に用いる．

一方，StudentのT検定とWelchのT検定はいずれも「対応のないデータ」を比較する場合に用いるが，両者の違いとして，比較先・比較元の両データについて「母集団の分散が等しいとする仮定（等分散性の仮定）を立てるか，立てないか」という点が異なる．

StudentのT検定は，両データの母集団の分散が等しいと仮定して行う検定であるのに対して，WelchのT検定は，母集団の分散が等しいという仮定を立てないで行う検定である．

結論を先に述べると，最近の統計学の流れでは，常にWelchのT検定を使うことが推奨される．実際にRでのt.testはデフォルトではWelchの検定を行う．

一昔前は，以下で述べる分散の比較（F検定）を行い，その結果に基づいてStudentのT検定かWelchのT検定かを選択することが一般的であった．すなわち，F検定の結果，「5%有意で分散が違っている」という結果が出ればWelchのT検定を使い，有意でなければStudentのT検定を行うということが行われていた．

しかし，こうした選択方法は最近では否定的に捉えられるようになっている．その理由として，「有意検定は『等しい』ことの保証にはならない」という点と，さらに本章の後の方で述べる検定の多重性問題がある．

特に1点目について，F検定の結果，「有意確率が5%を下回った」だった場合には「分散は等しくない」ということの強い根拠となる．しかしその逆である「有意確率が5%を上回った」という結果は，「分散は等しい」ということの強い根拠とはならない．あくまで「分散は等しい」という帰無仮説は棄却されなかったというだけに過ぎず，「本当は両者の間の分散は違っている」という可能性が残されているのである．このため，F検定の結果（5%を下回った）に基づいて，Wlechの検定を選択することは正しい選択であると言えるが，F検定の結果（5%を上回った）に基づいてStudentの検定を選択することは誤った選択である可能性が残されてしまう．

さらに，Welchの検定では，たとえ実際には分散が等しかった場合でも，かなり正確な結果を返してくることが様々なパターンでのシミュレーションの結果から明らかになっている[1]．

このため，最近では「分散が等しいかどうか」を気にせずにWelchの検定を使うことが推奨されるようになっている．

[1] Delacre, M., et al (2017). Why Psychologists Should by Default Use Welch’s t-test Instead of Student’s t-test. International Review of Social Psychology, 30(1), 92–101, DOI: https://doi.org/10.5334/irsp.82

1.5 t検定の結果の報告方法

t検定の結果を報告する際には棒グラフを用いるのが一般的である（平均なので箱ひげ図は用いない）．二つの群の平均と標準偏差を棒グラフで表現したうえで，角傘（ブラケットと呼ぶ）をつけてp値を記載する．

t検定は，あくまでそのデータの違いを見ているのではなく，そのデータから推定される母集団に違いがあるかを検証するものである．このため，t検定の結果を報告する際に用いる統計量は，推測統計量となる．このことはt検定だけに限らず，あらゆる分析にあてはまる．今後解説していくあらゆる分析は得られているデータを基に推定される母集団の比較をしたり，母集団間の関係性を分析したりするものである．このため，今後の分析の中で出てくる平均や標準偏差は全て基本的に推測統計量となる．

具体的にggplotを使って，先の2つの例題の結果を棒グラフで表記してみよう．まずは，こちらのコードを参考に棒グラフを作成する．まずは例題1．

library(ggplot2) # ggplotを使うためにはまずこのライブライを読み込んでおく必要がある

#まずは両群の平均と標準偏差を算出する．
mymean_new <- mean(mydata_ttest$新しい肥料を使用)
mymean_old <- mean(mydata_ttest$通常の肥料を使用)
mysd_new <- sd(mydata_ttest$新しい肥料を使用)
mysd_old <- sd(mydata_ttest$通常の肥料を使用)

mydf <- data.frame(
  条件=c("新しい肥料を使用", "通常の肥料を使用"),
  収穫量=c(mymean_new, mymean_old),
  標準偏差=c(mysd_new, mysd_old)
)

# ggplotで棒グラフを表示　
g <- ggplot(data=mydf)+ 
  geom_bar(aes(x=条件, y=収穫量, fill=条件), stat="identity")+
  geom_errorbar(aes(x=条件, ymin=収穫量-標準偏差, ymax=収穫量+標準偏差, width=0.1))
plot(g)

続いて，このグラフにブラケットをつけてp値を記載する．ブラケットをつけるためには，annotate()関数を用いて線分を描画する．この関数は第1引数として"segment"（線分の意味）を指定し，第2～４引数として始点と終点の座標を指定すると，その間に線分を描画する関数である．これを組み合わせてブラケットを描画する．x軸方向は「新しい肥料を使用」「通常の肥料を使用」というように数値ではなくラベルとなっているが，内部では座標として，それぞれ1，2とされているので，これを利用してブラケットを描画する．y座標については，エラーバーと干渉しないような値を設定してやる．以下ではy＝61, yend=64としているが，別にこの値でなくてもよい．

#エラーバーの上に実線を引く
g <- g + 
  annotate("segment",x=1, xend=1, y=61, yend=64)+
  annotate("segment",x=1, xend=2, y=64, yend=64)+
  annotate("segment",x=2, xend=2, y=61, yend=64)

plot(g)

さらに，これにp値を与える．これもannotate()関数を用いる．第1引数で"text"と指定して，x，y引数で座標（記載する文字列の中心位置）を指定し，label引数で表示する文字列を，size引数で文字の大きさを指定する．

# 各傘の上に
g <- g+annotate("text",x=1.5, y=66, label="p=0.0002922, ***", size=3)
plot(g)

p値にはアスタリスクをつけるのが一般的である．アスタリスク（*）の数はp値の大小に応じて変わる．p値が0.05未満であれば1つ，0.01未満であれば2つ，0.001未満であれば3つとする．

逆に有意でなかった場合には，P値を表記するだけで何も記載しないか，あるいは，n.s.と記載する．これはnot significantの略であり，「有意でない」という意味である（さらに多重比較を行うケースでは有意でないものについてはブラケットも含め何も表記しないのが一般的である）．

また，上記の例のようにP値が0.001を下回っている場合（多くの場合〇.〇e-7というように小数が指数表記される）には，桁を省略して，P<0.001と表記することもある．

例題2の方でも同様に棒グラフを作成してみよう．今度は描画関数をすべてまとめて記載する．

#まずは両群の平均と標準偏差を算出する．
mymean_new <- mean(mydata_ttest_paired$試作品)
mymean_old <- mean(mydata_ttest_paired$従来品)
mysd_new <- sd(mydata_ttest_paired$試作品)
mysd_old <- sd(mydata_ttest_paired$従来品)

# データフレームの作成
mydf <- data.frame(
  条件=c("試作品", "従来品"),
  おいしさ=c(mymean_new, mymean_old),
  標準偏差=c(mysd_new, mysd_old)
)

# ggplotで棒グラフを表示
g<- ggplot(data=mydf)+ 
  geom_bar(aes(x=条件, y=おいしさ, fill=条件), stat="identity")+
  geom_errorbar(aes(x=条件, ymin=おいしさ-標準偏差, ymax=おいしさ+標準偏差, width=0.1))+
  annotate("segment",x=1, xend=1, y=9, yend=10)+
  annotate("segment",x=1, xend=2, y=10, yend=10)+
  annotate("segment",x=2, xend=2, y=9, yend=10)+
  annotate("text",x=1.5, y=10.5, label="p=0.8123, n.s.", size=5)
plot(g)

なお，ロング型のデータの場合には，こちらで説明した通り，平均や標準偏差はaggregate()関数を用いて算出するとよい．以下これらを使って棒グラフを表示させる．

# 平均値の算出
平均<-aggregate(収穫量~肥料, long_data, mean)
# 標準偏差の算出
標準偏差<-aggregate(収穫量~肥料, long_data, sd)

#データフレームの作成
mydf <- data.frame(
  条件=平均$肥料,
  収穫量=平均$収穫量,
  標準偏差=標準偏差$収穫量
)

# ggplotで棒グラフを表示
g<- ggplot(data=mydf)+
  geom_bar(aes(x=条件, y=収穫量, fill=条件), stat="identity")+
  geom_errorbar(aes(x=条件, ymin=収穫量-標準偏差, ymax=収穫量+標準偏差, width=0.1))+
  annotate("segment",x=1, xend=1, y=61, yend=64)+
  annotate("segment",x=1, xend=2, y=64, yend=64)+
  annotate("segment",x=2, xend=2, y=61, yend=64)+
  annotate("text",x=1.5, y=66, label="p=0.0002922, ***", size=4)
plot(g)

ある自治体の高校生を対象に朝食摂取と成績の関係を明らかにするために調査をおこない，このようなデータを取得した．この調査では，調査日（抜き打ち実施）に朝食を食べているかどうかと，4限目の授業で実施した数学のテストの得点がデータとして取得されている．このデータについて，以下の課題を実施せよ．

このデータの冒頭6行分をhead()を用いて確認し，このデータがワイド形式かロング形式か確認せよ．
摂取群と欠食群のそれぞれの平均と標準偏差（こちらで述べている通り推測統計量でよい）を算出するとともに，それらを棒グラフで示せ．棒グラフにはエラーバーも付与すること．
朝食を食べている人と食べていない人とで数学のテストの得点に有意な差があるかどうかを検証せよ．
先に描いた棒グラフ上にt検定の結果を示すブラケットを書き入れよ．

あるショコラティエで，新しいチョコレートの試作に取り組んでいる．同じチョコレートでも，どのような形に成形するかで，味の評価が変わるのではないかと考え，店に来店したお客さん50人にそれぞれを食べてもらい，おいしさを10段階で評価してもらい，次のようなデータを取得した．このデータについて，以下の課題を実施せよ．

このデータの冒頭6行分をhead()を用いて確認し，このデータがワイド形式かロング形式か確認せよ．
形状Aと形状Bのそれぞれの平均と標準偏差を算出するとともに，それらを棒グラフで示せ．棒グラフにはエラーバーも付与すること．
形状Aと形状Bとでおいしさに有意な差があるかどうかを検証せよ．
先に描いた棒グラフ上にt検定の結果を示すブラケットとt検定の結果えられた有意確率を書き入れよ．

2 分散の比較

先ほどは2つのグループの間での平均の比較，すなわち「2つのグループの間で平均値の差が意味のある差なのか誤差の範囲に収まる差なのか」を検証する方法について説明した．次に，2つのグループの間で「ばらつき」，すなわち分散に違いがあるのかどうかを検定する方法（F検定）について説明する．

分散比

実際にはRの中での処理なので，意識することはないが，分散の場合には「差」をではなく，「比率」（すなわち割り算をする）で比較をしている．すなわち，2つのグループの分散の比率を求め，それが1（すなわち等しい）かどうかを検定することになる．

平均の場合には差であったものが，なぜ分散では比を取るのか．それは分散はそれぞれのサンプルの値と平均との差の2乗から算出されるため，同じ程度の有意性があったとしても，分散が大きいもの同士を比較している場合には差は大きくなるし，分散が小さいもの同士を比較していると差は小さくなる，というように差の大きさを絶対的な物差しで評価することができないためである．

分散の比較における帰無仮説は「2つの群の母集団の分散は等しい」となる．この帰無仮説が棄却された場合，すなわちP値が有意水準を下回った場合，「2つの群の母集団の分散は等しくない」と結論づけることができる．

分散の比較においては，「対応のあるデータ」についての検定は行わない．その理由として，対応のあるデータでは，「群の分散」よりも「それぞれのペアごとの差」の平均が0となるかどうかが重要なのであり，各群を別個に捉えて分散を比較することに意味はないからである．

2.1 F検定の方法（基本）

F検定は，Rのvar.test()関数を用いて行う．この関数の基本的な引数は以下の通りである．

第1引数：比較先のベクトルデータ
第2引数：比較元のベクトルデータ

第1引数と第2引数に比較したい2つのグループのデータを与える．この関数は，t検定と同様に「2つのグループは同じ母集団から抽出されたサンプルである」という仮説の元で，2つのグループの分散の比率が1であるかどうかを検定する．

ある製品の製造工程を変更したところ，製品の重量にばらつきが生じた．このばらつきが工程変更の影響によるものなのか，単なる偶然によるものなのかを検証するために，工程変更前と工程変更後の製品の重量を測定した．その結果，次のようなデータ（リンク）が得られた．このデータから工程変更によるばらつきの変化があるのかを検証せよ．

Rのスクリプトは以下の通りとなる．

mydata_vartest <- read.csv(
  "./practice/example_vartest.csv", 
  header=T, 
  fileEncoding="UTF8"
) # データ読み込み

変更前 <- mydata_vartest$工程変更前
変更後 <- mydata_vartest$工程変更後

# それぞれの分散を確認しておく
var(変更前)

## [1] 43.05095

var(変更後)

## [1] 123.1153

# F検定
var.test(変更後,変更前 )

## 
##  F test to compare two variances
## 
## data:  変更後 and 変更前
## F = 2.8598, num df = 19, denom df = 19, p-value = 0.02698
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  1.131927 7.225038
## sample estimates:
## ratio of variances 
##           2.859758

この結果から，P値が0.027となっており，5％の有意水準を下回っているため，工程変更によって製品の重量のばらつきに有意な変化が生じたと結論づけることができる．

2.2 ロング形式のデータでのF検定

ｔ検定の時と同様に，ロング形式のデータの場合には測定変数~要因変数, データの形でも検定を行うことができる．以下にその例を示す．

#ワイド形式からロング形式に変換
long_data_vartest <- pivot_longer(
  mydata_vartest,
  cols=colnames(mydata_vartest)[],
  names_to = "工程",
  values_to= "重量" #測定変数に与える列名
) 

# 工程を要因型（要因型）に変換
long_data_vartest$工程 <- as.factor(long_data_vartest$工程)

# F検定の実施
var.test(重量~工程, long_data_vartest)

## 
##  F test to compare two variances
## 
## data:  重量 by 工程
## F = 2.8598, num df = 19, denom df = 19, p-value = 0.02698
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  1.131927 7.225038
## sample estimates:
## ratio of variances 
##           2.859758

2.3 F検定の結果の報告方法

F検定の結果を報告する際には，一般には棒グラフを用いることは少ない．F検定の結果は，P値が有意水準を下回っているかどうかを確認するだけでよい．P値が有意水準を下回っている場合には「2つの群の分散は5%の有意水準のもとで有意な違いがある（F(df1, df2)=F-Value, P=value）」と記載する．df1とdf2はそれぞれの2つの群の自由度（各群のサンプル数―1）を示し，F-ValueはF値を示す．

例えば，上の例題の場合には，「工程変更前と工程変更後の製品の重量の分散に有意な違いがあることが示された（F(19, 19)=2.860, P=0.027）．」と記載する．

ある製品メーカでは公称重量70gの金属部品を製造しているが，製造工程の中でどうしても多少のばらつきが発生してしまうため，69.950g～70.049gまでを許容範囲とし，この範囲を超えるものについては不良品といて品質検査時に弾くことになっている．

今，製品の歩留まり率（検品の際に不良品として弾かれなかった製品の割合）の向上を狙って，製造工程を変更を計画している．この工程変更によって，製品の品質のバラツキが抑えられるかどうかを検証するために，変更前の製造された製品と，試験的に導入した新工程での製品をランダムに200個ずつ取り出してきたデータ．このデータを用いて，以下のことを行え．

このデータの冒頭6行分をhead()を用いて確認し，このデータの形式か確認せよ．
変更前と変更後のそれぞれの平均と標準偏差を算出するとともに，それらを棒グラフで示せ．棒グラフにはエラーバーも付与すること．
変更前と変更後とで品質検査の点数そのものに有意な差があるかどうかを検証せよ．
変更前と変更後とで品質検査の点数のバラツキに有意な差があるかどうかを検証せよ．

3 独立性の検定

続いて比率の比較を行う方法を説明する．比率の検定とは，例えば以下の表のような「クロス集計表」が作られている時に，列方向に並べた属性変数Aと行方向に並べた属性変数Bが互いに独立しているのか，何等かの関係性があるのかを検定する方法である．帰無仮説は「互いに独立している」であり，それが否定することによって「関係性がある」と結論づける統計方法である．

	Bあり	Bなし
Aあり	10	20
Aなし	20	40

3.1 「独立」とは？

独立とは，たとえば上記の表の場合，条件Aが「あり」の場合での条件Bの「あり」と「なし」の比率は1:2であるが，条件Aが「なし」の場合での条件Bの「あり」と「なし」の比率も2:1である．すなわち，条件Aが「あり」か「なし」かによって，条件Bの「あり」か「なし」かのなり易さ（生起確率）に違いはない．つまり，条件Bの「あり」「なし」のいずれになるかには条件Aは関わってこない．このような場合には，条件Aと条件Bが互いに独立していると呼ぶ．

逆に以下の表のような場合，条件Aが「あり」の場合での条件Bの「あり」と「なし」の比率は1:2であるが，条件Aが「なし」の場合での条件Bの「あり」と「なし」の比率は2:1である．すなわち，条件Aが「あり」か「なし」かによって，条件Bの「あり」か「なし」かの生起確率が異なっている．このような場合には，条件Aと条件Bが互いに独立していないと呼ぶ．

	Bあり	Bなし
Aあり	10	20
Aなし	40	20

上記の1つ目の例では両方とも1:2の比で完全に一致しているが，実際には，誤差が入り込む．では，実際に何らかのクロス集計表がえられたときに，比率にどの程度以上の差があれば「独立していない＝関連性がある」とみなすのか？その検定を行うのが独立性の検定である．

3.2 カイ2乗検定とフィッシャーの正確確率検定

独立性の検定を行う方法として，カイ2乗検定とフィッシャーの正確確率検定がある．以下，それぞれの方法のRでの行い方を説明する．

3.2.1 カイ2乗検定の方法

独立性の検定ではカイ2乗($\chi^2$)検定を用いるのが基本である．．カイ2乗検定は，Rのchisq.test()関数を用いて行う．この関数は，以下の2通りの方法で実行可能である．

両要因のベクトルデータを与える．

第1引数：要因Aについてのベクトルデータ
第2引数：要因Bについてのベクトルデータ
第3引数（correct）: FALSEにする．詳しくはこちら

クロス集計したテーブルデータを引数して与える．

第1引数：クロス集計したテーブルデータ
第2引数（correct）: FALSEにする．

通常，比率の検定を行う際にはその前段でクロス集計を行うのが一般的なので，後者の場合だとその結果をそのまま使うことができる．

クロス集計表の作成は，Rのtable()関数を用いて行う．この関数の引数は以下の通りである．

第1引数: 行方向（左側）に並べるベクトルデータ（要因A）
第2引数: 列方向（上側）に並べるベクトルデータ（要因B）

なお，要因A，要因Bにそれぞれ与えるベクトルデータは要因型である必要がある．

また，ここで用いられている要因A，要因Bの各ベクトルデータは互いに対応関係があるデータである．すなわち，例えば要因Aのベクトルの1つ目のデータと要因Bの1つ目のデータは同じサンプルが取られたデータである，ということである．このように対応関係がデータの並びによって取られていくので，分析を掛けるにあたっては，各ベクトルのデータの並びがきちんと対応関係がとれているかに注意しておく必要がある．

ある試験を受験した100人の受験者について，性別(男性, 女性）と合否（不合格, 合格）のデータ（リンク）がある．このデータから，クロス集計表を作成するとともに，性別と合否の間に関連性があるのかどうかを検証せよ．

Rのスクリプトは以下の通りとなる．

mydata_chisq <- read.csv(
  "./practice/example_chisq.csv", 
  header=T, 
  fileEncoding="UTF8"
) # データ読み込み

# 要因型に変換
mydata_chisq$合否 <- as.factor(mydata_chisq$合否)
mydata_chisq$性別 <- as.factor(mydata_chisq$性別)

# データの確認
head(mydata_chisq)

##   性別   合否
## 1 女性   合格
## 2 女性 不合格
## 3 男性 不合格
## 4 女性 不合格
## 5 女性   合格
## 6 女性   合格

summary(mydata_chisq)

##    性別        合否   
##  女性:55   合格  :46  
##  男性:45   不合格:54

#クロス集計表の作成 
# クロス集計表にはtable()関数を用いる.引数には行，列に割り振りたい要因を与える．
# 以下のように引数に名前を与えることもできる．
cross_table <- table("性別"=mydata_chisq$性別, "合否"=mydata_chisq$合否)
cross_table

##       合否
## 性別 合格 不合格
##   女性   27     28
##   男性   19     26

#カイ二乗検定
## 要因Aと要因Bをそれぞれベクトルデータとして与える方法
chisq.test(mydata_chisq$性別, mydata_chisq$合否, correct=FALSE)

## 
##  Pearson's Chi-squared test
## 
## data:  mydata_chisq$性別 and mydata_chisq$合否
## X-squared = 0.47008, df = 1, p-value = 0.493

## クロス集計表の結果を与える方法
chisq.test(cross_table, correct=FALSE)

## 
##  Pearson's Chi-squared test
## 
## data:  cross_table
## X-squared = 0.47008, df = 1, p-value = 0.493

この結果から，P値が0.493 となっており，5%の有意水準を下回っていないため，性別と合否の間に有意な関連性があるとは言えない．すなわち，性別は入試の合否に対して「独立である」（影響を与えていない）という可能性が否定出来ないという結果となった．

Yateの補正

カイ2乗検定の際には，クロス集計表のいずれかのセルが5以下の値の場合には，その値を補正することが一般的である．この補正をYateの補正と呼ぶ．Rのchisq.test()関数では，Yateの補正を行うのがデフォルトとなっている．

ただし最近では，セルが5以下の場合には以下で述べるFisherの正確確率検定を行う方が望ましいとされている．

さらに，サンプル数が十分にある場合にYateの補正を行うと，有意確率の値が本来の値よりも大きくなってしまい，本来は有意であるにもかかわらず，それを見逃すエラーが発生しうる．

このため，chisq.test()を行うときには，correct=FALSEのオプションを付けて，Yateの補正を行わず，純粋なカイ2乗検定の値を出力するようにすることが望ましい．

3.2.2 Fisherの正確確率検定の方法

カイ2乗検定はクロス集計した際のいずれかのセルが5以下の値の場合には正確な値は返してこない場合がある．そこでそうしたケースでも対応できるような方法としてFisherの正確確率検定というものがある．

Fisherの正確確率検定は，RのFisher.test()関数を用いて行う．引数の与え方はchisq.testと同様に2通りの方法が可能である．

ある喫茶店でコーヒーを注文したかどうか（しなかった，した）と，スイーツを注文したかどうか（しなかった，した）のデータ（リンク）がある．このデータから，クロス集計表を作成するとともに，コーヒーを注文するかどうかとスイーツを注文するかどうかの間に関連性があるのかどうかを検証せよ．

Rのスクリプトは以下の通りとなる．

mydata_fisher <- read.csv(
  "./practice/example_fisher.csv", 
  header=T, 
  fileEncoding="UTF8"
) # データ読み込み

mydata_fisher$珈琲 <- as.factor(mydata_fisher$珈琲)
mydata_fisher$甘味 <- as.factor(mydata_fisher$甘味)
summary(mydata_fisher)

##          珈琲            甘味   
##  した      :35   した      :26  
##  しなかった:15   しなかった:24

#クロス集計表の作成
cross_table <- table("珈琲"=mydata_fisher$珈琲, "甘味"=mydata_fisher$甘味)
cross_table

##             甘味
## 珈琲       した しなかった
##   した         23         12
##   しなかった    3         12

この結果から，スイーツもコーヒーも注文していない人が3人しかいないためカイ2乗検定では対応できないデータとなっている．このため，フィッシャーの正確確率検定を行う．

#Fisherの正確確率検定
## 要因Aと要因Bをそれぞれベクトルデータとして与える方法
fisher.test(mydata_fisher$珈琲, mydata_fisher$甘味)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  mydata_fisher$珈琲 and mydata_fisher$甘味
## p-value = 0.004889
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   1.571857 48.315800
## sample estimates:
## odds ratio 
##   7.334964

## クロス集計表の結果を与える方法
fisher.test(cross_table)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  cross_table
## p-value = 0.004889
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   1.571857 48.315800
## sample estimates:
## odds ratio 
##   7.334964

P値が0.0049となっており，1%の有意水準を下回っているため，コーヒーを注文するかしないかとスイーツを注文するかしないかの間に有意な関連性があると結論づけることができる．

3.3 カイ2乗検定とFisherの正確確率検定の使い分け

上記の通り，一般にいずれかのセルが5以下の小さな値になっている場合には，カイ2乗検定では正しい結果が導けなくなるので，Fisherの正確確率検定が良い，と述べた．であるならば，どのような場合でもFisherの正確確率検定を用いればよいではないか，という意見が出てくるだろう．

確かにFisherの正確確率検定を常に用いれば，カイ2乗検定の問題を回避できる．しかし，Fisherの正確確率検定は，カイ2乗検定に比べて計算量が多くなるというデメリットがある．この量はデータのサンプル数に依存しており，サンプル数が大きい場合には計算規模が大きくなりすぎてアプリケーションがフリーズしたり，あるいはRの側でエラーを返してくることもある．そのため，データの規模が大きい場合には，カイ2乗検定を用いることが一般的である．

実際にどれくらいのサンプル数であればカイ2乗検定の方が良いかの目安というものは特には存在しないので，サンプル数が5以上なのであれば，一旦はFisherの正確確率検定を実施してみて，もし何らかの問題が発生すればカイ2乗検定に切り替える，という方針が良いだろう．

3.4 結果の報告方法

クロス集計表の結果を報告する際には，まずはクロス集計表をしめす．その上で，カイ2乗検定を使った場合，その結果に表記されている値を引用しながら，「カイ2乗検定の結果2つの群の間に有意な関連性があることが確認された（$\chi^2$(df)=X-squared, P=p-value）」と文章の中で記載する．dfは自由度を示し，X-squaredはカイ2乗値を，p-valueはP値を示す．一方，フィッシャーの正確確率検定の場合には，「フィッシャーの正確確率検定の結果2つの群の間に有意な関連性があることが確認された（P=p-value）」と記載する．

3.5 3つ以上の水準を持つ要因での独立性の検定

上記の例では，2つの水準を持つ要因同士（つまりクロス集計表が$2 \times 2$となる）での独立性の検定を行ったが，3つ以上の水準を持つ要因での独立性の検定も同様に行うことができる．

ある大学の食堂の利用頻度が学部によって異なるのではないかという疑問が生じた．そこで，ある大学の学生に対して，学部（文学部，理学部，工学部）と食堂の利用頻度（毎日，週に3回以上，週に1回以上，月に1回以上，年に1回以上，利用しない）について尋ねるアンケートを行った．その結果，次のようなデータ（リンク）が得られた．このデータから，学部と食堂の利用頻度の間に関連性があるのかを検証せよ．

mydata_restaurant <- read.csv(
  "./practice/example_Restaurant.csv", 
  header=T, 
  fileEncoding="UTF8"
) # データ読み込み

mydata_restaurant$学部 <- as.factor(mydata_restaurant$学部)
mydata_restaurant$食堂利用頻度 <- as.factor(mydata_restaurant$食堂利用頻度)
summary(mydata_restaurant)

##        X              学部          食堂利用頻度
##  Min.   :   1.0   工学部:700   月に1回以上:139  
##  1st Qu.: 325.8   文学部:200   週に1回以上: 92  
##  Median : 650.5   理学部:400   週に3回以上:275  
##  Mean   : 650.5                年に1回以上: 65  
##  3rd Qu.: 975.2                毎日       :667  
##  Max.   :1300.0                利用しない : 62

#クロス集計表の作成
cross_table <- table(mydata_restaurant$学部, mydata_restaurant$食堂利用頻度)
cross_table

##         
##          月に1回以上 週に1回以上 週に3回以上 年に1回以上 毎日 利用しない
##   工学部          45          27         145          11  462         10
##   文学部          56          40          14          33   26         31
##   理学部          38          25         116          21  179         21

# カイ2乗検定
## 直接要因を与える場合
chisq.test(mydata_restaurant$学部, mydata_restaurant$食堂利用頻度)

## 
##  Pearson's Chi-squared test
## 
## data:  mydata_restaurant$学部 and mydata_restaurant$食堂利用頻度
## X-squared = 381.22, df = 10, p-value < 2.2e-16

## クロス集計表を与える場合
chisq.test(cross_table)

## 
##  Pearson's Chi-squared test
## 
## data:  cross_table
## X-squared = 381.22, df = 10, p-value < 2.2e-16

この結果から，P値が$2.2 \times 10^{-16}$未満という非常に小さな値となっており，有意水準を裕に下回っているため，学部と食堂の利用頻度の間に有意な関連性がある，すなわち，学部によって食堂利用頻度は異なると結論づけることができる．

このように3つ以上の水準を持つ要因での独立性も，chisq.test()で検証することができる．

ちなみ，このデータをfisher.testで行った場合には以下の通りエラーが出力される．

fisher.test(cross_table)

## Error in fisher.test(cross_table): FEXACT error 6 (f5xact).  LDKEY=614 is too small for this problem: kval=54241438.
## Try increasing the size of the workspace.

これは計算量が膨大になり，計算のためのメモリ領域(エラーメッセージのworkplace)が足りないことによって生じたエラーである．メッセージにあるとおり，以下の例の通り，workplaceに大きな数字を設定してやれば（以下の例では2e9=$2\times 10^9$という膨大な大きさのメモリを設定している），計算させることができるがかもしれないが，計算できたとして膨大な計算時間がかかってしまうだろう．

fisher.test(cross_table, workspace=2e9)

サンプル数が大きい場合にはカイ2乗検定でも十分に正確な結果が出力されるので，chisq.testで十分である．

ある大学3年生の学生が，最近後輩たちが使っているメッセンジャーアプリに違いが出てきているように感じている．具体的には，自分達はLINEを主なメッセンジャーツールとして使っているが，最近の後輩はInstagramのDMを主なツールとして使っているようだ．そこで実際に自分達の世代と後輩世代と使っているメッセンジャーアプリに違いがあるのかをゼミ研究として調査をしてみることにした．

調査は，大学1から4年生の学生約2000人に対して，LINEとInstagramのどちらを主なメッセンジャーツールとして使っているかを尋ねるアンケートをメールで送付した．結果，次のようなデータ（リンク）が得られた．このデータから，先輩群(3，4年生）と後輩群(1，2年生）の間で使っているメッセンジャーアプリに違いがあるのかを検証せよ．

ある大学のあるゼミの学生が，ふと「政治学のゼミにいる学生はそのほかのゼミの学生より．選挙に行く人が多いのだろうか」という疑問を抱いた．これを検証するために，実際に最近あった衆議院選挙に投票に行ったかどうかを尋ねるアンケートを行った．その結果，次のようなデータ（リンク）が得られた．このデータから，政治学ゼミの学生とその他のゼミの学生の間で投票に行く人の割合に違いがあるのかを検証せよ．

4 比率の差の検定

独立性の検定と類似したものとして「比率の差の検定」というものもある．独立性の検定が要因Aと要因Bのクロス集計表から，AとBの間に関連性があるかどうかを検証するのに対して，比率の差の検定は，要因Aの2つの水準（例えば，男性と女性）の間で，それぞれの中での何らかの事柄の占める割合（例えば，合格率＝受験者に占める合格者の割合）の差が有意であるかどうかを検証するものである．この場合にはZ検定を用いる．

比率の差の検定における帰無仮説は「2つの水準間で比率に差がない」である．

4.1 Z検定の方法

RでZ検定を行う場合にはprop.test()関数を用いて行う．この関数は，以下の引数を取る．

第1引数: 比較したい2つのグループの成功数（比率を計算するときの分子にくる数）を並べたベクトルデータ．並べ方はc(比較先, 比較元)とする．
第2引数: 比較したい2つのグループの試行総数（（比率を計算するときの分母にくる数）を並べたベクトルデータ．並べ方は同上．
第3引数correct: 「連続性の補正」を行うかどうかのオプション．カイ2乗検定におけるYateの補正と同等のもの．一般にサンプルサイズが5以下の場合にはTRUEにすることでより正確な値を返す．デフォルトではTRUEになっているが，サンプルサイズが十分に大きい場合にはこの補正は不要であるため，明示的にFALSEを与えるほうが良い．

先ほどの合否の例で男女での合格率に差があるかを検証してみよう．第1引数には，今回だと男女それぞれの合格者数をベクトルで与え，第2引数には，男女それぞれの受験者数をベクトルで与えることになる．

このためまずは男女別の受験者数と合格者数を求める必要がある． filter()関数とlength()（ベクトルデータの要素数を返す）やnrow()（データフレームの行数＝サンプル数を返す）を用いてデータを抽出することで求めることもできるが，Rにはsum()関数というものがある．この関数は元々は与えたベクトルデータの総和を返す関数であるが，以下のように検索条件の式を与えることによって，それにマッチしたデータの数を返してくれる．

mydata_chisq <- read.csv(
  "./practice/example_chisq.csv", 
  header=T, 
  fileEncoding="UTF8"
) # データ読み込み

mydata_chisq$合否 <- as.factor(mydata_chisq$合否)
mydata_chisq$性別 <- as.factor(mydata_chisq$性別)
summary(mydata_chisq)

##    性別        合否   
##  女性:55   合格  :46  
##  男性:45   不合格:54

# クロス集計表の作成
cross_table <- table("性別"=mydata_chisq$性別, "合否"=mydata_chisq$合否)
print(cross_table)

##       合否
## 性別 合格 不合格
##   女性   27     28
##   男性   19     26

この結果から，合格者は女性で27人，男性で19人．受験者数は女性で55人，男性で45人であることがわかる．

男女別の合格者数, 受験者数が求まったので，これらを使ってprop.test()関数を実行する．

# それぞれの成功数 (Success) と試行数 (Total) をベクトルにする
success <- c(27, 19) 
total <- c(55, 45)
# 確認のためそれぞれの合格率を算出
print(success/total)

## [1] 0.4909091 0.4222222

# Z検定
prop.test(success,total, correct=FALSE)

## 
##  2-sample test for equality of proportions without continuity correction
## 
## data:  success out of total
## X-squared = 0.47008, df = 1, p-value = 0.493
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.1269670  0.2643407
## sample estimates:
##    prop 1    prop 2 
## 0.4909091 0.4222222

結果として，男性と女性で合格率には有意な差はない，という結果となった．

なお有意確率を見ると，実は先の独立性の検定でのchisq.testの出力結果と同じ値となっている．これは，特に今回のように$2\times 2$のクロス集計が可能である場合，比率の差に有意な差があるかないかと，要因間の関連性があるかどうかという検定が，数学的には同じであるためである．

一方，3つ以上の水準を持つ要因を分析対象としている場合（すなわち，$3\times 2$以上のクロス集計となる場合）には独立性の検定と比率の差の検定は意味合いも，数学的にも，処理の面でも異なるものとなる．独立性の検定においては，先にも述べた通り，水準すべてを考慮したうえで，全体として独立か，関連性があるのかを検定することになる一方，比率の差の検定の場合は，「差」というものはそもそも2つのものの間での比較なので，3つ以上の水準を含んだ要因を分析対象にしていた場合でも，結局，その要因の中のどの水準とどの水準を比べるのか，という話になる．

例えば，ある要因でA, B, Cという群分けをして，合格率を比べる，という場合，全体として群分けと合格率に関係があるか，ということを調べたいのであれば，独立性の検定であり，こちらで述べた通りにカイ2乗検定を掛ける．一方で，各群の間での合格率の差が有意であるかどうかを調べたい，という場合にはA-B，B-C，C-Aという3つの組み合わせのそれぞれでZ検定を行う，ということになる（なお，この場合，次に述べる「多重性問題」への対応が必要になる）．

数値をコードで入力する方法

上記の例のように，直接数値を自分で手入力する方法でも全く問題はないが，手入力の場合，入力時に数値を誤ってしまう可能性がある．そうしたヒューマンエラーを避けるために，以下のようにコードで計算させる方法もある．

cross_table<-table(mydata_chisq$性別, mydata_chisq$合否)
print(cross_table)

##       
##        合格 不合格
##   女性   27     28
##   男性   19     26

success <- c(cross_table[1,1], cross_table[2,1]) # [行, 列]で指定する 
total <- c(cross_table[1,1]+cross_table[1,2], cross_table[2,1]+cross_table[2,2])
# 確認のためそれぞれの合格率を算出
print(success/total)

## [1] 0.4909091 0.4222222

# Z検定
prop.test(success,total)

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  success out of total
## X-squared = 0.23423, df = 1, p-value = 0.6284
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.1471690  0.2845427
## sample estimates:
##    prop 1    prop 2 
## 0.4909091 0.4222222

ある学部の就職支援担当の教員は，今年度の卒業予定者の9月末時点での内定率が例年より高いように感じた．今年度は卒業予定者153人中，内定を得ている学生は141名である．一方，過去10年間の状況を通算すると同時点に内定を得ていたのはのべ1732名中，1492名であった．単純に比較するならば，今年度は約92%の内定率であるのに対し，過去10年間の平均内定率は約86%である．この内定率の差は有意であるかどうかを検証せよ．（ヒント：この場合はすでに成功数と試行総数が端的に文章に与えられている！）

5 多重比較

例えば以下のようなケースではどのような分析をすればよいだろうか．

あるオンライン小売業者が，新しい広告キャンペーンを展開する際，3つの異なる広告バナーを作成した．それぞれの広告バナーには異なるキャッチコピーとデザインが使用されている．マーケティングチームは，どの広告がクリック率（広告が表示された回数に対する広告がクリックされた回数）を最も向上させるのかを調べるため，3つの広告キャンペーンを同時に2か月間実施し，毎日のA，B，Cのそれぞれの広告バナーのクリック率を調べた．その結果，次のようなデータが得られた．このデータをもとに，どの広告バナーが最もクリック率を向上させるのかを検証せよ（リンク）．

2つのバナーでの比較だった場合には，単純にt検定を実施すれば良いが，今回の例の場合には3つの比較となっている．このように3つ以上のデータで，各データ間の平均の比較を行うような場合には，どのようにすればよいだろうか？

5.1 検定の多重性問題

こういうケースでシンプルに浮かぶ方法としては，「AとB」「BとC」「CとA」という形で，1対1の組み合わせを作り，それぞれの組み合わせで平均値の比較を行う（このような比較のことを一対比較という）というものである．

例えば，平均値の値としてはA＜B＜Cとなっていたとして，上記の3つの組み合わせで，それぞれで差が有意であるかどうかを一対比較のt検定で調べた結果，「BとCの間には有意差はないが，AとB，CとAの間には有意差がある」という結果となった場合，「BとCの間にはそれほどの差はないが，AはB，Cに比べて広告としては弱い」という結論を出すことができる．

ただ，ここで注意しないといけないのは，検定の多重性という問題である

5.1.1 そもそもP値とは

検定では一般に5%を基準として，P値がこれを下回れば「有意である」と判定し，逆に上回っていた時には「測定時に生じた単なる誤差の可能性を否定できない」と判断する．

このP値の意味や見方はこちらで説明しているが，要するに「帰無仮説が真である確率」のことであり，言い換えれば「有意であるという判定結果が誤っている確率」である．つまり，P値が5％であった場合には，本当は「単なる誤差」であるのに「有意である」と誤って判定する確率が5%ある，という意味である（このことから「有意水準」という言葉の変わりに「危険率」という言葉が用いられることもある）．

このことから，「有意水準を下回っていれば有意とする」ということについて，なぜ「下回っていれば」なのかは，要するに「誤っている確率が5%未満であれば，誤る可能性は低いわけだから，結果を信頼してもよいだろう」と判断しているということである．ちなみに，「5%」を「十分に低い」の基準としているのは特段の理由があってのことではなく，統計学において経験上，設定された基準でしかない．

5.1.2 検定を多重化すると誤り率が上る！？

ある事柄Aについて，5%の有意水準で検定をしたときには，仮に「有意だ」と判定しても，それは5%の確率で誤っているということである．逆に言えば，95%の確率で正しいということである．

一方，5%水準で別の事柄Bを検定して「有意だ」と判定した場合，同じようにその判定結果は95%の確率で正しい．

では，「AとBの両方が正しい結果である」確率はいくらであろうか？

Aが正しい確率：95% (0.95)

Bが正しい確率：95% (0.95)

$\rarrow$ 両方がそろって正しい確率：0.95×0.95＝0.9025 (90.25%) となる．

つまり，AとBがそれぞれ個別には「95%の確率で正しい」となっていたとしても，AとB両方が同時に正しい結果である確率は90.25%となり，逆に「AとBのいずれか一方が間違っている」，もしくは「両方間違っている」確率は約10%になるということである．

より正確には，以下のような表となる．

	Bの判定は正しい	Bの判定は誤り
Aの判定は正しい	0.95×0.95＝0.9025	0.95×0.05＝0.0475
Aの判定は誤り	0.05×0.95＝0.0475	0.05×0.05＝0.025

この表から，AとBの両方が正しい確率は0.9025となる一方，AとBの少なくともどちらか一方が誤っている確率は0.0475×2＋0.025＝0.0975となる．

当然ながら判定は「両方正しい結果」であってほしい．しかしながら，個別に5%水準で判定をして，それを組み合わせるということをすると，「両方ともに正しい結果である確率」は90%ほどとなってしまう．

これは，さらにCという事柄についての検定を行うとすると，A，B，Cの3つがそろって正しい結果である確率は，

95%×95%×95%＝85.73%

同様に事柄の検定の数を増やすと，

4つ：95% × 95% × 95% × 95%＝81.45%
5つ：95% × 95% × 95% × 95% × 95%＝77.38%
6つ：95% × 95% × 95% × 95% × 95%×95%＝73.51%

となってしまう．

このように，検定を多重に行ってしまうと，個々の検定の結果が正しい確率が95%（誤り率が5%）でも，多重に行った検定が全て正しい確率は95%を下回ってしまう．これが検定の多重性問題である．

5.1.3 検定の多重性問題の解決方法

「検定の多重性問題を解決する」というのは，つまりは多重に行った検定の「全てが正しい」となる確率が95%を下回らないようにすることである．

そのための方法として，もっとも単純で簡単な方法は，個別の検定を行う際の有意水準である「5%」という値を検定を行う回数で割った値を個別の検定の有意水準とすることである．

より具体的には，例えば，

2回の検定を行うのであれば，5%÷2＝2.5%を個々の検定の有意水準とする．
3回の検定を行うのであれば，5%÷3≒1.67%を個々の検定の有意水準とする．（端数は切り上げ）
4回の検定を行うのであれば，5%÷4＝1.25%を個々の検定の有意水準とする．
5回の検定を行うのであれば，5%÷5＝1%を個々の検定の有意水準とする．
6回の検定を行うのであれば，5%÷6≒0.84%を個々の検定の有意水準とする．（端数は切り上げ）

ということである．

このように補正すると，多重に行った検定が全て正しい確率は以下の通りとなる．

検定回数	全て正しい確率
2回の場合	97.5% × 97.5％＝95.06%
3回の場合	98.33% × 98.33% × 98.33%＝95.07%
4回の場合	98.75% × 98.75% × 98.75% × 98.75%＝95.09%
5回の場合	99% × 99% × 99% × 99% × 99%＝95.10%
6回の場合	99.16% × 99.16% × 99.16% × 99.16% × 99.16% × 99.16%＝95.06%

このように，多重検定を行う際に実際に行う検定の数で有意水準を割った値を個々の有意性評価の基準とする方法をボンフェローニ（Bonferroni）補正と呼ぶ．

実際には有意水準を補正するのではなく，確率であるP値の方を補正し（検定の回数を掛ける），補正された有意確率が5%を上回るか下回るかで評価することが多い．すなわち，

2回の検定であれば，一対比較時のP値を2倍する
3回の検定であれば，一対比較時のP値を3倍する
4回の検定であれば，一対比較時のP値を4倍する
5回の検定であれば，一対比較時のP値を5倍する

・・・

とする．

5.1.4 補正の方法

補正方法には他にもいくつかの方法があるが，ボンフェローニ補正は最もシンプルで直感的な方法であるため，手計算を行う際にはよく使われる．

ただ，ボンフェローニ補正は，検定の数が増えるにつれて有意水準が厳しくなるという性質があるため，検定の数が多い場合には，ホルム法やシダック法，ダンスカー法などの補正方法を用いることもある．

以下で紹介するRでの多重比較の補正方法は，ホルム補正がデフォルトとなっている．

5.2 多重比較時のP値の補正方法

Rでは，p.adjust()という関数を用いて，上記のようなP値の補正を行うことができる．引数は以下の通りである．

第1引数：多重比較時の個々の一対比較の結果のP値をまとめたものベクトルで与える．
第2引数（省略可）：補正方法を指定する．デフォルトではholm補正を行うように設定されており，省略した場合にはholm補正となる．上記のボンフェローニ補正を行う場合には，method="bonferroni"と指定する．

出力結果は，補正されたP値を返す．この補正後のP値が5%を下回るかで，その組み合わせが有意であるかどうかを判断する．

5.2.1 平均の多重比較

先の例題では以下のようになる．

mydata_multcomp <- read.csv(
  "./practice/example_multcomp.csv", 
  header=T, 
  fileEncoding="UTF8"
) # データ読み込み

# データの確認
head(mydata_multcomp)

##         Date      A      B      C
## 1 2024-07-01 0.0772 0.0461 0.0552
## 2 2024-07-02 0.1040 0.0446 0.0580
## 3 2024-07-03 0.0650 0.0478 0.0594
## 4 2024-07-04 0.0685 0.0509 0.0606
## 5 2024-07-05 0.0747 0.0525 0.0505
## 6 2024-07-06 0.0709 0.0520 0.0583

# AとBの比較
res_A_B<-t.test(mydata_multcomp$A, mydata_multcomp$B)
print(res_A_B)

## 
##  Welch Two Sample t-test
## 
## data:  mydata_multcomp$A and mydata_multcomp$B
## t = 11.382, df = 63.999, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.02378236 0.03390797
## sample estimates:
##  mean of x  mean of y 
## 0.07918226 0.05033710

# BとCの比較
res_B_C<-t.test(mydata_multcomp$B, mydata_multcomp$C)
print(res_B_C)

## 
##  Welch Two Sample t-test
## 
## data:  mydata_multcomp$B and mydata_multcomp$C
## t = -14.834, df = 110.62, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.011330410 -0.008659913
## sample estimates:
##  mean of x  mean of y 
## 0.05033710 0.06033226

# CとAの比較
res_C_A<-t.test(mydata_multcomp$C, mydata_multcomp$A)
print(res_C_A)

## 
##  Welch Two Sample t-test
## 
## data:  mydata_multcomp$C and mydata_multcomp$A
## t = -7.3551, df = 66.822, p-value = 3.541e-10
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.02396572 -0.01373428
## sample estimates:
##  mean of x  mean of y 
## 0.06033226 0.07918226

# 各P値をベクトルにまとめる．
# それぞれのP値は$p.valueで取得できる
# わかりやすさのため，名前付きベクトルにしている．
P_values <- c(
  AとB=res_A_B$p.value, 
  BとC=res_B_C$p.value, 
  CとA=res_C_A$p.value
  )

print(P_values)

##         AとB         BとC         CとA 
## 5.030179e-17 4.600584e-28 3.541057e-10

# 多重比較のP値　補正方法を指定していないのでホルム補正が用いられる
p.adjust(P_values)

##         AとB         BとC         CとA 
## 1.006036e-16 1.380175e-27 3.541057e-10

# 多重比較のP値のボンフェローニ補正
p.adjust(P_values, method="bonferroni")

##         AとB         BとC         CとA 
## 1.509054e-16 1.380175e-27 1.062317e-09

これらの出力結果から，今回の場合だと補正後でも，いずれの組み合わせでも有意な差があることが分かった．

このp.adjust()を用いた補正方法は，平均の比較，分散の比較，比率の比較いずれの比較においても用いることができる．

5.2.2 分散の多重比較

分散の多重比較の例を示す．

ある製造工程において、製品の重量を測定した結果、以下のようなデータが得られた．このデータをもとに、各工程間の重量のばらつきに有意な差があるかどうかを検証せよ（リンク）．

この場合，A-B，B-C，C-Aでそれぞれ分散の一対比較を行った後，そのP値を補正することになる．

mydata_multcomp_var <- read.csv(
  "./practice/example_multcomp_var.csv",
  header=T,
  fileEncoding="UTF8"
) # データ読み込み

# データの確認
head(mydata_multcomp_var)

##       工程A    工程B     工程C
## 1 106.85479 108.4668 106.02999
## 2  97.17651 101.0935 103.19471
## 3 101.81564 109.1404 108.79082
## 4 103.16431 116.0734 101.36648
## 5 102.02134 119.4760  98.15859
## 6  99.46938 107.8477 107.16409

#分散の一対比較
res_AB <- var.test(mydata_multcomp_var$工程A, mydata_multcomp_var$工程B)
res_BC <- var.test(mydata_multcomp_var$工程B, mydata_multcomp_var$工程C)
res_CA <- var.test(mydata_multcomp_var$工程C, mydata_multcomp_var$工程A)

#P値をベクトルにまとめる
P_values_var <- c(
  AとB=res_AB$p.value, 
  BとC=res_BC$p.value, 
  CとA=res_CA$p.value
)
print(P_values_var) #P値の確認

##      AとB      BとC      CとA 
## 0.4709159 0.7281850 0.2868715

#P値の補正　補正方法は指定していないのでホルム補正が用いられる
p.adjust(P_values_var)

##      AとB      BとC      CとA 
## 0.9418318 0.9418318 0.8606146

#P値の補正　補正方法はボンフェローニ補正
p.adjust(P_values_var, method = "bonferroni")

##      AとB      BとC      CとA 
## 1.0000000 1.0000000 0.8606146

# 参考：補正前のP値
p.adjust(P_values_var, method = "none")

##      AとB      BとC      CとA 
## 0.4709159 0.7281850 0.2868715

この結果から、補正を行った結果、全ての組み合わせでP値が0.05を上回っていることが分かる．

5.2.3 比率の多重比較

比率の多重比較の例を示す．

ある製造工程A，B, Cにおいて、製品の不良品率を測定したするためにランダムに1000個のサンプルを取り出し、不良品かどうかを調べた（良・不良）．その結果、以下のようなデータが得られた．このデータをもとに、各工程間の不良品率に有意な差があるかどうかを検証せよ（リンク）．

この場合もこれまでの多重比較と同様に，A-B，B-C，C-Aでそれぞれ比率の検定を一対比較で行った後，そのP値を補正することになる．

実際に行ってみよう．まずは各工程の不良の数を取得する必要がある．

mydata_multcomp_prop <- read.csv(
  "./practice/example_multcomp_prop.csv",
  header=T,
  fileEncoding="UTF8"
) # データ読み込み

# データの確認
mydata_multcomp_prop$工程A <- as.factor(mydata_multcomp_prop$工程A)
mydata_multcomp_prop$工程B <- as.factor(mydata_multcomp_prop$工程B)
mydata_multcomp_prop$工程C <- as.factor(mydata_multcomp_prop$工程C)
head(mydata_multcomp_prop)

##   工程A 工程B 工程C
## 1    良    良  不良
## 2    良    良    良
## 3    良    良    良
## 4    良    良  不良
## 5    良    良    良
## 6    良    良    良

summary(mydata_multcomp_prop)

##   工程A      工程B      工程C    
##  不良: 43   不良:104   不良:121  
##  良  :957   良  :896   良  :879

summaryの結果から，各工程の不良数は，Aが43，Bが104, Cが121であることがわかった．総数は課題の設定から1000であることがわかっているので，これらを用いて，まずは一対比較での比率の検定を行う．後でprop.testの出力結果からp値を取り出したいので，結果はオブジェクトに保存しておく．

#比率の検定
res_AB <- prop.test(c(43, 104), c(1000, 1000))
res_BC <- prop.test(c(104, 121), c(1000, 1000))
res_CA <- prop.test(c(121, 43), c(1000, 1000))

これらの結果からp値を取り出して，p.adjust()関数を用いて補正を行う．

#P値をベクトルにまとめる
P_values_prop <- c(
  AとB=res_AB$p.value, 
  BとC=res_BC$p.value, 
  CとA=res_CA$p.value
)
print(P_values_prop) #P値の確認

##         AとB         BとC         CとA 
## 2.728969e-07 2.575269e-01 3.485408e-10

#P値の多重比較のための補正　補正方法：ホルム補正（デフォルト値）
p.adjust(P_values_prop)

##         AとB         BとC         CとA 
## 5.457937e-07 2.575269e-01 1.045622e-09

結果として，BとCの間には有意な差はないが，AとB，CとAの間には不良品率に有意な差がある，という結果となった．

なお，今回は工程A, B, Cはワイド形式のデータとして与えられたので，summary()関数で各工程の不良数を取得することができたが，ロング形式でデータが与えられている場合には，table()関数を用いてクロス集計表の作成することによって各工程の不良数を取得することができる．以下では，ロング形式に変換したうえで，table()関数を用いてクロス集計表を作成する例を示す．

## ロング形式に変換
library(tidyverse)
mydata_multcomp_prop_long <- pivot_longer(
  mydata_multcomp_prop,
  cols=colnames(mydata_multcomp_prop),
  names_to = "工程",
  values_to= "品質" #測定変数に与える列名
)
mydata_multcomp_prop_long$工程 <- as.factor(mydata_multcomp_prop_long$工程)
mydata_multcomp_prop_long$品質 <- as.factor(mydata_multcomp_prop_long$品質)

# クロス集計表の作成
cross_table <- table(mydata_multcomp_prop_long$工程, mydata_multcomp_prop_long$品質)
print(cross_table)

##        
##         不良  良
##   工程A   43 957
##   工程B  104 896
##   工程C  121 879

5.3 平均の多重比較を行う別の方法

さて、先ほどの例では、A-B、B-C、C-Aの3つの組み合わせでそれぞれの平均値の差を検定を行った後、p.adjust()関数を用いてそのP値を補正する、という方法をとったが，Rでは平均の多重比較に関してはpairwise.t.test()関数を用いることで、個別の検定を行うことなしに多重比較を行うことができる．

pairwise.t.test()関数の引数は以下の通りである．

第1引数 : 比較をしたい変数
第2引数 : グループ分けに用いる属性変数
第3引数pool.sd : FALSEを与える．デフォルトではTRUEとなっており，比較する各データの等分散性が仮定されているが，こちらで述べている通り，最近では等分散を仮定しないのが一般的となっているためFALSEにする．
第4引数p.adjust.method (省略可) : 補正方法を指定する．省略した場合，デフォルトとして設定されているholm補正が用いられる．

ただし，注意しなければならない点として、pairwise.t.test()で多重比較を行う場合にはデータはロング形式でなければならない．今回の例題のデータはワイド形式になっているので，以下の例ではまずはロング形式に変換している．

mydata_multcomp <- read.csv(
  "./practice/example_multcomp.csv", 
  header=T, 
  fileEncoding="UTF8"
) # データ読み込み

#ワイド形式からロング形式に変換
mydata_multcomp_long <- pivot_longer(
  mydata_multcomp,
  cols=colnames(mydata_multcomp)[2:4],
  names_to = "広告",
  values_to= "クリック率" #測定変数に与える列名
)

#平均と標準偏差の確認
aggregate(クリック率~広告, data=mydata_multcomp_long, FUN=mean)

##   広告 クリック率
## 1    A 0.07918226
## 2    B 0.05033710
## 3    C 0.06033226

aggregate(クリック率~広告, data=mydata_multcomp_long, FUN=sd)

##   広告  クリック率
## 1    A 0.019713926
## 2    B 0.003091990
## 3    C 0.004311471

#多重比較の実施　デフォルトではホルム補正
pairwise.t.test(mydata_multcomp_long$クリック率, mydata_multcomp_long$広告, pool.sd = FALSE)

## 
##  Pairwise comparisons using t tests with non-pooled SD 
## 
## data:  mydata_multcomp_long$クリック率 and mydata_multcomp_long$広告 
## 
##   A       B      
## B < 2e-16 -      
## C 3.5e-10 < 2e-16
## 
## P value adjustment method: holm

#多重比較の実施　ボンフェローニ補正
pairwise.t.test(mydata_multcomp_long$クリック率, mydata_multcomp_long$広告, pool.sd = FALSE, p.adjust.method="bonferroni")

## 
##  Pairwise comparisons using t tests with non-pooled SD 
## 
## data:  mydata_multcomp_long$クリック率 and mydata_multcomp_long$広告 
## 
##   A       B      
## B < 2e-16 -      
## C 1.1e-09 < 2e-16
## 
## P value adjustment method: bonferroni

結果の見方として，A-Bの比較のP値（補正済み）が< 2e-16となっている．これは，「$2.0 \times 10^{-16}$よりもさらに小さい」を意味しており，ほぼ0と見なせる値である．この結果から，AとBの間には有意な差があると結論づけることができる．同様に，A-Cの比較のP値（補正済み）は$1.1\times 10^{-9}$であり，有意な差であるといえる．さらに，B-Cについても< 2e-16であり，有意な差があるといえる．

これら結果から，広告バナーはAが最もクリック率を向上させる効果があり，次いでCとなり，Bがもっとも効果が低いと結論づけることができる．

なお，pairwise.t.test()関数はデフォルトではホルム補正を行う．ボンフェローニ補正を行いたい場合には，p.adjust.method引数に"bonferroni"を指定することで行うことができる．また補正前の値を確認したければ"none"を指定するとよい．

なお，補正の結果，P値が1を超えてしまうケースが出てくるが，その場合には「1」と表記される．

5.3.1 対応ありデータでの多重比較

対応ありデータでの多重比較を行う場合には，pairwise.t.test()関数の引数にpaired=TRUEのオプションを追加する．

あるパンメーカでは新しい食パンの開発に取り組んでおり，小麦粉の種類や使うバターの量などを変えた4種類のパンが市販化の最終候補として残った．そこで，この4つのパンのいずれが消費者に好まれるのかを調査するため，20人の消費者にモニターとして集まってもらい，それぞれのパンを試食してもらい，風味，食感，味わい，色味の4つの項目にそれぞれ5段階評価で得点をつけてもらった（データ）．その合計点から，いずれの食パンが最も好まれるのかを検証せよ．

この例題の場合，各モニターがA～Dの全てのパンに対して評価を行っており，対応関係のあるデータとなる．そのため，多重比較を行う際にはpaired=TRUEを指定する必要がある．

mydata_multcomp_paired <- read.csv(
  "./practice/example_multcomp_paired.csv",
  header=T,
  fileEncoding="UTF8"
) # データ読み込み

# データの確認
head(mydata_multcomp_paired)

##   協力者No.  パン 風味 味わい 食感 色味
## 1         1 パンA    4      3    4    3
## 2         1 パンB    4      4    3    3
## 3         1 パンC    4      4    3    4
## 4         1 パンD    3      4    4    3
## 5         2 パンA    3      3    3    4
## 6         2 パンB    3      4    3    3

# 総合評価の算出
mydata_multcomp_paired$総合評価 <- mydata_multcomp_paired$風味 + 
  mydata_multcomp_paired$味わい + 
  mydata_multcomp_paired$食感 + 
  mydata_multcomp_paired$色味

# 平均と標準偏差の確認
aggregate(総合評価~パン, mydata_multcomp_paired, mean)

##    パン 総合評価
## 1 パンA    14.05
## 2 パンB    13.60
## 3 パンC    15.40
## 4 パンD    12.85

aggregate(総合評価~パン, mydata_multcomp_paired, sd)

##    パン  総合評価
## 1 パンA 1.4317821
## 2 パンB 1.3138934
## 3 パンC 0.9947229
## 4 パンD 1.3484884

#多重比較の実施 (ホルム補正)
pairwise.t.test( mydata_multcomp_paired$総合評価,mydata_multcomp_paired$パン, pool.sd=FALSE, paired=TRUE)

## 
##  Pairwise comparisons using paired t tests 
## 
## data:  mydata_multcomp_paired$総合評価 and mydata_multcomp_paired$パン 
## 
##       パンA  パンB  パンC  
## パンB 0.2513 -      -      
## パンC 0.0342 0.0023 -      
## パンD 0.0766 0.2100 7.4e-06
## 
## P value adjustment method: holm

多重比較の結果を報告するときには，各組み合わせの補正済みのP値を示すことが一般的である．たとえば，上記の例の場合，「パンCはパンA，パンB，パンDのいずれとも有意な差がある（それぞれ有意確率は順に0.034, 0.002, 7.4$\times 10^{-6}$．いずれもholm補正済みの値）」と記述する．

また，平均の比較なのでグラフもあわせて示すことが多い．以下は上記の結果をグラフ表記したものである．

# 棒グラフをつくるために平均と標準偏差を収めたデータフレームを作成
平均=aggregate(総合評価~パン, mydata_multcomp_paired, mean)
標準偏差=aggregate(総合評価~パン, mydata_multcomp_paired, sd)

df <- data.frame(
  パン=平均$パン,
  平均=平均$総合評価, 
  標準偏差=標準偏差$総合評価)

# 内容確認
print(df)

##    パン  平均  標準偏差
## 1 パンA 14.05 1.4317821
## 2 パンB 13.60 1.3138934
## 3 パンC 15.40 0.9947229
## 4 パンD 12.85 1.3484884

# グラフ表示
g <- ggplot(df, aes(x=パン, y=平均, fill=パン)) +
  geom_bar(stat="identity") +
  geom_errorbar(aes(ymin=平均-標準偏差, ymax=平均+標準偏差), width=.1) +
  coord_cartesian(ylim = c(0, 22)) +
  labs(title="各パンの総合評価", x="パン", y="総合評価") +
  # A-Cの結果
  annotate("segment", x = 1, xend = 1, y = 17, yend = 20) +
  annotate("segment", x = 1, xend = 3, y = 20, yend = 20) +
  annotate("segment", x = 3, xend = 3, y = 20, yend = 17) +
  annotate("text", x = 2, y = 21, label = "p=0.0342, *") +
  # B-Cの結果
  annotate("segment", x = 2, xend = 2, y = 17, yend = 19.5) +
  annotate("segment", x = 2, xend = 2.9, y = 19.5, yend = 19.5) +
  annotate("segment", x = 2.9, xend = 2.9, y = 19.5, yend = 17) +
  annotate("text", x = 2.5, y = 18.5, label = "p=0.0023, **") +
  # C-Dの結果
  annotate("segment", x = 4, xend = 4, y = 17, yend = 20) +
  annotate("segment", x = 4, xend = 3.1, y = 20, yend = 20) +
  annotate("segment", x = 3.1, xend = 3.1, y = 20, yend = 17) +
  annotate("text", x = 3.5, y = 21, label = "p=7.4e-06, ***") 


print(g)

行ごとの和や平均の算出

データフレームに含まれている複数の変数について，サンプルごとの合計や平均を算出する場合には，もっとも簡単なのはそれぞれの変数を$で指定して計算させる方法である（上記のその通りに行っている）．ただ，この方法だと変数が多い場合には書く手間もかかるし，コードも長くなってしまって見づらくなる．

そこで，より短いコードで簡潔に書けるコードとして以下のような書き方がある．

rowSums, rowMeans関数を用いる方法
apply関数を用いる方法
dplyrパッケージ（tidyverseにふくまれている）のmutate()関数を用いる方法．

それぞれ以下で例を示すが，結論としてはmutate関数を用いる方法が最も簡潔でわかりやすいのでおすすめである．

rowSums, rowMeans関数を用いる方法

計算させたい列を指定する際には，c(3:6)のように列番号を指定するか，c("風味","味わい","食感","色味")のように列名を指定する．

mydata_multcomp_paired$総合評価_和 <- rowSums(mydata_multcomp_paired[c(3:6)]) # 合計
head(mydata_multcomp_paired)

##   協力者No.  パン 風味 味わい 食感 色味 総合評価 総合評価_和
## 1         1 パンA    4      3    4    3       14          14
## 2         1 パンB    4      4    3    3       14          14
## 3         1 パンC    4      4    3    4       15          15
## 4         1 パンD    3      4    4    3       14          14
## 5         2 パンA    3      3    3    4       13          13
## 6         2 パンB    3      4    3    3       13          13

mydata_multcomp_paired$総合評価_平均 <- rowSums(mydata_multcomp_paired[c("風味","味わい","食感","色味")]) # 合計
head(mydata_multcomp_paired)

##   協力者No.  パン 風味 味わい 食感 色味 総合評価 総合評価_和 総合評価_平均
## 1         1 パンA    4      3    4    3       14          14            14
## 2         1 パンB    4      4    3    3       14          14            14
## 3         1 パンC    4      4    3    4       15          15            15
## 4         1 パンD    3      4    4    3       14          14            14
## 5         2 パンA    3      3    3    4       13          13            13
## 6         2 パンB    3      4    3    3       13          13            13

mydata_multcomp_paired$総合評価_和 <- rowMeans(mydata_multcomp_paired[c(3:6)]) # 平均
head(mydata_multcomp_paired)

##   協力者No.  パン 風味 味わい 食感 色味 総合評価 総合評価_和 総合評価_平均
## 1         1 パンA    4      3    4    3       14        3.50            14
## 2         1 パンB    4      4    3    3       14        3.50            14
## 3         1 パンC    4      4    3    4       15        3.75            15
## 4         1 パンD    3      4    4    3       14        3.50            14
## 5         2 パンA    3      3    3    4       13        3.25            13
## 6         2 パンB    3      4    3    3       13        3.25            13

mydata_multcomp_paired$総合評価_平均 <- rowMeans(mydata_multcomp_paired[c("風味","味わい","食感","色味")]) # 平均
head(mydata_multcomp_paired)

##   協力者No.  パン 風味 味わい 食感 色味 総合評価 総合評価_和 総合評価_平均
## 1         1 パンA    4      3    4    3       14        3.50          3.50
## 2         1 パンB    4      4    3    3       14        3.50          3.50
## 3         1 パンC    4      4    3    4       15        3.75          3.75
## 4         1 パンD    3      4    4    3       14        3.50          3.50
## 5         2 パンA    3      3    3    4       13        3.25          3.25
## 6         2 パンB    3      4    3    3       13        3.25          3.25

apply関数を用いる方法

apply関数は，行列やデータフレームの各行や各列に対して指定した関数を適用するための関数である．引数は以下の通り3つある．

第1引数：処理対象の行列やデータフレーム
第2引数：処理を行う方向（1:行方向，2:列方向）
第3引数：適用する関数．例えば平均であればmean, 合計であればsum，標準偏差であればsdなど．

今回は行ごとに算出させたいので，第2引数は1となる．第1引数に与えるデータはrowMeansやrowSumsと同様に，列番号を指定するか，列名を指定するかで指定する．

mydata_multcomp_paired$総合評価_和 <- apply(mydata_multcomp_paired[c("風味","味わい","食感","色味")], 1, sum) # 合計
head(mydata_multcomp_paired)

##   協力者No.  パン 風味 味わい 食感 色味 総合評価 総合評価_和 総合評価_平均
## 1         1 パンA    4      3    4    3       14          14          3.50
## 2         1 パンB    4      4    3    3       14          14          3.50
## 3         1 パンC    4      4    3    4       15          15          3.75
## 4         1 パンD    3      4    4    3       14          14          3.50
## 5         2 パンA    3      3    3    4       13          13          3.25
## 6         2 パンB    3      4    3    3       13          13          3.25

mydata_multcomp_paired$総合評価_平均 <- apply(mydata_multcomp_paired[c("風味","味わい","食感","色味")], 1, mean) # 平均
head(mydata_multcomp_paired)

##   協力者No.  パン 風味 味わい 食感 色味 総合評価 総合評価_和 総合評価_平均
## 1         1 パンA    4      3    4    3       14          14          3.50
## 2         1 パンB    4      4    3    3       14          14          3.50
## 3         1 パンC    4      4    3    4       15          15          3.75
## 4         1 パンD    3      4    4    3       14          14          3.50
## 5         2 パンA    3      3    3    4       13          13          3.25
## 6         2 パンB    3      4    3    3       13          13          3.25

dplyrパッケージのmutate()関数を用いる方法

mutate()関数は，データフレームに含まれている変数を使って，新しい列を元のデータフレームに追加する関数である．引数は以下の通り．

第1引数：データフレーム
第2引数以降：新しく追加する列名と，その値をデータフレームに含まれる変数を使って算出する式．複数の列を一度に追加することも可能であり，その場合，同じ要領で順に第3引数，第4引数，第5引数…というように追加していく．

返り値は，元のデータフレームに新しい列が追加されたデータフレームである．したがって，その返り値を受けるオブジェクトは，rowMeansやrowSums，applyのように$を用いて特定の列を指定するのではなく，あくまでデータフレームでなければならない

以下の例では，総合評価_和と総合評価_平均という2つの列を一度に追加している．

なお，dplyrパッケージはtidyverseにふくまれているのでtidyverseを読み込んでいるならば，改めてインストールしたり読み込んだりする必要はない．

library(tidyverse)
mydata_multcomp_paired <- mutate(mydata_multcomp_paired, 
                                 総合評価_和 = 風味 + 味わい + 食感 + 色味,
                                 総合評価_平均= (風味 + 味わい + 食感 + 色味)/4
                                 )
head(mydata_multcomp_paired)

##   協力者No.  パン 風味 味わい 食感 色味 総合評価 総合評価_和 総合評価_平均
## 1         1 パンA    4      3    4    3       14          14          3.50
## 2         1 パンB    4      4    3    3       14          14          3.50
## 3         1 パンC    4      4    3    4       15          15          3.75
## 4         1 パンD    3      4    4    3       14          14          3.50
## 5         2 パンA    3      3    3    4       13          13          3.25
## 6         2 パンB    3      4    3    3       13          13          3.25

rowSumsやrowMeansはその名の通り，行毎に和や平均を算出する関数であり，最もシンプルである．ただ，あくまで和と平均しか算出できない．

apply関数の場合，所定の組み込み関数（meanやsum，sdなど）だけでなく，自分で定義した関数を適用することも可能であるため，汎用性という点ではapply関数が最も優れているだろう．ただ，もともとapply関数は行方向にも列方向にも適用できる関数であり，いずれの方向の計算をさせるかを第2引数で指定する（行方向:1，列方向:2）必要があるため，使い方が少し複雑である．また，自分で関数を定義するのは単純に統計分析に用いるよりも1段上のプログラミングの知識が必要になる．

mutate()関数は，applyに比べて直感的でわかりやすく，かつ複数の変数を一度に作ることもできるため，使いやすさの面でapplyよりも優れているので，データフレームの操作においてはmutate()関数を使うことが多い．

全国チェーンのあるスーパーが，割引クーポン、ポイント還元、限定セールのいずれが売り上げ向上に効果的かを検証するための社会実験を行いました．市場規模や顧客特性がほぼ同程度の3つの地域（A，B，C）の店舗を対象に，3日間限定で，A地域では割引クーポン、B地域ではポイント還元、C地域では限定セールを行った．それぞれの割引・還元率は同じとなってる．この3日間の各地域の各店舗の売上データ（各地域の店舗別の売り上げデータ）を用いて，以下の課題を実行せよ．

head()とsummary()を用いて与データがどういうデータかの確認を行え．
それぞれの地域ごとの売上の平均値と標準偏差を算出せよ．
A-B，B-C，C-Aの3つの組み合わせのそれぞれについて，一対比較の結果を示せ．
一対比較の結果を用いて，p.adjust()関数を用いて補正を行い，各組み合わせのP値を示せ（補正方法はholm補正でよい）．
pairwise.t.test()を用いた多重比較を行い，上の課題の結果と一致することを確認せよ（補正方法はholm補正でよい）．
この調査結果から導かれる結論を述べよ．

ある学生が卒業研究で，オンライン講義動画を視聴するときの視聴速度と内容の定着率に有意な差があるのかを調べたいと考えた．このことを検証するために，実験協力者36名全員に対して，1倍速， 1.5倍速, 2倍速の3つの速度で講義を視聴してもらい，それぞれ講義の内容に関するクイズに回答してもらってその成績で比較をする方法を取った．この実験の結果を示すデータがこちらである．このデータを用いて，以下の課題を実行せよ．

head()とsummary()を用いて与データがどういうデータかの確認を行え．
1倍速，1.5倍速，2倍速のそれぞれの視聴速度でのクイズの成績の平均値と標準偏差を算出せよ．
A-B，B-C，C-Aの3つの組み合わせのそれぞれについて，一対比較の結果を示せ．
一対比較の結果を用いて，p.adjust()関数を用いて補正を行い，各組み合わせのP値を示せ（補正方法はbonferroni補正を用いること）．
pairwise.t.test()を用いた多重比較を行い，上の課題の結果と一致することを確認せよ（補正方法は補正方法はbonferroni補正を用いること）．
この実験結果から導かれる結論を述べよ．

課題遂行には関係しないが，実際に上記のような実験を行う際に考慮しなければならない点について説明しているので，興味ああれば見てみてほしい．

参考　実験条件の組み合せ

このような実験を行う際には順序効果など，実験結果にバイアスを与える要因を相殺できるように実験を組み必要がある．

今回の例題の場合には，実験を行うにあたって，用いるオンライン講義が同じだと，視聴する毎に成績が良くなると考えられるので，用いるオンライン講義動画は内容の全く異なる3つ（A, B, Cとする）を用意する必要があるだろう．

ただそうすると，A，B, Cをどの順序でどの速度で視聴するのかも結果に影響を与えてしまう．この影響を相殺するためには，考えられるすべての組み合わせで実験を行うことである．

今回の例だと，1つ目に1倍速，2つ目に1.5倍速，3つ目に2倍速の動画を割り当てるケースを考えると，1つ目，2つ目，3つ目にA, B, Cどの動画を割り当てるかの割り当て方は6通りある．それがさらに，1つ目に1倍速，2つ目に2倍速，3つ目に1.5倍速の動画を割り当てるケースや，1つ目に1.5倍速，2つ目に1倍速，3つ目に2倍速の動画を割り当てるケースなど，倍速の割り当て方も6通りになる．各倍速の割り当て方に対して各動画の割り当て方が6通りあるので，全体で6x6=36通りの実験設計が必要となる．

具体的には以下の表のような組み合わせになる，ということである．

被験者	1つ目に視聴する動画	2つ目に視聴する動画	3つ目に視聴する動画
1	A・1倍速	B・1.5倍速	C・2倍速
2	A・1倍速	C・1.5倍速	B・2倍速
3	B・1倍速	A・1.5倍速	C・2倍速
4	B・1倍速	C・1.5倍速	A・2倍速
5	C・1倍速	A・1.5倍速	B・2倍速
6	C・1倍速	B・1.5倍速	A・2倍速
7	A・1倍速	B・2倍速	C・1.5倍速
8	A・1倍速	C・2倍速	B・1.5倍速
9	B・1倍速	A・2倍速	C・1.5倍速
10	B・1倍速	C・2倍速	A・1.5倍速
11	C・1倍速	A・2倍速	B・1.5倍速
12	C・1倍速	B・2倍速	A・1.5倍速
13	A・1.5倍速	B・1倍速	C・2倍速
14	A・1.5倍速	C・1倍速	B・2倍速
15	B・1.5倍速	A・1倍速	C・2倍速
16	B・1.5倍速	C・1倍速	A・2倍速
17	C・1.5倍速	A・1倍速	B・2倍速
18	C・1.5倍速	B・1倍速	A・2倍速
19	A・1.5倍速	B・2倍速	C・1倍速
20	A・1.5倍速	C・2倍速	B・1倍速
21	B・1.5倍速	A・2倍速	C・1倍速
22	B・1.5倍速	C・2倍速	A・1倍速
23	C・1.5倍速	A・2倍速	B・1倍速
24	C・1.5倍速	B・2倍速	A・1倍速
25	A・2倍速	B・1倍速	C・1.5倍速
26	A・2倍速	C・1倍速	B・1.5倍速
27	B・2倍速	A・1倍速	C・1.5倍速
28	B・2倍速	C・1倍速	A・1.5倍速
29	C・2倍速	A・1倍速	B・1.5倍速
30	C・2倍速	B・1倍速	A・1.5倍速
31	A・2倍速	B・1.5倍速	C・1倍速
32	A・2倍速	C・1.5倍速	B・1倍速
33	B・2倍速	A・1.5倍速	C・1倍速
34	B・2倍速	C・1.5倍速	A・1倍速
35	C・2倍速	A・1.5倍速	B・1倍速
36	C・2倍速	B・1.5倍速	A・1倍速

被験者内実験計画（一人の被験者が3つの実験条件（1倍速，1.5倍速，2倍速）をすべて実施する実験計画）の場合には，この表にある組み合わせに対して一人ずつ被験者を割り当てることによって，全体としてバイアスを除去した実験が可能となる．

ただ，このような実験計画は被験者が最低でも36人必要となる．一方で，一人があくまで1つの条件しか経験しない形（被験者間実験計画）の実験条件であれば，被験者が全体として同質なのであれば，バイアス要因の相殺は考慮に入れる必要はない（動画を3種類も用意しなくてもよいし，実験順序という概念も存在しない）．したがって，被験者内計画にするか，被験者間計画にするかは，その時々での実験を実施できる環境に応じて選択するとよい．

基本的な検定

藤野秀則,

初版日:2024-05-02 更新日:2025-02-02

1 平均の比較

1.1 比較の原理

1.2 t検定の方法（基本）

1.3 対応関係のあるデータでのt検定

1.4 ロング形式データでのt検定

1.4.1 式, data形式では対応ありデータの分析はできない！！

1.5 t検定の結果の報告方法

2 分散の比較

2.1 F検定の方法（基本）

2.2 ロング形式のデータでのF検定

2.3 F検定の結果の報告方法

3 独立性の検定

3.1 「独立」とは？

3.2 カイ2乗検定とフィッシャーの正確確率検定

3.2.1 カイ2乗検定の方法

3.2.2 Fisherの正確確率検定の方法

3.3 カイ2乗検定とFisherの正確確率検定の使い分け

3.4 結果の報告方法

3.5 3つ以上の水準を持つ要因での独立性の検定

4 比率の差の検定

4.1 Z検定の方法

5 多重比較

5.1 検定の多重性問題

5.1.1 そもそもP値とは

5.1.2 検定を多重化すると誤り率が上る！？

5.1.3 検定の多重性問題の解決方法

5.1.4 補正の方法

5.2 多重比較時のP値の補正方法

5.2.1 平均の多重比較

5.2.2 分散の多重比較

5.2.3 比率の多重比較

5.3 平均の多重比較を行う別の方法

5.3.1 対応ありデータでの多重比較