[R] Fitting models in a loop

[Bill.Venables at csiro.au]

This (below) also runs into trouble if you try to predict with new data
since you have no rule for re-constructing the formula.  Also, you can't
plot the term as a single contributor to the linear predictor with
termplot().

I'm sure given enough ingenuity you can get round these two, but why
avoid the language manipulation solution, when it does the lot?

Bill.


-----Original Message-----
From: Gabor Grothendieck [mailto:ggrothendieck at gmail.com] 
Sent: Wednesday, 2 August 2006 12:01 PM
To: Venables, Bill (CMIS, Cleveland)
Cc: Markus.Gesmann at lloyds.com; maj at waikato.ac.nz;
r-help at stat.math.ethz.ch
Subject: Re: [R] Fitting models in a loop

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