是否有可能绘制与ggplot2适合gamp平滑组件？ [英] Is it possible to plot the smooth components of a gam fit with ggplot2?
问题描述
我使用 mgcv 包中的 gam 拟合模型，并将结果存储在 model ，到目前为止，我一直在使用 plot（model）来查看流畅的组件。我最近开始使用ggplot2并且喜欢它的输出。所以我想知道，是否可以使用ggplot2绘制这些图？
以下是一个示例：
<$ p （1000）exp（x1）+ x2 ^ 2）
$ 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）
我对<$ c感兴趣$ c> s（x1，k = 10）和 s（x2，k = 20）不适合。
部分回答： 我深入了解 这将返回一个带有光滑组件的融化数据框，所以现在可以在上面的例子中使用 如果有人知道在一般情况下允许这样做的包，我会很感激。 您可以使用与plyr软件包结合的visreg软件包。 （$） 我们可以把整个东西放到一个函数中， （res = TRUE）： I am fitting a model using Here is an example: And I am interest in Partial answer: I dug deeper into This returns a "molten" data frame with the smooth components, so it is now possible to use 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. We can put the whole thing into a function, and add an option to show the residuals from the model (res = TRUE):
Colors are picked from http://colorbrewer2.org/. 这篇关于是否有可能绘制与ggplot2适合gamp平滑组件？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
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）
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）
ggplot.model< - function（model，type =conditional，res = FALSE，
col.line =＃7fc97f，col.point =＃beaed4，size.line = 1，size.point = 1）{
require（ （模型，类型=类型，绘图= FALSE）
smooths < - ldply（绘图数据，函数（部分）
） data.frame（变量=部分$ meta $ x，
x = part $ fit [[part $ meta $ x]]，
smooth = part $ fit $ visregFit，
lower = part $适合$ visregLwr，
upper = part $ fit $ visregUpr））
残差< - 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，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/ 中挑选。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?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) and s(x2, k=20) not in the fit.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
))
}))
}
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)
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")
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)
