使用 p 值逐步回归以删除 p 值不显着的变量

Posted 2023-03-12

【中文标题】使用 p 值逐步回归以删除 p 值不显着的变量

【英文标题】：Stepwise regression using p-values to drop variables with nonsignificant p-values

【发布时间】：2011-04-11 16:46:15

【问题描述】：

我想使用 p 值 作为选择标准来执行 逐步线性回归，例如：在每一步删除具有最高即最不显着 p 的变量-values，当所有值都由某个阈值 alpha 定义时停止。

我完全知道我应该使用 AIC（例如命令 step 或 stepAIC）或其他一些标准，但我的老板不懂统计并坚持关于使用 p 值。

如果有必要，我可以编写自己的例程，但我想知道是否有已经实现的版本。

【问题讨论】：

希望你的老板不看***。 ;-)

你可以考虑在这里问这个问题：stats.stackexchange.com

教你老板好的数据可能比让 R 做坏数据更容易。

随机选择三个变量 - 你可能会做的一样好逐步回归。

你的老板是否也告诉他的医生开什么药和他的机械师如何修理他的车？

【参考方案1】：

向你的老板展示以下内容：

set.seed(100)
x1 <- runif(100,0,1)
x2 <- as.factor(sample(letters[1:3],100,replace=T))

y <- x1+x1*(x2=="a")+2*(x2=="b")+rnorm(100)
summary(lm(y~x1*x2))

这给出了：

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.1525     0.3066  -0.498  0.61995    
x1            1.8693     0.6045   3.092  0.00261 ** 
x2b           2.5149     0.4334   5.802 8.77e-08 ***
x2c           0.3089     0.4475   0.690  0.49180    
x1:x2b       -1.1239     0.8022  -1.401  0.16451    
x1:x2c       -1.0497     0.7873  -1.333  0.18566    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

现在，根据 p 值，您会排除哪一个？ x2 同时是最显着和最不显着的。

---

编辑：澄清：这个例子不是最好的，如 cmets 所示。 Stata 和 SPSS 中的程序是 AFAIK，也不是基于系数 T 检验的 p 值，而是基于去除变量之一后的 F 检验。

我有一个功能可以做到这一点。这是对“p 值”的选择，而不是对系数或方差分析结果的 T 检验。好吧，如果它对您有用，请随意使用。

#####################################
# Automated model selection
# Author      : Joris Meys
# version     : 0.2
# date        : 12/01/09
#####################################
#CHANGE LOG
# 0.2   : check for empty scopevar vector
#####################################

# Function has.interaction checks whether x is part of a term in terms
# terms is a vector with names of terms from a model
has.interaction <- function(x,terms)
    out <- sapply(terms,function(i)
        sum(1-(strsplit(x,":")[[1]] %in% strsplit(i,":")[[1]]))==0
    )
    return(sum(out)>0)


# Function Model.select
# model is the lm object of the full model
# keep is a list of model terms to keep in the model at all times
# sig gives the significance for removal of a variable. Can be 0.1 too (see SPSS)
# verbose=T gives the F-tests, dropped var and resulting model after 
model.select <- function(model,keep,sig=0.05,verbose=F)
      counter=1
      # check input
      if(!is(model,"lm")) stop(paste(deparse(substitute(model)),"is not an lm object\n"))
      # calculate scope for drop1 function
      terms <- attr(model$terms,"term.labels")
      if(missing(keep)) # set scopevars to all terms
          scopevars <- terms
       else            # select the scopevars if keep is used
          index <- match(keep,terms)
          # check if all is specified correctly
          if(sum(is.na(index))>0)
              novar <- keep[is.na(index)]
              warning(paste(
                  c(novar,"cannot be found in the model",
                  "\nThese terms are ignored in the model selection."),
                  collapse=" "))
              index <- as.vector(na.omit(index))
          
          scopevars <- terms[-index]
      

      # Backward model selection : 

      while(T)
          # extract the test statistics from drop.
          test <- drop1(model, scope=scopevars,test="F")

          if(verbose)
              cat("-------------STEP ",counter,"-------------\n",
              "The drop statistics : \n")
              print(test)
          

          pval <- test[,dim(test)[2]]

          names(pval) <- rownames(test)
          pval <- sort(pval,decreasing=T)

          if(sum(is.na(pval))>0) stop(paste("Model",
              deparse(substitute(model)),"is invalid. Check if all coefficients are estimated."))

          # check if all significant
          if(pval[1]<sig) break # stops the loop if all remaining vars are sign.

          # select var to drop
          i=1
          while(T)
              dropvar <- names(pval)[i]
              check.terms <- terms[-match(dropvar,terms)]
              x <- has.interaction(dropvar,check.terms)
              if(x)i=i+1;next else break              
           # end while(T) drop var

          if(pval[i]<sig) break # stops the loop if var to remove is significant

          if(verbose)
             cat("\n--------\nTerm dropped in step",counter,":",dropvar,"\n--------\n\n")              
          

          #update terms, scopevars and model
          scopevars <- scopevars[-match(dropvar,scopevars)]
          terms <- terms[-match(dropvar,terms)]

          formul <- as.formula(paste(".~.-",dropvar))
          model <- update(model,formul)

          if(length(scopevars)==0) 
              warning("All variables are thrown out of the model.\n",
              "No model could be specified.")
              return()
          
          counter=counter+1
       # end while(T) main loop
      return(model)

