Obtain spline surface on R

Question

How do I generate a b-spline surface, let's say:

x=attitude$rating
y=attitude$complaints
z=attitude$privileges

would be x and y for the spline basis. z is the set of control points.

Arbela

It's a standard dataset of R – Arbela

+1 for providing a dataset. – Meggs

Answer 1

If I understand you, you have x,y, and z data and you want to use bivariate spline interpolation on x and y, using z for the control points. You can do this with interp(...) in the akima package.

library(akima)
spline <- interp(x,y,z,linear=FALSE)
# rotatable 3D plot of points and spline surface
library(rgl)
open3d(scale=c(1/diff(range(x)),1/diff(range(y)),1/diff(range(z))))
with(spline,surface3d(x,y,z,alpha=.2))
points3d(x,y,z)
title3d(xlab="rating",ylab="complaints",zlab="privileges")
axes3d()

The plot itself is fairly uninteresting with your dataset because x, y, and x are highly correlated.

EDIT response to OP's comment.

If you want a b-spline surface, try out mba.surf(...) in the unfortunately named MBA package.

library(MBA)
spline <- mba.surf(data.frame(x,y,z),100,100)

library(rgl)
open3d(scale=c(1/diff(range(x)),1/diff(range(y)),1/diff(range(z))))
with(spline$xyz,surface3d(x,y,z,alpha=.2))
points3d(x,y,z)
title3d(xlab="rating",ylab="complaints",zlab="privileges")
axes3d()

Meggs

Ok, I forgot to mention I was looking for a b-spline surface – Arbela

@user3083324: Seems like a very useful response even if jhoward did not read your mind. – Torus

from what I could deduce, the parameter to control degree is h – Arbela

Read the documentation on mba.surf(...). In the example, apline$xyz$z is a matrix of interpolated s-values. – Meggs

Answer 2

 require(rms)  # Harrell's gift to the R world.

 # Better to keep the original names and do so within a dataframe.
 att <- attitude[c('rating','complaints','privileges')]
 add <- datadist(att)  # records ranges and descriptive info on data
 options(datadist="add")  # need these for the rms functions

#  rms-`ols` function (ordinary least squares) is a version of `lm`
 mdl <- ols( privileges ~ rcs(rating,4)*rcs(complaints,4) ,data=att)
# Predict is an rms function that works with rms's particular classes
 pred <- Predict(mdl, 'rating','complaints')
# bplot calls lattice functions; levelplot by default; this gives a "3d" plot
 bplot(pred, yhat~rating+complaints, lfun=wireframe)

It's a crossed restricted-cubic spline model. If you have a favorite spline function you want to use instead, then by all means try it out. I've had good luck with the rcs- function.

This gives a more open mesh with fewer calculated points:

pred <- Predict(mdl, 'rating','complaints', np=25)
bplot(pred, yhat~rating+complaints, lfun=wireframe)
png()
bplot(pred, yhat~rating+complaints, lfun=wireframe)
dev.off()

You could use the rgl methods being illustrated by jhoward. The top of str(pred) looks like:

 str(pred)
Classes ‘Predict’ and 'data.frame': 625 obs. of  5 variables:
 $ rating    : num  43 44.6 46.2 47.8 49.4 ...
 $ complaints: num  45 45 45 45 45 ...
 $ yhat      : num  39.9 39.5 39.1 38.7 38.3 ...
 $ lower     : num  28 28.3 27.3 25 22 ...
 $ upper     : num  51.7 50.6 50.9 52.4 54.6 ...
snipped

library(rgl)
open3d()
with(pred, surface3d(unique(rating),unique(complaints),yhat,alpha=.2))
with(att, points3d(rating,complaints,privileges, col="red"))
title3d(xlab="rating",ylab="complaints",zlab="privileges")
axes3d()
aspect3d(x=1,z=.05)

Good illustration of the dangers of extrapolation once you realize there are no data out on the extremes of inappropriate extrapolations from that model. The rms-package has a perimeter function and the plotting functions have a perim argument to which perimeter-objects are passed.

Torus