我想推荐广义最小二乘法作为一种选择。同样，您可以为给定预测变量的数据均值创建模型，也可以为方差创建模型。nlme 包中的gls()函数允许我们这样做。

我将在下面演示如何：

set.seed(1984)
n <- 25
# Create heteroskedastic data, balanced design
dat <- data.frame(
  y = c(
    rnorm(n, 1, 1), rnorm(n, 0, 1), rnorm(n, .5, 1), 
    rnorm(n, 0, 3), rnorm(n, -.5, 4)),
  groups = sort(rep(LETTERS[1:5], n))
)
boxplot(y ~ groups, dat)

# Assuming all cases have identical variances
summary(aov(y ~ groups, dat))

             Df Sum Sq Mean Sq F value Pr(>F)  
groups        4   69.0  17.255   3.475 0.0101 *
Residuals   120  595.9   4.966                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Modeling the error variances
library(nlme)
# This is your standard linear model
coef(summary(fit.gls.baseline <- gls(y ~ groups, dat))))

                 Value Std.Error   t-value      p-value
(Intercept)  1.1784414 0.4457005  2.644021 0.0092875812
groupsB     -0.9748277 0.6303157 -1.546571 0.1246000669
groupsC     -0.8980932 0.6303157 -1.424831 0.1568017302
groupsD     -1.2903790 0.6303157 -2.047195 0.0428223361
groupsE     -2.3076197 0.6303157 -3.661054 0.0003747902

sigma(fit.gls.baseline) # get residual standard deviation
[1] 2.228502

# Next we fit a heteroskedastic model, hidden some output
# varIdent says the variances are identical for cases that have the same value of group
summary(fit.gls <- gls(
  y ~ groups, dat, weights = varIdent(form = ~ 1 | groups)))

Generalized least squares fit by REML
  Model: y ~ groups 
  Data: dat 
       AIC      BIC    logLik
  479.5354 507.4103 -229.7677

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | groups 
 Parameter estimates:
        A         B         C         D         E 
1.0000000 0.9741804 0.8698579 3.0062743 3.7867899 

Coefficients:
                 Value Std.Error   t-value p-value
(Intercept)  1.1784414 0.1951410  6.038923  0.0000
groupsB     -0.9748277 0.2724316 -3.578248  0.0005
groupsC     -0.8980932 0.2586375 -3.472401  0.0007
groupsD     -1.2903790 0.6182517 -2.087142  0.0390
groupsE     -2.3076197 0.7642898 -3.019299  0.0031

Residual standard error: 0.975705 

关于这个模型的一些注意事项。您会发现，与早期模型的相同估计值相比，异方差模型的残差标准偏差（误差）要小得多。现在看一下方差函数部分，您会看到 A 组的标准差为 1。这些值是相对标准差。因此，模型报告的 .976 是 A 组的残差是 B 组的值和 E × .976) 组的值要大得多残差标准差。$.976$$.974 \times .976$$(3 \times .976)$$(3.79 \times .976)$

另外，让我们获得F检验的异方差版本：

# marginal for type III sum of squares, only makes a difference for the
# test of the grouping variable if we have more than one grouping variable
anova(fit.gls, type = "marginal")
Denom. DF: 120 
            numDF  F-value p-value
(Intercept)     1 36.46859  <.0001
groups          4  5.51477   4e-04

足够的证据表明并非所有组都具有相同的平均值，。$F(4, 120)=5.51, p < .001$

现在将异方差模型与标准模型进行比较，注意系数大致相同，但是，除 D 组和 E 组外，异方差模型的标准误差要小得多。所以我们的p值也更小。模型比较的全局方法是：

anova(fit.gls, fit.gls.baseline) # Model fit improved
                 Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fit.gls              1 10 479.5354 507.4103 -229.7677                        
fit.gls.baseline     2  6 560.9588 577.6837 -274.4794 1 vs 2 89.42343  <.0001

这些结果表明异方差模型更好，AIC、BIC和统计学上显着的似然比检验更低。

现在，我们实际上无法访问数据生成过程，但我们总是可以绘制上面的 boxplt。假设我们想选择一个更简单的模型，我们可以假设 A、B 和 C 组具有相同的方差，而 D 和 E 组不同。然后我们可以尝试：

# Create a variable that classifies the groups according to variance
dat$var.groups <- (dat$groups %in% c("D", "E")) + 0
# then run:
summary(fit.gls.parsim <- gls(
  y ~ groups, dat, weights = varIdent(form = ~ 1 | var.groups)))

Generalized least squares fit by REML
  Model: y ~ groups 
  Data: dat 
       AIC      BIC    logLik
  475.3162 494.8287 -230.6581

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | var.groups 
 Parameter estimates:
       0        1 
