Может быть, не хорошая форма ответить собственный вопрос, но только обнаружили функция image.lm()
и contour.lm()
в пакете rsm
, который вычерчивает подгонку линейной модели как изображения или контур участок, который является то, что Я искал. Исходя из этого, я сделал следующую функцию, которая отображает результат (общей) линейной модели, как функцию двух объясняющих переменных в виде изображения или контура, вместе с фактическими точками данных, используя синтаксис, аналогичный plotPlane
в пакете rockchalk
(среднее значение используется для каких-либо переменных не в модели):
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)
Ха, это не плохо - большое спасибо за это! Одна вещь, которую можно было бы улучшить, - добавить тонкую черную линию вокруг точек (поскольку в противном случае подходящие точки не будут отображаться) и, возможно, использовать градиентную заливку для отдельных маленьких прямоугольников, возможно, используя http: // svitsrv25.epfl.ch/R-doc/library/plotrix/html/gradient.rect.html, чтобы сделать его более гладким. Как вы думаете? –
Ха и я также замечаем, что ваша функция работает только в том случае, если вы используете соотношение сторон 1: 1 - в моем окне сюжета R Studio я получаю белые вертикальные полосы, показывающие ... –
Каким будет способ создания сетки/сетки на самом деле? Я пробовал с breaks = 1000, но это не похоже на это. Есть предположения? –