This should have: Variables in Model IN P MSE R_2 R_2_adj CP AIC BIC Intercept # SOURCE: # http://ryouready.wordpress.com/2009/02/06/r-calculating-all-possible-linear-regression-models-for-a-given-set-of-predictors/ # identify your data dat <- dat # Identify your y y <- "medv" regressors <- names(dat)[names(dat) != dependent] # put transformation here y <- "log(medv)" #Now we want to construct a formula that contains the first and third regressor. vec_T <- rep(T, length(regressors)) vec <- vec_T i <- 2 vec[[i]] <- !vec[[i]] vec # #testing that paste works # paste(regressors[vec]) # # … and add the left side of the equation. The 1 in the formula models the intercept , 0 would be a model without intercept. # paste(c( paste(y, "~ 1"), regressors[vec]), collapse=" + ") # Now let’s make a formula out of it. as.formula(paste(c( paste(y, "~ 1"), regressors[vec]), collapse=" + ")) # So we can construct a formula from each row of a TRUE /FALSE matrix which determines # if a regressor is used or not. Now we need a TRUE / FALSE matrix of all the possible regressor combinations. # The expand.grid() function produces one (see ?expand.grid). # Total number of combinations are: n <- length(regressors) sum(sapply(seq(n), function(i) choose(n, i))) regMat <- expand.grid(c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE)) regMat <- regMat[1:(nrow(regMat)/2 -1), -(n+1)] regMat <- regMat[order(rowSums(regMat)),] # > regMat # Var1 Var2 Var3 Var4 # 1 TRUE TRUE TRUE TRUE # 2 FALSE TRUE TRUE TRUE # 3 TRUE FALSE TRUE TRUE # 4 FALSE FALSE TRUE TRUE # 5 TRUE TRUE FALSE TRUE # let's name the columns dimnames(regMat) <- list(1:nrow(regMat), regressors) head(regMat) # Now we can apply the above way of formula construction to each row of the matrix so # we get a list with all the possible models. allModelsList <- apply(regMat, 1, function(vec) as.formula(paste(c( paste(y, "~ 1"), regressors[vec]), collapse=" + ")) ) # > allModelsList # [[1]] # x1 ~ 1 + y1 + y2 + y3 + y4 # # [[2]] # x1 ~ 1 + y2 + y3 + y4 # # [[3]] # x1 ~ 1 + y1 + y3 + y4 #----------------------------------- # The last step is to use each list element for the calculation. allModelsResults <- lapply(allModelsList, function(x) lm(x, data=dat)) # So basically, here our computation work is done, but as in most cases a lot of # work follows to prepare the data in a nice way. # So now let’s get all the important information into one dataframe. # Let’s say we want a data frame like the following. # +-------+-----------------------------------------------------+ # | model | no. of | coefficients | se coef. | t-Val | etc. | # | | regressors | x1 | x2 ... | | | | # | | | | | | | # So we need to extract all the following information (coefficients, SE etc.) and cast them into one data frame. # x <- allModelsResults[[1]] # coef(x) # coef(summary(x))[, "Std. Error"] # ### ... and so on # This used to be one of the nasty tasks in R. Here Hadley Wickhams plyr package really helps a lot. ldply takes a list, applies a function and casts the results into ONE data frame (see ?ldply). As function return value it expects a data frame or a vector. The advantage to return data frames is that the ldply() function uses rbind.fill for combining the results when they are data frames. rbind.fill() allows a different number of columns in each data frame. Here this is the case as a different number of regressors are used each time. So we have to make sure that the function returns a data frame. Thus we use as.data.frame paying attention to the orientation of the data frame, using t() in case it is outputted as one column. library(plyr) dfCoefNum <- ldply(allModelsResults, function(x) as.data.frame( t(coef(x)))) dfStdErrors <- ldply(allModelsResults, function(x) as.data.frame( t(coef(summary(x))[, "Std. Error"]))) dftValues <- ldply(allModelsResults, function(x) as.data.frame( t(coef(summary(x))[, "t value"]))) dfpValues <- ldply(allModelsResults, function(x) as.data.frame( t(coef(summary(x))[, "Pr(>|t|)"]))) # rename DFs so we know what the column contains names(dfStdErrors) <- paste("se", names(dfStdErrors), sep=".") names(dftValues) <- paste("t", names(dftValues), sep=".") names(dfpValues) <- paste("p", names(dfpValues), sep=".") # p-value for overall model fit calcPval <- function(x){ fstat <- summary(x)$fstatistic pVal <- pf(fstat[1], fstat[2], fstat[3], lower.tail = FALSE) return(pVal) } # Before creating ONE data frame with all important entries, # we need to compute some more indices NoOfCoef <- unlist(apply(regMat, 1, sum)) R2 <- unlist(lapply(allModelsResults, function(x) summary(x)$r.squared)) adjR2 <- unlist(lapply(allModelsResults, function(x) summary(x)$adj.r.squared)) RMSE <- unlist(lapply(allModelsResults, function(x) summary(x)$sigma)) fstats <- unlist(lapply(allModelsResults, calcPval)) # now we can combine all the data into one data frame results <- data.frame( model = as.character(allModelsList), NoOfCoef = NoOfCoef, dfCoefNum, dfStdErrors, dftValues, dfpValues, R2 = R2, adjR2 = adjR2, RMSE = RMSE, pF = fstats ) # round the results results[,-c(1,2)] <- round(results[,-c(1,2)], 3) results$model <- as.character(results$model) rez <- results rez[, 63:66] <- round(rez[, 63:66], 3) names(rez)[names(rez) == "NoOfCoef"] <- "Vars" names(rez)[names(rez) == "RMSE"] <- "MSE" # Which columns will be displayed displayCols <- c(1:4, 63:66) varCols <- 5:17 varsIn <- c(NULL) #------ for (i in 1:5){ nexts <- findInRez(varsIn) nexts[, displayCols] (topRow <- nexts[1, varCols]) (varsIn <- names(topRow)[!is.na(topRow)]) } #--- repeat findInRez <- function(varsIn) { range <- rez$Vars==(length(varsIn) + 1) if (!is.null(varsIn)) range <- range & apply(!is.na(rez[,varsIn, drop=FALSE]), 1, all) current <- rez[range,] current[order(current$R2, decreasing=TRUE), ] }