Memprediksi respons dari kurva baru menggunakan paket fda di R

Pada dasarnya yang ingin saya lakukan adalah memprediksi respons skalar menggunakan beberapa kurva. Saya sudah sejauh melakukan regresi (menggunakan fRegress dari paket fda) tetapi tidak tahu bagaimana menerapkan hasil ke set kurva BARU (untuk prediksi).

Saya memiliki N = 536 kurva, dan 536 respons skalar. Inilah yang telah saya lakukan sejauh ini:

• Saya telah membuat dasar untuk kurva.
• Saya telah membuat objek fdPar untuk memperkenalkan penalti
• Saya telah membuat objek fd menggunakan smooth.basis untuk memuluskan kurva dengan penalti yang dipilih pada dasar yang ditentukan.
• Saya telah menjalankan regresi menggunakan fRegress (), mundur kurva pada respon skalar.

Sekarang, yang ingin saya lakukan adalah menggunakan regresi itu untuk menghasilkan prediksi untuk set data baru yang saya miliki. Sepertinya saya tidak dapat menemukan cara mudah untuk melakukan ini.

Bersulang

Saya tidak peduli dengan fdapenggunaan struktur objek Inception - like-list-dalam-daftar-dalam-daftar, tetapi tanggapan saya akan mematuhi sistem yang telah dibuat oleh penulis paket.

Saya pikir itu penting untuk pertama-tama berpikir tentang apa yang sedang kita lakukan sebenarnya. Berdasarkan uraian Anda tentang apa yang telah Anda lakukan sejauh ini, inilah yang saya yakin Anda lakukan (beri tahu saya jika saya salah mengartikan sesuatu). Saya akan terus menggunakan notasi dan, karena kurangnya data nyata, contoh dari Analisis Data Fungsional Ramsay dan Silverman dan Analisis Data Fungsional Ramsay, Hooker, dan Graves dengan R dan MATLAB (Beberapa persamaan dan kode berikut diangkat langsung dari buku-buku ini).

Kami memodelkan respons skalar melalui model linier fungsional, yaitu

$y_i = \beta_0 + \int_0^T X_i(s)\beta(s)ds + \epsilon_i$

Kami memperluas dengan dasar tertentu. Kami menggunakan, katakanlah, fungsi berbasisBegitu,$\beta$$K$

$\beta(s) = \sum\limits_{k=1}^K b_k\theta_k (s)$

Dalam notasi matriks, ini adalah .$\beta (s)=\boldsymbol{\theta}'(s)\mathbf{b}$

Kami juga memperluas fungsi kovariat dalam beberapa basis, juga (misalkan fungsi basis ). Begitu,$L$

$X_i(s) = \sum\limits_{k=1}^L c_{ik}\psi_k (s)$

Sekali lagi, dalam notasi matriks, ini adalah .$X(s)=\mathbf{C}\boldsymbol{\psi}(s)$

Dan dengan demikian, jika kita membiarkan , model kita dapat dinyatakan sebagai$\mathbf{J}=\int\boldsymbol{\psi}(s)\boldsymbol{\theta}'(s)ds$

$y = \beta_0 + \mathbf{CJb}$ .