【讨论】：

很好的例子——几个月前我就可以用这个了！

我认为很明显该决定并非完全基于 p 值。您首先要删除最不重要的最高阶交互作用，然后是潜在的二次项，然后再考虑要删除的任何解释变量。

当然不是。但是删除交互项，您会看到情况保持不变。所以你实际上应该去anova，但是你遇到了哪种类型的问题。这是一罐我不会打开的蠕虫。

没有机会在统计数据足够的水平上与老板争论 ;-) 尽管如此，我会尝试并在此之前“手动”进行逐步选择。奇怪的是，没有一个 R 函数，即使使用 p 值的方法来自计算 AIC 的计算成本太高......

我知道，这远非最好的例子。这只是为了提出一个观点。在任何情况下，“基于 p 值”的方法都是基于 F 检验，而不是基于系数的 t 检验。所以我添加的功能应该给 OP 他想要的。

【参考方案2】：

为什么不尝试使用step() 函数指定您的测试方法？

例如，对于向后消除，您只需键入一个命令：

step(FullModel, direction = "backward", test = "F")

对于逐步选择，只需：

step(FullModel, direction = "both", test = "F")

这可以显示 AIC 值以及 F 和 P 值。

【讨论】：

请注意，此过程将添加/删除具有最重要测试值的变量。但是继续或停止的选择仍然基于 AIC。我认为一旦没有重要的变量（就像 OP 想要的那样），就不可能停止该过程。

【参考方案3】：

这是一个例子。从最复杂的模型开始：这包括所有三个解释变量之间的交互作用。

model1 <-lm (ozone~temp*wind*rad)
summary(model1)

Coefficients:
Estimate Std.Error t value Pr(>t)
(Intercept) 5.683e+02 2.073e+02 2.741 0.00725 **
temp          -1.076e+01 4.303e+00 -2.501 0.01401 *
wind          -3.237e+01 1.173e+01 -2.760 0.00687 **
rad           -3.117e-01 5.585e-01 -0.558 0.57799
temp:wind      2.377e-01 1.367e-01 1.739 0.08519   
temp:rad       8.402e-03 7.512e-03 1.119 0.26602
wind:rad       2.054e-02 4.892e-02 0.420 0.47552
temp:wind:rad -4.324e-04 6.595e-04 -0.656 0.51358

三向交互显然不显着。这是您删除它的方式，开始模型简化过程：

model2 <- update(model1,~. - temp:wind:rad)
summary(model2)

根据结果，您可以继续简化模型：

model3 <- update(model2,~. - temp:rad)
summary(model3)
...

或者你可以使用自动模型简化功能step，看看 效果如何：

model_step <- step(model1)

【讨论】：

感谢您的回答。我刚刚意识到我没有足够清楚地说明我的问题：我正在寻找一个使用 p 值作为标准自动执行此操作的命令或包。命令“step”使用 AIC。当您建议或编写一些自动执行此操作的例程时，我也可以“手动”执行此操作，但已经实现的版本会非常方便。

【参考方案4】：

包rms: Regression Modeling Strategies 具有fastbw()，可以满足您的需求。甚至还有一个参数可以从 AIC 翻转到基于 p 值的消除。

【参考方案5】：

如果您只是想获得最好的预测模型，那么也许这并没有太大的关系，但对于其他任何事情，不要为这种模型选择而烦恼。错了。

使用收缩方法，例如岭回归（例如，在 MASS 包中的 lm.ridge() 中）、套索或弹性网（岭和套索约束的组合）。其中，只有 lasso 和弹性网络会做某种形式的模型选择，即将一些协变量的系数强制为零。

请参阅 CRAN 上 Machine Learning 任务视图的正则化和收缩部分。

【参考方案6】：

正如 Gavin Simpson 所述，rms 包中的函数 fastbw 可用于使用 p 值选择变量。 Bellow 是一个使用 George Dontas 给出的例子的例子。使用选项 rule='p' 选择 p 值标准。

require(rms)
model1 <- ols(Ozone ~ Temp * Wind * Solar.R, data=airquality)
model2 <- fastbw(fit=model1, rule="p", sls=0.05)
model2

【讨论】：

在文本中，您指的是 rms 包。还有一个 rsm 包，所以这可能会引起混淆。

感谢 rvl 的澄清。我修好了。