12. \(t\) 検定（2グループ）

• このセクションで使う R パッケージ一覧

library(tidyverse)
library(broom)
library(patchwork)
library(DT)
library(ggbeeswarm)
library(ggsignif)
library(rmarkdown)

1. \(t\) 検定

1.1 \(t\) 検定の目的

1) 標本平均が特定の値かどうか調べる

2) 二つの異なるグループの平均が異なるかどうか調べる

• 前セッションで紹介した中心極限定理によれば、母集団が正規分布していなくても、平均 \(mu\)，分散 \(σ^2\) の母集団からの無作為標本における標準誤差の分布は「標本サイズが十分に大きければ」近似的に正規分布に従う 
  
• しかし、標本サイズが十分に大きいケースばかりとは限らず、少ない標本サイズしか入手できないことが多々ある 
  
• この問題を解決したのが \(t\) 検定 を考え出した William Sealy Gosset という元ギネスビール社員 
  
• student's t という名前で論文を投稿したため、\(t\) 検定と呼ばれる 
  
• Gosset は「標本サイズが小さくても」母集団からの無作為標本における標準誤差の分布は \((n-1)\) の自由度で \(t\) 分布することを発見した 
• この発見のおかげで、小さい標本サイズの場合でも \(t\) 分布の特徴を使うことで統計的推計が可能となった 
• 標本をとり、それぞれの標本平均に明らかに違いがあったとする
  
• しかし、その違いは「標本における」違いに過ぎず、その結果を一般化することはできない
  
• ここでは、架空のデータを R 上で作り出して、Paired と Unpaired データの二種類のデータに関して \(t\) 検定の方法を演習する

1.2 \(t\) 検定結果の表記方法

• 一般的に \(t\) 検定の結果は次の 4 種類の示し方がある

結果の示し方	特徴
数値だけ	2 群の差が分かりにくい
箱ひげ図・バイオリン図	2 群の差が可視化されて見やすく、統計的有意性が明確
棒グラフ	2 群の差が可視化されて見やすく（平均値も表記）、統計的有意性が明確
「差分」を表示	2 群の差が可視化されて見やすく（平均値も表記）、95%有意水準が明確

参考：近年の実証研究では「棒グラフ」もしくは「差分」を表示する方法が使われることが多い

1.3 \(t\) 検定の種類

• 比較するデータの数によって、「t検定」と「分散分析」を使い分ける
  
• 比較するデータの数 = 2 →　「t 検定」
  
• 比較するデータの数 > 2 →　「分散分析」
  
• 対応のあるデータ (paired data) → t 検定： pairedデータによる平均値の検定
  
• 対応のないデータ (unpaired data) → t 検定： unpairedデータによる平均値の検定

• ここでは分散分析は取り扱わない

「対応のあるデータ」と「対応のないデータ」の違い

対応のあるデータ (paired data) →　比較する 2 つの変数の対象が同じ

事例：

1. モスとマックのハンバーガー、どちらが美味しいか調べたい
   →　10人がモスとマックのハンバーガーをどちらも食べてつけた点数を比較する

2. 補講を受けた効果を調べたい
   →　20人の学生の「補講を受ける前」と「補講を受けた後」のテストの点数を比較する

3. ヤクルト1000が睡眠を促進するかどうか調べたい
   →　10人と対象として、ヤクルト1000を飲む前と飲んだ後の睡眠時間を比較する

対応のないデータ (unpaired data) →　比較する 2 つの変数の対象が異なる

事例：

1. モスとマックのハンバーガー、どちらが美味しいか調べたい
   →　10人がモスのハンバーガーだけ、別の10人がマックのハンバーガーだけを食べてつけた点数を比較する

2. 補講を受けた効果を調べたい
   →　「補講を受けた10人の学生」と「補講を受けなかった10人の学生」のテストの点数を比較する

3. ヤクルト1000が睡眠を促進するかどうか調べたい
   →　ヤクルト1000を飲んでいない10人と、ヤクルト1000を飲んだ10人の睡眠時間を比較する

2. Paired データの \(t\) 検定

2.1 データの読み取り

• 対応のあるサンプルの \(t\) 検定 (Paired-samples t-test)
  
• 比較するサンプルに「対応がある」データ
  
• 10 人がモスバーガーとマクドナルドハンバーガーを食べて味を評価する（満点は100点）
  
• 10 人はモスバーガーとマクドナルドハンバーガーの両方を食べる (=> Paired data)
  
• サンプルの平均点：モスバーガー 80.5 点、マクドナルド 79.5 点 => モスの方がおいしい
  
• 「母集団でもそうなのか？」→ \(t\) 検定を実行して確かめる必要がある

データの読み取りと準備 (mos_mc_paired.csv)

Paired（対応あり）データ

• データ mos_mc_paired.csv を読み取る

df_mos_mc <- read_csv("data/mos_mc_paired.csv")

• データフレームを表示

DT::datatable(df_mos_mc)

2.2 データ形式の変換

• このようなデータ形式は「ワイド形式 (wide format）」と呼ばれ、R 上では分析しにくい
  → tidyr::pivot_longer() 関数を使って「ロング形式 (long format）」に変換する

df_long <- df_mos_mc %>% 
    tidyr::pivot_longer(mos:mc,
                 names_to = "burger", # バーガー店名を入力  
                 values_to = "score") # バーガーの評価点を入力  

• ロング形式に変更したデータフレームを確認

DT::datatable(df_long)

2.3 得票率の平均値の可視化

2.3.1 箱ひげ図

文字化けしないためのコマンド（マックユーザのみ）

theme_set(theme_bw(base_size = 14,base_family = "HiraKakuPro-W3"))

• 箱ひげ図を描く時にはロング形式のデータを使う
  
• 簡単な箱ひげ図を描いてみる

df_long %>% 
  ggplot() +
  geom_boxplot(aes(x = burger, 
                   y = score)) 

scale_x_discrete(labels = c())関数を使って変数名を付けることができるが、順番に注意

df_long %>% 
  mutate(burger = fct_inorder(burger)) %>% 
  ggplot(aes(x = burger, y = score)) +
  geom_boxplot() +
  scale_x_discrete(labels = c( "モスバーガー", "マクドナルド")) +
  labs(x = "店名", y = "評価")

2.3.2 バイオリン図

• 箱ひげ図によく似ている「バイオリン図」の方が t 検定ではよく使われる
  
• バイオリン図では、確率密度の分布状況が分かるという利点がある

df_long %>% 
  mutate(burger = fct_inorder(burger)) %>% 
  ggplot(aes(x = burger, y = score)) +
  geom_violin() + 
  scale_x_discrete(labels = c( "モスバーガー", "マクドナルド")) +
  labs(x = "店名", y = "評価")

2.3.3 箱ひげ図＋バイオリン図

• 箱ひげ図とバイオリン図を同時に示すことも可能
  
• 大きさを指定したり、平均値を付けたり、次元ごと（バーガーごと）に色づけすることも可能

df_long %>% 
  mutate(burger = fct_inorder(burger)) %>% 
  ggplot(aes(x = burger, y = score, color = burger)) +
  geom_violin() + 
  geom_boxplot(width = .1) + # 箱ひげ図の幅を 0.1 と指定
  stat_summary(fun.y = mean, geom = "point") + # 平均値を点で示す  
  scale_x_discrete(labels = c( "モスバーガー", "マクドナルド")) +
  labs(x = "店名", y = "評価")

• モスバーガーとマックの評価点の記述統計量を表示する

summary(df_mos_mc)

      mos              mc       
 Min.   :70.00   Min.   :70.00  
 1st Qu.:76.25   1st Qu.:75.00  
 Median :80.00   Median :80.00  
 Mean   :80.50   Mean   :79.50  
 3rd Qu.:85.00   3rd Qu.:83.75  
 Max.   :90.00   Max.   :90.00  

• この標本ではマクドナルドの評価 (79.5) よりモスバーガーへの評価 (80.5) の方が 1 点高い
  
• 我々が知りたいのは母集団でもモスバーガーの方が評価が高いのか否か、ということ
  
• ひょっとしたらこの標本は、偶然、モスバーガーの方が評価が高かっただけ、という可能性もある
  
• これが偶然に得られた結果なのか、それともそうでは無いのかということを \(t\) 検定で知ることができる
  
• ttest を実行してみる(paired ttest) 　

2.4 \(t\) 検定の手順（手計算）

帰無仮説: \(H_0\):「マクドナルドとモスバーガーの評価に違いはない（＝違いは 0）」

対抗仮説: \(H_1\):「マクドナルドとモスバーガーの評価に違いはある」

\(t\)値を計算する

帰無仮説を棄却する臨界値を特定する

\(t\)値が棄却域内かどうかを確かめる

結論

Pairedデータの \(t\)値の計算式

\[T = \frac{\bar{d} - d_0}{u_x / \sqrt{n}}\]

ただし、

\[\bar{d} = \frac{\sum (x_i - y_i)}{n}\]

• \(n\) : 標本サイズ（ここでは 10)
  

• \(u_x\) : 不偏標準偏差（= 標本標準偏差）
  

• \(\bar{d}\) : マクドナルドとモスバーガーへの評価差平均
  

• \(d_0\) : 母集団で推定したい値（ここでは 0)
  

