Berbagai cara memodelkan interaksi antara prediktor kontinu dan kategoris dalam GAM

Pertanyaan berikut didasarkan pada diskusi yang ditemukan di halaman ini . Diberikan variabel respons y, variabel penjelas kontinu xdan faktor fac, dimungkinkan untuk mendefinisikan Model Aditif Umum (GAM) dengan interaksi antara xdan facmenggunakan argumen by=. Menurut file bantuan ?gam.models dalam paket R mgcv, ini dapat dilakukan sebagai berikut:

gam1 <- gam(y ~ fac +s(x, by = fac), ...)

@ GavinSimpson di sini menyarankan pendekatan yang berbeda:

gam2 <- gam(y ~ fac +s(x) +s(x, by = fac, m=1), ...)

Saya telah bermain-main dengan model ketiga:

gam3 <- gam(y ~ s(x, by = fac), ...)

Pertanyaan utama saya adalah: apakah beberapa model ini salah, atau mereka berbeda? Dalam kasus terakhir, apa perbedaan mereka? Berdasarkan contoh yang akan saya bahas di bawah ini saya pikir saya bisa memahami beberapa perbedaan mereka, tetapi saya masih kehilangan sesuatu.

Sebagai contoh saya akan menggunakan dataset dengan spektrum warna untuk bunga dari dua spesies tanaman berbeda yang diukur di lokasi yang berbeda.

rm(list=ls())
# install.packages("RCurl")
library(RCurl) # allows accessing data from URL
df <- read.delim(text=getURL("https://raw.githubusercontent.com/marcoplebani85/datasets/master/flower_color_spectra.txt"))
library(mgcv)

Untuk kejelasan, setiap baris pada gambar di atas mewakili spektrum warna rata-rata yang diprediksi untuk setiap lokasi dengan bentuk GAM terpisah density~s(wl)berdasarkan sampel ~ 10 bunga. Daerah abu-abu mewakili 95% CI untuk setiap GAM.

Tujuan akhir saya adalah untuk memodelkan efek (potensial interaktif) dari Taxondan panjang gelombang wlpada reflektansi (disebut densitydalam kode dan dataset) sambil menghitung Localitysebagai efek acak dalam GAM efek campuran. Untuk saat ini saya tidak akan menambahkan bagian efek campuran ke piring saya, yang sudah cukup penuh dengan mencoba memahami cara memodelkan interaksi.

Saya akan mulai dengan yang paling sederhana dari tiga GAM interaktif:

gam.interaction0 <- gam(density ~ s(wl, by = Taxon), data = df) 
# common intercept, different slopes
plot(gam.interaction0, pages=1)

summary(gam.interaction0)

Menghasilkan:

Family: gaussian 
Link function: identity 

Formula:
density ~ s(wl, by = Taxon)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  28.3490     0.1693   167.4   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                      edf Ref.df     F p-value    
s(wl):TaxonSpeciesA 8.938  8.999 884.3  <2e-16 ***
s(wl):TaxonSpeciesB 8.838  8.992 325.5  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.523   Deviance explained = 52.4%
GCV = 284.96  Scale est. = 284.42    n = 9918

Bagian parametrik sama untuk kedua spesies, tetapi spline yang berbeda dipasang untuk masing-masing spesies. Agak membingungkan untuk memiliki bagian parametrik dalam ringkasan GAM, yang non-parametrik. @IsabellaGhement menjelaskan:

Jika Anda melihat plot perkiraan efek halus (smooths) yang sesuai dengan model pertama Anda, Anda akan melihat bahwa mereka berpusat sekitar nol. Jadi, Anda perlu 'menggeser' kelancaran ke atas (jika pencegatan diperkirakan positif) atau turun (jika pencegatan diperkirakan negatif) untuk mendapatkan fungsi kelancaran yang Anda pikir Anda perkirakan. Dengan kata lain, Anda perlu menambahkan perkiraan intersep ke smooth untuk mendapatkan apa yang Anda inginkan. Untuk model pertama Anda, 'shift' diasumsikan sama untuk kedua smooths.

