library(conflicted) # 関数名の衝突を警告
conflicts_prefer( # 優先的に使う関数を指定
dplyr::filter(),
dplyr::select(),
dplyr::lag(),
)
library(tidyverse)
#' 日本語を用いるので macOS ではフォントの設定を行う
if(Sys.info()["sysname"] == "Darwin") { # macOS か調べて日本語フォントを指定
theme_update(text = element_text(family = "HiraginoSans-W4"))
update_geom_defaults("text", list(family = theme_get()$text$family))
update_geom_defaults("label", list(family = theme_get()$text$family))}第12講 分散分析
複数のグループ間の違いを検証する
準備
以下で利用する共通パッケージを読み込む．
一元配置分散分析
問題
一元配置分散分析について確率シミュレーションを行いなさい．
解答欄
解答例
人工データにより一元配置分散分析の Monte-Carlo 実験を行う．
A,B,Cの3つの水準をもつデータを考える．
fact <- as.factor(rep(LETTERS[1:3], c(10,8,12))) # 因子(水準A,B,Cが10,8,12個)
toy_data0 <- # 平均が全て等しい場合 (帰無仮説が正しい例)
tibble(factor = fact, # 水準
value = rnorm(length(fact), mean = 3, sd = 2)) # 観測値データの状況を図示する．
toy_data0 |>
ggplot(aes(x = factor, y = value)) +
geom_boxplot(fill = "pink") +
geom_jitter(colour = "red", width = 0.2) +
labs(title = "帰無仮説が正しい場合")
水準ごとの平均にずれのある対立仮説の例を作成する．
alt <- rep(c(1,-1,1), c(10,8,12)) # 水準ごとの平均のずれ
toy_data1 <- # 平均が水準ごとに異なる場合(対立仮説が正しい例)
tibble(factor = fact, # 水準
value = rnorm(length(fact), mean = 3, sd = 2) + alt) # 観測値同様に図示する．
toy_data1 |>
ggplot(aes(x = factor, y = value)) +
geom_boxplot(fill = "lightblue") +
geom_jitter(colour = "blue", width = 0.2) +
labs(title = "帰無仮説が誤りの場合")
これらの2つのデータに対して分散分析を行う． 帰無仮説が正しい場合は以下のようになる．
aov(value ~ factor, data = toy_data0) Call:
aov(formula = value ~ factor, data = toy_data0)
Terms:
factor Residuals
Sum of Squares 1.21898 94.95877
Deg. of Freedom 2 27
Residual standard error: 1.875364
Estimated effects may be unbalancedaov(value ~ factor, data = toy_data0) |>
broom::tidy() # 分散分析表# A tibble: 2 × 6
term df sumsq meansq statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 factor 2 1.22 0.609 0.173 0.842
2 Residuals 27 95.0 3.52 NA NA 帰無仮説が誤り，すなわち対立仮説が正しい場合は以下のようになる．
aov(value ~ factor, data = toy_data1) Call:
aov(formula = value ~ factor, data = toy_data1)
Terms:
factor Residuals
Sum of Squares 40.98323 59.51545
Deg. of Freedom 2 27
Residual standard error: 1.48468
Estimated effects may be unbalancedaov(value ~ factor, data = toy_data1) |>
broom::tidy() # 分散分析表# A tibble: 2 × 6
term df sumsq meansq statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 factor 2 41.0 20.5 9.30 0.000848
2 Residuals 27 59.5 2.20 NA NA
Note
関数 aov() のようにデータが第1引数でない場合にパイプ演算子を使うには _ (placeholder) を使って位置を明示すればよい．
toy_data0 |> aov(value ~ factor, data = _) Call:
aov(formula = value ~ factor, data = toy_data0)
Terms:
factor Residuals
Sum of Squares 1.21898 94.95877
Deg. of Freedom 2 27
Residual standard error: 1.875364
Estimated effects may be unbalancedtoy_data0 |>
aov(value ~ factor, data = _) |>
broom::tidy()# A tibble: 2 × 6
term df sumsq meansq statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 factor 2 1.22 0.609 0.173 0.842
2 Residuals 27 95.0 3.52 NA NA Monte-Carlo 実験を行う． まず，試行のための関数を作成する．
mc_trial <- function(h0 = TRUE){
if(h0) { # 帰無仮説が正しい場合
tmp <-
tibble(factor = fact,
value = rnorm(length(fact), mean = 3, sd = 2)) |>
aov(value ~ factor, data = _) |> broom::tidy()
} else { # 対立仮説が正しい場合
tmp <-
tibble(factor = fact, # 水準ごとに平均が異なる
value = rnorm(length(fact), mean = 3, sd = 2) + alt) |>
aov(value ~ factor, data = _) |> broom::tidy()
}
#' 分散分析表における factor 行の統計量とp値を返す
return(c(statistic = tmp[1,5, drop = TRUE],
p.value = tmp[1,6, drop = TRUE]))
#' 以下のように書いても良い
#' return(tmp |> select(5:6) |> t() |> _[,1])
}帰無仮説が正しい場合のF統計量およびp値の分布を調べる． この場合，p値は一様分布になることに注意する．
replicate(2000, mc_trial()) |>
t() |> as_tibble() |> pivot_longer(everything()) |> # 描画用に変形
ggplot(aes(x = value, fill = name)) +
geom_histogram(show.legend = FALSE) +
facet_wrap(vars(name), ncol = 2, scales = "free") + # それぞれ個別の座標で描画
labs(title = "帰無仮説が正しい場合")
対立仮説が正しい場合のF統計量およびp値の分布を調べる． この場合，p値は小さな値に偏ることが確認できる．
replicate(2000, mc_trial(h0 = FALSE)) |>
t() |> as_tibble() |> pivot_longer(everything()) |>
ggplot(aes(x = value, fill = name)) +
geom_histogram(show.legend = FALSE) +
facet_wrap(vars(name), ncol = 2, scales = "free") +
labs(title = "帰無仮説が誤りの場合")
問題
東京の気候データ (tokyo_weather.csv) の 気温の項目について以下の問に答えよ．
- 曜日ごとの気温の平均と分散を求めよ．
- 曜日ごとに平均が異なるといえるかどうか分散分析を用いて検定しなさい．
Note
週日を因子(factor)とするには例えば以下のような項目を加えればよい．
tw_data <- read_csv("data/tokyo_weather.csv") |>
mutate(day_of_week = ordered(day_of_week, # 曜日を順序付き因子にする
levels = wday(1:7, label = TRUE))) 関数 lubridate::wday() の出力は使用言語に依存するので， 強制的に英語にするには locale を指定して wday(..., locale = "en_US") とすれば良い． 省略形を用いない場合は abbr を指定する． 以下はドイツ語の例．
wday(1:7, label = TRUE, abbr = FALSE, locale = "de_DE")利用可能な言語を調べるには stringi::stri_locale_list() を用いる．
解答欄
解答例
気候データにおいて 曜日ごとの気温に差があるか否かを分散分析を用いて調べる．
データを読み込み，整理する．
tw_data <- read_csv("data/tokyo_weather.csv") |>
mutate(day_of_week = ordered(day_of_week, # 曜日を順序付き因子にする
levels = wday(1:7, label = TRUE))) 曜日ごとの分布を見るとともに平均・分散(標準偏差)を確認する．
tw_data |>
ggplot(aes(x = day_of_week, y = temp)) +
geom_boxplot(fill = "lavender") +
labs(x = "曜日", y = "気温", title = "曜日ごとの気温の分布") 
tw_data |>
group_by(day_of_week) |>
summarize(across(temp, list(mean = mean,
var = var,
sd = sd)))# A tibble: 7 × 4
day_of_week temp_mean temp_var temp_sd
<ord> <dbl> <dbl> <dbl>
1 Sun 17.6 78.0 8.83
2 Mon 17.8 75.2 8.67
3 Tue 17.4 79.3 8.91
4 Wed 17.2 79.1 8.89
5 Thu 17.5 71.0 8.43
6 Fri 17.1 65.8 8.11
7 Sat 16.9 76.7 8.76曜日ごとの気温差に関する分散分析を行う．
aov(temp ~ day_of_week, data = tw_data) # aovによる分析Call:
aov(formula = temp ~ day_of_week, data = tw_data)
Terms:
day_of_week Residuals
Sum of Squares 29.139 26859.345
Deg. of Freedom 6 358
Residual standard error: 8.661761
Estimated effects may be unbalancedaov(temp ~ day_of_week, data = tw_data) |>
summary() # 分散分析表の表示(棄却されない) Df Sum Sq Mean Sq F value Pr(>F)
day_of_week 6 29 4.86 0.065 0.999
Residuals 358 26859 75.03 aov(temp ~ day_of_week, data = tw_data) |>
broom::tidy() # 分散分析表のデータフレーム (行・列の番号で必要な情報が取得できる)# A tibble: 2 × 6
term df sumsq meansq statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 day_of_week 6 29.1 4.86 0.0647 0.999
2 Residuals 358 26859. 75.0 NA NA aov(temp ~ day_of_week, data = tw_data) |>
model.tables(type = "means") # 水準(曜日)ごとの平均値Tables of means
Grand mean
17.36877
day_of_week
Sun Mon Tue Wed Thu Fri Sat
17.62 17.77 17.42 17.23 17.54 17.11 16.9
rep 52.00 52.00 52.00 53.00 52.00 52.00 52.0aov(temp ~ day_of_week, data = tw_data) |>
model.tables(type = "effects") # 水準(曜日)ごとの効果Tables of effects
day_of_week
Sun Mon Tue Wed Thu Fri Sat
0.2466 0.3985 0.05046 -0.1424 0.1755 -0.2553 -0.4707
rep 52.0000 52.0000 52.00000 53.0000 52.0000 52.0000 52.0000検定のみ実行する場合は，関数 oneway.test() を用いることができる．
級内の分散が等しいことを仮定する場合は以下のように指定する．
oneway.test(temp ~ day_of_week, data = tw_data, var.equal = TRUE)
One-way analysis of means
data: temp and day_of_week
F = 0.064732, num df = 6, denom df = 358, p-value = 0.9989級内の分散が等しいことが仮定できない場合はWelchの近似法による検定が行われる．
oneway.test(temp ~ day_of_week, data = tw_data)
One-way analysis of means (not assuming equal variances)
data: temp and day_of_week
F = 0.063645, num df = 6.00, denom df = 159.09, p-value = 0.999月ごとの気温に差があるか否かを分散分析で検証する． 図示した結果から棄却されることが予想される．
tw_data <- tw_data |> mutate(month = as_factor(month)) # 月を因子化
tw_data |>
ggplot(aes(x = month, y = temp)) +
geom_boxplot(fill = "lavender") +
labs(x = "月", y = "気温", title = "月ごとの気温の分布") 
aov(temp ~ month, data = tw_data) |>
summary() # 分散分析表の表示(棄却される) Df Sum Sq Mean Sq F value Pr(>F)
month 11 24132 2193.8 280.9 <2e-16 ***
Residuals 353 2757 7.8
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Note
曜日の因子化などは package::lubridate に含まれる関数が利用できる．
曜日の文字列を得る方法には以下のものがある．
wday(1:7, label = TRUE) # 省略名[1] Sun Mon Tue Wed Thu Fri Sat
Levels: Sun < Mon < Tue < Wed < Thu < Fri < Satwday(1:7, label = TRUE, abbr = FALSE) # 正式名[1] Sunday Monday Tuesday Wednesday Thursday Friday Saturday
7 Levels: Sunday < Monday < Tuesday < Wednesday < Thursday < ... < Saturdaywday(2:8, label = TRUE) # 月曜から始まる文字列を得る場合[1] Mon Tue Wed Thu Fri Sat <NA>
Levels: Sun < Mon < Tue < Wed < Thu < Fri < Satwday(1:7, label = TRUE, week_start = 1) # 順序を月曜から始める場合[1] Sun Mon Tue Wed Thu Fri Sat
Levels: Mon < Tue < Wed < Thu < Fri < Sat < Sun日付から曜日を得る方法は例えば以下のようにすれば良い．
tw_data |>
mutate(date = as_date(paste(year,month,day, sep = "-")), # 日付に変換
weekday = wday(date, label = TRUE), # 日付から曜日(順序付因子)を追加
.before = year) # year列の前に付加# A tibble: 365 × 15
date weekday year month day day_of_week temp rain solar snow
<date> <ord> <dbl> <fct> <dbl> <ord> <dbl> <dbl> <dbl> <dbl>
1 2025-01-01 Wed 2025 1 1 Wed 7.1 0 11.4 0
2 2025-01-02 Thu 2025 1 2 Thu 7.8 0 11.2 0
3 2025-01-03 Fri 2025 1 3 Fri 5.5 0 3.97 0
4 2025-01-04 Sat 2025 1 4 Sat 5.1 0 11.8 0
5 2025-01-05 Sun 2025 1 5 Sun 4.6 0 10.4 0
6 2025-01-06 Mon 2025 1 6 Mon 4.5 25 4.26 0
7 2025-01-07 Tue 2025 1 7 Tue 8.1 0 7.34 0
8 2025-01-08 Wed 2025 1 8 Wed 6.6 0 10.1 0
9 2025-01-09 Thu 2025 1 9 Thu 6.7 0 10.8 0
10 2025-01-10 Fri 2025 1 10 Fri 4.6 0 12.5 0
# ℹ 355 more rows
# ℹ 5 more variables: wdir <chr>, wind <dbl>, press <dbl>, humid <dbl>,
# cloud <dbl>二元配置分散分析
問題
package::datarium に含まれている jobsatisfaction データ (性別 (gender) と 学歴 (education_level) の違いによる 仕事の満足度 (score) を収集したデータ) について以下の問に答えよ．
- データを適当な方法で可視化しなさい．
- それぞれの因子で満足度の平均に違いがあるか，二元配置の分散分析を用いて検討しなさい．
Note
パッケージのインストールと読み込みは以下のように行うことができる．
install.packages("datarium") # インストール
library("datarium") # 読み込み解答欄
解答例
datarium::jobsatisfaction は 性別(gender)と学歴(education_level)による仕事の満足度(score)を集計したデータである． 二元配置分散分析を用いて満足度にそれぞれの因子の効果があるかを検証する．
必要であれば package::datarium をインストールし， パッケージの読み込みを行う．
#' install.packages("datarium") # 検定用のデータが集められている
library("datarium") # packageの読み込み
jobsatisfaction # データの内容を確認# A tibble: 58 × 4
id gender education_level score
<fct> <fct> <fct> <dbl>
1 1 male school 5.51
2 2 male school 5.65
3 3 male school 5.07
4 4 male school 5.51
5 5 male school 5.94
6 6 male school 5.8
7 7 male school 5.22
8 8 male school 5.36
9 9 male school 4.78
10 10 male college 6.01
# ℹ 48 more rowsデータの視覚化を行う．
jobsatisfaction |> # 学歴ごとに性別を比較
ggplot(aes(x = education_level, y = score, fill = gender)) +
geom_boxplot()
jobsatisfaction |> # 性別ごとに学歴を比較
ggplot(aes(x = gender, y = score, fill = education_level)) +
geom_boxplot()
二元配置分散分析を行い，結果をいくつかの方法で表示する．
aov(score ~ gender + education_level, data=jobsatisfaction) |>
summary() Df Sum Sq Mean Sq F value Pr(>F)
gender 1 0.54 0.54 1.447 0.234
education_level 2 113.68 56.84 152.172 <2e-16 ***
Residuals 54 20.17 0.37
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1aov(score ~ gender + education_level, data=jobsatisfaction) |>
model.tables(type = "means") Tables of means
Grand mean
6.963276
gender
male female
7.063 6.87
rep 28.000 30.00
education_level
school college university
5.594 6.351 8.846
rep 19.000 19.000 20.000aov(score ~ gender + education_level, data=jobsatisfaction) |>
model.tables(type = "effects")Tables of effects
gender
male female
0.09994 -0.09328
rep 28.00000 30.00000
education_level
school college university
-1.369 -0.612 1.882
rep 19.000 19.000 20.000