Koefisien tergantung waktu dalam R - bagaimana cara melakukannya?

Pembaruan : Maaf untuk pembaruan lain tetapi saya telah menemukan beberapa solusi yang mungkin dengan polinomial fraksional dan paket risiko yang bersaing yang saya perlukan bantuan.

Masalah

Saya tidak dapat menemukan cara mudah untuk melakukan analisis koefisien dependen waktu dalam R. Saya ingin dapat mengambil koefisien variabel saya dan melakukannya menjadi koefisien dependen waktu (bukan variabel) dan kemudian plot variasi terhadap waktu:

$\beta_{my\_variable}=\beta_0+\beta_1*t+\beta_2*t^2...$

Solusi yang memungkinkan

1) Memisahkan dataset

Saya telah melihat contoh ini (Se bagian 2 dari sesi lab) tetapi pembuatan dataset terpisah tampak rumit, komputasional mahal dan tidak terlalu intuitif ...

2) Model Reduced Rank - Paket coxvc

The paket coxvc menyediakan cara yang elegan untuk menangani masalah - di sini adalah panduan . Masalahnya adalah bahwa penulis tidak lagi mengembangkan paket (versi terakhir adalah sejak 23/5/2007), setelah beberapa percakapan email saya mendapatkan paket untuk bekerja tetapi satu kali berjalan membutuhkan 5 jam pada dataset saya (140 000 entri) dan memberikan perkiraan ekstrim pada akhir periode. Anda dapat menemukan paket yang sedikit diperbarui di sini - Saya baru saja memperbarui fungsi plot.

Mungkin hanya masalah penyesuaian tetapi karena perangkat lunak tidak mudah memberikan interval kepercayaan dan prosesnya sangat memakan waktu, saya sedang mencari solusi lain sekarang.

3) Paket penghitung waktu

Paket timereg yang mengesankan juga mengatasi masalah tetapi saya tidak yakin bagaimana menggunakannya dan itu tidak memberi saya alur yang mulus.

4) model Fractional Polynomial Time (FPT)

Saya menemukan disertasi Anika Buchholz yang sangat baik tentang "Penilaian efek jangka panjang yang bervariasi dari terapi dan faktor prognostik" yang melakukan pekerjaan luar biasa yang mencakup model yang berbeda. Dia menyimpulkan bahwa usulan FPT Sauerbrei et al tampaknya paling tepat untuk koefisien tergantung waktu:

FPT sangat baik dalam mendeteksi efek yang bervariasi waktu, sedangkan pendekatan Reduced Rank menghasilkan model yang terlalu rumit, karena tidak termasuk pemilihan efek yang bervariasi waktu.

Penelitian ini kelihatannya sangat lengkap tetapi sedikit di luar jangkauan saya. Saya juga sedikit bertanya-tanya karena dia kebetulan bekerja dengan Sauerbrei. Tampaknya terdengar baik dan saya kira analisis dapat dilakukan dengan paket mfp tapi saya tidak yakin bagaimana caranya.

5) Paket cmprsk

Saya sudah berpikir untuk melakukan analisis risiko saya yang bersaing, tetapi perhitungannya memakan waktu lama, jadi saya beralih ke regresi cox biasa. The crr memiliki thoug pilihan untuk waktu kovariat tergantung:

....
cov2        matrix of covariates that will be multiplied 
            by functions of time; if used, often these 
            covariates would also appear in cov1 to give 
            a prop hazards effect plus a time interaction
....

Ada contoh kuadrat tapi saya tidak cukup mengikuti di mana waktu sebenarnya muncul dan saya tidak yakin bagaimana menampilkannya. Saya juga sudah melihat file test.R tapi contohnya pada dasarnya sama ...

Contoh kode saya

Inilah contoh yang saya gunakan untuk menguji berbagai kemungkinan

library("survival")
library("timereg")
data(sTRACE)

# Basic cox regression    
surv <- with(sTRACE, Surv(time/365,status==9))
fit1 <- coxph(surv~age+sex+diabetes+chf+vf, data=sTRACE)
check <- cox.zph(fit1)
print(check)
plot(check, resid=F)
# vf seems to be the most time varying

######################################
# Do the analysis with the code from #
# the example that I've found        #
######################################

