R：在rgl中的三维散点图中将拟合平面添加到拟合的平面

Question

我正在通过rgl中的点云拟合一个半透明平面。 我得到尽可能

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) 

事情我想提高或增加有：

• 球和拟合平面
• 之间的虚线droplines添加一个黑色网状到合适的飞机
• 也根据Z值颜色代码拟合平面
• 添加Z轴标签“血压”（但旁边的Z轴旋转）

是否有人以任何机会我怎么能做到这一点知道吗？

编辑：响应于所述代码下面的答案我现在用下面的代码来显示（普通）的线性模型与实际数据点一起配合，使用Z轴颜色编码，和具有类似于plotPlane语法在包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

所以这四个问题。从`库（ “RGL”）

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)