[R] using poly() to predict

[Stanislav Aggerwal]

Thanks Prof Ripley.
For anybody else wondering about this, see:
http://stackoverflow.com/questions/26728289/extracting-orthogonal-polynomial-coefficients-from-rs-poly-function

=========================

The polynomials are defined recursively using the alpha and norm2
coefficients of the poly object you've created. Let's look at an example:

z <- poly(1:10, 3)
attributes(z)$coefs# $alpha# [1] 5.5 5.5 5.5# $norm2# [1]    1.0
10.0   82.5  528.0 3088.8

For notation, let's call a_d the element in index d of alpha and let's call
n_d the element in index d of norm2. F_d(x) will be the orthogonal
polynomial of degree d that is generated. For some base cases we have:

F_0(x) = 1 / sqrt(n_2)
F_1(x) = (x-a_1) / sqrt(n_3)

The rest of the polynomials are recursively defined:

F_d(x) = [(x-a_d) * sqrt(n_{d+1}) * F_{d-1}(x) - n_{d+1} / sqrt(n_d) *
F_{d-2}(x)] / sqrt(n_{d+2})

To confirm with x=2.1:

x <- 2.1
predict(z, newdata=x)#               1         2         3# [1,]
-0.3743277 0.1440493 0.1890351# ...

a <- attributes(z)$coefs$alpha
n <- attributes(z)$coefs$norm2
f0 <- 1 / sqrt(n[2])(f1 <- (x-a[1]) / sqrt(n[3]))# [1] -0.3743277(f2
<- ((x-a[2]) * sqrt(n[3]) * f1 - n[3] / sqrt(n[2]) * f0) /
sqrt(n[4]))# [1] 0.1440493(f3 <- ((x-a[3]) * sqrt(n[4]) * f2 - n[4] /
sqrt(n[3]) * f1) / sqrt(n[5]))# [1] 0.1890351

The most compact way to export your polynomials to your C++ code would
probably be to export attributes(z)$coefs$alpha and
attributes(z)$coefs$norm2 and then use the recursive formula in C++ to
evaluate your polynomials.


On Wed, Jan 14, 2015 at 2:38 PM, Prof Brian Ripley <ripley at stats.ox.ac.uk>
wrote:

> On 14/01/2015 14:20, Stanislav Aggerwal wrote:
>
>> This method of finding yhat as x %*% b works when I use raw polynomials:
>>
>> x<-1:8
>> y<- 1+ 1*x + .5*x^2
>> fit<-lm(y~poly(x,2,raw=T))
>> b<-coef(fit)
>> xfit<-seq(min(x),max(x),length=20)
>> yfit<-b[1] + poly(xfit,2,raw=T) %*% b[-1]
>> plot(x,y)
>> lines(xfit,yfit)
>>
>> But it doesn't work when I use orthogonal polynomials:
>>
>> fit<-lm(y~poly(x,2))
>> b<-coef(fit)
>> yfit<-b[1] + poly(xfit,2) %*% b[-1]
>> plot(x,y)
>> lines(xfit,yfit,col='red')
>>
>> I have a feeling that the second version needs to incorporate poly() coefs
>> (alpha and norm2) somehow. If so, please tell me how.
>>
>> I do know how to use predict() for this. I just want to understand how
>> poly() works.
>>
>
> What matters is how lm() and predict() use poly(): see ?makepredictcall
> and its code.
>
> str(fit) might also help.
>
>
>> Thanks very much for any help
>> Stan
>>
>>         [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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.
>>
>>  Please do, and do not send HTML.
>
> --
> Brian D. Ripley,                  ripley at stats.ox.ac.uk
> Emeritus Professor of Applied Statistics, University of Oxford
> 1 South Parks Road, Oxford OX1 3TG, UK
>

	[[alternative HTML version deleted]]