# Split the dataset according to the splitSurv() from prof. Wesley O. Johnson
# http://anson.ucdavis.edu/~johnson/st222/lab8/splitSurv.ssc
new_split_dataset = splitSuv(sTRACE$time/365, sTRACE$status==9, sTRACE[, grep("(age|sex|diabetes|chf|vf)", names(sTRACE))])

surv2 <- with(new_split_dataset, Surv(start, stop, event))
fit2 <- coxph(surv2~age+sex+diabetes+chf+I(pspline(stop)*vf), data=new_split_dataset)
print(fit2)

######################################
# Do the analysis by just straifying #
######################################
fit3 <- coxph(surv~age+sex+diabetes+chf+strata(vf), data=sTRACE)
print(fit3)

# High computational cost!
# The price for 259 events
sum((sTRACE$status==9)*1)
# ~240 times larger dataset!
NROW(new_split_dataset)/NROW(sTRACE)

########################################
# Do the analysis with the coxvc and   #
# the timecox from the timereg library #
########################################
Ft_1 <- cbind(rep(1,nrow(sTRACE)),bs(sTRACE$time/365,df=3))
fit_coxvc1 <- coxvc(surv~vf+sex, Ft_1, rank=2, data=sTRACE)

fit_coxvc2 <- coxvc(surv~vf+sex, Ft_1, rank=1, data=sTRACE)

Ft_3 <- cbind(rep(1,nrow(sTRACE)),bs(sTRACE$time/365,df=5))
fit_coxvc3 <- coxvc(surv~vf+sex, Ft_3, rank=2, data=sTRACE)

layout(matrix(1:3, ncol=1))
my_plotcoxvc <- function(fit, fun="effects"){
    plotcoxvc(fit,fun=fun,xlab='time in years', ylim=c(-1,1), legend_x=.010)
    abline(0,0, lty=2, col=rgb(.5,.5,.5,.5))
    title(paste("B-spline =", NCOL(fit$Ftime)-1, "df and rank =", fit$rank))
}
my_plotcoxvc(fit_coxvc1)
my_plotcoxvc(fit_coxvc2)
my_plotcoxvc(fit_coxvc3)

# Next group
my_plotcoxvc(fit_coxvc1)

fit_timecox1<-timecox(surv~sex + vf, data=sTRACE)
plot(fit_timecox1, xlab="time in years", specific.comps=c(2,3))

Hasil kode dalam grafik ini: Perbandingan pengaturan yang berbeda untuk coxvc dan coxvc dan plot timecox . Saya kira hasilnya ok tapi saya tidak berpikir saya akan bisa menjelaskan grafik timecox - sepertinya rumit ...

Pertanyaan saya (saat ini)

Bagaimana saya melakukan analisis FPT di R?
Bagaimana cara menggunakan kovariat waktu di cmprsk?
Bagaimana saya memplot hasilnya (lebih disukai dengan interval kepercayaan)?

r regression survival cox-model Max Gordon
sumber

y = x β_{m y}

$y=x\beta_{my}$

y = x β_{0} + x \cdot t β_{1} + x \cdot t^{2} β_{2}

$y=x\beta_0+x\cdot t\beta_1+x\cdot t^2\beta_2$ y~xy~x*(t+t^2)-ty~x+x:t+x:t^2

Saya memikirkan bagian kedua: "2. Model Kovariat yang bergantung pada Waktu untuk Memeriksa Asumsi PH" akan menjadi bagian yang berurusan dengan pertanyaan saya. Saya berharap untuk melakukan sesuatu dari rumus yang Anda gambarkan tetapi ketika saya mencobanya saya entah mendapat kesalahan atau variabel waktu yang terpisah ... Saya mendapat nilai-p rendah untuk waktu: -D

Max Gordon

@ max-gordon, apakah variabel respons Anda kuantitas, atau waktu berlalu hingga genap terjadi? Karena sebagian besar metode yang Anda kutip khusus untuk data waktu-ke-peristiwa.

f1r3br4nd

@ f1r3br4nd: Ini adalah kuantitas (Umur dalam penelitian saya) di mana bahayanya tidak proporsional, yaitu bervariasi dari waktu ke waktu dalam model waktu-ke-peristiwa saya. Pada akhirnya saya memutuskan untuk membagi menjadi dua kerangka waktu yang berbeda karena saya tidak senang dengan melakukan grafik 3-D - yang tidak akan pernah melewati pengulas ...

Max Gordon

