Regresión paso a paso usando valores p para descartar variables con valores de p no significativos

Quiero realizar una regresión lineal paso a paso utilizando p-values como criterio de selección, por ejemplo: en cada paso soltar las variables que tienen los p-values ​​más altos, es decir, los más insignificantes, deteniéndose cuando todos los valores son significativos definidos por algún umbral alfa .

Soy totalmente consciente de que debería usar el AIC (por ejemplo, paso de comando o paso AIC ) o algún otro criterio en su lugar, pero mi jefe no tiene conocimiento de las estadísticas e insiste en usar valores p.

Si es necesario, podría progtwigr mi propia rutina, pero me pregunto si ya hay una versión implementada de esto.

Muestra a tu jefe lo siguiente:

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

Lo que da :

  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 

Ahora, según los valores p, ¿excluiría cuál? x2 es más significativo y la mayoría no significativo al mismo tiempo.


Editar: Para aclarar: este ejemplo no es el mejor, como se indica en los comentarios. El procedimiento en Stata y SPSS es AFAIK y no se basa en los valores p de la prueba T en los coeficientes, sino en la prueba F después de la eliminación de una de las variables.

Tengo una función que hace exactamente eso. Esta es una selección en "el valor p", pero no de la prueba T en los coeficientes ni en los resultados de anova. Bueno, siéntase libre de usarlo si le parece útil.

 ##################################### # 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] 

¿Por qué no intentar usar la función step() especificando su método de prueba?

Por ejemplo, para la eliminación hacia atrás, escribe solo un comando:

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

y para la selección gradual, simplemente:

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

Esto puede mostrar tanto los valores AIC como los valores F y P.

Aquí hay un ejemplo. Comience con el modelo más complicado: esto incluye las interacciones entre las tres variables explicativas.

 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 

La interacción de tres vías claramente no es significativa. Así es como lo quitas, para comenzar el proceso de simplificación del modelo:

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

Dependiendo de los resultados, puede continuar simplificando su modelo:

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

Alternativamente, puede usar el step función de simplificación automática del modelo para ver qué tan bien funciona:

 model_step <- step(model1) 

Package rms: Regression Modeling Strategies tiene fastbw() que hace exactamente lo que necesita. Incluso hay un parámetro para cambiar de AIC a eliminación basada en p-value.

Si solo está tratando de obtener el mejor modelo predictivo, entonces tal vez no importe demasiado, pero para cualquier otra cosa, no se moleste con este tipo de selección de modelos. Está mal.

Utilice métodos de contracción como la regresión de cresta (en lm.ridge() en el paquete MASS, por ejemplo), o el lazo, o el elasticnet (una combinación de restricciones de cresta y lazo). De estos, solo el lazo y la red elástica harán alguna forma de selección del modelo, es decir, forzarán los coeficientes de algunas covariables a cero.

Consulte la sección Regularización y contracción de la vista de tareas de Aprendizaje automático en CRAN.