Rの復習

パッケージ

　通常、Rでのパッケージのインストールとアップーでとにはinstall.packages()関数、読み込みにはlibrary()、またはrequire()関数を使う。また、R公式レポジトリにないパッケージは{devtools}か{remote}パッケージを使う。これらの関数を使い分けることは面倒なので、本講義ではこれらの処理を統合した{pacman}パッケージを使用する。まずは、{pacman}パッケージをインストールする。

# NIIオンライン分析システムを利用する場合、導入済み
install.packages("pacman")

　パッケージを読み込む際、pacman::p_load(読み込むパッケージ名)を入力する。インストールされていない場合は、自動的にCRANからダウンロード&インストールした上で読み込んでくれるので便利だ¹。以下では本講義で使用するパッケージとして{tidyverse}、{summarytools}、{fastDummies}、{modelsummary}を読み込む。

pacman::p_load(tidyverse, summarytools, fastDummies,
               modelsummary, broom)

　CRANでなく、GitHub上で公開されているパッケージを使う場合はpacman::p_load_gh()を使用する。()の中には"ユーザー名/リポジトリ名"を入力。たとえば、{BalanceR}の作成者のGitHubアカウント名はJaehyunSongであり、{BalanceR}のリポジトリ名はBalanceRだから、以下のように入力する。

pacman::p_load_gh("JaehyunSong/BalanceR")

データの読み込みと確認

　.csv形式のデータを読み込むにはread_csv()関数を使用する。()内には読み込むファイルのパスを"で囲んで記入する。read_csv()関数はファイルの読み込みのみの機能しか持たない。現在の作業環境内に読み込んだデータを格納するためには代入演算子<-を使う。ここではdataフォルダー内のrct_data.csvを読み込み²、raw_dfという名のオブジェクトとしてく格納する。作業環境内のオブジェクトはRを再起動すると削除されるため、改めてパッケージ/データの読み込みが必要だ。

