[R] Predict polynomial problem

[Peter Ehlers]

Charles C. Berry wrote:
> 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.
> 
Right. It might be worth pointing out that the 'i' in
environment(terms(m2)) or in environment(terms(m3)),
which also has an 'i', is there even if you use poly(x,j).
It is (if I'm not mistaken) the number of 'variables' in
the formula: response plus predictor terms. Thus

j <- 5
m4 <- lm(bquote(y ~ sqrt(x) + poly(x, .(j))), data=d)
environment(terms(m4))$i
[1] 3

  -Peter Ehlers

> 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
> 
> 
> ------------------------------------------------------------------------
> 
> ______________________________________________
> R-help at r-project.org 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.

-- 
Peter Ehlers
University of Calgary
403.202.3921