带韦尔奇校正的方差分析重采样

Question

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

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

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

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

Answer 1

我相信你是对的，因为没有一种简单的方法可以通过重采样来进行韦尔奇校正方差分析，但应该可以将一些东西放在一起使其发挥作用。

require('Ecdat')

我将使用“Ecdat”包中的“Star”数据集，该数据集研究小 class 尺寸对标准化考试成绩的影响。

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))
## 
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) {
#
  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

带韦尔奇校正的方差分析重采样

ANOVA Resampling with Welch Correction

r

permutation

resampling

anova