带有 Welch 校正的 ANOVA 重采样

Question

我正在使用相同的数据进行一些探索，并试图突出显示组内方差与组间方差。现在我已经能够成功地显示组间方差非常强，但是，数据的性质应该显示组内方差较弱。 （即我的 Shapiro-Wilk 正态性检验显示了这一点）我相信如果我使用 welch 校正进行一些重新采样，可能就是这种情况。

我想知道是否有人知道在 R 中是否存在基于重新采样的带有 Welch 校正的方差分析。我看到置换测试的 R 实现但没有校正。如果不是，我将如何在使用此实现时直接编写测试代码。 http://finzi.psych.upenn.edu/library/lmPerm/html/aovp.html

这是我的基本组间方差分析的大纲：

fit <- lm(formula = data$Boys ~ data$GroupofBoys)
anova(fit)

Answer 1

我相信您是正确的，因为没有一种简单的方法可以通过重新采样来进行 welch 校正 anova，但是应该可以将一些东西放在一起以使其工作。

require('Ecdat')

我将使用“Ecdat”包中的“Star”数据集，该数据集着眼于小班教学对标准化考试成绩的影响。

star<-Star
attach(star)

head(star)

        tmathssk treadssk  classk      totexpk sex  freelunk race  schidkn
2       473      447       small.class       7 girl       no white      63
3       536      450       small.class      21 girl       no black      20
5       463      439 regular.with.aide       0  boy      yes black      19
11      559      448           regular      16  boy       no white      69
12      489      447       small.class       5  boy      yes white      79
13      454      431           regular       8  boy      yes white       5

一些探索性分析：

#bloxplots 
boxplot(treadssk ~ classk, ylab="Total Reading Scaled Score")
title("Reading Scores by Class Size")

#histograms
hist(treadssk, xlab="Total Reading Scaled Score")

运行常规方差分析

model1 = aov(treadssk ~ classk, data = star)
summary(model1)

              Df  Sum Sq Mean Sq F value   Pr(>F)    
classk         2   37201   18601   18.54 9.44e-09 ***
Residuals   5745 5764478    1003                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

查看方差分析残差

#qqplot
qqnorm(residuals(model1),ylab="Reading Scaled Score")
qqline(residuals(model1),ylab="Reading Scaled Score")

qqplot 显示方差分析残差偏离正常 qqline

#Fitted Y vs. Residuals
plot(fitted(model1), residuals(model1))

拟合的 Y 与残差显示残差的收敛趋势，可以通过 Shapiro-Wilk 检验进行测试以确保

shapiro.test(treadssk[1:5000]) #shapiro.test contrained to sample sizes between 3 and 5000

Shapiro-Wilk normality test

data:  treadssk[1:5000]
W = 0.92256, p-value < 2.2e-16

只是确认我们无法假设正态分布。

我们可以使用 bootstrap 来估计真实的 F-dist。

#Bootstrap version (with 10,000 iterations)
mean_read = mean(treadssk)
grpA = treadssk[classk=="regular"] - mean_read[1]
grpB = treadssk[classk=="small.class"]  - mean_read[2]
grpC = treadssk[classk=="regular.with.aide"]  - mean_read[3]
sim_classk <- classk
R = 10000
sim_Fstar = numeric(R)
for (i in 1:R) {
  groupA = sample(grpA, size=2000, replace=T)
  groupB = sample(grpB, size=1733, replace=T)
  groupC = sample(grpC, size=2015, replace=T)
  sim_score = c(groupA,groupB,groupC)
  sim_data = data.frame(sim_score,sim_classk)
}

现在我们需要得到组因子的唯一对集合

allPairs <- expand.grid(levels(sim_data$sim_classk), levels(sim_data$sim_classk))
## http://stackoverflow.com/questions/28574006/unique-combination-of-two-columns-in-r/28574136#28574136
allPairs <- unique(t(apply(allPairs, 1, sort)))
allPairs <- allPairs[ allPairs[,1] != allPairs[,2], ]

allPairs
     [,1]                [,2]               
[1,] "regular"           "small.class"      
[2,] "regular"           "regular.with.aide"
[3,] "regular.with.aide" "small.class" 

由于 oneway.test() 默认应用 Welch 校正，我们可以在模拟数据上使用它。

allResults <- apply(allPairs, 1, function(p) {
#http://stackoverflow.com/questions/28587498/post-hoc-tests-for-one-way-anova-with-welchs-correction-in-r
  dat <- sim_data[sim_data$sim_classk %in% p, ]
  ret <- oneway.test(sim_score ~ sim_classk, data = sim_data, na.action = na.omit)
  ret$sim_classk <- p
  ret
})

length(allResults)
[1] 3


allResults[[1]]

One-way analysis of means (not assuming equal variances)

data:  sim_score and sim_classk
F = 1.7741, num df = 2.0, denom df = 1305.9, p-value = 0.170