ezANOVA R 检查误差正态分布
ezANOVA R check error normally distributed
我正在使用 ezANOVA 对具有主题内变量和主题间变量的实验设计进行分析。我成功实施了 ezANOVA，如下所示：
structure(list(Sub = structure(c(3L, 3L, 3L, 4L, 4L, 4L, 1L,
1L, 1L, 2L, 2L, 2L), .Label = c("A7011", "A7022", "B13", "B14"
), class = "factor"), Depvariable = c(0.375, 0.066667, 0.15,
0.275, 0.025, 0.78333, 0.24167, 0.058333, 0.14167, 0.19167, 0.5,
0), Group = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L), .Label = c("A", "B"), class = "factor"), WithinFactor = c(0.6,
0, -0.3, 0.6, 0, -0.3, 0.6, 0, -0.3, 0.6, 0, -0.3)), .Names = c("Sub",
"Depvariable", "Group", "WithinFactor"), row.names = c(NA, 12L
), class = "data.frame")
mod.ez<-ezANOVA(data,
dv = .(Depvariable),
wid = .(Sub), # subject
within = .(WithinFactor),
between=.(Group),
type=3,
detailed=TRUE,
return_aov=TRUE)
我坚持检查残差正态分布的程序。
我试过以下方法：
shapiro.test(as.numeric(residuals(mod.ez$aov)))
但是我得到以下错误
shapiro.test(as.numeric(残差(mod.ez$aov))) 错误：
样本量必须介于 3 和 5000
之间
如果我调用 residuals(mod.ez$aov) 结果为 NULL。
我或者使用 lmer 来检查残差
plot(fitted(model_lmer), residuals(model_lmer))
但是，由于 ezANOVA 也实施了球形度的测试和校正，我想坚持下去并找到一种方法来检查残差的正态性假设。
非常感谢任何帮助
步骤：
完整示例
首先，您的代码的完整版本是：
library(ez)
data <- structure(list(Sub = structure(c(3L, 3L, 3L, 4L, 4L, 4L, 1L,
1L, 1L, 2L, 2L, 2L), .Label = c("A7011", "A7022", "B13", "B14"
), class = "factor"), Depvariable = c(0.375, 0.066667, 0.15,
0.275, 0.025, 0.78333, 0.24167, 0.058333, 0.14167, 0.19167, 0.5,
0), Group = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L), .Label = c("A", "B"), class = "factor"), WithinFactor = c(0.6,
0, -0.3, 0.6, 0, -0.3, 0.6, 0, -0.3, 0.6, 0, -0.3)), .Names = c("Sub",
"Depvariable", "Group", "WithinFactor"), row.names = c(NA, 12L
), class = "data.frame")
mod.ez <- ezANOVA(
data,
dv = .(Depvariable),
wid = .(Sub), # subject
within = .(WithinFactor),
between = .(Group),
type = 3,
detailed = TRUE,
return_aov = TRUE)
如何探索复杂的 R 结构
其次，如果您找不到残差（等），问题是：ezANOVA 的结果是否真的包含它们，或者它是否丢弃了信息？对于这类问题，我喜欢使用这个函数：
wtf_is <- function(x) {
# For when you have no idea what something is.
#
cat("1. typeof():\n")
print(typeof(x))
cat("\n2. class():\n")
print(class(x))
cat("\n3. mode():\n")
print(mode(x))
cat("\n4. names():\n")
print(names(x))
cat("\n5. slotNames():\n")
print(slotNames(x))
cat("\n6. attributes():\n")
print(attributes(x))
cat("\n7. str():\n")
print(str(x))
}
因此：
wtf_is(mod.ez)
寻找 ezANOVA 输出中的残差
输出很长。我们正在寻找长度为 12 的列表（因为您有 12 个数据点），或者看起来像残差或预测值的东西。部分输出为：
[...]
7. str():
List of 2
$ ANOVA:'data.frame': 3 obs. of 9 variables:
[...]
$ aov :List of 4
..$ (Intercept) :List of 9
[...]
..$ Sub :List of 9
[...]
.. ..$ residuals : Named num [1:3] 0.102 -0.116 0.164
.. .. ..- attr(*, "names")= chr [1:3] "2" "3" "4"
[...]
.. ..$ fitted.values: Named num [1:3] -1.39e-17 1.28e-01 9.03e-02
.. .. ..- attr(*, "names")= chr [1:3] "2" "3" "4"
..$ Sub:WithinFactor:List of 9
[...]
.. ..$ residuals : Named num [1:4] 0.00964 0.00964 0.23081 -0.23081
.. .. ..- attr(*, "names")= chr [1:4] "5" "6" "7" "8"
[...]
.. ..$ fitted.values: Named num [1:4] 0.0804 -0.0804 -0.0444 -0.0444
.. .. ..- attr(*, "names")= chr [1:4] "5" "6" "7" "8"
[...]
..$ Within :List of 6
[...]
.. ..$ residuals : num [1:4, 1] 0.3286 0.1098 -0.4969 0.0564
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:4] "9" "10" "11" "12"
.. .. .. ..$ : NULL
.. ..$ fitted.values: num [1:4, 1] 0 0 0 0
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:4] "9" "10" "11" "12"
.. .. .. ..$ : NULL
[...]
..- attr(*, "error.qr")=List of 5
.. ..$ qr : num [1:12, 1:8] -3.464 0.289 0.289 0.289 0.289 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:12] "1" "2" "3" "4" ...
.. .. .. ..$ : chr [1:8] "(Intercept)" "Sub1" "Sub2" "Sub3" ...
.. .. ..- attr(*, "assign")= int [1:8] 0 1 1 1 2 2 2 2
.. .. ..- attr(*, "contrasts")=List of 1
.. .. .. ..$ Sub: chr "contr.helmert"
[...]
... none 其中看起来对我很有帮助。所以答案很可能是"it's not there"，或者"not obviously there"，其他人也同意：
改用afex::aov_ez
因此您可以改用：
library(afex)
model2 <- aov_ez(
id = "Sub", # subject
dv = "Depvariable",
data = data,
between = c("Group"),
within = c("WithinFactor"),
type = "III" # or 3; type III sums of squares
)
anova(model2)
summary(model2)
residuals(model2$lm)
...这确实会给你残差。
然而，它也给出了不同的 F/p 值。
注意为什么 aov_ez 和 ezANOVA 在这里给出不同的答案
我们有：
> mod.ez
$ANOVA
Effect DFn DFd SSn SSd F p p<.05 ges
1 Group 1 2 0.024449088 0.05070517 0.96436277 0.4296328 0.134418588
2 WithinFactor 1 2 0.001296481 0.10673345 0.02429382 0.8904503 0.008167579
3 Group:WithinFactor 1 2 0.015557350 0.10673345 0.29151781 0.6433264 0.089928978
> anova(model2)
Anova Table (Type III tests)
Response: Depvariable
num Df den Df MSE F ges Pr(>F)
Group 1.0000 2.0000 0.025353 0.9644 0.07197 0.4296
WithinFactor 1.4681 2.9363 0.090093 0.2322 0.08876 0.7471
Group:WithinFactor 1.4681 2.9363 0.090093 1.5001 0.38628 0.3370
不同的结果。注意来自 mod.ez:
的警告消息
Warning: "WithinFactor" will be treated as numeric
... 即作为连续预测变量（协变量），而不是离散预测变量（因子）。所以我们应该看看 covariate 和 factorize 参数；参见 ?aov_ez。我必须说，我费了一番功夫才弄清楚如何在这里进行受试者内 ANCOVA。如果我正确阅读文档，factorize 部分仅适用于主体间预测变量，同样，covariate 仅适用于主体间协变量。
作为快速检查，如果您使用 ezANOVA 并强制它使用 WithinFactor 作为离散（非连续）预测变量，如下所示：
data$WithinCovariate <- data$WithinFactor # so the name is clearer!
data$WithinFactorDiscrete <- as.factor(data$WithinFactor)
mod.ez.discrete <- ezANOVA(
data,
dv = .(Depvariable),
wid = .(Sub), # subject
within = .(WithinFactorDiscrete),
between = .(Group),
type = 3,
detailed = TRUE,
return_aov = TRUE)
... 你得到 F/p 值匹配 aov_ez:
> mod.ez.discrete
$ANOVA
Effect DFn DFd SSn SSd F p p<.05 ges
1 (Intercept) 1 2 0.65723113 0.05070517 25.9236350 0.03647725 * 0.67583504
2 Group 1 2 0.02444909 0.05070517 0.9643628 0.42963280 0.07197457
3 WithinFactorDiscrete 2 4 0.03070651 0.26453641 0.2321534 0.80280844 0.08876045
4 Group:WithinFactorDiscrete 2 4 0.19841198 0.26453641 1.5000731 0.32651697 0.38627588
这样你就可以得到匹配的结果，Greenhouse-Geisser/Huynh-Feldt 更正和残差，除了受试者内的协变量。
终于...
用连续的受试者内预测变量检查球形度是什么意思？我对此完全不清楚；球形性与主体内因素不同水平的值对之间差异的方差同质性有关。如果预测变量是连续的，则没有对。
所以冒着犯错的风险，我要么
(a) 信任 ezANOVA 并放弃残差；
(b) 使用除了球形测试之外可以做任何事情的东西，像这样：
library(lme4)
library(lmerTest) # upgrades reports from lme4 to include p values! ;)
mod.lmer.wscov_interact <- lmer(
Depvariable ~
Group * WithinCovariate
+ (1 | Sub),
data = data
)
anova(mod.lmer.wscov_interact)
residuals(mod.lmer.wscov_interact)
mod.lmer.wscov_no_interact <- lmer(
Depvariable ~
Group + WithinCovariate
+ (1 | Sub),
data = data
)
anova(mod.lmer.wscov_no_interact)
mod.lmer.wsfac <- lmer(
Depvariable ~
Group * WithinFactorDiscrete
+ (1 | Sub),
data = data
)
anova(mod.lmer.wsfac)
给予
> anova(mod.lmer.wscov_interact)
Analysis of Variance Table of type III with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
Group 0.033586 0.033586 1 8 0.50936 0.4957
WithinCovariate 0.001296 0.001296 1 8 0.01966 0.8920
Group:WithinCovariate 0.015557 0.015557 1 8 0.23594 0.6402
> residuals(mod.lmer.wscov_interact)
1 2 3 4 5 6 7 8 9 10 11 12
0.130059250 -0.219344250 -0.156546500 0.030059250 -0.261011250 0.476783500 -0.009225679 -0.118156464 0.002383643 -0.059225679 0.323510536 -0.139286357
> anova(mod.lmer.wscov_no_interact)
Analysis of Variance Table of type III with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
Group 0.0244491 0.0244491 1 9 0.40519 0.5403
WithinCovariate 0.0012965 0.0012965 1 9 0.02149 0.8867
> anova(mod.lmer.wsfac)
Analysis of Variance Table of type III with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
Group 0.024449 0.024449 1 6 0.46534 0.5206
WithinFactorDiscrete 0.030707 0.015353 2 6 0.29222 0.7567
Group:WithinFactorDiscrete 0.198412 0.099206 2 6 1.88819 0.2312
我正在使用 ezANOVA 对具有主题内变量和主题间变量的实验设计进行分析。我成功实施了 ezANOVA，如下所示：
structure(list(Sub = structure(c(3L, 3L, 3L, 4L, 4L, 4L, 1L,
1L, 1L, 2L, 2L, 2L), .Label = c("A7011", "A7022", "B13", "B14"
), class = "factor"), Depvariable = c(0.375, 0.066667, 0.15,
0.275, 0.025, 0.78333, 0.24167, 0.058333, 0.14167, 0.19167, 0.5,
0), Group = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L), .Label = c("A", "B"), class = "factor"), WithinFactor = c(0.6,
0, -0.3, 0.6, 0, -0.3, 0.6, 0, -0.3, 0.6, 0, -0.3)), .Names = c("Sub",
"Depvariable", "Group", "WithinFactor"), row.names = c(NA, 12L
), class = "data.frame")
mod.ez<-ezANOVA(data,
dv = .(Depvariable),
wid = .(Sub), # subject
within = .(WithinFactor),
between=.(Group),
type=3,
detailed=TRUE,
return_aov=TRUE)
我坚持检查残差正态分布的程序。 我试过以下方法：
shapiro.test(as.numeric(residuals(mod.ez$aov)))
但是我得到以下错误
shapiro.test(as.numeric(残差(mod.ez$aov))) 错误： 样本量必须介于 3 和 5000
之间如果我调用 residuals(mod.ez$aov) 结果为 NULL。
我或者使用 lmer 来检查残差
plot(fitted(model_lmer), residuals(model_lmer))
但是，由于 ezANOVA 也实施了球形度的测试和校正，我想坚持下去并找到一种方法来检查残差的正态性假设。
非常感谢任何帮助
步骤：
完整示例
首先，您的代码的完整版本是：
library(ez)
data <- structure(list(Sub = structure(c(3L, 3L, 3L, 4L, 4L, 4L, 1L,
1L, 1L, 2L, 2L, 2L), .Label = c("A7011", "A7022", "B13", "B14"
), class = "factor"), Depvariable = c(0.375, 0.066667, 0.15,
0.275, 0.025, 0.78333, 0.24167, 0.058333, 0.14167, 0.19167, 0.5,
0), Group = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L,
1L, 1L), .Label = c("A", "B"), class = "factor"), WithinFactor = c(0.6,
0, -0.3, 0.6, 0, -0.3, 0.6, 0, -0.3, 0.6, 0, -0.3)), .Names = c("Sub",
"Depvariable", "Group", "WithinFactor"), row.names = c(NA, 12L
), class = "data.frame")
mod.ez <- ezANOVA(
data,
dv = .(Depvariable),
wid = .(Sub), # subject
within = .(WithinFactor),
between = .(Group),
type = 3,
detailed = TRUE,
return_aov = TRUE)
如何探索复杂的 R 结构
其次，如果您找不到残差（等），问题是：ezANOVA 的结果是否真的包含它们，或者它是否丢弃了信息？对于这类问题，我喜欢使用这个函数：
wtf_is <- function(x) {
# For when you have no idea what something is.
#
cat("1. typeof():\n")
print(typeof(x))
cat("\n2. class():\n")
print(class(x))
cat("\n3. mode():\n")
print(mode(x))
cat("\n4. names():\n")
print(names(x))
cat("\n5. slotNames():\n")
print(slotNames(x))
cat("\n6. attributes():\n")
print(attributes(x))
cat("\n7. str():\n")
print(str(x))
}
因此：
wtf_is(mod.ez)
寻找 ezANOVA 输出中的残差
输出很长。我们正在寻找长度为 12 的列表（因为您有 12 个数据点），或者看起来像残差或预测值的东西。部分输出为：
[...]
7. str():
List of 2
$ ANOVA:'data.frame': 3 obs. of 9 variables:
[...]
$ aov :List of 4
..$ (Intercept) :List of 9
[...]
..$ Sub :List of 9
[...]
.. ..$ residuals : Named num [1:3] 0.102 -0.116 0.164
.. .. ..- attr(*, "names")= chr [1:3] "2" "3" "4"
[...]
.. ..$ fitted.values: Named num [1:3] -1.39e-17 1.28e-01 9.03e-02
.. .. ..- attr(*, "names")= chr [1:3] "2" "3" "4"
..$ Sub:WithinFactor:List of 9
[...]
.. ..$ residuals : Named num [1:4] 0.00964 0.00964 0.23081 -0.23081
.. .. ..- attr(*, "names")= chr [1:4] "5" "6" "7" "8"
[...]
.. ..$ fitted.values: Named num [1:4] 0.0804 -0.0804 -0.0444 -0.0444
.. .. ..- attr(*, "names")= chr [1:4] "5" "6" "7" "8"
[...]
..$ Within :List of 6
[...]
.. ..$ residuals : num [1:4, 1] 0.3286 0.1098 -0.4969 0.0564
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:4] "9" "10" "11" "12"
.. .. .. ..$ : NULL
.. ..$ fitted.values: num [1:4, 1] 0 0 0 0
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:4] "9" "10" "11" "12"
.. .. .. ..$ : NULL
[...]
..- attr(*, "error.qr")=List of 5
.. ..$ qr : num [1:12, 1:8] -3.464 0.289 0.289 0.289 0.289 ...
.. .. ..- attr(*, "dimnames")=List of 2
.. .. .. ..$ : chr [1:12] "1" "2" "3" "4" ...
.. .. .. ..$ : chr [1:8] "(Intercept)" "Sub1" "Sub2" "Sub3" ...
.. .. ..- attr(*, "assign")= int [1:8] 0 1 1 1 2 2 2 2
.. .. ..- attr(*, "contrasts")=List of 1
.. .. .. ..$ Sub: chr "contr.helmert"
[...]
... none 其中看起来对我很有帮助。所以答案很可能是"it's not there"，或者"not obviously there"，其他人也同意：
改用afex::aov_ez
因此您可以改用：
library(afex)
model2 <- aov_ez(
id = "Sub", # subject
dv = "Depvariable",
data = data,
between = c("Group"),
within = c("WithinFactor"),
type = "III" # or 3; type III sums of squares
)
anova(model2)
summary(model2)
residuals(model2$lm)
...这确实会给你残差。
然而，它也给出了不同的 F/p 值。
注意为什么 aov_ez 和 ezANOVA 在这里给出不同的答案
我们有：
> mod.ez
$ANOVA
Effect DFn DFd SSn SSd F p p<.05 ges
1 Group 1 2 0.024449088 0.05070517 0.96436277 0.4296328 0.134418588
2 WithinFactor 1 2 0.001296481 0.10673345 0.02429382 0.8904503 0.008167579
3 Group:WithinFactor 1 2 0.015557350 0.10673345 0.29151781 0.6433264 0.089928978
> anova(model2)
Anova Table (Type III tests)
Response: Depvariable
num Df den Df MSE F ges Pr(>F)
Group 1.0000 2.0000 0.025353 0.9644 0.07197 0.4296
WithinFactor 1.4681 2.9363 0.090093 0.2322 0.08876 0.7471
Group:WithinFactor 1.4681 2.9363 0.090093 1.5001 0.38628 0.3370
不同的结果。注意来自 mod.ez:
的警告消息Warning: "WithinFactor" will be treated as numeric
... 即作为连续预测变量（协变量），而不是离散预测变量（因子）。所以我们应该看看 covariate 和 factorize 参数；参见 ?aov_ez。我必须说，我费了一番功夫才弄清楚如何在这里进行受试者内 ANCOVA。如果我正确阅读文档，factorize 部分仅适用于主体间预测变量，同样，covariate 仅适用于主体间协变量。
作为快速检查，如果您使用 ezANOVA 并强制它使用 WithinFactor 作为离散（非连续）预测变量，如下所示：
data$WithinCovariate <- data$WithinFactor # so the name is clearer!
data$WithinFactorDiscrete <- as.factor(data$WithinFactor)
mod.ez.discrete <- ezANOVA(
data,
dv = .(Depvariable),
wid = .(Sub), # subject
within = .(WithinFactorDiscrete),
between = .(Group),
type = 3,
detailed = TRUE,
return_aov = TRUE)
... 你得到 F/p 值匹配 aov_ez:
> mod.ez.discrete
$ANOVA
Effect DFn DFd SSn SSd F p p<.05 ges
1 (Intercept) 1 2 0.65723113 0.05070517 25.9236350 0.03647725 * 0.67583504
2 Group 1 2 0.02444909 0.05070517 0.9643628 0.42963280 0.07197457
3 WithinFactorDiscrete 2 4 0.03070651 0.26453641 0.2321534 0.80280844 0.08876045
4 Group:WithinFactorDiscrete 2 4 0.19841198 0.26453641 1.5000731 0.32651697 0.38627588
这样你就可以得到匹配的结果，Greenhouse-Geisser/Huynh-Feldt 更正和残差，除了受试者内的协变量。
终于...
用连续的受试者内预测变量检查球形度是什么意思？我对此完全不清楚；球形性与主体内因素不同水平的值对之间差异的方差同质性有关。如果预测变量是连续的，则没有对。
所以冒着犯错的风险，我要么 (a) 信任 ezANOVA 并放弃残差； (b) 使用除了球形测试之外可以做任何事情的东西，像这样：
library(lme4)
library(lmerTest) # upgrades reports from lme4 to include p values! ;)
mod.lmer.wscov_interact <- lmer(
Depvariable ~
Group * WithinCovariate
+ (1 | Sub),
data = data
)
anova(mod.lmer.wscov_interact)
residuals(mod.lmer.wscov_interact)
mod.lmer.wscov_no_interact <- lmer(
Depvariable ~
Group + WithinCovariate
+ (1 | Sub),
data = data
)
anova(mod.lmer.wscov_no_interact)
mod.lmer.wsfac <- lmer(
Depvariable ~
Group * WithinFactorDiscrete
+ (1 | Sub),
data = data
)
anova(mod.lmer.wsfac)
给予
> anova(mod.lmer.wscov_interact)
Analysis of Variance Table of type III with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
Group 0.033586 0.033586 1 8 0.50936 0.4957
WithinCovariate 0.001296 0.001296 1 8 0.01966 0.8920
Group:WithinCovariate 0.015557 0.015557 1 8 0.23594 0.6402
> residuals(mod.lmer.wscov_interact)
1 2 3 4 5 6 7 8 9 10 11 12
0.130059250 -0.219344250 -0.156546500 0.030059250 -0.261011250 0.476783500 -0.009225679 -0.118156464 0.002383643 -0.059225679 0.323510536 -0.139286357
> anova(mod.lmer.wscov_no_interact)
Analysis of Variance Table of type III with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
Group 0.0244491 0.0244491 1 9 0.40519 0.5403
WithinCovariate 0.0012965 0.0012965 1 9 0.02149 0.8867
> anova(mod.lmer.wsfac)
Analysis of Variance Table of type III with Satterthwaite
approximation for degrees of freedom
Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
Group 0.024449 0.024449 1 6 0.46534 0.5206
WithinFactorDiscrete 0.030707 0.015353 2 6 0.29222 0.7567
Group:WithinFactorDiscrete 0.198412 0.099206 2 6 1.88819 0.2312