Bergerak:

gam.interaction1 <- gam(density ~ Taxon +s(wl, by = Taxon, m=1), data = df)
plot(gam.interaction1,pages=1)

summary(gam.interaction1)

Memberi:

Family: gaussian 
Link function: identity 

Formula:
density ~ Taxon + s(wl, by = Taxon, m = 1)

Parametric coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    40.3132     0.1482   272.0   <2e-16 ***
TaxonSpeciesB -26.0221     0.2186  -119.1   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                      edf Ref.df    F p-value    
s(wl):TaxonSpeciesA 7.978      8 2390  <2e-16 ***
s(wl):TaxonSpeciesB 7.965      8  879  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.803   Deviance explained = 80.3%
GCV = 117.89  Scale est. = 117.68    n = 9918

Sekarang, setiap spesies juga memiliki perkiraan parametrik sendiri.

Model selanjutnya adalah model yang sulit saya pahami:

gam.interaction2 <- gam(density ~ Taxon + s(wl) + s(wl, by = Taxon,  m=1), data = df)
plot(gam.interaction2, pages=1)

Saya tidak memiliki gagasan yang jelas tentang apa yang digambarkan oleh grafik ini.

summary(gam.interaction2)

Memberi:

Family: gaussian 
Link function: identity 

Formula:
density ~ Taxon + s(wl) + s(wl, by = Taxon, m = 1)

Parametric coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    40.3132     0.1463   275.6   <2e-16 ***
TaxonSpeciesB -26.0221     0.2157  -120.6   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                      edf Ref.df     F p-value    
s(wl)               8.940  8.994 30.06  <2e-16 ***
s(wl):TaxonSpeciesA 8.001  8.000 11.61  <2e-16 ***
s(wl):TaxonSpeciesB 8.001  8.000 19.59  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.808   Deviance explained = 80.8%
GCV = 114.96  Scale est. = 114.65    n = 9918

Bagian parametrik gam.interaction2hampir sama dengan untuk gam.interaction1, tetapi sekarang ada tiga perkiraan untuk istilah yang lancar, yang tidak dapat saya artikan.

Terima kasih sebelumnya kepada siapa saja yang mau meluangkan waktu untuk membantu saya memahami perbedaan dalam ketiga model ini.

gam1dan gam2baik-baik saja; mereka adalah model yang berbeda, walaupun mereka mencoba untuk melakukan hal yang sama, yaitu model smooth khusus-kelompok.

The gam1bentuk

y ~ f + s(x, by = f)

melakukan ini dengan memperkirakan halus yang terpisah untuk setiap tingkat f(dengan asumsi itu fadalah faktor standar), dan memang, parameter kelancaran terpisah diperkirakan untuk setiap kelancaran juga.

The gam2bentuk

y ~ f + s(x) + s(x, by = f, m = 1)

