[R] Predict polynomial problem

Gavin Simpson gavin.simpson at ucl.ac.uk
Tue Jan 19 15:37:43 CET 2010


On Tue, 2010-01-19 at 08:27 +0000, Barry Rowlingson wrote:
> On Tue, Jan 19, 2010 at 1:36 AM, Charles C. Berry <cberry at tajo.ucsd.edu> wrote:
> 
> > Its the environment thing.
> >
> > I think you want something like this:
> >
> >        models[[i]]=lm( bquote( y ~ poly(x,.(i)) ), data=d)
> >
> > Use
> >        terms( mmn[[3]] )
> >
> > both with and without this change and
> >
> >
> >        ls( env = environment( formula( mmn[[3]] ) ) )
> >        get("i",env=environment(formula(mmn[[3]])))
> >        sapply(mmn,function(x) environment( formula( x ) ) )
> >
> >
> > to see what gives.
> 
>  Think I see it now. predict involves evaluating poly, and poly here
> needs 'i' for the order. If the right 'i' isn't gotten when predict is
> called then I get the error. Your fix sticks the right 'i' into the
> environment when predict is called.
> 
>  I haven't quite got my head round _how_ it does it, and I have no
> idea how I could have figured this out for myself. Oh well...

Perhaps this Programmer's Niche article by Bill Venables might also be
useful as it discuss how to manipulate formulas to automate model
fitting...?

Bill Venables. Programmer's Niche. R News, 2(2):24-26, June 2002.

http://cran.r-project.org/doc/Rnews/Rnews_2002-2.pdf

HTH

G

> 
>  The following lines are also illustrative:
> 
> d = data.frame(x=1:10,y=runif(10))
> 
> i=3
> #1 naive model:
> m1 = lm(y~poly(x,i),data=d)
> #2,3 bquote, without or with i-wrapping:
> m2 = lm(bquote(y~poly(x,i)),data=d)
> m3 = lm(bquote(y~poly(x,.(i))),data=d)
> 
> #1 works, gets 'i' from global i=3 above:
> predict(m1,newdata=data.frame(x=9:11))
> #2 fails - why?
> predict(m2,newdata=data.frame(x=9:11))
> #3 works, gets 'i' from within:
> predict(m3,newdata=data.frame(x=9:11))
> 
> rm(i)
> 
> #1 now fails because we removed 'i' from top level:
> predict(m1,newdata=data.frame(x=9:11))
> #2 still fails:
> predict(m2,newdata=data.frame(x=9:11))
> #3 still works:
> predict(m3,newdata=data.frame(x=9:11))
> 
> Thanks
> 

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%



More information about the R-help mailing list