• \(x_i\) : マクドナルドへの評価
  

• \(y_i\) : モスバーガーへの評価

• この式を使って、R で \(t\) 値を計算すると次のようになる
  

• df_mos_mc というデータフレーム内の変数 mos を x、mc を y と名付ける

x <- df_mos_mc$mos
y <- df_mos_mc$mc

• マクドナルドへの評価とモスバーガーへの評価の差をとり、d と名付ける

d <- x - y

• 上記の公式に従って \(t\)値を求める

t <- (mean(d) - 0) / (sd(d) / sqrt(10))

• 求めた \(t\) 値を表示する

t

[1] 0.557086

• ここで得られた 0.557086 という \(t\) 値を使って「マクドナルドとモスバーガーのどちらの方が評価が高いのか？」ということを標本サイズ 10 のサンプルを使って、母集団では「どちらの評価も同じ」であるという仮説が妥当かどうかを推定する
  
• t 分布表は次のとおり

• \(t\) 分布表内の数値は、検定を行う際の有意水準 (significance level) に対応した棄却限界値を示している
  
• 有意水準のことを \(\alpha\) と表記する
  
• 標本から得られた \(t\) 値が棄却限界値より大きい場合、帰無仮説を棄却する
  
• 横軸は片側検定における有意水準を示している
  
• 例えば、二列目の「確率95%」は両側検定における有意水準が 5% （\(α = 0.05\)）であることを示す
  
• \(t\) 検定では両側検定で 5% （\(α = 0.05\)）の有意水準を使うのが一般的
  
• 従って、 \(t\) 分布表二列目の「確率95%」を使う
  
• 分布表の一列目は自由度 (df: degree of freedom) を示している
  
• ここでは「標本サイズ 10 から 1 を引いた 9」 が自由度
  
• 縦軸の 9 と横軸の「0.025」交差する数字「2.262」が 5% の両側検定における棄却限界値
  
• これを図示すると次のようになる
  
• 自由度 9 の \(t\) 分布における 5% 有意水準（\(α = 0.05\)）の 2 つの棄却限界値 (critical value) -2.26 と 2.26 を示している

解釈 ・サンプルから得られた \(t\)値 (0.557086) が -2.26 と 2.26 の採択区間内にある
→帰無仮説 \(H_0\)「マクドナルドとモスバーガーへの評価は同じ」は棄却できない
・「得られた 10 の標本サイズでマクドナルドへの評価が 79.5、モスバーガーへの評価が 80.5 と、1 点だけモスバーガーが高く評価されたが、これは偶然の可能性が高い」

2.5 \(t\) 検定の手順（Rで計算）

2.5.1 ロング形式データの場合

df_long

# A tibble: 20 × 2
   burger score
   <chr>  <dbl>
 1 mos       80
 2 mc        75
 3 mos       75
 4 mc        70
 5 mos       80
 6 mc        80
 7 mos       90
 8 mc        85
 9 mos       85
10 mc        90
11 mos       80
12 mc        75
13 mos       75
14 mc        85
15 mos       85
16 mc        80
17 mos       85
18 mc        80
19 mos       70
20 mc        75

• 簡単な箱ひげ図を描いてみる

df_long %>% 
  ggplot() +
  geom_boxplot(aes(x = burger, 
                   y = score)) 

t.test(df_long$score[df_long$burger == "mos"],
       df_long$score[df_long$burger == "mc"],
       paired = TRUE)          # unpaired が default なので paired と指定

    Paired t-test

data:  df_long$score[df_long$burger == "mos"] and df_long$score[df_long$burger == "mc"]
t = 0.55709, df = 9, p-value = 0.5911
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -3.060696  5.060696
sample estimates:
mean difference 
              1 

解釈 ・帰無仮説:
「マクドナルドもモスバーガーも評価に差はない」
・4 行目には t = 0.55709 が示されている
・両側検定の p-value = 0.5911 が 0.05 よりも大きい
　→5% の有意水準（α = 0.05）で帰無仮説を棄却できない
　→「マクドナルドもモスバーガーも評価に差がない」
　 得られた標本ではマクドナルドへの評価が 79.5、モスバーガーへの評価が 80.5 と、1 点だけモスバーガーが高く評価されたが、これは偶然の可能性が高い

• 統計的有意性 (\(p\) 値) を加えた箱ひげ図を描くこともできる
  
• ggsignif()関数では \(t\) 検定は unpaired data がデフォルト
  
• ここでは paired の t 検定　→下から3行目：　test.args = list(paired = TRUE)と指定
  
• unpaired の t 検定なら → test.args = list(paired = TRUE)は不要
  
• もしくは test.args = list(paired = FALSE) と指定

df_long %>% 
  mutate(burger = fct_inorder(burger)) %>% 
  ggplot(aes(x = burger, y = score, color = burger)) +
  geom_violin() +
  geom_boxplot(width = .1) + # 箱ひげ図の幅を 0.1 と指定
  stat_summary(fun.y = mean, geom = "point") + # 平均値を点で示す 
  scale_x_discrete(labels = c( "モスバーガー", "マクドナルド")) +
  labs(x = "店名", y = "評価") +
    ggsignif::geom_signif(comparisons = combn(sort(unique(df_long$burger)), 2, FUN = list),
                          test = "t.test",　　　　　　　　# t 検定だと指定
                          test.args = list(paired = TRUE), # unpaired なら、paired = FALSE
                          na.rm = T,
                          step_increase = 0.1)

2.5.2 ワイド形式データの場合

df_mos_mc

# A tibble: 10 × 2
     mos    mc
   <dbl> <dbl>
 1    80    75
 2    75    70
 3    80    80
 4    90    85
 5    85    90
 6    80    75
 7    75    85
 8    85    80
 9    85    80
10    70    75

t.test(df_mos_mc$mos,
       df_mos_mc$mc, 
       paired = TRUE) # unpaired が default なので paired = TRUE と指定

    Paired t-test

data:  df_mos_mc$mos and df_mos_mc$mc
t = 0.55709, df = 9, p-value = 0.5911
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -3.060696  5.060696
sample estimates:
mean difference 
              1 

2.5.3 平均値の差の検定 = \(t\) 検定

One Sample t-test

• ここで行った ttest は次のようにしても同じ結果を求めることができる
  

• One Sample t-testとしても ttest できる
  

• モスバーガーとマクドナルドの平均値 diff を計算する
  

