ggplot2로 gam fit의 부드러운 구성 요소를 플로팅 할 수 있습니까?

Question

mgcv 패키지에서 gam을 사용하여 모델을 맞추고 결과를 model에 저장하고 지금까지 plot(model)을 사용하여 부드러운 구성 요소를 살펴 봤습니다. 나는 최근에 ggplot2를 사용하고 출력을 좋아하기 시작했다. 그래서 궁금, 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) 및 s(x2, k=20)하지 적합성에 관심을입니다.

부분 대답은 :

나는 plot.gam 및 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(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를 사용하는 것이 가능하다 : 사람이 이것을 할 수있는 패키지를 알고있는 경우

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) 

일반적인 경우 나는 매우 감사 할 것입니다.

Answer 1

plyr 패키지와 결합 된 visreg 패키지를 사용할 수 있습니다. visreg는 기본적으로 predict()를 사용할 수있는 모든 모델을 그립니다.

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

우리는 함수로 전체를 놓고, 모델 (고해상도 = TRUE)에서 잔류를 표시하는 옵션을 추가 할 수 있습니다 색상 http://colorbrewer2.org/에서 선택됩니다

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) 

.

Answer 2

참고로, visreg 직접 출력 gg 개체 수 있습니다 : 나는 geom_smooth을 알고있는 것처럼

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