mencapai tujuan yang sama dengan gam1(memodelkan hubungan yang mulus antara xdan yuntuk setiap level f) tetapi melakukannya dengan memperkirakan efek kelancaran global atau rata-rata xpada y( s(x)istilah) ditambah istilah perbedaan halus ( s(x, by = f, m = 1)istilah kedua ). Karena penalti di sini adalah pada turunan pertama ( m = 1) for this difference smoother, it is penalising departure from a flat line, which when added to the global or average smooth term (s (x) `) mencerminkan penyimpangan dari efek global atau rata-rata.

gam3 bentuk

y ~ s(x, by = f)

salah terlepas dari seberapa baik itu cocok dalam situasi tertentu. Alasan saya mengatakan itu salah adalah bahwa setiap kelancaran yang ditentukan oleh s(x, by = f)bagian berpusat tentang nol karena batasan jumlah-ke-nol yang dikenakan untuk pengidentifikasian model. Dengan demikian, tidak ada dalam model yang memperhitungkan rata-rata$Y$di setiap grup yang ditentukan oleh f. Hanya ada rata-rata keseluruhan yang diberikan oleh model intersep. Ini berarti bahwa lebih halus, yang berpusat sekitar nol dan yang memiliki fungsi basis datar dihapus dari perluasan basis x(seperti yang dikacaukan dengan model mencegat) sekarang bertanggung jawab untuk pemodelan kedua perbedaan dalam rata-rata$Y$untuk grup saat ini dan rata-rata keseluruhan (model mencegat), ditambah efek halus xpada$Y$.

Namun tidak satu pun dari model ini yang sesuai untuk data Anda; mengabaikan untuk saat ini distribusi yang salah untuk respons ( densitytidak boleh negatif dan ada masalah heterogenitas yang familyakan diperbaiki atau ditangani oleh non-Gaussian ), Anda belum memperhitungkan pengelompokan menurut bunga ( SampleIDdalam dataset Anda).

Jika tujuan Anda adalah memodelkan Taxonkurva tertentu, maka model formulir akan menjadi titik awal:

m1 <- gam(density ~ Taxon + s(wl, by = Taxon, k = 20) + s(SampleID, bs = 're'),
          data = df, method = 'REML')

di mana saya telah menambahkan efek acak untuk SampleIDdan meningkatkan ukuran ekspansi basis untuk Taxonsmooths tertentu.

Model ini m1, memodelkan pengamatan yang berasal dari wlefek halus tergantung pada spesies mana ( Taxon) pengamatan berasal ( Taxonistilah parametrik hanya menetapkan rata-rata densityuntuk setiap spesies dan diperlukan seperti yang dibahas di atas), ditambah intersep acak. Secara bersama-sama, kurva untuk bunga individu muncul dari versi bergeser dari Taxonkurva tertentu, dengan jumlah pergeseran yang diberikan oleh intersep acak. Model ini mengasumsikan bahwa semua individu memiliki bentuk smooth yang sama seperti yang diberikan oleh smooth untuk Taxonbunga yang berasal dari individu tersebut.

Versi lain dari model ini adalah gam2formulir dari atas tetapi dengan efek acak tambahan

m2 <- gam(density ~ Taxon + s(wl) + s(wl, by = Taxon, m = 1) + s(SampleID, bs = 're'),
          data = df, method = 'REML')

Model ini lebih cocok tetapi saya tidak berpikir itu menyelesaikan masalah sama sekali, lihat di bawah. Satu hal yang saya pikir disarankan adalah bahwa standarnya kberpotensi terlalu rendah untuk Taxonkurva tertentu dalam model ini . Masih ada banyak variasi smooth residual yang tidak kami modelkan jika Anda melihat plot diagnostik.

Model ini kemungkinan besar terlalu membatasi data Anda; beberapa kurva dalam plot smooths individual Anda tampaknya bukan versi bergeser sederhana dari Taxonkurva rata - rata. Model yang lebih kompleks akan memungkinkan untuk smooths individual juga. Model tersebut dapat diperkirakan menggunakan fsatau dasar interaksi faktor-halus . Kami masih menginginkan Taxonkurva tertentu tetapi kami juga ingin memiliki smooth yang terpisah untuk masing-masing SampleID, tetapi tidak seperti smooth by, saya akan menyarankan bahwa pada awalnya Anda ingin semua SampleIDkurva khusus tersebut memiliki kegigihan yang sama. Dalam arti yang sama dengan intersepsi acak yang kami sertakan sebelumnya, thefs basis menambahkan intersep acak, tetapi juga termasuk spline "acak" (Saya menggunakan kutipan menakut-nakuti seperti dalam interpretasi Bayesian tentang GAM, semua model ini hanya variasi pada efek acak).

Model ini cocok untuk data Anda sebagai

m3 <- gam(density ~ Taxon + s(wl, by = Taxon, k = 20) + s(wl, SampleID, bs = 'fs'), 
          data = df, method = 'REML')

Perhatikan bahwa saya memiliki peningkatan di ksini, kalau-kalau kita perlu lebih banyak kegoyahan dalam Taxonsmooths-spesifik. Kami masih membutuhkan Taxonefek parametrik untuk alasan yang dijelaskan di atas.

Model itu membutuhkan waktu lama untuk dipasangkan pada satu inti dengan gam()- bam()kemungkinan besar akan lebih baik dalam menyesuaikan model ini karena ada sejumlah besar efek acak di sini.

Jika kita membandingkan model-model ini dengan AIC versi koreksi parameter pemilihan kelancaran kita melihat betapa jauh lebih baik model yang terakhir ini m3, dibandingkan dengan dua lainnya meskipun menggunakan urutan besarnya lebih banyak derajat kebebasan

> AIC(m1, m2, m3)
          df      AIC
m1  190.7045 67264.24
m2  192.2335 67099.28
m3 1672.7410 31474.80

Jika kita melihat smooths model ini kita mendapatkan ide yang lebih baik tentang bagaimana hal itu sesuai dengan data:

(Catatan ini diproduksi menggunakan draw(m3)menggunakan draw()fungsi dari saya gratia paket. Warna-warna dalam plot kiri bawah tidak relevan dan tidak membantu di sini.)

Masing SampleID- masing kurva yang dipasang dibangun dari intersep atau istilah parametrik TaxonSpeciesBditambah satu dari dua Taxonsmooth-spesifik, tergantung pada milik Taxonmasing-masing SampleID, ditambah itu sendiri SampleID-specifc smooth.

Perhatikan bahwa semua model ini masih salah karena tidak memperhitungkan heterogenitas; model gamma atau Tweedie dengan tautan log akan menjadi pilihan saya untuk melangkah lebih jauh. Sesuatu seperti:

m4 <- gam(density ~ Taxon + s(wl, by = Taxon) + s(wl, SampleID, bs = 'fs'), 
          data = df, method = 'REML', family = tw())

Tapi saya mengalami masalah dengan pemasangan model ini saat ini, yang mungkin mengindikasikan terlalu rumit dengan banyak smooths yang wldisertakan.

Bentuk alternatif adalah dengan menggunakan pendekatan faktor yang dipesan, yang melakukan dekomposisi seperti ANOVA pada smooths:

• Taxon istilah parametrik dipertahankan
• s(wl)adalah smooth yang akan mewakili level referensi
• s(wl, by = Taxon)akan memiliki perbedaan terpisah halus untuk setiap level lainnya. Dalam kasus Anda, Anda hanya akan memiliki satu di antaranya.

Model ini dipasang seperti m3,

df <- transform(df, fTaxon = ordered(Taxon))
m3 <- gam(density ~ fTaxon + s(wl) + s(wl, by = fTaxon) +
            s(wl, SampleID, bs = 'fs'), 
          data = df, method = 'REML')

tetapi interpretasinya berbeda; yang pertama s(wl)akan merujuk TaxonAdan smooth yang tersirat s(wl, by = fTaxon)akan menjadi perbedaan yang halus antara smooth untuk TaxonAdan yang TaxonB.

Ini adalah apa yang Jacolien van Rij menulis di halaman tutorial nya:

Cara mengatur interaksi tergantung pada jenis prediktor pengelompokan:

• dengan faktor termasuk perbedaan intersepsi: Group + s(Time, by=Group)
• dengan faktor pesanan termasuk perbedaan intersep dan referensi halus: Group + s(Time) + s(Time, by=Group)
• dengan prediktor biner termasuk referensi smooth: s(Time) + s(Time, by=IsGroupChildren)

Variabel kategorikal harus ditentukan sebagai faktor, faktor terurut, atau faktor biner dengan fungsi R yang sesuai. Untuk memahami bagaimana menafsirkan output dan apa yang masing-masing model dapat dan tidak bisa memberi tahu kami, lihat halaman tutorial Jacolien van Rij secara langsung. Tutorialnya juga menjelaskan cara menyesuaikan GAM dengan efek campuran. Untuk memahami konsep interaksi dalam konteks GAM, halaman tutorial ini oleh Peter Laurinec juga berguna. Kedua halaman menyediakan banyak informasi lebih lanjut untuk menjalankan GAM dengan benar dalam berbagai skenario.