Rの復習
パッケージ
通常、Rでのパッケージのインストールとアップーでとにはinstall.packages()
関数、読み込みにはlibrary()
、またはrequire()
関数を使う。また、R公式レポジトリにないパッケージは{devtools}か{remote}パッケージを使う。これらの関数を使い分けることは面倒なので、本講義ではこれらの処理を統合した{pacman}パッケージを使用する。まずは、{pacman}パッケージをインストールする。
パッケージを読み込む際、pacman::p_load(読み込むパッケージ名)
を入力する。インストールされていない場合は、自動的にCRANからダウンロード&インストールした上で読み込んでくれるので便利だ1。以下では本講義で使用するパッケージとして{tidyverse}、{summarytools}、{fastDummies}、{modelsummary}を読み込む。
CRANでなく、GitHub上で公開されているパッケージを使う場合はpacman::p_load_gh()
を使用する。()
の中には"ユーザー名/リポジトリ名"
を入力。たとえば、{BalanceR}の作成者のGitHubアカウント名はJaehyunSong
であり、{BalanceR}のリポジトリ名はBalanceR
だから、以下のように入力する。
データの読み込みと確認
.csv
形式のデータを読み込むにはread_csv()
関数を使用する。()
内には読み込むファイルのパスを"
で囲んで記入する。read_csv()
関数はファイルの読み込みのみの機能しか持たない。現在の作業環境内に読み込んだデータを格納するためには代入演算子<-
を使う。ここではdata
フォルダー内のrct_data.csv
を読み込み2、raw_df
という名のオブジェクトとしてく格納する。作業環境内のオブジェクトはRを再起動すると削除されるため、改めてパッケージ/データの読み込みが必要だ。
オブジェクトの中身を出力するためにはオブジェクト名を入力する。
# A tibble: 5,000 × 8
treatment gender yob hh_size voted2000 voted2002 voted2004 voted2006
<chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr>
1 Hawthorne female 1955 1 no yes no no
2 Hawthorne male 1969 2 no yes no no
3 Control female 1944 2 no no no no
4 Control male 1974 2 no no yes no
5 Neighbors female 1969 2 no no yes no
6 Control male 1942 1 no yes no no
7 Civic Duty male 1969 2 yes no no no
8 Hawthorne male 1966 2 no yes no yes
9 Control male 1949 2 no yes yes no
10 Control female 1958 2 no yes no no
# ℹ 4,990 more rows
表形式データの大きさ(行の列の数)の確認にはdim()
関数を使う。
表形式データの場合、各列には名前が付いており、それぞれが一つの変数に該当する。これら変数名のみの出力にはnames()
関数を使う。今回のデータだと、列の数が少ないこともあり、一画面に全列が表示されるが、数百列のデータとなると画面に収まらないので、変数名を確認しておくことを推奨する。
データハンドリングとパイプ演算子
パイプ演算子には{magrittr}パッケージが提供する%>%
とR 4.1から提供されるネイティブパイプ演算子の|>
がある。現在の主流は古くから使われてきた%>%
であるが、今後、|>
が主流になると考えられるため、本講義では|>
を使用する。しかし、多くの場合、|>
の代わりに%>%
を使っても同じ結果が得られる。
パイプ演算子はパイプ前のオブジェクトを、パイプ後の関数の第一引数として渡す単純な演算子だ。たとえば、列名を変更する関数はrename()
であるが、使い方はrenames(データ名, 新しい列名 = 既存の列名, ...)
である。raw_df
のgender
列の名前をfemale
に変更する場合は以下のように書く。
# A tibble: 5,000 × 8
treatment female yob hh_size voted2000 voted2002 voted2004 voted2006
<chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr>
1 Hawthorne female 1955 1 no yes no no
2 Hawthorne male 1969 2 no yes no no
3 Control female 1944 2 no no no no
4 Control male 1974 2 no no yes no
5 Neighbors female 1969 2 no no yes no
6 Control male 1942 1 no yes no no
7 Civic Duty male 1969 2 yes no no no
8 Hawthorne male 1966 2 no yes no yes
9 Control male 1949 2 no yes yes no
10 Control female 1958 2 no yes no no
# ℹ 4,990 more rows
ここで第1引数がraw_df
だが、パイプ演算子を使うと以下のようになり、人間にとって読みやすいコードになる。
# A tibble: 5,000 × 8
treatment female yob hh_size voted2000 voted2002 voted2004 voted2006
<chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr>
1 Hawthorne female 1955 1 no yes no no
2 Hawthorne male 1969 2 no yes no no
3 Control female 1944 2 no no no no
4 Control male 1974 2 no no yes no
5 Neighbors female 1969 2 no no yes no
6 Control male 1942 1 no yes no no
7 Civic Duty male 1969 2 yes no no no
8 Hawthorne male 1966 2 no yes no yes
9 Control male 1949 2 no yes yes no
10 Control female 1958 2 no yes no no
# ℹ 4,990 more rows
要するに、X |> Y
は「X
(の結果)を使ってY
を行う」ことを意味する。
続いて、変数のリコーディングをしてみよう。x
の値が"A"
なら1、それ以外は0のように、戻り値が2種類の場合のリコーディングにはif_else()
を使用する。書き方は以下の通りだ。
たとえば、raw_df
のgender
列の値が"female"
なら1、それ以外なら0とし、その結果をfemale
列として追加するコードは以下の通り。同値を意味する演算子が=
でなく、==
であることに注意すること(=
は<-
と同じ代入演算子であるが、Rでは代入演算子として=
より<-
の使用を推奨している)。
# 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 Hawthorne female 1955 1 no yes no no 1
2 Hawthorne male 1969 2 no yes no no 0
3 Control female 1944 2 no no no no 1
4 Control male 1974 2 no no yes no 0
5 Neighbors female 1969 2 no no yes no 1
6 Control male 1942 1 no yes no no 0
7 Civic Du… male 1969 2 yes no no no 0
8 Hawthorne male 1966 2 no yes no yes 0
9 Control male 1949 2 no yes yes no 0
10 Control female 1958 2 no yes no no 1
# ℹ 4,990 more rows
mutate()
は指定された列に対して何らかの処理を行い、その結果を新しい列として追加するか、上書きする関数である。このmutate()
関数の第1引数もデータであるため、以下のようにパイプ演算子を使うこともできる。
# 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 Hawthorne female 1955 1 no yes no no 1
2 Hawthorne male 1969 2 no yes no no 0
3 Control female 1944 2 no no no no 1
4 Control male 1974 2 no no yes no 0
5 Neighbors female 1969 2 no no yes no 1
6 Control male 1942 1 no yes no no 0
7 Civic Du… male 1969 2 yes no no no 0
8 Hawthorne male 1966 2 no yes no yes 0
9 Control male 1949 2 no yes yes no 0
10 Control female 1958 2 no yes no 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 Hawthorne female 1955 1 0 1 0 0 1
2 Hawthorne male 1969 2 0 1 0 0 0
3 Control female 1944 2 0 0 0 0 1
4 Control male 1974 2 0 0 1 0 0
5 Neighbors female 1969 2 0 0 1 0 1
6 Control male 1942 1 0 1 0 0 0
7 Civic Du… male 1969 2 1 0 0 0 0
8 Hawthorne male 1966 2 0 1 0 1 0
9 Control male 1949 2 0 1 1 0 0
10 Control female 1958 2 0 1 0 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 Hawthorne 1 1955 1 0 1 0 0
2 Hawthorne 0 1969 2 0 1 0 0
3 Control 1 1944 2 0 0 0 0
4 Control 0 1974 2 0 0 1 0
5 Neighbors 1 1969 2 0 0 1 0
6 Control 0 1942 1 0 1 0 0
7 Civic Duty 0 1969 2 1 0 0 0
8 Hawthorne 0 1966 2 0 1 0 1
9 Control 0 1949 2 0 1 1 0
10 Control 1 1958 2 0 1 0 0
# ℹ 4,990 more rows
以上のコードはデータを加工し、その結果を出力するだけであって、その結果を保存しない。もう一度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 Hawthorne female 1955 1 no yes no no
2 Hawthorne male 1969 2 no yes no no
3 Control female 1944 2 no no no no
4 Control male 1974 2 no no yes no
5 Neighbors female 1969 2 no no yes no
6 Control male 1942 1 no yes no no
7 Civic Duty male 1969 2 yes no no no
8 Hawthorne male 1966 2 no yes no yes
9 Control male 1949 2 no yes yes no
10 Control female 1958 2 no yes no 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 Hawthorne 1 1955 1 0 1 0 0
2 Hawthorne 0 1969 2 0 1 0 0
3 Control 1 1944 2 0 0 0 0
4 Control 0 1974 2 0 0 1 0
5 Neighbors 1 1969 2 0 0 1 0
6 Control 0 1942 1 0 1 0 0
7 Civic Duty 0 1969 2 1 0 0 0
8 Hawthorne 0 1966 2 0 1 0 1
9 Control 0 1949 2 0 1 1 0
10 Control 1 1958 2 0 1 0 0
# ℹ 4,990 more rows
ちなみに、across()
関数とラムダ式(無名関数)を組み合わせると以上のコードをより効率的に書くこともできる。across()
は強力な関数だが、初心者にはやや難しいかも知れない。詳細は『私たちのR』の第13.1章を参照されたい。
記述統計量
記述統計量の計算には{summarytools}のdescr()
関数が便利だ。descr(データ名)
を入力するだけで各変数の記述統計量が出力される。実際にやってみると分かるが、情報量がかなり多い。しかし、実際の論文では各変数の歪度や尖度まで報告することはあまりないだろう。ここではstats
引数を追加して、論文などでよく使う平均値("mean"
)、標準偏差("sd"
)、最小値("min"
)、最大値("max"
)、有効ケース数("n.valid"
)のみ出力する。
# descr()の第一引数は表形式データのオブジェクト名であるため、
# オブジェクト名 |> descr()
# のような書き方でもOK
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.18 0.26 0.39 0.41 0.33 1955.75
Std.Dev 0.50 0.79 0.44 0.49 0.49 0.47 14.41
Min 0.00 1.00 0.00 0.00 0.00 0.00 1910.00
Max 1.00 6.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)であるため、表示されない3。各グループがサンプルの何割かを計算するためには、treatment
変数をダミー変数へ変換する必要がある。ダミー変数の作成は面倒な作業であるが、{fastDummies}パッケージのdummy_cols()
を使えば簡単にできる。dummy_cols()
の中にはselect_columns = "ダミー化する列名"
を入れれば、当該変数をダミー変数へ変換し、新しい列として追加してくれる。それではtreatment
列をダミー化&追加し、その結果をdf
に上書きしてみよう。
# A tibble: 5,000 × 13
treatment female yob hh_size voted2000 voted2002 voted2004 voted2006
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Hawthorne 1 1955 1 0 1 0 0
2 Hawthorne 0 1969 2 0 1 0 0
3 Control 1 1944 2 0 0 0 0
4 Control 0 1974 2 0 0 1 0
5 Neighbors 1 1969 2 0 0 1 0
6 Control 0 1942 1 0 1 0 0
7 Civic Duty 0 1969 2 1 0 0 0
8 Hawthorne 0 1966 2 0 1 0 1
9 Control 0 1949 2 0 1 1 0
10 Control 1 1958 2 0 1 0 0
# ℹ 4,990 more rows
# ℹ 5 more variables: `treatment_Civic Duty` <int>, treatment_Control <int>,
# treatment_Hawthorne <int>, treatment_Neighbors <int>, treatment_Self <int>
画面には表示されないが、出力結果の下段を見るとtreatment_
で始まるいくつかの変数が追加されたことが分かる。ここでは"tretmant"
で始まる列のみを抽出つして確認してみよう。
# A tibble: 5,000 × 6
treatment `treatment_Civic Duty` treatment_Control treatment_Hawthorne
<chr> <int> <int> <int>
1 Hawthorne 0 0 1
2 Hawthorne 0 0 1
3 Control 0 1 0
4 Control 0 1 0
5 Neighbors 0 0 0
6 Control 0 1 0
7 Civic Duty 1 0 0
8 Hawthorne 0 0 1
9 Control 0 1 0
10 Control 0 1 0
# ℹ 4,990 more rows
# ℹ 2 more variables: treatment_Neighbors <int>, treatment_Self <int>
select()
関数内には抽出する列名を入力するだけで良い。たとえば、female
とyob
列を抽出するならselect(female, yob)
である。また、female
からvoted2006
までの意味でfemale:voted2006
のような書き方もできる。他にも上の例のようにstarts_with()
やends_with()
、contain()
を使って特定の文字列で始まる(で終わる、を含む)列を指定することもできる。一部の列を除外する場合は変数名の前に!
か-
を付ける。
とにかく、問題なくダミー化されていることが分かる。もう一度記述統計量を出してみよう。descr()
は仕様上、出力される変数の順番はアルファベット順になるが、ここでは元の順番を維持するためにorder = "p"
を追加する。また、通常の記述統計表が、先ほど見たものとは違って、各行が変数を、列は記述統計量を表す場合が多い。このように行と列を交換するためにはtranspose = TRUE
を追加する4。
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 | 1955.75 | 14.41 | 1910.00 | 1986.00 | 5000.00 |
hh_size | 2.18 | 0.79 | 1.00 | 6.00 | 5000.00 |
voted2000 | 0.26 | 0.44 | 0.00 | 1.00 | 5000.00 |
voted2002 | 0.39 | 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.56 | 0.50 | 0.00 | 1.00 | 5000.00 |
treatment_Hawthorne | 0.11 | 0.32 | 0.00 | 1.00 | 5000.00 |
treatment_Neighbors | 0.11 | 0.31 | 0.00 | 1.00 | 5000.00 |
treatment_Self | 0.11 | 0.31 | 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] |
|
76 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.0.1 (R version 4.3.1)
2023-09-25
バランスチェック
バランスチェックの簡単な方法はグループごとに処置前変数(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.517 1957. 2.22 0.239 0.374 0.394
2 Control 0.491 1956. 2.17 0.261 0.390 0.411
3 Hawthorne 0.507 1956. 2.23 0.257 0.385 0.416
4 Neighbors 0.507 1955. 2.13 0.247 0.406 0.439
5 Self 0.526 1956. 2.19 0.253 0.383 0.401
それぞれの変数の平均値は非常に似ているため、無作為割当が成功したと考えられる。しかし、変数の単位によって判断が難しいかも知れない。たとえば、2つのグループがあり、年齢の平均値の差は3、世帯規模のそれは2だとする。これを見ると年齢の方がよりバランスが取れていないようにも見えるが、年齢の幅は数十であるに対し、世帯規模はせいぜい5〜6程度であろう。したがって、各変数のばらつきまで考慮した比較が適切であり、その方法の一つが標準化バイアス(=標準化差分)である。
標準化差分を計算する便利パッケージ、{BalanceR}を使ってみよう。第1引数はデータだから、パイプで渡せば良い。BalanceR()
内にはgroup
引数にグループ識別変数を、cov
には処置前変数のベクトルを入れる。
blc_chk <- df |>
BalanceR(group = treatment,
cov = c(female, yob, hh_size, voted2000, voted2002, voted2004))
blc_chk
Covariate Mean:Civic Duty SD:Civic Duty Mean:Control SD:Control
1 female 0.517 0.500 0.491 0.500
2 yob 1957.304 14.515 1955.555 14.453
3 hh_size 2.225 0.810 2.169 0.777
4 voted2000 0.239 0.427 0.261 0.439
5 voted2002 0.374 0.484 0.390 0.488
6 voted2004 0.394 0.489 0.411 0.492
Mean:Hawthorne SD:Hawthorne Mean:Neighbors SD:Neighbors Mean:Self SD:Self
1 0.507 0.500 0.507 0.500 0.526 0.500
2 1955.909 14.638 1954.976 14.510 1955.719 13.639
3 2.231 0.809 2.129 0.786 2.193 0.772
4 0.257 0.437 0.247 0.432 0.253 0.435
5 0.385 0.487 0.406 0.491 0.383 0.487
6 0.416 0.493 0.439 0.497 0.401 0.491
SB:Civic Duty-Control SB:Civic Duty-Hawthorne SB:Civic Duty-Neighbors
1 5.299 2.062 2.153
2 12.082 9.573 16.047
3 7.102 -0.723 12.018
4 -5.157 -4.225 -1.860
5 -3.353 -2.249 -6.546
6 -3.420 -4.405 -9.094
SB:Civic Duty-Self SB:Control-Hawthorne SB:Control-Neighbors SB:Control-Self
1 -1.746 -3.236 -3.145 -7.046
2 11.255 -2.437 3.997 -1.173
3 3.994 -7.847 5.060 -3.199
4 -3.260 0.931 3.296 1.896
5 -1.897 1.103 -3.191 1.456
6 -1.435 -0.985 -5.669 1.985
SB:Hawthorne-Neighbors SB:Hawthorne-Self SB:Neighbors-Self
1 0.090 -3.809 -3.899
2 6.404 1.341 -5.281
3 12.763 4.739 -8.258
4 2.365 0.965 -1.400
5 -4.295 0.353 4.648
6 -4.684 2.970 7.656
ちなみに、df
内にfemale
からvoted2004
は連続している(names(df)
で確認してみよう)。この場合は以下のように(female:voted2004
)書き換えることもできる。
標準化差分(標準化バイアス)を用いたバランスチェックはそれぞれのペアごとに計算を行うため、グループが多い場合は凡例が圧迫される場合が多い。しかし、重要なのは標準化差分の最大値だろう。ペア1、2、3でバランスが取れても、ペア4のバランスが取られていない場合は無意味だからだ。また、標準化差分の場合、符号の意味はなく、絶対値が重要だ。また、バランスチェックにおいてグループごとの平均値や標準偏差は不要である。ここでsummary()
関数を使うと、絶対値が最も大きい標準化差分のみ出力される。
Covariate Abs_Maximum_SB
1 female 7.046
2 yob 16.047
3 hh_size 12.763
4 voted2000 5.157
5 voted2002 6.546
6 voted2004 9.094
plot()
関数を使えば、これらの結果を可視化することもできる。
先ほど述べたようにバランスチェックで重要なのは絶対値が最も大きい標準化差分である。plot()
内にsimplify = TRUE
を指定すれば最大値のみ表示され、更にabs = TRUE
にすると絶対値へ変換される。また、垂直のガイドラインはvline
引数で変更できる。
処置効果の確認
グループごとの応答変数の平均値
処置効果を確認するためには各グループごとの応答変数(ここではvoted2006
)の平均値を計算し、処置群の平均値から統制群の平均値を引く必要がある。まずは、特定の変数の平均値を計算する方法について紹介する。データ内にある特定の変数の平均値を計算するためにはsummarise()
関数内に平均値を求めるmean()
関数を入れる。たとえば、df
のvoted2006
の平均値を計算するコードは以下の通りである。
# A tibble: 1 × 1
`mean(voted2006, na.rm = TRUE)`
<dbl>
1 0.325
na.rm = TRUE
は「欠損値があれば、それを除外する」を意味し、指定されていない場合(=既定値)はFALSE
になる。今回は欠損値がないものの、念の為に入れておく。
出力結果を見ると、平均値が表示される列の名前が`mean(voted2006, na.rm = TRUE)`
となっており、非常に見にくい。この場合、以下のようにmean()
の前に出力される列名を予め指定することもできる。
# A tibble: 1 × 1
Outcome
<dbl>
1 0.325
我々が知りたいのはvoted2006
の平均値でなく、グループごとの平均値だろう。被験者がどのグループに属しているかわ示す変数はtreatment
であるが、summarise()
にデータを渡す前にgroup_by()
変数を使うと、グループごとに計算を行い、その結果を返す。
# A tibble: 5 × 2
treatment Outcome
<chr> <dbl>
1 Civic Duty 0.317
2 Control 0.299
3 Hawthorne 0.339
4 Neighbors 0.421
5 Self 0.364
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.317
2 Control 0.299
3 Hawthorne 0.339
4 Neighbors 0.421
5 Self 0.364
ここで一つ注目したいのが、グループの表示順番である。変数のデータ型が文字型だと(Rコンソール上でclass(df$treatment)
を入力するか、df
の出力画面でtreatment
の下に<chr>
と表示されていることで確認できる)、今のようにアルファベット順で表示される。しかし、統制群は最初か最後に来るのが通例である。この順番をアルファベット順でなく、任意の順番にするためにはtreatment
変数をfactor型変数へ変換する必要がある。Factor型は「順序付きの文字型変数」だと理解しても良い5。列の追加・上書き(今回は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
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Hawthorne 1 1955 1 0 1 0 0
2 Hawthorne 0 1969 2 0 1 0 0
3 Control 1 1944 2 0 0 0 0
4 Control 0 1974 2 0 0 1 0
5 Neighbors 1 1969 2 0 0 1 0
6 Control 0 1942 1 0 1 0 0
7 Civic Duty 0 1969 2 1 0 0 0
8 Hawthorne 0 1966 2 0 1 0 1
9 Control 0 1949 2 0 1 1 0
10 Control 1 1958 2 0 1 0 0
# ℹ 4,990 more rows
# ℹ 5 more variables: `treatment_Civic Duty` <int>, treatment_Control <int>,
# treatment_Hawthorne <int>, treatment_Neighbors <int>, treatment_Self <int>
treatment
列名の下が<fct>
となっていることが分かる。これはtreatment
列のデータ型がfactor型であることを意味する。問題なく動くことが確認できたので、df
を上書きしよう。
それでは、改めてグループごとの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.299
2 Civic Duty 0.317
3 Self 0.364
4 Neighbors 0.421
5 Hawthorne 0.339
今回は統制群は最初に出力されていることが確認できる。
それではこの結果をグラフとして示してみよう。作図には{ggplot2}パッケージを使う。まずはout_mean_df
をggplot()
関数に渡す。ggplot()
関数以降は、+
演算子を使ってレイヤーを足していくこととなる。棒グラフのレイヤーはgeom_bar()
関数であり、その中にaes()
関数を入れる。aes()
の中には棒グラフの作図に必要な情報を入れる必要がある(これをマッピング(mapping)と呼ぶ)。棒グラフを作成するために必要な最低限の情報とは各棒の横軸上の位置(x
)と棒の高さ(y
)だ。今回は横軸がグループ名、縦軸が平均値となる棒グラフを作る。aes()
外側にはstat = "identity"
を忘れずに付けること。
続いて、このグラフの見た目を調整してみよう。
out_mean_df |>
ggplot() +
geom_bar(aes(x = Groups, y = Outcome), stat = "identity") +
# 縦軸(y軸)のラベルを変更する
labs(y = "Mean(Outcome)") +
# grayテーマ(デフォルトのテーマ)を使用し、フォントサイズは14
theme_gray(base_size = 14)
また、geom_label()
レイヤーを足すと、棒の上にラベルを付けることもできる。ラベルに必要な情報は各ラベルの横軸上の位置(x
)、縦軸上の位置(y
)、ラベルの表示内容(label
)だ。今回のラベルは平均値の具体的な数値を入れてみよう。
out_mean_df |>
ggplot() +
geom_bar(aes(x = Groups, y = Outcome), stat = "identity") +
geom_label(aes(x = Groups, y = Outcome, label = Outcome)) +
labs(y = "Mean(Outcome)") +
theme_gray(base_size = 14)
小数点が長すぎるので3桁まで表示としよう。ここではsprintf()
を使用する。使い方が簡単とは言えないが、覚える必要はなく、必要な時にググるか、本資料のコードをコピペすれば良い6。
out_mean_df |>
ggplot() +
geom_bar(aes(x = Groups, y = Outcome), stat = "identity") +
# 2桁までなら %.3f を %.2f に変更
geom_label(aes(x = Groups, y = Outcome, label = sprintf("%.3f", Outcome))) +
labs(y = "Mean(Outcome)") +
theme_gray(base_size = 14)
これで可視化ができた。ただし、以上のコードには改善の余地がある。geom_bar()
とgeom_label()
内のaes()
関数に注目して欲しい。よく見るとx
とy
と同じだろう。geom_*()
が共有するマッピングがあれば、ggplot()
内で指定することでコードを効率化することもできる。
統計的推定(単回帰分析)
これまでの作業はグループごとの応答変数の平均値であって、処置効果ではない。処置効果を計算するためには処置群の平均値から統制群の平均値を引く必要がある。たとえば、Civic Dutyはがき群の平均値は約0.317、統制群のそれは0.299であるため、Civic Dutyはがきの処置効果は約0.018である。しかし、これを各グループごとに計算することは面倒だし、何よりも得られた値が点推定値だという限界がある。得られた処置効果の不確実性は計算できない。
ここで有効なのが線形回帰分析である。回帰分析を行うことで処置効果の点推定値のみならず、不確実性の指標である標準誤差も計算され、区間推定や統計的仮説検定も可能となる。線形回帰分析の関数はlm()
だ。第1引数としては回帰式であり、応答変数 ~ 説明変数
と表記する。第2引数はdata
であり、回帰式で指定した変数が入っているデータ名を指定する。回帰分析の結果は名前を付けてオブジェクトとして格納し、summary()
関数を使うと、詳細が確認できる。
Call:
lm(formula = voted2006 ~ treatment, data = df)
Residuals:
Min 1Q Median 3Q Max
-0.4206 -0.3166 -0.2985 0.6357 0.7015
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.298524 0.008864 33.680 < 2e-16 ***
treatmentCivic Duty 0.018085 0.021355 0.847 0.3971
treatmentSelf 0.065789 0.022002 2.990 0.0028 **
treatmentNeighbors 0.122037 0.022053 5.534 3.3e-08 ***
treatmentHawthorne 0.040637 0.021447 1.895 0.0582 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4671 on 4995 degrees of freedom
Multiple R-squared: 0.007123, Adjusted R-squared: 0.006328
F-statistic: 8.959 on 4 and 4995 DF, p-value: 3.312e-07
ちなみに、これもパイプ演算子を使うことができる。ただし、第1引数として渡すパイプ演算子の特徴上、そのまま使うことはできない。なぜならlm()
関数の第1引数はデータでなく、回帰式(formula型)だから。この場合はプレースホルダー(place holder)を指定する必要がある。パイプ前のオブジェクトが入る位置を任意に指定することであり、_
を使う。%>%
演算子を使う場合は_
でなく、.
を使う。上記のコードと以下のコードは同じコードとなる。プレースホルダーは自分が使うパイプ演算子によって使い分けること。
Factor型、または文字型変数が説明変数の場合、自動的にダミー変数として処理され、Factor型の場合、最初の水準(ここでは"Control"
)がベースカテゴリとなる。説明変数が文字型ならアルファベット順で最初の水準がベースカテゴリとなり、今回の例だと"Civic Duty"
がベースカテゴリとなる。処置効果は「統制群に比べて〜」が重要となるので、数値型以外の説明変数は予めfactor化しておいた方が望ましい。
Civic Dutyの推定値は約0.018であり、これは統制群に比べ、Civic Duty群のvoted2006
の平均値は約0.018高いことを意味する。応答変数が0、1であるため、これを割合(=投票率)で換算すると、約1.8%p高いことを意味する。つまり、Civic Dutyのはがきをもらった被験者はそうでない被験者に比べて投票率が約1.8%p高いことを意味する。他の推定値も同じやり方で解釈すれば良い。
それではこれらの処置効果が統計的に有意なものかを確認してみよう。統計的有意か否かを判定するためには有意と非有意の境界線が必要である、これは通常、有意水準(significance level; \(\alpha\))と呼ばれる。この有意水準は分析者が決めるものであるが、社会科学で広く使われる基準は\(\alpha = 0.05\)、つまり5%だ。分析結果の画面にはPr(>|t|)
列が表示されているが、これが\(p\)値と呼ばれるもので、これが0.05を下回る場合、統計的に有意と判定する。もし、\(\alpha = 0.1\)を採用するなら、\(p < 0.1\)の場合において統計的に有意と判定する。Civic Dutyの\(p\)値は3.3e-08であり、これは\(3.3 \times 10^{-8}\)を意味する。\(10^{-1}\)は0.1、\(10^{-2}\)は0.01であることを考えると非常に小さい数値であり、統計的に有意であると考えられる。また、\(p\)値が一定値以下であれば< 2e-16
と表示される。今回の結果は、2つの処置群(SelfとNeighbors)において処置効果は統計的に有意であると判定できよう。
今回はサンプルサイズ34万から無作為抽出した\(n\)=5000のデータを使ったが、実際の論文ではいずれも統計的に有意である。なぜなら統計的有意性を判定する\(p\)値は処置効果の大きさ以外にも、推定値のばらつき具合(=標準誤差)の影響を受けるが、この標準誤差はサンプルサイズに反比例するからだ。ここで重要なことが分かる。それは「統計的に非有意な処置効果」は「処置効果がない」ことを意味しないことだ。たとえば、我々が日本人から男女3名ずつ抽出し、平均身長が同じだとしよう。ここから「母集団において男女の身長差はない」と解釈できるだろうか。それは全く誤りであろう。統計的に非有意な結果から分かるのは「現在のサンプル/手法/モデルにおいて男女間に身長差があるとは言えない」程度である。もう一回、男女3人ずつ抽出して比較してみれば差はあるかも知れない。また、今回はたまたま平均身長が同じだけであって、サンプルサイズが大きければ統計的に有意な差は得られるかも知れない。
一方、統計的に有意な差が得られたのであれば、「母集団において男女の身長差がある」と解釈できるようになる。これはもう一回、男女3人ずつ無作為抽出しても差が確認されることを意味する。むろん、100%ではない。もう一回やってみたところ、差がない可能性もある。ただし、その可能性は非常に小さく、その可能性が小さいことを、統計学では\(p\)値が小さいと言う。\(p\)値の厳密な意味は「帰無仮説が正しいと仮定した場合、今回得られた統計量と同じか、より極端な結果が得られる確率」である。身長差の例だと、帰無仮説は「男女間に身長差はない」であり、統計量とは\(t\)検定や単回帰分析の場合、「\(t\)統計量」を意味する。詳細はp値に関するアメリカ統計学会の声明や、矢内・SONG (2019)を参照すること。
話がややずれたが、Rの実習に戻そう。続いて、この結果を可視化してみよう。ここでも{ggplot2}パッケージを使って可視化をするが、{ggplot2}で使用可能なオブジェクトは表形式のデータである。Rコンソール上でclass(オブジェクト名)
を入力すると、データのクラスが出力されるが、このクラスに"data.frame"
があれば、{ggplot2}で使用できる。たとえば、fit1
オブジェクトのクラスは"lm"
であるため、そのまま{ggplot2}で使うことはできない。
推定結果を表形式に変換するためには{broom}パッケージのtidy()
関数が便利だ。使い方は簡単でtidy()
内に回帰分析の推定結果が格納されたオブジェクトを入れるだけである。ただし、デフォルトの設定では95%信頼区間が表示されないため、中にはconf.int = TRUE
を追加しておく必要がある。
# 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.299 0.00886 33.7 2.80e-224 0.281 0.316
2 treatmentCivic Duty 0.0181 0.0214 0.847 3.97e- 1 -0.0238 0.0599
3 treatmentSelf 0.0658 0.0220 2.99 2.80e- 3 0.0227 0.109
4 treatmentNeighbors 0.122 0.0221 5.53 3.30e- 8 0.0788 0.165
5 treatmentHawthorne 0.0406 0.0214 1.89 5.82e- 2 -0.00141 0.0827
[1] "tbl_df" "tbl" "data.frame"
fit1_coef
のクラスに"data.frame"
が含まれているので、これを使って作図することができる。
作図する前に、fit1_coef
の加工しておきたい。それぞれの係数(estimate
列)は処置効果を表しているが、切片("(Intercept)"
)の推定値は処置効果とは無関係である。したがって、予め切片の行を除外しておきたい。特定の行を残したり、除外する関数はfilter()
である。今回はterm
列の値が"(Intercept)"
ではない行を残したいので、同値演算子(==
)の否定を意味する!=
演算子を使用する。
# 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.0181 0.0214 0.847 3.97e-1 -0.0238 0.0599
2 treatmentSelf 0.0658 0.0220 2.99 2.80e-3 0.0227 0.109
3 treatmentNeighbors 0.122 0.0221 5.53 3.30e-8 0.0788 0.165
4 treatmentHawthorne 0.0406 0.0214 1.89 5.82e-2 -0.00141 0.0827
それでは作図に入ろう。処置効果を示す場合は、点推定値以外にもその不確実性を示すのは一般的である。不確実性の指標として幅広く使われるのは標準誤差(standard error; 標準偏差ではない)であるが、可視化の際にはこの標準誤差に基づき計算した信頼区間を示すのが一般的だ。有意水準が5%であれば、95%信頼区間を示し、10%なら90%信頼区間を用いる。
点と区間を同時に示すプロットが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のところに水平線が描ける。
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)
まだ気になる点がある。それは横軸の目盛りラベルにtreatment
という不要な情報がある点だ。これは作図の時点で修正することも可能だが、まずはdf
のterm
変数の値を修正する方法を紹介する。変数の値を修正する時にはrecode()
関数を使用する。第1引数はリコーディングする変数名であり、引き続き"元の値" = "新しい値"
を指定すれば良い。スペルミスに注意すること。
fit1_coef <- fit1_coef |>
mutate(term = recode(term,
"treatmentCivic Duty" = "Civic Duty",
"treatmentHawthorne" = "Hawthorne",
"treatmentNeighbors" = "Neighbors",
"treatmentSelf" = "Self"))
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.0181 0.0214 0.847 0.397 -0.0238 0.0599
2 Self 0.0658 0.0220 2.99 0.00280 0.0227 0.109
3 Neighbors 0.122 0.0221 5.53 0.0000000330 0.0788 0.165
4 Hawthorne 0.0406 0.0214 1.89 0.0582 -0.00141 0.0827
以上の作業はterm
列の各値から"treatment"
文字を""
に置換することなので、文字列を置換する関数であるstr_replace()
を使えば、より短くすることができる。
fit1_coef
も修正できたので、 図 9 と同じコードでもう一度作図してみよう。
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()
を使ってみよう。
# 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.0181 0.0214 0.847 0.397 -0.0238 0.0599
2 Self 0.0658 0.0220 2.99 0.00280 0.0227 0.109
3 Neighbors 0.122 0.0221 5.53 0.0000000330 0.0788 0.165
4 Hawthorne 0.0406 0.0214 1.89 0.0582 -0.00141 0.0827
それでは、 図 10 と同じコードでもう一度作図してみよう。
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)
これで処置効果の可視化もバッチリだ。
多重比較の問題
グループが2つ、つまり統制群と統制群のみが存在する場合、我々が比較を行う回数は1回のみである(統制群 - 処置群)。しかし、今回のデータの場合、処置群は4つである。これは比較を4回行うことを意味する。具体的には「統制群 - 処置群1」、「統制群 - 処置群2」、「統制群 - 処置群3」、「統制群 - 処置群4」だ。比較を繰り返すほど、統計的に有意な結果が得られる可能性は高い。極端な話、1000回程度検定を繰り返せば、本当は効果がなくてもたまたま統計的に有意な結果が何回かは得られるだろう。これが多重検定(multiple testing)の問題である。したがって、比較の回数が多くなるにつれ、統計的有意性検定にも何らかのペナルティーを課す必要がある。
多重比較におけるペナルティーの付け方はいくつかあるが、ここでは最も保守的な(=研究者にとって都合の悪い)補正法であるボンフェローニ補正(Bonferroni correction)を紹介する。これは非常に単純で、\(p\)値や信頼区間を計算する際、「統計的有意」と判定されるハードルを上げる方法である。予め決めておいた有意水準(\(\alpha\))が0.05で、比較の回数が4回であれば、\(p\)値が\(0.05 \times \frac{1}{4} = 0.0125\)を下回る場合において「5%水準で有意である」と判定する。信頼区間でいえば通常の95%信頼区間(1 - 0.05)でなく、98.75%信頼区間(1 - 0.0125)を使うこととなる。この結果、統計的に有意な結果が得られたら「1.25%水準で〜」と解釈するのではなく、「5%水準で〜」と解釈する必要がある。
95%以外の信頼区間を求めるのは簡単で、tidy()
関数内にconf.level
を修正すれば良い。指定されていない場合はデフォルトで0.95が割り当てられているが、これを0.9875と修正する。
# 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.299 0.00886 33.7 2.80e-224 0.276 0.321
2 treatmentCivic Duty 0.0181 0.0214 0.847 3.97e- 1 -0.0353 0.0714
3 treatmentSelf 0.0658 0.0220 2.99 2.80e- 3 0.0108 0.121
4 treatmentNeighbors 0.122 0.0221 5.53 3.30e- 8 0.0669 0.177
5 treatmentHawthorne 0.0406 0.0214 1.89 5.82e- 2 -0.0130 0.0942
それでは 図 11 と同じ図を作ってみよう。まず、切片の行を除外するが、ここでは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
と指定すれば良い。
# 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.0181 0.0214 0.847 3.97e-1 -0.0353 0.0714
2 treatmentSelf 0.0658 0.0220 2.99 2.80e-3 0.0108 0.121
3 treatmentNeighbors 0.122 0.0221 5.53 3.30e-8 0.0669 0.177
4 treatmentHawthorne 0.0406 0.0214 1.89 5.82e-2 -0.0130 0.0942
続いて、term
変数の値から"treatment"
の文字を除去し、fit1_coef
での出力順番でterm
をfactor化する。
fit1_coef <- fit1_coef |>
mutate(term = recode(term,
"treatmentCivic Duty" = "Civic Duty",
"treatmentHawthorne" = "Hawthorne",
"treatmentNeighbors" = "Neighbors",
"treatmentSelf" = "Self"),
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.0181 0.0214 0.847 0.397 -0.0353 0.0714
2 Self 0.0658 0.0220 2.99 0.00280 0.0108 0.121
3 Neighbors 0.122 0.0221 5.53 0.0000000330 0.0669 0.177
4 Hawthorne 0.0406 0.0214 1.89 0.0582 -0.0130 0.0942
最後に 図 11 と同じコードで作図する。
統計的推定(重回帰分析)
今回の例は無作為割当が成功しており、処置前変数の偏りは見られない。しかし、何らかの理由で処置前変数の偏りが生じる場合がある。その「何らかの理由」が応答変数にまで影響を与えるのであれば、それは交絡変数(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.7573 -0.3447 -0.1953 0.5184 0.9506
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.0886118 0.9164461 6.644 3.39e-11 ***
treatmentCivic Duty 0.0307108 0.0205474 1.495 0.135074
treatmentSelf 0.0703356 0.0211562 3.325 0.000892 ***
treatmentNeighbors 0.1147107 0.0212088 5.409 6.65e-08 ***
treatmentHawthorne 0.0423327 0.0206241 2.053 0.040165 *
female -0.0240666 0.0127276 -1.891 0.058696 .
yob -0.0030326 0.0004705 -6.445 1.27e-10 ***
hh_size 0.0025215 0.0085056 0.296 0.766896
voted2000 0.0666981 0.0147950 4.508 6.69e-06 ***
voted2002 0.1644528 0.0133278 12.339 < 2e-16 ***
voted2004 0.1583835 0.0130339 12.152 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.449 on 4989 degrees of freedom
Multiple R-squared: 0.08363, Adjusted R-squared: 0.0818
F-statistic: 45.53 on 10 and 4989 DF, p-value: < 2.2e-16
{modelsummary}パッケージのmodelsummary()
関数を使えば、推定結果がより見やすくなる。
(1) | |
---|---|
(Intercept) | 6.089 |
(0.916) | |
treatmentCivic Duty | 0.031 |
(0.021) | |
treatmentSelf | 0.070 |
(0.021) | |
treatmentNeighbors | 0.115 |
(0.021) | |
treatmentHawthorne | 0.042 |
(0.021) | |
female | −0.024 |
(0.013) | |
yob | −0.003 |
(0.000) | |
hh_size | 0.003 |
(0.009) | |
voted2000 | 0.067 |
(0.015) | |
voted2002 | 0.164 |
(0.013) | |
voted2004 | 0.158 |
(0.013) | |
Num.Obs. | 5000 |
R2 | 0.084 |
R2 Adj. | 0.082 |
AIC | 6195.0 |
BIC | 6273.2 |
Log.Lik. | −3085.506 |
F | 45.533 |
RMSE | 0.45 |
また、複数のモデルをlist()
関数でまとめると、モデル間比較もできる。
w/o Covariates | w/ Covariates | |
---|---|---|
(Intercept) | 0.299 | 6.089 |
(0.009) | (0.916) | |
treatmentCivic Duty | 0.018 | 0.031 |
(0.021) | (0.021) | |
treatmentSelf | 0.066 | 0.070 |
(0.022) | (0.021) | |
treatmentNeighbors | 0.122 | 0.115 |
(0.022) | (0.021) | |
treatmentHawthorne | 0.041 | 0.042 |
(0.021) | (0.021) | |
female | −0.024 | |
(0.013) | ||
yob | −0.003 | |
(0.000) | ||
hh_size | 0.003 | |
(0.009) | ||
voted2000 | 0.067 | |
(0.015) | ||
voted2002 | 0.164 | |
(0.013) | ||
voted2004 | 0.158 | |
(0.013) | ||
Num.Obs. | 5000 | 5000 |
R2 | 0.007 | 0.084 |
R2 Adj. | 0.006 | 0.082 |
AIC | 6584.0 | 6195.0 |
BIC | 6623.1 | 6273.2 |
Log.Lik. | −3285.982 | −3085.506 |
F | 8.959 | 45.533 |
RMSE | 0.47 | 0.45 |
modelsummary()
は推定値と標準誤差(カッコ内)が別々の行として出力する。これを一行でまとめるためには、以下のようにコードを修正する。
modelsummary(list("w/o Covariates" = fit1, "w/ Covariates" = fit2),
estimate = "{estimate} ({std.error})",
statistic = NULL)
w/o Covariates | w/ Covariates | |
---|---|---|
(Intercept) | 0.299 (0.009) | 6.089 (0.916) |
treatmentCivic Duty | 0.018 (0.021) | 0.031 (0.021) |
treatmentSelf | 0.066 (0.022) | 0.070 (0.021) |
treatmentNeighbors | 0.122 (0.022) | 0.115 (0.021) |
treatmentHawthorne | 0.041 (0.021) | 0.042 (0.021) |
female | −0.024 (0.013) | |
yob | −0.003 (0.000) | |
hh_size | 0.003 (0.009) | |
voted2000 | 0.067 (0.015) | |
voted2002 | 0.164 (0.013) | |
voted2004 | 0.158 (0.013) | |
Num.Obs. | 5000 | 5000 |
R2 | 0.007 | 0.084 |
R2 Adj. | 0.006 | 0.082 |
AIC | 6584.0 | 6195.0 |
BIC | 6623.1 | 6273.2 |
Log.Lik. | −3285.982 | −3085.506 |
F | 8.959 | 45.533 |
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",
"treatmentSelf" = "Self",
"treatmentNeighbors" = "Neighbors",
"treatmentHawthorne" = "Hawthorne",
"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.299 (0.009) | 6.089 (0.916) |
Civic Duty | 0.018 (0.021) | 0.031 (0.021) |
Self | 0.066 (0.022) | 0.070 (0.021) |
Neighbors | 0.122 (0.022) | 0.115 (0.021) |
Hawthorne | 0.041 (0.021) | 0.042 (0.021) |
Female | −0.024 (0.013) | |
Year of Birth | −0.003 (0.000) | |
Household Size | 0.003 (0.009) | |
Voted (2000) | 0.067 (0.015) | |
Voted (2002) | 0.164 (0.013) | |
Voted (2004) | 0.158 (0.013) | |
Num.Obs. | 5000 | 5000 |
R2 | 0.007 | 0.084 |
R2 Adj. | 0.006 | 0.082 |
AIC | 6584.0 | 6195.0 |
BIC | 6623.1 | 6273.2 |
Log.Lik. | −3285.982 | −3085.506 |
F | 8.959 | 45.533 |
RMSE | 0.47 | 0.45 |
処置効果に注目すると、共変量の有無が推定結果に影響をほぼ与えないことが分かる。これは無作為割当に成功したことを意味する。
番外編
modelsummary()
を使えば、複数のモデルの推定結果を一つの表としてまとめられる。しかし、図の場合はどうだろう。共変量なしモデルとありモデルを 図 12 のように一つにまとめることはできるだろうか。もちろん出来る。
まず、重回帰分析を行った結果(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",
"treatmentNeighbors" = "Neighbors",
"treatmentSelf" = "Self"),
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.0307 0.0205 1.49 0.135 -0.0206 0.0821
2 Self 0.0703 0.0212 3.32 0.000892 0.0175 0.123
3 Neighbors 0.115 0.0212 5.41 0.0000000665 0.0617 0.168
4 Hawthorne 0.0423 0.0206 2.05 0.0402 -0.00920 0.0939
処置効果の推定値や標準誤差などが異なるが、構造としては同じである。続いて、bind_rows()
を用い、この2つのデータを一つの表として結合する。2つの表はlist()
関数でまとめるが、それぞれ"モデル名" = データ名
と指定する。最後に、.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.0181 0.0214 0.847 3.97e-1 -0.0353 0.0714
2 Model 1 Self 0.0658 0.0220 2.99 2.80e-3 0.0108 0.121
3 Model 1 Neighbors 0.122 0.0221 5.53 3.30e-8 0.0669 0.177
4 Model 1 Hawthorne 0.0406 0.0214 1.89 5.82e-2 -0.0130 0.0942
5 Model 2 Civic Duty 0.0307 0.0205 1.49 1.35e-1 -0.0206 0.0821
6 Model 2 Self 0.0703 0.0212 3.32 8.92e-4 0.0175 0.123
7 Model 2 Neighbors 0.115 0.0212 5.41 6.65e-8 0.0617 0.168
8 Model 2 Hawthorne 0.0423 0.0206 2.05 4.02e-2 -0.00920 0.0939
2つの表が1つとなり、Model
という列が追加される(これは.id
で指定した名前)。そして、fit1_coef
だった行は"Model 1"
、fit2_coef
だった行は"Model 2"
が付く。ただし、これだけだと表が結合されて出力されるだけなので、fit_coef
という名のオブジェクトとして作業環境内に格納しておく。
それではfit_coef
を使って、作図をしてみよう。コードは 図 12 と同じであるが、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")