• diff の平均値（サンプル平均）を計算する
  

• この統計量を使ってパラメタを推定する
  

• ワイド形式データ (df_mos_mc) を使って計算してみる

• 帰無仮説：diff の平均値は 0 = モスとマックの味に差はない

diff <- df_mos_mc$mos - df_mos_mc$mc
diff

 [1]   5   5   0   5  -5   5 -10   5   5  -5

• diff の平均値は

mean(diff)

[1] 1

t.test(diff)

    One Sample t-test

data:  diff
t = 0.55709, df = 9, p-value = 0.5911
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -3.060696  5.060696
sample estimates:
mean of x 
        1 

3. Unpaired データの \(t\) 検定

• 前のセクションでは、Paired データに関して \(t\) 検定の方法を演習した
  
• ここでは Unpaired データを使って演習する
  
• 対応のないサンプルの \(t\) 検定 (Unpaired-samples t-test)
  
• 比較するサンプルに「対応がない」データ
  
• 「マクドナルドとモスバーガーを 0 点から 100 点で評価して下さい」という世論調査結果を人工的に生成
  
• 標本サイズ：20 人
  
• 10 人がマクドナルドだけを評価し、残りの 10 人がモスバーガーだけを評価して点数をつける
  ＝ Unpaired data （対応のないサンプル）

データの種類	解説
Paired（対応あり）	バーガーを評価した人は 10 人、両方のバーガーを食べた
Unpaired（対応なし）	バーガーを評価した人は 20 人、片方のバーガーだけ食べた

Unpaired（対応なし）データ

• 標本をとり、それぞれの標本平均に違いがあったとする
  
• しかし、その違いは「標本における」違いに過ぎず、その結果を一般化することはできない

3.1 データの読み取り

• データ mos_mc_paired.csv を読み取る

df_mos_mc <- read_csv("data/mos_mc_paired.csv")

• データフレームを表示

DT::datatable(df_mos_mc)

3.2 データ形式の変換

