13 実験的手法

佐久間智広

sakuma@port.kobe-u.ac.jp

神戸大学大学院経営学研究科

2024/07/08

0.1 準備

0.1.1 パッケージの読み込み

Code

1if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse,magrittr,ggpubr,
               estimatr,car,modelsummary,ggrepel,
               patchwork,broom,psych,datarium,
2               sjPlot,dae)

1: 複数のパッケージを一度に読み込めるpacmanパッケージが入っていない場合は，インストールする
2: パッケージの読み込み

pacman: 複数のパッケージを読み込めるパッケージ。そのパッケージの中のp_loadコマンドを使うと，（）内に指示したパッケージを読み込める。さらにそのパッケージがそもそもインストールされていない場合にはインストールした後に読み込んでくれる
tidyverse: Rのコードを楽にするたくさんのパッケージをパッケージにしたもの
magrittr: tidyverseに含まれる%>%などのコードをさらに拡張するもの。これがあると%$%とか%<>%が使える

1 はじめに

ここまでしてきた話は
- 分析の前提となる確率・統計の話
  - 確率→統計的仮設検定
- 基礎的な分析の話
  - 平均値の差の検定（t検定）
  - 回帰分析
今日は，具体的な調査方法として，最もシンプルで，なおかつ（因果関係を明らかにすると言う意味では）一番厳密な実験を扱います。

2 実験的手法とは？

調査者が何らかの操作をして行う調査・研究

2.1 例1：治験

参加者をランダムに2つ（以上）のグループに割り振って，片方には新しい薬を，片方には薬ではないもの（プラセボ）を投与して，新薬の効果を検証する

2.2 例2：肥料の実験

同じ場所に並べた複数の植木のうち，ランダムに選んだいくつかには肥料を与え，残りには与えないことで，肥料の効果を検証する

2.3 例3：ABテスト

アクセスしてくる人のうち一定割合に別のデザインの画面を表示し，クリックする先やクリック率等を比較する。（どちらのデザインの方が良いか？）

出所：https://unique1.co.jp/column/ads/4324/

2.4 例4：従業員の業績評価方法

企業の人事制度の違いが効果をもたらすのか明らかにするため，ランダムに選択された部署のみ業績評価の方法を変えてみる。評価の方法を変えなかった部署と行動の変化を比較する
企業トップのメッセージが現場の従業員の行動に与える影響を見るため，ランダムに選ばれた店舗だけにトップのメッセージ動画を定期的に配信する。メッセージを受け取った店舗と受け取らなかった店舗の行動の違いを比較する (Cronin, Erkens, Schloetzer, and Tinsley 2021)

2.5 実験的手法の特徴

実験的手法は，調査者が何らかの介入をして，その介入効果を見るようなもの。そのような特徴から

因果関係を検証できる
実験で検証したい要因以外が調整されているので，堅い証拠が得られる（内的妥当性が高い）

元々理系では実験的手法が中心だったけれど，上記の特徴から社会科学領域でも実験的研究が重視されてきている。

2019年のノーベル経済学賞をとったBanerjee, Duflo, Kremerも「世界の貧困を改善するための実験的アプローチに関する功績」

一方で，

実験環境以外でも同じような結果が得られるかはわからない（外的妥当性が低い）
実験で検証できる要因（≒調査者が操作できる要因）しか実験できない
- 「職場における仲の良さが，仕事に対するモチベーションに与える影響」みたいなのは研究しづらい。
- 「株式市場の値動き」みたいなのも研究しづらい

2.6 実験研究と因果推論

実験研究と関連して，統計的手法を用いた因果関係の推論（統計的因果推論）が注目されつつある
- 2021年のノーベル経済学賞をとったCard, Angrist, Imbensの受賞理由は「因果関係の分析への方法論的貢献」
回帰分析は，yを結果，xを原因とみなした因果関係を検証しているように見える（し，多くの場合因果関係を暗に想定している）が，実際に算出しているのは相関関係
- 逆因果の可能性も否定できない。
実験的な調査方法と組み合わせることで，回帰分析の結果は因果関係と見なせる
- 他の調査方法で因果関係を証明するためには，統計的に行動な技術が必要
  - 過激派は実験以外では因果関係は証明不可能と主張

3 実験におけるデータ分析例

3.1 データ

Rにサンプルデータとして入っているheadacheデータを使用

ある製薬会社が片頭痛患者を対象に3種類の治療法を試験した。実験には72人の参加者が登録された。その目的は、片頭痛エピソードに関連する痛みのスコアを下げる新しいクラスの治療法の可能性を調べることである。

参加者は男性36名、女性36名である。男性と女性はさらに（等しく）片頭痛のリスクが低いか高いかに細分化された。

以下の変数を含む：

gender性別：「男性」と「女性」；

riskリスク：低 “と”高”

treatment治療の種類：3つのカテゴリーがある： X”、“Y”、“Z”の3つのカテゴリーがある。

Code

data("headache")
headache %<>% 
  filter(treatment != "Y") %>% 
  mutate(treatment = as.character(treatment),
         treatment = as.factor(treatment))

