使用显著性检验分析对于一个具有多个不平等水平的因素和二元依赖变量的影响。

Question

在现实生活中，我试图根据众多故障代码来预测车辆的健康状况。

有些故障码（因素的级别）出现得非常频繁（1000+），而其他的只出现两三次。通常情况下，那些很少出现的故障码是健康状态的“完美预测器”（0或1）。我试图找到一种可靠的统计方法来确定哪些因素级别是好的预测变量（显著的）。这样，那些很少出现但是是好的预测变量的故障码就不会仅仅因为它们的罕见性而被舍弃。

数据创建

library(tidyverse)
n_small =   4
n_big   = 100

set.seed(567)
df_big_1 <- data.frame(class = rep("A", n_big),
                       health = rbinom(n = n_big, size = 1, prob = .4))
df_small_1 <- data.frame(class = rep("B", n_small),
                         health = rbinom(n = n_small, size = 1, prob = 1))
df_small_2 <- data.frame(class = rep("C", n_small), 
                         health = rbinom(n = n_small, size = 1, prob = 1))
df_big_2 <- data.frame(class = rep("D", n_big),
                       health = rbinom(n = n_big, size = 1, prob = .4))
df_big_3 <- data.frame(class = rep("E", n_big),
                       health = rbinom(n = n_big, size = 1, prob = .4))
df_data <- rbind(df_big_1 ,df_small_1, df_big_2, df_small_2, df_big_3)
df_data <- df_data %>% mutate(class = factor(class))

数据检查

df_data %>%
  group_by(class) %>%
  summarise(N_health = sum(health), Mean = mean(health))
# A tibble: 5 × 3
  class N_health  Mean
  <fct>    <int> <dbl>
1 A           36  0.36
2 B            4  1   
3 C            4  1   
4 D           40  0.4 
5 E           40  0.4 

（二元）逻辑回归

在这个简化数据集上运行二元逻辑回归时，我未能获取罕见但“完美”的预测变量：

