[R-sig-ME] plotting gams in mgcv

[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:
> 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 		 	   		  
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

--

_____________________________________________________________________

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