Zeitabhängige Koeffizienten in R - wie geht das?

Update : Es tut mir leid für ein weiteres Update, aber ich habe einige mögliche Lösungen mit gebrochenen Polynomen und dem konkurrierenden Risikopaket gefunden, bei denen ich Hilfe benötige.

Das Problem

Ich kann keine einfache Möglichkeit finden, eine zeitabhängige Koeffizientenanalyse in R durchzuführen. Ich möchte in der Lage sein, meinen Variablenkoeffizienten in einen zeitabhängigen Koeffizienten (nicht variabel) umzuwandeln und dann die Variation gegen die Zeit aufzuzeichnen:

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

Mögliche Lösungen

1) Aufteilen des Datensatzes

Ich habe mir dieses Beispiel angesehen (siehe Teil 2 der Laborsitzung), aber die Erstellung eines separaten Datensatzes scheint kompliziert, rechenintensiv und nicht sehr intuitiv zu sein ...

2) Modelle mit reduziertem Rang - Das coxvc-Paket

Das coxvc-Paket bietet eine elegante Möglichkeit, mit dem Problem umzugehen - hier ist ein Handbuch . Das Problem ist, dass der Autor das Paket nicht mehr entwickelt (letzte Version ist seit dem 23.05.2007). Nach einigen E-Mail-Gesprächen habe ich das Paket zum Laufen gebracht, aber ein Durchlauf hat 5 Stunden für meinen Datensatz (140 000) gedauert Einträge) und gibt extreme Schätzungen am Ende des Zeitraums. Ein leicht aktualisiertes Paket finden Sie hier - ich habe gerade die Plot-Funktion aktualisiert.

Es könnte nur eine Frage der Optimierung sein, aber da die Software nicht leicht Konfidenzintervalle bereitstellt und der Prozess so zeitaufwendig ist, suche ich gerade nach anderen Lösungen.

3) Das timereg-Paket

Das beeindruckende Timereg-Paket behebt auch das Problem, aber ich bin mir nicht sicher, wie ich es verwenden soll, und es gibt mir keine glatte Handlung.

4) Fractional Polynomial Time (FPT) -Modell

Ich fand Anika Buchholz 'exzellente Dissertation "Bewertung zeitlich variierender Langzeiteffekte von Therapien und prognostischen Faktoren" , die hervorragende Arbeit leistet und verschiedene Modelle abdeckt. Sie kommt zu dem Schluss, dass das von Sauerbrei et al. Vorgeschlagene FPT für zeitabhängige Koeffizienten am besten geeignet zu sein scheint:

FPT kann sehr gut zeitvariable Effekte erkennen, während der Reduced Rank-Ansatz zu viel zu komplexen Modellen führt, da er keine Auswahl zeitvariabler Effekte enthält.

Die Nachforschungen scheinen sehr vollständig zu sein, sind für mich jedoch leicht unerreichbar. Ich wundere mich auch ein bisschen, dass sie zufällig mit Sauerbrei arbeitet. Es scheint zwar vernünftig und ich denke, die Analyse könnte mit dem MFP-Paket durchgeführt werden, aber ich bin nicht sicher, wie.

5) Das cmprsk-Paket

Ich habe darüber nachgedacht, meine konkurrierende Risikoanalyse durchzuführen, aber die Berechnungen waren zu zeitaufwändig, sodass ich auf die reguläre Cox-Regression umgestiegen bin. Die CRR bietet eine Option für zeitabhängige Kovariaten:

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

Es gibt das quadratische Beispiel, aber ich verfolge nicht genau, wo die Uhrzeit tatsächlich angezeigt wird, und ich bin nicht sicher, wie ich sie anzeigen soll. Ich habe mir auch die test.R-Datei angeschaut, aber das Beispiel dort ist im Grunde das gleiche ...

Mein Beispielcode

Hier ist ein Beispiel, anhand dessen ich die verschiedenen Möglichkeiten teste

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

Der Code führt zu folgenden Diagrammen: Vergleich der verschiedenen Einstellungen für coxvc sowie für coxvc- und timecox- Diagramme. Ich denke, die Ergebnisse sind in Ordnung, aber ich glaube nicht, dass ich das Timecox-Diagramm erklären kann - es scheint zu komplex ...

Meine (aktuellen) Fragen

• Wie mache ich die FPT-Analyse in R?
• Wie verwende ich die Zeitkovariate in cmprsk?
• Wie zeichne ich das Ergebnis (vorzugsweise mit Konfidenzintervallen)?

