krokowej regresji za pomocą wartości p spadek zmiennymi nieistotna wartości p

Question

że chce wykonać krokowej regresji liniowej pomocą wartości p jako kryterium wyboru, na przykład: na każdym etapie spada zmienne, które mają najwyższy IE najbardziej nieistotne wartości p, zatrzymujące się, gdy wszystkie wartości są znaczące określone przez pewien próg alfa.

Jestem całkowicie świadomy, że powinienem użyć AIC (np polecenia krok lub stepAIC) lub innego kryterium w zamian, ale mój szef ma pojęcie o statystyce i nalegać na użyciu wartości p.

Jeśli to konieczne, mogę zaprogramować własną procedurę, ale zastanawiam się, czy istnieje już wdrożona wersja tego.

Answer 1

Pokaż szef następujące:

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

Co daje:

  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 

Teraz, na podstawie wartości p byś wykluczają który? x2 jest najbardziej znaczący i najbardziej nieistotny w tym samym czasie.

---

Edytuj: Aby wyjaśnić: To exaxmple nie jest najlepsze, jak wskazano w komentarzach. Procedura w Stata i SPSS to AFAIK również nie bazująca na wartościach p testu T na współczynnikach, ale na teście F po usunięciu jednej ze zmiennych.

Mam funkcję, która dokładnie to robi. Jest to wybór "wartości p", ale nie testu T na współczynnikach lub wyników anova. Możesz go używać, jeśli ci się przyda.

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

Oto przykład. Zacznij od najbardziej skomplikowanego modelu: obejmuje to interakcje między wszystkimi trzema zmiennymi objaśniającymi.

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 

Interakcja trójstronna najwyraźniej nie ma znaczenia. W ten sposób można usunąć go, aby rozpocząć proces modelu Uproszczenie:

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

zależności od wyników, można kontynuować uproszczenie modelu:

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

Alternatywnie można skorzystać z funkcji automatycznego modelu uproszczenia step, aby zobaczyć jak dobrze robi:

model_step <- step(model1) 

Answer 3

Jeśli jesteś po prostu staramy się uzyskać najlepszą modelu prognostycznego, to może to nie ma znaczenia zbyt wiele, ale na nic innego, nie przejmuj się tym rodzajem wyboru modelu. To jest złe.

Użyj metod skurczu, takich jak regresja grzbietu (na przykład w przypadku lm.ridge() w pakiecie MASS) lub lasso, lub elasticnet (połączenie ograniczeń grzbietów i lasso). Spośród nich tylko lasso i elastyczna sieć dokonają pewnej selekcji modelu, tj. Wymuszą zerowanie współczynników niektórych współzmiennych.

Zobacz sekcję Unieważniania i kurczenie się w widoku zadań na CRAN w widoku Machine Learning.

Answer 4

Pakiet rms: Regression Modeling Strategies ma fastbw(), który spełnia dokładnie to, czego potrzebujesz. Istnieje nawet parametr do przechodzenia z AIC do eliminacji opartej na wartości p.

Answer 5

Dlaczego nie spróbować użyć funkcji step(), określając swoją metodę testowania?

Na przykład, dla wstecznej eliminacji wpisać tylko polecenie:

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

oraz dla stopniowej selekcji, po prostu:

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

ten może wyświetlać zarówno wartości AIC oraz F i wartości P.