• このようなデータ形式は「ワイド形式 (wide format）」と呼ばれ、R 上では分析しにくい
  → tidyr::pivot_longer() 関数を使って「ロング形式 (long format）」に変換する

df_long <- df_mos_mc %>% 
    tidyr::pivot_longer(mos:mc,
                 names_to = "burger", # バーガー店名を入力  
                 values_to = "score") # バーガーの評価点を入力  

• ロング形式に変更したデータフレームを確認

DT::datatable(df_long)

3.3 箱ひげ図を描く

文字化けしないためのコマンド（マックユーザのみ）

theme_set(theme_bw(base_size = 14,base_family = "HiraKakuPro-W3"))

• 箱ひげ図を描く時にはロング形式のデータを使う
  
• 簡単な箱ひげ図を描いてみる

df_long %>% 
  ggplot() +
  geom_boxplot(aes(x = burger, 
                   y = score)) 

df_long %>% 
  mutate(burger = fct_inorder(burger)) %>% 
  ggplot(aes(x = burger, y = score, color = burger)) +
  geom_violin() +
  geom_boxplot(width = .1) + # 箱ひげ図の幅を 0.1 と指定
  stat_summary(fun.y = mean, geom = "point") + # 平均値を点で示す 
  scale_x_discrete(labels = c( "モスバーガー", "マクドナルド")) +
  labs(x = "店名", y = "評価")

• これはセッション 2.3 で描いた箱ひげ図と同じもの
  
• 箱ひげ図をみると、この標本ではマクドナルドの評価 (79.5) よりモスバーガーへの評価 (80.5) の方が 1 点高い
  
• 我々が知りたいのは母集団でもモスバーガーの方が評価が高いのか否か、ということ
  
• ひょっとしたらこの標本は、偶然、モスバーガーの方が評価が高かっただけ、という可能性もある
  
• これが偶然に得られた結果なのか、それともそうでは無いのかということを \(t\) 検定で知ることができる
  
• ttest を実行してみる(unpaired ttest) 　

3.4 \(t\) 検定の手順（Rで計算）

3.4.1 ロング形式データの場合

df_long

# A tibble: 20 × 2
   burger score
   <chr>  <dbl>
 1 mos       80
 2 mc        75
 3 mos       75
 4 mc        70
 5 mos       80
 6 mc        80
 7 mos       90
 8 mc        85
 9 mos       85
10 mc        90
11 mos       80
12 mc        75
13 mos       75
14 mc        85
15 mos       85
16 mc        80
17 mos       85
18 mc        80
19 mos       70
20 mc        75

t.test(df_long$score[df_long$burger == "mos"],
       df_long$score[df_long$burger == "mc"]) # unpaired が default 

    Welch Two Sample t-test

data:  df_long$score[df_long$burger == "mos"] and df_long$score[df_long$burger == "mc"]
t = 0.37354, df = 18, p-value = 0.7131
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -4.624301  6.624301
sample estimates:
mean of x mean of y 
     80.5      79.5 

解釈 ・帰無仮説:
「マクドナルドもモスバーガーも評価に差はない」
・4 行目には t = 0.37354 が示されている
・両側検定の p-value = 0.7131 が 0.05 よりも大きい
　→5% の有意水準（α = 0.05）で帰無仮説を棄却できない
　→「マクドナルドもモスバーガーも評価に差がない」
　 　得られた標本ではマクドナルドへの評価が 79.5、モスバーガーへの評価が 80.5 と、1 点だけモスバーガーが高く評価されたが、これは偶然の可能性が高い

• 簡単な箱ひげ図を描いてみる

df_long %>% 
  ggplot() +
  geom_boxplot(aes(x = burger, 
                   y = score)) 

• 観測値と統計的有意性 (\(p\) 値) を加えた箱ひげ図を描くこともできる
  
• ggsignif()関数における \(t\) 検定は unpaired data がデフォルト
  
• ここでは unpaired の t 検定 → test.args = list(paired = TRUE)は不要
• もしくは test.args = list(paired = FALSE) と指定

df_long %>% 
  mutate(burger = fct_inorder(burger)) %>% 
  ggplot(aes(x = burger, y = score, fill = burger)) +
  geom_violin() + 
  geom_boxplot(width = .1) + # 箱ひげ図の幅を 0.1 と指定
  stat_summary(fun.y = mean, geom = "point") + # 平均値を点で示す 
    ggbeeswarm::geom_beeswarm() +
  scale_x_discrete(labels = c( "モスバーガー", "マクドナルド")) +
  labs(x = "店名", y = "評価") +
    ggsignif::geom_signif(comparisons = combn(sort(unique(df_long$burger)), 2, FUN = list),
                          test = "t.test",
                          na.rm = T,
                          step_increase = 0.1)

3.4.2 ワイド形式データの場合

df_mos_mc

# A tibble: 10 × 2
     mos    mc
   <dbl> <dbl>
 1    80    75
 2    75    70
 3    80    80
 4    90    85
 5    85    90
 6    80    75
 7    75    85
 8    85    80
 9    85    80
10    70    75

t.test(df_mos_mc$mos,
       df_mos_mc$mc) # unpaired が default 

    Welch Two Sample t-test

data:  df_mos_mc$mos and df_mos_mc$mc
t = 0.37354, df = 18, p-value = 0.7131
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -4.624301  6.624301
sample estimates:
mean of x mean of y 
     80.5      79.5 

4.総選挙データを使った \(t\) 検定

得票率の比較（自民と立憲）

4.1 データの準備 (hr96-21.csv)

4.1.1 データのダウンロード方法

予めダウンロード先を指定する方法

• RProject フォルダ内に data という名称のフォルダを作成する
  
• 下のコマンドを実行すると、csv ファイルがパソコンにダウンロードされ、data 内に自動的に保存される

download.file(url = "https://asanoucla.github.io/hr96-21.csv",
              destfile = "data/hr96-21.csv")

注意：一度ダウンロードを完了すれば、このコマンドを実行する必要はありません

ダウンロード先を指定しない方法

• hr96-21.csv をクリックしてデータをパソコンにダウンロード　　

• RProject フォルダ内に data という名称のフォルダを作成する
  

• ダウンロードした hr96-21.csv を手動でRProject フォルダ内にある data フォルダに入れる

4.1.2 選挙データの読み取り方法

• 次のいずれかの方法で hr96-21.csv を読み取る

読み取り方法 1

• na = "."というコマンドは「欠損値をドットで置き換える」という意味
  
• 欠損値を空欄のまま残すと、本来「数値 (numeric)」型のデータが「」文字型 (character)」として認識されるなど、エラーの原因になるため、読み取る時点で事前に対処する

hr <- read_csv("data/hr96-21.csv",
               na = ".")  

読み取り方法 2

• 読み取った値の日本語が文字化けする場合
  
• locale()関数を使って日本語エンコーディング (cp932) を指定する

hr <- read_csv("data/hr96-21.csv",
               na = ".",
               locale = locale(encoding = "cp932"))

読み取り方法 3

hr <- read.csv("data/hr96-21.csv",
               na = ".")  

4.1.3 読み取った選挙データを確認

• hr96_21.csv は1996年に衆院選挙に小選挙区が導入されて以来実施された 9 回の衆議院選挙（1996, 2000, 2003, 2005, 2009, 2012, 2014, 2017, 2021）の結果のデータ
  
• hr に含まれる変数名を表示させる

names(hr)

 [1] "year"          "pref"          "ku"            "kun"          
 [5] "wl"            "rank"          "nocand"        "seito"        
 [9] "j_name"        "gender"        "name"          "previous"     
[13] "age"           "exp"           "status"        "vote"         
[17] "voteshare"     "eligible"      "turnout"       "seshu_dummy"  
[21] "jiban_seshu"   "nojiban_seshu"

• df1 には 22 個の変数が入っている

変数名	詳細
year	選挙年 (1996-2017)
pref	都道府県名
ku	小選挙区名
kun	小選挙区
rank	当選順位
wl	選挙の当落: 1 = 小選挙区当選、2 = 復活当選、0 = 落選
nocand	立候補者数
seito	候補者の所属政党
j_name	候補者の氏名（日本語）
name	候補者の氏名（ローマ字）
previous	これまでの当選回数（当該総選挙結果は含まない）
gender	立候補者の性別: “male”, “female”
age	立候補者の年齢
exp	立候補者が使った選挙費用（総務省届け出）
status	候補者のステータス: 0 = 非現職、1 現職、2 = 元職
vote	得票数
voteshare	得票率 (%)
eligible	小選挙区の有権者数
turnout	小選挙区の投票率 (%)
seshu_dummy	世襲候補者ダミー: 1 = 世襲、0 = 非世襲（地盤世襲 or 非世襲）
jiban_seshu	地盤の受け継ぎ元の政治家の氏名と関係
nojiban_seshu	世襲元の政治家の氏名と関係

• データの型をチェック

str(hr)

spc_tbl_ [9,660 × 22] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ year         : num [1:9660] 1996 1996 1996 1996 1996 ...
 $ pref         : chr [1:9660] "愛知" "愛知" "愛知" "愛知" ...
 $ ku           : chr [1:9660] "aichi" "aichi" "aichi" "aichi" ...
 $ kun          : num [1:9660] 1 1 1 1 1 1 1 2 2 2 ...
 $ wl           : num [1:9660] 1 0 0 0 0 0 0 1 0 2 ...
 $ rank         : num [1:9660] 1 2 3 4 5 6 7 1 2 3 ...
 $ nocand       : num [1:9660] 7 7 7 7 7 7 7 8 8 8 ...
 $ seito        : chr [1:9660] "新進" "自民" "民主" "共産" ...
 $ j_name       : chr [1:9660] "河村たかし" "今枝敬雄" "佐藤泰介" "岩中美保子" ...
 $ gender       : chr [1:9660] "male" "male" "male" "female" ...
 $ name         : chr [1:9660] "KAWAMURA, TAKASHI" "IMAEDA, NORIO" "SATO, TAISUKE" "IWANAKA, MIHOKO" ...
 $ previous     : num [1:9660] 2 2 2 0 0 0 0 2 0 0 ...
 $ age          : num [1:9660] 47 72 53 43 51 51 45 51 71 30 ...
 $ exp          : num [1:9660] 9828097 9311555 9231284 2177203 NA ...
 $ status       : num [1:9660] 1 2 1 0 0 0 0 1 2 0 ...
 $ vote         : num [1:9660] 66876 42969 33503 22209 616 ...
 $ voteshare    : num [1:9660] 40 25.7 20.1 13.3 0.4 0.3 0.2 32.9 26.4 25.7 ...
 $ eligible     : num [1:9660] 346774 346774 346774 346774 346774 ...
 $ turnout      : num [1:9660] 49.2 49.2 49.2 49.2 49.2 49.2 49.2 51.8 51.8 51.8 ...
 $ seshu_dummy  : num [1:9660] 0 0 0 0 0 0 0 0 1 0 ...
 $ jiban_seshu  : chr [1:9660] NA NA NA NA ...
 $ nojiban_seshu: chr [1:9660] NA NA NA NA ...
 - attr(*, "spec")=
  .. cols(
  ..   year = col_double(),
  ..   pref = col_character(),
  ..   ku = col_character(),
  ..   kun = col_double(),
  ..   wl = col_double(),
  ..   rank = col_double(),
  ..   nocand = col_double(),
  ..   seito = col_character(),
  ..   j_name = col_character(),
  ..   gender = col_character(),
  ..   name = col_character(),
  ..   previous = col_double(),
  ..   age = col_double(),
  ..   exp = col_double(),
  ..   status = col_double(),
  ..   vote = col_double(),
  ..   voteshare = col_double(),
  ..   eligible = col_double(),
  ..   turnout = col_double(),
  ..   seshu_dummy = col_double(),
  ..   jiban_seshu = col_character(),
  ..   nojiban_seshu = col_character()
  .. )
 - attr(*, "problems")=<externalptr> 

• 数値は numeric 文字は character として認識されていることがわかる

4.2 得票率の \(t\) 検定（2 つの政党間）

• 2017年総選挙での自民党と立憲民主党が得た得票率の箱ひげ図を描いてみよう

• 次の 2 つの条件を指定して hr からデータを抜き出す
  

1. 2017年の総選挙だけ
   
2. 自民党と立憲民主党の候補者だけ
   

• 抜き出したデータフレームに ldp_cdp17 と何前をつける

ldp_cdp17 <- hr %>% 
  dplyr::filter(year == 2017) %>%                          #2017年に実施された総選挙だけ  
  dplyr::filter(seito == "自民" | seito == "立憲") %>% #自民党と立憲民主党だけ
  dplyr::select(year, voteshare, seito)    #三つの変数だけ    

• 変数の中身を確認する

DT::datatable(ldp_cdp17)

• 箱ひげ図を描いてみる

ldp_cdp17 %>% 
  ggplot(aes(x = seito, y = voteshare, color = seito)) +
  geom_violin() +
  geom_boxplot(width = .1) + # 箱ひげ図の幅を 0.1 と指定
  stat_summary(fun.y = mean, geom = "point") + # 平均値を点で示す 
  labs(x = "政党名", y = "得票率 (%)")

• 二つの政党が得た得票率を箱ひげ図で表示すれば、二つの政党が得た得票率の分布がわかる
  
• しかし、両者の間に「差がある」かどうかを確かめるためには \(t\) 検定を実行する必要がある
  
• 自民党と立憲民主党の候補者が得た得票率は「対応のない \(t\) 検定」に該当
  →　unpaired ttest を実行
  →　R ではunpaired ttest がデフォルト
  →　paired = FALSE という指定は不要

t.test(ldp_cdp17$voteshare[ldp_cdp17$seito == "自民"], 
       ldp_cdp17$voteshare[ldp_cdp17$seito == "立憲"])   

    Welch Two Sample t-test

data:  ldp_cdp17$voteshare[ldp_cdp17$seito == "自民"] and ldp_cdp17$voteshare[ldp_cdp17$seito == "立憲"]
t = 9.4581, df = 86.8, p-value = 5.287e-15
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 11.77808 18.04572
sample estimates:
mean of x mean of y 
 50.39634  35.48444 

• 箱ひげ図をみると、この標本では立憲民主党の得票率平均 (35.5%) より自民党の得票率平均 (50.4%) の方が 14.9 %高い
  
• 我々が知りたいのは母集団でも自民党の方が評価が高いのか否か、ということ
  
• ひょっとしたらこの標本は、偶然、自民党の方が得票率が高かっただけ、という可能性もある
  
• これが偶然に得られた結果なのか、それともそうでは無いのかということを \(t\) 検定で知ることができる

解釈 ・帰無仮説:
「立憲民主党も自民党も得票率に差はない」
・4 行目には t = 9.4518 が示されている
・両側検定の p-value = 5.54e-15 が 0.01 よりも遙かに小さい
　→1% の有意水準（α = 0.01）で帰無仮説を棄却
　→「立憲民主党も自民党も得票率に差がない」わけではない
・5 行目に alternative hypothesis: true difference in means is not equal to 0 と対抗仮説が明示されている
・サンプルでの得票率が高い自民党の方が、母集団においても得票率が高いという対抗仮説を採択する（自民党　50.4%、立憲民主党　35.5%）

立憲民主党より自民党の方が、母集団においても得票率が高い

箱ひげ図上に \(t\) 検定の結果 (p値) を表示させたい場合

library(ggsignif)

• 最もシンプルな方法
  
• ここでは unpaired の t 検定 → test.args = list(paired = TRUE)は不要
• もしくは test.args = list(paired = FALSE) と指定
  （ここでは設定を省きtest.args = list(paired = TRUE)を書き入れていない）

ldp_cdp17 %>% 
    ggplot(aes(seito, voteshare)) + 
  geom_boxplot() +
  geom_signif(comparisons = list(c("立憲", "自民")), # 分類の指定  
              test = "t.test")     　　　　　　　　　# t 検定だと指定

- バイオリンプロットを加えて変数名もカスタマイズしたいとき

ldp_cdp17 %>% 
    ggplot(aes(seito, voteshare, fill = seito)) + 
    geom_violin() + 
  　geom_boxplot(width = .1) +     # 箱ひげ図の幅を 0.1 と指定
  　stat_summary(fun.y = mean, 
  　             geom = "point") + # 平均値を点で示す 
    geom_signif(comparisons = list(c("立憲", "自民")),  # 分類の指定  
                test = "t.test",   # t 検定だと指定
                na.rm = T) +       # 欠損値を省く設定
  labs(x = "政党名", y = "得票率 (%)") +
  theme(legend.position = "none")

4.3 得票率の \(t\) 検定（3 つの政党間）

• 前のセクションでは自民党と立憲民主党の得票率の平均値を \(t\) 検定した
  
• ここでは 3 つの政党間の得票率の平均値をペアごとに \(t\) 検定する方法を紹介する
  
• 次の 2 つの条件を指定して hr からデータを抜き出す
• 2017年の総選挙だけ
• 自民党と立憲民主党と希望の党の候補者だけ
• 抜き出したデータフレームに ldp_cdp_poh17 と何前をつける

ldp_cdp_poh17 <- hr %>% 
  filter(year == 2017) %>% 
  filter(seito == "自民" | seito == "立憲" |  seito == "希望") %>% 
  select(voteshare, seito)

• データフレーム ldp_cdp_poh17 が含む変数 seito の中身を表示させる

unique(ldp_cdp_poh17$seito)

[1] "自民" "立憲" "希望"

• 自民党と立憲民主党の間で得票率の差を確認する

library(ggsignif)

• ここでは unpaired の t 検定 → test.args = list(paired = TRUE)は不要
• もしくは test.args = list(paired = FALSE) と指定
  （ここでは設定を省きtest.args = list(paired = TRUE)を書き入れていない）

ldp_cdp_poh17 %>% 
    ggplot(aes(seito, voteshare, fill = seito)) + 
    geom_violin() + 
    geom_boxplot(width = .1) + # 箱ひげ図の幅を 0.1 と指定
    stat_summary(fun.y = mean, geom = "point") + # 平均値を点で示す 
    geom_signif(comparisons = list(c("立憲", "自民")),
                test = "t.test",
                na.rm = T,
                step_increase = 0.1)　+
  theme(legend.position = "none")

• 自民党と立憲民主党、自民党と希望の党の間で得票率の差を確認する

ldp_cdp_poh17 %>% 
    ggplot(aes(seito, voteshare, fill = seito)) + 
    geom_violin() + 
    geom_boxplot(width = .1) + # 箱ひげ図の幅を 0.1 と指定
    stat_summary(fun.y = mean, geom = "point") + # 平均値を点で示す 
    geom_signif(comparisons = list(c("立憲", "自民"),
                                   c("希望", "自民")),
                test = "t.test",
                na.rm = T,
                step_increase = 0.1)　+
  theme(legend.position = "none")

• 自民党と立憲民主党、自民党と希望の党、希望の党と立憲民主党の間で得票率の差を確認する

ldp_cdp_poh17 %>% 
    ggplot(aes(seito, voteshare, fill = seito)) + 
    geom_violin() + 
    geom_boxplot(width = .1) + # 箱ひげ図の幅を 0.1 と指定
    stat_summary(fun.y = mean, geom = "point") + # 平均値を点で示す 
    geom_signif(comparisons = list(c("立憲", "自民"),
                                   c("希望", "自民"),
                                   c("希望", "立憲")),
                test = "t.test",
                na.rm = T,
                step_increase = 0.1)　+
  theme(legend.position = "none")

解釈 ・三つの政党間それぞれにおいて得票率の平均に差はある
・民主党が二つ（希望の党と立憲民主党）に分党したが、立憲民主党の方が得票率が高い（1 %で統計的に有意）

二群同士の \(t\) 検定で結果を確かめる（政党間得票率の比較）

• 自民党、立憲民主党、希望の党、共産党の間で得票率の差は母集団が異なる unpaired data なので「対応のない \(t\) 検定」を実行する

自民党 vs 立憲民主党

• p value = 5.5e-15

# unpaired が default（つまり unpaired がデフォルト） 
t.test(ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "自民"],
       ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "立憲"]) 

    Welch Two Sample t-test

