是否有可能绘制与ggplot2适合gamp平滑组件？

本文介绍了是否有可能绘制与ggplot2适合gamp平滑组件？的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！ 问题描述 我使用 mgcv 包中的 gam 拟合模型，并将结果存储在 model ，到目前为止，我一直在使用 plot（model）来查看流畅的组件。我最近开始使用ggplot2并且喜欢它的输出。所以我想知道，是否可以使用ggplot2绘制这些图？ 以下是一个示例： $ b $ = $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ b 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 code>和 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 （数据[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） 名称（范围）= c（平滑$ term，平滑$ by） mat = PredictMat（平滑，范围） par =平滑$ first.para：平滑$ last .para y = mat％*％model $系数[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 ，值= y ，se = se ））}））} 这将返回一个带有光滑组件的融化数据框，所以现在可以在上面的例子中使用 ggplot ： smooths = EvaluateSmooths（模型） 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） 如果有人知道在一般情况下允许这样做的包，我会很感激。解决方案您可以使用与plyr软件包结合的visreg软件包。 （$） library（mgcv） library（visreg） 库（plyr）库（ggplot2） ＃估算gam模型： x1 = rnorm（1000） x2 = rnorm（1000）$ （x1，k = 10）+ s（x2，k = 20），family =poisson ） ＃使用plot = FALSE从visreg获得绘图数据而不绘制 plotdata< - visreg（model，type =contrast，plot = FALSE） ＃visreg的输出是一个长度与'x'变量数量相同的列表，＃因此我们使用ldply从每个列表部分选择我们想要的对象并创建一个数据框： smooths< - ldply（plot data，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）） ＃ggplot： ggplot（smooths，aes（x，smooth））+ geom_line（）+ geom_line（aes（y = lower），linetype =dashed）+ geom_line（aes（y = upper） =dashed）+ facet_grid（。 〜变量，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（ （模型，类型=类型，绘图= FALSE） smooths ） data.frame（变量=部分$ meta $ x，x = part $ fit [[part $ meta $ x]]， smooth = part $ fit $ visregFit， lower = part $适合$ visregLwr， upper = part $ fit $ visregUpr））残差 data.frame（Variable = part $ meta $ x，x = part $ res [[part $ meta $ x]]，y = part $ res $ visregRes）） if（res） ggplot（smooths，aes（x，smo oth））+ 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（。 〜变量，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（模型） ggplot。 model（model，res = TRUE） 颜色可从 http://colorbrewer2.org/ 中挑选。 I am fitting a model using gam from the mgcv package and store the result in model and so far I have been looking at the smooth components using plot(model). I have recently started using ggplot2 and like its output. So I am wondering, is it possible to plot these graphs using ggplot2?Here is an example: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)And I am interest in s(x1, k=10) and s(x2, k=20) not in the fit.Partial answer:I dug deeper into plot.gam and mgcv:::plot.mgcv.smooth and built my own function which extracts the predicted effects and standard errors from the smooth components. It doesn't handle all options and cases of plot.gam so I only consider it a partial solution, but it works well for me.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 )) }))}This returns a "molten" data frame with the smooth components, so it is now possible to use ggplot with the example above :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)If anyone knows a package which allows this in the general case I would be very grateful. 解决方案 You can use the visreg package combined with the plyr package. visreg basically plots any model that you can use predict() on.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 plottingplotdata <- 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")We can put the whole thing into a function, and add an option to show the residuals from the model (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)Colors are picked from http://colorbrewer2.org/. 这篇关于是否有可能绘制与ggplot2适合gamp平滑组件？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！