Tracciamento della FC stimata dall'oggetto coxph con coefficiente e spline dipendenti dal tempo

Question

Voglio rappresentare il rapporto di rischio stimato in funzione del tempo nel caso di un modello coxph con un coefficiente tempo dipendente che si basa su un termine spline. Ho creato il coefficiente dipendente dal tempo utilizzando funzione tt, analogo a questo esempio che arriva direttamente da ?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)) 

Calling survfit(cox) risultati in un errore che survfit non comprende modelli con durata tt (as described in 2011 by Terry Therneau).

È possibile estrarre il predittore lineare utilizzando cox$linear.predictors, ma in qualche modo avrei bisogno di estrarre le età e meno banalmente, le volte in cui andare con ciascuna. Poiché tt divide il set di dati sugli orari degli eventi, non posso semplicemente abbinare le colonne del frame di dati di input con l'output coxph. Inoltre, mi piacerebbe molto tracciare la funzione stimata stessa, non solo le previsioni per i punti di dati osservati.

C'è a related question che coinvolge spline qui, ma non coinvolge tt.

Edit (7/7)

sto ancora bloccato su questo. Ho cercato in profondità a questo oggetto:

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

così, ho i nodi, ma non sono ancora sicuro di come combinare le coxph coefficienti con i nodi al fine di tracciare in realtà la funzione. Qualsiasi vantaggio molto apprezzato.

Answer 1

Credo ciò è necessario può essere generato da generare una matrice di ingresso utilizzando pspline e matrice moltiplicando questo per i coefficienti rilevanti dalla coxph uscita. Per ottenere l'HR, devi quindi prendere l'esponente.

cioè

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

nota l'età qui è l'età al momento di interesse (vale a dire la somma delle età al basale e il tempo trascorso dal momento dell'ingresso nello studio). Deve utilizzare l'intervallo specificato per corrispondere ai parametri nel modello originale. Potrebbe anche essere calcolato utilizzando il x uscita da usare x = TRUE come mostrato:

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