R help - Jan 2014 - Safe prediction does not work for bivariate polynomial terms?

Hi everyone,

R documents says the safe prediction is implemented, when basis
functions are used, such as poly(), bs(), ns()

This works for univariate basis, but fails in my bivariate polynomial setting.
Can anyone explain the reason?


The following is a small example.

set.seed(731)
x<-runif(300)
y<-runif(300)

f <- function(x, y) {  0.5+2*x+13*x^2-14*y^2+12*x*y+y }

z <- f(x,y)+rnorm(length(x))*0.2

# default orthogonal polynomials basis
mod <- lm (z ~ poly(x,y,degree = 2))

# raw polynomials basis
mod1 <- lm (z ~ poly(x,y,degree = 2, raw = T))

# data points to evaluate, just the first five data
new <- data.frame(x=x[1:5],y= y[1:5])

z[1:5]
[1]  9.796620 10.397851  1.280832  4.028284  4.811709

# two predicted values differ dramatically, and the orthogonal
polynomials basis fails
predict(mod, new)
        1         2         3         4         5
121.46776  40.85002  18.67273  31.82417  20.81673
predict(mod1, new)
        1         2         3         4         5
 9.981091 10.418628  1.223148  4.031664  4.837099

# However, the fitted.values are the same
mod$fitted.values[1:5]
        1         2         3         4         5
 9.981091 10.418628  1.223148  4.031664  4.837099
mod1$fitted.values[1:5]
        1         2         3         4         5
 9.981091 10.418628  1.223148  4.031664  4.837099


Thanks in advance
Xing

On 24/01/2014 02:09, Xing Zhao wrote:
> Hi everyone,
>
> R documents says the safe prediction is implemented, when basis
> functions are used, such as poly(), bs(), ns()
>
> This works for univariate basis, but fails in my bivariate polynomial
setting.
> Can anyone explain the reason?
Because there is a makepredictcall() method for class "poly", and 
bivariate polynomials are not of that class (see ?poly).  The 
documentation only says safe prediction is available for 'polynomial' 
fits, and 'bivariate polynomial' is not conventionally included in 
'polynomial'.

Your call is really to polym(), not to poly(), and it may be better to 
call polym() explicitly to remind yourself.
>
>
> The following is a small example.
>
> set.seed(731)
> x<-runif(300)
> y<-runif(300)
>
> f <- function(x, y) {  0.5+2*x+13*x^2-14*y^2+12*x*y+y }
>
> z <- f(x,y)+rnorm(length(x))*0.2
>
> # default orthogonal polynomials basis
> mod <- lm (z ~ poly(x,y,degree = 2))
>
> # raw polynomials basis
> mod1 <- lm (z ~ poly(x,y,degree = 2, raw = T))
>
> # data points to evaluate, just the first five data
> new <- data.frame(x=x[1:5],y= y[1:5])
>
> z[1:5]
> [1]  9.796620 10.397851  1.280832  4.028284  4.811709
>
> # two predicted values differ dramatically, and the orthogonal
> polynomials basis fails
> predict(mod, new)
>          1         2         3         4         5
> 121.46776  40.85002  18.67273  31.82417  20.81673
> predict(mod1, new)
>          1         2         3         4         5
>   9.981091 10.418628  1.223148  4.031664  4.837099
>
> # However, the fitted.values are the same
> mod$fitted.values[1:5]
>          1         2         3         4         5
>   9.981091 10.418628  1.223148  4.031664  4.837099
> mod1$fitted.values[1:5]
>          1         2         3         4         5
>   9.981091 10.418628  1.223148  4.031664  4.837099
>
>
> Thanks in advance
> Xing
>
> ______________________________________________
> 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.
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595