raw_df <- read_csv("data/rct_data.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 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()関数を使う。

dim(raw_df)

[1] 5000    8

　表形式データの場合、各列には名前が付いており、それぞれが一つの変数に該当する。これら変数名のみの出力にはnames()関数を使う。今回のデータだと、列の数が少ないこともあり、一画面に全列が表示されるが、数百列のデータとなると画面に収まらないので、変数名を確認しておくことを推奨する。

names(raw_df)

[1] "treatment" "gender"    "yob"       "hh_size"   "voted2000" "voted2002"
[7] "voted2004" "voted2006"

データハンドリングとパイプ演算子

　パイプ演算子には{magrittr}パッケージが提供する%>%とR 4.1から提供されるネイティブパイプ演算子の|>がある。現在の主流は古くから使われてきた%>%であるが、今後、|>が主流になると考えられるため、本講義では|>を使用する。しかし、多くの場合、|>の代わりに%>%を使っても同じ結果が得られる。

　パイプ演算子はパイプ前のオブジェクトを、パイプ後の関数の第一引数として渡す単純な演算子だ。たとえば、列名を変更する関数は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 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だが、パイプ演算子を使うと以下のようになり、人間にとって読みやすいコードになる。

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 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()を使用する。書き方は以下の通りだ。

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 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引数もデータであるため、以下のようにパイプ演算子を使うこともできる。

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 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を出力してみても、これまでのデータ加工内容は反映されていないことが分かる。

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章を参照されたい。

df <- raw_df |>
  rename(female = gender) |>
  mutate(female = if_else(female == "female", 1, 0),
         # 第1引数: votedで始まる変数を対象に処理を行う
         # 第2引数: 当該変数の値が"yes"なら1、それ以外なら0を割り当てる無名関数
         #          無名関数は「~」で始まり、変数が入る箇所は.xと表記する
         #          引数が当該変数のみであれば、「~」を付けずに関数のみでもOK
         across(starts_with("voted"), ~if_else(.x == "yes", 1, 0)))

記述統計量

　記述統計量の計算には{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）であるため、表示されない³。各グループがサンプルの何割かを計算するためには、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
   <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"で始まる列のみを抽出つして確認してみよう。

df |>
  select(starts_with("treatment"))

# 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を追加する⁴。

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]	1. Civic Duty 2. Control 3. Hawthorne 4. Neighbors 5. Self	578 ( 11.6% ) 2777 ( 55.5% ) 572 ( 11.4% ) 535 ( 10.7% ) 538 ( 10.8% )		5000 (100.0%)	0 (0.0%)
2	female [numeric]	Min : 0 Mean : 0.5 Max : 1	0 : 2494 ( 49.9% ) 1 : 2506 ( 50.1% )		5000 (100.0%)	0 (0.0%)
3	yob [numeric]	Mean (sd) : 1955.8 (14.4) min ≤ med ≤ max: 1910 ≤ 1956 ≤ 1986 IQR (CV) : 18 (0)	76 distinct values		5000 (100.0%)	0 (0.0%)
4	hh_size [numeric]	Mean (sd) : 2.2 (0.8) min ≤ med ≤ max: 1 ≤ 2 ≤ 6 IQR (CV) : 0 (0.4)	1 : 699 ( 14.0% ) 2 : 3112 ( 62.2% ) 3 : 831 ( 16.6% ) 4 : 308 ( 6.2% ) 5 : 45 ( 0.9% ) 6 : 5 ( 0.1% )		5000 (100.0%)	0 (0.0%)
5	voted2000 [numeric]	Min : 0 Mean : 0.3 Max : 1	0 : 3722 ( 74.4% ) 1 : 1278 ( 25.6% )		5000 (100.0%)	0 (0.0%)
6	voted2002 [numeric]	Min : 0 Mean : 0.4 Max : 1	0 : 3058 ( 61.2% ) 1 : 1942 ( 38.8% )		5000 (100.0%)	0 (0.0%)
7	voted2004 [numeric]	Min : 0 Mean : 0.4 Max : 1	0 : 2941 ( 58.8% ) 1 : 2059 ( 41.2% )		5000 (100.0%)	0 (0.0%)
8	voted2006 [numeric]	Min : 0 Mean : 0.3 Max : 1	0 : 3373 ( 67.5% ) 1 : 1627 ( 32.5% )		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）書き換えることもできる。

blc_chk <- df |>
  BalanceR(group = treatment,
           cov   = female:voted2004)

blc_chk

　標準化差分（標準化バイアス）を用いたバランスチェックはそれぞれのペアごとに計算を行うため、グループが多い場合は凡例が圧迫される場合が多い。しかし、重要なのは標準化差分の最大値だろう。ペア1、2、3でバランスが取れても、ペア4のバランスが取られていない場合は無意味だからだ。また、標準化差分の場合、符号の意味はなく、絶対値が重要だ。また、バランスチェックにおいてグループごとの平均値や標準偏差は不要である。ここでsummary()関数を使うと、絶対値が最も大きい標準化差分のみ出力される。

summary(blc_chk)

  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(blc_chk)

図 1: ?(caption)

　先ほど述べたようにバランスチェックで重要なのは絶対値が最も大きい標準化差分である。plot()内にsimplify = TRUEを指定すれば最大値のみ表示され、更にabs = TRUEにすると絶対値へ変換される。また、垂直のガイドラインはvline引数で変更できる。

# plot() の第1引数は blc_chk なのでパイプの使える
blc_chk |>
  plot(vline = c(5, 10), simplify = TRUE, abs = TRUE)

図 2: ?(caption)

処置効果の確認

グループごとの応答変数の平均値

　処置効果を確認するためには各グループごとの応答変数（ここではvoted2006）の平均値を計算し、処置群の平均値から統制群の平均値を引く必要がある。まずは、特定の変数の平均値を計算する方法について紹介する。データ内にある特定の変数の平均値を計算するためにはsummarise()関数内に平均値を求めるmean()関数を入れる。たとえば、dfのvoted2006の平均値を計算するコードは以下の通りである。

df |>
  summarise(mean(voted2006, na.rm = TRUE))

# A tibble: 1 × 1
  `mean(voted2006, na.rm = TRUE)`
                            <dbl>
1                           0.325

　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.325

　我々が知りたいのは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.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型は「順序付きの文字型変数」だと理解しても良い⁵。列の追加・上書き（今回は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を上書きしよう。

df <- df |>
  mutate(treatment = factor(treatment,
                            levels = c("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.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")

図 3: ?(caption)

　続いて、このグラフの見た目を調整してみよう。

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)

図 4: ?(caption)

　また、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)

図 5: ?(caption)

　小数点が長すぎるので3桁まで表示としよう。ここではsprintf()を使用する。使い方が簡単とは言えないが、覚える必要はなく、必要な時にググるか、本資料のコードをコピペすれば良い⁶。

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)

図 6: ?(caption)

　これで可視化ができた。ただし、以上のコードには改善の余地がある。geom_bar()とgeom_label()内のaes()関数に注目して欲しい。よく見るとxとyと同じだろう。geom_*()が共有するマッピングがあれば、ggplot()内で指定することでコードを効率化することもできる。

