9

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

NicE
  • 21,165
  • 3
  • 51
  • 68
Lina Moreno
  • 101
  • 1
  • 1
  • 3

2 Answers2

15

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.

Adam Kimberley
  • 879
  • 5
  • 12
  • Thanks for your suggestion, Adam! I run it and apparently it worked, but I get no estimate values for the factor, other than that, I got the AIC, R2 values, etc. for all the possible models... Im not sure if that is correct... – Lina Moreno Feb 20 '15 at 15:09
  • I've edited what I think should help get the estimates for factor variables you're after. I'd suggest reading the linked helpfile if you want more details on what exactly is produced in the dredge output. – Adam Kimberley Feb 20 '15 at 15:33
  • Thank you again Adam, I totally forgot to ask for the str(combination) so to see what I can extract.. Now the problem is solve, thanks! – Lina Moreno Feb 23 '15 at 20:54
  • 6
    Error in dredge(globalmodel) : 'global.model''s 'na.action' argument is not set and options('na.action') is "na.omit" – Mox Jul 13 '18 at 21:57
  • With this line `coefTable(combinations)`, do you know how to also add p-values or t-stats to the coefTable? Thank you – Jeremy K. Jul 08 '19 at 16:58
1

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)
Mox
  • 511
  • 5
  • 15
  • 8
    This is not perfect, as it changes an R global option. It would be better to include this in the model function, such as `lm(..., na.action = "na.fail")`. See https://stackoverflow.com/questions/25281739/dredge-function-error-r-package-mumln – JASC Jan 10 '19 at 04:37