1.000000 3.600023 

Coefficients:
                 Value Std.Error   t-value p-value
(Intercept)  1.1784414 0.1853218  6.358893  0.0000
groupsB     -0.9748277 0.2620846 -3.719516  0.0003
groupsC     -0.8980932 0.2620846 -3.426730  0.0008
groupsD     -1.2903790 0.6924235 -1.863569  0.0648
groupsE     -2.3076197 0.6924235 -3.332671  0.0011

# A model comparison between our previous best and this new model
anova(fit.gls, fit.gls.parsim)
               Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fit.gls            1 10 479.5354 507.4103 -229.7677                        
fit.gls.parsim     2  7 475.3162 494.8287 -230.6581 1 vs 2 1.780859  0.6191

似然比检验无法区分这个参数较少的更简单模型和我们早期的异方差模型。AIC 和 BIC 也倾向于这种模式。所以我们可能会选择这种模式，因为我们不知道真相。理论也可能暗示有关群体差异差异的某些事情。可以看出，异方差这一特定问题本身如何具有实质性的意义。

F -test 再次将是：

anova(fit.gls.parsim, type = "marginal")
Denom. DF: 120 
            numDF  F-value p-value
(Intercept)     1 40.43552  <.0001
groups          4  6.04201   2e-04

现在，下一件事可能是对比。我在工作中不太关心这个，所以我在这里做最简单的事情。可能有更好的方法来处理对比，请参阅 emmeans 包文档以获取更多选项。我只使用默认值。

library(emmeans)
pairs(emmeans(fit.gls.parsim, "groups"))
 contrast   estimate        SE  df t.ratio p.value
 A - B     0.9748277 0.2620846 120   3.720  0.0028
 A - C     0.8980932 0.2620846 120   3.427  0.0073
 A - D     1.2903790 0.6924235 120   1.864  0.3427
 A - E     2.3076197 0.6924235 120   3.333  0.0099
 B - C    -0.0767345 0.2620846 120  -0.293  0.9984
 B - D     0.3155513 0.6924235 120   0.456  0.9910
 B - E     1.3327920 0.6924235 120   1.925  0.3099
 C - D     0.3922858 0.6924235 120   0.567  0.9796
 C - E     1.4095265 0.6924235 120   2.036  0.2555
 D - E     1.0172407 0.9435106 120   1.078  0.8174

P value adjustment: tukey method for comparing a family of 5 estimates 

这些对比使用来自异方差模型的信息，因为如果您使用该模型，结果会有所不同fit.gls.baseline。

我希望上面的演示展示了在执行方差分析时如何使用广义最小二乘来解释异方差性。与许多其他方法相比，这种方法并不简单地将异方差视为一种麻烦。也可以包括控制变量，如 ANCOVA。

顺便说一句，如果有人还怀疑 ANCOVA 中的控制变量也存在异方差，那么glmmTMB()glmmTMB 包中的函数可以处理这种情况。您只需使用参数指定方差的回归方程dispformula =。

几种可能的方法之一是oneway.test在 R 中使用。这是一个包含三组的示例，每组有十个观察值：

x1 = rnorm(10, 100, 10);  x2 = rnorm(10, 95, 15);  x3 = rnorm(10, 90, 5)
x = c(x1, x2, x3); group = rep(1:3, each=10)
boxplot(x ~ group)

你可以看到我的假数据在三组中的每一组中都有不同的标准差。箱线图显示了这种异方差性。Bartlett 检验确认显着性（P 值 0.004）。

sd(x1); sd(x2); sd(x3)
[1] 10.88923
[1] 16.35099
[1] 4.656118

bartlett.test(x, group)

        Bartlett test of homogeneity of variances

data:  x and group
Bartlett's K-squared = 11.103, df = 2, p-value = 0.003881

该oneway.test过程允许不同的方差，其方式与 Welch 2 样本 t 检验的方式有些相同。它表明并非所有组群均值都相等（P 值 < 5%）。请注意，分母 DF假设方差相等的标准方差分析的分母 DF [我相信这是@SalMangiafico 建议的“韦尔奇方差分析”。]$\approx 15;$$= 27.$

mean(x1); mean(x2); mean(x3)
[1] 98.00458
[1] 105.3806
[1] 92.0077

oneway.test(x ~ group)

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

data:  x and group
F = 3.818, num df = 2.000, denom df = 14.536, p-value = 0.04649

您可以使用 Welch 2 样本 t 检验来探索配对比较，可能使用 Bonferroni 家族错误率。

参考： 这个问答提到了多种替代方法。