このシミュレーションの目的は一つのグループに何人くらい配分したら全体の平均値に近づくのかを経験的に確認するためのものです。
趣味でしたものなので授業とは全く関係ありません。Rを覚える必要もありません。理解しなくてもいいです。Rを触ったことがないと、何をやってるのか分からないと思うので読まなくてもいいです。
それでは、さっそく
作図のためにggplot2パッケージの呼び出す。もし、ggplot2がR上にない場合、まずinstall.packages("ggplot2")でインストールしておく。
library(ggplot2)次に実験群の数n.groupを設定する。ここでは統制群, 実験群1, 実験群2ということで3に設定
n.group <- 3続いて、シミュレーション結果を保存するデータフレームの作成する。
result.df <- data.frame(n = rep(NA, 996),
sex.mean = rep(NA, 996),
age.mean = rep(NA, 996),
trust.mean = rep(NA, 996),
n1 = rep(NA, 996),
n2 = rep(NA, 996),
n3 = rep(NA, 996),
sex.mean1 = rep(NA, 996),
sex.mean2 = rep(NA, 996),
sex.mean3 = rep(NA, 996),
age.mean1 = rep(NA, 996),
age.mean2 = rep(NA, 996),
age.mean3 = rep(NA, 996),
trust.mean1 = rep(NA, 996),
trust.mean2 = rep(NA, 996),
trust.mean3 = rep(NA, 996))本格的にシミュレーションコードを書いてみる。一つの実験群に理論的に5から1000まで入るようにする。つまり、サンプルサイズは15から3000まで3ずつ増やしていく。
for (i in 5:1000){
# 架空のデータセットの生成する
# 各グループに大体x人が属し、３つのグループで分けるなら
# xの3倍のサンプルサイズが必要
sex <- sample(0:1, i * n.group, replace = TRUE)
age <- sample(20:70, i * n.group, replace = TRUE)
trust <- sample(1:5, i * n.group, replace = TRUE)
# グループを分ける
group.id <- sample(1:n.group, i * n.group, replace = TRUE)
# 架空データセットをデータフレームとして保存
temp.df <- data.frame(sex, age, trust, group.id)
# サンプルサイズ
sample.size <- i * n.group
# サンプル全体のsex, age, trustの平均値
sex.mean <- mean(temp.df$sex )
age.mean <- mean(temp.df$age )
trust.mean <- mean(temp.df$trust)
# 各グループにいくつのケースが属されているかを計算
n1 <- nrow(subset(temp.df, group.id == 1))
n2 <- nrow(subset(temp.df, group.id == 2))
n3 <- nrow(subset(temp.df, group.id == 3))
# 各グループごとのsex, age, trustの平均値を計算する
sex.mean1 <- mean(subset(temp.df, group.id == 1)$sex )
sex.mean2 <- mean(subset(temp.df, group.id == 2)$sex )
sex.mean3 <- mean(subset(temp.df, group.id == 3)$sex )
age.mean1 <- mean(subset(temp.df, group.id == 1)$age )
age.mean2 <- mean(subset(temp.df, group.id == 2)$age )
age.mean3 <- mean(subset(temp.df, group.id == 3)$age )
trust.mean1 <- mean(subset(temp.df, group.id == 1)$trust)
trust.mean2 <- mean(subset(temp.df, group.id == 2)$trust)
trust.mean3 <- mean(subset(temp.df, group.id == 3)$trust)
# 結果を保存
result.df[i-4, ] <- c(sample.size, sex.mean, age.mean, trust.mean,
n1, n2, n3,
sex.mean1, sex.mean2, sex.mean3,
age.mean1, age.mean2, age.mean3,
trust.mean1, trust.mean2, trust.mean3)
}結果をグラフで確認してみる。
まず、各グループにちゃんと3分の1ずつ割り当てられているかをみる。赤の実線は3分の1を表し、この実線に近いということは理想的に割り当てられていることを意味する。
# 各グループにちゃんと1/3ずつ入ったかを確認。
# 赤の実線に近いほど理想的な割合で配分されたということ
ggplot(data = result.df) +
geom_line(aes(x = n, y = n1/n, color = "red"), alpha = 0.4) +
geom_line(aes(x = n, y = n2/n, color = "blue"), alpha = 0.4) +
geom_line(aes(x = n, y = n3/n, color = "green"), alpha = 0.4) +
geom_hline(yintercept = 1/3, color = "red") +
geom_vline(xintercept = 750, color = "red") +
scale_color_discrete(name = "Group",
labels = c("Group1", "Group3", "Group2")) +
labs(x = "Sample size", y = "Sample size in each group") +
theme(legend.position = "bottom")
だいたいサンプルサイズが500~750くらいになると3つのグループは概ね3分の1ずつ割り当てられることが確認できる。
次は、割り当てられたグループの性別の平均値が均一か否かをみる。紫を除く3つの線が集まっていると3つのグループの性別の平均値は概ね一緒ということを意味する。赤い実線は理論値である0.5を表すが、そもそもサンプルの平均が0.5になっているという保障はないので注意が必要。
# 各グループのsex, age, trustの平均値。
# 紫(?)は全体の平均値
# 赤い実線は理論値
ggplot(data = result.df) +
geom_line(aes(x = n, y = sex.mean, color = "black"), alpha = 0.4) +
geom_line(aes(x = n, y = sex.mean1, color = "red"), alpha = 0.4) +
geom_line(aes(x = n, y = sex.mean2, color = "blue"), alpha = 0.4) +
geom_line(aes(x = n, y = sex.mean3, color = "green"), alpha = 0.4) +
geom_hline(yintercept = 0.5, color = "red") +
geom_vline(xintercept = 750, color = "red") +
scale_color_discrete(name = "Group",
labels = c("Group1", "Group3", "Group2", "Entire")) +
labs(x = "Sample size", y = "Mean of Sex") +
theme(legend.position = "bottom")
サンプルサイズが大きくなると紫を除く3つの線は赤い実線に集まっていくことが確認できる。ただし、サンプルサイズが15(最小値)から600に変化したときと1000から3000に変化した時の差をみると明らかに前者の方の差が大きい。つまり、3つのグループで、0/1変数の場合はサンプルサイズがだいたい500、つまり一つのグループに170人くらいいるとグループ間の均質性は保証され、そこからは劇的にサンプルサイズを増やさないかぎり、大きな変化はないと考えられる。
同じ手順でageとtrustも見てみよう。
ggplot(data = result.df) +
geom_line(aes(x = n, y = age.mean, color = "black"), alpha = 0.4) +
geom_line(aes(x = n, y = age.mean1, color = "red"), alpha = 0.4) +
geom_line(aes(x = n, y = age.mean2, color = "blue"), alpha = 0.4) +
geom_line(aes(x = n, y = age.mean3, color = "green"), alpha = 0.4) +
geom_hline(yintercept = 45, color = "red") +
geom_vline(xintercept = 750, color = "red") +
scale_color_discrete(name = "Group",
labels = c("Group1", "Group3", "Group2", "Entire")) +
labs(x = "Sample size", y = "Mean of Age") +
theme(legend.position = "bottom")
ggplot(data = result.df) +
geom_line(aes(x = n, y = trust.mean, color = "black"), alpha = 0.4) +
geom_line(aes(x = n, y = trust.mean1, color = "red"), alpha = 0.4) +
geom_line(aes(x = n, y = trust.mean2, color = "blue"), alpha = 0.4) +
geom_line(aes(x = n, y = trust.mean3, color = "green"), alpha = 0.4) +
geom_hline(yintercept = 3, color = "red") +
geom_vline(xintercept = 750, color = "red") +
scale_color_discrete(name = "Group",
labels = c("Group1", "Group3", "Group2", "Entire")) +
labs(x = "Sample size", y = "Mean of Trust") +
theme(legend.position = "bottom")
結果的にみると一つのグループに250人くらい割り当てられるサンプルサイズ(ここではグループは3つだから750人くらい)で、各グループのサンプルサイズと各変数の期待値がある程度安定することが確認できる。
もしかしたら、測定質問の尺度(つまり、分散)と関係があるかも知れない。二項変数 (sex) やカテゴリ変数 (trust) は確かにN = 500、つまりグループごとに166人くらいで安定するように見えるが、連続変数 (age) は250人くらいだと思う。尺度の最小値の最大値の幅が広いと分散が大きくなるせいではないかとも思う。これに関しては尺度の分散をいじりながらシミュレーションしてみる価値があるかも。