Created
March 7, 2017 23:56
-
-
Save mbedward/f5b1c14b54d9d8d6b158fc0460618cb3 to your computer and use it in GitHub Desktop.
This function retrieves data for a specified smooth term from a fitted gam object (package mgcv) which is handy when you want to graph the smoother with ggplot.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(ggplot2) | |
library(mgcv) | |
# Toy data set | |
dat <- data.frame( | |
x = seq(0, 2*pi, length.out = 100), | |
y = sin(3*x) / sin(x/2) | |
) | |
# Randomly thin the data | |
dat <- dat[sample(1:nrow(dat), nrow(dat)/3), ] | |
# Fit the model | |
m <- gam(y ~ s(x), data = dat) | |
# Plot the smooth term with plot.gam | |
plot(m, shade = TRUE, residuals = TRUE, cex = 5) | |
# Plot the smooth term with ggplot | |
smooth.dat <- getSmoother(m, 1) | |
ggplot(data = smooth.dat, aes(x = x)) + | |
# smooth term bounds and line | |
geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "#ccccff") + | |
geom_line(aes(y = fit), colour = "blue") + | |
# data points | |
geom_point(data = dat, aes(x = x, y = y)) + | |
theme_bw() |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(mgcv) | |
getSmoother <- function(model, term = 1) { | |
# mgcv::plot.gam invisibly returns the data used to construct the plot | |
dat <- plot(model, select = 0)[[term]] | |
ndim <- sum(c("x", "y") %in% names(dat)) | |
fit <- dat$fit[,1] | |
lwr <- fit - 2 * dat$se | |
upr <- fit + 2 * dat$se | |
if (ndim == 1) { | |
out <- data.frame(x = dat$x, fit, lwr, upr) | |
} | |
else if (ndim == 2) { | |
xy <- expand.grid(dat$x, dat$y) | |
out <- data.frame(x = xy[,1], y = xy[,2], fit, lwr, upr) | |
} | |
else | |
stop("Only 1D and 2D smoothers are supported") | |
out | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment