Menafsirkan hasil spline

20

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.)

Eric
sumber
7
Saya sarankan, jika mungkin, untuk menghindari memasukkan kode yang menghapus semua variabel ( 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 disebut x, y, dfatau spline1) 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.
Glen_b -Reinstate Monica

Jawaban:

25

Anda bisa merekayasa balik rumus spline tanpa harus masuk ke Rkode. Cukup untuk mengetahui hal itu

  • Spline adalah fungsi polinomial piecewise.

  • Polinomial derajat ditentukan oleh nilainya pada titik .dd+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+1xxdd=34×4=16d+1=4xjika 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.64RR

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)RR

Plot R

Plot Excel

(Garis-garis abu-abu vertikal dalam Rversi menunjukkan di mana simpul internal berada.)


Ini Rkode lengkapnya . Ini adalah hack yang tidak canggih, bergantung sepenuhnya pada pastefungsi untuk menyelesaikan manipulasi string. (Cara yang lebih baik adalah membuat templat rumus dan mengisinya menggunakan perintah pencocokan string dan substitusi.)

#
# Create and display a spline basis.
#
x <- 1:1000
n <- ns(x, knots=c(200, 500, 800))

colors <- c("Orange", "Gray", "tomato2", "deepskyblue3")
plot(range(x), range(n), type="n", main="R Version",
     xlab="x", ylab="Spline value")
for (k in attr(n, "knots")) abline(v=k, col="Gray", lty=2)
for (j in 1:ncol(n)) {
  lines(x, n[,j], col=colors[j], lwd=2)
}
#
# Export this basis in Excel-readable format.
#
ns.formula <- function(n, ref="A1") {
  ref.p <- paste("I(", ref, sep="")
  knots <- sort(c(attr(n, "Boundary.knots"), attr(n, "knots")))
  d <- attr(n, "degree")
  f <- sapply(2:length(knots), function(i) {
    s.pre <- paste("IF(AND(", knots[i-1], "<=", ref, ", ", ref, "<", knots[i], "), ", 
                   sep="")
    x <- seq(knots[i-1], knots[i], length.out=d+1)
    y <- predict(n, x)
    apply(y, 2, function(z) {
      s.f <- paste("z ~ x+", paste("I(x", 2:d, sep="^", collapse=")+"), ")", sep="")
      f <- as.formula(s.f)
      b.hat <- coef(lm(f))
      s <- paste(c(b.hat[1], 
            sapply(1:d, function(j) paste(b.hat[j+1], "*", ref, "^", j, sep=""))), 
            collapse=" + ")
      paste(s.pre, s, ", 0)", sep="")
    })
  })
  apply(f, 1, function(s) paste(s, collapse=" + "))
}
ns.formula(n) # Each line of this output is one basis formula: paste into Excel

Rumus keluaran spline pertama (dari empat yang diproduksi di sini) adalah

"IF(AND(1<=A1, A1<200), -1.26037447288906e-08 + 3.78112341937071e-08*A1^1 + -3.78112341940948e-08*A1^2 + 1.26037447313669e-08*A1^3, 0) + IF(AND(200<=A1, A1<500), 0.278894459758071 + -0.00418337927419299*A1^1 + 2.08792741929417e-05*A1^2 + -2.22580643138594e-08*A1^3, 0) + IF(AND(500<=A1, A1<800), -5.28222778473101 + 0.0291833541927414*A1^1 + -4.58541927409268e-05*A1^2 + 2.22309136420529e-08*A1^3, 0) + IF(AND(800<=A1, A1<1000), 12.500000000002 + -0.0375000000000067*A1^1 + 3.75000000000076e-05*A1^2 + -1.25000000000028e-08*A1^3, 0)"

Rxx

Cuplikan Excel

