Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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?

  • This post has 4 questions total, highlighted below.
  • I'm seeking to evaluate the output of a poly()-output object (orthogonal or raw polynomials). Such would give me the ability to use polynomials to generate a row similar to that of my model matrix that I may use to evaluate a result (i.e., I seek to push multivariate test data values through a poly() call so that it can be evaluated similarly to a row of my regression's method matrix).
  • My background: I'm relatively new to R, R's poly(), and R's regression routines.
  • I've tried several approaches and would appreciate help with each:

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

like image 234
blake Avatar asked Nov 09 '22 10:11

blake


1 Answers

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)
like image 141
Benjamin Christoffersen Avatar answered Nov 15 '22 06:11

Benjamin Christoffersen