[R] Fitting models in a loop

[Gabor Grothendieck]

Here is another attempt.  This one allows general prediction
yet its actually shorter and does not use any advanced
language constructs (although to understand why it works
one must understand that formulas have environments
and the environment of the formula corresponding to each
component of mod is the environment within the anonymous
function instance that created it):

# test data - as before
set.seed(1)
x <- 1:10
y <- x^3 + rnorm(10)

mod <- lapply(1:3, function(i) lm(y ~ poly(x,i)))
print(mod)

# test - each component of mod remembers its 'i'
# This returns 1, 2 and 3 as required.
for (j in 1:3) print(environment(formula(mod[[j]]))$i)

# following two lines give same answer
# showing prediction works
predict(mod[[2]], list(x = 1:10))
fitted(lm(y ~ poly(x,2)))

On 8/1/06, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:
> Actually in thinking about this some more that still gets you
> into a mess if you want to do prediction at anything other
> than the original points.
>
> On 8/1/06, Gabor Grothendieck <ggrothendieck at gmail.com> wrote:
> > A simple way around this is to pass it as a data frame.
> > In the code below the only change we made was to change
> > the formula from y ~ poly(x, i) to y ~ . and pass poly(x,i)
> > in a data frame as argument 2 of lm:
> >
> > # test data
> > set.seed(1)
> > x <- 1:10
> > y <- x^3 + rnorm(10)
> >
> > # run same code except change the lm call
> > mod <- list()
> > for (i in 1:3) {
> >        mod[[i]] <- lm(y ~., data.frame(poly(x, i)))
> >        print(summary(mod[[i]]))
> > }
> >
> > After running the above we can test that it works:
> >
> > > for(i in 1:3) print(formula(mod[[i]]))
> > y ~ X1
> > y ~ X1 + X2
> > y ~ X1 + X2 + X3
> >
> > On 8/1/06, Bill.Venables at csiro.au <Bill.Venables at csiro.au> wrote:
> > >
> > > Markus Gesmann writes:
> > >
> > > > Murray,
> > > >
> > > > How about creating an empty list and filling it during your loop:
> > > >
> > > >  mod <- list()
> > > >  for (i in 1:6) {
> > > >         mod[[i]] <- lm(y ~ poly(x,i))
> > > >         print(summary(mod[[i]]))
> > > >         }
> > > >
> > > > All your models are than stored in one object and you can use lapply
> > > to
> > > > do something on them, like:
> > > >  lapply(mod, summary) or lapply(mod, coef)
> > >
> > > I think it is important to see why this deceptively simple
> > > solution does not achieve the result that Murray wanted.
> > >
> > > Take any fitted model object, say mod[[4]].  For this object the
> > > formula component of the call will be, literally,  y ~ poly(x, i),
> > > and not y ~ poly(x, 4), as would be required to use the object,
> > > e.g. for prediction.  In fact all objects have the same formula.
> > >
> > > You could, of course, re-create i and some things would be OK,
> > > but getting pretty messy.
> > >
> > > You would still have a problem if you wanted to plot the fit with
> > > termplot(), for example, as it would try to do a two-dimensional
> > > plot of the component if both arguments to poly were variables.
> > >
> > > >
> > > > -----Original Message-----
> > > > From: r-help-bounces at stat.math.ethz.ch
> > > > [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of
> > > > Bill.Venables at csiro.au
> > > > Sent: 01 August 2006 06:16
> > > > To: maj at waikato.ac.nz; r-help at stat.math.ethz.ch
> > > > Subject: Re: [R] Fitting models in a loop
> > > >
> > > >
> > > > Murray,
> > > >
> > > > Here is a general paradigm I tend to use for such problems.  It
> > > extends
> > > > to fairly general model sequences, including different responses, &c
> > > >
> > > > First a couple of tiny, tricky but useful functions:
> > > >
> > > > subst <- function(Command, ...) do.call("substitute", list(Command,
> > > > list(...)))
> > > >
> > > > abut <- function(...)  ## jam things tightly together
> > > >   do.call("paste", c(lapply(list(...), as.character), sep = ""))
> > > >
> > > > Name <- function(...) as.name(do.call("abut", list(...)))
> > > >
> > > > Now the gist.
> > > >
> > > > fitCommand <- quote({
> > > >         MODELi <- lm(y ~ poly(x, degree = i), theData)
> > > >         print(summary(MODELi))
> > > > })
> > > > for(i in 1:6) {
> > > >         thisCommand <- subst(fitCommand, MODELi = Name("model_", i), i
> > > =
> > > > i)
> > > >         print(thisCommand)  ## only as a check
> > > >         eval(thisCommand)
> > > > }
> > > >
> > > > At this point you should have the results and
> > > >
> > > > objects(pat = "^model_")
> > > >
> > > > should list the fitted model objects, all of which can be updated,
> > > > summarised, plotted, &c, because the information on their construction
> > > > is all embedded in the call.
> > > >
> > > > Bill.
> > > >
> > > > -----Original Message-----
> > > > From: r-help-bounces at stat.math.ethz.ch
> > > > [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Murray
> > > Jorgensen
> > > > Sent: Tuesday, 1 August 2006 2:09 PM
> > > > To: r-help at stat.math.ethz.ch
> > > > Subject: [R] Fitting models in a loop
> > > >
> > > > If I want to display a few polynomial regression fits I can do
> > > something
> > > >
> > > > like
> > > >
> > > > for (i in 1:6) {
> > > >         mod <- lm(y ~ poly(x,i))
> > > >         print(summary(mod))
> > > >         }
> > > >
> > > > Suppose that I don't want to over-write the fitted model objects,
> > > > though. How do I create a list of blank fitted model objects for later
> > >
> > > > use in a loop?
> > > >
> > > > Murray Jorgensen
> > > > --
> > >
> > > ______________________________________________
> > > R-help at stat.math.ethz.ch mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> > > and provide commented, minimal, self-contained, reproducible code.
> > >
> >
>