data:  ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "自民"] and ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "立憲"]
t = 9.4581, df = 86.8, p-value = 5.287e-15
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 11.77808 18.04572
sample estimates:
mean of x mean of y 
 50.39634  35.48444 

自民党 vs 希望の党

• p value = 2.22e-16

# unpaired が default  
t.test(ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "自民"],
       ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "希望"]) 

    Welch Two Sample t-test

data:  ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "自民"] and ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "希望"]
t = 19.823, df = 397.34, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 18.50649 22.58137
sample estimates:
mean of x mean of y 
 50.39634  29.85241 

立憲民主党 vs 希望の党

• p value = 0.00087

# unpaired が default  
t.test(ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "立憲"],
       ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "希望"]) 

    Welch Two Sample t-test

data:  ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "立憲"] and ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "希望"]
t = 3.3819, df = 105.41, p-value = 0.001011
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 2.330125 8.933940
sample estimates:
mean of x mean of y 
 35.48444  29.85241 

5. 有意水準に関するシミュレーション

5.1 有意水準 5% の意味を考える

• 「有意水準と有意確率」の関係に関して「有意水準 5 %（α = 0.05）」が意味するのは、帰無仮説が正しい時、誤って帰無仮説をたまたま棄却してしまう危険性が 5%（100回のうち 5回）存在すること
  
