... and 5 cents more.

I was solving similar problem some time ago.
Here is my script.
I had a data frame, containing a response and several other variables,
which were assumed predictors.
I was trying to choose the best linear approximation.
This approach now seems to me useless, please, don't blame me for that.
However, the script might be useful to you.

<code>
library(forward)

# dfr is a data.frame, that contains everything.
# The response variable is named med5x
# The following lines construct linear models for all possibe formulas
# of the form 
# med5x~T+a+height
# med5x~a+height+RH
# T, a, RH, etc are the names of possible predictors

inputs<-names(dfr)[c(10:30,1)]  # dfr was a very large data frame,
containing lot of variables.
# here we have chosen only a subset of them.

for(nc in 11:length(inputs)){ # the linear models were assumed to have at
least 11 terms
# now we are generating character vectors containing formulas.

  formulas<-paste("med5x",sep="~",
                 
fwd.combn(inputs,nc,fun=function(x){paste(x,collapse="+")}))

# and then, are trying to fit every

  for(f in formulas){
    lms<-lm(eval(parse(text=f)),data=dfr)

   
cat(file="linear_models.txt",f,sum(residuals(lms)^2),"\n",sep="\t",append=TRUE)
  }
}
</code>

Hmm, looking back, I see that this is rather inefficient script.
For example, the inner cycle can easily be replaced with the apply
function.