whuber
sumber
2
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 ..
geotheory
Ini mungkin pertanyaan bodoh: tetapi apakah Anda merencanakan 4 splines, atau 4 basis dari satu spline?
Erosennin
@ Erosennin saya tergantung pada apa yang Anda maksud dengan "satu spline." Keempat kurva ini adalah dasar untuk spline yang masing-masing berbentuk kubik dalam empat interval dan terus menerus dibedakan kedua pada tiga titik di mana interval tersebut bertemu, seperti yang dijelaskan oleh tiga poin poin yang memperkenalkan jawaban saya.
whuber
Terima kasih! Saya tidak bermaksud menjadi nitpicking, Hanya terlihat karena ada empat splines (dari jawaban), dan bukan empat kurva yang merupakan basis. Sekali lagi, saya hanya di sini mencoba memahami ...
Erosennin
1
@ Erosennin Tidak masalah. Mungkin ini akan membantu: "spline" adalah kombinasi linear dari keempat kurva ini yang ditentukan oleh proses pemasangan regresi. Cara lain untuk menggambarkannya: spline terdiri dari ruang vektor kurva yang dapat dibuat dengan mengambil kombinasi linear dari empat kurva ini.
whuber
4

Anda sudah melakukan yang berikut:

> rm(list=ls())
> 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))
> library(splines)
> spline1 <- glm(y~ns(x,knots=c(500)),data=df,family=Gamma(link="log"))
> 

Sekarang saya akan menunjukkan kepada Anda bagaimana memprediksi (respons) untuk x = 12 dalam dua cara berbeda: Pertama menggunakan fungsi prediksi (cara mudah!)

> new.dat=data.frame(x=12)
> predict(spline1,new.dat,type="response")
       1 
68.78721 

Cara ke-2 didasarkan pada matriks model secara langsung. Catatan yang saya gunakan expkarena fungsi tautan yang digunakan adalah log.

> m=model.matrix( ~ ns(df$x,knots=c(500))) 
> prd=exp(coefficients(spline1) %*% t(m)) 
> prd[12]
[1] 68.78721

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

> predict(spline1, newdata=data.frame(x=1100),type="response")
       1 
366.3483 
Stat
sumber
Terima kasih atas tanggapan Anda! Tapi, saya masih bingung: /. Saya tidak yakin saya tahu apa yang harus dilakukan dengan matriks ini. Sebagai contoh, jika saya memiliki x = 12, maka perkirakan kata y = 68.78721, tetapi mencari 12 dari matriks itu saya dapatkan 0,016816392. Intersep dan koefisien asli untuk x <500 adalah masing-masing 4.174603 dan 3.830416. exp (4.174603 + 3.8304116 * 0.016816392) <> 68.78721. Plus, bagaimana saya mendapatkan nilai untuk x jika x tidak ada dalam set pelatihan?
Eric
Saya mengubah jawaban saya.
Stat
Saya menambahkan kode untuk case ketika x tidak ada dalam set pelatihan.
Stat
2
Apakah ada cara untuk mendapatkan 366.3483 untuk x = 1100 tanpa menggunakan fungsi prediksi?
Eric
4

Anda mungkin lebih mudah menggunakan basis daya terpotong untuk splines regresi kubik, menggunakan rmspaket R. Setelah Anda cocok dengan model, Anda dapat mengambil representasi aljabar dari fungsi spline yang pas menggunakan fungsi Functionatau latexdi rms.

Frank Harrell
sumber
Terima kasih. Saya benar-benar membaca tanggapan Anda di sini stats.stackexchange.com/questions/67607/… sebelum memposting. Saya kira saya hanya perlu pemahaman yang lebih baik tentang apa yang dapat saya lakukan dengan rms.
Eric
Dokumentasi untuk Function()tidak benar-benar mengatakan apa yang dilakukannya. Dalam kasus saya (lihat rincian tentang Rpubs rpubs.com/EmilOWK/rms_splines ), saya mendapatkan function(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.7787nilai adalah koefisien pertama dalam model, yang 245.72672kedua, dan koefisien terakhir -873.0223tidak terlihat di mana saja persamaan. Hal yang sama berlaku untuk output dari latex().
Menghapus
Functionbekerja dengan Glm()ketika Anda menggunakan rcsfungsi 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 .
Frank Harrell