机器算法验证 - 异方差数据的单向方差分析的替代方案 - 吾爱随笔录

异方差数据的单向方差分析的替代方案

机器算法验证 r 方差分析数据转换异方差

2022-01-30 10:24:33

我有来自 3 组藻类生物量（、、）的数据，其中包含不同的样本量（、、），我想比较这些组是否来自同一种群。 $A$ $B$ $C$ $n_A=15$ $n_B=13$ $n_C=12$

单向方差分析肯定是要走的路，但是在对我的数据进行正态性测试时，异方差似乎是主要问题。我的原始数据，未经任何转换，产生的方差比（）远高于临界值（），因此我无法执行单向方差分析. $F_{\max} = 19.1$ $F_{\rm crit} = 4.16$

我还尝试转换以规范化我的数据。即使在尝试了各种变换（对数、平方根、平方）之后，使用变换后产生的最低相比仍然更高。 $F_{\max}$ $\log_{10}$ $7.16$ $F_{\rm crit}$

这里的任何人都可以告诉我从这里去哪里吗？我想不出其他转换方法来按数据进行标准化。是否有任何替代单向方差分析的方法？

PS：我的原始数据如下：

A: 0.178 0.195 0.225 0.294 0.315 0.341 0.36  0.363 0.371 0.398 0.407 0.409 0.432 
   0.494 0.719
B: 0.11  0.111 0.204 0.416 0.417 0.441 0.492 0.965 1.113 1.19  1.233 1.505 1.897
C: 0.106 0.114 0.143 0.435 0.448 0.51  0.576 0.588 0.608 0.64  0.658 0.788 0.958

1个回答

处理异方差数据时有许多可用的选项。不幸的是，它们都不能保证始终有效。以下是我熟悉的一些选项：

转变
韦尔奇方差分析
加权最小二乘
稳健回归
异方差一致标准误
引导程序
Kruskal-Wallis 检验
序数逻辑回归

更新：当您具有方差异方差/异质性时，这里演示 R 了一些拟合线性模型（即方差分析或回归）的方法。

让我们先看看您的数据。为方便起见，我将它们加载到两个数据框中，称为my.data（其结构类似于上面，每组一列）和stacked.data（有两列：values数字和ind组指示符）。

我们可以用 Levene 的检验正式检验异方差性：

library(car)
leveneTest(values~ind, stacked.data)
# Levene's Test for Homogeneity of Variance (center = median)
#       Df F value   Pr(>F)   
# group  2  8.1269 0.001153 **
#       38

果然，你有异方差性。我们将检查组的方差是什么。一个经验法则是，只要最大方差不超过，线性模型对方差的异质性就相当稳健，因此我们也会找到该比率： $4\!\times$

apply(my.data, 2, function(x){ var(x, na.rm=T) })
#          A          B          C 
# 0.01734578 0.33182844 0.06673060 
var(my.data$B, na.rm=T) / var(my.data$A, na.rm=T)
# [1] 19.13021

您的方差差异很大，最大的B是最小的。这是一个有问题的异方差水平。 $19\!\times$ A

您曾考虑使用对数或平方根等转换来稳定方差。这在某些情况下会起作用，但是Box-Cox类型的变换通过不对称地压缩数据来稳定方差，要么向下压缩它们，同时压缩最多的最高数据，要么向上压缩它们，同时压缩最多的最低数据。因此，您需要数据的方差随着平均值的变化而变化，以使其发挥最佳作用。您的数据的方差差异很大，但均值和中位数之间的差异相对较小，即分布大多重叠。作为教学练习，我们可以parallel.universe.data通过将添加到所有值和来创建一些 $2.7$ B $.7$ C来展示它是如何工作的：

parallel.universe.data = with(my.data, data.frame(A=A, B=B+2.7, C=C+.7))
apply(parallel.universe.data,       2, function(x){ var(x, na.rm=T) })
#          A          B          C 
# 0.01734578 0.33182844 0.06673060 
apply(log(parallel.universe.data),  2, function(x){ var(x, na.rm=T) })
#          A          B          C 
# 0.12750634 0.02631383 0.05240742 
apply(sqrt(parallel.universe.data), 2, function(x){ var(x, na.rm=T) })
#          A          B          C 
# 0.01120956 0.02325107 0.01461479 
var(sqrt(parallel.universe.data$B), na.rm=T) / 
var(sqrt(parallel.universe.data$A), na.rm=T)
# [1] 2.074217

使用平方根变换可以很好地稳定这些数据。您可以在这里看到平行宇宙数据的改进：

在此处输入图像描述

