Nous allons modéliser les données tongue du package survival dont voici les premiers individus.

library(survival)
library(KMsurv)
data(tongue)
head(tongue)
##   type time delta
## 1    1    1     1
## 2    1    3     1
## 3    1    3     1
## 4    1    4     1
## 5    1   10     1
## 6    1   13     1

Soit \(T_i\) la durée de survie du \(i\)-ème individu. Nous allons ici considérer le modèle statistique suivant \[ T_i \sim \text{Weibull}(\kappa, \lambda_i), \qquad \lambda_i = \exp(\beta_0 + \beta_1 X_i), \qquad i = 1, \ldots, n, \] où \(X_i = 1_{\{\text{$i$-ème individu a un cancer de type 2}\}}\).

Pour cela nous allons ajuster par maximum de vraisemblance ce modèle en créant notre propre fonction R puis nous verrons comment faire de même via le pacakge survival.

Création de notre fonction R

Après avoir relu votre cours, le code ci dessous devrait vous être limpide. Sinon, retournez à votre cours !!!

fitweibull <- function(obs, delta, covariates){
  
  ## la fonction model.matrix cree la matrice de design via une formule R
  dsgn.mat <- model.matrix(~ covariates)
  nllik <- function(par){
    kappa <- par[1]
    lambda <- exp(dsgn.mat %*% par[-1])
    
    if (kappa <= 0)
      return(1e6 * (1 - kappa))
    
    -sum(dweibull(obs[delta==1], kappa, lambda[delta==1], log = TRUE)) -
      sum(pweibull(obs[delta==0], kappa, lambda[delta==0], log = TRUE, lower.tail = FALSE))
  }
  
  init <- rep(0, 1 + ncol(dsgn.mat))##kappa puis les beta_i
  init[1:2] <- c(1, log(mean(obs)))## estimation par methode des moments d'une Exp(lambda), lambda constant
  opt <- nlm(nllik, init, hessian = TRUE, iterlim = 10^4)
  opt$ihessian <- solve(opt$hessian)
  opt$std.err <- sqrt(diag(opt$ihessian))
  
  return(opt)
}

Nous pouvons exécuter notre code sur les données tongue, i.e.,

fit <- fitweibull(tongue$time, tongue$delta, tongue$type)
fit
## $minimum
## [1] 298.9146
## 
## $estimate
## [1]  0.8059018  5.6404967 -0.6689499
## 
## $gradient
## [1] -7.730705e-06  6.348972e-07  5.627498e-06
## 
## $hessian
##            [,1]      [,2]      [,3]
## [1,] 119.588435  6.386282  3.718498
## [2,]   6.386282 34.406754 48.696674
## [3,]   3.718498 48.696674 77.276940
## 
## $code
## [1] 1
## 
## $iterations
## [1] 13
## 
## $ihessian
##              [,1]         [,2]         [,3]
## [1,]  0.008694476 -0.009449086  0.005536046
## [2,] -0.009449086  0.279074768 -0.175406487
## [3,]  0.005536046 -0.175406487  0.123207865
## 
## $std.err
## [1] 0.09324417 0.52827528 0.35100978

Attention, ici chaque élement de la sortie précédente devrait être parfaitement compris !!!

Nous pouvons ainsi prédire facilement des quantiles voire tracer la fonction de survie estimée…

prob <- c(0.25, 0.5, 0.75)

fittedQuantileFct <- function(prob, type = 1)
  qweibull(prob, fit$estimate[1], exp(fit$estimate[2] + fit$estimate[3] * type))

fittedQuantileFct(prob)
## [1]  30.74039  91.53859 216.33963
fittedQuantileFct(prob, type = 2)
## [1]  15.74665  46.89029 110.81915
## Faisons une estimation non parametrique via Kaplan Meier
plot(survfit(Surv(time,delta)~type, data = tongue), col = 1:2)

## Puis comparons cette derniere a la modelisation paramétrique
prob <- seq(0, 1, length = 500)
lines(fittedQuantileFct(prob), 1 - prob)
lines(fittedQuantileFct(prob, type = 2), 1 - prob, col = 2)

Nous pouvons voir que la modélisation paramétrique et celle via Kaplan Meier sont sensiblement les mêmes là où nous avaons des observations. L’approche paramétrique ayant l’avantage de pouvoir extrapoler de manière plutôt raisonnable, i.e., estimer des probabilités de survie en dehors de notre domaine d’observation comme par exemple \(S(350)\) là où la version non paramétrique fournit une estimation constante.

Modélisation via le package survival

Voilà comment utiliser le package survival pour faire ce que nous avons fait précédement.

fit2 <- survreg(Surv(time,delta)~type, data = tongue)
fit2
## Call:
## survreg(formula = Surv(time, delta) ~ type, data = tongue)
## 
## Coefficients:
## (Intercept)        type 
##   5.6405171  -0.6689624 
## 
## Scale= 1.240846 
## 
## Loglik(model)= -298.9   Loglik(intercept only)= -300.7
##  Chisq= 3.58 on 1 degrees of freedom, p= 0.059 
## n= 80

Nous voyons que nous obtenons les mêmes résultats à la différence près d’une paramétrisation différente (et un optimiseur numérique différent). En effet, remarquons que

coef(fit2) - fit$estimate[-1]
##   (Intercept)          type 
##  2.046241e-05 -1.246215e-05
1/fit2$scale - fit$estimate[1]
## [1] -1.378384e-07

indiquant bien que notre fonction R fait bien la chose attendue. (Autrement dit, nous sommes trop forts !!!!)

Quant à la prédiction de quantiles, c’est bien évidemment possible

predict(fit2, type = "quantile", p = c(0.15, 0.9, 0.95), newdata = data.frame(type = c(1,2)))
##           [,1]     [,2]     [,3]
## [1,] 15.134636 406.0437 562.8410
## [2,]  7.752563 207.9917 288.3095
## Un graphique de la fonction de survie (K.-M. et Weibull)
prob <- seq(0, 1, length = 500)
plot(survfit(Surv(time,delta)~type, data = tongue), col = 1:2)
lines(predict(fit2, newdata = data.frame(type = 1), type = "quantile", p = prob), 1-prob)
lines(predict(fit2, newdata = data.frame(type = 2), type = "quantile", p = prob), 1-prob, col = 2)

Au fait… Le type de cellules cancéreuses joue-t-il un rôle ?

La variable type est elle statistiquement significative ? Pour rappel nous pouvons utiliser le test du log-rang pour cela

survdiff(Surv(time,delta)~type, data = tongue)
## Call:
## survdiff(formula = Surv(time, delta) ~ type, data = tongue)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## type=1 52       31     36.6     0.843      2.79
## type=2 28       22     16.4     1.873      2.79
## 
##  Chisq= 2.8  on 1 degrees of freedom, p= 0.0949

Ici puisque nous avons utilisé une approche via la vraisemblance nous pouvons faire par exemple un test du rapport de vraisemblance pour tester \[ H_0\colon \beta_1 = 0 \qquad H_1\colon \beta_1 \neq 0. \] A la main le test du rapport de vraisemblance pourrait se faire ainsi

fit3 <- survreg(Surv(time,delta)~1, data = tongue)
Tobs <- as.numeric(2 * (logLik(fit2) - logLik(fit3)))
Tobs
## [1] 3.579104
pval <- pchisq(Tobs, 1, lower.tail = FALSE)
print(pval)
## [1] 0.0585107

Ici la décision du test (du moins à 5%) est quelque peu ambigue puisque nous sommes assez proche des 5%. Néanmoins si nous devions tout de même conclure nous pourrions dire que nous sommes pas en capacité de rejeter \(H_0\) au profit de \(H_1\), i.e., que le type de cellule cancéreuse n’aurait pas effet sur la durée de survie du cancer de la langue. (ceci est en accord avec le test du log-rang précédent d’ailleurs)

Le résultat de ce test du rapport de vraisemblance était déjà présent dans la sortie de notre objet fit2 (tout à la fin)

fit2
## Call:
## survreg(formula = Surv(time, delta) ~ type, data = tongue)
## 
## Coefficients:
## (Intercept)        type 
##   5.6405171  -0.6689624 
## 
## Scale= 1.240846 
## 
## Loglik(model)= -298.9   Loglik(intercept only)= -300.7
##  Chisq= 3.58 on 1 degrees of freedom, p= 0.059 
## n= 80

On aurait pu faire un test de Wald à la main également via

Tobs <- coef(fit2)["type"] / sqrt(diag(fit2$var)["type"])
Tobs
##      type 
## -1.906323
pval <- 2 * pnorm(abs(Tobs), lower.tail = FALSE)
pval
##      type 
## 0.0566083

On obtient, comme souvent, à peu près les mêmes p-valeurs. Via le package survival on le fait en analysant la ligne type de la sortie suivante

summary(fit2)
## 
## Call:
## survreg(formula = Surv(time, delta) ~ type, data = tongue)
##              Value Std. Error     z        p
## (Intercept)  5.641      0.528 10.68 1.24e-26
## type        -0.669      0.351 -1.91 5.66e-02
## Log(scale)   0.216      0.116  1.87 6.21e-02
## 
## Scale= 1.24 
## 
## Weibull distribution
## Loglik(model)= -298.9   Loglik(intercept only)= -300.7
##  Chisq= 3.58 on 1 degrees of freedom, p= 0.059 
## Number of Newton-Raphson Iterations: 5 
## n= 80

A vous de jouer

  1. A quoi correspond la ligne log(scale) de la dernière sortie ?
  2. Tester si un modèle exponentiel ne serait pas plus adapté ici…
  3. Refaites une analyse en considérant un autre modèle paramétrique.