政策分析の量的基礎(入門)：補足 実験データの分析
授業計画
※太字はこのレジュメと関係する回
- 授業概要の説明
2. リサーチデザイン・データ収集の方法
- Rのインストール・記述統計の出し方
4. 二変数の関係（クロス表・散布図）
5. 統計的検定
回帰分析（1）
回帰分析（2）
回帰分析（3）
パネルデータ分析
計量分析を用いた文献を読む
最終発表に関する実習・相談
学生による分析結果発表（1）
学生による分析結果発表（2）
学生による分析結果発表（3）
フィードバック
I. 実験の設計とデータの収集
1. 実験データで測定するもの
潜在的帰結（potential outcome）にもとづく因果推論の立場に立って、処置を与える処置群（介入群）、与えない対照群（統制群）との間で結果の違いを比較する（ランダム化比較試験（Randomized Controlled Trial: RCT））。
⇒結果において生じた差を、処置や介入の平均因果効果（Average Treatment Effetc: ATE）とみなす。
※結果のの違いについての効果量は、いつもATEとして求められるとは限らない。例えば、他に…
局所因果効果(Local Average Treatment Effect: LATE)
平均限界効果（Average Marginal Effect: AME）
平均要素限界効果（Average Component Marginal Effect: ACME）
処置群における平均処置効果（Average Treatment effect on the Treated: ATT）\(\leftrightarrow\) 統制群における平均処置効果（Average Treatment effect on the Untreated: ATU）
このように多様な平均効果が存在し、リサーチ・デザイン、分析手法によって厳密に分けられるので注意しよう。
2. 実験の設計、何を使って実験のための調査票を作るのか？
ここでは、研究設計の実例を示す。
研究の設計：企業の社会的責任（corporate social responsibility: CSR）に関する情報は投資行動に影響するか、さらに投資家の党派性は、そのメカニズムにどのように作用するかを検証
研究の背景: アメリカをはじめとする欧米諸国で問題となる党派的分極化のもと、投資行動をはじめとする経済行動に党派性が与える影響に関心が集まっている。アメリカでは、投資をめぐる意思決定において、民主党支持者の左派が企業の環境配慮といった社会的に責任に関する情報をもとに投資行動を決める場合があるのに対して、共和党支持者の右派にその傾向は認められないことが明らかになっている。こうした傾向は、他の民主主義体制においても発現しているのか。
研究の問い: 日本においても、左派の野党支持者がCSR情報をもとに投資行動を決め、右派の与党支持者にそのような傾向が認められないという対比があるか？党派性が投資行動を左右するのか？
仮説の設定：
仮説1：左派の野党支持者ほど、CSRにもとづいて投資行動を行う
仮説2：右派の与党支持者ほど、CSRにもとづいて投資選択を行わず、企業の業績に従って投資選択を行う
変数の設定
結果変数（outcome variables)：投資意欲／投資意欲金額
処置変数（treatment variables）：CSRに関する肯定的情報／CSRに関する否定的情報／企業業績に関する肯定的情報／企業業績に関する否定的情報＋無情報（統制群） ※処置の詳細については動画内の紹介を参照
調整変数(moderators): 党派性
共変量（covariates）：性別、年齢、教育歴、所得、投資額、NISA利用状況
実験をめぐるフロー
初期説明、最終意思確認、倫理審査、pre-regsitrationについての戦略・検討
初期説明：
- 政治経済に対する立ち入った質問を行うことへの断り
- 身体的拘束はないが不快を感じたら調査をやめれるように伝達
- 個人情報保護
最終意思確認・debriefing:
- 実験の処置についての解説を行う。虚偽情報、誇張情報、架空情報があったことを伝達
- 不快感をもたらす恐れがあった点について丁寧に説明
- 調査結果を奏しない選択があることを伝達
- 無断での情報の転用・転載は法的な処罰の対象となることを伝達
倫理審査(Institutional Review Board: IRB)： 実験計画を立案し調査票を作り終えたら、所属機関による倫理審査（IRB）を受けねばならない。京都大学法学部・法学研究科の場合、法政策共同研究センターにおいて、研究倫理審査が実施される。本調査に関わっての実際の倫理審査の文書はこちら。
実験の事前登録(pre-registration)：
3. Qualtricsを用いた調査票・画面の作成
今日、政治学においてRCTを行う際に、最も利用可能性が高い手法はオンライン上でのサーベイ実験（online survey experiment）である。以下では、上記の研究設計の実例を示しながら、オンライン上でのサーベイ実験を行う手順を概観し、得られたデータの分析方法について説明する。
オンライン上でのサーベイ実験において、もっと強力なツールとなっているのがQualtrics。京都大学法学部生であれば、機関契約のもと多様なツールを無料で使用できる（Qualtricsについての詳細な日本語マニュアルについては、五十嵐祐研究室HPを参照してほしい）。
さて、Qualtricsでアカウントを作り、プロジェクトを作った画面は以下の例となる（プロジェクト作成についての詳細は動画を参照）。示しているのは、調査の初期画面である。私もかかわった共同研究であるが、調査前のインフォームド・コンセント（調査の初期説明と承諾の受領プロセス）がいかに重要かわかってもらえるだろう。
調査プロジェクトは、次のように1つ1つのブロックから構成される。ブロックを構成しないとRCTができない。ブロックをもとに類似質問をまとめる習慣をつけておこう。
このブロックをもとに、処置群と統制群を設定していく。以下の例は、「株価行政に関する肯定的情報群」である（他群の概観は動画を参照）。
以下図から、ブロックを配することによって、ランダム化を設定していることが分かってもらえるだろう（ランダム化の詳細については動画を参照）。
II. 実験データを使った分析
1. 変数の設定とデータフレームへの追加
収集したデータに対して、分析のためにいくつもの加工が必要になる。以下のように、パイプ関数を使って1本のコードとして、効率的に変数設定を行ってみよう。
【主要変数】
党派性変数：仮説にもとづくと、党派性・政党支持が重要な変数となるので、党派性変数を作成
処置変数：閲覧した文章に従って、処置群と統制群を整理
結果変数：結果として、①株の購入意欲(intention)、②株の購入希望額(price)
まずは、分析で必要となるパッケージの導入。
#==========================================
# load library
#==========================================
library(ggeffects)
library(ggfortify)
library(ggpubr)
library(ggrepel)
library(ggsignif)
library(gridExtra)
library(MASS)
library(Rmisc)
library(scales)
library(stargazer)
library(stringr)
library(tidyverse)#==========================================
# setting variables
#==========================================
df_csr <- read.csv("csr.csv") %>%
mutate(
# —————————————— partisanship —————————————— #
psu_rul = case_when( # ruling party support
Q2.1 %in% c(1,4) ~ 1,
Q2.1 %in% c(13,14) ~ NA_real_,
TRUE ~ 0
),
psu_op = case_when( # opposite party support
Q2.1 %in% c(2,3,5:11) ~ 1,
Q2.1 %in% c(13,14) ~ NA_real_,
TRUE ~ 0
),
psu_indep = case_when( # independentn(無党派)
Q2.1 == 12 ~ 1,
Q2.1 %in% c(13,14) ~ NA_real_,
TRUE ~ 0
),
# —————————————— group setting —————————————— #
group = coalesce( #Treatment group setting
if_else(Q10.2 > 0, "Performance positive", NA_character_),
if_else(Q11.2 > 0, "Performance negative", NA_character_),
if_else(Q12.2 > 0, "CSR positive", NA_character_),
if_else(Q13.2 > 0, "CSR negative", NA_character_),
if_else(Q14.1_1 > 0, "Control", NA_character_)
),
# —————————————— group dummies —————————————— #
perf_pos = dplyr::if_else(!is.na(Q10.2) & Q10.2 > 0, 1L, 0L),
perf_neg = dplyr::if_else(!is.na(Q11.2) & Q11.2 > 0, 1L, 0L),
csr_pos = dplyr::if_else(!is.na(Q12.2) & Q12.2 > 0, 1L, 0L),
csr_neg = dplyr::if_else(!is.na(Q13.2) & Q13.2 > 0, 1L, 0L),
control_group = dplyr::if_else(!is.na(Q14.1_1) & Q14.1_1 > 0, 1L, 0L),
# —————————————— Outcome 1: stock‐purchase intention ——————————— #
intention = coalesce(Q10.2, Q11.2, Q12.2, Q13.2, Q14.1_1),
outcome_intention = case_when(
intention == 1 ~ 1,
intention == 2 ~ 0,
intention %in% c(3,4) ~ NA_real_,
TRUE ~ NA_real_
),
# —————————————— Outcome 2: rescale price —————————————— #
buy_pospr = scales::rescale(Q10.3_1, to = c(0,1), from = c(1800,2000)),
buy_negpr = scales::rescale(Q11.3_1, to = c(0,1), from = c( 300, 700)),
buy_poscsr = scales::rescale(Q12.3_1, to = c(0,1), from = c(1800,2000)),
buy_negcsr = scales::rescale(Q13.3_1, to = c(0,1), from = c(1800,2000)),
buy_control = scales::rescale(Q14.1_1, to = c(0,1), from = c(1000,1400)),
outcome_buy = coalesce(buy_pospr, buy_negpr, buy_poscsr, buy_negcsr, buy_control),
# —————————————— Outcome 3: cabinet approval —————————————— #
cabap = case_when( # ruling party support
Q15.1 %in% c(1) ~ 2,
Q15.1 %in% c(3) ~ 1,
Q15.1 %in% c(2) ~ 0,
Q15.1 %in% c(4) ~ NA_real_
),
# —————————————— Covariates —————————————— #
gender = case_when(
Q17.2 == 1 ~ 1, # male
Q17.2 == 2 ~ 0, # female
Q17.2 %in% c(3,4) ~ NA_real_
),
age = Q17.3, # age
education = case_when( # 1-3: no university degree 4: higher than undergraduate degree
Q17.4 %in% 1:4 ~ Q17.4,
Q17.4 == 5 ~ NA_real_,
TRUE ~ NA_real_
),
income = Q17.5_1, # income
investment = Q17.7_2 + Q17.7_4, # stock
savings = Q17.7_1,
nisa = case_when( # NISA dummy
Q17.9_1 == 1 ~ 1,
Q17.9_1 == 2 ~ 0,
Q17.9_1 == 3 ~ NA_real_,
TRUE ~ NA_real_
)
)2. 分散分析、平均の差のプロット(mean-difference plot)による分析
2群間の差であれば、ここまでも見てきた\(t\)検定を用いればよい。しかし、3群以上、すなわち多群間の検定に際しては、まず分散分析(Analysis of Variance: ANOVA)を用いる。
ANOVAにおける帰無仮説と対立仮説は以下の通りである。
\(H_0\):すべての群平均が等しい(\(\tau_1=\tau_2...=0\))
\(H_1\):少なくとも1群の平均が異なる
ANOVAでは、以下のモデルにおいて分散を群間変動(Between Sum of Squares; SSB)と群内変動（Within Sum of Squares: SSW）とに分解する。
\[ Y_ij = \mu + \tau_i+\epsilon_{ij}, \]
ここで、\(Y_{ij}\)はグループ\(_i\)の\(j\)番目の観測値、\(\mu\)は全体平均、\(\tau_{i}\)はグループ\(_{i}\)の効果（群平均―全体平均）, \(\epsilon_{ij}\)は誤差項であり、互いに独立で分散\(\sigma^2\)に従う。この設定もと、全変動(Total Sum of Squares; SST)を、
\[ SST=SSB+SSW \]
として分解して、自由度のもとに平均平方和を算出し、\(F\)検定を実施。群間平方和の平均平方\(MSB=\frac{SSB}{df_B}\)、群内平方和の平均平方が\(MSW=\frac{SSW}{df_W}\)とすると、
\[ F-statistics=\frac{MSB}{MSW} \]
となる。ここで「群間変動が群内変動に比べて大きい=平均差が有意」かを確かめる。
この結果をもとに、群間変動が群内変動に比して大きく、少なくとも1群の平均が他と異なることがわかれば、どこに差があるのかを個別のペアワイズのもとに検証していくことになる。そこで用いられるのが、Tuley-HSD検定。以下では、ANOVAからTuley-HSD検定までを行う。
#==========================================
# Mean difference and Tukey-HSD test
#==========================================
library(ggpubr)
library(multcompView)
# ------------------
# 1) summarise data
# ------------------
df_csrc <- df_csr %>%
dplyr::filter(!is.na(group), !is.na(outcome_intention))
df_summary <- df_csrc %>%
dplyr::group_by(group) %>%
dplyr::summarise(
mean_int = mean(outcome_intention),
sd = sd(outcome_intention),
n = length(outcome_intention),
se = sd / sqrt(n),
ci = qt(0.975, df = n - 1) * se,
.groups = "drop"
)
# --------------------------------
# 2) Analysis of Variance (ANOVA)
# --------------------------------
fit_aov <- aov(outcome_intention ~ group, data = df_csrc)
summary(fit_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## group 3 27.31 9.104 56.01 <2e-16 ***
## Residuals 898 145.98 0.163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ---------------------------------------------------
# 3) Tukey-HSD test: Between group pairwised test
# ---------------------------------------------------
tuk_df <- TukeyHSD(fit_aov, conf.level = 0.95)$group
tuk_df## diff lwr upr
## CSR positive-CSR negative 0.07299065 -0.02907325 0.17505455
## Performance negative-CSR negative -0.08902834 -0.18774357 0.00968689
## Performance positive-CSR negative 0.36112033 0.26185682 0.46038385
## Performance negative-CSR positive -0.16201899 -0.25893354 -0.06510445
## Performance positive-CSR positive 0.28812968 0.19065672 0.38560264
## Performance positive-Performance negative 0.45014867 0.35618785 0.54410949
## p adj
## CSR positive-CSR negative 2.549434e-01
## Performance negative-CSR negative 9.396315e-02
## Performance positive-CSR negative 2.619016e-13
## Performance negative-CSR positive 1.099703e-04
## Performance positive-CSR positive 6.703527e-13
## Performance positive-Performance negative 2.212674e-13まずANOVAから\(F=56.01, \quad Pr(>F)<0.0000\)との結果から、群間変動は群内変動よりも大きく、少なくとも1つ以上の群の平均が他と異なることが分かる。
そして、4群ペア間の結果の比較から、「CSR positive-CSR negative」間に有意差はないが、他の組み合わせでは有意差があることが分かる。このTukey-HSD検定の結果を反映させ、群間の平均値の差を描画してみよう。
# ----------------------------------------
# 3) Make the table of Tukey-HSD results
# ----------------------------------------
tuk_df <- as.data.frame(tuk_df) %>%
tibble::rownames_to_column("pair") %>%
tidyr::separate(pair, into = c("group1", "group2"), sep = "-") %>%
dplyr::rename(p.adj = `p adj`) %>%
mutate(
y.position = max(df_summary$mean_int + df_summary$ci) +
seq(0.02, 0.08, length.out = nrow(.)),
p.adj.signif = stats::symnum(
p.adj,
corr = FALSE,
na = FALSE,
cutpoints = c(0, .001, .01, .05, 1),
symbols = c("***", "**", "*", "ns")
)
)
# ------------------------------------------------------
# 4) Plot mean difference and reflect Tukey-HSD result
# ------------------------------------------------------
ggplot(df_summary, aes(x = group, y = mean_int)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = mean_int - ci, ymax = mean_int + ci),
width = 0.2) +
geom_label(aes(label = paste0(round(mean_int, 2), "\n[",
round(mean_int - ci, 2), ", ",
round(mean_int + ci, 2), "]")),
fill = "white",
label.r = unit(0.15, "lines"),
label.size = 0.25,
size = 3,
hjust = 0,
vjust = 1,
nudge_x = 0.1,
nudge_y = -0.05
) +
ggpubr::stat_pvalue_manual(
tuk_df,
label = "p.adj.signif",
xmin = "group1",
xmax = "group2",
y.position = "y.position",
tip.length = 0.01,
step.increase = 0.02
) +
scale_x_discrete(guide = guide_axis(angle = 45)) +
theme_bw() +
labs(
title = "Mean intention by group (mean[lower, upper])",
x = "Group",
y = "Intention"
)
次に党派性によって、株式の購入意欲が変化するかを確かめたいので、党派性によって結果をサブグループに分ける。ここでは、党派性を調整変数(moderator)とする分析を行っていることになる。調整変数という用語も理解しておくようにしよう。
#==========================================
# mean difference: partisanship*intention
#==========================================
df_summary_rul <- df_csr %>%
dplyr::filter(!is.na(group), !is.na(outcome_intention), !is.na(psu_rul)) %>%
mutate(
group = factor(group),
psu_rul = factor(psu_rul,
levels = c(0,1),
labels = c("Others","In‑partisans"))
) %>%
# summarise mean and sd by group
dplyr::group_by(group, psu_rul) %>%
dplyr::summarise(
mean_int = mean(outcome_intention, na.rm = TRUE),
sd = sd(outcome_intention, na.rm = TRUE),
n = length(outcome_intention), se = sd / sqrt(n),
ci = qt(0.975, df = n - 1) * se,
.groups = "drop"
)
p_rul <- ggplot(df_summary_rul,
aes(x = group, y = mean_int, shape = psu_rul)) +
geom_errorbar(aes(ymin = mean_int - ci, ymax = mean_int + ci),
width = 0.2,
position = position_dodge(width = 0.3),
color = "grey40") +
geom_point(position = position_dodge(width = 0.3),
size = 3,
color = "black") +
scale_shape_manual(values = c(1, 16)) +
scale_x_discrete(guide = guide_axis(angle = 45)) +
theme_bw(base_size = 14) +
theme(legend.position = "bottom",
legend.title = element_blank()) +
labs(title = "Mean Intention by Group\n
(In-partisans' effect)",
x = "Group",
y = "Mean Intention")
# ----------------------
# 2) Out-partisan case
# ----------------------
df_summary_op <- df_csr %>%
dplyr::filter(!is.na(group), !is.na(outcome_intention), !is.na(psu_op)) %>%
mutate(
group = factor(group),
psu_op = factor(psu_op,
levels = c(0,1),
labels = c("Others","Out-partisans"))
) %>%
dplyr::group_by(group, psu_op) %>%
dplyr::summarise(
mean_int = mean(outcome_intention, na.rm = TRUE),
sd = sd(outcome_intention, na.rm = TRUE),
n = length(outcome_intention),
se = sd / sqrt(n),
ci = qt(0.975, df = n - 1) * se,
.groups = "drop"
)
p_op <- ggplot(df_summary_op,
aes(x = group, y = mean_int, shape = psu_op)) +
geom_errorbar(aes(ymin = mean_int - ci, ymax = mean_int + ci),
width = 0.2,
position = position_dodge(width = 0.3),
color = "grey40") +
geom_point(position = position_dodge(width = 0.3),
size = 3,
color = "black") +
scale_shape_manual(values = c(1, 16)) +
scale_x_discrete(guide = guide_axis(angle = 45)) +
theme_bw(base_size = 14) +
theme(legend.position = "bottom",
legend.title = element_blank()) +
labs(title = "Mean Intention by Group\n
(Out-partisans' effect)",
x = "Group",
y = "Mean Intention")
# ----------------------
# 3) Arrange two plots
# ----------------------
grid.arrange(p_rul, p_op, ncol = 2)
WORK
- 結果変数をoutcome_intentionから、Outcome_buyに変えてみよう。
## Df Sum Sq Mean Sq F value Pr(>F)
## group 4 1.28 0.3206 5.044 5e-04 ***
## Residuals 1001 63.63 0.0636
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## diff lwr upr
## CSR negative-Control 0.040916416 -0.029879544 0.11171238
## CSR positive-Control 0.035795013 -0.035461474 0.10705150
## Performance negative-Control 0.099786061 0.030014414 0.16955771
## Performance positive-Control 0.084335884 0.016652160 0.15201961
## CSR positive-CSR negative -0.005121403 -0.075725339 0.06548253
## Performance negative-CSR negative 0.058869645 -0.010235431 0.12797472
## Performance positive-CSR negative 0.043419468 -0.023576915 0.11041585
## Performance negative-CSR positive 0.063991048 -0.005585747 0.13356784
## Performance positive-CSR positive 0.048540871 -0.018941972 0.11602371
## Performance positive-Performance negative -0.015450177 -0.081363226 0.05046287
## p adj
## CSR negative-Control 0.5109995480
## CSR positive-Control 0.6453233669
## Performance negative-Control 0.0009397778
## Performance positive-Control 0.0061763058
## CSR positive-CSR negative 0.9996558785
## Performance negative-CSR negative 0.1368665639
## Performance positive-CSR negative 0.3911415939
## Performance negative-CSR positive 0.0884476907
## Performance positive-CSR positive 0.2836684699
## Performance positive-Performance negative 0.9683706954
3. 線形確率モデル(Linear Probability Model: LPM)による分析
1つ目の結果変数、購入意欲は\({0,1}\)からなる二値変数である。OLSにおいて紹介したガウス・マルコフ定理において、二値変数の推定は、誤差の不均一分散性、線形性の侵をもたらす。そこで、二値変数に対して、OLS推定は一切許されないかというとそうではない。二値変数であったとしても、OLS推定を使うことにより、Xの増分に対するYの増分をとらえることができる。二値変数を従属変数とする線形モデルを、線形確率モデル(Linear Probability Model: LPM)と呼び、実験データの分析に際してはよく利用されるので覚えておこう。
仮説に即して推定したい1つ目のモデルは、処置効果を直接的に測るベースラインのモデルである；
\[ Y_i=\beta_0 +\beta_1D_i+X_i'\gamma+\epsilon_i \]
ここで、\(Y_i\)は被験者\(_i\)の購入意欲、\(D_i\)は処置ダミー、\(X_i\)は共変量ベクトルで、本分析では、年齢(\(age\))、性別(\(gender\))、教育歴(\(education\))、所得(\(income\))、株式などの投資額(\(investment\))、NISA利用の有無(\(nisa\))である。
次に条件付けの仮説に即して、調整変数として党派性を加えるモデルは以下の通りである；
\[ Y_i=\beta_0 +\beta_1D_i+\beta_2Z_i+\beta_3(D_i\times Z_i{})X_i'\gamma+\epsilon_i \]
2つのモデルに、データを当てはめた結果が以下の通りである。
#=================================================
# Average Treatment Effect with the OLS estimator
#=================================================
# ----------------------
# 1) Base-line model
# ----------------------
model_lpm <- df_csr %>%
dplyr::select(outcome_intention, group, perf_pos, perf_neg, csr_pos, csr_neg, control_group, psu_rul, psu_op, group, age, gender, education, income, investment, nisa) %>%
na.omit() %>%
with(lm(outcome_intention ~ group + age + gender + education + income + investment + nisa))
# --------------------------
# 2) Interaction term model
# --------------------------
model_lpm_int <- df_csr %>%
dplyr::select(outcome_intention, group, perf_pos, perf_neg, csr_pos, csr_neg, control_group, psu_rul, psu_op, group, age, gender, education, income, investment, nisa)%>%
na.omit() %>%
with(lm(outcome_intention ~ group + psu_rul + psu_op + group*psu_rul + group*psu_op + age + gender + education + income + investment + nisa))
stargazer(model_lpm, model_lpm_int,
type = "text",
title = "OLS estimation results (LPM)",
dep.var.labels = "Outcome: intention",
out.header = FALSE)##
## OLS estimation results (LPM)
## ================================================================================
## Dependent variable:
## ----------------------------------------------
## Outcome: intention
## (1) (2)
## --------------------------------------------------------------------------------
## groupCSR positive 0.064 0.311*
## (0.089) (0.158)
##
## groupPerformance negative -0.202** 0.006
## (0.084) (0.152)
##
## groupPerformance positive 0.356*** 0.420***
## (0.084) (0.147)
##
## psu_rul 0.351**
## (0.169)
##
## psu_op 0.381**
## (0.149)
##
## age -0.002 -0.002
## (0.002) (0.002)
##
## gender 0.086 0.057
## (0.064) (0.065)
##
## education 0.032 0.039
## (0.033) (0.035)
##
## income 0.0001 0.00004
## (0.0001) (0.0001)
##
## investment -0.00001 -0.00001
## (0.00001) (0.00001)
##
## nisa 0.149** 0.140**
## (0.060) (0.060)
##
## groupCSR positive:psu_rul -0.340
## (0.229)
##
## groupPerformance negative:psu_rul -0.311
## (0.220)
##
## groupPerformance positive:psu_rul -0.108
## (0.216)
##
## groupCSR positive:psu_op -0.410*
## (0.212)
##
## groupPerformance negative:psu_op -0.352*
## (0.202)
##
## groupPerformance positive:psu_op -0.141
## (0.200)
##
## Constant 0.100 -0.123
## (0.139) (0.169)
##
## --------------------------------------------------------------------------------
## Observations 236 236
## R2 0.251 0.289
## Adjusted R2 0.221 0.233
## Residual Std. Error 0.432 (df = 226) 0.429 (df = 218)
## F Statistic 8.397*** (df = 9; 226) 5.208*** (df = 17; 218)
## ================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01III. 質的変数（離散型確率変数）を従属変数とする分析
質的変数とは何か？：ダミー変数（二値変数、二項選択型変数）、名義変数、順序変数を含み、値と値との距離が任意に定まるもの。元のデータは、ほぼ文字情報であり、ここでの\(\texttt{outcome_intention}\)変数は典型的な質的変数で、
線形モデルと非線形モデルの対照：線形モデルでは量的変動を説明・予測するのに対して、非線形モデルでは発生確率を説明・予測することに主眼が置かれる。例示によって、対比してみよう。
1. ロジスティック回帰分析：従属変数の形態がダミー変数（二値変数、二項選択型変数）の場合
\(Y\)が0または1からなる、ロジスティック回帰モデルは以下の通り；
\[log(\frac{p(Y=1)}{1-p(Y=1)})=\beta_{o}+\beta_{1}x_{1}+\beta_{2}x_{2}+...+\beta_{k}x_{k}, i=1,...,n.\]
\(Y=\{1,0\}\)のもと、\(Y=1\)になる確率\(p\)の説明が目的。
#====================
# Logistic regression
#====================
# ---------------------------------------
# 1) Base-line logistic regression model
# ---------------------------------------
model_logit <- df_csr %>%
dplyr::select(outcome_intention, group, perf_pos, perf_neg, csr_pos, csr_neg, control_group, psu_rul, psu_op, group, age, gender, education, income, investment, nisa) %>%
na.omit() %>%
with(glm(
outcome_intention ~ perf_pos + perf_neg + csr_pos + age + gender + education + income + investment + nisa,
family = binomial(link = "logit")
)
)
# ----------------------------------------------------
# 2) Logistic regression model with interaction terms
# ----------------------------------------------------
model_logit_int <- df_csr %>%
dplyr::select(outcome_intention, group, perf_pos, perf_neg, csr_pos, csr_neg, control_group, psu_rul, psu_op, group, age, gender, education, income, investment, nisa) %>%
na.omit() %>%
with(glm(
outcome_intention ~ perf_pos + perf_neg + csr_pos +
psu_rul +
perf_pos:psu_rul + perf_neg:psu_rul + csr_pos:psu_rul +
age + gender + education + income + investment + nisa,
family = binomial(link = "logit")
)
)
stargazer(
model_logit,
model_logit_int,
type = "text",
title = "Logistic Regression Results",
dep.var.labels = "Outcome: intention",
out.header = FALSE
)##
## Logistic Regression Results
## ==============================================
## Dependent variable:
## ----------------------------
## Outcome: intention
## (1) (2)
## ----------------------------------------------
## perf_pos 1.583*** 1.584***
## (0.437) (0.520)
##
## perf_neg -1.274** -1.133*
## (0.505) (0.607)
##
## csr_pos 0.270 0.430
## (0.449) (0.548)
##
## psu_rul 0.664
## (0.758)
##
## age -0.011 -0.009
## (0.010) (0.010)
##
## gender 0.494 0.486
## (0.359) (0.361)
##
## education 0.195 0.159
## (0.188) (0.190)
##
## income 0.0003 0.0003
## (0.0003) (0.0003)
##
## investment -0.00004 -0.00004
## (0.0001) (0.0001)
##
## nisa 0.831** 0.820**
## (0.334) (0.339)
##
## perf_pos:psu_rul -0.120
## (0.972)
##
## perf_neg:psu_rul -0.513
## (1.085)
##
## csr_pos:psu_rul -0.569
## (0.978)
##
## Constant -2.064*** -2.137***
## (0.790) (0.820)
##
## ----------------------------------------------
## Observations 236 236
## Log Likelihood -125.809 -125.005
## Akaike Inf. Crit. 271.617 278.010
## ==============================================
## Note: *p<0.1; **p<0.05; ***p<0.01これらの推定結果からでも、どの変数が統計的に有意かを確かめることができる。では、リンク関数であるロジスティック関数に当てはめるとはどういうことだろうか？
group: Performance positiveの係数値は「1.583**」。これは先ほどのロジスティック関数の式に当てはめると…
\[log(\frac{p(Y=1)}{1-p(Y=1)})=-2.064+1.583x_{per.good}\]
\[\Leftrightarrow p=\frac{exp(-2.064+1.583x_{per.good})}{1+exp(-2.064+1.583x_{per.good})}\]
上式を利用してモデル化された曲線を作図すると、以下のようになる。
# ------------------------------------------
# 1) setting the range of perf_pos variable
# ------------------------------------------
x.new <- seq(min(df_csr$perf_pos, na.rm = TRUE),
max(df_csr$perf_pos, na.rm = TRUE),
length.out = 100)
# ---------------------------------------------------------------
# 2) create newdata with mean and median values of all variables
# ---------------------------------------------------------------
newdata <- df_csr %>%
summarise(
perf_neg = mean(perf_neg, na.rm = TRUE),
csr_pos = mean(csr_pos, na.rm = TRUE),
csr_neg = mean(csr_neg, na.rm = TRUE),
control_group = 0,
psu_rul = 0,
psu_op = 0,
age = mean(age, na.rm = TRUE),
gender = 0, # 男性を0, 女性を1でコーディングしているなら
education = mean(education, na.rm = TRUE),
income = mean(income, na.rm = TRUE),
investment = mean(investment, na.rm = TRUE),
nisa = 0
) %>%
slice(rep(1, length(x.new))) %>%
mutate(perf_pos = x.new)
# -------------------------------------------------------------------
# 3) Get link‐scale fit and SE, then transform to probability and CI #--------------------------------------------------------------------
pred <- predict(model_logit,
newdata = newdata,
type = "link",
se.fit = TRUE)
newdata <- newdata %>%
mutate(
fit_link = pred$fit,
se_link = pred$se.fit,
logit_low = fit_link - 1.96 * se_link,
logit_high = fit_link + 1.96 * se_link,
pred_prob = plogis(fit_link),
prob_low = plogis(logit_low),
prob_high = plogis(logit_high)
)
# --------
# 4) Plot
#---------
ggplot(newdata, aes(x = perf_pos, y = pred_prob)) +
# CI ribbon
geom_ribbon(aes(ymin = prob_low, ymax = prob_high),
fill = "steelblue", alpha = 0.2) +
# prediction line
geom_line(color = "steelblue", size = 1.2) +
# jittered raw points
geom_jitter(
data = df_csr,
aes(x = perf_pos, y = outcome_intention),
inherit.aes = FALSE,
width = 0.1,
height = 0.02,
alpha = 0.2,
color = "grey40",
size = 1.5
) +
labs(
x = "Treatment (Performance(+))",
y = "Predicted Probability of Purchase Intention",
title = "Logistic Predictions with 95% CI"
) +
theme_bw(base_size = 14) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(hjust = 0.5)
)
さらに知りたいことは、党派性をある値に固定し、共変量を平均や中央値にセットしたときの処置効果である。上式\(\beta_0+\beta_1x_1+\beta_2x_2+...\beta_kx_k,\quad i=1,...,n\)をもとにして、
\[ \eta = \beta_0+\beta_1x_1+\beta_2x_2+...+\beta_kx_{k} \]
とする。この時の\(\eta\)を線形予測子と呼び、この線形予測子をもとに、予測確率、
\[ \hat{p}=\mathrm{Pr}(y=1|\mathbf{X})=\frac{\exp(\eta)}{1+exp(\eta)} \]
を、それぞれの要因について計算することにより予測確率を求める。最も関心があるのは、党派性のもとでの処置によるアウトカムの発生確率であり、以下によって計算できる。
# ----------------------------------------------------
# 1) Logistic regression model with interaction terms
# ----------------------------------------------------
model_logit_int <- df_csr %>%
dplyr::select(outcome_intention, group, perf_pos, perf_neg, csr_pos, csr_neg, control_group, psu_rul, psu_op, group, age, gender, education, income, investment, nisa) %>%
na.omit() %>%
with(glm(
outcome_intention ~ group +
psu_rul +
group:psu_rul +
age + gender + education + income + investment + nisa,
family = binomial(link = "logit")
)
)
# -------------------------
# 2) generate predictions
# -------------------------
preds <- ggpredict(model_logit_int, terms = c("group", "psu_rul"))
# -------------------------
# 3) plot results
# -------------------------
p <- plot(preds) +
labs(
x = "Group",
y = "Predicted Probability of Intention",
colour = "In-partisanship or not",
fill = "In-partisanship or not"
) +
theme_bw(base_size = 14) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "bottom"
)
print(p)
2. プロビットモデルを使った分析
もう1つの二値変数を従属変数とする推定はプロビット・モデル(probit model)。プロビットモデルでは、推定式の誤差項\(\mu_{i}\)が正規分布に従うと仮定し、確率を推定する。そのための確率密度関数は、以下の通り；
\[f(x)=\frac{1}{\sqrt{2\pi\sigma}}exp\{-\frac{(x-\mu)^2}{2\sigma^2} \}.\]
推定式は上式を積分したもので求められる；
\[P(Y_{i}=1)=\frac{1}{\sqrt{2\pi\sigma}}\int^{\beta_{0}+\beta_{1}X_{i}}_{-\infty}exp\{-\frac{(x-\mu)^2}{2\sigma^2}\}dx\]
Rでのプロビット分析の実行:分析結果は、上記のロジット・モデルの場合と統計的に有意性において大きく変わらないことを確かめてみよう。
#===============================================
# Comparison of Logistic model and Probit model
#===============================================
# -----------------------------
# 1) Logistic regression model
# -----------------------------
model_logit <- df_csr %>%
dplyr::select(outcome_intention, group, perf_pos, perf_neg, csr_pos, csr_neg, control_group, psu_rul, psu_op, group, age, gender, education, income, investment, nisa) %>%
na.omit() %>%
with(glm(
outcome_intention ~ perf_pos + perf_neg + csr_pos + age + gender + education + income + investment + nisa,
family = binomial(link = "logit")
)
)
# ----------------------------------------------------
# 2) Probit model
# ----------------------------------------------------
model_probit <- df_csr %>%
dplyr::select(outcome_intention, group, perf_pos, perf_neg, csr_pos, csr_neg, control_group, psu_rul, psu_op, group, age, gender, education, income, investment, nisa) %>%
na.omit() %>%
with(glm(
outcome_intention ~ perf_pos + perf_neg + csr_pos + age + gender + education + income + investment + nisa,
family = binomial(link = "probit")
)
)
stargazer(
model_logit,
model_probit,
type = "text",
title = "Comparison of logistic regression and probit model results",
dep.var.labels = "Outcome: intention",
out.header = FALSE
)##
## Comparison of logistic regression and probit model results
## ==============================================
## Dependent variable:
## ----------------------------
## Outcome: intention
## logistic probit
## (1) (2)
## ----------------------------------------------
## perf_pos 1.583*** 0.978***
## (0.437) (0.261)
##
## perf_neg -1.274** -0.692**
## (0.505) (0.285)
##
## csr_pos 0.270 0.188
## (0.449) (0.271)
##
## age -0.011 -0.007
## (0.010) (0.006)
##
## gender 0.494 0.260
## (0.359) (0.209)
##
## education 0.195 0.105
## (0.188) (0.110)
##
## income 0.0003 0.0002
## (0.0003) (0.0002)
##
## investment -0.00004 -0.00002
## (0.0001) (0.00003)
##
## nisa 0.831** 0.468**
## (0.334) (0.195)
##
## Constant -2.064*** -1.157**
## (0.790) (0.459)
##
## ----------------------------------------------
## Observations 236 236
## Log Likelihood -125.809 -126.239
## Akaike Inf. Crit. 271.617 272.479
## ==============================================
## Note: *p<0.1; **p<0.05; ***p<0.01WORK
- 党派性を考慮した上での処置のもとで、購入意欲の予測確率がどのようになるかを検討してほしい。以下のような図が得られただろうか。
3. 順序ロジット・モデルを使った分析
実証モデルにおいて、しばしば従属変数に順序値を使う場合がある。この実験においても、処置の後に内閣への支持を聞いており、
支持する＝2
どちらともいえない = 1
支持しない = 0
として順序値を設定し、これを結果変数とする分析が可能である。順序値の特性は、単純な離散型の変数というだけでなく、連続的な推移の中に上記であれば2個の「閾値（仕切り）・thresholds」があるとも理解できるように、連続性の性質も有する変数である。
但し、推定に際しては、0～1、1～2への切り替わりに際して、どの程度の効果があるのかを測定するため、上記のロジット・モデルまたはプロビット・モデルと組み合わせて、順序ロジット・モデル、順序プロビット・モデルとして推定することが多い。
順序ロジット・モデルは、結果変数、\(Y_i\)が複数カテゴリ(\(1,2,..,J\))を撮るとき、切片（閾値）を\(\alpha_1<...\alpha_{J-1}\)とするとき、次式によって表される。
\[ \Pr\bigl(Y_i \le j \mid \mathbf{x}_i\bigr) = \frac{1}{1 + \exp\bigl[-\bigl(\alpha_j - \mathbf{x}_i^\top \boldsymbol{\beta}\bigr)\bigr]} = \mathrm{logit}^{-1}\!\bigl(\alpha_j - \mathbf{x}_i^\top \boldsymbol{\beta}\bigr), \quad j = 1, \dots, J-1. \]
順序ロジット推定は、以下のように\(\texttt{MASS}\) Package内のpolr関数によって行う（他にも可能な関数は多数）。以下では、係数だけでは効果が判別し難いため、オッズ比(odds-ratio)を示す。1.0倍から考えて、どの程度当該要因によって帰結が起こりやすくなるのか否かを表す比率である。
library(MASS)
# -----------------------
# 1) Ordered logit model
# -----------------------
df_csr <- df_csr %>%
mutate(cabap = factor(cabap, ordered = TRUE))
model_ologit <- df_csr %>%
polr(
formula = cabap ~
perf_pos +
perf_neg +
csr_pos +
age +
gender +
education +
income +
investment +
nisa,
data = .,
Hess = TRUE
)
# ----------------------------------------------
# 2) Extract coefficients, SEs, and odds‑ratios
# ----------------------------------------------
coefs <- coef(model_ologit) #係数の抽出
ses <- sqrt(diag(vcov(model_ologit))) #seの抽出
ors <- exp(coefs) #オッズ比の算出
# ----------------------------------------------
# 3) Output tables with coeffients and odds-ratio
# ----------------------------------------------
stargazer(
list(model_ologit, model_ologit),
coef = list(coefs, ors),
se = list(ses, rep(NA_real_, length(ses))),
column.labels = c("Estimate", "Odds Ratio"),
dep.var.labels = "CABAP (ordered logit)",
covariate.labels = c(
"Threshold 0 | 1",
"Threshold 1 | 2",
"Performance (+)",
"Performance (–)",
"CSR (+)",
"Age",
"Gender (Male=1)",
"Education",
"Income",
"Investment",
"NISA Investor"
),
omit.stat = c("LL", "AIC"),
digits = 3,
type = "text",
intercept.bottom = FALSE,
intercept.top = TRUE,
out.header = FALSE
)##
## ============================================
## Dependent variable:
## ----------------------------
## CABAP (ordered logit)
## Estimate Odds Ratio
## (1) (2)
## --------------------------------------------
## Threshold 0 | 1 0.354** 1.425
## (0.174)
##
## Threshold 1 | 2 0.068 1.070
## (0.160)
##
## Performance (+) 0.001 1.001
## (0.157)
##
## Performance (–) 0.010** 1.010
## (0.005)
##
## CSR (+) -0.029 0.971
## (0.195)
##
## Age 0.064 1.066
## (0.072)
##
## Gender (Male=1) -0.0003 1.000
## (0.0003)
##
## Education 0.0001 1.000
## (0.0001)
##
## Income 0.107 1.112
## (0.194)
##
## --------------------------------------------
## Observations 443 443
## ============================================
## Note: *p<0.1; **p<0.05; ***p<0.01# ----------------------------------------------
# 1) Ordered logit model with interaction terms
# ----------------------------------------------
model_ologit_int <- df_csr %>%
polr(
formula = cabap ~
psu_rul +
perf_pos +
perf_neg +
csr_pos +
perf_pos * psu_rul +
perf_neg * psu_rul +
csr_pos * psu_rul +
age +
gender +
education +
income +
investment +
nisa,
data = .,
Hess = TRUE
)
# ----------------------------------------------
# 2) Extract coefficients, SEs, and odds‑ratios
# ----------------------------------------------
coefs_int <- c(model_ologit_int$zeta, model_ologit_int$coefficients)
ses_int <- sqrt(diag(vcov(model_ologit_int)))
ors_int <- exp(coefs_int)
# ----------------------------------------------
# 3) Output tables with coeffients and odds-ratio
# ----------------------------------------------
stargazer(
list(model_ologit_int, model_ologit_int),
coef = list(coefs_int, ors_int),
se = list(ses_int, rep(NA_real_, length(ses_int))),
column.labels = c("Estimate", "Odds Ratio"),
dep.var.labels= "CABAP (ordered logit with interactions)",
covariate.labels = c(
"Threshold 0 | 1",
"Threshold 1 | 2",
# main effects
"In‑party",
"Performance (+)",
"Performance (–)",
"CSR (+)",
"Age",
"Gender (Male=1)",
"Education",
"Income",
"Investment",
"NISA Investor",
# interaction terms
"Performance (+) × In‑party",
"Performance (–) × In‑party",
"CSR (+) × In‑party"
),
intercept.bottom = FALSE,
intercept.top = TRUE,
omit.stat = c("LL","AIC"),
digits = 3,
type = "text",
out.header = FALSE
)##
## ===================================================================
## Dependent variable:
## ----------------------------------------
## CABAP (ordered logit with interactions)
## Estimate Odds Ratio
## (1) (2)
## -------------------------------------------------------------------
## Threshold 0 | 1 1.629*** 5.101
## (0.174)
##
## Threshold 1 | 2 0.356** 1.427
## (0.153)
##
## In‑party 0.176*** 1.192
## (0.146)
##
## Performance (+) 0.161 1.175
## (0.147)
##
## Performance (–) 0.011 1.011
## (0.005)
##
## CSR (+) -0.173*** 0.841
## (0.204)
##
## Age 0.0002** 1.000
## (0.079)
##
## Gender (Male=1) -0.0003 1.000
## (0.0003)
##
## Education 0.00005 1.000
## (0.0001)
##
## Income 0.201*** 1.222
## (0.200)
##
## Investment -0.644*** 0.525
## (0.077)
##
## NISA Investor -0.941*** 0.390
## (0.074)
##
## Performance (+) × In‑party -0.936*** 0.392
## (0.080)
##
## -------------------------------------------------------------------
## Observations 416 416
## ===================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01「1．ロジスティック回帰分析」でも行ったように、順序ロジット・モデルの推定結果をもとに、交差項のもとでの各選択肢についての選択の予測確率を求め、図示する。ロジット分析の場合と異なるのは、順序の各水準ごとに予測確率が異なることである。facet_wrap()関数を使って、選択肢ごとの図を折りたたむ。
library(ggeffects)
# ----------------------------------------------------
# 1) prediction with ggpredict() in ggeffects package
# ----------------------------------------------------
pred_pp <- ggpredict(
model_ologit_int,
terms = c("perf_pos [0,1]", "psu_rul [0,1]")
)
# -------------------------------------
# 2) tidy up the data.frame and relabel
# -------------------------------------
pred_pp_df <- as.data.frame(pred_pp) %>%
mutate(
cabap = factor(response.level,
levels = c("0","1","2"),
labels = c("Disapprove(0)","Neither \n dis-/approve(1)","Approve(2)")),
# PSU ruling party labels
psu_rul = factor(group,
levels = c("0","1"),
labels = c("Out-party","In-party"))
)
# 3) Dot‐and‐error‐bar plot, dodged by psu_rul, faceted by cabap
ggplot(pred_pp_df,
aes(x = factor(x),
y = predicted,
colour = psu_rul)) +
geom_point(position = position_dodge(width = 0.6),
size = 3) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
width = 0.2,
position = position_dodge(width = 0.6)) +
facet_wrap(~ cabap, nrow = 1) +
scale_x_discrete(name = "Treatment (Performance: positive)") +
scale_y_continuous(name = "Predicted probability", labels = scales::percent_format(accuracy = 1)) +
scale_colour_manual(name = "Partisanship",
values = c("Out-party" = "#E69F00",
"In-party" = "#56B4E9")) +
labs(title = "Predicted probability of cabinet support \n (Performance evaluation: positive)") +
theme_bw(base_size = 14) +
theme(
legend.position = "bottom",
axis.text.x = element_text(angle = 45, hjust = 1),
strip.background = element_rect(fill = "grey90"),
strip.text = element_text(face = "bold")
)
他の処置についても、上記の図と同じようにpredicted probabilityを算出してみよう。下のような図を示せただろうか。
他の処置についてのコードも含めて、コードの効率化を試みてほしい。
