[R] Predict polynomial problem

[Barry Rowlingson]

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

 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

-- 
blog: http://geospaced.blogspot.com/
web: http://www.maths.lancs.ac.uk/~rowlings
web: http://www.rowlingson.com/
twitter: http://twitter.com/geospacedman
pics: http://www.flickr.com/photos/spacedman