[R] Predict polynomial problem

[Charles C. Berry]

On Tue, 19 Jan 2010, Charles C. Berry wrote:

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

And I should have mentioned that environment(terms(m2)) happens to have an 
object 'i' in it regardless of whether poly() is used.

Chuck


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

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