R: adding droplines to fitted plane in 3D scatterplot in rgl

Question

I was making a plot with a semi-transparent plane fitted through a cloud of points in rgl. I got as far as

library(rgl)
BLOOD_PRESSURE=c(132,143,153,162,154,168,137,149,159,128,166)
AGE=c(52,59,67,73,64,74,54,61,65,46,72)
WEIGHT=c(78,83,87,95,88,99,85,85,93,75,98)
fit=lm(BLOOD_PRESSURE~AGE+WEIGHT)
npp=10
plot3d(x=AGE, y=WEIGHT, z=BLOOD_PRESSURE, type="s", col=rainbow(length(BLOOD_PRESSURE))[rank(BLOOD_PRESSURE)], radius=1, zlab="")
grd <- expand.grid(AGE=seq(min(AGE),max(AGE),length.out=npp),
                   WEIGHT=seq(min(WEIGHT),max(WEIGHT),length.out=npp) )
grd$pred <-predict(fit, newdata=grd)
persp3d(x=unique(grd[[1]]), y=unique(grd[[2]]), 
        z=matrix(grd[[3]],npp,npp),  color="lightgrey",alpha=0.7, lit=T, back="lines", add=TRUE)

Things I would like to improve or add are:

• dashed droplines between the spheres and the fitted plane
• add a black mesh to the fitted plane
• also colour code the fitted plane according to the Z value
• add a Z axis label "Blood pressure" (but rotated next to the Z axis)

Does anybody know by any chance how I could achieve this?

EDIT: in response to the answer below the code I am now using the following code to display a (general) linear model fit together with the actual data points, using Z axis colour coding, and with a syntax similar to plotPlane in package rockchalk:

plotPlaneFancy=function(model=NULL,plotx1=NULL,plotx2=NULL,plotPoints=T,plotDroplines=T,npp=50,x1lab=NULL,x2lab=NULL,ylab=NULL,x1lim=NULL,x2lim=NULL,cex=0.5,col.palette=NULL,segcol="darkgrey",interval="none",confcol="lightgrey",confalpha=0.4,lit=T,outfile="graph.png",aspect=c(1,1,0.7),zoom=1,userMatrix=matrix(c(0.80,-0.60,0.022,0,0.23,0.34,0.91,0,-0.55,-0.72,0.41,0,0,0,0,1),ncol=4,byrow=T),windowRect=c(0,29,1920,1032)) { # or library(colorRamps);col.palette <- matlab.like(1000)
  library(rockchalk)
  library(rgl)
  library(colorRamps)
  mf=model.frame(model);emf=rockchalk::model.data(model)
  if (is.null(x1lab)) x1lab=plotx1
  if (is.null(x2lab)) x2lab=plotx2
  if (is.null(ylab)) ylab=names(mf)[[1]]
  if (is.null(col.palette)) col.palette=rev(colorRampPalette(rainbow(13,s=0.9,v=0.8),bias=0.6,interpolate ="spline")(1000))
  x1=emf[,plotx1]
  x2=emf[,plotx2]
  y=mf[,1]
  if (is.null(x1lim)) x1lim=c(min(x1),max(x1))
  if (is.null(x2lim)) x2lim=c(min(x2),max(x2))
  preds=predictOMatic(model,predVals=c(plotx1,plotx2),n=npp,divider="seq",interval=interval)
  ylim=c(min(c(preds$fit,y)),max(c(preds$fit,y)))
  open3d(zoom=zoom,userMatrix=userMatrix,windowRect=windowRect)
  if (plotPoints) plot3d(x=x1,y=x2,z=y,type="s",col=col.palette[(y-min(y))*999/diff(range(y))+1],radius=cex,aspect=aspect,xlab=x1lab,ylab=x2lab,zlab=ylab,lit=lit)
  if (!plotPoints) plot3d(x=x1,y=x2,z=y,type="n",col=col.palette[(y-min(y))*999/diff(range(y))+1],radius=cex,aspect=aspect,xlab=x1lab,ylab=x2lab,zlab=ylab)
  if ("lwr" %in% names(preds)) persp3d(x=unique(preds[,plotx1]),y=unique(preds[,plotx2]),z=matrix(preds[,"lwr"],npp,npp),color=confcol, alpha=confalpha, lit=lit, back="lines",add=TRUE)
  ypred=matrix(preds[,"fit"],npp,npp)
  cols=col.palette[(ypred-min(ypred))*999/diff(range(ypred))+1]
  persp3d(x=unique(preds[,plotx1]),y=unique(preds[,plotx2]),z=ypred,color=cols, alpha=0.7, lit=lit, back="lines",add=TRUE)
  if ("upr" %in% names(preds)) persp3d(x=unique(preds[,plotx1]),y=unique(preds[,plotx2]),z=matrix(preds[,"upr"],npp,npp),color=confcol, alpha=confalpha, lit=lit, back="lines",add=TRUE)
  if (plotDroplines) segments3d(x=rep(x1,each=2),y=rep(x2,each=2),z=matrix(t(cbind(y,fitted(model))),nc=1),col=segcol,lty=2)
  if (!is.null(outfile)) rgl.snapshot(outfile, fmt="png", top=TRUE)
}

# simulate some data
n=10000
age=rnorm(n,mean=40,sd=5)
height=rnorm(n,mean=180,sd=7)
weight=-85+0.8*age+0.004*height^2+rnorm(n,mean=0,sd=7)
bmi=weight/((height/100)^2)
sbp=33+1.8*age+2.1*bmi-0.035*age*bmi+rnorm(n,mean=0,sd=5)
mydata=data.frame(cbind(age,height,weight,bmi,sbp))


fit1=lm(sbp~age*bmi,data=mydata)
plotPlaneFancy(fit1, plotx1 = "age", plotx2 = "bmi",cex=0.6)
plotPlaneFancy(fit1, plotx1 = "age", plotx2 = "bmi",cex=0.5,interval="confidence")
plotPlaneFancy(fit1, plotx1 = "age", plotx2 = "bmi",cex=0.5,interval="prediction")

Answer 1

So this is four questions. Here is an answer for two of them (starting right after you finish defining grd):

z <- matrix(grd[[3]],npp,npp)
col.palette <- rainbow(100)
colors <- col.palette[(z-min(z))*99/diff(range(z))+1]
persp3d(x=unique(grd[[1]]), y=unique(grd[[2]]), z=z, 
        color=colors, alpha=0.7, lit=T, back="lines", add=TRUE)

segments3d(x=rep(AGE,each=2),
           y=rep(WEIGHT,each=2),
           z=matrix(t(cbind(BLOOD_PRESSURE,predict(fit))),nc=1),
           col="red", lty=2)