除了尝试不同的变换，更系统的方法是优化 Box-Cox 参数（尽管通常建议将其四舍五入到最接近的可解释变换）。在您的情况下，平方根或对数都是可以接受的，尽管它们实际上都不起作用。对于平行宇宙数据，平方根是最好的： $\lambda$ $\lambda = .5$ $\lambda = 0$

boxcox(values~ind, data=stacked.data,    na.action=na.omit)
boxcox(values~ind, data=stacked.pu.data, na.action=na.omit)

在此处输入图像描述

由于这种情况是 ANOVA（即没有连续变量），处理异质性的一种方法是使用Welch 校正检验中的分母自由度（nb, , 一个小数值，而不是）： $F$ df = 19.445df = 38

oneway.test(values~ind, data=stacked.data, na.action=na.omit, var.equal=FALSE)
#  One-way analysis of means (not assuming equal variances)
# 
# data:  values and ind
# F = 4.1769, num df = 2.000, denom df = 19.445, p-value = 0.03097

更通用的方法是使用加权最小二乘法。由于某些组 ( B) 分散得更多，因此这些组中的数据提供的关于均值位置的信息比其他组中的数据少。我们可以通过为每个数据点提供权重来让模型结合这一点。一个常见的系统是使用组方差的倒数作为权重：

wl = 1 / apply(my.data, 2, function(x){ var(x, na.rm=T) })
stacked.data$w = with(stacked.data, ifelse(ind=="A", wl[1], 
                                           ifelse(ind=="B", wl[2], wl[3])))
w.mod = lm(values~ind, stacked.data, na.action=na.omit, weights=w)
anova(w.mod)
# Response: values
#           Df Sum Sq Mean Sq F value  Pr(>F)  
# ind        2   8.64  4.3201  4.3201 0.02039 *
# Residuals 38  38.00  1.0000

和值与未加权 ANOVA ( , )略有不同，但它很好地解决了异质性问题： $F$ $p$ 4.50890.01749

在此处输入图像描述

然而，加权最小二乘法并不是灵丹妙药。一个令人不安的事实是，只有权重恰到好处，这意味着它们是先验已知的。它也没有解决非正态性（例如偏斜）或异常值。但是，使用根据您的数据估计的权重通常可以正常工作，特别是如果您有足够的数据来以合理的精度估计方差（这类似于使用表而不是表，当您有或 $z$ $t$ $50$ $100$ 自由度），您的数据足够正常，并且您似乎没有任何异常值。不幸的是，您的数据相对较少（每组 13 或 15 个），有些偏差，可能还有一些异常值。我不确定这些是否足够糟糕，是否可以大肆宣传，但您可以将加权最小二乘法与稳健的方法混合使用。您可以使用四分位间距的倒数（不受每组中高达 50% 的异常值的影响），而不是使用方差作为散布的度量（这对异常值很敏感，尤其是在低然后可以使用不同的损失函数（如 Tukey 的双方）将这些权重与稳健回归相结合： $N$

1 / apply(my.data,       2, function(x){ var(x, na.rm=T) })
#         A         B         C 
# 57.650907  3.013606 14.985628 
1 / apply(my.data,       2, function(x){ IQR(x, na.rm=T) })
#        A        B        C 
# 9.661836 1.291990 4.878049 

rw = 1 / apply(my.data, 2, function(x){ IQR(x, na.rm=T) })
stacked.data$rw = with(stacked.data, ifelse(ind=="A", rw[1], 
                                            ifelse(ind=="B", rw[2], rw[3])))
library(robustbase)
w.r.mod = lmrob(values~ind, stacked.data, na.action=na.omit, weights=rw)
anova(w.r.mod, lmrob(values~1,stacked.data,na.action=na.omit,weights=rw), test="Wald")
# Robust Wald Test Table
# 
# Model 1: values ~ ind
# Model 2: values ~ 1
# Largest model fitted by lmrob(), i.e. SM
# 
#   pseudoDf Test.Stat Df Pr(>chisq)  
# 1       38                          
# 2       40    6.6016  2    0.03685 *

这里的重量没有那么极端。预测的组均值略有不同（A：WLS 0.36673，稳健0.35722；B：WLS 0.77646，稳健0.70433；C：WLS 0.50554，稳健0.51845），均值B和C受极值影响较小。

在计量经济学中，Huber-White（“三明治”）标准误非常流行。与 Welch 校正一样，这不需要您先验地知道方差，也不需要您根据数据和/或可能不正确的模型来估计权重。另一方面，我不知道如何将其与 ANOVA 结合，这意味着您只能将它们用于单个虚拟代码的测试，这让我觉得在这种情况下不太有用，但无论如何我都会展示它们：

