install.packages("pacman")Rの復習
1 パッケージ
通常、Rでのパッケージのインストールとアップデートはinstall.packages()関数、読み込みはlibrary()、またはrequire()関数を使う。一つ注意すべき点はinstall.packages()の場合、R公式レポジトリであるCRANに登録されているパッケージのみが対象となっている点だ。しかし、今はCRANでなくGitHub上で公開されているパッケージも非常に多い。これらパッケージは{devtools}か{remote}パッケージを使う。これらの関数を使い分けることは面倒なので、本講義ではこれらの処理を統合した{pacman}パッケージを使用する。まずは、{pacman}パッケージをインストールする。
まず、CRANに登録されているパッケージを読み込む際は、pacman::p_load(読み込むパッケージ名)を入力する1。インストールされていない場合は、自動的にCRANからダウンロード&インストールした上で読み込んでくれるので便利だ2。以下では本講義で使用するパッケージとして{tidyverse}、{summarytools}、{fastDummies}、{modelsummary}、{broom}を読み込む。
# pacman::p_load(tidyverse, summarytools, fastDummies, modelsummary, broom) もOK
pacman::p_load(tidyverse,
summarytools,
fastDummies,
modelsummary,
broom,
cobalt)
ノート参考）CRANでなく、GitHubで配布しているパッケージのインストール
CRANでなく、GitHub上で公開されているパッケージを使う場合はpacman::p_load_gh()を使用する。()の中には"ユーザー名/リポジトリ名"を入力する。たとえば、{BalanceR}の作成者のGitHubアカウント名はJaehyunSongであり、{BalanceR}のリポジトリ名はBalanceRだから、以下のように入力する。p_load()とは違って、文字列は"で囲む必要があることに注意しよう。
pacman::p_load_gh("JaehyunSong/BalanceR")2 データの読み込み
.csv形式のデータを読み込むにはread_csv()関数を使用する3。()内には読み込むファイルのパスを"で囲んで記入する。read_csv()関数はファイルの読み込みのみの機能しか持たない。現在の作業環境内に読み込んだデータを格納するためには代入演算子<-を使う。ここではdataフォルダー内のrct_data.csvを読み込み4、raw_dfという名のオブジェクトとしてく格納する。作業環境内のオブジェクトはRを再起動すると削除されるため、改めてパッケージ・データの読み込みが必要だ。
raw_df <- read_csv("_data/Gerber_et_al_APSR_2008.csv")オブジェクトの中身を出力するためにはオブジェクト名を入力する。
raw_df# A tibble: 5,000 × 8
treatment gender yob hh_size voted2000 voted2002 voted2004 voted2006
<chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr>
1 Control male 1981 3 yes no no no
2 Control female 1957 2 no yes no no
3 Control male 1975 2 no no yes no
4 Hawthorne male 1982 4 no no no no
5 Control female 1949 2 no no yes yes
6 Control female 1937 2 yes yes yes yes
7 Control male 1967 2 no no yes no
8 Control male 1960 2 no no no no
9 Civic Duty male 1952 2 no no yes no
10 Control female 1953 2 no no yes no
# ℹ 4,990 more rows 表形式データの大きさ（行の列の数）の確認にはdim()関数を使う。長さ2のnumeric（数値）型ベクトルが返され、それぞれデータセットの行と列の数を意味する。
dim(raw_df)[1] 5000 8 表形式データの場合、各列には名前が付いており、それぞれが一つの変数に該当する。これら変数名のみの出力にはnames()関数を使う。今回のデータだと、列の数が少ないこともあり、一画面に全列が表示されるが、数百列のデータとなると画面に収まらないので、変数名を確認しておくことを推奨する。
names(raw_df)[1] "treatment" "gender" "yob" "hh_size" "voted2000" "voted2002" "voted2004" "voted2006"3 データハンドリング
パイプ演算子には{magrittr}パッケージが提供する%>%とR 4.1から提供されるネイティブパイプ演算子の|>がある。現在の主流は古くから使われてきた%>%であるが、今後、|>が主流になると考えられるため、本講義では|>を使用する。しかし、多くの場合、|>の代わりに%>%を使っても同じ結果が得られる。
ヒントパイプ演算子のショートカットキー
Rのnativeパイプ演算子（|>）は意外と打ちにくい。ショートカットキーを活用しよう。
- macOS：Cmd（⌘） + Shift + m
- Windows：Ctrl（Control） + Shift + m
パイプ演算子はパイプ前のオブジェクトを、パイプ後の関数の第一引数として渡す単純な演算子だ。たとえば、列名を変更する関数はrename()であるが、使い方はrenames(データ名, 新しい列名 = 既存の列名, ...)である。raw_dfのgender列の名前をfemaleに変更する場合は以下のように書く。
rename(raw_df, female = gender)# A tibble: 5,000 × 8
treatment female yob hh_size voted2000 voted2002 voted2004 voted2006
<chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr>
1 Control male 1981 3 yes no no no
2 Control female 1957 2 no yes no no
3 Control male 1975 2 no no yes no
4 Hawthorne male 1982 4 no no no no
5 Control female 1949 2 no no yes yes
6 Control female 1937 2 yes yes yes yes
7 Control male 1967 2 no no yes no
8 Control male 1960 2 no no no no
9 Civic Duty male 1952 2 no no yes no
10 Control female 1953 2 no no yes no
# ℹ 4,990 more rows ここで第1引数がraw_dfだが、パイプ演算子を使うと以下のようになり、人間にとって読みやすいコードになる。
raw_df |>
rename(female = gender)# A tibble: 5,000 × 8
treatment female yob hh_size voted2000 voted2002 voted2004 voted2006
<chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr>
1 Control male 1981 3 yes no no no
2 Control female 1957 2 no yes no no
3 Control male 1975 2 no no yes no
4 Hawthorne male 1982 4 no no no no
5 Control female 1949 2 no no yes yes
6 Control female 1937 2 yes yes yes yes
7 Control male 1967 2 no no yes no
8 Control male 1960 2 no no no no
9 Civic Duty male 1952 2 no no yes no
10 Control female 1953 2 no no yes no
# ℹ 4,990 more rows 要するに、X |> Yは「X（の結果）を使ってYを行う」ことを意味する。より詳しいパイプ演算子の解説は『私たちのR』の「データハンドリング [抽出]」を参照されたい。
続いて、変数のリコーディングをしてみよう。xの値が"A"なら1、それ以外は0のように、戻り値が2種類の場合、if_else()関数でリコーディングする。書き方は以下の通りだ。
if_else(条件式, 条件が満たされる場合の戻り値, 条件が満たされない場合の戻り値) たとえば、raw_dfのgender列の値が"female"なら1、それ以外なら0とし、その結果をfemale列として追加するコードは以下の通り。同値を意味する演算子が=でなく、==であることに注意すること（=は<-と同じ代入演算子であるが、Rでは代入演算子として=より<-の使用を推奨している）。
mutate(raw_df,
female = if_else(gender == "female", 1, 0))# A tibble: 5,000 × 9
treatment gender yob hh_size voted2000 voted2002 voted2004 voted2006 female
<chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr> <dbl>
1 Control male 1981 3 yes no no no 0
2 Control female 1957 2 no yes no no 1
3 Control male 1975 2 no no yes no 0
4 Hawthorne male 1982 4 no no no no 0
5 Control female 1949 2 no no yes yes 1
6 Control female 1937 2 yes yes yes yes 1
7 Control male 1967 2 no no yes no 0
8 Control male 1960 2 no no no no 0
9 Civic Duty male 1952 2 no no yes no 0
10 Control female 1953 2 no no yes no 1
# ℹ 4,990 more rows mutate()は指定された列に対して何らかの処理を行い、その結果を新しい列として追加するか、上書きする関数である。このmutate()関数の第1引数もデータであるため、以下のようにパイプ演算子を使うこともできる。
raw_df |>
mutate(female = if_else(gender == "female", 1, 0))# A tibble: 5,000 × 9
treatment gender yob hh_size voted2000 voted2002 voted2004 voted2006 female
<chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr> <dbl>
1 Control male 1981 3 yes no no no 0
2 Control female 1957 2 no yes no no 1
3 Control male 1975 2 no no yes no 0
4 Hawthorne male 1982 4 no no no no 0
5 Control female 1949 2 no no yes yes 1
6 Control female 1937 2 yes yes yes yes 1
7 Control male 1967 2 no no yes no 0
8 Control male 1960 2 no no no no 0
9 Civic Duty male 1952 2 no no yes no 0
10 Control female 1953 2 no no yes no 1
# ℹ 4,990 more rows また、mutate()内には複数のコードを書くこともできる。voted2000列からvoted2006列までそれぞれの値が"yes"であれば、1を、それ以外の場合は0にリコーディングしてみよう。
raw_df |>
mutate(female = if_else(gender == "female", 1, 0),
voted2000 = if_else(voted2000 == "yes", 1, 0),
voted2002 = if_else(voted2002 == "yes", 1, 0),
voted2004 = if_else(voted2004 == "yes", 1, 0),
voted2006 = if_else(voted2006 == "yes", 1, 0))# A tibble: 5,000 × 9
treatment gender yob hh_size voted2000 voted2002 voted2004 voted2006 female
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Control male 1981 3 1 0 0 0 0
2 Control female 1957 2 0 1 0 0 1
3 Control male 1975 2 0 0 1 0 0
4 Hawthorne male 1982 4 0 0 0 0 0
5 Control female 1949 2 0 0 1 1 1
6 Control female 1937 2 1 1 1 1 1
7 Control male 1967 2 0 0 1 0 0
8 Control male 1960 2 0 0 0 0 0
9 Civic Duty male 1952 2 0 0 1 0 0
10 Control female 1953 2 0 0 1 0 1
# ℹ 4,990 more rows また、パイプ演算子は2つ以上使うこともできる。たとえば、rename()を使ってgender列をfemaleに変更し、mutate()でリコーディングを行う場合、以下のように書く。これはraw_dfを使ってrename()の処理を行い、その結果をmutate()関数のデータとして渡すことを意味する。
raw_df |>
rename(female = gender) |>
mutate(female = if_else(female == "female", 1, 0),
voted2000 = if_else(voted2000 == "yes", 1, 0),
voted2002 = if_else(voted2002 == "yes", 1, 0),
voted2004 = if_else(voted2004 == "yes", 1, 0),
voted2006 = if_else(voted2006 == "yes", 1, 0))# A tibble: 5,000 × 8
treatment female yob hh_size voted2000 voted2002 voted2004 voted2006
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Control 0 1981 3 1 0 0 0
2 Control 1 1957 2 0 1 0 0
3 Control 0 1975 2 0 0 1 0
4 Hawthorne 0 1982 4 0 0 0 0
5 Control 1 1949 2 0 0 1 1
6 Control 1 1937 2 1 1 1 1
7 Control 0 1967 2 0 0 1 0
8 Control 0 1960 2 0 0 0 0
9 Civic Duty 0 1952 2 0 0 1 0
10 Control 1 1953 2 0 0 1 0
# ℹ 4,990 more rows 以上のコードはデータを加工し、その結果を出力するだけであって、その結果を保存しない。もう一度raw_dfを出力してみても、これまでのデータ加工内容は反映されていないことが分かる。
raw_df# A tibble: 5,000 × 8
treatment gender yob hh_size voted2000 voted2002 voted2004 voted2006
<chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr>
1 Control male 1981 3 yes no no no
2 Control female 1957 2 no yes no no
3 Control male 1975 2 no no yes no
4 Hawthorne male 1982 4 no no no no
5 Control female 1949 2 no no yes yes
6 Control female 1937 2 yes yes yes yes
7 Control male 1967 2 no no yes no
8 Control male 1960 2 no no no no
9 Civic Duty male 1952 2 no no yes no
10 Control female 1953 2 no no yes no
# ℹ 4,990 more rows このように頑張ってデータを加工したもののその結果が全く反映されていない。加工したデータを引き続き使っていくためには、加工結果を作業環境内に保存する必要がある。作業環境内にオブジェクトを保存するためには代入演算子（<-）を使い、名前を付けて作業空間内に保存する（ファイルとして保存されるわけではない）必要がある。今回は加工の結果をdfという名で保存する。raw_dfに上書きしても問題はないが、生データはとりあえず作業空間内に残しておくことを推奨する（Rに慣れれば上書きしても良い）。
df <- raw_df |>
rename(female = gender) |>
mutate(female = if_else(female == "female", 1, 0),
voted2000 = if_else(voted2000 == "yes", 1, 0),
voted2002 = if_else(voted2002 == "yes", 1, 0),
voted2004 = if_else(voted2004 == "yes", 1, 0),
voted2006 = if_else(voted2006 == "yes", 1, 0))
df# A tibble: 5,000 × 8
treatment female yob hh_size voted2000 voted2002 voted2004 voted2006
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Control 0 1981 3 1 0 0 0
2 Control 1 1957 2 0 1 0 0
3 Control 0 1975 2 0 0 1 0
4 Hawthorne 0 1982 4 0 0 0 0
5 Control 1 1949 2 0 0 1 1
6 Control 1 1937 2 1 1 1 1
7 Control 0 1967 2 0 0 1 0
8 Control 0 1960 2 0 0 0 0
9 Civic Duty 0 1952 2 0 0 1 0
10 Control 1 1953 2 0 0 1 0
# ℹ 4,990 more rows
ノート参考）より効果的な書き方
ちなみに、across()関数とラムダ式（無名関数）を組み合わせると以上のコードをより効率的に書くこともできる。across()は強力な関数だが、初心者にはやや難しいかも知れない。詳細は『私たちのR』の「データハンドリング[要約]」を参照されたい。
df <- raw_df |>
rename(female = gender) |>
mutate(female = if_else(female == "female", 1, 0),
across(starts_with("voted"), \(x) if_else(x == "yes", 1, 0)))4 記述統計量
記述統計量の計算には{summarytools}のdescr()関数が便利だ。descr(データ名)を入力するだけで各変数の記述統計量が出力される。実際にやってみると分かるが、情報量がかなり多い。しかし、実際の論文では各変数の歪度や尖度まで報告することはあまりないだろう。ここではstats引数を追加して、論文などでよく使う平均値（"mean"）、標準偏差（"sd"）、最小値（"min"）、最大値（"max"）、有効ケース数（"n.valid"）のみ出力する。
df |>
descr(stats = c("mean", "sd", "min", "max", "n.valid"))Descriptive Statistics
df
N: 5000
female hh_size voted2000 voted2002 voted2004 voted2006 yob
------------- --------- --------- ----------- ----------- ----------- ----------- ---------
Mean 0.50 2.19 0.25 0.40 0.41 0.33 1956.48
Std.Dev 0.50 0.79 0.43 0.49 0.49 0.47 14.58
Min 0.00 1.00 0.00 0.00 0.00 0.00 1911.00
Max 1.00 7.00 1.00 1.00 1.00 1.00 1986.00
N.Valid 5000.00 5000.00 5000.00 5000.00 5000.00 5000.00 5000.00 ただし、descr()を使うと数値型（numeric）変数の記述統計量のみ表示される。dfだと、treatment列は文字型（character）であるため、表示されない5。各グループがサンプルの何割かを計算するためには、treatment変数をダミー変数へ変換する必要がある。ダミー変数の作成は面倒な作業であるが、{fastDummies}パッケージのdummy_cols()を使えば簡単にできる。dummy_cols()の中にはselect_columns = "ダミー化する列名"を入れれば、当該変数をダミー変数へ変換し、新しい列として追加してくれる。それではtreatment列をダミー化&追加し、その結果をdfに上書きしてみよう。
df <- df |>
dummy_cols(select_columns = "treatment")
df# A tibble: 5,000 × 13
treatment female yob hh_size voted2000 voted2002 voted2004 voted2006 `treatment_Civic Duty`
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 Control 0 1981 3 1 0 0 0 0
2 Control 1 1957 2 0 1 0 0 0
3 Control 0 1975 2 0 0 1 0 0
4 Hawthorne 0 1982 4 0 0 0 0 0
5 Control 1 1949 2 0 0 1 1 0
6 Control 1 1937 2 1 1 1 1 0
7 Control 0 1967 2 0 0 1 0 0
8 Control 0 1960 2 0 0 0 0 0
9 Civic Duty 0 1952 2 0 0 1 0 1
10 Control 1 1953 2 0 0 1 0 0
# ℹ 4,990 more rows
# ℹ 4 more variables: treatment_Control <int>, treatment_Hawthorne <int>,
# treatment_Neighbors <int>, treatment_Self <int> 画面には表示されないが、出力結果の下段を見るとtreatment_で始まるいくつかの変数が追加されたことが分かる。ここでは"tretmant"で始まる列のみを抽出つして確認してみよう。
df |>
select(starts_with("treatment"))# A tibble: 5,000 × 6
treatment `treatment_Civic Duty` treatment_Control treatment_Hawthorne treatment_Neighbors
<chr> <int> <int> <int> <int>
1 Control 0 1 0 0
2 Control 0 1 0 0
3 Control 0 1 0 0
4 Hawthorne 0 0 1 0
5 Control 0 1 0 0
6 Control 0 1 0 0
7 Control 0 1 0 0
8 Control 0 1 0 0
9 Civic Duty 1 0 0 0
10 Control 0 1 0 0
# ℹ 4,990 more rows
# ℹ 1 more variable: treatment_Self <int> select()関数内には抽出する列名を入力するだけで良い。たとえば、femaleとyob列を抽出するならselect(female, yob)である。また、femaleからvoted2006までの意味でfemale:voted2006のような書き方もできる。他にも上の例のようにstarts_with()やends_with()、contain()を使って特定の文字列で始まる（で終わる、を含む）列を指定することもできる。一部の列を除外する場合は変数名の前に!か-を付ける。
とにかく、問題なくダミー化されていることが分かる。もう一度記述統計量を出してみよう。descr()は仕様上、出力される変数の順番はアルファベット順になるが、ここでは元の順番を維持するためにorder = "p"を追加する。また、通常の記述統計表が、先ほど見たものとは違って、各行が変数を、列は記述統計量を表す場合が多い。このように行と列を交換するためにはtranspose = TRUEを追加する6。
df |>
descr(stats = c("mean", "sd", "min", "max", "n.valid"),
order = "p", transpose = TRUE, headings = FALSE)| Mean | Std.Dev | Min | Max | N.Valid | |
|---|---|---|---|---|---|
| female | 0.50 | 0.50 | 0.00 | 1.00 | 5000.00 |
| yob | 1956.48 | 14.58 | 1911.00 | 1986.00 | 5000.00 |
| hh_size | 2.19 | 0.79 | 1.00 | 7.00 | 5000.00 |
| voted2000 | 0.25 | 0.43 | 0.00 | 1.00 | 5000.00 |
| voted2002 | 0.40 | 0.49 | 0.00 | 1.00 | 5000.00 |
| voted2004 | 0.41 | 0.49 | 0.00 | 1.00 | 5000.00 |
| voted2006 | 0.33 | 0.47 | 0.00 | 1.00 | 5000.00 |
| treatment_Civic Duty | 0.12 | 0.32 | 0.00 | 1.00 | 5000.00 |
| treatment_Control | 0.55 | 0.50 | 0.00 | 1.00 | 5000.00 |
| treatment_Hawthorne | 0.11 | 0.31 | 0.00 | 1.00 | 5000.00 |
| treatment_Neighbors | 0.11 | 0.31 | 0.00 | 1.00 | 5000.00 |
| treatment_Self | 0.12 | 0.32 | 0.00 | 1.00 | 5000.00 |
他にも以下のようにdfSummary()関数を使えば、綺麗な表としてまとめてくれる。しかも文字型、factor型変数の場合も度数分布表を作成してくれるので非常に便利だ。これも{summarytools}パッケージに含まれた機能なので、別途、パッケージを読み込む必要はない。
df |>
select(-starts_with("treatment_")) |>
dfSummary(headings = FALSE) |>
print(method = "render", round.digits = 3)| No | Variable | Stats / Values | Freqs (% of Valid) | Graph | Valid | Missing | ||||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | treatment [character] |
|
|
 |
