Is there an R library that estimates a multivariate natural cubic spline (or similar) function?

note: originally posted on Cross Validated (stats SE) on 07-26-2011, with no correct answers to date.

Background

I have a model, f, where Y=f(X)

X is an n x m matrix of samples from m parameters and Y is the n x 1 vector of model outputs.

f is computationally intensive, so I would like to approximate f using a multivariate cubic spline through (X,Y) points, so that I can evaluate Y at a larger number of points.

Question

Is there an R function that will calculate an arbitrary relationship between X and Y?

Specifically, I am looking for a multivariate version of the splinefun function, which generates a spline function for the univariate case.

e.g. this is how splinefun works for the univariate case

x <- 1:100
y <- runif(100)
foo <- splinefun(x,y, method = "monoH.FC")
foo(x)开发者_StackOverflow #returns y, as example

The test that the function interpolates exactly through the points is successful:

all(y == foo(1:100))
## TRUE

What I have tried

I have reviewed the mda package, and it seems that the following should work:

library(mda)
x   <- data.frame(a = 1:100, b = 1:100/2, c = 1:100*2)
y   <- runif(100)
foo <- mars(x,y)
predict(foo, x) #all the same value

however the function does not interpolate exactly through the design points:

all(y == predict(foo,x))
## FALSE

I also could not find a way to implement a cubic-spline in either the gam, marss, or earth packages.

Actually several packages can do it. The one I use is the "rms" package which has rcs, but the survival package also has pspline and the splines package has the ns function {}. "Natural splines" (constructed with ns) are also cubic splines. You will need to form multivariate fitting function with the '*' operator in the multivariate formula creating "crossed" spline terms. that the example you offered was not sufficiently rich.

I guess I am confused that you want exact fits. R is a statistical package. Approximate estimation is the goal. Generally exact fits are more of a problem because they lead to multicollinearity.

Have a look at the DiceKriging package which was developed to undertake tasks like this. http://cran.r-project.org/web/packages/DiceKriging/index.html

I've provided an example application at https://stats.stackexchange.com/questions/13510/fitting-multivariate-natural-cubic-spline/65012#65012

I'm not sure if this is precisely what you are looking for, but you could try Tps() in the R package fields. It's meant for doing thin-plate splines interpolations (2D equivalent of cubic splines) for spatial data, but will take up to four covariates, although it will expect them to be euclidean x,y,z + time, so you need to be clear that you are selecting the correct options for your particular case. If you want to interpolate, set the smoothing parameter lambda to zero. You might also try the function polymars() in the R package polspline.