@mpiktas kam dem Angebot eines realisierbaren Modells sehr nahe, aber der Begriff, der für das quadratische in time = t verwendet werden muss, wäre I(t^2)). Dies ist so, weil in R die Formelinterpretation von "^" Wechselwirkungen erzeugt und keine Exponentiation durchführt, so dass die Wechselwirkung von "t" mit "t" nur "t" ist. (Sollte dies nicht mit einem [r] -Tag nach SO migriert werden?)

Alternativen zu diesem Prozess, der für Inferenzzwecke etwas zweifelhaft erscheint und der wahrscheinlich zu Ihrem Interesse an der Verwendung der unterstützenden Funktionen in Harrells rms / Hmisc-Paketen passt, finden Sie unter Harrells "Regressionsmodellierungsstrategien". Er erwähnt (aber nur beiläufig, obwohl er einige seiner eigenen Veröffentlichungen zitiert), dass die Konstruktion von Splines der Zeitskala entspricht, um badewannenförmige Gefahren zu modellieren. In seinem Kapitel über parametrische Überlebensmodelle werden verschiedene Darstellungstechniken beschrieben, mit denen Annahmen über proportionale Gefahren überprüft und die Linearität der geschätzten logarithmischen Gefährdungseffekte auf der Zeitskala untersucht werden können.

Bearbeiten: Eine zusätzliche Option ist die Verwendung coxphdes Parameters 'tt', der als "optionale Liste von Zeittransformationsfunktionen" beschrieben wird.

Ich habe die Antwort auf diese Frage geändert, da weder @ DWins noch @ Zachs Antworten vollständig darauf antworten, wie man zeitvariable Koeffizienten modelliert. Ich habe kürzlich einen Beitrag darüber geschrieben. Hier ist der Kern davon.

$h(t)$

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

$time_0$$S(t)$

Wenn wir zulassen, dass Personen zu anderen Zeitpunkten eintreten, müssen wir Survvon Surv(time, status)auf ändern Surv(start_time, end_time, status). Während das end_timestark mit dem Ergebnis korreliert, start_timeist das nun als Interaktionsbegriff verfügbar (wie in den ursprünglichen Vorschlägen angedeutet). In einer normalen Einstellung start_timeist der Wert 0, mit Ausnahme einiger später erscheinender Themen. Wenn wir jedoch jede Beobachtung in mehrere Zeiträume aufteilen, haben wir plötzlich viele Startzeiten ungleich Null. Der einzige Unterschied besteht darin, dass jede Beobachtung mehrmals auftritt, wobei alle außer der letzten Beobachtung die Option eines nicht zensierten Ergebnisses haben.

Zeitteilung in der Praxis

Ich habe gerade auf CRAN den veröffentlichten Greg Paket , das dies macht die zeitlich geteilte einfach. Zunächst beginnen wir mit einigen theoretischen Beobachtungen:

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

Wir können dies grafisch darstellen, wobei * ein Indikator für das Ereignis ist:

Wenn wir timeSplitterfolgendes anwenden :

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

Wir bekommen folgendes:

Wie Sie sehen, wurde jedes Objekt in mehrere Ereignisse aufgeteilt, wobei die letzte Zeitspanne den tatsächlichen Ereignisstatus enthält. Auf diese Weise können wir nun Modelle mit einfachen :Interaktionstermen erstellen (verwenden Sie das nicht, *wenn dies erweitert wird, time + var + time:varund wir interessieren uns nicht für die Zeit an sich). Die I()Funktion muss nicht verwendet werden. Wenn Sie jedoch die Nichtlinearität mit der Zeit überprüfen möchten, erstelle ich häufig eine separate Zeit-Interaktions-Variable, der ich einen Spline hinzufüge und dann mithilfe von anzeige rms::contrast. Wie auch immer, Ihr Regressionsaufruf sollte jetzt so aussehen:

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

Verwendung der ttFunktion des Überlebenspakets

Mit der ttFunktion können auch zeitabhängige Koeffizienten direkt im Überlebenspaket modelliert werden . Prof. Therneau gibt in seiner Vignette eine ausführliche Einführung in das Thema . Leider ist dies bei großen Datenmengen aufgrund von Speicherbeschränkungen schwierig. Es scheint, dass die ttFunktion die Zeit in sehr feine Teile aufteilt und dabei eine riesige Matrix erzeugt.

Sie können die Verwendung apply.rolling Funktion in PerformanceAnalytics eine lineare Regression durch ein Rollfenster ausführen, die Ihre Koeffizienten Zeit variieren ermöglicht.

Beispielsweise:

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

Dies funktioniert auch mit anderen Funktionen.