R パッケージ
一覧library(tidyverse)
library(broom)
library(patchwork)
library(DT)
library(ggbeeswarm)
library(ggsignif)
library(rcompanion)
library(rmarkdown)
1) 標本平均が特定の値かどうか調べる
2) 二つの異なるグループの平均が異なるかどうか調べる
student's t
という名前で論文を投稿したため、\(t\) 検定と呼ばれる Paired
と Unpaired
データの二種類のデータに関して \(t\) 検定の方法を演習する
結果の示し方 | 特徴 |
---|---|
数値だけ | 2 群の差が分かりにくい |
箱ひげ図・バイオリン図 | 2 群の差が可視化されて見やすく、統計的有意性が明確 |
棒グラフ | 2 群の差が可視化されて見やすく(平均値も表記)、統計的有意性が明確 |
「差分」を表示 | 2 群の差が可視化されて見やすく(平均値も表記)、95%有意水準が明確 |
参考:近年の実証研究では「棒グラフ」もしくは「差分」を表示する方法が使われることが多い
Paired-samples t-test
)Paired data
)データの読み取りと準備 (mos_mc.csv
)
<- read_csv("data/mos_mc.csv") df_mos_mc
::datatable(df_mos_mc) DT
wide format
)」と呼ばれ、R 上では分析しにくいtidyr::pivot_longer()
関数を使って「ロング形式 (long format
)」に変換する<- df_mos_mc %>%
df_long ::pivot_longer(mos:mc,
tidyrnames_to = "burger", # バーガー店名を入力
values_to = "score") # バーガーの評価点を入力
::datatable(df_long) DT
文字化けしないためのコマンド(マックユーザのみ)
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 = "評価")
%>%
df_long mutate(burger = fct_inorder(burger)) %>%
ggplot(aes(x = burger, y = score)) +
geom_violin() +
scale_x_discrete(labels = c( "モスバーガー", "マクドナルド")) +
labs(x = "店名", y = "評価")
%>%
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)
ID mos mc
Min. : 1.00 Min. :70.00 Min. :70.00
1st Qu.: 3.25 1st Qu.:76.25 1st Qu.:75.00
Median : 5.50 Median :80.00 Median :80.00
Mean : 5.50 Mean :80.50 Mean :79.50
3rd Qu.: 7.75 3rd Qu.:85.00 3rd Qu.:83.75
Max. :10.00 Max. :90.00 Max. :90.00
ttest
を実行してみる(paired ttest
) 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 と名付ける
<- df_mos_mc$mos
x <- df_mos_mc$mc y
<- x - y d
<- (mean(d) - 0) / (sd(d) / sqrt(10)) t
t
[1] 0.557086
df_long
# A tibble: 20 x 3
ID burger score
<dbl> <chr> <dbl>
1 1 mos 80
2 1 mc 75
3 2 mos 75
4 2 mc 70
5 3 mos 80
6 3 mc 80
7 4 mos 90
8 4 mc 85
9 5 mos 85
10 5 mc 90
11 6 mos 80
12 6 mc 75
13 7 mos 75
14 7 mc 85
15 8 mos 85
16 8 mc 80
17 9 mos 85
18 9 mc 80
19 10 mos 70
20 10 mc 75
%>%
df_long ggplot() +
geom_boxplot(aes(x = burger,
y = score))
t.test(df_long$score[df_long$burger == "mos"],
$score[df_long$burger == "mc"],
df_longpaired = 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 difference in means is not equal to 0
95 percent confidence interval:
-3.060696 5.060696
sample estimates:
mean of the differences
1
t = 0.55709
が示されているp-value = 0.5911
が 0.05 よりも大きいggsignif()
関数では \(t\) 検定は unpaired data
がデフォルトtest.args = list(paired = TRUE)
と指定することunpaired data
を使って \(t\) 検定するときには、この指定は不要%>%
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 = "評価") +
::geom_signif(comparisons = combn(sort(unique(df_long$burger)), 2, FUN = list),
ggsigniftest = "t.test",
test.args = list(paired = TRUE), # paired data の $t$ 検定だと指定
na.rm = T,
step_increase = 0.1)
df_mos_mc
# A tibble: 10 x 3
ID mos mc
<dbl> <dbl> <dbl>
1 1 80 75
2 2 75 70
3 3 80 80
4 4 90 85
5 5 85 90
6 6 80 75
7 7 75 85
8 8 85 80
9 9 85 80
10 10 70 75
t.test(df_mos_mc$mos,
$mc,
df_mos_mcpaired = 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 difference in means is not equal to 0
95 percent confidence interval:
-3.060696 5.060696
sample estimates:
mean of the differences
1
One Sample t-test
ここで行った ttest
は次のようにしても同じ結果を求めることができる
One Sample t-test
としても ttest
できる
モスバーガーとマクドナルドの平均値 diff
を計算する
diff
の平均値(サンプル平均)を計算する
この統計量を使ってパラメタを推定する
ワイド形式データ (df_mos_mc
) を使って計算してみる
帰無仮説:diff
の平均値は 0 = モスとマックの味に差はない
<- df_mos_mc$mos - df_mos_mc$mc
diff 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
Paired
データに関して \(t\) 検定の方法を演習したUnpaired
データを使って演習するUnpaired-samples t-test
)Unpaired data
(対応のないサンプル)データの種類 | 解説 |
---|---|
Paired (対応あり) |
バーガーを評価した人は 10 人、両方のバーガーを食べた |
Unpaired (対応なし) |
バーガーを評価した人は 20 人、片方のバーガーだけ食べた |
<- read_csv("data/mos_mc.csv") df_mos_mc
::datatable(df_mos_mc) DT
wide format
)」と呼ばれ、R 上では分析しにくいtidyr::pivot_longer()
関数を使って「ロング形式 (long format
)」に変換する<- df_mos_mc %>%
df_long ::pivot_longer(mos:mc,
tidyrnames_to = "burger", # バーガー店名を入力
values_to = "score") # バーガーの評価点を入力
::datatable(df_long) DT
文字化けしないためのコマンド(マックユーザのみ)
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 = "評価")
ttest
を実行してみる(unpaired ttest
) df_long
# A tibble: 20 x 3
ID burger score
<dbl> <chr> <dbl>
1 1 mos 80
2 1 mc 75
3 2 mos 75
4 2 mc 70
5 3 mos 80
6 3 mc 80
7 4 mos 90
8 4 mc 85
9 5 mos 85
10 5 mc 90
11 6 mos 80
12 6 mc 75
13 7 mos 75
14 7 mc 85
15 8 mos 85
16 8 mc 80
17 9 mos 85
18 9 mc 80
19 10 mos 70
20 10 mc 75
t.test(df_long$score[df_long$burger == "mos"],
$score[df_long$burger == "mc"]) # unpaired が default df_long
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
t = 0.37354
が示されているp-value = 0.7131
が 0.05 よりも大きい%>%
df_long ggplot() +
geom_boxplot(aes(x = burger,
y = score))
ggsignif()
関数における \(t\) 検定は unpaired data
がデフォルト%>%
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") + # 平均値を点で示す
::geom_beeswarm() +
ggbeeswarmscale_x_discrete(labels = c( "モスバーガー", "マクドナルド")) +
labs(x = "店名", y = "評価") +
::geom_signif(comparisons = combn(sort(unique(df_long$burger)), 2, FUN = list),
ggsigniftest = "t.test",
na.rm = T,
step_increase = 0.1)
df_mos_mc
# A tibble: 10 x 3
ID mos mc
<dbl> <dbl> <dbl>
1 1 80 75
2 2 75 70
3 3 80 80
4 4 90 85
5 5 85 90
6 6 80 75
7 7 75 85
8 8 85 80
9 9 85 80
10 10 70 75
t.test(df_mos_mc$mos,
$mc) # unpaired が default df_mos_mc
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
hr96-21.csv
)RProject
フォルダ内に data という名称のフォルダを作成するcsv
ファイルがパソコンにダウンロードされ、data 内に自動的に保存されるdownload.file(url = "http://www.ner.takushoku-u.ac.jp/masano/class_material/waseda/keiryo/Data/hr96-21.csv",
destfile = "data/hr96-21.csv")
注意:一度ダウンロードを完了すれば、このコマンドを実行する必要はありません
hr96-21.csv をクリックしてデータをパソコンにダウンロード
RProject
フォルダ内に data という名称のフォルダを作成する
ダウンロードした hr96-21.csv
を手動でRProject
フォルダ内にある data フォルダに入れる
hr96-21.csv
を読み取るna = "."
というコマンドは「欠損値をドットで置き換える」という意味numeric
)」型のデータが「」文字型 (character
)」として認識されるなど、エラーの原因になるため、読み取る時点で事前に対処する<- read_csv("data/hr96-21.csv",
hr na = ".")
locale()
関数を使って日本語エンコーディング (cp932
) を指定する<- read_csv("data/hr96-21.csv",
hr na = ".",
locale = locale(encoding = "cp932"))
<- read.csv("data/hr96-21.csv",
hr na = ".")
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)
spec_tbl_df [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()
.. )
numeric
文字は character
として認識されていることがわかる2017年総選挙での自民党と立憲民主党が得た得票率の箱ひげ図を描いてみよう
次の 2 つの条件を指定して hr
からデータを抜き出す
ldp_cdp17
と何前をつける<- hr %>%
ldp_cdp17 ::filter(year == 2017) %>% #2017年に実施された総選挙だけ
dplyr::filter(seito == "自民" | seito == "立憲") %>% #自民党と立憲民主党だけ
dplyr::select(year, voteshare, seito) #三つの変数だけ dplyr
::datatable(ldp_cdp17) DT
%>%
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 = "得票率 (%)")
unpaired ttest
を実行unpaired ttest
がデフォルトpaired = FALSE
という指定は不要t.test(ldp_cdp17$voteshare[ldp_cdp17$seito == "自民"],
$voteshare[ldp_cdp17$seito == "立憲"]) ldp_cdp17
Welch Two Sample t-test
data: ldp_cdp17$voteshare[ldp_cdp17$seito == "自民"] and ldp_cdp17$voteshare[ldp_cdp17$seito == "立憲"]
t = 9.3267, df = 88.101, p-value = 8.633e-15
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
11.61689 17.90779
sample estimates:
mean of x mean of y
50.24679 35.48444
解釈 ・帰無仮説:
「立憲民主党も自民党も得票率に差はない」
・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%)
%>%
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 = combn(sort(unique(ldp_cdp17$seito)), 2, FUN = list),
ggsigniftest = "t.test", na.rm = T,
step_increase = 0.1)+
labs(x = "政党名", y = "得票率 (%)")
hr
からデータを抜き出すldp_cdp_poh17
と何前をつける<- hr %>%
ldp_cdp_poh17 filter(year == 2017) %>%
filter(seito == "自民" | seito == "立憲" | seito == "希望") %>%
select(voteshare, seito)
ldp_cdp_poh17
が含む変数 seito
の中身を表示させるunique(ldp_cdp_poh17$seito)
[1] "自民" "立憲" "希望"
%>%
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 = combn(sort(unique(ldp_cdp_poh17$seito)), 2, FUN = list),
ggsigniftest = "t.test", na.rm = T,
step_increase = 0.1)
unpaired data
なので「対応のない \(t\) 検定」を実行するp value = 5.5e-15
# unpaired が default(つまり unpaired がデフォルト)
t.test(ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "自民"],
$voteshare[ldp_cdp_poh17$seito == "立憲"]) ldp_cdp_poh17
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.3267, df = 88.101, p-value = 8.633e-15
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
11.61689 17.90779
sample estimates:
mean of x mean of y
50.24679 35.48444
p value = 2.22e-16
# unpaired が default
t.test(ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "自民"],
$voteshare[ldp_cdp_poh17$seito == "希望"]) ldp_cdp_poh17
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.502, df = 404.26, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
18.33860 22.45015
sample estimates:
mean of x mean of y
50.24679 29.85241
p value = 0.00087
# unpaired が default
t.test(ldp_cdp_poh17$voteshare[ldp_cdp_poh17$seito == "立憲"],
$voteshare[ldp_cdp_poh17$seito == "希望"]) ldp_cdp_poh17
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
<- function(alpha = .05, n = 50, trials = 100) {
sig_sim ## Arguments:
## alpha = significance level(有意水準)、defaultで 5% (α = 0.05)に設定
## n = 標本サイズ、defaultで 50 に設定
## trials = シミュレーションの試行数, defaultで 100 に設定
## Return:
## 帰無仮説を棄却する割合 (rejection rate) を表示させる
## t 分布もヒストグラムとして表示させる
## vector to save the result
<- rep(NA, trials)
res <-rep(NA, trials)
T_vec
## critical value
<- abs(qt(alpha / 2, df = n -1))
cv
for (i in 1:trials) {
<- rnorm(50, mean = 50, sd = 5) # 標本 x は母平均 50、母標準偏差 5 から 50 個無作為抽出
x <- rnorm(50, mean = 50, sd = 5) # 標本 y は母平均 50、母標準偏差 5 から 50 個無作為抽出
y <- x - y # x と y の差を d とする
d <- (mean(d) - 0) / (sd(d) / sqrt(n)) # サンプルから得られる $t$値を計算
T <- T
T_vec[i] <- abs(T) > cv # $t$値の絶対値が cv(棄却限界値)を超える場合を選択
res[i]
}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.0527
上の数値は、帰無仮説を棄却する割合 (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.0973
rejection rate
)が約 10% (0.1に近い値) まで増えていることがわかる
近年、様々な学会誌 (Academic Journal
) では \(t\) 検定の結果が様々な方法で可視化
ここでは棒グラフを使った \(t\) 検定の可視化方法を紹介
ここではモスバーガーとマクドナルドそれぞれの店で「チーズバーガー」「フライドポテト」「テリヤキバーガー」「シェイク」を食べてもらい点数をつけたという架空データを使う
回答者はモスバーガーもしくはマクドナルドのいずれかで4種類食べたとする
ここで知りたいこと
モスとマックの「フライドポテト」はどちらがおいしいか
df_menu
と名前を付ける(架空データ)<- read_csv("data/menu.csv") df_menu
::datatable(df_menu) DT
unpaired data
なので対応のない \(t\) 検定 = Welch の \(t\) 検定を実行<- t.test(fried_potato ~ mosmc, data = df_menu)
potato_ttest 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 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%で統計的に有意)
関数とデータの準備
平均値と95%信頼区間 mean_ci()
を定義する関数を作成する
<- function(data, by, vari){
mean_ci <- function(x) sqrt(var(x)/length(x))
se <- data %>%
meanci 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)
}
<- df_menu %>%
potato_mean mean_ci(mosmc, fried_potato)
potato_mean
# A tibble: 2 x 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
) を使っている
<- t.test(fried_potato ~ mosmc, data = df_menu) potato_ttest
<- tidy(potato_ttest) %>%
potato_tidy select(estimate, p.value, conf.low, conf.high) %>%
mutate(
p_label = case_when( # p値のラベル指定
<= 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:統計的に有意ではない) と表示
p.value
) )
potato_mean
と potato_tidy
を bind_cols()
関数を使って結合するdf_potato
と名前をつけるmosmc
という新たな変数を作るp_label
という新たな変数を作る<- bind_cols(potato_mean, potato_tidy) %>%
df_potato 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 x 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 4.49e-4 3.23 9.17
2 モスバー… 74.7 72.4 77 74.7 6.2 4.49e-4 3.23 9.17
# … with 2 more variables: p_label <chr>, menu <chr>
<- df_potato %>%
pl_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
<- 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 x 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
# … with 4 more variables: p_label <chr>, menu <chr>, diff_x <chr>,
# diff_label <chr>
<- df_potato %>%
pl_mean_poteto 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())
<- df_potato %>%
pl_diff_poteto 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_poteto + pl_diff_poteto + plot_layout(widths = c(3, 1))
pl_mean_diff
pl_mean_diff
%>%
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 = "評価") +
::geom_signif(comparisons = combn(sort(unique(df_menu$mosmc)), 2, FUN = list),
ggsigniftest = "t.test", na.rm = T,
step_increase = 0.1)
総選挙データ (hr96-21.csv
) を使って、2021年総選挙(小選挙区)における自民党と公明党の得票率の平均値に差があるかどうかを知りたい
Q1: この検定における帰無仮説を書きなさい
Q2: この検定における対立仮説を書きなさい
Q3: t.test()
関数を使って検定結果を出力し、その結果をわかりやすい言葉で説明しなさい
Q4: 統計的有意性を含む分析結果をgeom_signif() function
を作って箱ひげ図+バイオリン図で示しなさい
Q5: \(t\) 検定の結果を棒グラフを使って示しなさい
Q6: \(t\) 検定の結果を「差分」を使って示しなさい
✔ 分析で使うcsvファイル:hr96-21.csv
「5. 有意水準に関するシミュレーション
」を参照して有意水準を 1% (α = 0.01
) に変更すると、帰無仮説を棄却する割合がどう変化するかグラフで示しなさい
モスバーガーとマクドナルドそれぞれの店で「チーズバーガー」「フライドポテト」「テリヤキバーガー」「シェイク」を食べてもらい点数をつけたという架空データがある
Q1: この検定における帰無仮説を書きなさい
Q2: この検定における対立仮説を書きなさい
Q3: t.test()
関数を使って検定結果を出力し、その結果をわかりやすい言葉で説明しなさい
Q4: 統計的有意性を含む分析結果をgeom_signif() function
を作って箱ひげ図+バイオリン図で示しなさい
Q5: \(t\) 検定の結果を棒グラフを使って示しなさい
Q6: \(t\) 検定の結果を「差分」を使って示しなさい
✔ 分析で使うcsvファイル:menu.csv
次のデータは「計量分析(政治)」の受講生30名に対して行った試験結果(架空データ)
Q1: この検定における帰無仮説を書きなさい
Q2: この検定における対立仮説を書きなさい
Q3: t.test()
関数を使って検定結果を出力し、その結果をわかりやすい言葉で説明しなさい
Q4: 統計的有意性を含む分析結果をgeom_signif() function
を作って箱ひげ図+バイオリン図で示しなさい
Q5: \(t\) 検定の結果を棒グラフを使って示しなさい
Q6: \(t\) 検定の結果を「差分」を使って示しなさい
✔ 分析で使うcsvファイル:test_score.csv
2021年総選挙において「自民党」と「立憲民主党」の候補者の年齢 (age
) に差があるかどうか調べたい
Q1: この検定における帰無仮説を書きなさい
Q2: この検定における対立仮説を書きなさい
Q3: t.test()
関数を使って検定結果を出力し、その結果をわかりやすい言葉で説明しなさい
Q4: 統計的有意性を含む分析結果をgeom_signif() function
を作って箱ひげ図+バイオリン図で示しなさい
Q5: \(t\) 検定の結果を棒グラフを使って示しなさい
Q6: \(t\) 検定の結果を「差分」を使って示しなさい
✔ 分析で使うcsvファイル:hr96-21.csv
参考文献