• ここでは R を使って母数を指定した架空の母集団から標本を採取するシミュレーションを 10,000 回行い、上記のようなことが実際に起こるのかどうか確かめてみる
  
• まず、同一の母集団から標本サイズ 50、標本数 2 のサンプルを採取して二つの標本平均を比較する
  
• 二つの標本は同一の母集団から採取されているのだから、標本平均は同じはず
  
• しかし、現実に標本を採取してみると、標本平均は異なる
  
• 有意水準 5% (α = 0.05) に設定したシミュレーションを行えば、理論的には 100回に 5回は「母集団の平均値は同じ」という帰無仮説を誤って棄却するはずである

sig_sim <- function(alpha = .05, n = 50, trials = 100) {
    ## Arguments:
    ##     alpha = significance level（有意水準）、defaultで 5% (α = 0.05)に設定
    ##     n = 標本サイズ、defaultで 50 に設定
    ##     trials = シミュレーションの試行数, defaultで 100 に設定
    ## Return:
    ##     帰無仮説を棄却する割合 (rejection rate) を表示させる
    ## 　　t 分布もヒストグラムとして表示させる
    
    ## vector to save the result
    res <- rep(NA, trials)
    T_vec <-rep(NA, trials)
    
    ## critical value
    cv <- abs(qt(alpha / 2, df = n -1))
    
    for (i in 1:trials) {
        x <- rnorm(50, mean = 50, sd = 5) # 標本 x は母平均 50、母標準偏差 5 から 50 個無作為抽出
        y <- rnorm(50, mean = 50, sd = 5) # 標本 y は母平均 50、母標準偏差 5 から 50 個無作為抽出
        d <- x - y                        # x と y の差を d とする
        T <- (mean(d) - 0) / (sd(d) / sqrt(n)) # サンプルから得られる $t$値を計算
        T_vec[i] <- T
        res[i] <- abs(T) > cv　　　　# $t$値の絶対値が cv(棄却限界値)を超える場合を選択
    }
    hist(T_vec, freq = FALSE, col = "gray", # 計算された $t$値をヒストグラムとして表示
         xlab = "test statistic",
         main = "Distribution of the test statistic")
    abline(v = c(-cv, cv), col = "red")     # cv(棄却限界値)の上限と下限それぞれに赤線
    return(mean(res))                       # 帰無仮説を棄却する割合 (rejection rate) を表示
}
sig_sim(trials = 10000)  # シミュレーションの試行回数を指定

[1] 0.0482

• 上の数値は、帰無仮説を棄却する割合 (rejection rate)
  

• ここではシミュレーションを行っているので、試行するたびにこの数値は異なる
  

• 上のヒストグラムは 10000 回の標本採取シミュレーションの結果得られた t 分布を示す
  

• -2 と 2 に縦に引かれている赤線は 5% 有意水準 (α = 0.05) の棄却限界値を示している
  

• 例えば、ここで得られた数値が 0.0486 だとする（試行するたびに異なることに注意）
  

• -2 より左側と 2 より右側が帰無仮説の棄却域であり、10000回の標本試行において、帰無仮説が棄却された割合が 0.0486（つまり約 5%）となる
  

• ここでは母数を予め指定した架空の母集団を作り、同じ母集団から二つの標本を採取した
  

• 母集団を推定するために設定した帰無仮説は「母平均は同じである」である（これは「事実」である）

• つまり、同一の母集団から標本を採取した結果、予想通り、100 回のうち約 5 回は、帰無仮説「母平均は同じである」を棄却し、誤った結果を得たことになる

次に、有意水準を 10% (α = 0.10) に変更して、同様のシミュレーションを行ってみる

#有意水準を 10% (α = 0.10) 、シミュレーションの試行回数は 10000回
sig_sim(alpha = .1, trials = 10000)  

[1] 0.0947

結果 ・予想どおり、t 分布における二本の棄却限界値 を示す赤線はどちらも中央に寄り、帰無仮説を棄却する割合(rejection rate)が約 10% (0.1に近い値) まで増えていることがわかる

6. \(t\) 検定の可視化（棒グラフ）

• 近年、様々な学会誌 (Academic Journal) では \(t\) 検定の結果が様々な方法で可視化
  

• ここでは棒グラフを使った \(t\) 検定の可視化方法を紹介
  

• ここではモスバーガーとマクドナルドそれぞれの店で「チーズバーガー」「フライドポテト」「テリヤキバーガー」「シェイク」を食べてもらい点数をつけたという架空データを使う
  

• 回答者はモスバーガーもしくはマクドナルドのいずれかで4種類食べたとする

• ここで知りたいこと

モスとマックの「フライドポテト」はどちらがおいしいか

6.1 データの読み取りと準備

• menu.csv を読み込み df_menu と名前を付ける（架空データ）

df_menu <- read_csv("data/menu.csv")

DT::datatable(df_menu)

6.2 mos と mc ポテトの \(t\) 検定

• mos ポテトと mc ポテトの点数の \(t\) 検定を実行
  
• unpaired data なので対応のない \(t\) 検定 = Welch の \(t\) 検定を実行

potato_ttest <- t.test(fried_potato ~ mosmc, data = df_menu)
potato_ttest

    Welch Two Sample t-test

data:  fried_potato by mosmc
t = 4.4463, df = 15.32, p-value = 0.0004489
alternative hypothesis: true difference in means between group mc and group mos is not equal to 0
95 percent confidence interval:
 3.233228 9.166772
sample estimates:
 mean in group mc mean in group mos 
             80.9              74.7 

解釈 ・帰無仮説:
「マクドナルドのポテトとモスバーガーのポテトの味の評価に差はない」
・3 行目には t = 4.4463 が示されている
・両側検定の p-value = 0.0004489 が 0.01 よりも遙かに小さい
　→1% の有意水準（α = 0.01）で帰無仮説を棄却
　→「マクドナルドのポテトとモスバーガーのポテトの味の評価に差はない」わけではない
・4 行目に alternative hypothesis: true difference in means is not equal to 0 と対抗仮説が明示されている
・サンプルでの評価が高いマクドナルドの方が、母集団においても評価が高いという対抗仮説を採択する（マクドナルド　80.9点、モスバーガー　74.4点）
ポテトの評価はモスバーガーよりマクドナルドの方が高い（1%で統計的に有意）

6.3 棒グラフによる可視化

統計的有意性が非表示の場合

• ここで使うデータセット

library(DT)
library(tidyverse)
library(ggsignif)
library(Rcpp)

DT::datatable(df_menu)

• 棒グラフで表示する次の変数を作成する

