1

我正在寻找 R 中的一个函数，它将双变量（一般）线性模型拟合 ( fit=lm(z~x+y)) 显示为 alevelplot或contourplot，其中拟合值以及实际数据点是彩色编码的。

我要寻找的最终结果是这样的（在 Mathematica 中制作了这个，但我现在正在寻找 R 解决方案）

在此处输入图像描述

有谁知道某处是否已经有一个功能可以做这样的事情？

编辑：与此同时，我找到了一个解决方案 - 见下文！

4

3 回答 3

1

也许不是回答自己问题的好形式，但只是发现了函数image.lm()和contour.lm()包rsm，它将线性模型的拟合绘制为图像或等高线图，这正是我正在寻找的。在此基础上，我制作了以下函数，该函数使用类似于plotPlanein package 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)

在此处输入图像描述

于 2014-10-10T09:50:01.110 回答
0

对于不规则间隔的点，我写了 colPoints。如果有人知道一个好的最终包装，请告诉我！

install.packages("berryFunctions")
library(berryFunctions)
?colPoints

i <- c( 22,  40,  80,  45,  60,  63,  30,  70,  55,  48,  32,  48,  70,  40)
j <- c(  5,  33,  12,  56,  20,  40,  45,  45,  30,  36,  23,  15,  30,  10)
k <- c(175, 174, 120, 105, 132, 130, 190, 110, 131, 160, 183, 163, 117, 168)

mod <- lm(k~i+j)
modcoord <- expand.grid(i=seq(20,80, 0.1), j=seq(0,60,0.1))
modvals <- predict(mod, newdata=modcoord)

colPoints(modcoord$i, modcoord$j, modvals, pch=15, add=FALSE)
colPoints(i,j,k, cex=1.5)
points(i,j, cex=1.5)

在此处输入图像描述

于 2014-10-09T21:56:05.693 回答
0

如果将 z 的模型拟合为 x 和 y 的多项式函数，看起来会更好：

 yourdata <- data.frame(
          "x" = c(22, 40, 80, 45, 60,  63, 30, 70, 55, 48, 32, 48, 70, 40),
          "y" = c(5, 33, 12, 56, 20, 40, 45, 45, 30, 36, 23, 15, 30, 10),
          "z" = c(175, 174, 120, 105, 132, 130, 190, 110, 131, 160, 183, 163, 117, 168) )

fit <- lm(z ~ poly(x, y, degree = 2), data = yourdata)

然后绘制：

library(rsm)
image(fit, y ~ x)
contour(fit, y ~ x)
persp(fit, y ~ x, zlab = "z")

结果将显示曲面而不是平面。

使用 ggplot：

# ggplots
library(rsm)
SurfMod <- contour(fit, y ~ x)

# extract matrix values from rsm contour
Xvals <- SurfMod$`x ~ y`[1]
Yvals <- SurfMod$`x ~ y`[2]
Zvals <- SurfMod$`x ~ y`[3]

# form matrix with col and row names 
SurfMatrix <- Zvals$z
colnames(SurfMatrix) <- Yvals$y
rownames(SurfMatrix) <- Xvals$x

# Convert matrix to data frame
library(reshape2)
SurfDF <- melt(SurfMatrix)

library(ggplot2)
library(directlabels)
gg <- ggplot(data = SurfDF) +
  geom_tile(data = SurfDF, aes(Var1, Var2,z = value, fill = value)) +
  stat_contour(data = SurfDF, aes(Var1, Var2, z = value, color = ..level..)) +
  scale_colour_gradient(low = "brown", high = "red") +
  geom_point(data = testdata, aes(x, y, z = z, color = z)) +
  geom_text(data = testdata, aes(x, y,label=z),hjust=0, vjust=0) +
  xlab("x") +
  ylab("y")
direct.label.ggplot(gg, "angled.endpoints")

有关更多 direct.label 绘图方法，请访问http://directlabels.r-forge.r-project.org/docs/index.html

于 2016-01-14T01:51:35.460 回答