Saya mencoba menyesuaikan spline untuk GLM menggunakan R. Setelah saya cocok dengan spline, saya ingin dapat mengambil model yang saya hasilkan dan membuat file pemodelan dalam buku kerja Excel.
Sebagai contoh, katakanlah saya memiliki kumpulan data di mana y adalah fungsi acak x dan kemiringan berubah secara tiba-tiba pada titik tertentu (dalam hal ini @ x = 500).
set.seed(1066)
x<- 1:1000
y<- rep(0,1000)
y[1:500]<- pmax(x[1:500]+(runif(500)-.5)*67*500/pmax(x[1:500],100),0.01)
y[501:1000]<-500+x[501:1000]^1.05*(runif(500)-.5)/7.5
df<-as.data.frame(cbind(x,y))
plot(df)
Saya sekarang cocok menggunakan ini
library(splines)
spline1 <- glm(y~ns(x,knots=c(500)),data=df,family=Gamma(link="log"))
dan hasil saya menunjukkan
summary(spline1)
Call:
glm(formula = y ~ ns(x, knots = c(500)), family = Gamma(link = "log"),
data = df)
Deviance Residuals:
Min 1Q Median 3Q Max
-4.0849 -0.1124 -0.0111 0.0988 1.1346
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.17460 0.02994 139.43 <2e-16 ***
ns(x, knots = c(500))1 3.83042 0.06700 57.17 <2e-16 ***
ns(x, knots = c(500))2 0.71388 0.03644 19.59 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for Gamma family taken to be 0.1108924)
Null deviance: 916.12 on 999 degrees of freedom
Residual deviance: 621.29 on 997 degrees of freedom
AIC: 13423
Number of Fisher Scoring iterations: 9
Pada titik ini, saya dapat menggunakan fungsi prediksi dalam r dan mendapatkan jawaban yang bisa diterima. Masalahnya adalah saya ingin menggunakan hasil model untuk membangun buku kerja di Excel.
Pemahaman saya tentang fungsi prediksi adalah bahwa diberi nilai "x" baru, r memasukkan x baru itu ke fungsi spline yang sesuai (baik fungsi untuk nilai di atas 500 atau yang untuk nilai di bawah 500), maka dibutuhkan hasil dan dikalikan dengan koefisien yang sesuai dan sejak saat itu memperlakukannya seperti istilah model lainnya. Bagaimana cara mendapatkan fungsi spline ini?
(Catatan: Saya menyadari bahwa GLM gamma terkait-log mungkin tidak sesuai untuk set data yang disediakan. Saya tidak bertanya tentang bagaimana atau kapan agar sesuai dengan GLM. Saya menyediakan set itu sebagai contoh untuk tujuan reproduktifitas.)
rm(list=ls())
), terutama bukan tanpa peringatan. Seseorang mungkin copy-paste kode Anda ke sesi terbuka R di mana mereka memiliki beberapa variabel yang sudah (tapi tidak disebutx
,y
,df
atauspline1
) dan miss bahwa kode Anda menghapuskan pekerjaan mereka. Apakah agak bodoh bagi mereka untuk melakukan itu? Iya nih. Tetapi masih sopan untuk membiarkan mereka memutuskan kapan harus menghapus variabel mereka sendiri.Jawaban:
Anda bisa merekayasa balik rumus spline tanpa harus masuk ke
R
kode. Cukup untuk mengetahui hal ituSpline adalah fungsi polinomial piecewise.
Polinomial derajat ditentukan oleh nilainya pada titik .d d+ 1
Koefisien polinomial dapat diperoleh melalui regresi linier.
Dengan demikian, Anda hanya perlu membuat titik ditempatkan di antara setiap pasangan simpul yang berurutan (termasuk titik akhir implisit dari rentang data), memprediksi nilai spline, dan mundur prediksi terhadap kekuatan hingga . Akan ada formula terpisah untuk setiap elemen basis spline dalam setiap simpul "tempat sampah" tersebut. Misalnya, dalam contoh di bawah ini ada tiga simpul internal (untuk empat simpul simpul) dan splines kubik ( ) digunakan, menghasilkan polinomial kubik, masing-masing dengan koefisien. Karena kekuatan relatif tinggid+ 1 x xd d= 3 4 × 4 = 16 d+ 1 = 4 x jika terlibat, sangat penting untuk menjaga semua presisi dalam koefisien. Seperti yang Anda bayangkan, rumus lengkap untuk elemen basis spline bisa jadi cukup panjang!
Seperti yang saya sebutkan beberapa waktu lalu , dapat menggunakan output dari satu program sebagai input dari program lain (tanpa intervensi manual, yang dapat menyebabkan kesalahan yang tidak dapat direproduksi) adalah keterampilan komunikasi statistik yang berguna. Pertanyaan ini memberikan contoh yang bagus tentang bagaimana prinsip itu berlaku: alih-alih menyalin koefisien enam belas digit secara manual, kita bisa meretas bersama cara untuk mengubah splines yang dikomputasi menjadi formula yang dapat dipahami Excel. Yang perlu kita lakukan adalah mengekstrak koefisien spline dari seperti yang dijelaskan di atas, minta itu memformatnya kembali menjadi formula seperti Excel, dan menyalin dan menempelkannya ke dalam Excel.64
R
R
Metode ini akan bekerja dengan perangkat lunak statistik apa pun, bahkan perangkat lunak berpemilik yang tidak berdokumen yang kode sumbernya tidak tersedia.
Berikut adalah contoh yang diambil dari pertanyaan, tetapi dimodifikasi untuk memiliki simpul pada tiga titik internal ( ) serta di titik akhir . Plot menunjukkan versi diikuti oleh rendering Excel. Sangat sedikit kustomisasi yang dilakukan di kedua lingkungan (selain dari menentukan warna untuk mencocokkan warna default Excel sekitar).200 , 500 , 800 ( 1 , 1000 )
R
R
(Garis-garis abu-abu vertikal dalam
R
versi menunjukkan di mana simpul internal berada.)Ini
R
kode lengkapnya . Ini adalah hack yang tidak canggih, bergantung sepenuhnya padapaste
fungsi untuk menyelesaikan manipulasi string. (Cara yang lebih baik adalah membuat templat rumus dan mengisinya menggunakan perintah pencocokan string dan substitusi.)Rumus keluaran spline pertama (dari empat yang diproduksi di sini) adalah
R
sumber
ns.formula
.. apakah Anda berpikir dalam R? Serius meskipun metode Anda terlihat sangat berguna tetapi tampaknya ironis harus meretas hack untuk mendapatkan parameter ini. Akan sangat berguna untuk menampilkan tabel ..Anda sudah melakukan yang berikut:
Sekarang saya akan menunjukkan kepada Anda bagaimana memprediksi (respons) untuk x = 12 dalam dua cara berbeda: Pertama menggunakan fungsi prediksi (cara mudah!)
Cara ke-2 didasarkan pada matriks model secara langsung. Catatan yang saya gunakan
exp
karena fungsi tautan yang digunakan adalah log.Perhatikan bahwa di atas saya mengekstrak elemen ke-12, karena itu sesuai dengan x = 12. Jika Anda ingin memprediksi untuk x di luar set pelatihan, maka cukup Anda dapat kembali menggunakan fungsi prediksi. Katakanlah kita ingin menemukan nilai respons yang diprediksi untuk x = 1100
sumber
Anda mungkin lebih mudah menggunakan basis daya terpotong untuk splines regresi kubik, menggunakan
rms
paket R. Setelah Anda cocok dengan model, Anda dapat mengambil representasi aljabar dari fungsi spline yang pas menggunakan fungsiFunction
ataulatex
dirms
.sumber
Function()
tidak benar-benar mengatakan apa yang dilakukannya. Dalam kasus saya (lihat rincian tentang Rpubs rpubs.com/EmilOWK/rms_splines ), saya mendapatkanfunction(x = NA) {-2863.7787+245.72672* x-0.1391794*pmax(x-10.9,0)^3+0.27835881*pmax(x-50.5,0)^3-0.1391794*pmax(x-90.1,0)^3 } <environment: 0x556156e80db8>
The-2863.7787
nilai adalah koefisien pertama dalam model, yang245.72672
kedua, dan koefisien terakhir-873.0223
tidak terlihat di mana saja persamaan. Hal yang sama berlaku untuk output darilatex()
.Function
bekerja denganGlm()
ketika Anda menggunakanrcs
fungsi spline. Outputnya adalah pengubahan ulang spline dalam bentuk paling sederhana dengan menulis seolah-olah pembatasan ekor linier tidak ada (tetapi mereka) sebagaimana dirinci dalam catatan kursus RMS saya .