I want to perform model comparison according to several criteria using R. My dataframe's name is df


   Y  X1  X2   X3   X4
1 18 307 130 3504 12.0
2 15 350 165 3693 11.5
3 18 318 150 3436 11.0
4 16 304 150 3433 12.0
5 17 302 140 3449 10.5
6 15 429 198 4341 10.0

I want to compare all possible combinations with Y as the dependent variable. I'll give you what I've done so far to get the R.squared for all possible models. I would like your comments for better coding, especially on getting all those model formulas

# First get all possible combinations of the 4 inependent variables
comb_list <- lapply(1:4,function(i) combn(4,i))

# now create the formulas
comb_list_forms <- lapply(unlist(
  sapply(1:length(comb_list), function(i)
  sapply(1:dim(comb_list[[i]])[2], function(x) {nam <- names(df[1+comb_list[[i]][,x]])
  formul <- "Y~"
 sapply(1:length(nam), function(y) formul <<- paste(formul,nam[y],sep="+"))
 formul <- sub("~+", "~", formul,fixed = T)

# finally get r.squared
models.r.sq <- sapply(comb_list_forms, function(i) summary(lm(i))$r.squared)

Thank you


I think package leaps will help you. It is designed specifically to select a "best" subset of variables for linear regression.

As for the formula creation, you can use foreach package like this:

 #Create the list where each list element contains one combination
  xcomb <- foreach(i=1:4, .combine=c) %do% { 
       combn(paste("X", 1:4, sep=""), i, simplify=FALSE)
 #Convert list elements to formula   
 formlist <- lapply(xcomb, function(l) formula(paste("Y", paste(l, collapse="+"), sep="~")))

Ok, here is another way (which is certainly not as elegant or concise than the other ones, but it has the merit of not requiring Emacs to check parenthesis matches :-)

Let's say we have a vector of predictors of interest like this

> Xs <- paste("X", 1:4, sep="")

Then, we can just use

> allXs <- lapply(seq(along=Xs), 
                  function(x) combn(Xs, x, simplify=FALSE, 
                                    FUN=paste, collapse="+"))


> unlist(allXs)

gives you all 15 combinations of the X's.

Another option is to just change the right-hand side of a formula, say

 fm <- as.formula(paste("Y ~ ", paste(Xs, collapse= "+")))

so as to reflect the different combinations that you enumerate in your comb_list object. This can be done using the Formula package:

> fm <- as.formula(paste("Y ~ ", paste(Xs, collapse= "|")))
> fm2 <- Formula(fm)
> foo <- function(x) formula(fm2, rhs=x, collapse=TRUE)
> foo(1:2)
Y ~ X1 + X2
<environment: 0x102593040>
> foo(c(1,3,4))
Y ~ X1 + X3 + X4
<environment: 0x1021d02a0>
  • $\begingroup$ +1. Did not know about this Formula feature. Thanks $\endgroup$ – mpiktas Dec 16 '10 at 15:16

This is another way to get the formulas

comb_list_forms <- unlist(lapply(1:length(comb_list), function(i)
  lapply(1:dim(comb_list[[i]])[2], function(x) 
    as.formula(paste("Y ~ ", paste(names(df[1+comb_list[[i]][,x]]), collapse= "+"))))))

You could consider using step on the full model. This uses AIC rather than $R^2$ as the criterion for selection as this is considered to be superior.

In any case, the functions add1 and drop1 that step uses may be of use.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.