multiple - regresion lineal ejemplos

Regresión por pasos utilizando valores de p para eliminar variables con valores de p no significativos

Question

Quiero realizar una Regresión lineal por pasos utilizando los valores de p como criterio de selección, por ejemplo: en cada paso, se eliminan las variables que tienen los valores p 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, el paso de comando o el paso AIC ) o algún otro criterio, pero mi jefe no entiende las estadísticas e insiste en usar valores p.

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

Answer 1

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

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

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

y para la selección por pasos, simplemente:

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

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

Answer 2

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 se elimina, para comenzar el proceso de simplificación del modelo:

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

Dependiendo de los resultados, puedes continuar simplificando tu modelo:

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

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

model_step <- step(model1)

Answer 3

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, basado en los valores p, ¿excluirías cuál? x2 es más significativo y más no significativo al mismo tiempo.

Edición: para aclarar: este ejemplo no es el mejor, como se indica en los comentarios. El procedimiento en Stata y SPSS para AFAIK tampoco 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 o en los resultados del anova. Bueno, siéntete libre de usarlo si te 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]<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)
}

Answer 4

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

Answer 5

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 modelo. Está mal.

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

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