[R] Fitting models in a loop

[Bill.Venables at csiro.au]

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