研究概要
大学時代に実験経済学で行った実験結果のデータがUSBに入っていたので、振り返って分析などをしてみたいと思います。
研究目的
ピア効果に関して、競争相手が自分よりも秀でた人がいいのか劣った人がいいのかを確かめる。
実験方法
・1分間で100マス計算を2セット解いてもらう。（めちゃ速い人には3枚目も渡した）
・実験開始後、実験対象のクラスによって、途中で「平均的なクラスは○○マスまで進んでいます！」とアナウンスします。アナウンスすることで、競争相手のレベルを知り、焦るなり余裕を感じるなりしてもらおうという計画です。
なお、対照群はアナウンスをしていません。アナウンス内容は「平均告知（１８秒）、上告知（１５秒）、超上告知（１２秒）、下告知（２０秒）」と4パターンとなります。
・計算が間違っているものは加点しません。
実験対象
某国立大学の経済学部生の1~2年の必修科目履修者217名（先生に交渉して授業の開始5分を頂いて実験を行いました。）
内部進学やスポ専などがない分、計算能力的にある程度近い集団ではないかと思われます。
検証方法
アナウンスごとに100マス計算の点数の水準が変わりうるのかを回帰分析などで判断。
データ可視化
以下、実験カテゴリごとの略記です。
下告知（２０秒）・・・slow_20
上告知（１５秒）・・・fast_15
超上告知（１２秒）・・・fastest_12
平均告知（１８秒）・・・average_18
対照群・・・baseline
データ構造の確認です。
|
1 2 3 4 5 6 |
> str(dataset) 'data.frame': 217 obs. of 4 variables: $ categories : chr "baseline" "baseline" "baseline" "baseline" ... $ points : int 72 79 81 81 98 99 100 101 102 104 ... $ errors : int 0 0 0 0 2 1 0 0 1 0 ... $ genuine_points: int 72 79 81 81 96 98 100 101 101 104 ... |
平均値、中央値、標準偏差、サンプルサイズを出してみます。
|
1 2 3 4 5 6 7 8 9 10 11 12 |
> dataset %>% group_by(categories) %>% summarise(average=mean(genuine_points), + median=median(genuine_points), + stdev=sd(genuine_points), + sample=n()) # A tibble: 5 × 5 categories average median stdev sample <chr> <dbl> <dbl> <dbl> <int> 1 average_18 120.8605 117 28.08455 43 2 baseline 121.0000 121 22.19109 46 3 fast_15 123.0357 126 24.32355 56 4 fastest_12 123.6774 125 24.09756 31 5 slow_20 126.3902 126 23.20547 41 |
中央値で見てみると、baselineに対してわずかですが点数に違いがありそうに見えます。
実験種別で点数に関するヒストグラムと確率密度関数を確認してみます。
|
1 2 3 4 5 6 7 8 |
library(ggplot2) g <- ggplot(data = dataset, aes(x = genuine_points, y = ..density..)) + geom_histogram(alpha = 0.5,position = "identity") + geom_density(alpha = 0) g + facet_wrap(~categories,nrow=5) |
baselineが多峰性がありそうなのが気になります。average_18は低そうに見えますね。
RStanで重回帰
『StanとRでベイズ統計モデリング』にあるコードを参考にしています。正規分布を事前分布にした線形回帰モデルです。
被説明変数が点数、説明変数が実験種別のダミー変数だけからなります。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 |
library(rstan) library(dummies) dummies <- dummy.data.frame(dataset, sep = "_", names = c("categories")) analytical_dataset <- dummies %>% select(categories_average_18, categories_fast_15, categories_fastest_12, categories_slow_20, genuine_points) data <- list(N=nrow(analytical_dataset), average_18=analytical_dataset$categories_average_18, fast_15=analytical_dataset$categories_fast_15, fastest_12=analytical_dataset$categories_fastest_12, slow_20=analytical_dataset$categories_slow_20, genuine_points=analytical_dataset$genuine_points) stan_code <- " data{ int N; //the number of student int<lower=0> genuine_points[N]; real<lower=0, upper=1> average_18[N]; real<lower=0, upper=1> fast_15[N]; real<lower=0, upper=1> fastest_12[N]; real<lower=0, upper=1> slow_20[N]; } parameters{ real b1; real b2; real b3; real b4; real b5; real<lower=0> sigma; } transformed parameters{ real mu[N]; for(n in 1:N) mu[n] = b1 + b2*average_18[n] + b3*fast_15[n] + b4*fastest_12[n] + b5*slow_20[n]; } model{ for(n in 1:N) genuine_points[n] ~ normal(mu[n], sigma); } " fit <- stan(model_code =stan_code, data=data, seed=1234) fit.summary <-data.frame(summary(fit)$summary) head(fit.summary,6) |
結果
traceplot(fit)でMCMCのサンプリング結果を確認する。
収束しているように見えます。
以下推定結果ですが、残念ながらベイズ予測区間において符号の逆転が起きていないものはなかったので、アナウンスによる効果があるとは言えないようです。ただ、slow_20の係数がおしいですね。少なくとも他の実験種別よりも、アナウンス効果があるかもしれないという考察に止まりそうです。
|
1 2 3 4 5 6 7 8 |
> head(fit.summary,6) mean se_mean sd X2.5. X25. X50. X75. X97.5. n_eff Rhat b1 121.0253094 0.1034836 3.460792 114.214104 118.793988 121.0140383 123.213054 127.902759 1118.429 0.9998152 b2 -0.1195318 0.1326887 4.957153 -10.009827 -3.372795 -0.1304577 3.117195 9.694184 1395.715 1.0004034 b3 1.9440496 0.1248433 4.661571 -7.494857 -1.043060 1.9756929 5.007503 10.827640 1394.228 1.0000830 b4 2.5460694 0.1436941 5.447930 -8.026616 -1.128587 2.4997953 6.220891 13.044467 1437.426 1.0008688 b5 5.3029965 0.1351662 4.988323 -4.479865 2.004280 5.3108476 8.543292 15.458379 1361.988 1.0004238 sigma 24.6158717 0.0211252 1.150481 22.501612 23.811858 24.5900106 25.368057 26.970197 2965.905 1.0017414 |
分布でも見てみます。アナウンス効果が0を確実に超えているとは言えないですね。
|
1 2 3 4 5 6 7 8 9 |
library (reshape) library(dplyr) library(ggplot2) post <- extract (fit, permuted = F) m.post <- melt (post) m.post <- m.post %>% filter(parameters %in% c("b1","b2","b3","b4","b5")) graph <- ggplot (m.post, aes(x = value)) graph <- graph + geom_density () + facet_grid(. ~ parameters, scales = "free") + theme_bw() plot (graph) |
結局、学部時代のレポートと結論は変わらないのですが、係数が0よりも大きい確率という観点で結果に向き合えたのは良かったと思います。
参考文献