5000 (100.0%) | 0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||
| 2 | female [numeric] |
|
|
 |
5000 (100.0%) | 0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||
| 3 | yob [numeric] |
|
74 distinct values |  |
5000 (100.0%) | 0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||
| 4 | hh_size [numeric] |
|
|
 |
5000 (100.0%) | 0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||
| 5 | voted2000 [numeric] |
|
|
 |
5000 (100.0%) | 0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||
| 6 | voted2002 [numeric] |
|
|
 |
5000 (100.0%) | 0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||
| 7 | voted2004 [numeric] |
|
|
 |
5000 (100.0%) | 0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||
| 8 | voted2006 [numeric] |
|
|
 |
5000 (100.0%) | 0 (0.0%) |
Generated by summarytools 1.1.5 (R version 4.6.1)
2026-08-19
5 バランスチェック
バランスチェックの簡単な方法はグループごとに処置前変数（pre-treatment variables）の平均値を比較することである。無作為割当が成功しているのであれば、処置前に測定された変数の平均値は近似するはずである。ここではグループ（treatment）ごとに性別、誕生年、世帯規模、2000〜2004年の投票参加の平均値を比較してみる。
df |>
group_by(treatment) |>
summarise(female = mean(female, na.rm = TRUE),
yob = mean(yob, na.rm = TRUE),
hh_size = mean(hh_size, na.rm = TRUE),
voted2000 = mean(voted2000, na.rm = TRUE),
voted2002 = mean(voted2002, na.rm = TRUE),
voted2004 = mean(voted2004, na.rm = TRUE))# A tibble: 5 × 7
treatment female yob hh_size voted2000 voted2002 voted2004
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Civic Duty 0.490 1957. 2.19 0.241 0.355 0.395
2 Control 0.505 1956. 2.20 0.251 0.401 0.423
3 Hawthorne 0.472 1957. 2.21 0.241 0.418 0.412
4 Neighbors 0.530 1957. 2.17 0.207 0.386 0.401
5 Self 0.506 1957. 2.18 0.272 0.397 0.398それぞれの変数の平均値は群間でかなり近いため、無作為割当が成功したと考えられる。しかし、変数の単位によって判断が難しいかも知れない。たとえば、2つのグループがあり、年齢の平均値の差は3、世帯規模のそれは2だとする。これを見ると年齢の方がよりバランスが取れていないようにも見えるが、年齢の幅は数十であるに対し、世帯規模はせいぜい5〜6程度であろう。したがって、各変数のばらつきまで考慮した比較が適切であり、その方法の一つが標準化バイアス（=標準化平均差）である。
標準化平均差を計算する便利パッケージ、{coblat}を使ってみよう。第1引数はformulaで処置変数 ~ 共変量1 + 共変量2 + ... + 共変量kのように書く。第2引数はdataでこれらの変数が格納されているデータフレームのオブジェクト名を明記する。続いて、s.d.denom引数だが、多くの場合は、"pooled"で良い（実はs.d.denomの既定値が"pooled"なので、省略しても良い）。
blc_chk <- bal.tab(treatment ~ female + yob + hh_size + voted2000 + voted2002 + voted2004,
data = df, s.d.denom = "pooled")
blc_chkBalance summary across all treatment pairs
Type Max.Diff.Un
female Binary 0.0584
yob Contin. 0.0353
hh_size Contin. 0.0456
voted2000 Binary 0.0650
voted2002 Binary 0.0627
voted2004 Binary 0.0283
Sample sizes
Civic Duty Control Hawthorne Neighbors Self
All 580 2763 536 526 595 ここで注目することはMax.Diff.Un列だ。標準化平均差はペアで計算されるため、今回は各共変量ごとに計10個7の標準化平均差が計算される。一つでも標準化平均差が0から大きく離れているとバランスがとれていないことになるため、10個の標準化平均差はすべて見る必要はなく、最大値だけで十分だろう。ここに出力されているMax.Diff.Unがその標準化平均差の絶対値の最大値である。これが予め決めておいた閾値を超えたらアンバランスと判定し、その基準としてよく使われるのは0.1である。
love.plot()関数を使えば、これらの結果を可視化することもできる。threholds = 0.1を追加しておくと、0.1にガイドラインが追加され、どの共変量がアンバランスかがすぐに確認できるようになる。
love.plot(blc_chk, thresholds = 0.1)
グラフは{ggplot2}で作図されたものだから、{ggplot2}の知識があるとレイヤーを追加して図をカスタマイズすることもできる。たとえば、凡例をなくすためには以下のように書く。実はlove.plot()内部のposition引数で変更可能だが、ここでは{ggplot2}のレイヤーをそのまま足せることを見せておこう。以降のコードではlove.plot()内部で凡例の位置を変更する。
love.plot(blc_chk, thresholds = 0.1) +
theme(legend.position = "none")
以上の図でバランスチェックの可視化は終わりだが、一つ注意すべきものがある。それは変数のタイプによって計算式が異なることだ。通常の連続変数であれば、標準化後にグループ間の差分を計算するが、二値変数は標準化をしない。名目変数もバランスチェックではダミー変数に変換されるため、標準化を施さずに群間比較をすることとなる。誤解を避けるためにも、できれば標準化をしなかった変数には何らかの目印を付けた方が良いだろう。love.plot()内にstars = "raw"を付けると、標準化を施さず、生（raw）の値から群間比較をした変数にはアスタリスク（*）が付くようになる。
love.plot(blc_chk, thresholds = 0.1, stars = "raw", position = "none")
この図を日本語の論文に掲載するためには、中身の言語も日本語にすべきだろう。以下はコードはタイトルや値のラベルを日本語に変換した例である。
love.plot(blc_chk, thresholds = 0.1, stars = "raw", position = "none",
var.names = c("female" = "女性",
"yob" = "誕生年",
"hh_size" = "世帯規模",
"voted2000" = "投票有無 (2000)",
"voted2002" = "投票有無 (2002)",
"voted2004" = "投票有無 (2004)")) +
labs(x = "標準化平均差の絶対値", title = "共変量のバランスチェック", subtitle = "ペアごと比較の最大値",
caption = "注：「*」は標準化を施さなかった二値変数を意味する")
6 処置効果の推定
6.1 グループごとの応答変数の平均値
処置効果を確認するためには各グループごとの応答変数（ここではvoted2006）の平均値を計算し、処置群の平均値から統制群の平均値を引く必要がある。まずは、特定の変数の平均値を計算する方法について紹介する。データ内にある特定の変数の平均値を計算するためにはsummarise()関数内に平均値を求めるmean()関数を入れる。たとえば、dfのvoted2006の平均値を計算するコードは以下の通りである。
df |>
summarise(mean(voted2006, na.rm = TRUE))# A tibble: 1 × 1
`mean(voted2006, na.rm = TRUE)`
<dbl>
1 0.326 na.rm = TRUEは「欠損値があれば、それを除外する」を意味し、指定されていない場合（=既定値）はFALSEになる。今回は欠損値がないものの、念の為に入れておく。
出力結果を見ると、平均値が表示される列の名前が`mean(voted2006, na.rm = TRUE)`となっており、非常に見にくい。この場合、以下のようにmean()の前に出力される列名を予め指定することもできる。
df |>
# voted2006の平均値が表示される列名を Outcome にする。
summarise(Outcome = mean(voted2006, na.rm = TRUE))# A tibble: 1 × 1
Outcome
<dbl>
1 0.326 我々が知りたいのはvoted2006の平均値でなく、グループごとの平均値だろう。被験者がどのグループに属しているかわ示す変数はtreatmentであるが、summarise()にデータを渡す前にgroup_by()変数を使うと、グループごとに計算を行い、その結果を返す。
df |>
group_by(treatment) |>
summarise(Outcome = mean(voted2006, na.rm = TRUE))# A tibble: 5 × 2
treatment Outcome
<chr> <dbl>
1 Civic Duty 0.331
2 Control 0.303
3 Hawthorne 0.360
4 Neighbors 0.405
5 Self 0.328 group_by()内でも=演算子を使うと、グループ名が出力される列名を変更することができる。
df |>
# グループ名が表示される列名を Group にする。
group_by(Groups = treatment) |>
summarise(Outcome = mean(voted2006, na.rm = TRUE))# A tibble: 5 × 2
Groups Outcome
<chr> <dbl>
1 Civic Duty 0.331
2 Control 0.303
3 Hawthorne 0.360
4 Neighbors 0.405
5 Self 0.328
ヒントTip）
group_by()を使わないやり方
group_by()を使わず、以下のようにsummarise()内に.by引数を使っても良い。
df |>
summarise(Outcome = mean(voted2006, na.rm = TRUE), .by = treatment) ここで一つ注目したいのが、グループの表示順番である。変数のデータ型が文字型だと（Rコンソール上でclass(df$treatment)を入力するか、dfの出力画面でtreatmentの下に<chr>と表示されていることで確認できる）、今のようにアルファベット順で表示される。しかし、統制群は最初か最後に来るのが通例である。この順番をアルファベット順でなく、任意の順番にするためにはtreatment変数をfactor型変数へ変換する必要がある。Factor型は「順序付きの文字型変数」だと理解しても良い8。列の追加・上書き（今回はtreatment列の上書き）の処理が必要なのでmutate()関数を使う。変数をfactor型に変換する関数はfactor()関数で、第1引数としてはfactor型へ変換する変数名を指定する。第2引数はlevelsであり、出力したい順番の文字型ベクトルを指定する。スペルミスには注意すること。
df |>
mutate(treatment = factor(treatment,
levels = c("Control", "Civic Duty",
"Self", "Neighbors", "Hawthorne")))# A tibble: 5,000 × 13
treatment female yob hh_size voted2000 voted2002 voted2004 voted2006 `treatment_Civic Duty`
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 Control 0 1981 3 1 0 0 0 0
2 Control 1 1957 2 0 1 0 0 0
3 Control 0 1975 2 0 0 1 0 0
4 Hawthorne 0 1982 4 0 0 0 0 0
5 Control 1 1949 2 0 0 1 1 0
6 Control 1 1937 2 1 1 1 1 0
7 Control 0 1967 2 0 0 1 0 0
8 Control 0 1960 2 0 0 0 0 0
9 Civic Duty 0 1952 2 0 0 1 0 1
10 Control 1 1953 2 0 0 1 0 0
# ℹ 4,990 more rows
# ℹ 4 more variables: treatment_Control <int>, treatment_Hawthorne <int>,
# treatment_Neighbors <int>, treatment_Self <int> treatment列名の下が<fct>となっていることが分かる。これはtreatment列のデータ型がfactor型であることを意味する。問題なく動くことが確認できたので、dfを上書きしよう。
df <- df |>
mutate(treatment = factor(treatment,
levels = c("Control", "Civic Duty", "Hawthorne",
"Self", "Neighbors")))
ヒントTip）より簡単なfactor化方法
fct_relevel()を使えば、より簡単にfactor型に変換できる。factor化する変数名を最初に書き、続いて要素の順番を指定するだけだ。fct_relevel()を使うには予め{forcats}パッケージを読み込んでおく必要があるが、{tidyverse}を読み込むと{forcats}も同時に読み込まれる。
df <- df |>
mutate(treatment = fct_relevel(treatment, "Control", "Civic Duty",
"Self", "Neighbors", "Hawthorne")) それでは、改めてグループごとのvoted2006の平均値を計算してみよう。今回は計算結果をout_mean_dfという名のオブジェクトとして格納する。
out_mean_df <- df |>
group_by(Groups = treatment) |>
summarise(Outcome = mean(voted2006, na.rm = TRUE))
out_mean_df# A tibble: 5 × 2
Groups Outcome
<fct> <dbl>
1 Control 0.303
2 Civic Duty 0.331
3 Hawthorne 0.360
4 Self 0.328
5 Neighbors 0.405今回は統制群は最初に出力されていることが確認できる。
それではこの結果をグラフとして示してみよう。作図には{ggplot2}パッケージを使う。まずはout_mean_dfをggplot()関数に渡す。ggplot()関数以降は、+演算子を使ってレイヤーを足していくこととなる。棒グラフのレイヤーはgeom_col()関数であり、その中にaes()関数を入れる。aes()の中には棒グラフの作図に必要な情報を入れる必要がある（これをマッピング（mapping）と呼ぶ）。棒グラフを作成するために必要な最低限の情報とは各棒の横軸上の位置（x）と棒の高さ（y）だ。今回は横軸がグループ名、縦軸が平均値となる棒グラフを作る。`
out_mean_df |>
ggplot() +
geom_col(aes(x = Groups, y = Outcome))
続いて、このグラフの見た目を調整してみよう。
out_mean_df |>
ggplot() +
geom_col(aes(x = Groups, y = Outcome)) +
# 縦軸（y軸）のラベルを変更する
labs(y = "Mean(Outcome)") +
# grayテーマ（デフォルトのテーマ）を使用し、フォントサイズは14
theme_gray(base_size = 12)
また、geom_label()レイヤーを足すと、棒の上にラベルを付けることもできる。ラベルに必要な情報は各ラベルの横軸上の位置（x）、縦軸上の位置（y）、ラベルの表示内容（label）だ。今回のラベルは平均値の具体的な数値を入れてみよう。
out_mean_df |>
ggplot() +
geom_col(aes(x = Groups, y = Outcome)) +
geom_label(aes(x = Groups, y = Outcome, label = Outcome)) +
labs(y = "Mean(Outcome)") +
theme_gray(base_size = 12)
小数点が長すぎるので3桁まで表示としよう。ここではsprintf()を使用する。使い方が簡単とは言えないが、覚える必要はなく、必要な時にググるか、本資料のコードをコピペすれば良い9。
out_mean_df |>
ggplot() +
geom_col(aes(x = Groups, y = Outcome)) +
# 2桁までなら %.3f を %.2f に変更
geom_label(aes(x = Groups, y = Outcome, label = sprintf("%.3f", Outcome))) +
labs(y = "Mean(Outcome)") +
theme_gray(base_size = 12)
これで可視化ができた。ただし、以上のコードには改善の余地がある。geom_bar()とgeom_label()内のaes()関数に注目して欲しい。よく見るとxとyと同じだろう。geom_*()が共有するマッピングがあれば、ggplot()内で指定することでコードを効率化することもできる。
out_mean_df |>
ggplot(aes(x = Groups, y = Outcome)) +
geom_col() +
geom_label(aes(label = sprintf("%.3f", Outcome))) +
labs(y = "Mean(Outcome)") +
theme_gray(base_size = 12)
ggplot()内でまとめる
6.2 統計的推定（単回帰分析）
これまでの作業はグループごとの応答変数の平均値であって、処置効果ではない。処置効果を計算するためには処置群の平均値から統制群の平均値を引く必要がある。たとえば、Civic Dutyはがき群の平均値は約0.331、統制群のそれは0.303であるため、Civic Dutyはがきの処置効果は約0.028である。しかし、これを各グループごとに計算することは面倒だし、何よりも得られた値が点推定値だという限界がある。得られた処置効果の不確実性は計算できない。
ここで有効なのが線形回帰分析である。回帰分析を行うことで処置効果の点推定値のみならず、不確実性の指標である標準誤差も計算され、区間推定や統計的仮説検定も可能となる。線形回帰分析の関数はlm()だ。第1引数としては回帰式であり、応答変数 ~ 説明変数と表記する。第2引数はdataであり、回帰式で指定した変数が入っているデータ名を指定する。回帰分析の結果は名前を付けてオブジェクトとして格納し、summary()関数を使うと、詳細が確認できる。
fit1 <- lm(voted2006 ~ treatment, data = df)
summary(fit1)
Call:
lm(formula = voted2006 ~ treatment, data = df)
Residuals:
Min 1Q Median 3Q Max
-0.4049 -0.3277 -0.3026 0.6399 0.6974
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.302570 0.008899 34.002 < 2e-16 ***
treatmentCivic Duty 0.028465 0.021364 1.332 0.18279
treatmentHawthorne 0.057505 0.022076 2.605 0.00922 **
treatmentSelf 0.025161 0.021140 1.190 0.23401
treatmentNeighbors 0.102373 0.022251 4.601 4.31e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4677 on 4995 degrees of freedom
Multiple R-squared: 0.004947, Adjusted R-squared: 0.00415
F-statistic: 6.209 on 4 and 4995 DF, p-value: 5.575e-05 ちなみに、これもパイプ演算子を使うことができる。ただし、第1引数として渡すパイプ演算子の特徴上、そのまま使うことはできない。なぜならlm()関数の第1引数はデータでなく、formulだからだ。この場合はプレースホルダー（place holder）を指定する必要がある。パイプ前のオブジェクトが入る位置を任意に指定することであり、_を使う10。
fit1 <- df |>
lm(voted2006 ~ treatment, data = _) Factor型、または文字型変数が説明変数の場合、自動的にダミー変数として処理され、factor型の場合、最初の水準（ここでは"Control"）がベースカテゴリとなる。もし説明変数が文字型なら、アルファベット順で最初の水準がベースカテゴリとなり、今回の例だと"Civic Duty"がベースカテゴリとなる。処置効果は「統制群に比べて〜」が重要となるので、数値型以外の説明変数は予めfactor化しておいた方が望ましい。
Civic Dutyの推定値は約0.028であり、これは統制群（Control）に比べ、Civic Duty群のvoted2006の平均値は約0.028高いことを意味する。応答変数が0、1であるため、これを割合（= 投票率）で換算すると、約2.8pp11高いことを意味する。つまり、Civic Dutyのはがきをもらった被験者はそうでない被験者に比べて投票率が約2.8pp高いことを意味する。他の推定値も同じやり方で解釈すれば良い。
続いて、これらの処置効果が統計的に有意なものかを確認してみよう。統計的有意か否かを判定するためには有意と非有意の境界線が必要である、これは通常、有意水準（significance level; \(\alpha\)）と呼ばれる。この有意水準は分析者が決めるものではあるが、社会科学で広く使われる基準は\(\alpha\) = 0.05、つまり5%だ。分析結果の画面にはPr(>|t|)列が表示されているが、これが p 値と呼ばれるもので、これが0.05以下る場合、処置効果は統計的に有意だと判定する。もし、\(\alpha\) = 0.1を採用するなら、 p \(\leq\) 0.1の場合において統計的に有意と判定する。Civic Dutyの p 値は約0.183であり、0.05より大きいことから統計的に有意な処置効果と判定できない。一方、Neighborsはがきの処置効果の p 値は0.234で、0.05以下であることからNeighborsが投票率に与える影響は統計的に有意であると判定できる。
ノート「5.85e-12」のように表示されたら…
Rでは0に非常に近い値を表示する際、「5.85e-12」のような表記を採用する。たとえば5.85e-12は、5.85 \(\times\) 10-12を意味する。10-1は0.1、10-2は0.01であることを考えると、0に限りなく近い値である。こういう表記は特に p 値においてよく登場し、 p 値が一定値以下であれば< 2e-16と表示される。このように非常に小さい p 値を論文等に掲載する際は、< 0.001と表記するケースが多い。
続いて、この結果を可視化してみよう。ここでも{ggplot2}パッケージを使って可視化をするが、{ggplot2}で使用可能なオブジェクトは表形式のデータである。Rコンソール上でclass(オブジェクト名)を入力すると、データのクラスが出力されるが、このクラスに"data.frame"があれば、{ggplot2}で使用できる。たとえば、fit1オブジェクトのクラスは"lm"であるため、そのまま{ggplot2}で使うことはできない。
class(fit1)[1] "lm" 推定結果を表形式に変換するためには{broom}パッケージのtidy()関数が便利だ。使い方は簡単でtidy()内に回帰分析の推定結果が格納されたオブジェクトを入れるだけである。ただし、デフォルトの設定では95%信頼区間が表示されないため、中にはconf.int = TRUEを追加しておく必要がある。
# 90%信頼区間を使うのであれば conf.int = 0.9 を追加（デフォルトは0.95）
fit1_coef <- tidy(fit1, conf.int = TRUE)
fit1_coef# A tibble: 5 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.303 0.00890 34.0 3.90e-228 0.285 0.320
2 treatmentCivic Duty 0.0285 0.0214 1.33 1.83e- 1 -0.0134 0.0703
3 treatmentHawthorne 0.0575 0.0221 2.60 9.22e- 3 0.0142 0.101
4 treatmentSelf 0.0252 0.0211 1.19 2.34e- 1 -0.0163 0.0666
5 treatmentNeighbors 0.102 0.0223 4.60 4.31e- 6 0.0588 0.146 class(fit1_coef)[1] "tbl_df" "tbl" "data.frame" fit1_coefのクラスに"data.frame"が含まれているので、これを使って作図することができる。
作図する前に、fit1_coefの加工しておきたい。それぞれの係数（estimate列）は処置効果を表しているが、切片（"(Intercept)"）の推定値は処置効果とは無関係である。したがって、予め切片の行を除外しておきたい。特定の行を残したり、除外する関数はfilter()である。今回はterm列の値が"(Intercept)"ではない行を残したいので、同値演算子（==）の否定を意味する!=演算子を使用する。
fit1_coef <- fit1_coef |>
filter(term != "(Intercept)")
fit1_coef# A tibble: 4 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 treatmentCivic Duty 0.0285 0.0214 1.33 0.183 -0.0134 0.0703
2 treatmentHawthorne 0.0575 0.0221 2.60 0.00922 0.0142 0.101
3 treatmentSelf 0.0252 0.0211 1.19 0.234 -0.0163 0.0666
4 treatmentNeighbors 0.102 0.0223 4.60 0.00000431 0.0588 0.146 それでは作図に入ろう。処置効果を示す場合は、点推定値以外にもその不確実性を示すのは一般的である。不確実性の指標として幅広く使われるのは標準誤差（standard error; 標準偏差ではない）であるが、可視化の際にはこの標準誤差に基づき計算した信頼区間を示すのが一般的だ。有意水準が5%（\(\alpha\) = 0.05）であれば、95%信頼区間を示し、10%（\(\alpha\) = 0.1）なら90%信頼区間を用いる。tidy()で得られたデータの場合、信頼区間の下限と上限はそれぞれconf.lowとconf.highという名の列に格納されている（conf.int = TRUEを指定しておかないと、信頼区間は計算されない）。
点と区間を同時に示すプロットがpoint-rangeプロットであり、{ggplot2}ではgeom_pointrange()レイヤーを使う。必要な情報はpoint-rangeの横軸上の位置（x）、点の縦軸上の位置（y）、区間の上限（ymax）と下限（ymin）である。これらの情報は全てfit1_coefに入っているため、fit1_coefをそのままggplot()関数に渡して作図することができる。
fit1_coef |>
ggplot() +
geom_pointrange(aes(x = term, y = estimate,
ymin = conf.low, ymax = conf.high))
それでは図をカスタマイズしてみよう。図内の様々なラベルを修正するlabs()レイヤーでラベルを修正する。テーマはデフォルトのtheme_gray()の代わりに白黒テーマ（theme_bw()）を使用し、フォントサイズは12とする。また、y = 0の水平線を追加する。95%信頼区間内に0が含まれる場合、「5%水準で統計的に有意でない」と判断できる。水平線を描くにはgeom_hline()レイヤーを追加し、yintercept = 0を指定することで、0のところに水平線が表示される。この水平線は通常、実線であるがlinetype = "dashed"を追加することで破線に変更することができる。
fit1_coef |>
ggplot() +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_pointrange(aes(x = term, y = estimate,
ymin = conf.low, ymax = conf.high)) +
labs(x = "Treatments", y = "Average Treatment Effects") +
theme_bw(base_size = 12)
まだ気になる点がある。それは横軸の目盛りラベルにtreatmentという不要な情報がある点だ。これは作図の時点で修正することも可能だが、まずはdfのterm変数の値を修正する方法を紹介する。変数の値を修正する時にはrecode()関数を使用する。第1引数はリコーディングする変数名であり、引き続き"元の値" = "新しい値"を指定すれば良い。スペルミスに注意すること。
fit1_coef <- fit1_coef |>
mutate(term = recode(term,
"treatmentCivic Duty" = "Civic Duty",
"treatmentHawthorne" = "Hawthorne",
"treatmentSelf" = "Self",
"treatmentNeighbors" = "Neighbors"))
fit1_coef# A tibble: 4 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Civic Duty 0.0285 0.0214 1.33 0.183 -0.0134 0.0703
2 Hawthorne 0.0575 0.0221 2.60 0.00922 0.0142 0.101
3 Self 0.0252 0.0211 1.19 0.234 -0.0163 0.0666
4 Neighbors 0.102 0.0223 4.60 0.00000431 0.0588 0.146 以上の作業はterm列の各値から"treatment"文字を""に置換することなので、文字列を置換する関数であるstr_replace()を使えば、より短くすることができる[^str-remove]。
ヒントTip）文字列から特定の文字列を削除する方法
str_remove(変数名, "削除する文字")で指定した変数から特定の文字列を削除することができる。今回の例でいうと、term列の文字列から”treatment”を削除すれば、recode()は不要であろう。
fit1_coef <- fit1_coef |>
mutate(term = str_remove(term, "treatment")) 類似した関数としてstr_replace(変数名, "文字列1", "文字列2")がある。これは指定した変数の文字列1を文字列2に置換する関数だ。str_remove()は”treatment”を”“に置換したものだから、実はstr_replace()はstr_remove()の上位互換とも言えよう。
fit1_coef <- fit1_coef |>
mutate(term = str_replace(term, "treatment", "")) fit1_coefも修正できたので、 図 11 と同じコードでもう一度作図してみよう。
fit1_coef |>
ggplot() +
geom_hline(yintercept = 0) +
geom_pointrange(aes(x = term, y = estimate,
ymin = conf.low, ymax = conf.high)) +
labs(x = "Treatments", y = "Average Treatment Effects") +
theme_bw(base_size = 12)
最後に横軸の順番を修正してみよう。fit1_coefのterm列は文字型変数であるため、アルファベット順になる。これをdfのtreatment列と同様、Civic Duty、Self、Neighbors、Hawthorneの順にしたい。この場合fit1_coefのterm列をfactor化すれば良い。factor()関数を使っても良いが、ここではまた便利な技を紹介しよう。それはfct_inorder()関数だ。これは表示されている順番をfactorの順番とする関数だ。実際、fit1_coefの中身を見ると、表示順番はCivic Duty、Self、Neighbors、Hawthorneだ。非常に嬉しい状況なので、fct_inorder()を使ってみよう。
fit1_coef <- fit1_coef |>
mutate(term = fct_inorder(term))
fit1_coef# A tibble: 4 × 7
term estimate std.error statistic p.value conf.low conf.high
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Civic Duty 0.0285 0.0214 1.33 0.183 -0.0134 0.0703
2 Hawthorne 0.0575 0.0221 2.60 0.00922 0.0142 0.101
3 Self 0.0252 0.0211 1.19 0.234 -0.0163 0.0666
4 Neighbors 0.102 0.0223 4.60 0.00000431 0.0588 0.146 それでは、 図 12 と同じコードでもう一度作図してみよう。
fit1_coef |>
ggplot() +
geom_hline(yintercept = 0) +
geom_pointrange(aes(x = term, y = estimate,
ymin = conf.low, ymax = conf.high)) +
labs(x = "Treatments", y = "Average Treatment Effects") +
theme_bw(base_size = 12)
これで処置効果の可視化もバッチリだ。
6.3 多重比較の問題
グループが2つ、つまり統制群と統制群のみが存在する場合、我々が比較を行う回数は1回のみである（統制群 - 処置群）。しかし、今回のデータの場合、処置群は4つである。これは比較を4回行うことを意味する。具体的には「統制群 - 処置群1」、「統制群 - 処置群2」、「統制群 - 処置群3」、「統制群 - 処置群4」だ。比較を繰り返すほど、統計的に有意な結果が得られる可能性は高い。極端な話、1000回程度検定を繰り返せば、本当は効果がなくてもたまたま統計的に有意な結果が何回かは得られるだろう。これが多重比較の問題（または、多重検定の問題）である。したがって、比較の回数が多くなるにつれ、統計的有意性検定にも何らかのペナルティーを課す必要がある。詳細は「ランダム化比較試験」の講義資料を参照されたい。
多重比較におけるペナルティーの付け方はいくつかあるが、ここでは最も保守的な（=研究者にとって都合の悪い）補正法であるボンフェローニ補正（Bonferroni correction）を紹介する。これは非常に単純で、 p 値や信頼区間を計算する際、「統計的有意」と判定されるハードルを上げる方法である。予め決めておいた有意水準（\(\alpha\)）が0.05で、比較の回数が4回であれば、 p 値が0.05 \(\times \frac{\sf 1}{\sf 4}\) = 0.0125を下回る場合において「5%水準で有意である」と判定する。信頼区間でいえば通常の95%信頼区間（1 - 0.05）でなく、98.75%信頼区間（1 - 0.0125）を使うこととなる。補正を施してもなお統計的に有意な結果が得られたら「1.25%水準で〜」と解釈するのではなく、「5%水準で〜」と解釈する必要がある。
95%以外の信頼区間を求めるのは簡単で、tidy()関数内にconf.levelを修正すれば良い。指定されていない場合はデフォルトで0.95が割り当てられているが、これを0.9875と修正する。
fit1_coef <- tidy(fit1, conf.int = TRUE, conf.level = 0.9875)
fit1_coef# A tibble: 5 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 0.303 0.00890 34.0 3.90e-228 0.280 0.325
2 treatmentCivic Duty 0.0285 0.0214 1.33 1.83e- 1 -0.0249 0.0818
3 treatmentHawthorne 0.0575 0.0221 2.60 9.22e- 3 0.00234 0.113
4 treatmentSelf 0.0252 0.0211 1.19 2.34e- 1 -0.0277 0.0780
5 treatmentNeighbors 0.102 0.0223 4.60 4.31e- 6 0.0468 0.158 それでは 図 13 と同じ図を作ってみよう。まず、切片の行を除外するが、ここではfilter()を使わず、slice()の使った方法を紹介する。slice()は()内に指定した行を残す関数だ。たとえば、slice(fit1_coef, 2)ならfit1_coefの2行目のみを残す。fit1_coefはslice()の第1引数だから、パイプ演算子を使うことも可能で、こちらの方を推奨する。そうすれば()内には残す行のみの指定で済む。slice(2)のみなら2行目を残し、slice(1, 3, 5)なら1、3、5行目を残す。:を使うと「〜行目から〜行目まで」の指定ができる。処置効果の係数はfit1_coefの2行目から5行目までなので、2:5と指定すれば良い。
fit1_coef <- fit1_coef |>
slice(2:5)
fit1_coef# A tibble: 4 × 7
term estimate std.error statistic p.value conf.low conf.high
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 treatmentCivic Duty 0.0285 0.0214 1.33 0.183 -0.0249 0.0818
2 treatmentHawthorne 0.0575 0.0221 2.60 0.00922 0.00234 0.113
3 treatmentSelf 0.0252 0.0211 1.19 0.234 -0.0277 0.0780
4 treatmentNeighbors 0.102 0.0223 4.60 0.00000431 0.0468 0.158 続いて、term変数の値から"treatment"の文字を除去し、fit1_coefでの出力順番でtermをfactor化する。
fit1_coef <- fit1_coef |>
mutate(term = recode(term,
"treatmentCivic Duty" = "Civic Duty",
"treatmentHawthorne" = "Hawthorne",
"treatmentSelf" = "Self",
"treatmentNeighbors" = "Neighbors"),
term = fct_inorder(term))
fit1_coef# A tibble: 4 × 7
term estimate std.error statistic p.value conf.low conf.high
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Civic Duty 0.0285 0.0214 1.33 0.183 -0.0249 0.0818
2 Hawthorne 0.0575 0.0221 2.60 0.00922 0.00234 0.113
3 Self 0.0252 0.0211 1.19 0.234 -0.0277 0.0780
4 Neighbors 0.102 0.0223 4.60 0.00000431 0.0468 0.158 最後に 図 13 と同じコードで作図する。縦軸のタイトルがかなり長くなったので\nで改行する。
fit1_coef |>
ggplot() +
geom_hline(yintercept = 0) +
geom_pointrange(aes(x = term, y = estimate,
ymin = conf.low, ymax = conf.high)) +
labs(x = "Treatments",
y = "Average Treatment Effects\n(w/ 98.75% CI)") +
# シンプルなclassicテーマを使用し、文字サイズを16に
theme_classic(base_size = 16)
全体的に信頼区間が広くなり、統計的有意な結果は得られにくくなったもののHawhtrone、Neighborsの処置効果はいずれも統計的に有意な結果が得られた。
ヒント{multcomp}パッケージを使った補正
pacman::p_load(tidyverse)
multicomp <- glht(fit1, linfct = mcp(treatment = "Tukey"))
summary(multicomp)
summary(multicomp, test = adjusted(type = "bonferroni"))
summary(multicomp, test = adjusted(type = "holm"))
glht(fit1, linfct = mcp(treatment = "Dunnett")) |>
summary() |> modelsummary()6.4 統計的推定（重回帰分析）
今回の例は無作為割当が失敗したという直接的な証拠もなく、共変量（処置前変数）の大きな偏りは見られなかった。しかし、何らかの理由で処置前変数の偏りが生じる場合がある。その何らかの理由でバランスが崩れた共変量が応答変数にまで影響を与えるのであれば、それは交絡変数（confounder）となり、バイアスの原因となる。たとえば、今回の例でいうと、なんらかの理由で年齢のバランスが崩れたとしよう。年齢と投票率の間には強い関連があると言われており、今まで通りに処置群のみの回帰分析を行うとバイアスが発生してしまう。この場合、偏りが生じている処置前変数を統制（control）することによってバイアスを小さくすることができる。今回はバランスチェックで問題は見られなかったから不要であるが、性別や誕生年などの共変量を統制した推定をしてみよう。
やり方は簡単で、lm()内の回帰式を応答変数 ~ 説明変数1 + 説明変数2 + ...のように説明変数を+で足していけば良い。
fit2 <-lm(voted2006 ~ treatment + female + yob + hh_size +
voted2000 + voted2002 + voted2004, data = df)
summary(fit2)
Call:
lm(formula = voted2006 ~ treatment + female + yob + hh_size +
voted2000 + voted2002 + voted2004, data = df)
Residuals:
Min 1Q Median 3Q Max
-0.7524 -0.3422 -0.1955 0.5219 0.9552
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.3762557 0.9028102 7.063 1.86e-12 ***
treatmentCivic Duty 0.0420120 0.0205106 2.048 0.04058 *
treatmentHawthorne 0.0583915 0.0211860 2.756 0.00587 **
treatmentSelf 0.0293015 0.0202883 1.444 0.14873
treatmentNeighbors 0.1140621 0.0213632 5.339 9.75e-08 ***
female -0.0204459 0.0127110 -1.609 0.10778
yob -0.0031784 0.0004634 -6.858 7.81e-12 ***
hh_size 0.0004358 0.0084609 0.052 0.95892
voted2000 0.1014383 0.0149395 6.790 1.25e-11 ***
voted2002 0.1379055 0.0132289 10.425 < 2e-16 ***
voted2004 0.1720424 0.0129749 13.260 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4487 on 4989 degrees of freedom
Multiple R-squared: 0.08528, Adjusted R-squared: 0.08345
F-statistic: 46.51 on 10 and 4989 DF, p-value: < 2.2e-16 {modelsummary}パッケージのmodelsummary()関数を使えば、推定結果がより見やすくなる。
modelsummary(fit2)| (1) | |
|---|---|
| (Intercept) | 6.376 |
| (0.903) | |
| treatmentCivic Duty | 0.042 |
| (0.021) | |
| treatmentHawthorne | 0.058 |
| (0.021) | |
| treatmentSelf | 0.029 |
| (0.020) | |
| treatmentNeighbors | 0.114 |
| (0.021) | |
| female | -0.020 |
| (0.013) | |
| yob | -0.003 |
| (0.000) | |
| hh_size | 0.000 |
| (0.008) | |
| voted2000 | 0.101 |
| (0.015) | |
| voted2002 | 0.138 |
| (0.013) | |
| voted2004 | 0.172 |
| (0.013) | |
| Num.Obs. | 5000 |
| R2 | 0.085 |
| R2 Adj. | 0.083 |
| AIC | 6189.2 |
| BIC | 6267.4 |
| Log.Lik. | -3082.598 |
| F | 46.513 |
| RMSE | 0.45 |
また、複数のモデルをlist()関数でまとめると、モデル間比較もできる。
modelsummary(list("w/o Covariates" = fit1, "w/ Covariates" = fit2))| w/o Covariates | w/ Covariates | |
|---|---|---|
| (Intercept) | 0.303 | 6.376 |
| (0.009) | (0.903) | |
| treatmentCivic Duty | 0.028 | 0.042 |
| (0.021) | (0.021) | |
| treatmentHawthorne | 0.058 | 0.058 |
| (0.022) | (0.021) | |
| treatmentSelf | 0.025 | 0.029 |
| (0.021) | (0.020) | |
| treatmentNeighbors | 0.102 | 0.114 |
| (0.022) | (0.021) | |
| female | -0.020 | |
| (0.013) | ||
| yob | -0.003 | |
| (0.000) | ||
| hh_size | 0.000 | |
| (0.008) | ||
| voted2000 | 0.101 | |
| (0.015) | ||
| voted2002 | 0.138 | |
| (0.013) | ||
| voted2004 | 0.172 | |
| (0.013) | ||
| Num.Obs. | 5000 | 5000 |
| R2 | 0.005 | 0.085 |
| R2 Adj. | 0.004 | 0.083 |
| AIC | 6598.1 | 6189.2 |
| BIC | 6637.2 | 6267.4 |
| Log.Lik. | -3293.044 | -3082.598 |
| F | 6.209 | 46.513 |
| RMSE | 0.47 | 0.45 |
modelsummary()は推定値と標準誤差（カッコ内）が別々の行として出力する。これを一行でまとめるためには、以下のようにコードを修正する。
# listをmodelsummary()の外側へ移し、パイプで繋いでみた
list("w/o Covariates" = fit1,
"w/ Covariates" = fit2) |>
modelsummary(estimate = "{estimate} ({std.error})",
statistic = NULL)| w/o Covariates | w/ Covariates | |
|---|---|---|
| (Intercept) | 0.303 (0.009) | 6.376 (0.903) |
| treatmentCivic Duty | 0.028 (0.021) | 0.042 (0.021) |
| treatmentHawthorne | 0.058 (0.022) | 0.058 (0.021) |
| treatmentSelf | 0.025 (0.021) | 0.029 (0.020) |
| treatmentNeighbors | 0.102 (0.022) | 0.114 (0.021) |
| female | -0.020 (0.013) | |
| yob | -0.003 (0.000) | |
| hh_size | 0.000 (0.008) | |
| voted2000 | 0.101 (0.015) | |
| voted2002 | 0.138 (0.013) | |
| voted2004 | 0.172 (0.013) | |
| Num.Obs. | 5000 | 5000 |
| R2 | 0.005 | 0.085 |
| R2 Adj. | 0.004 | 0.083 |
| AIC | 6598.1 | 6189.2 |
| BIC | 6637.2 | 6267.4 |
| Log.Lik. | -3293.044 | -3082.598 |
| F | 6.209 | 46.513 |
| RMSE | 0.47 | 0.45 |
また、alignで各列を左寄せや右寄せに（文字列は左寄せ、数値は右寄せが一般的）、coef_rename引数で表示される変数名を変更することもできる。
modelsummary(list("w/o Covariates" = fit1, "w/ Covariates" = fit2),
estimate = "{estimate} ({std.error})",
statistic = NULL,
align = "lrr", # 1列は左寄せ、2列は右寄せ、3列は右寄せ
coef_rename = c("treatmentCivic Duty" = "Civic Duty",
"treatmentHawthorne" = "Hawthorne",
"treatmentSelf" = "Self",
"treatmentNeighbors" = "Neighbors",
"female" = "Female",
"yob" = "Year of Birth",
"hh_size" = "Household Size",
"voted2000" = "Voted (2000)",
"voted2002" = "Voted (2002)",
"voted2004" = "Voted (2004)"))| w/o Covariates | w/ Covariates | |
|---|---|---|
| (Intercept) | 0.303 (0.009) | 6.376 (0.903) |
| Civic Duty | 0.028 (0.021) | 0.042 (0.021) |
| Hawthorne | 0.058 (0.022) | 0.058 (0.021) |
| Self | 0.025 (0.021) | 0.029 (0.020) |
| Neighbors | 0.102 (0.022) | 0.114 (0.021) |
| Female | -0.020 (0.013) | |
| Year of Birth | -0.003 (0.000) | |
| Household Size | 0.000 (0.008) | |
| Voted (2000) | 0.101 (0.015) | |
| Voted (2002) | 0.138 (0.013) | |
| Voted (2004) | 0.172 (0.013) | |
| Num.Obs. | 5000 | 5000 |
| R2 | 0.005 | 0.085 |
| R2 Adj. | 0.004 | 0.083 |
| AIC | 6598.1 | 6189.2 |
| BIC | 6637.2 | 6267.4 |
| Log.Lik. | -3293.044 | -3082.598 |
| F | 6.209 | 46.513 |
| RMSE | 0.47 | 0.45 |
処置効果に注目すると、Civic Dutyは処置効果の点推定値が2.8ppから4.2ppへ1.5倍になり、（多重比較の補正をしなければ）統計的に有意な結果が得られた。ランダム化が完璧に成功すると、つまり標準化平均差がすべて0であれば、これら2つのモデルの点推定値は変わらないはずであるが、 図 3 を見ると、今回は若干のズレはある。この場合は、今回のように2つのモデルの推定結果は一致しない。どのモデルを採用するかは分析する側が決める問題であるが、宋の個人的な見解としては、バランスが完全に取れていないのであれば、やはり重回帰分析で共変量のしっかり統制したモデルを採用すべきだと考えている。しかし、今回の例だと、多重比較の問題まで考慮すれば、統計的に有意な処置効果が確認できた処置はHawthorneとNeighborsのみという点で変化はなかったとも言えよう（この後で紹介する可視化で確認できる）。
6.5 複数モデルの可視化
modelsummary()を使えば、複数のモデルの推定結果を一つの表としてまとめられる。しかし、図の場合はどうだろう。共変量なしモデルとありモデルを 図 14 のように一つにまとめることはできるだろうか。もちろん出来る。
まず、重回帰分析を行った結果（fit2）から処置効果の推定値情報を抽出し、fit1_coefと同じ構造のデータとしてまとめる。多重比較の問題に対応するために、ここでもボンフェローニ補正を施す。
fit2_coef <- tidy(fit2, conf.int = TRUE, conf.level = 0.9875)
fit2_coef <- fit2_coef |>
slice(2:5) |>
mutate(term = recode(term,
"treatmentCivic Duty" = "Civic Duty",
"treatmentHawthorne" = "Hawthorne",
"treatmentSelf" = "Self",
"treatmentNeighbors" = "Neighbors"),
term = fct_inorder(term))
fit2_coef# A tibble: 4 × 7
term estimate std.error statistic p.value conf.low conf.high
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Civic Duty 0.0420 0.0205 2.05 0.0406 -0.00924 0.0933
2 Hawthorne 0.0584 0.0212 2.76 0.00587 0.00546 0.111
3 Self 0.0293 0.0203 1.44 0.149 -0.0214 0.0800
4 Neighbors 0.114 0.0214 5.34 0.0000000975 0.0607 0.167 処置効果の推定値や標準誤差などが異なるが、構造としては同じである。続いて、bind_rows()を用い、この2つのデータを一つの表として結合する。2つの表はlist()関数でまとめるが、それぞれ"モデル名" = データ名と指定する。最後に、.id = "Model"を追加する。
bind_rows(list("Model 1" = fit1_coef,
"Model 2" = fit2_coef),
.id = "Model")# A tibble: 8 × 8
Model term estimate std.error statistic p.value conf.low conf.high
<chr> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Model 1 Civic Duty 0.0285 0.0214 1.33 0.183 -0.0249 0.0818
2 Model 1 Hawthorne 0.0575 0.0221 2.60 0.00922 0.00234 0.113
3 Model 1 Self 0.0252 0.0211 1.19 0.234 -0.0277 0.0780
4 Model 1 Neighbors 0.102 0.0223 4.60 0.00000431 0.0468 0.158
5 Model 2 Civic Duty 0.0420 0.0205 2.05 0.0406 -0.00924 0.0933
6 Model 2 Hawthorne 0.0584 0.0212 2.76 0.00587 0.00546 0.111
7 Model 2 Self 0.0293 0.0203 1.44 0.149 -0.0214 0.0800
8 Model 2 Neighbors 0.114 0.0214 5.34 0.0000000975 0.0607 0.167 2つの表が1つとなり、Modelという列が追加される（これは.idで指定した名前）。そして、fit1_coefだった行は"Model 1"、fit2_coefだった行は"Model 2"が付く。ただし、これだけだと表が結合されて出力されるだけなので、fit_coefという名のオブジェクトとして作業環境内に格納しておく。
fit_coef <- bind_rows(list("Model 1" = fit1_coef,
"Model 2" = fit2_coef),
.id = "Model") それではfit_coefを使って、作図をしてみよう。コードは 図 14 と同じであるが、facet_wrap()レイヤーを追加する。これはグラフのファセット（facet）分割を意味し、ファセットとは「面」を意味する。()内には~分割の基準となる変数名を入れる。2つのモデルがあり、fit_coefだとModel列がどのモデルの推定値かを示している。
fit_coef |>
ggplot() +
geom_hline(yintercept = 0) +
geom_pointrange(aes(x = term, y = estimate,
ymin = conf.low, ymax = conf.high)) +
labs(x = "Treatments", y = "Average Treatment Effects") +
facet_wrap(~ Model) +
theme_bw(base_size = 12)
今回の結果だとモデル1もモデル2も推定値がほぼ同じである。ファセット分割の場合、小さい差の比較が難しいというデメリットがある。この場合、ファセット分割をせず、一つのファセットにpoint-rangeの色分けした方が読みやすくなる。point-rangeをModelの値に応じて色分けする場合、aes()内にcolor = Modelを追加する。
fit_coef |>
ggplot() +
geom_hline(yintercept = 0) +
geom_pointrange(aes(x = term, y = estimate,
ymin = conf.low, ymax = conf.high,
color = Model)) +
labs(x = "Treatments", y = "Average Treatment Effects") +
theme_bw(base_size = 12)
何かおかしい。point-rangeの横軸上の位置が同じということから重なってしまい、モデル1のpoint-rangeがよく見えない。これをずらすためにaes()の外側にposition = position_dodge2(1/2)を追加する。
fit_coef |>
ggplot() +
geom_hline(yintercept = 0) +
geom_pointrange(aes(x = term, y = estimate,
ymin = conf.low, ymax = conf.high,
color = Model),
position = position_dodge2(1/2)) +
labs(x = "Treatments", y = "Average Treatment Effects") +
theme_bw(base_size = 12)
これで図は完成だが、少し修正してみよう。{ggplot2}の場合、凡例は右側に表示されるが、これを下側へ移動させるためにはtheme()レイヤーを追加し、legend.position = "bottom"を指定する。また、モデル1とモデル2が具体的に何を意味するのかを明確に示したい。これはfit_coefのModel列を修正しても良いが、今回はscale_color_discrete()レイヤーで修正する例を紹介する。
fit_coef |>
ggplot() +
geom_hline(yintercept = 0) +
geom_pointrange(aes(x = term, y = estimate,
ymin = conf.low, ymax = conf.high,
color = Model),
position = position_dodge2(1/2)) +
labs(x = "Treatments", y = "Average Treatment Effects") +
scale_color_discrete(labels = c("Model 1" = "w/o Covariates",
"Model 2" = "w/ Covariates")) +
theme_bw(base_size = 12) +
theme(legend.position = "bottom")
脚注
pacman::p_load()は「{pacman}パッケージのp_load()関数」を意味する。このような書き方をすると、パッケージを読み込まず、関数を使うことができる。むろん、library(pacman)で予め{pacman}パッケージを読み込んでおくと、pacman::は省略し、p_load()だけでも問題ない。ただし、最初の1、2回程度しか使わないパッケージをわざわざ読み込んでおくのは目盛りの無駄遣いなので、このようにパッケージから直接呼び出したほうが効率が良い。↩︎パッケージのアップデートには
pacman::p_update()またはpacman::p_up()関数を使う。()内に何も入力しない場合、全パッケージがアップデートされる。↩︎Microsoft Excel形式（
.xls、.xlsx）は{readxl}パッケージのread_excel()関数を、Stata形式（.dta）は{heaven}のread_dta()（またはread_stata()）、SPSS形式（.sav）は{haven}のread_sav()（またはread_spss()）を使用する。ただし、データ分析界隈の標準は.csvフォーマットである。↩︎実習用データ（
rct_data.csv）はLMSから入手可能↩︎変数のデータ型はデータを出力する際、列名の下段に表示される。
<chr>は文字型、<dbl>と<int>は数値型、<fct>はfactor型である。他にもいくつかのデータ型がある。詳細は『私たちのR』の第8章を参照すること。↩︎RMarkdown内に埋め込むなら更に
style = "rmarkdown"を追加してみよう。ただし、Chunkオプションにresults = "asis"（Quartoなら#| results: "asis"）を付けること。↩︎A-B / A-C / A-D / A-E / B-C / B-D / B-E / C-D / C-E / D-E↩︎
中身が1、2、3、…であってもfactor型であれば1、2、3、…は数字でなく文字として認識される。↩︎
もっと使いやすい
round()があるが、round()の場合、丸めた結果が1.100なら1.1としか表記されない。表示される桁数を固定するためにはsprintf()を使う。↩︎従来のパイプ演算子（
%>%）を使う場合は_でなく、.を使う。↩︎%と%の差分の単位は%でなく、パーセントポイント（percent point; pp）を使用する。例えば35%から30%を引いたら、その結果は5%でなく、5%ポイントである。日本語では単純に「ポイント」と表記することが多いが、本講義ではppと表記する。他にもp.p、ppt、%p、%ptといった書き方がある。↩︎