在 R 中使用 Anova 计算 I 型错误时出现问题答案

【问题标题】：Troubles computing Type I error with Anova in R在 R 中使用 Anova 计算 I 型错误时出现问题
【发布时间】：2017-03-10 10:23:10
【问题描述】：

我正在尝试为简单的 Anova 测试计算 I 类错误，但我得到了奇怪的结果。

假设我想检测剂量 (DOSE) 对 4 家不同医院 (HOP) 中观察变量 (obs) 的影响。在药物不会影响观察到的变量但医院可以生成这样的数据集的假设下：

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 次并且我将 pvalue 小于或等于 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%。

你能帮我解决这个问题吗？

最好，西蒙

【问题讨论】：

您应该在开头添加set.seed(number)，以便我们得到相同的答案。例如，我得到了 2 个。
感谢您的回复，让我们使用 set.seed(1234) 给出 1

标签： r statistics anova

【解决方案1】：

好的，我终于设法解决了这个问题。这是由于不同的错误：

1) 你需要在方差分析中使用因子(DOSE)

2）我模拟数据集的方式不正确：

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) 您需要考虑 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)

感谢您的宝贵时间，西蒙

【讨论】：