CODEWITHSUNDEEP
rregressionpolynomialsorthogonal
blake
blake

Reputation: 119

How do you make R poly() evaluate (or "predict") multivariate new data (orthogonal or raw)?

With the poly function in R, how do I evaluate a multivariate polynomial?

(A): Direct approach with predict

This method FAILED, apparently due to some unexpected class of the inputs. I'm aware that these particular x1 & x2 values, being collinear, are not ideal for a general fit (I'm simply trying to get the predict machinery operational). The use of predict was inspired by this SO post. (Q1) Is it possible to call the predict method directly to evaluate this polynomial?

> x1 = seq(1,  10,  by=0.2)
> x2 = seq(1.1,10.1,by=0.2)
> t = poly(cbind(x1,x2),degree=2,raw=T)
> predict(t,newdata=data.frame(x1=2.03,x2=2.03))
Error in UseMethod("predict") : 
  no applicable method for 'predict' applied to an object of class "c('matrix', 'double', 'numeric')"

(B) Direct evaluation only works for raw polynomials (not orthogonal)

Due to (A), I tried a workaround with a direct call to poly(). For raw polynomials, I could get it to work, but I had to repeat data for each respective variable. Following shows (1st) the failure with single data point, (2nd) success with repeating the value. (Q2) Is there any way to avoid the redundant repeat of data in the 2nd listing to make raw poly() evaluate properly?

> poly(cbind(x1=c(2.03),x2=c(2.13)),degree=2,raw=T)
Error in `colnames<-`(`*tmp*`, value = apply(z, 1L, function(x) paste(x,  : 
  attempt to set 'colnames' on an object with less than two dimensions

> poly(cbind(x1=c(2.03,2.03),x2=c(2.13,2.13)),degree=3,raw=T)
      1.0    2.0      3.0  0.1    1.1      2.1    0.2      1.2      0.3
[1,] 2.03 4.1209 8.365427 2.13 4.3239 8.777517 4.5369 9.209907 9.663597
[2,] 2.03 4.1209 8.365427 2.13 4.3239 8.777517 4.5369 9.209907 9.663597
attr(,"degree")
[1] 1 2 3 1 2 3 2 3 3

If I try a similar redundantly-listed-data approach with orthogonal polynomials, I incur the "hey, your data's redundant!" error (which I also incur if I only list each variable's value once). (Q3) Is it possible to evaluate multivariate orthogonal polynomials via a direct call to poly()? 

> poly(cbind(x1=c(2.03, 2.03),x2=c(2.13, 2.13)),degree=2)
Error in poly(dots[[1L]], degree, raw = raw) : 
  'degree' must be less than number of unique points

(C) Inability to extract alpha & norm coefficients from multivariate orthogonal polynomials Finally, I'm aware that there is a coefs input variable to predict.poly. I understand the coefs are to be the alpha and norm values output from an orthogonal polynomial fit. However, I am only able to extract those from a univariate polynomial fit ... when I fit a multivariate orthogonal (or raw), the return value from poly does NOT have the coefs. (Q4) Is it possible to extract alpha and norm coefficients from calling poly() for an orthogonal polynomial fit to multivariate data?

> t = poly(cbind(x1),degree=2)   # univariate orthog poly --> WORKS
> attributes(t)$coefs
$alpha
[1] 5.5 5.5

$norm2
[1]    1.000   46.000  324.300 1826.458


> t = poly(cbind(x1,x2),degree=2)  # multivariate orthog poly --> DOES NOT WORK
> attributes(t)$coefs
NULL

Please let me know if I can clarify. Thank you kindly for any help you can provide.

Upvotes: 5

Views: 2960

Answers (2)

Ben2018
Ben2018

Reputation: 635

Use ModelMatrixModel package. ModelMatrixModel() in is similar to model.matrix(), but it save the transforming parameters, which can be applied new data.

library(ModelMatrixModel) 
traindf=data.frame(x1 = seq(1,  10,  by=0.2),
                   x2 = seq(1.1,10.1,by=0.2))

mm=ModelMatrixModel(~poly(x1,2)+poly(x2,3),traindf,sparse=F)
mm$x[1:2,] #output matrix
##   poly_x1__2_1 poly_x1__2_2 poly_x2__3_1 poly_x2__3_2 poly_x2__3_3
## 1   -0.2498843    0.3088653   -0.2498843    0.3088653   -0.3423492
## 2   -0.2387784    0.2676833   -0.2387784    0.2676833   -0.2510561
predict(mm,traindf[1:2,])$x
##   poly_x1__2_1 poly_x1__2_2 poly_x2__3_1 poly_x2__3_2 poly_x2__3_3
## 1   -0.2498843    0.3088653   -0.2498843    0.3088653   -0.3423492
## 2   -0.2387784    0.2676833   -0.2387784    0.2676833   -0.2510561

Upvotes: 0

Benjamin Christoffersen
Benjamin Christoffersen

Reputation: 4841

For the record, it seems that this has been fixed

> x1 = seq(1,  10,  by=0.2)
> x2 = seq(1.1,10.1,by=0.2)
> t = poly(cbind(x1,x2),degree=2,raw=T)
> 
> class(t) # has a class now
[1] "poly"   "matrix"
> 
> # does not throw error
> predict(t, newdata = cbind(x1,x2)[1:2, ])                                                     
     1.0  2.0 0.1  1.1  0.2
[1,] 1.0 1.00 1.1 1.10 1.21
[2,] 1.2 1.44 1.3 1.56 1.69
attr(,"degree")
[1] 1 2 1 2 2
attr(,"class")
[1] "poly"   "matrix"
> 
> # and gives the same
> t[1:2, ]
     1.0  2.0 0.1  1.1  0.2
[1,] 1.0 1.00 1.1 1.10 1.21
[2,] 1.2 1.44 1.3 1.56 1.69
> 
> sessionInfo()
R version 3.4.1 (2017-06-30)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Upvotes: 1

Related Questions