分析例を簡単にするため，treatment Yを落とした。なので，Xを治療なし，Zを薬による治療と解釈して話を進める

3.2 記述統計量

記述統計を見ると，gender(男女)とrisk(高い低い)は各条件に均等に配分されている。

Code

datasummary(pain_score + gender + risk ~ 
              treatment * (mean + sd+N + Percent()),
            data = headache)

tinytable_pwnpcvx4q8bivqd00bfq

		X				Z
		mean	sd	N	Percent	mean	sd	N	Percent
	pain_score	80.45	8.57	24	50.00	76.24	5.88	24	50.00
gender	male			12	25.00			12	25.00
	female			12	25.00			12	25.00
risk	high			12	25.00			12	25.00
	low			12	25.00			12	25.00

Code

ggplot(headache, aes(x = treatment, y = pain_score))+
  geom_boxplot()

3.3 分析方法(分散分析)

主な分析は分散分析（ANOVA）。aov(結果 ~ 要因, data = データ名)コマンドで実行可能

Code

aov(pain_score ~ treatment, data = headache) |>  
  summary()

            Df Sum Sq Mean Sq F value Pr(>F)  
treatment    1  213.2  213.21   3.946  0.053 .
Residuals   46 2485.7   54.04                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

仮説検定での帰無仮説は治療の効果はない。治療の違いの効果は10%水準で有意（水準を10%としたら，関係ないと言う帰無仮説は棄却される）

ANOVAの原理はダミー変数を使った回帰分析と一緒。なので以下のような単回帰分析をしても同じ

Code

lm(pain_score ~ treatment, data = headache) |> 
  summary()


Call:
lm(formula = pain_score ~ treatment, data = headache)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.093  -5.620  -1.165   4.403  19.547 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   80.453      1.501  53.617   <2e-16 ***
treatmentZ    -4.215      2.122  -1.986    0.053 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.351 on 46 degrees of freedom
Multiple R-squared:  0.079, Adjusted R-squared:  0.05898 
F-statistic: 3.946 on 1 and 46 DF,  p-value: 0.05297

更に，性別やそもそもの頭痛リスクを踏まえた分析をしてみる

3.3.1 リスクを含めた分析

Code

aov(pain_score ~ treatment * risk, data = headache) |> 
  summary()

               Df Sum Sq Mean Sq F value   Pr(>F)    
treatment       1  213.2   213.2   6.741   0.0128 *  
risk            1 1076.0  1076.0  34.019 5.94e-07 ***
treatment:risk  1   18.1    18.1   0.573   0.4531    
Residuals      44 1391.7    31.6                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

リスク高低で分けると，治療の効果は5％水準で有意に。
リスクも有意。そもそもリスクが高い人は頭痛度合いが強い。
効果は似た感じ（交互作用はない）。

また，ggline()コマンドを使うと，効果を図に表せる

Code

ggline(data = headache, #データを指定
       x = "treatment", #x軸
       y = "pain_score", #y軸
       color = "risk", #色分けしたい条件
       add = c("mean") #それぞれの条件での平均値を載せる
       ) +
  theme_minimal()

3.3.2 性別を含めた分析

男女で効き目が違うか

Code

aov(pain_score ~ treatment * gender, data = headache) |>  
  summary()

                 Df Sum Sq Mean Sq F value Pr(>F)  
treatment         1  213.2  213.21   4.475 0.0401 *
gender            1  273.4  273.36   5.738 0.0209 *
treatment:gender  1  116.2  116.20   2.439 0.1255  
Residuals        44 2096.2   47.64                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

headache %$% 
  interaction.plot(treatment, gender, pain_score, 
                   type = "b",  pch=c(1,2))

男性の方が効果が大きい。ただし，交互作用が有意とまではいかない（p < 0.123）

3.3.3 3要因を含んだ分析

Code

aov(pain_score ~ treatment+risk+gender, 
    data = headache) |>  
  summary()

            Df Sum Sq Mean Sq F value   Pr(>F)    
treatment    1  213.2   213.2   8.255  0.00623 ** 
risk         1 1076.0  1076.0  41.659 7.24e-08 ***
gender       1  273.4   273.4  10.584  0.00219 ** 
Residuals   44 1136.4    25.8                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

treatment, gender, riskはそれぞれ頭痛の度合いと有意に関係

Code

aov(pain_score ~ treatment*risk*gender, 
    data = headache) |> 
  summary()

                      Df Sum Sq Mean Sq F value   Pr(>F)    
treatment              1  213.2   213.2  11.646 0.001485 ** 
risk                   1 1076.0  1076.0  58.770 2.25e-09 ***
gender                 1  273.4   273.4  14.931 0.000399 ***
treatment:risk         1   18.1    18.1   0.990 0.325773    
treatment:gender       1  116.2   116.2   6.347 0.015854 *  
risk:gender            1   26.5    26.5   1.449 0.235751    
treatment:risk:gender  1  243.2   243.2  13.286 0.000762 ***
Residuals             40  732.3    18.3                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

