Plot geschätzte HR von Coxph-Objekt mit zeitabhängigem Koeffizienten und Splines

Question

Ich möchte die geschätzte Hazard Ratio als Funktion der Zeit im Fall eines coxph Modells mit einem zeitabhängigen Koeffizienten, der auf einem Spline-Begriff basiert, darstellen. Ich habe die zeitabhängige Koeffizient Funktion tt, analog zu diesem Beispiel, das gerade kommt aus ?coxph:

# Fit a time transform model using current age 
cox = coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung, 
    tt=function(x,t,...) pspline(x + t/365.25)) 

Aufruf survfit(cox) führt zu einem Fehler, dass survfit nicht verstehen Modelle mit tt Begriff (as described in 2011 by Terry Therneau).

Sie können den linearen Prädiktor unter Verwendung cox$linear.predictors extrahieren, aber ich würde irgendwie Zeiten und weniger trivial extrahieren müssen, Zeiten, um mit jedem zu gehen. Da tt die Datenmenge nach Ereigniszeiten aufteilt, kann ich die Spalten des Eingabedatenrahmens nicht einfach mit der Ausgabe coxph verknüpfen. Zusätzlich möchte ich die geschätzte Funktion selbst darstellen, nicht nur die Vorhersagen für die beobachteten Datenpunkte.

Es gibt a related question mit Splines hier, aber es handelt sich nicht tt.

Bearbeiten (7/7)

Ich bin auf diese noch fest. Ich habe in diesem Objekt in der Tiefe gesucht:

spline.obj = pspline(lung$age) 
str(spline.obj) 

# something that looks very useful, but I am not sure what it is 
# cbase appears to be the cardinal knots 
attr(spline.obj, "printfun") 

function (coef, var, var2, df, history, cbase = c(43.3, 47.6, 
51.9, 56.2, 60.5, 64.8, 69.1, 73.4, 77.7, 82, 86.3, 90.6)) 
{ 
    test1 <- coxph.wtest(var, coef)$test 
    xmat <- cbind(1, cbase) 
    xsig <- coxph.wtest(var, xmat)$solve 
    cmat <- coxph.wtest(t(xmat) %*% xsig, t(xsig))$solve[2, ] 
    linear <- sum(cmat * coef) 
    lvar1 <- c(cmat %*% var %*% cmat) 
    lvar2 <- c(cmat %*% var2 %*% cmat) 
    test2 <- linear^2/lvar1 
    cmat <- rbind(c(linear, sqrt(lvar1), sqrt(lvar2), test2, 
     1, 1 - pchisq(test2, 1)), c(NA, NA, NA, test1 - test2, 
     df - 1, 1 - pchisq(test1 - test2, max(0.5, df - 1)))) 
    dimnames(cmat) <- list(c("linear", "nonlin"), NULL) 
    nn <- nrow(history$thetas) 
    if (length(nn)) 
     theta <- history$thetas[nn, 1] 
    else theta <- history$theta 
    list(coef = cmat, history = paste("Theta=", format(theta))) 
} 

Also habe ich den Knoten haben, aber ich bin immer noch nicht sicher, wie die coxph Koeffizienten bei den Knoten zu kombinieren, um zu tatsächlich die Funktion grafisch darzustellen. Jeder führt sehr geschätzt.

Answer 1

Ich denke, was Sie brauchen, kann generiert werden, indem eine Eingangsmatrix mit pspline und Matrix-Multiplikation mit den entsprechenden Koeffizienten aus dem coxph Ausgang erzeugt wird. Um die HR zu erhalten, müssen Sie dann den Exponenten nehmen.

heißt

output <- data.frame(Age = seq(min(lung$age) + min(lung$time)/365.25, 
           max(lung$age + lung$time/365.25), 
           0.01)) 
output$HR <- exp(pspline(output$Age) %*% cox$coefficients[-1] - 
       sum(cox$means[-1] * cox$coefficients[-1])) 
library("ggplot2") 
ggplot(output, aes(x = Age, y = HR)) + geom_line() 

Hinweis das Alter hier ist das Alter zum Zeitpunkt von Interesse (das heißt die Summe der Basislinie Alter und die verstrichene Zeit seit dem Eintritt in die Studie). Es muss den angegebenen Bereich verwenden, um mit den Parametern im ursprünglichen Modell übereinzustimmen.

cox <- coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung, 
      tt=function(x,t,...) pspline(x + t/365.25), x = TRUE) 
index <- as.numeric(unlist(lapply(strsplit(rownames(cox$x), "\\."), "[", 1))) 
ages <- lung$age[index] 
output2 <- data.frame(Age = ages + cox$y[, 1]/365.25, 
         HR = exp(cox$x[, -1] %*% cox$coefficients[-1] - 
           sum(cox$means[-1] * cox$coefficients[-1])))