Skip to content

plotting gams in mgcv

2 messages · john benson, Ulf Köther

Original
john benson
# 
Hi,
I am running a GAM that has a smooth x factor interaction accomplished with the "by" function.
The model looks like this where X1 is a continuous smooth term and X2 is a categorical factor variable with 2 levels:
Gam1<-gam(y~s(X1, by = X2) + X2, data=datum)summary(Gam1)plot (Gam1)
The plot of this model gives separate figures for the separate smooth functions for each level (n=2) of the factor term.  I would like the rug plot for each figure to show only the locations of the covariate values along the x-axis that are specific to each level of the factor term.  
Instead, the gam plot gives me identical rug plots (with all covariate locations across both factor levels) for both of the 2 plots of the separate smooth functions.
Does anyone know if it is possible to only display the rug plot for covariate values associated with correct level of the factor variable for each of these 2 figures? 
Many thanks for any help with this!
Best,
John
Ulf Köther
# 
Dear John,

because no one answered your question yet, I give it a try. I am no
expert, so in case someone knows an answer via just tweaking the
plot.gam-function, please come forward..!

Here is a solution making use of the predict.gam-function and ggplot2:

# Generate data, example from help files of mgcv 1.8-4, function "s":
library(mgcv)

set.seed(0)
n <- 200
sig2 <- 4
x1 <- runif(n, 0, 1)
x2 <- runif(n, 0, 1)
x3 <- runif(n, 0, 1)
fac <- c(rep(1, n/2), rep(2, n/2)) # create factor
fac.1 <- rep(0, n) + (fac == 1)
fac.2 <- 1 - fac.1 # and dummy variables
fac <- as.factor(fac)
f1 <- exp(2 * x1) - 3.75887
f2 <- 0.2 * x1^11 * (10 * (1 - x1))^6 + 10 * (10 * x1)^3 * (1 - x1)^10
f <- f1 * fac.1 + f2 * fac.2 + x2
e <- rnorm(n, 0, sqrt(abs(sig2)))
y <- f + e

# Putting everything in a dataframe:

data <- data.frame(cbind(fac, x1, x2, y))
colnames(data) <- c("fac", "x1", "x2", "y")
data$fac <- factor(data$fac)

# Housekeeping:
rm(e, f, f1, f2, fac, fac.1, fac.2, n, sig2, x1, x2, x3, y)

# Model fitting:
b <- gam(y ~ fac + s(x1, by = fac) + x2, data = data) 
summary(b)
plot(b, pages = 1, shade = TRUE)

# Plotting with predict and ggplot2:
library(ggplot2)

pred <- expand.grid(fac = factor(c(1,2)),
                    x2 = mean(data$x2),
                    x1 = seq(from = min(data$x1), 
                             to = max(data$x1),
                             length = 100))
# Ordering works better in regard of the split of the predictions later on
pred <- pred[order(pred$fac),] 

# using the predict.gam-function with type = "terms" is essential:
predNew <- predict(b, newdata = pred, type = "terms", se.fit = TRUE)

# Notice the manual split of the fit and the se to get both factor levels
# in one variable:
predPlot <- data.frame(cbind(pred.fit = c(predNew$fit[1:100, "s(x1):fac1"],
                                          predNew$fit[101:200, "s(x1):fac2"]),
                             pred.se = c(predNew$se.fit[1:100, "s(x1):fac1"],
                                         predNew$se.fit[101:200, "s(x1):fac2"]),
                             pred))

# CIs:
predPlot$CI_Low <- predPlot$pred.fit - qnorm(0.975) * predPlot$pred.se
predPlot$CI_High <- predPlot$pred.fit + qnorm(0.975) * predPlot$pred.se

# Plotting with a coloured rug to highlight the difference:
ggplot(predPlot, aes(y = pred.fit, x = x1)) + geom_line() + 
    facet_grid(. ~ fac) + theme_bw() +
    geom_ribbon(aes(x = x1, ymin = CI_Low, ymax = CI_High), alpha = 0.3) +
    geom_rug(data = data, aes(y = y, x = x1, colour = fac),
             size = 0.6, alpha = 0.5, sides = "b") + 
    theme(legend.position = "none")


Good luck!


Am 25.01.2015 um 22:36 schrieb john benson:
--

_____________________________________________________________________

Universit?tsklinikum Hamburg-Eppendorf; K?rperschaft des ?ffentlichen Rechts; Gerichtsstand: Hamburg | www.uke.de
Vorstandsmitglieder: Prof. Dr. Burkhard G?ke (Vorsitzender), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Pr?l?, Rainer Schoppik
_____________________________________________________________________

SAVE PAPER - THINK BEFORE PRINTING