Add points and colored surface to locfit perspective plot

Question

I have a perspective plot of a locfit model and I wish to add two things to it

1. Predictor variables as points in the 3D space
2. Color the surface according to the Z axis value

For the first, I have tried to use the trans3d function. But I get the following error even though my variables are in vector format:

Error in cbind(x, y, z, 1) %*% pmat : requires numeric/complex matrix/vector arguments

Here is a snippet of my code

library(locfit)    
    X <- as.matrix(loc1[,1:2])
    Y <- as.matrix(loc1[,3])  
    zz <- locfit(Y~X,kern="bisq")
    pmat <- plot(zz,type="persp",zlab="Amount",xlab="",ylab="",main="Plains",
         phi = 30, theta = 30, ticktype="detailed")

    x1 <- as.vector(X[,1])
    x2 <- as.vector(X[,2])
    Y <- as.vector(Y)
    points(trans3d(x1,x2,Y,pmat))

My "loc1" data can be found here - https://www.dropbox.com/s/0kdpd5hxsywnvu2/loc1_amountfreq.txt?dl=0

Answer 1

TL,DR: not really in plot.locfit, but you can reconstruct it.

I don't think plot.locfit has good support for this sort of customisation. Supposedly get.data=T in your plot call will plot the original data points (point 1), and it does seem to do so, except if type="persp". So no luck there. Alternatively you can points(trans3d(...)) as you have done, except you need the perspective matrix returned by persp, and plot.locfit.3d does not return it. So again, no luck.

For colouring, typically you make a colour scale (http://r.789695.n4.nabble.com/colour-by-z-value-persp-in-raster-package-td4428254.html) and assign each z facet the colour that goes with it. However, you need the z-values of the surface (not the z-values of your original data) for this, and plot.locfit does not appear to return this either.

So to do what you want, you'll essentially be recoding plot.locfit yourself (not hard, though just cludgy).

You could put this into a function so you can reuse it. We:

1. make a uniform grid of x-y points
2. calculate the value of the fit at each point
3. use these to draw the surface (with a colour scale), saving the perspective matrix so that we can
4. plot your original data

so:

# make a grid of x and y coords, calculate the fit at those points
n.x <- 20 # number of x points in the x-y grid
n.y <- 30 # number of y points in the x-y grid
zz <- locfit(Total ~ Mex_Freq + Cal_Freq, data=loc1, kern="bisq")
xs <- with(loc1, seq(min(Mex_Freq), max(Mex_Freq), length.out=20))
ys <- with(loc1, seq(min(Cal_Freq), max(Cal_Freq), length.out=30))
xys <- expand.grid(Mex_Freq=xs, Cal_Freq=ys)
zs <- matrix(predict(zz, xys), nrow=length(xs))

# generate a colour scale
n.cols <- 100 # number of colours
palette <- colorRampPalette(c('blue', 'green'))(n.cols) # from blue to green
# palette <- colorRampPalette(c(rgb(0,0,1,.8), rgb(0,1,0,.8)), alpha=T)(n.cols) # if you want transparency for example
# work out which colour each z-value should be in by splitting it
#  up into n.cols bins
facetcol <- cut(zs, n.cols)

# draw surface, with colours (col=...)
pmat <- persp(x=xs, y=ys, zs, theta=30, phi=30, ticktype='detailed', main="plains", xlab="", ylab="", zlab="Amount", col=palette[facetcol])

# draw your original data
with(loc1, points(trans3d(Mex_Freq,Cal_Freq,Total,pmat), pch=20))

Note - doesn't look that pretty! might want to adjust say your colour scale colours, or the transparency of the facets, etc. Re: adding legend, there are some other questions that deal with that.

(PS: what a shame ggplot doesn't do 3D scatter plots.)