Повторная выборка ANOVA с коррекцией Welch

Я полагаю, что вы правы в том, что не существует простого способа сделать аномалию с поправкой на валлий с повторной выборкой, но должна быть возможность совместить несколько вещей, чтобы это заработало.

require('Ecdat')

Я буду использовать набор данных «Звезда» из пакета «Экдат», в котором рассматривается влияние небольшого размера класса на результаты стандартизированных тестов.

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 показывает, что остатки ANOVA отклоняются от нормальной qqline

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

Подобранная Y по сравнению с остатками показывает сходящуюся тенденцию в остатках, можно проверить с помощью теста Шапиро-Уилка, чтобы быть уверенным

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

Просто подтверждает, что мы не сможем предположить нормальное распределение.

Мы можем использовать начальную загрузку для оценки истинного F-расстояния.

#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() по умолчанию применяет поправку Велча, мы можем использовать ее для наших смоделированных данных.

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