Examples

packages

Vous avez besoin des packages suivants.

library(lubridate)

Attachement du package : ‘lubridate’

The following object is masked from ‘package:base’:

    date

Cox model

Stanford Heart Transplant data

Survival of patients on the waiting list for the Stanford heart transplant program.

• accept.dt: acceptance into program

• fustat: dead or alive

• surgery: prior bypass surgery

• age: age (in years)

• futime: followup time

• wait.time: time before transplant

• transplant: transplant indicator

On sélectionne les variables accept.dt, fustat, surgery,age,futime,wait.time,transplant

On transforme la variable accept.dt en gardant seulement l’année.

jasa = dplyr::mutate(jasa, accept.yr = year(ymd(accept.dt)))
jasa = dplyr::select(jasa,   fustat, surgery,age,futime,wait.time,transplant,accept.yr)

head(jasa)

Un premier modèle avec les variables non-dépendantes du temps.

coxph(Surv(futime,fustat) ~ accept.yr + surgery + age, data = jasa)

Call:
coxph(formula = Surv(futime, fustat) ~ accept.yr + surgery + 
    age, data = jasa)

             coef exp(coef) se(coef)     z     p
accept.yr -0.1320    0.8764   0.0681 -1.94 0.053
surgery   -0.6427    0.5259   0.3673 -1.75 0.080
age        0.0276    1.0280   0.0134  2.06 0.039

Likelihood ratio test=14.5  on 3 df, p=0.00226
n= 103, number of events= 75 

Avec les variables dépendantes du temps.

Mauvaise paramétrisation avec la variable transplant

autoplot(survfit(Surv(futime,fustat) ~transplant , data = jasa))

coxph(Surv(futime,fustat) ~ surgery + transplant + age , data = jasa)

Call:
coxph(formula = Surv(futime, fustat) ~ surgery + transplant + 
    age, data = jasa)

              coef exp(coef) se(coef)     z       p
surgery    -0.4190    0.6577   0.3712 -1.13    0.26
transplant -1.7171    0.1796   0.2785 -6.16 7.1e-10
age         0.0589    1.0607   0.0150  3.91 9.1e-05

Likelihood ratio test=45.9  on 3 df, p=6.11e-10
n= 103, number of events= 75 

Un nouveau format pour gérer les variables “time-dependent”

Pour un individu sans transplantation

jasa[1,]

jasa1[1,]

Pour un individu avec transplantation

jasa[4,]

jasa1[jasa1$id==4,]

coxph(Surv(start, stop, event) ~ age + surgery + transplant + year ,jasa1)

Call:
coxph(formula = Surv(start, stop, event) ~ age + surgery + transplant + 
    year, data = jasa1)

              coef exp(coef) se(coef)     z     p
age         0.0272    1.0276   0.0137  1.98 0.047
surgery    -0.6371    0.5288   0.3672 -1.73 0.083
transplant -0.0129    0.9872   0.3133 -0.04 0.967
year       -0.1464    0.8638   0.0705 -2.08 0.038

Likelihood ratio test=15.1  on 4 df, p=0.00447
n= 170, number of events= 75 

Quand \(p\) grandit : on régularise, par exemple via l’elasticnet

data("nki70")

model_matrix = model.matrix( ~ as.factor(Grade)  + . - Grade - 1   , data = nki70[3:77])


X = model_matrix[,-1]


elasticnet_solution = cv.glmnet(X,Surv(nki70$time, nki70$event), family = "cox" , alpha = 0.5,penalty.factor = c(rep(0,6),rep(1,70)))

coef(elasticnet_solution)

Quand \(n\) grandit : les problèmes se compliquent

Un exemple sur simulations

n = 100
p = 20
X = matrix(rnorm(n*p,0,0.1),n,p)
beta = runif(p)
Xbeta = X%*%beta

Dans un modèle logistique

Y_logistic = rbinom(n=n,size=1,prob = exp(Xbeta)/(1+exp(Xbeta)))

Dans un modèle de Cox

Y_cox = rexp(n=n, exp(-Xbeta))

Timings

n_s = c(100,1000,10000,100000)
# times_logistic = matrix(0,5,50)
# times_cox = matrix(0,5,50)
# for (m in 1:50){
# k=1
# for (n in n_s){
#   X = matrix(rnorm(n*p,0,0.1),n,p)
#   beta = runif(p)
#   Xbeta = X%*%beta
#   Y_logistic = rbinom(n=n,size=1,prob = exp(Xbeta)/(1+exp(Xbeta)))
#   Y_cox = rexp(n=n, exp(-Xbeta))
#   # Z = rexp(n=n, 1)
#   # Y = min(Y_cox,Z)
#   # delta = (Y <= Z)
#   start.time <- Sys.time()
#   cv.glmnet(X,Y_logistic, family = "binomial" , alpha = 0.8)
#   end.time <- Sys.time()
#   times_logistic[k,m] <- end.time - start.time
#   start.time <- Sys.time()
#   cv.glmnet(X,Surv(Y_cox,rep(1,n)), family = "cox" , alpha = 0.8)
#   end.time <- Sys.time()
#   times_cox[k,m] <- end.time - start.time
#   k = k+1 
# }}
# write.table(times_logistic,"times_logistic")
# write.table(times_cox,"times_cox")

