2014-10-09 1 views
1

Я ищу функцию в R, которая будет отображать двумерную (общую) линейную модель (fit=lm(z~x+y)) как levelplot или contourplot, в соответствии с которой установленные значения, а также фактические точки данных кодируются цветом.R: дисплей (g) lm подходит как цветное изображение/контурный график

Конечный результат я бы искал бы что-то вроде этого (сделал это один в Mathematica, но я искал R решение сейчас)

enter image description here

Кто-нибудь есть идеи, ли там уже есть какая-то функция, которая делает что-то вроде этого?

EDIT: тем временем я нашел решение - см. Ниже!

ответ

1

Может быть, не хорошая форма ответить собственный вопрос, но только обнаружили функция 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) 

enter image description here

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) 

enter image description here

+0

Ха, это не плохо - большое спасибо за это! Одна вещь, которую можно было бы улучшить, - добавить тонкую черную линию вокруг точек (поскольку в противном случае подходящие точки не будут отображаться) и, возможно, использовать градиентную заливку для отдельных маленьких прямоугольников, возможно, используя http: // svitsrv25.epfl.ch/R-doc/library/plotrix/html/gradient.rect.html, чтобы сделать его более гладким. Как вы думаете? –

+0

Ха и я также замечаем, что ваша функция работает только в том случае, если вы используете соотношение сторон 1: 1 - в моем окне сюжета R Studio я получаю белые вертикальные полосы, показывающие ... –

+0

Каким будет способ создания сетки/сетки на самом деле? Я пробовал с breaks = 1000, но это не похоже на это. Есть предположения? –

0

бы выглядеть лучше, если вы подходите модель г как полиномиальная функция от х и у:

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

Смежные вопросы