¿Cómo obtener los valores utilizados en plot.gam en mgcv?

Question

Me gustaría saber los valores (x, y) utilizado en el trazado plot(b, seWithMean=TRUE) en mgcv paquete. ¿Alguien sabe cómo puedo extraer o calcular estos valores?

He aquí un ejemplo:

library(mgcv) 
set.seed(0)
dat <- gamSim(1, n=400, dist="normal", scale=2) 
b   <- gam(y~s(x0), data=dat) 
plot(b, seWithMean=TRUE)

Answer 1

A partir de mgcv 1.8-6, plot.gam devuelve de forma invisible los datos que utiliza para generar los gráficos, es decir, hace

pd <- plot(<some gam() model>)

le da una lista con los datos de trazado en pd .

---

RESPUESTA A CONTINUACIÓN PARA mgcv <= 1.8-5:

He maldecido repetidamente el hecho de que las funciones de trama para mgcv no devuelvan lo que están tramando lo que sigue es feo pero funciona:

library(mgcv) 
set.seed(0)
dat <- gamSim(1, n = 400, dist = "normal", scale = 2)
b <- gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat)

plotData <- list()
trace(mgcv:::plot.gam, at = list(c(27, 1)), 
  ## tested for mgcv_1.8-4. other versions may need different at-argument.
  quote({
    message("ooh, so dirty -- assigning into globalenv()'s plotData...")
    plotData <<- pd
    }))
mgcv::plot.gam(b, seWithMean = TRUE, pages = 1)

par(mfrow = c(2, 2))
for (i in 1:4) {
  plot(plotData[[i]]$x, plotData[[i]]$fit, type = "l", xlim = plotData[[i]]$xlim,
    ylim = range(plotData[[i]]$fit + plotData[[i]]$se, plotData[[i]]$fit -
      plotData[[i]]$se))
  matlines(plotData[[i]]$x, cbind(plotData[[i]]$fit + plotData[[i]]$se, 
    plotData[[i]]$fit - plotData[[i]]$se), lty = 2, col = 1)
  rug(plotData[[i]]$raw)  
}

Answer 2

Paquete visreg puede hacer gráficos de efectos similares a los de GAM (¿pero quizás no idénticos?) y también da los componentes del gráfico como salida, formateados como una lista. Usando plyr se puede hacer un marco de datos de la salida. Ejemplo:

plot <- visreg(model, type = "contrast")
smooths <- ldply(plot, function(part)   
  data.frame(x=part$x$xx, smooth=part$y$fit, lower=part$y$lwr, upper=part$y$upr))

Answer 3

Esta no será una respuesta completa. Todo el trazado para gam objetos se hace con la función plot.gam . Puede consultar su código simplemente escribiendo

> plot.gam

en la consola de R. Como verás el código es enorme. Lo que he deducido de él, es que todo el trazado se hace recogiendo información relevante en pd que es una lista. Así que una de las posibles soluciones sería editar plot.gam , utilizando edit por ejemplo, para que devuelva ese objeto. Añadir pd antes de la última } será suficiente. Yo aconsejaría añadir invisible(pd) para que este objeto sea devuelto sólo si se lo pide:

> pd <- plot(b,seWithMean = TRUE)

A continuación, inspeccione este objeto y busque en el código de plot.gam para las líneas con plot y lines . A continuación, verá cuál de las x y y aparecen en el gráfico.

Answer 4

## And this is the code for multiple variables!
require(mgcv)
n      = 100
N      = n
tt     = 1:n
arfun  = c(rep(.7,round(n/3)),rep(.3,round(n/3)),rep(-.3,ceiling(n/3)))
arfun2 = c(rep(.8,round(n/3)),rep(.3,round(n/3)),rep(-.3,ceiling(n/3)))
int    = .1*(tt-mean(tt))/max(tt)-.1*((tt-mean(tt))/(max(tt)/10))^2
y      = rep(NA,n)
s.sample <- N
x        <- 10*rnorm(s.sample)
z        <- 10*rnorm(s.sample)
for(j in 1:n){
  y[j]=int[j]+x[j]*arfun[j]+z[j]*arfun2[j]+rnorm(1)  
}

mod = gam(y ~ s(tt) + s(tt, by=x) + s(tt, by=z)) 
## getting the data out of the plot
plotData <- list()
trace(mgcv:::plot.gam, at=list(c(25,3,3,3)),
      # this gets you to the location where plot.gam calls 
      #    plot.mgcv.smooth (see ?trace)
      # plot.mgcv.smooth is the function that does the actual plotting and
      # we simply assign its main argument into the global workspace
      # so we can work with it later.....

      quote({
        # browser()
        print(pd)
        plotData <<- c(plotData, pd)
      }))

# test: 
mgcv::plot.gam(mod, seWithMean=TRUE)

# see if it succeeded
slct = 3
plot(plotData[[slct]]$x, plotData[[slct]]$fit, type="l", xlim=plotData$xlim, 
     ylim=range(plotData[[slct]]$fit + plotData[[slct]]$se, plotData[[slct]]$fit - 
                plotData[[slct]]$se))
matlines(plotData[[slct]]$x, 
         cbind(plotData[[slct]]$fit + plotData[[slct]]$se, 
               plotData[[slct]]$fit - plotData[[slct]]$se), lty=2, col=1)
rug(plotData[[slct]]$raw)