この前、神戸大学のワークショップでStanford大学のIyenger先生は「実験データで回帰分析するな」と言われた覚えがある。そりゃ実験のデータセットは実験のために設計されたもんだし、別に疑問もなかったが、もうちょっと考えたら「別に回帰分析やってもいいんじゃない?」と思うようになった。

それでシミュレーション


ラーメン二郎の写真を見せることによってラーメン好き度がどう変化するかというリサーチクエスチョンを設定する。サンプルサイズ(\(N\))は500で、実験群と統制群を無作為配分する。ラーメン好き度は性別、年齢、教育水準という3つの変数の関数だとする。パラメータは以下のように設定する

\[\beta_0 = 1\] \[\beta_1 = 2\] \[\beta_2 = 0.5\] \[\beta_3 = 0.8\]

以上はそれぞれ切片、性別の係数、年齢の係数、教育水準の係数である。また、ラーメン好き度は、

\[y \sim N(\mu, \sigma = 2)\] \[\mu_i = \beta_0 + \beta_1 x_{1, i} + \beta_2 x_{2, i} + \beta_3 x_{3, i}\]

とする。このデータジェネレーティングプロセスをRで表現すると、

b0 <- 1; b1 <- 2; b2 <- 0.5; b3 <- 0.8
x1 <- sample(0:1, 500, replace = TRUE)
x2 <- sample(20:90, 500, replace = TRUE)
x3 <- sample(c(6, 9, 12, 16, 18, 21), 500, replace = TRUE)

mu <- b0 + b1 * x1 + b2 * x2 + b3 * x3
y <- rnorm(500, mean = mu, sd = 2)

続いて、実験群と統制群を割り当てる。sが0だと統制群、1だと実験群だとする。

s <- sample(0:1, 500, replace = TRUE)

仮に刺激sの効果がラーメン好き度を約2くらい上げるとする。このsの効果もまた、

\[{coef}_s \sim N(\mu = 2, \sigma = 2)\]

だとする。そうなると最終的に出来るデータセットは、

sim.df <- data.frame(y, x1, x2, x3, s)
head(sim.df)
##          y x1 x2 x3 s
## 1 47.52974  1 69 12 1
## 2 43.17594  0 79  6 0
## 3 52.39953  1 90 12 1
## 4 49.26891  0 70 18 0
## 5 40.12571  0 57 16 1
## 6 40.91571  1 52 16 0
for(i in 1:nrow(sim.df)){
        if(sim.df$s[i] == 1){
                sim.df$y[i] <- sim.df$y[i] + rnorm(1, 2, 2)
        }else{}
}
head(sim.df)
##          y x1 x2 x3 s
## 1 52.98304  1 69 12 1
## 2 43.17594  0 79  6 0
## 3 53.47549  1 90 12 1
## 4 49.26891  0 70 18 0
## 5 42.07727  0 57 16 1
## 6 40.91571  1 52 16 0

次に、実験群の割り当てを確認してみる。

# 統制群と実験群のサイズ
table(sim.df$s)
## 
##   0   1 
## 263 237
# 性別に差はあるか
t.test(sim.df[which(sim.df$s == 0), 2], sim.df[which(sim.df$s == 1), 2])
## 
##  Welch Two Sample t-test
## 
## data:  sim.df[which(sim.df$s == 0), 2] and sim.df[which(sim.df$s == 1), 2]
## t = -0.20189, df = 492.79, p-value = 0.8401
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.09676212  0.07872936
## sample estimates:
## mean of x mean of y 
## 0.5437262 0.5527426
# 年齢に差はあるか
t.test(sim.df[which(sim.df$s == 0), 3], sim.df[which(sim.df$s == 1), 3])
## 
##  Welch Two Sample t-test
## 
## data:  sim.df[which(sim.df$s == 0), 3] and sim.df[which(sim.df$s == 1), 3]
## t = 0.31185, df = 496.95, p-value = 0.7553
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.029235  4.172262
## sample estimates:
## mean of x mean of y 
##  54.82890  54.25738
# 教育水準に差はあるか
t.test(sim.df[which(sim.df$s == 0), 4], sim.df[which(sim.df$s == 1), 4])
## 
##  Welch Two Sample t-test
## 
## data:  sim.df[which(sim.df$s == 0), 4] and sim.df[which(sim.df$s == 1), 4]
## t = -0.36698, df = 492.94, p-value = 0.7138
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.0923799  0.7485381
## sample estimates:
## mean of x mean of y 
##  13.18251  13.35443

確認してみたら、かなりいい感じで割り当てられている(seedを設定していないので、完成した頃にはあまり良くないかも知れない)。それでは、本格的にt検定で差を見てみる。

t.test(sim.df[which(sim.df$s == 0), 1], sim.df[which(sim.df$s == 1), 1])
## 
##  Welch Two Sample t-test
## 
## data:  sim.df[which(sim.df$s == 0), 1] and sim.df[which(sim.df$s == 1), 1]
## t = -1.8449, df = 483.7, p-value = 0.06566
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.9258140  0.1236121
## sample estimates:
## mean of x mean of y 
##  40.23013  42.13123

t検定をやってみると、2つのグループ間の差は統計的には有意でない。それなら回帰分析をやってみたらどうだろうか。

lm.result <- lm(y ~ x1 + x2 + x3 + s, data = sim.df)
summary(lm.result)
## 
## Call:
## lm(formula = y ~ x1 + x2 + x3 + s, data = sim.df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9617 -1.7365  0.0396  1.6700  8.3591 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.126402   0.467903   2.407   0.0164
## x1          1.991548   0.235246   8.466 2.91e-16
## x2          0.498884   0.005712  87.344  < 2e-16
## x3          0.809224   0.022184  36.477  < 2e-16
## s           2.029141   0.231817   8.753  < 2e-16
## 
## Residual standard error: 2.588 on 495 degrees of freedom
## Multiple R-squared:  0.9497, Adjusted R-squared:  0.9493 
## F-statistic:  2337 on 4 and 495 DF,  p-value: < 2.2e-16

ちゃんと2くらいの係数が得られた。

ここまでの結果はあくまでも乱数を生成したもので、場合によってはt検定でも有意な結果が十分得られる。ただし、ラーメン二郎の写真の効果には人それぞれ違うと考えた方が妥当かも知れない。もし、その分散が大きい場合はt検定では結果が得られにくくなる可能性があると思う。つまり、刺激の効果の分散以外の分散も分析に含まれることによって、標準誤差が大きくなり、きちんとした結果が出ないかも知れない。

また、データ生成過程を上とは違って、

\[\mu_i = \beta_0 + \beta_1 x_{1, i} + \beta_2 x_{2, i} + \beta_3 x_{3, i} + \beta_4 s_i\] \[y \sim N(\mu_i, \sigma)\]

にしても同じであろう。ただし、この場合は分散を共有(?)することもあって、効果量自体にバイアスがかかるかも知れない(やってみてないからよく分からないが…)。また、もう一つのありうる問題は、刺激の効果(係数)が何かの関数である場合がある。この場合はもっとややこしくなると考えられる。

もちろん、その意味でt検定でも有意な結果が得られるということは、ロバストであることを意味するのかも知れない。ただし、未だよく分からないことは刺激の効果の生成過程がよく分からないことである。その意味で、今回のシミュレーションは設計段階から間違っているのかも知れない。