library(sandwich)
mod = lm(values~ind, stacked.data, na.action=na.omit)
sqrt(diag(vcovHC(mod)))
# (Intercept)        indB        indC 
#  0.03519921  0.16997457  0.08246131 
2*(1-pt(coef(mod) / sqrt(diag(vcovHC(mod))), df=38))
#  (Intercept)         indB         indC 
# 1.078249e-12 2.087484e-02 1.005212e-01

该函数vcovHC为您的 beta（您的虚拟代码）计算异方差一致的方差-协方差矩阵，这就是函数调用中的字母所代表的含义。要获得标准误差，您需要提取主对角线并取平方根。要为您的 beta 获得检验，您将系数估计值除以 SE，并将结果与适当的分布（即，具有剩余自由度 $t$ $t$ $t$

特别是对于R用户，@TomWenseleers 在下面的评论中指出，包中的?Anova函数car可以接受一个white.adjust参数，以使用异方差一致错误获取因子的值。 $p$

Anova(mod, white.adjust=TRUE)
# Analysis of Deviance Table (Type II tests)
# 
# Response: values
#           Df      F  Pr(>F)  
# ind        2 3.9946 0.02663 *
# Residuals 38                 
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

您可以尝试通过bootstrapping对测试统计的实际抽样分布进行经验估计。首先，通过使所有组均值完全相等来创建一个真正的空值。然后，您使用替换重新采样并计算每个引导样本的检验统计量 (在空值下的采样分布的经验估计，无论它们在正态性或同质性方面的状态如何。与您观察到的检验统计量一样极端或更极端的抽样分布的比例是值： $F$ $F$ $p$

mod    = lm(values~ind, stacked.data, na.action=na.omit)
F.stat = anova(mod)[1,4]
  # create null version of the data
nullA  = my.data$A - mean(my.data$A)
nullB  = my.data$B - mean(my.data$B, na.rm=T)
nullC  = my.data$C - mean(my.data$C, na.rm=T)
set.seed(1)
F.vect = vector(length=10000)
for(i in 1:10000){
  A = sample(na.omit(nullA), 15, replace=T)
  B = sample(na.omit(nullB), 13, replace=T)
  C = sample(na.omit(nullC), 13, replace=T)
  boot.dat  = stack(list(A=A, B=B, C=C))
  boot.mod  = lm(values~ind, boot.dat)
  F.vect[i] = anova(boot.mod)[1,4]
}
1-mean(F.stat>F.vect)
# [1] 0.0485

在某些方面，自举是进行参数分析（例如均值）的最终简化假设方法，但它确实假设您的数据很好地代表了总体，这意味着您有一个合理的样本量。由于您的很小，因此可能不太值得信赖。可能对非正态性和异质性的最终保护是使用非参数检验。ANOVA 的基本非参数版本是Kruskal-Wallis 检验： $n$

kruskal.test(values~ind, stacked.data, na.action=na.omit)
#  Kruskal-Wallis rank sum test
# 
# data:  values by ind
# Kruskal-Wallis chi-squared = 5.7705, df = 2, p-value = 0.05584

尽管 Kruskal-Wallis 检验绝对是防止 I 类错误的最佳方法，但它只能与单个分类变量一起使用（即，没有连续预测变量或因子设计），并且在所有讨论的策略中它的功效最小。另一种非参数方法是使用序数逻辑回归。这对很多人来说似乎很奇怪，但您只需要假设您的响应数据包含合法的序数信息，他们肯定会这样做，否则上述所有其他策略也是无效的：

library(rms)
olr.mod = orm(values~ind, stacked.data)
olr.mod
# Model Likelihood          Discrimination          Rank Discrim.    
# Ratio Test                 Indexes                Indexes       
# Obs            41    LR chi2      6.63    R2                  0.149    rho     0.365
# Unique Y       41    d.f.            2    g                   0.829
# Median Y    0.432    Pr(> chi2) 0.0363    gr                  2.292
# max |deriv| 2e-04    Score chi2   6.48    |Pr(Y>=median)-0.5| 0.179
#                      Pr(> chi2) 0.0391

从输出中可能不清楚，但模型作为一个整体的测试，在这种情况下是你的组的测试，是chi2under Discrimination Indexes。列出了两个版本，似然比测试和分数测试。似然比检验通常被认为是最好的。它产生一个值。 $p$ 0.0363

其它你可能感兴趣的问题

上一篇如何避免在 R 图中重叠标签？下一篇对于样本量为 1 的总体均值，我们能说什么？