[R] Fitting models in a loop

[Gabor Grothendieck]

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.
>