Re: [R] Fitting models in a loop

From: Gabor Grothendieck <ggrothendieck_at_gmail.com>
Date: Wed 02 Aug 2006 - 13:57:37 EST

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@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@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@csiro.au <Bill.Venables@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@stat.math.ethz.ch
> > > > [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of
> > > > Bill.Venables@csiro.au
> > > > Sent: 01 August 2006 06:16
> > > > To: maj@waikato.ac.nz; r-help@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@stat.math.ethz.ch
> > > > [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of Murray
> > > Jorgensen
> > > > Sent: Tuesday, 1 August 2006 2:09 PM
> > > > To: r-help@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@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.
> > >
> >
>



R-help@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. Received on Wed Aug 02 14:01:27 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Wed 02 Aug 2006 - 16:19:21 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.