社会科学における因果推論

5/ Lab session

宋財泫（関西大学）

1 本日の内容

本日の内容

  1. Rの導入（省略）
  2. データハンドリング
  3. 可視化
  4. 記述統計量の計算
  5. 平均処置効果の推定
    • 線形回帰分析
    • 多重比較
  6. 推定結果の可視化

2 データ操作

データ操作の手順

  1. 必要なパッケージの読み込み
  2. データの読み込み
  3. 列の整理
  4. 欠損値処理
  5. 変数の計算
  6. 列の整理

パイプ演算子

パイプ演算子（|>）：オブジェクトを関数の第一引数として渡す演算子

  • func(x, y)はx |> func(y)と同じ
    • パイプなし：print(mean(x, na.rm = TRUE))
    • パイプあり：x |> mean(na.rm = TRUE) |> print()
  • コーディングを考える順番（or コードを読む順番）と書く順番が一致
    • \(\therefore\) コーディングミスが少なく、メンテナンスしやすい

パイプ演算子の例（1）

x <- c(2, NA, 5, 7, 11)
print(mean(x, na.rm = TRUE))
## [1] 6.25

パイプ演算子の例（2）

x <- c(2, NA, 5, 7, 11)
# 必須ではないが、パイプ演算子の後は改行が推奨される。
x |> 
  mean(na.rm = TRUE) |> 
  print()
## [1] 6.25

列の整理

  • 列の抽出：select()
    • 列名の変更：rename()
    • 列の位置変更：relocate()
  • 行の抽出：filter()

欠損値の処理

欠損値（missing value）：本来あるべき値が記録されていないこと

  • 調査の無回答、観測機器の不具合、入力ミス、政治的な理由などで発生する
  • 分析に使用する変数が1つ以上欠損した行（ケース）は分析から除外される1ため、予め除外しておく2

欠損値処理の手順

同じデータセットでも欠損値の表記方法は複数あり得るため、自分が生産したデータでない場合はコードブックを参照すること

  1. Step 1：999や-9などの数値、""（空欄）が欠損値を意味する場合、これらの値をRの欠損値（NA）に置換
  2. Step 2：drop_na()やfilter()で欠損値を含む行を除去

値の置換

特定の値をNAに置換するreplace_with_na()関数（要{naniar}パッケージ）

  • 例）ageは999、educは-9と-10が欠損値（femaleは欠損値が既にNAになっている）
dfの中身
id female age educ
1 0 18 4
2 1 55 3
3 1 24 -9
4 1 999 -10
5 0 24 5
6 NA 39 4 
library(naniar)
df |> 
  replace_with_na(list(age  = 999,
                       educ = c(-9, -10)))
## # A tibble: 6 × 4
##      id female   age  educ
##   <int>  <dbl> <dbl> <dbl>
## 1     1      0    18     4
## 2     2      1    55     3
## 3     3      1    24    NA
## 4     4      1    NA    NA
## 5     5      0    24     5
## 6     6     NA    39     4

欠損値を含む行の削除

  1. filter()関数
    • a
  2. drop_na()関数
    • a

3 可視化

{ggplot2}の考え方

4 記述統計

記述統計量

必ず分析に使用するデータの記述統計量を報告すること

  • 分析に使われるケース（個体）、変数のみで十分
  • 分析に使われないケース（個体）、クリーニング前の変数は載せない。
  • よく使う記述統計量は平均値、標準偏差、最小値、最大値
    • 中央値も入れるとなおさら良い
  • 実験データの場合、グループごとに掲載する場合もあるが、別途バランスチェックを行うなら不要

記述統計表の例

変数名 平均値 中央値 標準偏差 最小値 最大値
女性 0.503 1 0.500 0 1
年齢 47.340 47 15.628 18 75
投票有無
　投票 0.736 1 0.441 0 1
　棄権 0.229 0 0.420 0 1
　参政権なし 0.035 0 0.184 0 1
感情温度
　自民党 41.130 50 28.015 0 100
　立憲民主党 34.248 40 25.947 0 100

記述統計量の計算方法

  • 方法1：R内蔵関数を利用する
    • mean()、median()、sd()、min()、max()
    • {dplyr}のsummarise()との組み合わせが効率的
    • 複数の変数の記述統計量を計算する場合、更にacross()と組み合わせるか、別途のパッケージ（{summarytools}など）を使った方が効率的
  • 方法2：{summarytools}のdescr()を利用する（推奨）

バランスチェック

{cobalt}パッケージのbal.tab()を利用した標準化平均差の計算

## Balance Measures
##                Type Diff.Un     M.Threshold.Un
## age         Contin. -0.2419 Not Balanced, >0.1
## educ        Contin.  0.0448     Balanced, <0.1
## race_black   Binary  0.6404 Not Balanced, >0.1
## race_hispan  Binary -0.0827     Balanced, <0.1
## race_white   Binary -0.5577 Not Balanced, >0.1
## married      Binary -0.3236 Not Balanced, >0.1
## re74        Contin. -0.5958 Not Balanced, >0.1
## re75        Contin. -0.2870 Not Balanced, >0.1
## 
## Balance tally for mean differences
##                    count
## Balanced, <0.1         2
## Not Balanced, >0.1     6
## 
## Variable with the greatest mean difference
##    Variable Diff.Un     M.Threshold.Un
##  race_black  0.6404 Not Balanced, >0.1
## 
## Sample sizes
##     Control Treated
## All     429     185

バランスチェック（可視化）

bal.tab() + love.plot()を利用した可視化

4 推定

t検定

分散分析

線形回帰モデル

lm(formula, data, weights)
  • formula：結果変数 ~ 説明変数1 + 説明変数2 + ...
  • data：結果変数および説明変数が格納されているデータフレーム名
    • パイプ演算子を使用する場合、データフレーム名 |> lm(formula, data = _)
  • weights：重み付け回帰分析の場合、重み変数の列名

多重比較の補正

処置効果

lm(結果変数 ~ 処置変数 + 統制変数1 + 統制変数2 + ..., data = データフレーム名)

予測値

predictions()

限界効果

slopes()

5 可視化

{ggplot2}による作図

棒グラフ

  • 幾何オブジェクト：geom_bar()、またはgeom_col()
  • マッピング
    • 必須：x、y
    • その他：fillなど

散布図

ヒストグラム

点—範囲グラフ

  • 幾何オブジェクト：geom_pointrange()
  • マッピング
    • 必須：x、y、ymin、ymax（または、x、xmin、xmax、y）
    • その他：color、alpha、shpaeなど