変数名	内容
mean_mosmc	モスとマックのバーガーの味の点数
sd_mosmc	標準偏差 (standard deviation)
se_mosmc	標準誤差 (standard error
count	ケース数

df_eb <- df_menu |> 
  mutate(potato = fried_potato) |> # fried_potato を potato に変更
  select(mosmc, potato)

names(df_eb)

[1] "mosmc"  "potato"

df_eb <- df_eb |> 
  group_by(mosmc) |> 
  summarise(min = min(potato),
            max = max(potato),
            mean_potato = mean(potato),
            sd_potato = sd(potato),
            count = n(),
            se_potato = (sd_potato)/(sqrt(count)))

• 計算した値をチェック

df_eb

# A tibble: 2 × 7
  mosmc   min   max mean_potato sd_potato count se_potato
  <chr> <dbl> <dbl>       <dbl>     <dbl> <int>     <dbl>
1 mc       76    83        80.9      2.38    10     0.752
2 mos      70    81        74.7      3.71    10     1.17 

plt1 <- ggplot(df_eb, aes(x = mosmc,
             y = mean_potato, fill = mosmc)) +
  geom_bar(stat = "identity",
           position = position_dodge(width = 0.8)) +
  geom_errorbar(aes(ymin = mean_potato - se_potato,
                    ymax = mean_potato + se_potato),
                    width = 0.2) +
  scale_fill_manual(values = c("skyblue", "pink")) +
  labs(x = "バーガー店",
       y = "得点") +
  ggtitle("ポテトの点数比較：エラーバーは標準誤差") +
  theme_bw(base_size = 14,base_family = "HiraKakuPro-W3") +
  theme(legend.position = "none") +
  geom_label(aes(label = mean_potato),
             size = 7.5, position = position_stack(vjust = 0.5),
             show.legend = F, fill = "white")
  
plt1

エラーバーに「標準偏差」を表示させる場合 ・「標準誤差」ではなく「標準偏差」を表示させることも可能
・この場合、ymin と ymax の指定を次のように変更する ymin = mean_potato - sd_potato,
ymax = mean_potato + sd_potato

plt2 <- ggplot(df_eb, aes(x = mosmc,
             y = mean_potato, fill = mosmc)) +
  geom_bar(stat = "identity",
           position = position_dodge(width = 0.8)) +
  geom_errorbar(aes(ymin = mean_potato - sd_potato,
                    ymax = mean_potato + sd_potato),
                    width = 0.2) +
  scale_fill_manual(values = c("skyblue", "pink")) +
  labs(x = "バーガー店",
       y = "得点") +
  ggtitle("ポテトの点数比較：エラーバーは標準偏差") +
  theme_bw(base_size = 14,base_family = "HiraKakuPro-W3") +
  theme(legend.position = "none") +
  geom_label(aes(label = mean_potato),
             size = 7.5, position = position_stack(vjust = 0.5),
             show.legend = F, fill = "white")
  
plt2

エラーバーに「最小値・最大値」を表示させる場合 ・「標準誤差」ではなく「最小値・最大値」を表示させることも可能
・この場合、ymin と ymax の指定を次のように変更する:
ymin = min,
ymax = max

plt3 <- ggplot(df_eb, aes(x = mosmc,
             y = mean_potato, fill = mosmc)) +
  geom_bar(stat = "identity",
           position = position_dodge(width = 0.8)) +
  geom_errorbar(aes(ymin = min,
                    ymax = max),
                    width = 0.2) +
  scale_fill_manual(values = c("skyblue", "pink")) +
  labs(x = "バーガー店",
       y = "得点") +
  ggtitle("ポテトの点数比較：エラーバーは最大値・最小値") +
  theme_bw(base_size = 14,base_family = "HiraKakuPro-W3") +
  theme(legend.position = "none") +
  geom_label(aes(label = mean_potato),
             size = 7.5, position = position_stack(vjust = 0.5),
             show.legend = F, fill = "white")
  
plt3

統計的有意性表示する場合

• 関数とデータの準備

• 平均値と95％信頼区間　mean_ci()　を定義する関数を作成する

mean_ci <- function(data, by, vari){
    se <- function(x) sqrt(var(x)/length(x))
    meanci <- data %>% 
        group_by({{by}}) %>%
        summarise(n = n(),
                  mean_out = mean({{vari}}),
                  se_out = se({{vari}}),
                  .group = "drop"
                  ) %>%
        mutate(
            lwr = mean_out - 1.96 * se_out,
            upr = mean_out + 1.96 * se_out
            ) %>%
        mutate(across(where(is.double), round, 1)) %>%
        mutate(mean_label = format(round(mean_out, 1), nsmall = 1)) %>% 
        select({{by}}, mean_out, lwr, upr, mean_label) %>% 
        mutate(across(.cols = {{by}}, as.factor))
    return(meanci)
}

• 上で作成した関数を使って、ポテトの点数の平均と95%信頼区間を計算してみる

potato_mean <- df_menu %>% 
  mean_ci(mosmc, fried_potato)

potato_mean

# A tibble: 2 × 5
  mosmc mean_out   lwr   upr mean_label
  <fct>    <dbl> <dbl> <dbl> <chr>     
1 mc        80.9  79.4  82.4 80.9      
2 mos       74.7  72.4  77   74.7      

• lwr は 95% 信頼区間の下側
  

• upr は 95% 信頼区間の上側
  

• mean_label はグラフに貼り付けるラベル

• ここで得られた \(t\) 検定結果を broom::tidy() 関数を使って tibble 形式にする
  

• 可視化する際には \(p\) 値が必要　→　　\(p\) 値を可視化する条件を指定
  

• 再度、t 検定を実行
  

• ここではロング型データ (df_menu) を使っている

potato_ttest <- t.test(fried_potato ~ mosmc, data = df_menu)

potato_tidy <- tidy(potato_ttest) %>% 
  select(estimate, p.value, conf.low, conf.high) %>%  
  mutate(
    p_label = case_when(                            # p値のラベル指定  
      p.value <= 0.01 ~ "p < .01",                  # p値が <= 0.01→ < .01 と表示
      p.value > 0.01 & p.value <= 0.05 ~ "p < .05", # p値が0.01と0.05の間→ p < .05 と表示（5%で有意）
      p.value > 0.05 & p.value <= 0.1 ~ "p < .1",   # p値が0.05と0.1の間→ p < .1 と表示（10%で有意）
      p.value > 0.1 ~ "N.S"                         # p値が > 0.1→ N.S (Not significant:統計的に有意ではない) と表示
      )
    )

• 上で作成した ポテト点数の平均と95%信頼区間データ potato_mean と potato_tidy を bind_cols()関数を使って結合する
  →　df_potatoと名前をつける
  
• mosmc という新たな変数を作る
  
• p_label という新たな変数を作る

df_potato <- bind_cols(potato_mean, potato_tidy) %>% 
  mutate(
    mosmc = as.factor(if_else(mosmc == "mc", 
                              "マクドナルド", "モスバーガー")),
    p_label = if_else(mosmc == "マクドナルド", p_label, NA_character_),
    menu = "ポテト"
    )

• df_potato の中を表示する

df_potato

# A tibble: 2 × 11
  mosmc     mean_out   lwr   upr mean_label estimate  p.value conf.low conf.high
  <fct>        <dbl> <dbl> <dbl> <chr>         <dbl>    <dbl>    <dbl>     <dbl>
1 マクドナ…     80.9  79.4  82.4 80.9            6.2 0.000449     3.23      9.17
2 モスバー…     74.7  72.4  77   74.7            6.2 0.000449     3.23      9.17
# ℹ 2 more variables: p_label <chr>, menu <chr>

• 得られた結果を棒グラフで表示させる

pl_potato <- df_potato %>% 
  ggplot(aes(x = mosmc, y = mean_out, fill = mosmc)) +
  geom_bar(stat = "identity") +
  geom_errorbar(aes(ymin = lwr, ymax = upr, width = 0.3)) +
  geom_label(aes(label = mean_label),
             size = 7.5, position = position_stack(vjust = 0.5),
             show.legend = F, fill = "white") +
  geom_segment(aes(x = 1, y = 90, xend = 1, yend = 95)) +
  geom_segment(aes(x = 1, y = 95, xend = 2, yend = 95)) +
  geom_segment(aes(x = 2, y = 90, xend = 2, yend = 95)) +
  geom_text(aes(x = 1.5, y = 100, label = p_label), 
            size = 4.5, family = "Times New Roman", inherit.aes = FALSE) +
  scale_fill_manual(values = c("red", "green4")) +
  scale_y_continuous(expand = c(0, 0),
                     limits = c(0, 105)) +
  labs(x = NULL, y = "評価点の平均値",
       title = "ポテトの評価点の平均値の比較") +
  theme(legend.position = "none",
        plot.title = element_text(size = 12, hjust = 0.5),
        axis.title = element_text(size = 13),
        axis.text = element_text(size = 13))

pl_potato

6.4 「差分」を表示する方法

• 2 群間の平均値の差を比較をする際、両者の「差分」の信頼区間を示すことが重要
  
• それを可能にする可視化が「差分」を表示する方法
  
• 「差分」を可視化するために必要な変数を計算する

df_potato <- df_potato %>% 
  mutate(across(where(is.double), ~ round(.x, 1))) %>% 
  mutate(
    diff_x = "差分", # 差分を表す変数 `diff_x` を作る  
    diff_label = format(round(estimate, 1), # 小数点1位だけ表記  
                        nsmall = 1)
    )

• 計算した結果を表示させる　　　

df_potato

# A tibble: 2 × 13
  mosmc     mean_out   lwr   upr mean_label estimate p.value conf.low conf.high
  <fct>        <dbl> <dbl> <dbl> <chr>         <dbl>   <dbl>    <dbl>     <dbl>
1 マクドナ…     80.9  79.4  82.4 80.9            6.2       0      3.2       9.2
2 モスバー…     74.7  72.4  77   74.7            6.2       0      3.2       9.2
# ℹ 4 more variables: p_label <chr>, menu <chr>, diff_x <chr>, diff_label <chr>

• モスバーガーとマクドナルドの平均点を可視化する
  
• その際、95％信頼区間も表示する

pl_mean_poteto <- df_potato %>% 
  ggplot(aes(x = mosmc, 
             y = mean_out, 
             ymin = lwr, 
             ymax = upr)) +
  geom_pointrange(size = 1) +
  geom_text(aes(label = mean_label), 
            size = 6.5, 
            nudge_x = .13) + # ドット（●） と表記された数値との距離 
  ylim(70, 100) +
  labs(x = NULL, y = NULL, title = "評価点の平均値") +
  theme(plot.title  = element_text(hjust = 0.5, size = 16),
        axis.text = element_text(size = 17),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.background = element_blank(),
        strip.text.y = element_blank())

• モスバーガーとマクドナルドの平均点の「差分」を可視化する

pl_diff_poteto <- df_potato %>% 
  ggplot(aes(x = diff_x, y = estimate)) +
  geom_hline(yintercept = 0, col = "red") +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), size = 1) +
  geom_text(aes(label = diff_label), 
            size = 6.5, 
            nudge_x = .19) +　# ドット（●） と表記された数値との距離 
  labs(x = NULL, y = NULL, title = "差分") +
  theme(plot.title  = element_text(hjust = 0.5, size = 16),
        axis.text = element_text(size = 17),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        strip.text.y = element_text(size = 17),
        strip.background = element_blank())

