Comment obtenir les valeurs utilisées dans plot.gam en mgcv?

10

Je voudrais connaître les valeurs (x, y)utilisées dans le traçage plot(b, seWithMean=TRUE)dans le package mgcv . Est-ce que quelqu'un sait comment extraire ou calculer ces valeurs?

Voici un exemple:

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 - Réintégrer Monica
la source
Je ne connais pas les gammodèles, mais avez-vous examiné les différents attributs de cet objet? Vous pouvez regarder les noms des objets avec names(b). Je suppose que tous les détails que vous recherchez seront conservés quelque part dans cet objet.
Chase

Réponses:

19

À partir de mgcv1.8-6, plot.gamretourne de manière invisible les données qu'il utilise pour générer les tracés, c'est-à-dire faire

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

vous donne une liste avec les données de traçage dans pd.


RÉPONSE CI-DESSOUS POUR mgcv<= 1,8-5:

J'ai maudit à plusieurs reprises le fait que les fonctions de tracé mgcvne renvoient pas les éléments qu'ils tracent - ce qui suit est laid mais cela fonctionne:

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)  
}
fabians
la source
Merci beaucoup pour votre aide. Lorsque je reproduis votre code jusqu'à plotData <<- c(plotData, pd[[i]])})) , le message suivant se produit Error in fBody[[i]] : no such index at level 3. Des idées pourquoi cela ne fonctionne pas?
L'astuce "trace" fonctionnait pour moi. Cependant, il m'a récemment échoué. Je soupçonne que cela a à voir avec une nouvelle version du package mgcv (j'utilise actuellement la v 1.8-3), qui peut nécessiter un argument "at" différent dans la fonction de trace. Quelqu'un peut-il m'aider à obtenir le bon vecteur pour l'argument "at" de la fonction trace? Merci d'avance!
@Pepijn voir ma modification.
fabians
4

Le package visregpeut créer des tracés d'effets similaires à GAM (mais peut-être pas identiques?) Et donne également les composants de tracé en sortie, formatés sous forme de liste. En utilisant plyr on peut créer une trame de données de la sortie. Exemple:

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
la source
3

Ce ne sera pas une réponse complète. Tout le traçage des gamobjets se fait avec la fonction plot.gam. Vous pouvez consulter son code en tapant simplement

> plot.gam

dans la console R. Comme vous le verrez, le code est énorme. Ce que j'en ai glané, c'est que tout le traçage se fait en rassemblant des informations pertinentes dans l' pdobjet qui est une liste. Ainsi, l'une des solutions possibles serait de modifier plot.gam, en utilisant editpar exemple, pour qu'il renvoie cet objet. Ajouter pdavant-dernier }sera suffisant. Je conseillerais d'ajouter invisible(pd), pour que cet objet ne soit retourné que si vous le demandez:

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

Inspectez ensuite cet objet et recherchez dans le code de plot.gamles lignes avec plotet lines. Ensuite , vous verrez que des pertinentes xet des yvaleurs apparaissent dans l'intrigue.

mpiktas
la source
oups, je n'ai pas vu la vôtre quand j'ai posté ma réponse. Eh bien, c'est un peu plus détaillé de toute façon ....
fabians
@fabians, pas de soucis, je n'aurais pas posté le mien si j'avais vu le vôtre. J'ai esquissé l'idée générale, vous avez fourni le code. Puisque la question demande du code, votre réponse est meilleure.
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
la source