R-code question: model selection based on individual significance in regression? I'm trying to generate an R function that keeps relevant variables based on their absolute t-value (or p, whichever is easier in code).
Basically what I want is to run one regression (1), retain all variables that are significant (based on the t-value, or p). Then run another regression (2), retain all variables that are significant (based on t-value, or p). Then take all retained variables from the regressions (1 & 2) above and combine them in a final regression (3).
Currently I am struggling with the step of retaining variables based on their t-value, I can extract the t-values of the variables but don't know how to code that they should be "retained" and kept for the later regression.
What I have currently (example edited for better understanding):
summaryI1 <- summary(lm(y~x1+x2+x3+x4)) #x1-x4 are the vars included in regression 1
(coefI1 <- coef(summaryI1))
tI1 <- coefI1[, "t value"] #extracts the t-values for all the x's

cutoff <- which(tI2>cut, arr.in=TRUE) #checks that t-stats of vars are above the cutoff level

Suppose x1 and x2 are significant and I would like to keep them "aside" for the next regression. How would I code that r should create a new matrix of x1 and x2 based on x1 and x2 having t-values higher than my cutoff?
I'm not that experienced in r programming yet but learning.
Can anyone help?
 A: Maybe some simulations will help you understand why these types of approaches can go bad.
First we will create some data, in this case x1 and x2 are correlated with each other (positive correlation) and y is a function of both of them (but more specifically their difference).  First we look at the regressions of each of the x's individually, then the regression that combines them:
set.seed(1)
x1 <- rnorm(25)
x2 <- rnorm(25, x1)
y <- x1-x2 + rnorm(25)
pairs( cbind(y,x1,x2) )
cor( cbind(y,x1,x2) )
summary(lm(y~x1))
summary(lm(y~x2))
summary(lm(y~x1+x2))

This shows that if you use some type of pre-screening on individual (or groups) of predictor variables is done, and only the "significant" predictors are then used in further analyses, then very important variables can be screened out.
Now lets go the other way and simulate some data where there is no relationship between the response and any of the predictors, but we first screen the potential predictors and put just the "likely" predictors into a new model:
set.seed(1)
y <- rnorm(100)
x <- matrix( rnorm(50*100), ncol=50 )
mydat <- data.frame( y,x )

fit1 <- lm( y~., data=mydat )
fit0 <- lm( y~1, data=mydat )
(out <- add1( fit0, fit1, test='F' ))

form <- as.formula( paste('y ~', 
    paste0('X', which(out$`Pr(>F)` <= 0.1)-1, collapse=' + ') ) )

fit2 <- lm( form, data=mydat )
summary(fit2)
summary(fit1)

The overall f-test on fit1 shows that there is nothing going on, which is the correct answer since y was simulated without any influence from the x's.  But the f-test on fit2 is significant (the wrong answer) because it does not take into account the screening done.
If simple simulations can show cases where stepwise type methods can leave out important variables and include useless ones, just think of the risks when using more complicated real world data.
A: So maybe it's not a good idea, but here's a start: 
    #first, some simulated data
    y = seq(from = 1, to = 10, by = 0.5)
    x1 = y + rnorm(length(y)) #not a significant variable
    x2 = y * runif(length(y), min = -5, max = 5) #a significant variable
    x3 = y + runif(length(y), min = -5, max = 5) #another insignificant variable

    model = lm(y ~ x1 + x2 + x3)
    summary(model)

    p_value = coef(summary(model))[,4]

    #create a new data frame to store y and only the significant regressors
    alpha = 0.05
    new_x = data.frame(model$model)
    for(i in 2:length(coef(model))){
       if(p_value[i] > alpha){
         new_x[,i] = NA
       }
    }

    #function to remove the NA columns
    remove.column.NA = function(data){
       i = 2
       while(i <= ncol(data)){
          if(is.na(data[1,i])){
             data = data[,-i]
             i = i-1
          } 
          else{
             i = i+1
          }
       }
       return(data)
    }

    new_x = remove.column.NA(new_x)

Now you have a dataframe (new_x) containing only the columns of significant regressors. I was trying to turn this into a complete function, but the names of the variables are lost in the lm() function.