• patchwork パッケージを使って作成した図を並べて表示させる

pl_mean_diff <- pl_mean_poteto + pl_diff_poteto + plot_layout(widths = c(3, 1))

pl_mean_diff

6.5 箱ひげ図＋バイオリン図（参考）

df_menu %>% 
    ggplot(aes(mosmc, fried_potato, fill = mosmc)) + 
    geom_violin() +
  geom_boxplot(width = .1) + # 箱ひげ図の幅を 0.1 と指定
    stat_summary(fun.y = mean, geom = "point") + # 平均値を点で示す 
  scale_x_discrete(labels = c( "モスバーガー", "マクドナルド")) +
  labs(x = "店名", y = "評価") +
    ggsignif::geom_signif(comparisons = combn(sort(unique(df_menu$mosmc)), 2, FUN = list),
                          test = "t.test", na.rm = T,
                          step_increase = 0.1)

7. Exercise

Exercise 7.1

総選挙データ (hr96-21.csv) を使って、2021年総選挙（小選挙区）における自民党と公明党の得票率の平均値に差があるかどうかを知りたい

Q1: この検定における帰無仮説を書きなさい
Q2: この検定における対立仮説を書きなさい
Q3: t.test()関数を使って検定結果を出力し、その結果をわかりやすい言葉で説明しなさい
Q4: 統計的有意性を含む分析結果をgeom_signif() function を作って箱ひげ図＋バイオリン図で示しなさい
Q5: \(t\) 検定の結果を棒グラフを使って示しなさい
Q6: \(t\) 検定の結果を「差分」を使って示しなさい

分析で使うcsvファイル：hr96-21.csv

Exercise 7.2

「5. 有意水準に関するシミュレーション」を参照して有意水準を 1% (α = 0.01) に変更すると、帰無仮説を棄却する割合がどう変化するかグラフで示しなさい

Exercise 7.3

モスバーガーとマクドナルドそれぞれの店で「チーズバーガー」「フライドポテト」「テリヤキバーガー」「シェイク」を食べてもらい点数をつけたという架空データがある

• 20人の回答者はモスバーガーもしくはマクドナルドのいずれかで 4 種類食べたとする
  
• ここで知りたいのは「チーズバーガー」は Mos と Mc どちらがおいしいかということ

Q1: この検定における帰無仮説を書きなさい
Q2: この検定における対立仮説を書きなさい
Q3: t.test()関数を使って検定結果を出力し、その結果をわかりやすい言葉で説明しなさい
Q4: 統計的有意性を含む分析結果をgeom_signif() function を作って箱ひげ図＋バイオリン図で示しなさい
Q5: \(t\) 検定の結果を棒グラフを使って示しなさい
Q6: \(t\) 検定の結果を「差分」を使って示しなさい

分析で使うcsvファイル：menu.csv

Exercise 7.4

次のデータは「計量分析（政治）」の受講生30名に対して行った試験結果（架空データ）

• 授業の初日に行った試験結果は「before」、授業最終日に行った試験結果は「after」に示している
  
• ここでは \(t\) 検定を行うことによって「計量政治学」の授業を受けたことによって、計量政治学の試験スコアーが良くなったかどうか、ということを知りたい

Q1: この検定における帰無仮説を書きなさい
Q2: この検定における対立仮説を書きなさい
Q3: t.test()関数を使って検定結果を出力し、その結果をわかりやすい言葉で説明しなさい
Q4: 統計的有意性を含む分析結果をgeom_signif() function を作って箱ひげ図＋バイオリン図で示しなさい
Q5: \(t\) 検定の結果を棒グラフを使って示しなさい
Q6: \(t\) 検定の結果を「差分」を使って示しなさい

分析で使うcsvファイル：test_score.csv

Exercise 7.5

2021年総選挙において「自民党」と「立憲民主党」の候補者の年齢 (age) に差があるかどうか調べたい

Q1: この検定における帰無仮説を書きなさい
Q2: この検定における対立仮説を書きなさい
Q3: t.test()関数を使って検定結果を出力し、その結果をわかりやすい言葉で説明しなさい
Q4: 統計的有意性を含む分析結果をgeom_signif() function を作って箱ひげ図＋バイオリン図で示しなさい
Q5: \(t\) 検定の結果を棒グラフを使って示しなさい
Q6: \(t\) 検定の結果を「差分」を使って示しなさい

分析で使うcsvファイル：hr96-21.csv

DT::datatable(df_long)

DT::datatable(df_menu)

参考文献

宋財泫 (Jaehyun Song)- 矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」

グラフ作成に関しては遠藤勇哉氏（東北大学大学院情報科学研究科博士後期課程）の助言を参考にしています

土井翔平（北海道大学公共政策大学院）「Rで計量政治学入門」

矢内勇生（高知工科大学）授業一覧

浅野正彦, 矢内勇生.『Rによる計量政治学』オーム社、2018年

浅野正彦, 中村公亮.『初めてのRStudio』オーム社、2018年

Winston Chang, R Graphics Cookbook, O’Reilly Media, 2012.

Kieran Healy, DATA VISUALIZATION, Princeton, 2019

Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017