Représentation graphique

times_cox = read.table("times_cox")
times_logistic = read.table("times_logistic")
plot(n_s,rowMeans(times_cox)[1:4],type="l",ylab="temps de calcul")
lines(n_s,rowMeans(times_logistic)[1:4],col="red")

legend(x =0, y=15,lty = c(1,1),legend = c("Cox","logistic"), col = c("black","red"))

rowMeans(times_cox)[1:4]

         1          2          3          4 
 0.1096696  0.2761387  1.8763669 22.2677560 

rowMeans(times_logistic)[1:4]

         1          2          3          4 
 0.1209730  0.2581053  1.4016656 15.3722102 

Régression de Poisson et de Poisson conditionelle

Récupérer les données simulées

Elles contiennent des observations pour 10000 individus sur 5 ans avec une mesure par semaine. J’ai reproduit des données de santé de type pharmacovigilance.

On s’intéresse à la survenue d’une pathologie. On suspecte qu’un traitement augmente ce risque, proportionnellement à la dose cumulée reçue par l’individu.

Chaque individu commence (à un moment aléatoire sur les 5 ans) le traitement.

Les autres variables mesurées sont le sexe de l’individu et une variable continue (qui dans la simulation n’a pas d’influence sur le risque).

Pour chaque individu et chaque semaine soit \(10 000*5*52 = 2 600 000\) lignes de données, ont été enregistrés

• l’identifiant de l’individu (ind)

• le numéro de la semaine (week)

• la dose cumulée (cum_dose)

• le sexe (sex)

• la variable continue (var)

• le nombre d’évènements survenus dans la semaine (nevent)

• une variable d’identification de la ligne (id)

paste(c(yourpath,"data_example_poisson.Rdata"),collapse = "")

[1] "~/Dropbox/Exposes/Formation IA/Rexamples/DATA/data_example_poisson.Rdata"

Seuls 26 individus ont développé (au moins une fois) la pathologie.

data_sum = data %>%  group_by(ind) %>% summarise (neventS = sum(nevent))
sum(data_sum$neventS >0)  

[1] 26

En régression de Poisson simple

start.time <- Sys.time()
outcomes  = select(data,id,ind,week,nevent)
#head(outcomes)
colnames(outcomes) = c("rowId","stratumId","time","y")
#head(outcomes)
covariates = select(data,id,ind,cum_dose,var,sex)
covariates = gather(covariates , "covariateId","covariateValue" ,2:4)
covariates$covariateId = 1*(covariates$covariateId =="cum_dose") + 2 * ((covariates$covariateId == "var")) +  3 * ((covariates$covariateId == "sex"))
colnames(covariates) = c("rowId","stratumId","covariateId","covariateValue")
cyclopsData <- convertToCyclopsData(outcomes, covariates, modelType = "pr",addIntercept = FALSE)
fit <- fitCyclopsModel(cyclopsData, prior = createPrior("laplace"))
coefs = coef(fit)
names(coefs) = c("cum_dose","var","sex")
print(coefs)

    cum_dose          var          sex 
 -0.27576300 -53.41001632  -0.07756173 

end.time <- Sys.time()
time.taken <- end.time - start.time
print(time.taken)

Time difference of 17.34703 secs

En régression de Poisson conditionelle

start.time <- Sys.time()
data_cond = filter(data,ind %in% data_sum$ind[data_sum$neventS >0])
outcomes  = select(data_cond,id,ind,week,nevent)
colnames(outcomes) = c("rowId","stratumId","time","y")
covariates = select(data_cond,id,ind,cum_dose,var)
covariates = gather(covariates , "covariateId","covariateValue" ,3:4)
covariates$covariateId = 1*(covariates$covariateId =="cum_dose") + 2 * ((covariates$covariateId == "var"))
colnames(covariates) = c("rowId","stratumId","covariateId","covariateValue")
cyclopsData <- convertToCyclopsData(outcomes, covariates, modelType = "cpr",addIntercept = FALSE)
fit <- fitCyclopsModel(cyclopsData, prior = createPrior("laplace"))
coefs = coef(fit)
names(coefs) = c("cum_dose","sex")
print(coefs)

  cum_dose        sex 
-0.7105996  0.0000000 

end.time <- Sys.time()
time.taken <- end.time - start.time
print(time.taken)

Time difference of 0.328902 secs