全嵌套方差分析

Question

我试图在NIST的下面一页的底部复制这个示例：

http://www.itl.nist.gov/div898/handbook/ppc/section2/ppc233.htm

使用Minitab和SPSS，这不是一个问题，但对于R，我不明白。

我试过：

mydata<-read.csv("mydata.csv",header=T,sep=";")
summary(aov(Measure~factor(Machine)/factor(Operator),data=mydata))

但F值不正确。

谢谢你的帮忙

PS:这是NIST (mydata.csv)的数据集

Operator;Machine;Replicate;Measure
1;1;1;0.125
1;1;2;0.127
1;1;3;0.125
1;1;4;0.126
1;1;5;0.128
1;2;1;0.118
1;2;2;0.122
1;2;3;0.12
1;2;4;0.124
1;2;5;0.119
1;3;1;0.123
1;3;2;0.125
1;3;3;0.125
1;3;4;0.124
1;3;5;0.126
1;4;1;0.126
1;4;2;0.128
1;4;3;0.126
1;4;4;0.127
1;4;5;0.129
1;5;1;0.118
1;5;2;0.129
1;5;3;0.127
1;5;4;0.12
1;5;5;0.121
2;1;1;0.124
2;1;2;0.128
2;1;3;0.127
2;1;4;0.126
2;1;5;0.129
2;2;1;0.116
2;2;2;0.125
2;2;3;0.119
2;2;4;0.125
2;2;5;0.12
2;3;1;0.122
2;3;2;0.121
2;3;3;0.124
2;3;4;0.126
2;3;5;0.125
2;4;1;0.126
2;4;2;0.129
2;4;3;0.125
2;4;4;0.13
2;4;5;0.124
2;5;1;0.125
2;5;2;0.123
2;5;3;0.114
2;5;4;0.124
2;5;5;0.117

Answer 1

您可以指定带有随机项的模型，以及如何使用公式中的Error项来划分F比(我将机器和操作符更改为数据中的因子)：

> summary(aov(Measure~Machine + Error(interaction(Machine,Operator)),data=mydat))

Error: interaction(Machine, Operator)
          Df    Sum Sq   Mean Sq F value  Pr(>F)   
Machine    4 0.0003033 7.583e-05   20.38 0.00269 **
Residuals  5 0.0000186 3.720e-06                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Within
          Df   Sum Sq  Mean Sq F value Pr(>F)
Residuals 40 0.000346 8.65e-06            

Answer 2

aov函数的摘要：

> summary(aov(Measure~factor(Machine)/factor(Operator),data=mydata))
                                 Df    Sum Sq   Mean Sq F value   Pr(>F)    
factor(Machine)                   4 0.0003033 7.583e-05   8.766 3.52e-05 ***
factor(Machine):factor(Operator)  5 0.0000186 3.720e-06   0.430    0.825    
Residuals                        40 0.0003460 8.650e-06                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

正如您在问题中所说的，8.766不同于您所期望的20.38。

但是，这不是函数的错误。这是因为aov函数是如何工作的。正如文档所述(?aov)：

通过调用lm对各层进行方差模型的拟合。

请注意上线中的“每一个”。这意味着aov函数总是根据variable MSE / Residuals MSE计算F统计量(F值)。如您所见，8.766是factor(Machine)和Residuals (即75.83 / 8.65 = 8.766 )之间划分的结果。

在源代码中也可以看到这一点，如果您想挖得更深一点(在控制台上键入summary.aov )，其中写着：

if (rdf > 0L) {
        TT <- ms/ms[nt]
        TP <- pf(TT, df, rdf, lower.tail = FALSE)
        TT[nt] <- TP[nt] <- NA
        x$"F value" <- TT
        x$"Pr(>F)" <- TP
    }

这也是为什么F value for factor(Machine):factor(Operator)是正确的。

如果您想要找到F-统计量和相关的概率，那么factor(Machine)和factor(Machine):factor(Operator)之间的除法只需计算如下：

a <- summary(aov(Measure~factor(Machine)/factor(Operator),data=mydata))

和

> a[[1]]$'Mean Sq'[1] / a[[1]]$'Mean Sq'[2] #which is what you want 
[1] 20.38441

对应概率Pr(>F)

> pf(20.38441, df1=4, df2=5, lower.tail = FALSE) 
[1] 0.00269263