Bagaimana cara mendapatkan nilai yang digunakan dalam plot.gam dalam mgcv?

10

Saya ingin mengetahui nilai-nilai yang (x, y)digunakan dalam merencanakan plot(b, seWithMean=TRUE)di mgcv paket. Adakah yang tahu bagaimana saya bisa mengekstraksi atau menghitung nilai-nilai ini?

Berikut ini sebuah contoh: 

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)
r time-series smoothing mgcv gung - Pasang kembali Monica
sumber
Saya tidak terbiasa dengan gammodel, tetapi apakah Anda sudah memeriksa atribut yang berbeda dari objek itu? Anda dapat melihat nama-nama objek dengan names(b). Saya menduga detail apa pun yang Anda cari akan disimpan dalam objek itu di suatu tempat.
Chase

Jawaban:

19

Dimulai dengan mgcv1.8-6, secara plot.gamtidak kasat mata mengembalikan data yang digunakannya untuk menghasilkan plot, yaitu melakukan

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

memberi Anda daftar dengan data plot di pd.

JAWABAN DI BAWAH UNTUK mgcv<= 1.8-5:

Saya telah berulang kali mengutuk fakta bahwa plot berfungsi untuk mgcvtidak mengembalikan hal-hal yang mereka rencanakan - yang berikut jelek tapi berhasil:

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)  
}
orang fabian
sumber
Terima kasih banyak atas bantuanmu. Ketika saya mereproduksi kode Anda hingga plotData <<- c(plotData, pd[[i]])})) , pesan berikut terjadi Error in fBody[[i]] : no such index at level 3. Adakah ide mengapa itu tidak berhasil?
Trik "jejak" dulu berhasil bagi saya. Namun, baru-baru ini saya gagal. Saya menduga itu ada hubungannya dengan versi baru dari paket mgcv (saya saat ini menggunakan v 1.8-3), yang mungkin memerlukan argumen "at" yang berbeda dalam fungsi penelusuran. Adakah yang bisa membantu saya tentang cara mendapatkan vektor yang benar untuk argumen "at" fungsi trace? Banyak terima kasih sebelumnya!
@ Paten lihat hasil edit saya.
Fabian
4

Paket visregdapat membuat plot efek yang mirip dengan GAM (tetapi mungkin tidak identik?) Dan apakah memberikan komponen plot sebagai output juga, diformat sebagai daftar. Menggunakan plyr, seseorang dapat membuat dataframe dari output. Contoh:

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))
pengguna13380
sumber
3

Ini bukan jawaban yang lengkap. Semua merencanakan gamobjek dilakukan dengan fungsi plot.gam. Anda dapat melihat kodenya hanya dengan mengetik

> plot.gam

di konsol R. Seperti yang akan Anda lihat, kode ini sangat besar. Apa yang saya dapatkan darinya, bahwa semua plot dilakukan dengan mengumpulkan informasi yang relevan di pdobjek yang merupakan daftar. Jadi salah satu solusi yang mungkin adalah mengedit plot.gam, menggunakan editmisalnya, sehingga mengembalikan objek itu. Menambahkan pdsebelum terakhir }sudah cukup. Saya akan menyarankan menambahkan invisible(pd), sehingga objek ini dikembalikan hanya jika Anda memintanya:

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

Kemudian periksa objek ini dan cari kode plot.gamuntuk baris dengan plotdan lines. Kemudian Anda akan melihat mana yang relevan xdan ynilai muncul di plot.

mpiktas
sumber
oops, saya tidak melihat jawaban Anda ketika saya mengirimkan jawaban saya. Yah, ini sedikit lebih detail sih ....
Fabians
@fabians, jangan khawatir, saya tidak akan memposting milik saya jika saya melihat milikmu. Saya menguraikan ide umum, Anda telah memberikan kode. Karena pertanyaannya meminta kode, jawaban Anda lebih baik.
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
sumber