完全嵌套方差分析

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

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):

Fit an analysis of variance model by a call to lm for each stratum.

注意上面一行中的'each'。这意味着 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 对应 factor(Machine):factor(Operator) 的原因。

如果您想找到 factor(Machine) 和 factor(Machine):factor(Operator) 之间除法的 F 统计量和相关概率，只需自行计算即可：

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

Answer 2

您可以使用公式中的 Error 项指定具有随机项的模型以及如何划分以获得 F 比率（我将 Machine 和 Operator 更改为数据中的因子）：

> 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

完全嵌套方差分析

Fully Nested ANOVA

r

statistics

nested

anova