regmod_01 <- glm(health ~ class, family = binomial, data = df_data)
summary(regmod_01

Call:
glm(formula = health ~ class, family = binomial, data = df_data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.0108  -1.0108  -0.9448   1.3537   1.4294  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)   
(Intercept)   -0.5754     0.2083  -2.762  0.00575 **
classB        17.1414  1199.7724   0.014  0.98860   
classC        17.1414  1199.7724   0.014  0.98860   
classD         0.1699     0.2917   0.583  0.56022   
classE         0.1699     0.2917   0.583  0.56022   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 415.22  on 307  degrees of freedom
Residual deviance: 399.89  on 303  degrees of freedom
AIC: 409.89

Number of Fisher Scoring iterations: 15

我还有其他的方法可以尝试，来区分好的和坏的预测变量并且包含那些不经常出现的吗？

- KoenV

Answer 1

这是统计学中所谓的Hauck-Donner效应的经典例子。

来自麦克马斯特大学，感谢@BenBolker：

Hauck-Donner效应发生在参数估计极端的情况下（例如完全或接近完全分离的情况下），此时二次近似非常糟糕：其标志是较大的参数估计值（例如|βˆ|>10）和非常大的置信区间（导致小的Z统计量和大的p值）。

类别B和类别C的预测变量不频繁并不是它们的p值比您预期的要高得多的原因。如果您将n_small改为90而不是5，再次运行相同的代码，类别B和类别C仍然具有0.98的高p值。

n_small =   90
...
regmod_01 <- glm(health ~ class, family = binomial, data = df_data)
summary(regmod_01)


Coefficients:
             Estimate Std. Error z value Pr(>|z|)   
(Intercept)   -0.5754     0.2083  -2.762  0.00575 **
classB        20.1414  1133.5725   0.018  0.98582   
classC        20.1414  1133.5725   0.018  0.98582   
classD         0.1699     0.2917   0.583  0.56022   
classE         0.1699     0.2917   0.583  0.56022   

然而，如果您将B和C的概率从完美的1略微降低到0.97，并保持n_small为90，则p值会从几乎为1.00急剧变化为极显著（p < 0.001），因为Wald检验重新开始正常工作。

n_small =   90
...
df_small_1 <- data.frame(class = rep("B", n_small),
                         health = rbinom(n = n_small, size = 1, prob = 0.97))
df_small_2 <- data.frame(class = rep("C", n_small), 
                         health = rbinom(n = n_small, size = 1, prob = 0.97))
...
regmod_01 <- glm(health ~ class, family = binomial, data = df_data)
summary(regmod_01)

              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -5.754e-01  2.083e-01  -2.762  0.00575 ** 
classB       5.064e+00  1.027e+00   4.931 8.18e-07 ***
classC       3.943e+00  6.231e-01   6.328 2.49e-10 ***
classD       3.779e-15  2.946e-01   0.000  1.00000    
classE       4.315e-02  2.938e-01   0.147  0.88323   

标准glm输出中得到的p值是通过Wald检验获得的，这在极端分离（完美预测器）的情况下被视为极其不可靠。对于二项回归，似然比检验（LRT）是更合适的替代方法。

目前，glm没有内置的Hauck-Donner效应检测支持。但是，如果我们使用VGAM :: vglm重复此测试，我们会在底部看到一个适当的警告信息。

library(VGAM)

regmod_01 <- VGAM::vglm(health ~ class, family = binomialff, data = df_data)
summary(regmod_01)

(output)

Warning: Hauck-Donner effect detected in the following estimate(s):
'classB', 'classC'

回答你的问题，你可以选择具有显著p值的因素，以及任何系数/估计值超过10或低于0.1的因素，这些因素可能在Wald测试中不被视为显著，但实际上非常显著。

请记住，每个因素水平的p值是参考虚拟对比编码方案中的截距级别，因此根据您选择的参考点，显著性和几率比可能会大相径庭。最好使用似然比检验来确定回归模型中整个变量的显着性，而不是单个水平。

- dcsuka

---

非常有趣。然而，当存在许多实际上不重要的小因子水平时，我想知道这种方法如何有用。例如，想象一个因子由100个水平组成，每个水平都有2个随机生成的观测值。你会发现其中一半似乎是完美的预测因子，尽管我们知道这并不是真实情况。 - Iyar Lin

@IyarLin，一般来说，这样的数据并不适合进行统计推断，而且在这种情况下，无论使用哪种方法，系数的意义都不是很明显。也许应该添加自定义的样本大小截止值或标准。但再次强调，Wald的p值和系数都是相对于参考水平的，并且可以用更稳健的测试方法（如似然估计）进行替换。我只是提供了一个经验法则，它在许多情况下可以被常识所取代。 - dcsuka

Answer 2

我认为你对p值的理解可能有些误解。该指标告诉您结果的可靠性。您的变量可能未能解释任何关于结果变量的内容，但仍然显著。假设该变量基本上没有影响，但其p值为0.01。这告诉您，如果有人使用与您的数据相同的样本重新进行分析，他们中有99个人会得到与您相同的结果...即变量对结果没有任何影响。

我不确定您对该分析方法或常见误解的p值有多忠诚。我认为您应该看看生存集成（包括时间因素）或基于决策的分类器。

例如，如果我将您的结果变量设置为因子，并使用caret和randomForest，则告诉我B和C类在您的数据集中最具影响力（importance）。varImp（）提供此信息。

例如（使用以下结果），B类的平均准确度降低为16.306。如果排除此类，则准确度将下降超过16％。E和D类具有负面结果。这表明，如果您不包括这些类别，则准确性将提高。

为使此方法有效，您需要安装randomForest软件包。不过，您不需要调用该库。

看看吧：

library(caret) 
    
df_data$health <- factor(df_data$health)

set.seed(3253)
tr <- createDataPartition(df_data$health, p = .7, list = F)

cfit <- train(health~., data = df_data[tr, ], importance = T, keep.forest = T)

varImp(cfit, scale = F)  # mean decrease in accuracy
# rf variable importance
# 
#        Importance
# classB     16.306
# classC     11.287
# classE     -0.967
# classD     -2.036 

varImp(cfit, type = 2)  # mean increase in purity
# rf variable importance
# 
#        Overall
# classB  100.00
# classC   57.21
# classD    0.60
# classE    0.00 

有很多基于决策树的分析方法。随机森林的一个好处是几乎没有任何正式的统计假设和过度拟合的能力。同样的，许多基于决策树的方法也可以这样说，但并非全部。

- Kat

---

完美/接近完美的预测器应该具有非常低甚至为零的p值。尽管存在完美分离，但高p值是由于Hauck-Donner效应造成的，默认的Wald检验的统计误差。 - dcsuka

@Kat，我在哪里可以了解更多关于p值与变量重要性的知识？正是由于对p值的困惑，我正在转向更为稳健的方法，如分位数回归。 - alejandro_hagan

p值与重要性完全无关。如果你找到了比较它们的东西，我会非常惊讶。有关p值的一些信息，请参见此处：https://towardsdatascience.com/how-to-understand-p-value-in-layman-terms-80a5cc206ec2。有关变量重要性和类似指标的信息，请查看此处：https://topepo.github.io/caret/variable-importance.html，参见第15.1节。 - Kat

Answer 3

一般来说，当考虑给定因素水平与预测变量之间的关系时，可以测量两个经常相互矛盾的指标：

1. 水平“预测准确性” 2. 频率水平

类别水平越小，它似乎仅凭偶然就具有更高的预测价值。例如，考虑100个类别水平，每个水平只有2个随机生成的观察结果。你会发现大约一半的水平“完美地”预测了结果（有时是纯0，有时是纯1），即使没有真正的关系。这就是为什么在回归中也会看到高的p值。

决定哪些类别水平具有较高的预测价值以及推广到哪些变量具有较高的预测价值的一种非常自然的方法是使用条件推断树。

这种类型的树使用置换检验，不仅考虑预测准确性，还考虑每个类别水平的大小。您可以在链接的vignette中了解更多有关此类树在这方面与经典CART不同的信息。

在您的数据集中，我们得到：

library(partykit)
treemod <- ctree(health ~ class, data = df_data)
treemod

Fitted party:
[1] root
|   [2] class in A, D, E: 0.387 (n = 300, err = 71.1)
|   [3] class in B, C: 1.000 (n = 8, err = 0.0)

Number of inner nodes:    1
Number of terminal nodes: 2

你可以看到，树确实找到了2个“不常见”的层级，就像你预期的那样。从树结构中还可以看出，A、D和E层级之间的差异不足以被视为单独的层级。

为了了解这种方法的有用性，让我们来看一个具有100个类别层级和每个层级随机生成的2个观测值的例子：

df_data2 <- data.frame(class = factor(rep(1:100, each = 2)), 
                      health = rbinom(n = 200, size = 1, prob = 0.5))

接下来我们再次运行决策树，可以看到它没有进行任何分裂（也就是没有重要因素），这正是我们所期望的：

treemod2 <- ctree(health ~ class, data = df_data2)
treemod2

Fitted party:
[1] root: 0.530 (n = 200, err = 49.8) 

Number of inner nodes:    0
Number of terminal nodes: 1

- Iyar Lin