何故か 5-4 から再開。つうか
良書ですね。Python で云々とか緑な本とか読んだりしたいです。
とりあえず、何となく 5-4-8 以降をもごもごしつつ控えを。
データを用意して統計モデルを作る
以下?
d2 <- data.frame(
length = c(2, 4, 10, 12, 6, 8),
medicine = c("A", "A", "B", "B", "C", "C")
)
modelANOVA <- lm(length ~ medicine, data=d2)
modelNull <- lm(length ~ 1, data=d2)
標本から F 比を計算
凄い、
anova(modelANOVA)
から F value
だけ取り出せるのか。
> anova(modelANOVA)$"F value"[1]
[1] 16
R は配列参照が 1 origin なのですね。
シミュレーションデータ作成
こう、なのかな。
> simData <- cbind(
+ simulate(modelNull, 1),
+ d2$medicine)
で、以下なのかな。
> colnames(simData) <- c("length", "medicine")
> simData
length medicine
1 7.073384 A
2 9.022781 A
3 1.486284 B
4 2.100310 B
5 10.555131 C
6 5.317685 C
ええと、値が変わらないのですが良いのかどうか。simulate
すると値は変わりますね。
F 比を多数計算する
合わせると以下なのかな。
Nsim <- 10000
sim <- simulate(modelNull, Nsim)
simFValue <- numeric()
for(i in 1:Nsim){
simData <- cbind(
sim[i],
c("A", "A", "B", "B", "C", "C")
)
colnames(simData) <- c("length", "medicine")
model <- lm(length ~ medicine, data=simData)
simFValue[i] <- anova(model)$"F value"[1]
}
で、値を見ても意味分からないので、以下で図示。
hist(
subset(simFValue, simFValue < 20),
main="F 比のヒストグラム"
)
スクリプトがバグッててアレ。繰返しにブラケット付けるの忘れてました。
p 値の算出
16 以上の F 比の数。
> sum(simFValue>=16)
[1] 262
Nsim で割れば p 値が出る模様。
> sum(simFValue>=16)/Nsim
0.0262
正規分布から派生した確率分布
正規分布偉い、らしい。とは言え、「母集団分布が推定できるわけではない」らしい。
つうかこの 5-5 章はちゃんと確認した方が良さげ。何度か目を通して別途必要があればエントリ入れます。