Współczynniki zależne od czasu w R - jak to zrobić?

Aktualizacja : Przepraszam za kolejną aktualizację, ale znalazłem kilka możliwych rozwiązań z ułamkami wielomianów i konkurencyjnym pakietem ryzyka, z którym potrzebuję pomocy.

Problem

Nie mogę znaleźć łatwego sposobu na wykonanie analizy współczynnika zależnego od czasu w R. Chcę mieć mój współczynnik zmiennych i zrobić to w zależności od czasu (nie zmienny), a następnie wykreślić zmienność względem czasu:

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

Możliwe rozwiązania

1) Podział zestawu danych

Patrzyłem na ten przykład (Se część 2 sesji laboratoryjnej), ale utworzenie osobnego zestawu danych wydaje się skomplikowane, kosztowne obliczeniowo i niezbyt intuicyjne ...

2) Modele z obniżoną rangą - Pakiet Coxvc

Pakiet coxvc zapewnia elegancki sposób radzenia sobie z tym problemem - oto instrukcja . Problem polega na tym, że autor nie pracuje już nad pakietem (ostatnia wersja jest od 23.05.2007), po kilku rozmowach e-mailowych udało mi się uruchomić pakiet, ale jedno uruchomienie zajęło 5 godzin na moim zestawie danych (140 000 wpisy) i podaje skrajne szacunki na koniec okresu. Możesz znaleźć nieco zaktualizowany pakiet tutaj - głównie właśnie zaktualizowałem funkcję drukowania.

Może to być tylko kwestia ulepszenia, ale ponieważ oprogramowanie nie zapewnia łatwych przedziałów ufności, a proces jest tak czasochłonny, że szukam teraz innych rozwiązań.

3) Pakiet timereg

Imponujący pakiet timereg również rozwiązuje ten problem, ale nie jestem pewien, jak go używać i nie zapewnia mi to płynnego działania.

4) Model ułamkowego czasu wielomianowego (FPT)

Znalazłem doskonałą rozprawę Aniki Buchholz na temat „Oceny zmieniających się w czasie długoterminowych efektów terapii i czynników prognostycznych” , która doskonale sprawdza się w różnych modelach. Dochodzi do wniosku, że FPT zaproponowane przez Sauerbrei i wsp. Wydaje się najbardziej odpowiednie dla współczynników zależnych od czasu:

FPT jest bardzo dobry w wykrywaniu efektów zmieniających się w czasie, podczas gdy podejście z obniżoną rangą daje w wyniku zbyt złożone modele, ponieważ nie obejmuje wyboru efektów zmieniających się w czasie.

Badania wydają się bardzo kompletne, ale dla mnie są nieco poza zasięgiem. Zastanawiam się również, skoro akurat współpracuje z Sauerbrei. Wydaje się, że jest to poprawne i myślę, że analiza mogłaby zostać przeprowadzona z pakietem mfp, ale nie jestem pewien, jak to zrobić .

5) Pakiet cmprsk

Zastanawiałem się nad analizą ryzyka konkurencji, ale obliczenia były czasochłonne, więc przełączyłem się na zwykłą regresję Coxa. CRR ma thoug opcję współzmiennych zależnych czasowych:

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

Jest kwadratowy przykład, ale nie do końca śledzę, gdzie faktycznie pojawia się czas i nie jestem pewien, jak go wyświetlić. Obejrzałem również plik test.R, ale przykład jest w zasadzie taki sam ...

Mój przykładowy kod

Oto przykład, którego używam do testowania różnych możliwości

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

Kod daje następujące wykresy: Porównanie różnych ustawień dla Coxvc i Coxvc i wykresów timecox . Myślę, że wyniki są w porządku, ale nie sądzę, że będę w stanie wyjaśnić wykres timecox - wydaje się skomplikowany ...

Moje (aktualne) pytania

• Jak przeprowadzić analizę FPT w języku R?
• Jak korzystać z zmiennej czasowej w cmprsk?
• Jak wykreślić wynik (najlepiej z przedziałami ufności)?

