Wie erhalte ich die in plot.gam verwendeten Werte in mgcv?

9

Ich möchte die Werte herausfinden, die (x, y)beim Plotten plot(b, seWithMean=TRUE)im mgcv- Paket verwendet werden. Weiß jemand, wie ich diese Werte extrahieren oder berechnen kann?

Hier ist ein Beispiel:

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)
gung - Monica wieder einsetzen
quelle
Ich bin mit gamModellen nicht vertraut , aber haben Sie die verschiedenen Attribute dieses Objekts untersucht? Sie können die Namen der Objekte mit anzeigen names(b). Ich vermute, dass alle Details, nach denen Sie suchen, irgendwo in diesem Objekt gespeichert werden.
Chase

Antworten:

19

Ab mgcv1.8-6 werden plot.gamdie Daten, die zum Generieren der Diagramme verwendet werden , unsichtbar zurückgegeben

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

gibt Ihnen eine Liste mit den Plotdaten in pd.


mgcvANTWORT UNTEN FÜR <= 1,8-5:

Ich habe wiederholt die Tatsache verflucht, dass die Plotfunktionen für das zurückgegebene mgcvMaterial nicht zurückgeben - was folgt, ist hässlich, aber es funktioniert:

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)  
}
Fabianer
quelle
Vielen Dank für deine Hilfe. Wenn ich Ihren Code bis zu reproduziere plotData <<- c(plotData, pd[[i]])})) , wird die folgende Meldung angezeigt Error in fBody[[i]] : no such index at level 3. Irgendwelche Ideen, warum es nicht funktioniert?
Der "Trace" -Trick hat bei mir funktioniert. Allerdings hat es mich kürzlich gescheitert. Ich vermute, es hat mit einer neuen Version des mgcv-Pakets zu tun (ich verwende derzeit Version 1.8-3), für die möglicherweise ein anderes "at" -Argument in der Trace-Funktion erforderlich ist. Kann mir jemand helfen, wie man den richtigen Vektor für das "at" -Argument der Trace-Funktion erhält? Vielen Dank im Voraus!
@Pepijn siehe meine Bearbeitung.
Fabians
4

Das Paket visregkann Effektdiagramme ähnlich wie GAM erstellen (aber möglicherweise nicht identisch?) Und gibt die Diagrammkomponenten auch als Ausgabe an, formatiert als Liste. Mit plyr kann man einen Datenrahmen der Ausgabe erstellen. Beispiel:

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))
user13380
quelle
3

Dies wird keine vollständige Antwort sein. Das gesamte Plotten für gamObjekte erfolgt mit Funktion plot.gam. Sie können den Code durch einfaches Eingeben anzeigen

> plot.gam

in der R-Konsole. Wie Sie sehen werden, ist der Code riesig. Was ich daraus gelernt habe, ist, dass das gesamte Zeichnen durch Sammeln relevanter Informationen in einem pdObjekt erfolgt, das eine Liste ist. Eine der möglichen Lösungen wäre also, beispielsweise zu bearbeiten plot.gam, editdass das Objekt zurückgegeben wird. pdVorletztes Hinzufügen }reicht aus. Ich würde empfehlen, hinzuzufügen invisible(pd), damit dieses Objekt nur zurückgegeben wird, wenn Sie danach fragen:

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

Untersuchen Sie dann dieses Objekt und suchen Sie im Code von plot.gamnach den Zeilen mit plotund lines. Dann werden Sie die von der entsprechenden sehen xund yWerte erscheinen in der Handlung.

mpiktas
quelle
Hoppla, ich habe deine nicht gesehen, als ich meine Antwort gepostet habe. Nun, es ist sowieso ein bisschen detaillierter ...
Fabians
@fabians, keine Sorge, ich hätte meine nicht gepostet, wenn ich deine gesehen hätte. Ich skizzierte allgemeine Idee, Sie haben den Code zur Verfügung gestellt. Da die Frage nach Code fragt, ist Ihre Antwort besser.
mpiktas
0
## 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)
Lauie
quelle