R 中分层样本的单向方差分析

我有一个包含三组（“a”、“b”、“c”）的分层样本，这些样本是从较大的总体 N 中抽取的。所有组都有 30 个观察值，但它们在 N 中的比例不相等，因此它们的采样权重不同。

我用surveyR 中的包来计算汇总统计数据和线性回归模型，并且想知道如何计算单向方差分析来校正调查设计（如果需要）。

我的假设是，如果我错了，请纠正我，对于权重较小的总体，方差的标准误差通常应该更高，因此不考虑调查设计的简单方差分析应该不可靠。

这是一个例子。任何帮助，将不胜感激。

## Oneway- ANOVA tests in R for surveys with stratified sampling-design
library("survey")
# create test data
test.df<-data.frame(
  id=1:90,
  variable=c(rnorm(n = 30,mean=150,sd=10),
             rnorm(n = 30,mean=150,sd=10),
             rnorm(n = 30,mean=140,sd=10)),
  groups=c(rep("a",30),
  rep("b",30),
  rep("c",30)),
  weights=c(rep(1,30), # undersampled
  rep(1,30),
  rep(100,30))) # oversampled


# correct for survey design
test.df.survey<-svydesign(id=~id,
                           strata=~groups,
                           weights=~weights,
                           data=test.df)

## descriptive statistics
# boxplot
svyboxplot(~variable~groups,test.df.survey)
# means
svyby(~variable,~groups,test.df.survey,svymean)
# variances
svyby(~variable,~groups,test.df.survey,svyvar)


### ANOVA ###
## One-way ANOVA without correcting for survey design
summary(aov(formula = variable~groups,data = test.df))

嗯，这是一个有趣的问题，据我所知，很难在单向方差分析中考虑权重。因此，我决定向您展示解决这个问题的方法。

我将使用双向方差分析，然后使用 soem port hoc 测试。

首先，让我们根据您的数据构建一个线性模型并检查它的外观。

library(car)
library(agricolae)
model.lm = lm(variable ~ groups * weights, data = test.df)
shapiro.test(resid(model.lm))

Shapiro-Wilk normality test

data:  resid(model.lm)
W = 0.98238, p-value = 0.263

leveneTest(variable ~ groups * factor(weights), data = test.df)
Levene's Test for Homogeneity of Variance (center = median)
Df F value  Pr(>F)  
group  2  2.6422 0.07692 .
      87                  
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

分布接近正态，组间方差不同，因此方差不是同质的 - 应该用于参数检验 - 方差分析。不管怎样，我们还是进行一下测试吧。

检查我们的数据是否适合此测试的几个图：

hist(resid(model.lm))
plot(model.lm)

Here https://stats.stackexchange.com/questions/58141/interpreting-plot-lm/65864是情节的解释，它们实际上看起来并不糟糕。

让我们运行双向方差分析：

anova(model.lm)
Analysis of Variance Table

Response: variable
          Df Sum Sq Mean Sq F value    Pr(>F)    
groups     2 2267.8 1133.88  9.9566 0.0001277 ***
Residuals 87 9907.8  113.88                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

如您所见，结果非常接近您的结果。一些事后测试：

(result.hsd = HSD.test(model.lm, list('groups', 'weights')))
$statistics
   MSerror Df     Mean     CV      MSD
  113.8831 87 147.8164 7.2195 6.570186

$parameters
   test         name.t ntr StudentizedRange alpha
  Tukey groups:weights   3         3.372163  0.05

$means
      variable       std  r      Min      Max      Q25      Q50      Q75
a:1   150.8601 11.571185 30 113.3240 173.0429 145.2710 151.9689 157.8051
b:1   151.8486  8.330029 30 137.1907 176.9833 147.8404 150.3161 154.7321
c:100 140.7404 11.762979 30 118.0823 163.9753 131.6112 141.1810 147.8231

$comparison
NULL

$groups
      variable groups
b:1   151.8486      a
a:1   150.8601      a
c:100 140.7404      b

attr(,"class")
[1] "group"

也许还有一些不同的方式：

aov_cont<- aov(test.df$variable ~ test.df$groups * test.df$weights)
summary(aov_cont)
               Df Sum Sq Mean Sq F value   Pr(>F)    
test.df$groups  2   2268  1133.9   9.957 0.000128 ***
Residuals      87   9908   113.9                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(TukeyHSD(aov_cont))
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = test.df$variable ~ test.df$groups * test.df$weights)

$`test.df$groups`
           diff        lwr       upr     p adj
b-a   0.9884608  -5.581725  7.558647 0.9315792
c-a -10.1197048 -16.689891 -3.549519 0.0011934
c-b -11.1081657 -17.678352 -4.537980 0.0003461

总结一下，结果非常接近您的结果。就我个人而言，我将运行双向方差分析(*)符号或(+)当您确定变量是独立的时 - 加性模型。

Group c体重较大的群体不同a and b基本上。