@mpiktas zbliżył się do oferowania wykonalnego modelu, jednak termin, który należy zastosować dla kwadratyki w czasie = t, byłby I(t^2)). Jest tak, ponieważ w R interpretacja wzoru „^” tworzy interakcje i nie wykonuje potęgowania, więc interakcja „t” z „t” jest po prostu „t”. (Czy nie należy tego migrować do SO za pomocą znacznika [r]?)

Aby zapoznać się z alternatywami dla tego procesu, który wydaje mi się nieco wątpliwy dla celów wnioskowania, i który prawdopodobnie odpowiada twojemu zainteresowaniu korzystaniem z funkcji pomocniczych w pakietach rms / Hmisc Harrella, zobacz „Strategie modelowania regresji” Harrella. Wspomina (ale tylko mimochodem, choć cytuje niektóre swoje własne prace), tworząc pasowanie splajnu do skali czasu, aby modelować zagrożenia w kształcie wanny. W rozdziale dotyczącym parametrycznych modeli przeżycia opisano różnorodne techniki kreślenia służące do sprawdzania założeń dotyczących proporcjonalnych zagrożeń i badania liniowości szacowanych skutków logarytmicznych w skali czasu.

Edycja: Dodatkową opcją jest użycie coxphparametru „tt” opisanego jako „opcjonalna lista funkcji przekształcania czasu”.

Zmieniłem odpowiedź na to, ponieważ ani @ DWin, ani @ Zach nie w pełni odpowiadają, w jaki sposób modelować zmienne w czasie współczynniki. Niedawno napisałem o tym post . Oto sedno tego.

$h(t)$

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

$time_0$$S(t)$

Zezwalając podmiotom na wejście w innym punkcie czasowym, musimy zmienić Survz Surv(time, status)na Surv(start_time, end_time, status). Chociaż end_timejest silnie skorelowany z wynikiem, start_timejest teraz dostępny jako termin interakcji (tak jak wskazano w oryginalnych sugestiach). W normalnym ustawieniu start_timewynosi 0, z wyjątkiem kilku obiektów, które pojawią się później, ale jeśli podzielimy każdą obserwację na kilka okresów, nagle mamy wiele początkowych czasów, które nie są zerowe. Jedyną różnicą jest to, że każda obserwacja występuje wielokrotnie, a wszystkie obserwacje oprócz ostatniej mają opcję nieocenzurowanego wyniku.

Podział czasu w praktyce

Właśnie opublikowałem w CRAN pakiet Grega , który ułatwia ten podział czasu . Najpierw zacznijmy od kilku obserwacji teoretycznych:

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

Możemy to pokazać graficznie z * będącym wskaźnikiem zdarzenia:

Jeśli zastosujemy, timeSplitterco następuje:

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

Otrzymujemy następujące:

Jak widać, każdy obiekt został podzielony na wiele zdarzeń, w których ostatni przedział czasu zawiera rzeczywisty stan zdarzenia. To pozwala nam teraz budować modele, które mają proste :warunki interakcji (nie używaj, *ponieważ to się rozwija time + var + time:vari nie jesteśmy zainteresowani czasem per se). Nie ma potrzeby korzystania z I()funkcji, chociaż jeśli chcesz sprawdzić nieliniowość z czasem, często tworzę osobną zmienną interakcji czasowej, do której dodam splajn, a następnie wyświetlam ją rms::contrast. W każdym razie połączenie regresji powinno teraz wyglądać następująco:

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

Korzystanie z ttfunkcji pakietu przetrwania

Istnieje również sposób na modelowanie współczynników zależnych od czasu bezpośrednio w pakiecie przeżycia za pomocą ttfunkcji. Prof. Therneau zapewnia dokładne wprowadzenie do tematu w swojej winiecie . Niestety w dużych zestawach danych jest to trudne ze względu na ograniczenia pamięci. Wydaje się, że ttfunkcja dzieli czas na bardzo drobne kawałki, generując w tym procesie ogromną matrycę.

Możesz skorzystać z funkcji Apply.rolling w PerformanceAnalytics, aby uruchomić regresję liniową za pomocą ruchomego okna, co pozwoli zmieniać współczynniki w czasie.

Na przykład:

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

Działa to również z innymi funkcjami.