是否可以绘制与ggplot2适合的gam平滑组件？

Question

我正在装配包装中的gam，并将结果存储在model中，到目前为止，我一直在使用plot(model)查看光滑组件。我最近开始使用ggplot2并且喜欢它的输出。所以我想知道，是否可以使用ggplot2绘制这些图？

下面是一个例子：

x1 = rnorm(1000) 
x2 = rnorm(1000) 
n = rpois(1000, exp(x1) + x2^2) 

model = gam(n ~ s(x1, k=10) + s(x2, k=20), family="poisson") 
plot(model, rug=FALSE, select=1) 
plot(model, rug=FALSE, select=2) 

我在s(x1, k=10)和s(x2, k=20)不配合我的兴趣。

部分答案：

我更深入地挖掘plot.gam和mgcv:::plot.mgcv.smooth，并建立了自己的功能，提取从平滑部件的预测效果和标准误差。它不处理所有的选项和plot.gam的情况，所以我只认为它是一个部分的解决方案，但它适用于我。

EvaluateSmooths = function(model, select=NULL, x=NULL, n=100) { 
    if (is.null(select)) { 
    select = 1:length(model$smooth) 
    } 
    do.call(rbind, lapply(select, function(i) { 
    smooth = model$smooth[[i]] 
    data = model$model 

    if (is.null(x)) { 
     min = min(data[smooth$term]) 
     max = max(data[smooth$term]) 
     x = seq(min, max, length=n) 
    } 
    if (smooth$by == "NA") { 
     by.level = "NA" 
    } else { 
     by.level = smooth$by.level 
    } 
    range = data.frame(x=x, by=by.level) 
    names(range) = c(smooth$term, smooth$by) 

    mat = PredictMat(smooth, range) 
    par = smooth$first.para:smooth$last.para 

    y = mat %*% model$coefficients[par] 

    se = sqrt(rowSums(
     (mat %*% model$Vp[par, par, drop = FALSE]) * mat 
    )) 

    return(data.frame(
     label=smooth$label 
     , x.var=smooth$term 
     , x.val=x 
     , by.var=smooth$by 
     , by.val=by.level 
     , value = y 
     , se = se 
    )) 
    })) 
} 

这将返回与光滑部件的“熔融”的数据帧，因此现在可以使用ggplot与上面的例子：如果有人知道的软件包，其允许该在

smooths = EvaluateSmooths(model) 

ggplot(smooths, aes(x.val, value)) + 
    geom_line() + 
    geom_line(aes(y=value + 2*se), linetype="dashed") + 
    geom_line(aes(y=value - 2*se), linetype="dashed") + 
    facet_grid(. ~ x.var) 

一般情况下我会很感激。

Answer 1

您可以将visreg软件包与plyr软件包结合使用。 visreg基本上绘制你可以使用预测（）的任何模型。

library(mgcv) 
library(visreg) 
library(plyr) 
library(ggplot2) 

# Estimating gam model: 
x1 = rnorm(1000) 
x2 = rnorm(1000) 
n = rpois(1000, exp(x1) + x2^2) 
model = gam(n ~ s(x1, k=10) + s(x2, k=20), family="poisson") 

# use plot = FALSE to get plot data from visreg without plotting 
plotdata <- visreg(model, type = "contrast", plot = FALSE) 

# The output from visreg is a list of the same length as the number of 'x' variables, 
# so we use ldply to pick the objects we want from the each list part and make a dataframe: 
smooths <- ldply(plotdata, function(part) 
    data.frame(Variable = part$meta$x, 
      x=part$fit[[part$meta$x]], 
      smooth=part$fit$visregFit, 
      lower=part$fit$visregLwr, 
      upper=part$fit$visregUpr)) 

# The ggplot: 
ggplot(smooths, aes(x, smooth)) + geom_line() + 
    geom_line(aes(y=lower), linetype="dashed") + 
    geom_line(aes(y=upper), linetype="dashed") + 
    facet_grid(. ~ Variable, scales = "free_x") 

我们可以把整个事情到一个函数，并添加一个选项，从模型（RES = TRUE）显示残差：

ggplot.model <- function(model, type="conditional", res=FALSE, 
         col.line="#7fc97f", col.point="#beaed4", size.line=1, size.point=1) { 
    require(visreg) 
    require(plyr) 
    plotdata <- visreg(model, type = type, plot = FALSE) 
    smooths <- ldply(plotdata, function(part) 
    data.frame(Variable = part$meta$x, 
      x=part$fit[[part$meta$x]], 
      smooth=part$fit$visregFit, 
      lower=part$fit$visregLwr, 
      upper=part$fit$visregUpr)) 
    residuals <- ldply(plotdata, function(part) 
    data.frame(Variable = part$meta$x, 
       x=part$res[[part$meta$x]], 
       y=part$res$visregRes)) 
    if (res) 
    ggplot(smooths, aes(x, smooth)) + geom_line(col=col.line, size=size.line) + 
     geom_line(aes(y=lower), linetype="dashed", col=col.line, size=size.line) + 
     geom_line(aes(y=upper), linetype="dashed", col=col.line, size=size.line) + 
     geom_point(data = residuals, aes(x, y), col=col.point, size=size.point) + 
     facet_grid(. ~ Variable, scales = "free_x") 
    else 
    ggplot(smooths, aes(x, smooth)) + geom_line(col=col.line, size=size.line) + 
     geom_line(aes(y=lower), linetype="dashed", col=col.line, size=size.line) + 
     geom_line(aes(y=upper), linetype="dashed", col=col.line, size=size.line) + 
     facet_grid(. ~ Variable, scales = "free_x") 
    } 

ggplot.model(model) 
ggplot.model(model, res=TRUE) 

颜色从http://colorbrewer2.org/挑。

Answer 2

FYI，visreg可以直接输出一个gg对象：

visreg(model, "x1", gg=TRUE)