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

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

Aquí hay 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)

Comenzando con mgcv1.8-6, plot.gaminvisiblemente devuelve los datos que usa para generar los gráficos, es decir

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:

Repetidamente he maldecido el hecho de que las funciones de trazado mgcvno devuelven 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)  
}

El paquete visregpuede hacer gráficos de efectos similares a GAM (¿pero quizás no idénticos?) Y también proporciona los componentes de gráficos como salida, formateados como una lista. El uso de plyr one 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))

Esta no será una respuesta completa. Todo el trazado de los gamobjetos se realiza con función plot.gam. Puede ver su código simplemente escribiendo

> plot.gam

en la consola R. Como verás, el código es enorme. Lo que deduje de esto es que todo el trazado se realiza mediante la recopilación de información relevante en el pdobjeto, que es una lista. Entonces, una de las posibles soluciones sería editar plot.gam, usando, editpor ejemplo, para que devuelva ese objeto. Agregar pdantes del último }será suficiente. Aconsejaría agregar invisible(pd), para que este objeto se devuelva solo si lo solicita:

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

Luego inspeccione este objeto y busque en el código de plot.gamlas líneas con ploty lines. Luego verá cuáles de los valores relevantes xy yaparecen en la gráfica.

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