Ada perbedaan antara prediktor dependen / bervariasi waktu dan interaksi waktu. Sebagian besar variabel tergantung waktu (jenis kelamin adalah pengecualian). Jika Anda memiliki satu pengamatan per orang, maka Anda akan memiliki sedikit atau tidak ada kesempatan untuk melakukan analisis tergantung / bervariasi waktu. Metode Anderson-Gill adalah yang paling sering digunakan untuk analisis survival tergantung waktu. Keuntungan dari metode tergantung waktu adalah bahwa nilai-nilai selama tindak lanjut mungkin lebih prediktif dari pengalaman bertahan hidup daripada nilai-nilai dasar. Situasi kedua, prediktor yang berinteraksi waktu hanyalah uji asumsi PH.

Adam Robinsson

Jawaban:

@mpiktas nyaris menawarkan model yang layak, namun istilah yang perlu digunakan untuk kuadratik dalam waktu = t adalah I(t^2)). Ini karena di R rumus interpretasi "^" menciptakan interaksi dan tidak melakukan eksponensial, sehingga interaksi "t" dengan "t" hanyalah "t". (Bukankah ini harus dimigrasi ke SO dengan tag [r]?)

Untuk alternatif dari proses ini, yang bagi saya agak meragukan untuk tujuan inferensi, dan yang mungkin sesuai dengan minat Anda dalam menggunakan fungsi-fungsi pendukung dalam paket rm / Hmisc Harrell, lihat "Strategi Regresi Modeling" Harrell. Dia menyebutkan (tetapi hanya secara sepintas meskipun dia mengutip beberapa makalahnya sendiri) membangun spline sesuai dengan skala waktu untuk memodelkan bahaya berbentuk bak mandi. Babnya tentang model survival parametrik menjelaskan berbagai teknik merencanakan untuk memeriksa asumsi bahaya proporsional dan untuk memeriksa linearitas dari perkiraan dampak bahaya log pada skala waktu.

Sunting: Opsi tambahan adalah menggunakan coxphparameter 'tt' yang dijelaskan sebagai "daftar opsional fungsi transformasi waktu."

DWIN
sumber

Saya setuju bahwa ini mungkin harus dipindahkan ke tag SO [r].

Zach

