R: display (g)lm fit as colour coded image / contour plot

Question

I am on the lookout for a function in R that would display a bivariate (general) linear model fit (fit=lm(z~x+y)) as a levelplot or contourplot, whereby the fitted values as well as the actual data points are colour coded.

The final result I would be looking for would be something like this (made this one in Mathematica, but I am looking for an R solution now)

Does anybody have any idea whether there is already a function out there somewhere that does something like this?

EDIT: in the meantime I found a solution - see below!

Answer 1

Maybe not good form to answer one's own question, but just discovered function image.lm() and contour.lm() in package rsm which plots the fit of a linear model as an image or contour plot, which is what I was looking for. Building on that, I made the following function which plots the result of a (general) linear model fit as a function of two explanatory variables as an image or contour plot, together with the actual data points, using a syntax similar to plotPlane in package rockchalk (the mean value is used for any variables not in the model):

plotImage=function(model=NULL,plotx=NULL,ploty=NULL,plotPoints=T,plotContours=T,plotLegend=F,npp=1000,xlab=NULL,ylab=NULL,zlab=NULL,xlim=NULL,ylim=NULL,pch=16,cex=1.2,lwd=0.1,col.palette=NULL) {
      library(rockchalk)
      library(aqfig)
      library(colorRamps)
      mf=model.frame(model);emf=rockchalk::model.data(model)
      if (is.null(xlab)) xlab=plotx
      if (is.null(ylab)) ylab=ploty
      if (is.null(zlab)) zlab=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))
      x=emf[,plotx];y=emf[,ploty];z=mf[,1]
      if (is.null(xlim)) xlim=c(min(x)*0.95,max(x)*1.05)
      if (is.null(ylim)) ylim=c(min(y)*0.95,max(y)*1.05)
      preds=predictOMatic(model,predVals=c(plotx,ploty),n=npp,divider="seq")
      zpred=matrix(preds[,"fit"],npp,npp)
      zlim=c(min(c(preds$fit,z)),max(c(preds$fit,z)))
      par(mai=c(1.2,1.2,0.5,1.2),fin=c(6.5,6))
      graphics::image(x=seq(xlim[1],xlim[2],len=npp),y=seq(ylim[1],ylim[2],len=npp),z=zpred,xlab=xlab,ylab=ylab,col=col.palette,useRaster=T,xaxs="i",yaxs="i")
      if (plotContours) graphics::contour(x=seq(xlim[1],xlim[2],len=npp),y=seq(ylim[1],ylim[2],len=npp),z=zpred,xlab=xlab,ylab=ylab,add=T,method="edge")
      if (plotPoints) {cols1=col.palette[(z-zlim[1])*999/diff(zlim)+1]
                       pch1=rep(pch,length(n))
                       cols2=adjustcolor(cols1,offset=c(-0.3,-0.3,-0.3,1))
                       pch2=pch-15
                       points(c(rbind(x,x)),c(rbind(y,y)), cex=cex,col=c(rbind(cols1,cols2)),pch=c(rbind(pch1,pch2)),lwd=lwd) }
      box()
      if (plotLegend) vertical.image.legend(zlim=zlim,col=col.palette) # TO DO: add z axis label, maybe make legend a bit smaller?
    }

# 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)
plotImage(fit1,plotx="age",ploty="bmi",plotContours=F,plotLegend=T)