TroublesコンピューティングでのタイプIエラーAnova with R

Question

単純なAnovaテストでタイプIエラーを計算しようとしていますが、私は奇妙な結果を得ています。

4つの異なる病院（HOP）の観察された変数（obs）に対する投与量（DOSE）の影響を検出したいと仮定します。仮説の下で は、薬物が観測変数に影響を与えないことが、病院は、私はそのようなデータセットを生成することができます：私は、テストを観察し、変数やエキスに投与量の効果をANOVAを使用することができます

data.frame(
obs=c(rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1), 
rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1), 
rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1), 
rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1)), 
HOP=rep(1:4,100), 
DOSE=rep(c(0,15,30,50),each=100))->data 

p値：

summary(aov(data$obs~data$DOSE))[[1]][[5]][1]->pvalue 

私はそれを100回行うと、私はp値が小さいか等しい回数を合計した場合0.05私はタイプIエラーが発生しますと、この値は5に等しいはずである。

Allp<-NULL 
for (i in 1:100){ 
data.frame(
obs=c(rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1), 
rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1), 
rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1), 
rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1)), 
HOP=rep(1:4,100), 
DOSE=rep(c(0,15,30,50),each=100))->data 
summary(aov(data$obs~data$DOSE))[[1]][[5]][1]->pvalue 
Allp<-rbind(Allp,pvalue)} 
sum(Allp<=0.05) 

ただし、0または1！

私は病院の影響を想定していない試してみました：

Allp<-NULL 
for (i in 1:100){ 
data.frame(
obs=c(rnorm(25,0,1),rnorm(25,0,1),rnorm(25,0,1),rnorm(25,0,1), 
rnorm(25,0,1),rnorm(25,0,1),rnorm(25,0,1),rnorm(25,0,1), 
rnorm(25,0,1),rnorm(25,0,1),rnorm(25,0,1),rnorm(25,0,1), 
rnorm(25,0,1),rnorm(25,0,1),rnorm(25,0,1),rnorm(25,0,1)), 
HOP=rep(1:4,100), 
DOSE=rep(c(0,15,30,50),each=100))->data 
summary(aov(data$obs~data$DOSE))[[1]][[5]][1]->pvalue 
Allp<-rbind(Allp,pvalue)} 
sum(Allp<=0.05) 

そして、ここで私が期待される5％を取得します。

この問題を解決するお手伝いをしてください。

ベスト

、 サイモン

Answer 1

[OK]を、私は最終的に問題を解決するために管理しました。

1）あなたはANOVA

2）私は、データセットが間違っていたシミュレートされた方法で因子（DOSE）を使用する必要があります：これは違いによるミスをした

data.frame(
obs=c(rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1), 
rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1), 
rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1), 
rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1)), 
HOP=rep(rep(1:4,each=25),4), 
DOSE=rep(c(0,15,30,50),each=100))->data 

3）あなたをタイプ1エラーの適切な推定値を得るために考慮HOPの影響を取る必要があります：あなたの時間のために

set.seed(1234) 
Allp<-NULL 
for (i in 1:100){ 
data.frame(
obs=c(rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1), 
rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1), 
rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1), 
rnorm(25,0,1),rnorm(25,1,1),rnorm(25,2,1),rnorm(25,3,1)), 
HOP=rep(rep(1:4,each=25),4), 
DOSE=rep(c(0,15,30,50),each=100))->data 
summary(aov(data$obs~as.factor(data$DOSE)+as.factor(data$HOP)))[[1]][[5]][1]->pvalue 
Allp<-rbind(Allp,pvalue)} 
sum(Allp<=0.05) 

おかげで、 Simon