out_mean_df |>
  ggplot(aes(x = Groups, y = Outcome)) +
  geom_bar(stat = "identity") +
  geom_label(aes(label = sprintf("%.3f", Outcome))) +
  labs(y = "Mean(Outcome)") +
  theme_gray(base_size = 14)

図 7: ?(caption)

統計的推定（単回帰分析）

　これまでの作業はグループごとの応答変数の平均値であって、処置効果ではない。処置効果を計算するためには処置群の平均値から統制群の平均値を引く必要がある。たとえば、Civic Dutyはがき群の平均値は約0.317、統制群のそれは0.299であるため、Civic Dutyはがきの処置効果は約0.018である。しかし、これを各グループごとに計算することは面倒だし、何よりも得られた値が点推定値だという限界がある。得られた処置効果の不確実性は計算できない。

　ここで有効なのが線形回帰分析である。回帰分析を行うことで処置効果の点推定値のみならず、不確実性の指標である標準誤差も計算され、区間推定や統計的仮説検定も可能となる。線形回帰分析の関数は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.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）を指定する必要がある。パイプ前のオブジェクトが入る位置を任意に指定することであり、_を使う。%>%演算子を使う場合は_でなく、.を使う。上記のコードと以下のコードは同じコードとなる。プレースホルダーは自分が使うパイプ演算子によって使い分けること。

fit1 <- df |> # |> パイプを使う場合
  lm(voted2006 ~ treatment, data = _)

fit1 <- df %>% # %>% パイプを使う場合
  lm(voted2006 ~ treatment, data = .)

　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}で使うことはできない。

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.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

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.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))

図 8: ?(caption)

　それでは図をカスタマイズしてみよう。図内の様々なラベルを修正する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)

図 9: ?(caption)

　まだ気になる点がある。それは横軸の目盛りラベルに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 <- fit1_coef |>
  mutate(term = str_replace(term, "treatment", ""))

　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)

図 10: ?(caption)

　最後に横軸の順番を修正してみよう。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.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)

図 11: ?(caption)

　これで処置効果の可視化もバッチリだ。

多重比較の問題

　グループが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と修正する。

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.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と指定すれば良い。

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.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 と同じコードで作図する。

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 (w/ 98.75% CI)") +
  theme_bw(base_size = 12)

図 12: ?(caption)

統計的推定（重回帰分析）

　今回の例は無作為割当が成功しており、処置前変数の偏りは見られない。しかし、何らかの理由で処置前変数の偏りが生じる場合がある。その「何らかの理由」が応答変数にまで影響を与えるのであれば、それは交絡変数（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()関数を使えば、推定結果がより見やすくなる。

modelsummary(fit2)

	(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()関数でまとめると、モデル間比較もできる。

modelsummary(list("w/o Covariates" = fit1, "w/ Covariates" = fit2))

	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"を追加する。

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.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 <- bind_rows(list("Model 1" = fit1_coef, 
                           "Model 2" = fit2_coef),
                      .id = "Model")

　それでは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)

図 13: ?(caption)

　今回の結果だとモデル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)

図 14: ?(caption)

　何かおかしい。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)

図 15: ?(caption)

　これで図は完成だが、少し修正してみよう。{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")

図 16: ?(caption)

脚注

1. パッケージのアップデートにはpacman::p_update()またはpacman::p_up()関数を使う。()内に何も入力しない場合、全パッケージがアップデートされる。

2. 使用するデータはGerber, Green, and Larimer (2010)のデータである。実際のサンプルサイズは約34万であるが、ここでは無作為に5000個を抽出したサンプルを用いる。抽出前のデータはYale ISPSから入手可能である。

3. 変数のデータ型はデータを出力する際、列名の下段に表示される。<chr>は文字型、<dbl>と<int>は数値型、<fct>はfactor型である。他にもいくつかのデータ型がある。詳細は『私たちのR』の第8章を参照すること。

4. RMarkdown内に埋め込むなら更にstyle = "rmarkdown"を追加してみよう。ただし、Chunkオプションにresults = "asis"（Quartoなら#| results: "asis"）を付けること。

5. 中身が1、2、3、…であってもfactor型であれば1、2、3、…は数字でなく文字として認識される。

6. もっと使いやすいround()があるが、round()の場合、丸めた結果が1.100なら1.1としか表記されない。表示される桁数を固定するためにはsprintf()を使う。