[R] Fitting models in a loop

[Gabor Grothendieck]

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