CODEWITHSUNDEEP
rregressionlinear-regressionlm
mchen
mchen

Reputation: 10166

How do I use lm() from inside a function?

It seems that calling lm() from within a function or via lapply screws up the $call associated with a fit. Minimal working example:

> library(MASS)
> dat <- data.frame(x = 1:100, y=1:100)
> dat <- within(dat, z <- x + log(y) + rnorm(100))
> fits <- lapply(list(z ~ x + y, z ~ x + log(y)), lm, dat)
> stepAIC(fits[[1]])               # <-- error when I try to use the fit in other functions
Error in eval(expr, envir, enclos) : could not find function "FUN"

> fits[[1]]$call
FUN(formula = X[[i]], data = ..1)  # Aha -- this must be why -- $call is screwed up

How do I resolve this problem and prevent the above error?

Upvotes: 1

Views: 1208

Answers (3)

G. Grothendieck
G. Grothendieck

Reputation: 270298

Try using this as the function in the lapply. This produces nice looking formulas in fits showing the actual formula and stepAIC works:

fun <- function(fo) do.call("lm", list(fo, quote(dat)))
fits <- lapply(list(z ~ x + y, z ~ x + log(y)), fun)

giving:

> fits[[1]]

Call:
lm(formula = z ~ x + y, data = dat)

Coefficients:
(Intercept)            x            y  
      2.154        1.031           NA  


> stepAIC(fits[[1]])
Start:  AIC=-3.34
z ~ x + y


Step:  AIC=-3.34
z ~ x

       Df Sum of Sq   RSS    AIC
<none>                 93  -3.34
- x     1     88600 88693 680.78

Call:
lm(formula = z ~ x, data = dat)

Coefficients:
(Intercept)            x  
      2.154        1.031

Upvotes: 3

coffeinjunky
coffeinjunky

Reputation: 11514

An alternative yould be to apply the stepAIC inside lapply directly:

 AICs <- lapply(list(z ~ x + y, z ~ x + log(y)), 
                function(x) stepAIC(lm(x,dat)))

This gives you a list of the output of stepAIC for all models.

Upvotes: 2

Roland
Roland

Reputation: 132989

Sometimes it's better to supply an anonymous function to lapply:

fits <- lapply(list(z ~ x + y, z ~ x + log(y)), 
                   function(f) lm(f, data = dat))
stepAIC(fits[[1]])
#works

Note that this (usually my preferred way to make scoping explicit) doesn't work because DF is not found by stepAIC:

fits <- lapply(list(z ~ x + y, z ~ x + log(y)), 
                   function(f, DF) lm(f, data = DF), DF = dat)

Note that stepwise regression is a bad method anyway.

Upvotes: 6

Related Questions