How to run lm models using all possible combinations of several variables and a factor

Question

this is not my subject so I am sorry if my question is badly asked or if the data is incomplete. I am trying to run 31 lineal models which have a single response variable (VELOC), and as predictor variables have a factor (TRAT, with 2 levels, A and B) and five continuous variables.

I have a loop I used for gls but only with continuous predictor variables, so I thought it could work. But it did not and I believe the problem must be a silly thing.

I don't know how to include the data, but it looks something like:

   TRAT    VELOC        l        b        h        t        m

1     A  0.02490 -0.05283  0.06752  0.03435 -0.03343  0.10088

2     A  0.01196 -0.01126  0.10604 -0.01440 -0.08675  0.18547

3     A -0.06381  0.00804  0.06248 -0.04467 -0.04058 -0.04890

4     A  0.07440  0.04800  0.05250 -0.01867 -0.08034  0.08049

5     A  0.07695  0.06373 -0.00365 -0.07319 -0.02579  0.06989

6     B -0.03860 -0.01909  0.04960  0.09184 -0.06948  0.17950

7     B  0.00187 -0.02076 -0.05899 -0.12245  0.12391 -0.25616

8     B -0.07032 -0.02354 -0.05741  0.03189  0.05967 -0.06380

9     B -0.09047 -0.06176 -0.17759  0.15136  0.13997  0.09663

10    B -0.01787  0.01665 -0.08228 -0.02875  0.07486 -0.14252

now, the script I used is:

pred.vars = c("TRAT","l", "b", "h","t","m") #define predictor variables

m.mat = permutations(n = 2, r = 6, v = c(F, T), repeats.allowed = T)# I run        all possible combinations of pred.vars
models = apply(cbind(T, m.mat), 1, function(xrow) {paste(c("1", pred.vars)
[xrow], collapse = "+")})# fill the models 
models = paste("VELOC", models, sep = "~")#fill the left side
all.aic = rep(NA, length(models))# AIC of models
m.mat = cbind(1, m.mat)# Which predictors are estimated in the models beside
#the intercept

colnames(m.mat) = c("(Intercept)", pred.vars)

n.par = 2 + apply(m.mat,1, sum)# number of parameters estimated in the Models
coefs=m.mat# define an object to store the coefficients 

for (k in 1:length(models)) {
   res = try(lm(as.formula(models[k]), data = xdata))
   if (class(res) != "try-error") {
   all.aic[k] = -2 * logLik(res)[1] + 2 * n.par[k]
   xx = coefficients(res)
   coefs[k, match(names(xx), colnames(m.mat))] = xx
    }
}

And I get this error:"Error in coefs[k, match(names(xx), colnames(m.mat))] = xx : NAs are not allowed in subscripted assignments"

Thanks in advance for your help. I'll appreciate any corrections about how to post properly questions. Lina

Answer 1

I suspect the dredge function in the MuMIn package would help you. You specify a "full" model with all parameters you want to include and then run dredge(fullmodel) to get all combinations nested within the full model.

You should then be able to get the coefficients and AIC values from the results of this.

Something like:

require(MuMIn)
data(iris)

globalmodel <- lm(Sepal.Length ~ Petal.Length + Petal.Width + Species, data = iris)

combinations <- dredge(globalmodel)

print(combinations)

to get the parameter estimates for all models (a bit messy) you can then use

coefTable(combinations)

or to get the coefficients for a particular model you can index that using the row number in the dredge object, e.g.

coefTable(combinations)[1]

to get the coefficients in the model at row 1. This should also print coefficients for factor levels.

See the MuMIn helpfile for more details and ways to extract information.

Hope that helps.

Answer 2

To deal with:

'global.model''s 'na.action' argument is not set and options('na.action') is "na.omit"

require(MuMIn)
data(iris)
options(na.action = "na.fail") # change the default "na.omit" to prevent models
# from being fitted to different datasets in
# case of missing values.
globalmodel <- lm(Sepal.Length ~ Petal.Length + Petal.Width + Species, data = iris)
combinations <- dredge(globalmodel)
print(combinations)