使用 lmrob 的稳健回归 p 值答案

【问题标题】：Robust regression p values with lmrob使用 lmrob 的稳健回归 p 值
【发布时间】：2022-06-10 16:39:36
【问题描述】：

我正在尝试对一个数据集进行稳健的多元回归，其中一些异常值不允许我通过通常的线性模型看到底层模式。我正在使用包robustbase 中的函数lmrob，我对我发现的重要关系的数量感到惊讶。我决定尝试使用随机数据的方法，代码如下：

library(robustbase)
set.seed(4)
ax<-data.frame(a1=rnorm(20,3),
               a2=rnorm(20,5),
               a3=rnorm(20,4),
               a4=rnorm(20,6),
               a5=rnorm(20,2))
 
axm<-lmrob(a1~a2*a3*a4*a5,data=ax)
summary(axm)

还有输出：

        Call:
lmrob(formula = a1 ~ a2 * a3 * a4 * a5, data = ax)
 \--> method = "MM"
Residuals:
         1          2          3          4          5          6          7          8          9         10         11         12         13 
-34.740270  -0.049493  -0.044379   0.002770   0.219825   0.041285   0.156152  -0.072825   0.034824  -0.014757  -0.088263  -0.185045  -0.079679 
        14         15         16         17         18         19         20 
 -0.045121  -0.007576   0.008813   0.010451   0.015716   0.060781   0.040187 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1160.5907    94.0095  -12.35 0.000247 ***
a2            205.6910    15.8689   12.96 0.000204 ***
a3            327.9787    24.2161   13.54 0.000172 ***
a4            193.2384    15.7300   12.29 0.000252 ***
a5            734.2203    49.8960   14.71 0.000124 ***
a2:a3         -57.6229     4.0533  -14.22 0.000142 ***
a2:a4         -33.5644     2.6130  -12.85 0.000212 ***
a3:a4         -54.1622     4.0438  -13.39 0.000180 ***
a2:a5        -138.8395     9.2697  -14.98 0.000116 ***
a3:a5        -198.4961    12.3168  -16.12 8.67e-05 ***
a4:a5        -123.0895     8.2792  -14.87 0.000119 ***
a2:a3:a4        9.3344     0.6659   14.02 0.000150 ***
a2:a3:a5       37.1371     2.2502   16.50 7.89e-05 ***
a2:a4:a5       23.0014     1.5152   15.18 0.000110 ***
a3:a4:a5       32.9766     2.0388   16.18 8.55e-05 ***
a2:a3:a4:a5    -6.0817     0.3660  -16.62 7.68e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Robust residual standard error: 0.4039 
Multiple R-squared:  0.9861,    Adjusted R-squared:  0.934 
Convergence in 5 IRWLS iterations

Robustness weights: 
 observation 1 is an outlier with |weight| = 0 ( < 0.005); 
 9 weights are ~= 1. The remaining 10 ones are
     2      3      5      7      8     11     12     13     14     19 
0.9986 0.9989 0.9732 0.9864 0.9970 0.9957 0.9810 0.9965 0.9989 0.9979 
Algorithmic parameters: 
       tuning.chi                bb        tuning.psi        refine.tol           rel.tol         scale.tol         solve.tol       eps.outlier 
        1.548e+00         5.000e-01         4.685e+00         1.000e-07         1.000e-07         1.000e-10         1.000e-07         5.000e-03 
            eps.x warn.limit.reject warn.limit.meanrw 
        1.150e-09         5.000e-01         5.000e-01 
     nResample         max.it       best.r.s       k.fast.s          k.max    maxit.scale      trace.lev            mts     compute.rd 
           500             50              2              1            200            200              0           1000              0 
fast.s.large.n 
          2000 
                  psi           subsampling                   cov compute.outlier.stats 
           "bisquare"         "nonsingular"         ".vcov.avar1"                  "SM" 
seed : int(0)

据此，我了解到其他随机变量与第一个相关，并且对它有很高的预测能力，这是没有意义的。

这里发生了什么？我做回归错了吗？

编辑：我放置了一个 p 值极低的种子以实现可复制性。

【问题讨论】：

再运行几次你会看到你会得到非常不同的结果。这纯粹是过拟合（包含所有可能的交互对您没有帮助）。
@user2974951 是的，我注意到它发生了变化，通常会选择 2 或 3 个不同的选项。在使用我的数据的真实模型中，我只测试了 1 个变量与其余变量的交互，但结果看起来还是过拟合。这个函数对了解哪些变量有显着影响有用吗？
我应该提一下，您向我们展示的模型（以及模型结果）来自常规线性模型，而不是来自 lmrob。
您的代码以如此低的 p 值给出输出的概率非常低。
@user2974951 哦，我的错，代码错误，因为我后来用 lm() 尝试并复制了错误的，但输出实际上来自 lmrob()。使用 lm，p 值要高得多，大约为 0.7

标签： r regression robust

【解决方案1】：

我想我可能已经找到了对如此高的 p 值的解释：结果表明 MM 估计样本量较小（如我的示例中为 20）is prone to type 1 error。 robustbase 软件包的作者之一发布了an article，为这种情况提出了另一种估计，但恐怕就我的数据而言，它的效果并没有那么好。

【讨论】：