Dan jika kita membiarkan dan , model kami adalah$\mathbf{Z} = [\mathbf{1}\quad \mathbf{CJ}]$$\boldsymbol{\xi} = [\beta_0\quad \mathbf{b}']'$

$\mathbf{y} = \mathbf{Z}\boldsymbol{\xi}$

Dan ini terlihat jauh lebih akrab bagi kita.

Sekarang saya melihat Anda menambahkan semacam regularisasi. The fdapaket bekerja dengan hukuman kekasaran bentuk

$P = \lambda\int[L\beta(s)]^2 ds$

untuk beberapa operator linear diferensial . Sekarang dapat ditunjukkan (detail ditinggalkan di sini - benar-benar tidak sulit untuk menunjukkan ini) bahwa jika kita mendefinisikan matriks penalti sebagai$L$$\mathbf{R}$

$\mathbf{R}= \lambda\left(\begin{array}{cccc} 0 & 0 & \cdots & 0 \\ 0 & R_1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & R_K \end{array}\right)$

di mana adalah dalam hal perluasan dasar , maka kami meminimalkan jumlah kotak yang terkena penalti:$R_i$$\beta_i$

$\left(\mathbf{y}-\mathbf{Z}\boldsymbol{\xi}\right)'\left(\mathbf{y}-\mathbf{Z}\boldsymbol{\xi}\right) + \lambda\boldsymbol{\xi}'\mathbf{R}\boldsymbol{\xi}$ ,

dan masalah kita hanyalah regresi ridge dengan solusi:

$\hat{\boldsymbol{\xi}}=(\mathbf{Z}'\mathbf{Z}+\lambda\mathbf{R})^{-1}\mathbf{Z}'\mathbf{y}$ .

Saya berjalan di atas karena, (1) Saya pikir penting untuk memahami apa yang kami lakukan, dan (2) beberapa di atas diperlukan untuk memahami beberapa kode yang akan saya gunakan nanti. Aktif ke kode ...

Berikut ini contoh data dengan kode R. Saya menggunakan dataset cuaca Kanada yang disediakan dalam fdapaket. Kami akan memodelkan curah hujan tahunan log untuk sejumlah stasiun cuaca melalui model linier fungsional dan kami akan menggunakan profil suhu (suhu dicatat sekali sehari selama 365 hari) dari masing-masing stasiun sebagai kovariat fungsional. Kami akan melanjutkan dengan cara yang Anda uraikan dalam situasi Anda. Data direkam di 35 stasiun. Saya akan memecah dataset menjadi 34 stasiun, yang akan digunakan sebagai data saya, dan stasiun terakhir, yang akan menjadi dataset "baru" saya.

Saya melanjutkan melalui kode R dan komentar (saya berasumsi Anda cukup akrab dengan fdapaket sehingga tidak ada yang berikut ini yang mengejutkan - jika ini tidak terjadi, tolong beri tahu saya):

# pick out data and 'new data'
dailydat <- daily$precav[,2:35]
dailytemp <- daily$tempav[,2:35]
dailydatNew <- daily$precav[,1]
dailytempNew <- daily$tempav[,1]

# set up response variable
annualprec <- log10(apply(dailydat,2,sum))

# create basis objects for and smooth covariate functions
tempbasis <- create.fourier.basis(c(0,365),65)
tempSmooth <- smooth.basis(day.5,dailytemp,tempbasis)
tempfd <- tempSmooth$fd

# create design matrix object
templist <- vector("list",2)
templist[[1]] <- rep(1,34)
templist[[2]] <- tempfd

# create constant basis (for intercept) and
# fourier basis objects for remaining betas
conbasis <- create.constant.basis(c(0,365))
betabasis <- create.fourier.basis(c(0,365),35)
betalist <- vector("list",2)
betalist[[1]] <- conbasis
betalist[[2]] <- betabasis

# set roughness penalty for betas 
Lcoef <- c(0,(2*pi/365)^2,0)
harmaccelLfd <- vec2Lfd(Lcoef, c(0,365))
lambda <- 10^12.5
betafdPar <- fdPar(betabasis, harmaccelLfd, lambda)
betalist[[2]] <- betafdPar

# regress
annPrecTemp <- fRegress(annualprec, templist, betalist)

Sekarang ketika saya pertama kali diajarkan tentang data fungsional setahun yang lalu, saya bermain-main dengan paket ini. Saya juga tidak bisa predict.fRegressmemberikan apa yang saya inginkan. Melihat ke belakang sekarang, saya masih tidak tahu bagaimana membuatnya berperilaku. Jadi, kita hanya perlu mendapatkan prediksi semi-manual. Saya akan menggunakan potongan yang saya tarik langsung dari kode fRegress(). Sekali lagi, saya melanjutkan melalui kode dan komentar.

Pertama, pengaturannya:

# create basis objects for and smooth covariate functions for new data
tempSmoothNew <- smooth.basis(day.5,dailytempNew,tempbasis)
tempfdNew <- tempSmoothNew$fd

# create design matrix object for new data
templistNew <- vector("list",2)
templistNew[[1]] <- rep(1,1)
templistNew[[2]] <- tempfdNew

# convert the intercept into an fd object
onebasis <- create.constant.basis(c(0,365))
templistNew[[1]] <- fd(matrix(templistNew[[1]],1,1), onebasis)

Sekarang untuk mendapatkan prediksi

$\hat{\mathbf{y}}_{\mathrm{new}}=\mathbf{Z}_{\mathrm{new}}\hat{\boldsymbol{\xi}}$

Saya hanya mengambil kode yang fRegressdigunakan untuk menghitung yhatfdobjdan mengeditnya sedikit. fRegressmenghitung yhatfdobjdengan memperkirakan integral melalui aturan trapesium (dengan dan diperluas di basis masing-masing). $\int_0^T X_i (s) \beta(s)$$X_i$$\beta$

Biasanya, fRegresshitung nilai-nilai yang dipasang dengan mengulang melalui kovariat yang disimpan annPrecTemp$xfdlist. Jadi untuk masalah kami, kami mengganti daftar kovariat ini dengan yang sesuai dalam daftar kovariat baru kami, yaitu templistNew,. Berikut kodenya (identik dengan kode yang ditemukan fRegressdengan dua pengeditan, beberapa penghapusan kode yang tidak dibutuhkan, dan beberapa komentar ditambahkan):

# set up yhat matrix (in our case it's 1x1)
yhatmat <- matrix(0,1,1)

# loop through covariates
p <- length(templistNew)
for(j in 1:p){
    xfdj       <- templistNew[[j]]
    xbasis     <- xfdj$basis
    xnbasis    <- xbasis$nbasis
    xrng       <- xbasis$rangeval
    nfine      <- max(501,10*xnbasis+1)
    tfine      <- seq(xrng[1], xrng[2], len=nfine)
    deltat     <- tfine[2]-tfine[1]
    xmat       <- eval.fd(tfine, xfdj)
    betafdParj <- annPrecTemp$betaestlist[[j]]
    betafdj    <- betafdParj$fd
    betamat    <- eval.fd(tfine, betafdj)
    # estimate int(x*beta) via trapezoid rule
    fitj       <- deltat*(crossprod(xmat,betamat) - 
                      0.5*(outer(xmat[1,],betamat[1,]) +
              outer(xmat[nfine,],betamat[nfine,])))
    yhatmat    <- yhatmat + fitj
}

(catatan: jika Anda melihat kode chunk dan sekitarnya ini fRegress, Anda akan melihat langkah-langkah yang saya uraikan di atas).

Saya menguji kode dengan menjalankan kembali contoh cuaca menggunakan semua 35 stasiun sebagai data kami dan membandingkan output dari loop di atas annPrecTemp$yhatfdobjdan semuanya cocok. Saya juga menjalankannya beberapa kali menggunakan stasiun yang berbeda sebagai data "baru" saya dan semuanya tampak masuk akal.

Beri tahu saya jika ada di atas yang tidak jelas atau ada yang tidak berfungsi dengan benar. Maaf atas tanggapan yang terlalu rinci. Saya tidak bisa menahan diri :) Dan jika Anda belum memilikinya, lihat dua buku yang saya gunakan untuk menulis tanggapan ini. Buku-buku itu sangat bagus.