Ist es möglich, die glatten Komponenten eines Gam Fit mit ggplot2 zu plotten?

Question

Ich bin ein Modell mit gam aus dem mgcv Paket und speichern Sie das Ergebnis in model und bis jetzt habe ich auf die glatten Komponenten mit plot(model). Ich habe kürzlich begonnen, ggplot2 zu verwenden und mag seine Ausgabe. Ich frage mich, ob es möglich ist, diese Graphen mit ggplot2 zu plotten? Hier

ein Beispiel:

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) 

Und ich bin Interesse an s(x1, k=10) und s(x2, k=20) nicht in der Passform.

Teil Antwort:

Ich grub tiefer in plot.gam und mgcv:::plot.mgcv.smooth und baute meine eigene Funktion, die die vorhergesagten Effekte und Standardfehler von den glatten Komponenten extrahiert. Es behandelt nicht alle Optionen und Fälle von plot.gam, also halte ich es nur für eine Teillösung, aber es funktioniert gut für mich.

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 
    )) 
    })) 
} 

Dies liefert einen „geschmolzenen“ Datenrahmen mit den glatten Komponenten, so ist es nun möglich, ggplot mit dem Beispiel zu verwenden, oben:

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) 

Wenn jemand ein Paket kennt, welche diese in die ermöglicht Allgemeiner Fall wäre ich sehr dankbar.

Answer 1

Sie können das visreg-Paket in Kombination mit dem plyr-Paket verwenden. visreg zeichnet im Prinzip jedes Modell auf, das Sie mit pred() on verwenden können.

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") 

Wir können das Ganze in eine Funktion setzen, und eine Option hinzufügen, die Residuen aus dem Modell zu zeigen (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) 

Farben von http://colorbrewer2.org/ gepflückt .

Answer 2

FYI, können visreg direkt Ausgang ein gg Objekt:

visreg(model, "x1", gg=TRUE)