R语言 使用p值的逐步回归法删除具有非显著p值的变量

vm0i2vca  于 2023-01-10  发布在  其他
关注(0)|答案(7)|浏览(304)

我想使用p值作为选择标准执行逐步线性回归,例如:在每一步丢弃具有最高即最不显著p值的变量,当由某个阈值 * α * 定义的所有值都显著时停止。
我完全知道我应该使用AIC(例如,命令 * 步骤 * 或 * 步骤AIC *)或其他一些标准,但我的老板不懂统计学,坚持使用p值。
如果有必要,我可以编写自己的例程,但我想知道是否有一个已经实现的版本。

voase2hg

voase2hg1#

向老板展示以下内容:

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同时最显著和最不显著。
编辑:澄清:这个例子不是最好的,正如评论中所指出的。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)
}
js81xvg6

js81xvg62#

为什么不尝试使用step()函数指定测试方法呢?
例如,对于向后消除,只需键入一个命令:

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

对于逐步选择,只需:

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

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

ui7jx7zq

ui7jx7zq3#

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

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)
vfwfrxfs

vfwfrxfs4#

软件包rms: Regression Modeling Strategies中的fastbw()完全满足您的需求,甚至还有一个参数可以从AIC转换为基于p值的消除。

7qhs6swi

7qhs6swi5#

如果你只是想得到最好的预测模型,那么也许这并不太重要,但对于其他任何事情,不要为这种模型选择而烦恼,这是错误的。
使用收缩方法,如岭回归(例如MASS软件包中的lm.ridge())、套索或弹性网络(岭约束和套索约束的组合)。在这些方法中,只有套索和弹性网络会执行某种形式的模型选择,即强制某些协变量的系数为零。
请参见CRAN上Machine Learning任务视图的正则化和收缩部分。

mec1mxoz

mec1mxoz6#

如Gavin Simpson所述,rms程序包中的函数fastbw可用于选择使用p值的变量。Bellow是使用乔治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
xu3bshqb

xu3bshqb7#

OLSRR封装可能是有用。
您可以定义pent(进入模型的p值)和prem(移除的p值)
输出提供了您需要的所有指标,甚至更多。

相关问题