Le document a été compilé avec le magique kernel jupyter nbstata programmé par Tim Huegerich. Comparé à l’ancienne solution statamarkdown qui exécutait Stata en mode batch, il n’y a pas photo au niveau du runtime, avec pratiquement aucune latence une fois le noyau chargé (attendre une quinzaine de seconde à la première compilation).
Le document n’a pas été compilé en PDF, seule cette version html est disponible.
Ouverture de la base
webuseset"https://raw.githubusercontent.com//mthevenin/analyse_duree/master/bases/"webuse"transplantation_m", clearwebusesetlist id year age surgery stime died in 1/10
Duree pour differents quantiles de la fonction de survie
Definition des bornes Stata-ltable
S(t)=0.90: t= .
S(t)=0.75: t= 7.623
S(t)=0.50: t= 74.729
S(t)=0.25: t= 849.325
S(t)=0.10: t= .
Avec la définition des bornes des intervalles de Sas
qlt, sas
Duree pour differents quantiles de la fonction de survie
Definition des bornes Sas-lifetest
S(t)=0.90: t= 13.977
S(t)=0.75: t= 37.623
S(t)=0.50: t= 104.729
S(t)=0.25: t= 906.993
S(t)=0.10: t= .
16.1.2 Méthode Kaplan-Meier
Mode analyse des durées: stset
Les données doivent être mises en mode analyse de durée avec la commande stset (help stset).
A minima la commande stset entre la variable de durée en argument principal et la variable de censure/évènement avec failure(nom_var) en option.
stset stime, f(died)list id stime died _st _d _t _t0 in 1/10
Survival-time data settings
Failure event: died!=0 & died<.
Observed time interval: (0, stime]
Exit on or before: failure
--------------------------------------------------------------------------
103 total observations
0 exclusions
--------------------------------------------------------------------------
103 observations remaining, representing
75 failures in single-record/single-failure data
31,938 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 1,799
+-----------------------------------------+
| id stime died _st _d _t _t0 |
|-----------------------------------------|
1. | 15 1 1 1 1 1 0 |
2. | 43 2 1 1 1 2 0 |
3. | 61 2 1 1 1 2 0 |
4. | 75 2 1 1 1 2 0 |
5. | 6 3 1 1 1 3 0 |
|-----------------------------------------|
6. | 42 3 1 1 1 3 0 |
7. | 54 3 1 1 1 3 0 |
8. | 38 5 1 1 1 5 0 |
9. | 85 5 1 1 1 5 0 |
10. | 2 6 1 1 1 6 0 |
+-----------------------------------------+
Estimation de la fonction de survie
La commande sts list permet d’afficher le tableau des estimations à chaque moment d’évènement
Note
Le tableau des estimateurs n’est reporté que pour le format htl du support
On récupére les quantiles de la durée avec stci. Par défaut la commande affiche la durée médiane. On peut choisir le quantile avec l’arguement p(valeur) en option
stcistci, p(75)stci, p(25)
Failure _d: died
Analysis time _t: stime
| Number of
| subjects 50% Std. err. [95% conf. interval]
-------------+-------------------------------------------------------------
Total | 103 100 38.64425 69 219
Failure _d: died
Analysis time _t: stime
| Number of
| subjects 75% Std. err. [95% conf. interval]
-------------+-------------------------------------------------------------
Total | 103 979 144.4728 340 .
Failure _d: died
Analysis time _t: stime
| Number of
| subjects 25% Std. err. [95% conf. interval]
-------------+-------------------------------------------------------------
Total | 103 36 9.033439 16 51
Pour afficher le graphique des estimateurs KM, on utilise la commande sts graph. On ajoute ici l’affichage des intervalles de confiance avec l’option ci
stsgraph, ci
Failure _d: died
Analysis time _t: stime
Les commandes sts list, stci et sts graph permettent avec l’option by(nom_variable) d’obtenir des comparaisons entre groupes.
16.1.3 Comparaison des fonctions de survie (KM)
On va comparer les fonctions de survie pour la variable surgery.
Graphique
On ajoute l’option by(surgery) à la commande sts graph
stsgraph, by(surgery) ci ci1opts(fc(stc1%40)) ci2opts(fc(stc2%40)) legend(pos(6))
Failure _d: died
Analysis time _t: stime
Tests du logrank
Fonction sts test. On affichera ici plusieurs variantes du test.
localtest`""l""w""tw""p""'foreach test2 of local test {sts test surgery, `test2'}
Failure _d: died
Analysis time _t: stime
Equality of survivor functions
Log-rank test
| Observed Expected
surgery | events events
--------+-------------------------
0 | 69 60.34
1 | 6 14.66
--------+-------------------------
Total | 75 75.00
chi2(1) = 6.59
Pr>chi2 = 0.0103
Failure _d: died
Analysis time _t: stime
Equality of survivor functions
Wilcoxon–Breslow–Gehan test
| Observed Expected Sum of
surgery | events events ranks
--------+--------------------------------------
0 | 69 60.34 623
1 | 6 14.66 -623
--------+--------------------------------------
Total | 75 75.00 0
chi2(1) = 8.99
Pr>chi2 = 0.0027
Failure _d: died
Analysis time _t: stime
Equality of survivor functions
Tarone–Ware test
| Observed Expected Sum of
surgery | events events ranks
--------+--------------------------------------
0 | 69 60.34 73.105398
1 | 6 14.66 -73.105398
--------+--------------------------------------
Total | 75 75.00 0
chi2(1) = 8.46
Pr>chi2 = 0.0036
Failure _d: died
Analysis time _t: stime
Equality of survivor functions
Peto–Peto–Prentice test
| Observed Expected Sum of
surgery | events events ranks
--------+--------------------------------------
0 | 69 60.34 6.0505875
1 | 6 14.66 -6.0505875
--------+--------------------------------------
Total | 75 75.00 0
chi2(1) = 8.66
Pr>chi2 = 0.0032
Comparaison des rmst
Installation de la commande strmst2
ssc install strmst2
Comparaison deux à deux (c’est pas plus mal comme cela). Je conseille de mettre les variables sous forme d’indicatrice.
Les labels de la variable ne sont jamais reportés. Ils sont remplacés par les modalités arm0 et arm1.
On peut changer la valeurs maximale de la durée avec l’option tau(valeur).
arm1 = Opération
arm0 = pas d’opération
strmst2 surgery
Number of observations for analysis = 103
The truncation time, tau, was not specified. Thus, the default tau (the minimum
> of the
largest observed event time within each group), 995.000, is used.
Restricted Mean Survival Time (RMST) by arm
-----------------------------------------------------------
Group | Estimate Std. Err. [95% Conf. Interval]
---------+-------------------------------------------------
arm 1 | 734.758 104.187 530.554 938.961
arm 0 | 310.168 42.481 226.907 393.429
-----------------------------------------------------------
Between-group contrast (arm 1 versus arm 0)
------------------------------------------------------------------------
Contrast | Estimate [95% Conf. Interval] P>|z|
---------------------+--------------------------------------------------
RMST (arm 1 - arm 0) | 424.590 204.064 645.115 0.000
RMST (arm 1 / arm 0) | 2.369 1.610 3.486 0.000
------------------------------------------------------------------------
16.2 Modèles à risques proportionnels
16.2.1 Modèle de Cox
16.2.1.1 Estimation du modèle
Commande stcox
Avec la correction d’Efron (conseillé)
HR:
stcoxyear age surgery, nolog noshow efron
Cox regression with Efron method for ties
No. of subjects = 103 Number of obs = 103
No. of failures = 75
Time at risk = 31,938
LR chi2(3) = 17.63
Log likelihood = -289.30639 Prob > chi2 = 0.0005
------------------------------------------------------------------------------
_t | Haz. ratio Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
year | 0.887 0.060 -1.78 0.076 0.778 1.012
age | 1.030 0.014 2.19 0.029 1.003 1.058
surgery | 0.373 0.163 -2.26 0.024 0.158 0.876
------------------------------------------------------------------------------
log(HR):
stcoxyear age surgery, nolog noshow efron nohr
Cox regression with Efron method for ties
No. of subjects = 103 Number of obs = 103
No. of failures = 75
Time at risk = 31,938
LR chi2(3) = 17.63
Log likelihood = -289.30639 Prob > chi2 = 0.0005
------------------------------------------------------------------------------
_t | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
year | -0.120 0.067 -1.78 0.076 -0.252 0.012
age | 0.030 0.014 2.19 0.029 0.003 0.056
surgery | -0.987 0.436 -2.26 0.024 -1.842 -0.132
------------------------------------------------------------------------------
16.2.1.2 Test de l’hypothèse PH
Test Grambsch-Therneau sur les résidus de Schoenfeld
Le test OLS est exécuté avec la commande estat phtest. Je conseille d’ajouté l’option detail pour récupéré les test à un degré de liberté, le test omnibus n’étant pas fiable.
Par défaut \(f(t)=t\).
* f(t)=t - par défaut estat phtest, detail* f(t)= 1-km - solution par défaut de Restat phtest, detail km
Test of proportional-hazards assumption
Time function: Analysis time
--------------------------------------------------------
| rho chi2 df Prob>chi2
-------------+------------------------------------------
year | 0.10162 0.80 1 0.3720
age | 0.12937 1.61 1 0.2043
surgery | 0.29664 5.54 1 0.0186
-------------+------------------------------------------
Global test | 8.76 3 0.0327
--------------------------------------------------------
Test of proportional-hazards assumption
Time function: 1 - Kaplan–Meier estimate
--------------------------------------------------------
| rho chi2 df Prob>chi2
-------------+------------------------------------------
year | 0.15920 1.96 1 0.1620
age | 0.10907 1.15 1 0.2845
surgery | 0.25096 3.96 1 0.0465
-------------+------------------------------------------
Global test | 7.99 3 0.0462
--------------------------------------------------------
Intéraction avec une fonction de la durée
Avec \(f(t)=t\)
L’intéraction s’ajoute directement en option en indiquant la ou les variables concernées (une seule de préférence) avec l’option tvc(nom_variable) et la transformation de la durée (variable t) avec l’option texp(expression). Ici on a choisi \(ft(t)\) soit l’expression dans l’option texp(_t)
stcoxyear age surgery, nolog noshow efron nohr tvc(surgery) texp(_t)
Cox regression with Efron method for ties
No. of subjects = 103 Number of obs = 103
No. of failures = 75
Time at risk = 31,938
LR chi2(4) = 21.58
Log likelihood = -287.32903 Prob > chi2 = 0.0002
------------------------------------------------------------------------------
_t | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
main |
year | -0.123 0.067 -1.84 0.066 -0.254 0.008
age | 0.029 0.013 2.15 0.032 0.003 0.055
surgery | -1.755 0.674 -2.60 0.009 -3.077 -0.433
-------------+----------------------------------------------------------------
tvc |
surgery | 0.002 0.001 2.02 0.043 0.000 0.004
------------------------------------------------------------------------------
Note: Variables in tvc equation interacted with _t.
16.2.1.3 Introduction d’une variable dynamique (binaire)
Transformation de la base en format long aux temps d’évènement
Etape 1
stset stime, f(died) id(id)
Survival-time data settings
ID variable: id
Failure event: died!=0 & died<.
Observed time interval: (stime[_n-1], stime]
Exit on or before: failure
--------------------------------------------------------------------------
103 total observations
0 exclusions
--------------------------------------------------------------------------
103 observations remaining, representing
103 subjects
75 failures in single-failure-per-subject data
31,938 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 1,799
Etape 2
stsplit, at(failure)stset stime, f(died) id(id)
(62 failure times)
(3,470 observations (episodes) created)
Survival-time data settings
ID variable: id
Failure event: died!=0 & died<.
Observed time interval: (stime[_n-1], stime]
Exit on or before: failure
--------------------------------------------------------------------------
3,573 total observations
0 exclusions
--------------------------------------------------------------------------
3,573 observations remaining, representing
103 subjects
75 failures in single-failure-per-subject data
31,938 total analysis time at risk and under observation
At risk from t = 0
Earliest observed entry t = 0
Last observed exit t = 1,799
Etape 3
gen tvc = transplant==1 & wait<=_tsort id _tlist id transplant wait tvc _d _t _t0 if id==10 , noobs
gen t2=t^2gen t3=t^3quilogite t , nologestaticquilogite t t2 , nologestaticquilogite t2 t3 , nologestatic
Akaike's information criterion and Bayesian information criterion
-----------------------------------------------------------------------------
Model | N ll(null) ll(model) df AIC BIC
-------------+---------------------------------------------------------------
. | 1,127 -275.6841 -250.2606 2 504.5212 514.5758
-----------------------------------------------------------------------------
Note: BIC uses N = number of observations. See [R] IC note.
Akaike's information criterion and Bayesian information criterion
-----------------------------------------------------------------------------
Model | N ll(null) ll(model) df AIC BIC
-------------+---------------------------------------------------------------
. | 1,127 -275.6841 -243.0576 3 492.1152 507.1972
-----------------------------------------------------------------------------
Note: BIC uses N = number of observations. See [R] IC note.
Akaike's information criterion and Bayesian information criterion
-----------------------------------------------------------------------------
Model | N ll(null) ll(model) df AIC BIC
-------------+---------------------------------------------------------------
. | 1,127 -275.6841 -254.3782 3 514.7563 529.8383
-----------------------------------------------------------------------------
Note: BIC uses N = number of observations. See [R] IC note.
Estimation du modèle
quisumyearquigen year2 = year - `r(mean)'quisum agequigen age2 = age - `r(mean)'logite t t2 t3 year2 age2 surgery, nologor
Le risque d’intérêt est compet=1, le risque concurrent est compet=2. La commande externe stcompet permet d’estimer les CIF mais ne propose pas d’output des estimateurs 1. On utilise alors la commande stcomlist en précisant le ou les risques concurrent (ici compet=2)
Rappel: Le modèle de Cox cause-specific s’excute facilement. Il suffit de passer la ou les causes concurrente en censure à droite. De même le modèle Fine-Gray (critiqué) est estimable avec la commande usine stcrreg
Attention le modèle modèle multinomial n’est pas directement relié aux CIF. Il permet néanmoins d’estimer des probabilités conditionnelles qui tiennent pleinement compte de concurrence entre les différentes issues. Pour la variable de durée on utilise la variable mois
qui expand moisquibysort id: gen t=_nquigen t2=t*tquisumyearquigen year2 = year - `r(mean)'quisum agequigen age2 = age - `r(mean)'gene = competreplacee=0 if t<moismlogite t t2 year2 age2 surgery, rrr noheader
on doit cependant l’installer pour utiliser la commande suivante↩︎
Code source
---jupyter: nbstata---# Stata::: callout-note* Le document a été compilé avec le magique kernel jupyter [**nbstata**](https://github.com/hugetim/nbstata) programmé par **Tim Huegerich**. Comparé à l'ancienne solution **statamarkdown** qui exécutait Stata en mode batch, il n'y a pas photo au niveau du runtime, avec pratiquement aucune latence une fois le noyau chargé (attendre une quinzaine de seconde à la première compilation).* Le document n'a pas été compilé en PDF, seule cette version html est disponible.:::Ouverture de la base```{stata}*| message : false*| warning : falsewebuse set "https://raw.githubusercontent.com//mthevenin/analyse_duree/master/bases/"webuse "transplantation_m", clearwebuse setlist id year age surgery stime died in 1/10```## Analyse non paramétrique### Méthode actuarielleContrairement à la formation, l'estimation sera faite sur des intervalles de 30 jours```{stata}ltable stime died, interval(30) graph ci```**Récupération des quartiles de la durée** Installation de la commande `qlt````{stata}*| eval : falsenet install qlt, from("https://raw.githubusercontent.com/mthevenin/analyse_duree/master/ado/qlt/") replace* help qlt``````{stata}qui ltable stime died, interval(30) saving(base, replace)use base, clear``````{stata}qlt```Avec la définition des bornes des intervalles de Sas```{stata}qlt, sas```### Méthode Kaplan-Meier**Mode analyse des durées: stset**Les données doivent être mises en mode analyse de durée avec la commande `stset` (`help stset`).```{stata}*| include: false*| echo : false*| message : false*| warning : falsewebuse set "https://raw.githubusercontent.com//mthevenin/analyse_duree/master/bases/"webuse "transplantation_m", clearwebuse set```A minima la commande `stset` entre la variable de durée en argument principal et la variable de censure/évènement avec `failure(nom_var)` en option. ```{stata}stset stime, f(died)list id stime died _st _d _t _t0 in 1/10```**Estimation de la fonction de survie** La commande `sts list` permet d'afficher le tableau des estimations à chaque moment d'évènement::: callout-noteLe tableau des estimateurs n'est reporté que pour le format htl du support:::::: {.content-visible when-format="html"}```{stata}sts list```:::On récupére les quantiles de la durée avec `stci`. Par défaut la commande affiche la durée médiane. On peut choisir le quantile avec l'arguement p(valeur) en option```{stata}stcistci, p(75)stci, p(25)```Pour afficher le graphique des estimateurs KM, on utilise la commande **`sts graph`**. On ajoute ici l'affichage des intervalles de confiance avec l'option `ci````{stata}sts graph, ci```Les commandes `sts list`, `stci` et `sts graph` permettent avec l'option `by(nom_variable)` d'obtenir des comparaisons entre groupes.### Comparaison des fonctions de survie (KM)On va comparer les fonctions de survie pour la variable *surgery*. **Graphique**On ajoute l'option `by(surgery)` à la commande `sts graph````{stata}sts graph, by(surgery) ci ci1opts(fc(stc1%40)) ci2opts(fc(stc2%40)) legend(pos(6))```**Tests du logrank** Fonction **`sts test`**. On affichera ici plusieurs variantes du test.```{stata}local test `" "l" "w" "tw" "p" "'foreach test2 of local test {sts test surgery, `test2'}```*Comparaison des rmst* Installation de la commande **`strmst2`**```{stata}*| eval : falsessc install strmst2```* Comparaison deux à deux (c'est pas plus mal comme cela). Je conseille de mettre les variables sous forme d'indicatrice.* Les labels de la variable ne sont jamais reportés. Ils sont remplacés par les modalités arm0 et arm1.* On peut changer la valeurs maximale de la durée avec l'option `tau(valeur)`.arm1 = Opération arm0 = pas d'opération ```{stata}strmst2 surgery```## Modèles à risques proportionnels### Modèle de Cox#### Estimation du modèleCommande **`stcox`**Avec la correction d'Efron (conseillé)**HR**:```{stata}stcox year age surgery, nolog noshow efron```**log(HR)**:```{stata}stcox year age surgery, nolog noshow efron nohr```#### Test de l'hypothèse PH**Test Grambsch-Therneau sur les résidus de Schoenfeld**Le test OLS est exécuté avec la commande **`estat phtest`**. Je conseille d'ajouté l'option detail pour récupéré les test à un degré de liberté, le test omnibus n'étant pas fiable.Par défaut $f(t)=t$. ```{stata}* f(t)=t - par défaut estat phtest, detail* f(t)= 1-km - solution par défaut de Restat phtest, detail km```**Intéraction avec une fonction de la durée**Avec $f(t)=t$L'intéraction s'ajoute directement en option en indiquant la ou les variables concernées (une seule de préférence) avec l'option `tvc(nom_variable)` et la transformation de la durée (variable *t*) avec l'option `texp(expression)`. Ici on a choisi $ft(t)$ soit l'expression dans l'option `texp(_t)````{stata}stcox year age surgery, nolog noshow efron nohr tvc(surgery) texp(_t)```#### Introduction d'une variable dynamique (binaire)**Transformation de la base en format long aux temps d'évènement** *Etape 1*```{stata}stset stime, f(died) id(id)```*Etape 2*```{stata}stsplit, at(failure)stset stime, f(died) id(id)```*Etape 3*```{stata}gen tvc = transplant==1 & wait<=_tsort id _tlist id transplant wait tvc _d _t _t0 if id==10 , noobs```**Estimation du modèle** ```{stata}stcox year age surgery tvc, nolog noshow efron nohr```### Modèle à durée discrèteVariable de durée = *mois*#### Mise en forme de la base```{stata}*| include: false*| echo : false*| message : false*| warning : falsewebuse set "https://raw.githubusercontent.com//mthevenin/analyse_duree/master/bases/"webuse "transplantation_m", clearwebuse set``````{stata}expand moisbysort id: gen t=_ngen e = diedreplace e=0 if t<moislist id year age surgery mois t e in 1/31```#### Paramétrisation avec durée quantitativeLes critères d'information ```{stata}gen t2=t^2gen t3=t^3qui logit e t , nolog estat icqui logit e t t2 , nolog estat icqui logit e t2 t3 , nolog estat ic```**Estimation du modèle**```{stata}qui sum yearqui gen year2 = year - `r(mean)'qui sum agequi gen age2 = age - `r(mean)'logit e t t2 t3 year2 age2 surgery, nolog or```#### Paramétrisation avec durée discrètePour l'exemple seulement, on procède à un regroupement des intervalles découpés sur les quartiles de la durée```{stata}*| include: false*| echo : false*| message : false*| warning : falsewebuse set "https://raw.githubusercontent.com//mthevenin/analyse_duree/master/bases/"webuse "transplantation_m", clearwebuse setexpand moisbysort id: gen t=_ngen e = diedreplace e=0 if t<mois``````{stata}xtile ct4=t, n(4)bysort id ct4: keep if _n==_Ntab ct4 e``````{stata}logit e i.ct4 year age surgery, nolog or```## Modèles paramétriquesJuste un exemple pour le modèle de **weibull**Commande **`streg`**Par défaut, le modèle de Weibull est exécuté sous paramétrisation PH. Pour une paramétrisation type AFT, ajouter l'option *time*.```{stata}*| include: false*| echo : false*| message : false*| warning : falsewebuse set "https://raw.githubusercontent.com//mthevenin/analyse_duree/master/bases/"webuse "transplantation_m", clearwebuse setqui stset stime, f(died)```**PH**: défaut```{stata}*qui stset stime, f(died) // penser à mettre la base en mode analyse de duréestreg year age surgery , dist(weibull) nolog noheader```**AFT**```{stata}streg year age surgery , dist(weibull) time nolog noheader```## Risques concurrents### Non paramétrique: estimation des ICInstaller les commandes `stcompet` et `stcomlist````{stata}*| eval: FALSEssc install stcompetssc install stcomlist```Le risque d'intérêt est compet=1, le risque concurrent est compet=2. La commande externe **`stcompet`** permet d'estimer les CIF mais ne propose pas d'output des estimateurs ^[on doit cependant l'installer pour utiliser la commande suivante]. On utilise alors la commande **`stcomlist`** en précisant le ou les risques concurrent (ici compet=2)```{stata}stset stime, failure(compet==1)stcomlist, compet1(2)```### Modèle multinomial à durée discrèteRappel: Le modèle de Cox *cause-specific* s'excute facilement. Il suffit de passer la ou les causes concurrente en censure à droite. De même le modèle Fine-Gray (critiqué) est estimable avec la commande usine **`stcrreg`**Attention le modèle modèle multinomial n'est pas directement relié aux CIF. Il permet néanmoins d'estimer des probabilités conditionnelles qui tiennent pleinement compte de concurrence entre les différentes issues.Pour la variable de durée on utilise la variable *mois*```{stata}qui expand moisqui bysort id: gen t=_nqui gen t2=t*tqui sum yearqui gen year2 = year - `r(mean)'qui sum agequi gen age2 = age - `r(mean)'gen e = competreplace e=0 if t<moismlogit e t t2 year2 age2 surgery, rrr noheader```