treatmentとriskとgenderの3要因に交互作用

これがどう言う意味かというと…次のページの図から，「男性で，尚且つリスクが高い人」にとって，治療Zは効果がある

Code

interaction.ABC.plot(pain_score, 
                     x.factor = treatment,
                     groups.factor = gender, 
                     trace.factor = risk,
                     data = headache,
                     title = "Effect")

回帰分析の結果を並べるとこんな感じ

Code

list(
  m1 <- lm(pain_score ~ treatment ,data = headache),
  m2 <- lm(pain_score ~ treatment + gender,data = headache),
  m3 <- lm(pain_score ~ treatment * gender,data = headache),
  m4 <- lm(pain_score ~ treatment + risk,data = headache),
  m5 <- lm(pain_score ~ treatment * risk,data = headache),
  m6 <- lm(pain_score ~ treatment + gender + risk ,data = headache),
  m7 <- lm(pain_score ~ treatment * gender * risk ,data = headache)
  ) %>% 
  msummary(.,stars = TRUE,
             gof_omit = "Log.Lik.|AIC|BIC|RMSE|F",
             statistic = NULL
           )

tinytable_ksjlt0u6u5ms5hdhu69t

	(1)	(2)	(3)	(4)	(5)	(6)	(7)
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
(Intercept)	80.453***	82.839***	84.395***	85.188***	85.802***	87.574***	92.739***
treatmentZ	-4.215+	-4.215*	-7.327*	-4.215*	-5.444*	-4.215**	-13.058***
genderfemale		-4.773*	-7.885**			-4.773**	-13.874***
treatmentZ × genderfemale			6.224				15.228***
risklow				-9.469***	-10.698***	-9.469***	-16.687***
treatmentZ × risklow					2.458		11.462**
genderfemale × risklow							11.978**
treatmentZ × genderfemale × risklow							-18.009***
Num.Obs.	48	48	48	48	48	48	48
R2	0.079	0.180	0.223	0.478	0.484	0.579	0.729
R2 Adj.	0.059	0.144	0.170	0.454	0.449	0.550	0.681

4 まとめ

実験研究の特徴と具体的な分散方法について扱いました。
実験は，厳密な因果関係の推論ができる反面，社会調査において実験できる機会や状況は限られている
実験では主に分散分析が使われるけれど，本質は回帰分析と同じ

5 課題

ハンバーガー統計学より引用

ある小学校で、算数の分数の計算を教えるためのマンガを使った新しい教材を開発した。この教材の効果は従来のものと比べ、効果があるようだということはわかっているが、さらにその効果が子どもの算数に対する好みによってどう違うのかを調べたいと思う。

この効果を調べるために、あるクラスでは従来の教材で教え（統制群）、別のクラスではマンガ教材で教えた（マンガ群）。一日おいて、分数の計算テストをした。このとき、各クラスで、算数が好きか嫌いかというアンケートをあらかじめ取っておき、算数が好きな子ども10人と嫌いな子ども10人とで比較することにした。

テストの点数データ（10点満点）は次のようになった。これを分散分析したい。

	統制群	マンガ群
算数が好き	7, 8, 6, 8, 10, 7, 8, 8, 9, 7	8, 9, 10, 10, 8, 8, 9, 7, 10, 8
算数が嫌い	4, 6, 5, 4, 3, 7, 5, 6, 4, 5	8, 7, 8, 6, 9, 7, 8, 8, 10, 8

(1)　この検定での帰無仮説を言いなさい。

(2)　この検定での対立仮説を言いなさい。

(3)　4つの条件におけるそれぞれの平均と標準偏差を計算しなさい（コードと結果）。

(4)　分散分析を行いなさい（コード）。

(5)　有意水準を1%としたとき、この分散分析表から言えることを書きなさい。

(6)　以上の検定の結果を、わかりやすいことばで説明しなさい。

次ページにデータ作成用コマンドがあります。

データ作成

Code

score <- c(7, 8, 6, 8, 10, 7, 8, 8, 9, 7,
           8, 9, 10, 10, 8, 8, 9, 7, 10, 8,
           4, 6, 5, 4, 3, 7, 5, 6, 4, 5,
1           8, 7, 8, 6, 9, 7, 8, 8, 10, 8)
like <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
          1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
2          0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
manga <-c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
          1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
          0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
3          1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
4kadai <- tibble(score,like,manga) %>%
  mutate(like = factor(like, levels = c(1, 0), labels = c("like", "dislike")),
         manga = factor(manga, levels = c(1, 0), labels = c("manga", "normal")))

1: 表にある得点を入力
2: 算数が好きなら1，そうでないなら0
3: 漫画が好きなら1
4: 上で作った３つのベクトルをまとめて一つのデータ (tibble)に。続けて0,1変数をダミー変数と認識させて（factor），ラベルをつけている

6 参考文献

Cronin, M., D. H. Erkens, J. D. Schloetzer, and C. H. Tinsley. 2021. How Controlling Failure Perceptions Affects Performance: Evidence from a Field Experiment. The Accounting Review 96 (2): 205–230.