[R] Predict polynomial problem

[Charles C. Berry]

and the values in those 
places are different:

On Tue, 19 Jan 2010, 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...
>

Per ?bquote, "bquote quotes its argument except that terms wrapped in 
'.()' are evaluated in the specified 'where' environment.

(which by default is the parent.frame)

Note:

> i <- 20
> bquote(y ~ poly(x,.(i)))
y ~ poly(x, 20)
>

So, now 'i' is irrelevant as the expression returned by bquote has '20' as 
the 'degree' arg.



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

Well, the terms() objects are the same:

> all.equal(terms(m1),terms(m2))
[1] TRUE
>

but they will look in different places for 'i':


> environment(terms(m2))
<environment: 0x01b7c178>
> environment(terms(m1))
<environment: R_GlobalEnv>
>

and the values in those places are different:

> environment(terms(m2))$i
[1] 2
> environment(terms(m1))$i
[1] 3
>


> #3 works, gets 'i' from within:
> predict(m3,newdata=data.frame(x=9:11))

It doesn't need 'i', because the i was evaluated and substituted by 
bquote. That is, it doesn't get("i") as the expression returned by bquote 
has no 'i' in it.

HTH,

Chuck

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

Charles C. Berry                            (858) 534-2098
                                             Dept of Family/Preventive Medicine
E mailto:cberry at tajo.ucsd.edu	            UC San Diego
http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901