Coefficients dépendant du temps dans R - comment faire?

Mise à jour : Désolé pour une autre mise à jour, mais j'ai trouvé des solutions possibles avec des polynômes fractionnaires et le package de risques concurrents avec lesquels j'ai besoin d'aide.

Le problème

Je ne trouve pas de moyen facile de faire une analyse de coefficient dépendant du temps dans R. Je veux pouvoir prendre mon coefficient de variables et le faire dans un coefficient dépendant du temps (non variable) et puis tracer la variation en fonction du temps:

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

Solutions possibles

1) Fractionnement de l'ensemble de données

J'ai regardé cet exemple (voir la partie 2 de la session de laboratoire) mais la création d'un ensemble de données séparé semble compliqué, coûteux en calcul et peu intuitif ...

2) Modèles à rang réduit - Le package coxvc

Le package coxvc fournit une manière élégante de traiter le problème - voici un manuel . Le problème est que l'auteur ne développe plus le package (la dernière version est depuis le 23/05/2007), après une conversation par e-mail, j'ai fait fonctionner le package mais une exécution a pris 5 heures sur mon jeu de données (140 000 entrées) et donne des estimations extrêmes à la fin de la période. Vous pouvez trouver un package légèrement mis à jour ici - je viens surtout de mettre à jour la fonction de tracé.

Il pourrait être juste une question de peaufinage, mais comme le logiciel ne fournit pas facilement des intervalles de confiance et que le processus prend tellement de temps, je regarde en ce moment d'autres solutions.

3) Le paquet timereg

L'impressionnant package timereg résout également le problème, mais je ne sais pas comment l'utiliser et il ne me donne pas une intrigue fluide.

4) Modèle de temps polynomial fractionné (FPT)

J'ai trouvé l'excellente dissertation d'Anika Buchholz sur «L'évaluation des effets à long terme variables des thérapies et des facteurs pronostiques» qui fait un excellent travail couvrant différents modèles. Elle conclut que le FPT proposé par Sauerbrei et al semble être le plus approprié pour les coefficients dépendants du temps:

Le FPT est très bon pour détecter les effets variant dans le temps, tandis que l'approche à rang réduit donne des modèles beaucoup trop complexes, car elle n'inclut pas la sélection d'effets variant dans le temps.

La recherche semble très complète mais elle est légèrement hors de portée pour moi. Je me demande aussi un peu puisqu'elle travaille avec Sauerbrei. Cela semble solide et je suppose que l'analyse pourrait être effectuée avec le package mfp mais je ne sais pas comment.

5) Le package cmprsk

J'ai pensé à faire mon analyse des risques concurrents, mais les calculs ont pris trop de temps, alors je suis passé à la régression cox régulière. Le CRR a une option pour les covariables dépendantes du temps:

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

Il y a l'exemple quadratique mais je ne suis pas tout à fait suivre où l'heure apparaît réellement et je ne sais pas comment l'afficher. J'ai également regardé le fichier test.R mais l'exemple il est fondamentalement le même ...

Mon exemple de code

Voici un exemple que j'utilise pour tester les différentes possibilités

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

Le code donne ces graphiques: Comparaison des différents paramètres pour coxvc et des graphiques coxvc et timecox . Je suppose que les résultats sont corrects, mais je ne pense pas que je serai en mesure d'expliquer le graphique de timecox - il semble complexe ...

Mes questions (actuelles)

• Comment faire l'analyse FPT dans R?
• Comment utiliser la covariable de temps dans cmprsk?
• Comment puis-je tracer le résultat (de préférence avec des intervalles de confiance)?

@mpiktas est venu près en offrant un modèle réalisable, mais le terme qui doit être utilisé pour le quadratique dans time = t serait I(t^2)). En effet, dans R, l'interprétation de la formule de "^" crée des interactions et n'effectue pas d'exponentiation, donc l'interaction de "t" avec "t" est simplement "t". (Cela ne devrait-il pas être migré vers SO avec une balise [r]?)

Pour des alternatives à ce processus, qui me semble quelque peu douteux à des fins d'inférence, et qui correspond probablement à votre intérêt à utiliser les fonctions de support dans les packages rms / Hmisc de Harrell, voir Harrell's "Regression Modeling Strategies". Il mentionne (mais seulement en passant bien qu'il cite certains de ses propres articles) la construction d'ajustements de spline à l'échelle de temps pour modéliser des dangers en forme de baignoire. Son chapitre sur les modèles de survie paramétriques décrit une variété de techniques de traçage pour vérifier les hypothèses de risques proportionnels et pour examiner la linéarité des effets log-aléatoires estimés sur l'échelle de temps.

Edit: Une option supplémentaire consiste à utiliser coxphle paramètre 's' tt 'décrit comme une "liste facultative de fonctions de transformation temporelle".

J'ai changé la réponse à cela car ni les réponses de @ DWin ni de @ Zach ne répondent pleinement à la façon de modéliser les coefficients variant dans le temps. J'ai récemment écrit un article à ce sujet. Voici l'essentiel.

$h(t)$

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

$time_0$$S(t)$

Lorsque vous autorisez des sujets à entrer à d'autres moments, nous devons changer le Survde Surv(time, status)vers Surv(start_time, end_time, status). Bien que le end_timesoit fortement corrélé au résultat, il start_timeest désormais disponible en tant que terme d'interaction (tout comme l'indiquaient les suggestions originales). Dans un cadre normal, le start_timeest 0, sauf pour quelques sujets qui apparaissent plus tard, mais si nous divisons chaque observation en plusieurs périodes, nous avons soudainement de nombreuses heures de début qui ne sont pas nulles. La seule différence est que chaque observation se produit plusieurs fois où toutes, sauf la dernière observation, ont la possibilité d'un résultat non censuré.

Le temps partagé en pratique

Je viens de publier sur CRAN le package Greg qui facilite ce time-split . Commençons d'abord par quelques observations théoriques:

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

Nous pouvons le montrer graphiquement avec * comme indicateur d'événement:

Si nous appliquons ce timeSplitterqui suit:

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

Nous obtenons ce qui suit:

Comme vous pouvez le voir, chaque objet a été divisé en plusieurs événements dont la dernière période contient l'état réel de l'événement. Cela nous permet maintenant de construire des modèles qui ont des :termes d'interaction simples (n'utilisez pas le *comme cela se développe time + var + time:varet nous ne sommes pas intéressés par le temps en soi). Il n'est pas nécessaire d'utiliser la I()fonction, bien que si vous souhaitez vérifier la non-linéarité avec le temps, je crée souvent une variable d'interaction temporelle distincte à laquelle j'ajoute une spline, puis affiche à l'aide de rms::contrast. Quoi qu'il en soit, votre appel de régression devrait maintenant ressembler à ceci:

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

Utilisation de la ttfonction du package de survie

Il existe également un moyen de modéliser les coefficients dépendants du temps directement dans le package de survie à l'aide de la ttfonction. Le professeur Therneau fournit une introduction approfondie au sujet dans sa vignette . Malheureusement, dans les grands ensembles de données, cela est difficile à faire en raison des limitations de mémoire. Il semble que la ttfonction divise le temps en morceaux très fins générant dans le processus une énorme matrice.

Vous pouvez utiliser la fonction apply.rolling dans PerformanceAnalytics pour exécuter une régression linéaire à travers une fenêtre mobile, ce qui permettra à vos coefficients de varier dans le temps.

Par exemple:

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

Cela fonctionne également avec d'autres fonctions.