+1 untuk jawaban Anda, saya tidak sadar bahwa ini akan menjadi jawaban yang sulit. Sepertinya masalah umum, mungkin pertanyaannya lebih merupakan pertanyaan tentang pengkodean daripada dan Anda mungkin benar tentang SO menjadi pilihan yang lebih baik. Saya mencoba formula Anda tampaknya vf + I (vf log (waktu)) memiliki kesesuaian yang sangat baik, saya mencoba hanya vf waktu dan vf * waktu ^ 2 tetapi log memberikan tarif p-value terendah. Saya mencoba menjalankannya dengan fungsi cph () untuk mendapatkan AIC tetapi memberikan kesalahan :( Apakah Anda punya ide tentang bagaimana melakukan plot pada estimasi?

Max Gordon

Saya berpikir bahwa check <- cox.zph(fit1); print(check); plot(check, resid=F)seperti dalam pengaturan Anda memberi plot informatif waktu "efek". Apakah maksud Anda cph () yang berasal dari paket rms atau coxph dari survival?

DWin

Ya, residu Schoenfeld memang memberikan ide yang bagus tentang variasi waktu tetapi saya pikir orang mungkin akan kesulitan memahaminya. Plot memberi seperti yang saya mengerti variasi residu tidak dijelaskan oleh model saya. Saya ingin plot di mana saya memiliki efek variabel lengkap pada sumbu y dan waktu pada sumbu x, saya percaya bahwa ini akan lebih mudah untuk ditafsirkan karena Anda tidak perlu melihat tabel dan plot untuk mendapatkan bahaya pada titik waktu tertentu ... Ya saya maksudkan cph () dan bukan coxph () karena orang itu tidak bekerja dengan AIC ()

Max Gordon

Saya juga sedikit bingung mengapa saya menemukan semua metode kompleks yang dijelaskan dalam pertanyaan saya sementara saya ini (variabel * waktu) tampaknya sangat mudah dan intuitif - sebagai non-ahli statistik yang saya pikirkan - apa yang telah saya lewatkan ?

Max Gordon

Saya telah mengubah jawaban untuk ini karena jawaban @ DWin atau @ Zach tidak sepenuhnya menjawab bagaimana memodelkan koefisien yang bervariasi waktu. Saya baru saja menulis posting tentang ini. Inilah intinya.

$h(t)$

h (t) = \frac{f (t)}{S (t)}

$h(t) = \frac{f(t)}{S(t)}$

$f(t)$ $S(t)$ $0$ $\infty$

$time_0$ $S(t)$

Saat mengizinkan subjek masuk pada titik waktu lain, kita harus mengubah Survdari Surv(time, status)menjadi Surv(start_time, end_time, status). Sementara end_timesangat berkorelasi dengan hasil, start_timesekarang tersedia sebagai istilah interaksi (seperti yang disarankan dalam saran asli). Dalam pengaturan reguler, start_timeis 0 kecuali untuk beberapa subjek yang muncul kemudian tetapi kami membagi setiap pengamatan menjadi beberapa periode, kami tiba-tiba memiliki banyak waktu mulai yang tidak nol. Satu-satunya perbedaan adalah bahwa setiap pengamatan terjadi beberapa kali di mana semua kecuali pengamatan terakhir memiliki opsi untuk hasil yang tidak disensor.

Pemecahan waktu dalam praktik

Saya baru saja menerbitkan di CRAN paket Greg yang membuat pembagian waktu ini mudah. Pertama kita mulai dengan beberapa pengamatan teoretis:

library(Greg)
test_data <- data.frame(
  id = 1:4,
  time = c(4, 3.5, 1, 5),
  event = c("censored", "dead", "alive", "dead"),
  age = c(62.2, 55.3, 73.7, 46.3),
  date = as.Date(
    c("2003-01-01", 
      "2010-04-01", 
      "2013-09-20",
      "2002-02-23"))
)

Kami dapat menampilkan ini secara grafis dengan * menjadi indikator acara:

Jika kami menerapkan timeSplittersebagai berikut:

library(dplyr)
split_data <- 
  test_data %>% 
  select(id, event, time, age, date) %>% 
  timeSplitter(by = 2, # The time that we want to split by
               event_var = "event",
               time_var = "time",
               event_start_status = "alive",
               time_related_vars = c("age", "date"))

Kami mendapatkan yang berikut ini:

Seperti yang Anda lihat, setiap objek telah dipecah menjadi beberapa peristiwa di mana rentang waktu terakhir berisi status peristiwa aktual. Ini memungkinkan kami untuk sekarang membangun model yang memiliki :istilah interaksi sederhana (jangan gunakan *saat ini diperluas time + var + time:vardan kami tidak tertarik pada waktu per se). Tidak perlu menggunakan I()fungsi meskipun jika Anda ingin memeriksa nonlinier dengan waktu saya sering membuat variabel interaksi waktu yang terpisah yang saya tambahkan spline dan kemudian menampilkan menggunakan rms::contrast. Bagaimanapun, panggilan regresi Anda sekarang akan terlihat seperti ini:

coxp(Surv(start_time, end_time, event) ~ var1 + var2 + var2:time, 
     data = time_split_data)

Menggunakan ttfungsi paket survival

Ada juga cara untuk memodelkan koefisien tergantung waktu secara langsung dalam paket survival menggunakan ttfungsi. Prof. Therneau memberikan pengantar menyeluruh untuk subjek dalam sketsa nya . Sayangnya dalam dataset besar hal ini sulit dilakukan karena keterbatasan memori. Tampaknya ttfungsi tersebut membagi waktu menjadi bagian-bagian yang sangat halus menghasilkan dalam proses matriks besar.

Max Gordon
sumber

Anda dapat menggunakan fungsi apply.rolling di PerformanceAnalytics untuk menjalankan regresi linier melalui jendela bergulir, yang akan memungkinkan koefisien Anda bervariasi dari waktu ke waktu.

Sebagai contoh:

library(PerformanceAnalytics)
library(quantmod)
getSymbols(c('AAPL','SPY'), from='01-01-1900')
chart.RollingRegression(Cl(AAPL),Cl(SPY), width=252, attribute='Beta')
#Note: Alpha=y-intercept, Beta=regression coeffient

Ini berfungsi dengan fungsi lain juga.

Zach
sumber

Terima kasih atas jawaban Anda, saya kira jendela waktu yang bergerak harus bekerja sama baiknya dengan pendekatan saya. Saya tidak bisa mendapatkan contoh Anda untuk dijalankan, bisakah Anda memberikan contoh berdasarkan contoh sTRACE saya sehingga saya tahu persis bagaimana menerapkannya?

Max Gordon