Inferencia estatística con datos censurados

Traballo Fin de Grao Inferencia estatística con datos censurados Álvaro Añón Dosil Xullo, 2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Inferencia estatística con datos censurados Álvaro Añón Dosil Xullo, 2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA iii iv Traballo proposto Área de Coñecemento: Estatística e Investigación Operativa Título: Inferencia estatística con datos censurados Breve descrición do contido No estudo dos tempos de vida, que poden ser duracións dunha enfermidade, dun artículo de consumo (coches, teléfonos, ordenadores, etc.) ou calquera outro tempo entre dous eventos, a obtención dunha observación adoita requerir dun seguimento do individuo ou obxecto en cuestión. Se ese seguimento se interrompe, só coñeceremos que o tempo durou polo menos ata o momento da perda do seguimento. Nese caso diremos que o tempo foi censurado. Neste traballo estudaranse modelos de distribución propios de tempos de vida, como o exponencial ou o Weibull, xunto con estimadores dos seus parámetros con datos censurados. Tamén se estudará o estimador de Kaplan-Meier, que é o estimador non paramétrico da función de distribución con datos censurados. Expoñeranse os métodos existentes, estudaranse as súas propiedades mediante simulacións e ilustraranse con datos reais. Recomendacións Ter un coñecemento básico do programa estatístico . Outras observacións Índice Resumo viii Introdución xi 1. Funcións relevantes na análise de supervivencia 1 1.1. Funcióndesupervivencia................................ 1 1.2. Funciónderisco .................................... 3 1.3. Función de esperanza de vida residual . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4. Relacións entre as funcións . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2. Métodos de estimación paramétricos 9 2.1. Introdución ....................................... 9 2.2. Función de verosimilitude baixo censura . . . . . . . . . . . . . . . . . . . . . . . 10 2.3. Modeloexponencial................................... 11 2.4. ModeloWeibull..................................... 13 2.5. Estudodesimulación.................................. 16 2.5.1. Propiedades dos estimadores . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.5.2. Propiedades dos intervalos de conanza . . . . . . . . . . . . . . . . . . . 17 3. Métodos de estimación non paramétricos 21 3.1. Táboasdevida..................................... 21 v vi ÍNDICE 3.2. Estimador de Kaplan-Meier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.3. O paquete survival de . Aplicación aos datos pharmacoSmoking ................................... 29 3.4. Probas de bondade de axuste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.4.1. Modeloexponencial .............................. 35 3.4.2. ModeloWeibull................................. 36 A. Código de empregado 39 A.1.Figurasdaintrodución................................. 39 A.1.1.Figuras1e2................................... 39 A.1.2.Figura3..................................... 40 A.2.FigurasdoCapítulo1 ................................. 41 A.2.1.Figura1.1.................................... 41 A.2.2.Figura1.2.................................... 41 A.3.Capítulo2........................................ 42 A.3.1.Figura2.1.................................... 42 A.3.2.Figura2.2.................................... 43 A.3.3. Estudo de simulación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 A.4.Capítulo3........................................ 48 A.4.1.Exemplo3.2................................... 48 A.4.2. Estudo dos datos pharmacoSmoking . . . . . . . . . . . . . . . . . . . . . 49 Bibliografía 51 xiv INTRODUCIÓN ● ● ● A B C D E F 0 1 2 3 4 5 6 7 8 9 10 11 12 13 Tempo (meses) Pacientes Censura ●Censurado Observado Tempo de remisión de pacientes tratados de leucemia aguda Figura 3: Exemplo de censura aleatoria Debido a presenza deste tipo de observacións, precisamos construír novos estimadores que teñan en conta a censura para así coñecer mellor e estimar de forma máis precisa a distribución dos datos. Estes estimadores, os cales podemos clasicar en paramétricos e non paramétricos , serán os que nos permitirán facer inferencia e intentar sacar conclusións a partir das mostras. Ao longo deste traballo veremos como construír os distintos tipos de estimadores, discutiremos e compararemos as súas propiedades e as ilustraremos con exemplos empregando o software estatístico . Capítulo 1 Funcións relevantes na análise de supervivencia Como xa se mencionou, a variable aleatoria a estudar vai ser o tempo de supervivencia, a cal denotaremos por T . Supoñemos que T é unha variable aleatoria continua e non negativa, e como toda variable aleatoria absolutamente continua ten funcións de distribución e densidade que a caracterizan. Sen embargo, no marco da Análise de Supervivencia resulta máis interesante considerar outras funcións que, dende outro punto de vista, tamén caracterizan a mesma variable aleatoria. Ditas funcións son a función de supervivencia , a función de risco e a función de esperanza de vida residual . Neste capítulo introduciremos estas funcións e veremos que efectivamente caracterizan a T , é dicir, a partir de calquera delas podemos determinar de forma única as restantes. Antes de continuar, denimos as funcións de distribución e densidade, que serán empregadas de aquí en diante. Denición 1.1. Sexa T unha variable aleatoria absolutamente continua. Defínense as funcións de densidade, f , e distribución, F , de T como: f(t) = l´ım ∆t→0 P(t<T≤t+ ∆t) ∆t (1.1) F(t) = P(T≤t) = Zt 0 f(t)dt (1.2) 1.1. Función de supervivencia A función de supervivencia, a cal denotaremos por S(t) , é a máis empregada para analizar e describir fenómenos que involucran tempos de supervivencia. Defínese como a probabilidade de que un individuo sobreviva (experimente o evento de interés) máis alá dun tempo t . Formalmente: 1 2 1. Funcións relevantes na análise de supervivencia S(t) = P(T > t) con 0⩽t < ∞ Polo tanto, pódese relacionar coa función de distribución F e a función de densidade f da seguinte maneira: S(t) = P(T > t) = 1 −P(T⩽t)=1−F(t)∀t∈[0,∞) (1.3) S(t) = P(T > t) = Z∞ t f(t)dt ⇒f(t) = −S′(t)∀t∈[0,∞) (1.4) A representación gráca de S(t) coñécese como curva de supervivencia. Existen moitos tipos de curvas de supervivencia, pero todas comparten unha serie de propiedades: S(t) é monótona decrecente, é dicir, S(t1)≥S(t2) se t1< t2 . S(0) = 1 , polo que a probabilidade de sobrevivir nun instante inicial t= 0 é 1. l´ım t→∞ S(t)=0 , é dicir, co paso do tempo a probabilidade de supervivencia converxe a cero. Na gura 1.1 móstrase o exemplo da curva de supervivencia dunha Exponencial(1), no que se poden apreciar claramente as propiedades mencionadas, herdadas das propiedades da función de distribución. 02468 0.0 0.2 0.4 0.6 0.8 1.0 Distribución t F(t) 02468 0.0 0.2 0.4 0.6 0.8 1.0 Curva de supervivencia t S(t) Figura 1.1: Función de distribución e supervivencia dunha Exp(1) No caso de ter unha mostra completa, é dicir, sen datos censurados, a función de supervivencia nun instante t pode estimarse como a proporción de individuos sobrevivindo máis que t , o cal se coñece como función de supervivencia empírica . 1.2. Función de risco 3 ˆ Sn(t) = 1 n n X i=1 I(Ti> t) (1.5) Sen embargo, máis adiante mostraremos que en presencia de datos censurados esta estimación non é boa, de feito non está ben denida, pois para as observacións censuradas non coñecemos se sobreviviron ou non máis aló dun tempo superior ao tempo de censura. Por este motivo construiremos outros estimadores máis apropiados que teñan en conta a censura. 1.2. Función de risco A función de risco, h(t) , é a que determina o que se coñece como taxa de fallo instantáneo , que se dene como a 'tendencia' a experimentar o evento de interese (fallo, morte...) nun intervalo de tempo moi pequeno supoñendo que o individuo sobreviviu ao comezo do mesmo. Formalmente: h(t) = l´ım ∆t→0 P(t<T≤t+ ∆t|T > t) ∆t (1.6) Podemos ver que se ∆t≈0 entón h(t)∆t≈ P(t<T≤t+∆t|T > t) , que é aproximadamente a probabilidade de fallo no intervalo (t, t+∆t] . Polo tanto, podemos interpretar a función de risco como unha medida de propensión ao fallo nun instante t . Ademais, aplicando denicións básicas de probabilidade condicionada podemos relacionala coa función de densidade e supervivencia da seguinte forma: h(t) = l´ım ∆t→0 P(t<T≤t+ ∆t|T > t) ∆t = l´ım ∆t→0 P(t<T≤t+ ∆t) ∆t 1 P(T > t)=f(t) S(t) (1.7) A función de risco pode ter moitas formas, coa restricción de que ten que ser non negativa, i.e., h(t)≥0 . En moitas investigacións pode ser máis doado ter unha idea de como é o risco de fallo nunha poboación, polo que esta función pode servir de guía para propoñer un modelo paramétrico para a distribución de T . Na gura 1.2 podemos ver varios exemplos de curvas de risco. 4 1. Funcións relevantes na análise de supervivencia Risco exponencial Taxa de fallo Tempo Risco Weibull Taxa de fallo Tempo Risco Weibull Taxa de fallo Tempo Risco Lognormal Taxa de fallo Tempo Figura 1.2: Exemplos de curvas de risco correspondentes coas distribucións Exp(1), Weibull(3,1) (vermella), Weibull( 1 2 ,1) (azul) e Lognormal(2,1). Nestes exemplos podemos observar diferentes comportamentos da taxa de fallo. No primeiro caso, Exponencial( λ= 1 ), vemos que a curva de risco permanece constante, o cal quere dicir que o fallo ten a mesma probabilidade de producirse ao longo do tempo. No seguinte caso, Weibull(3,1), a función de risco é monótona crecente, o que indica que o fallo faise máis probable co paso do tempo. Polo tanto, esta curva de risco aporta unha idea de envellecemento ou desgaste na poboación de estudo. O caso contrario é o exemplo da Weibull( 1 2 ,1), no que a curva de risco é monótona decrecente, e polo tanto dá unha idea de `rexuvenecemento' na poboación. A última curva de risco da gura, procedente dunha Lognormal(2,1), mostra o caso dunha poboación na que a probabiliade de fallar crece moi rápido nos primeiros instantes de vida, alcanzando un máximo a partir do cal este risco descende progresivamente. Ademais, na gura 1.3 móstrase unha curva de risco frecuente en moitos estudos e que describe o funcionamento de moitos procesos, incluída a vida humana, denominada ás veces pola súa forma curva de bañeira . Vemos que é unha `combinación' dos casos anteriores. No inicio, o risco de fallo é alto (mortalidade infantil), pero co paso do tempo este risco redúcese sendo constante ata certo punto (mortes accidentales ou aleatorias), a partir do cal aumenta de novo (envellecemento). 1.3. Función de esperanza de vida residual 5 Figura 1.3: Curva de bañeira Para rematar, de forma análoga a como se dene a función de distribución dunha variable aleatoria en función da súa densidade, podemos denir a función de risco acumulado , H(t) , da seguinte forma: H(t) = Zt 0 h(x)dx (1.8) 1.3. Función de esperanza de vida residual A última función de interese na análise de supervivencia que imos ver neste capítulo é a función de esperanza de vida residual, mrl(t) polas súas siglas en inglés, que mide a esperanza de vida restante para un individuo en tempo t . Polo tanto: mrl(t) = E(T−t|T > t) =Z∞ 0 (x−t)f(x|x>t)dx =Z∞ 0 (x−t)f(x) S(t)I(t,∞)(x)dx =R∞ t(x−t)f(x)dx S(t) (1.9) 6 1. Funcións relevantes na análise de supervivencia Agora, empregando (1.4) e integrando por partes chegamos a que: mrl(t) = R∞ t(x−t)f(x)dx S(t)=R∞ tS(x)dx S(t) (1.10) Vendo isto, podemos deducir que a esperanza de vida, µ=mrl(0) , correspóndese coa área total baixo a curva de supervivencia. µ=E(T) = Z∞ 0 tf(t)dt =Z∞ 0 S(t)dt =mrl(0) (1.11) Traballando no contexto da análise de supervivencia ás veces vai resultar máis interesante calcular a mediana en lugar da media (esperanza de vida), xa que nos vai aportar máis información, sobre todo cando se traballa con distribucións moi nesgadas. Polo tanto, para calcular o cuantil p , bastaría con resolver a ecuación S(tp)=1−p , sendo a mediana o cuantil t0,5 . 1.4. Relacións entre as funcións Para rematar, veremos as relacións entre as funcións introducidas, que nos permitirán obter unhas a través das outras, e polo tanto caracterizarán a variable T . Se consideramos unha función de supervivencia dada, S(t) , poderiamos obter a densidade correspondente empregando (1.4). Agora, tendo en conta isto e a expresión (1.7) tamén podemos obter a función de risco da seguinte forma: h(t) = f(t) S(t)=−S′(t) S(t)=−dlog S(t) dt (1.12) Por último, a ecuación (1.10) móstranos como obter a función de esperanza de vida residual a partir da supervivencia. Feito isto, vexamos agora como obter a función de supervivencia a partir das outras e desta forma teremos cumprido o noso obxectivo. Se coñecemos a función de risco e integramos en (1.12) podemos concluir o seguinte: h(t) = −dlog S(t) dt ⇒S(t) = exp[−Zt 0 h(x)dx] = exp[−H(t)] (1.13) Esta expresión permítenos tamén relacionar o risco e a densidade: f(t) = h(t)S(t)⇒f(t) = h(t)exp[−Zt 0 h(x)dx] (1.14) 1.4. Relacións entre as funcións 7 Para nalizar, en [9] podemos ver a seguinte expresión que relaciona S(t) con mrl(t) ; S(t) = mrl(0) mrl(t)exp[−Zt 0 dx mrl(x)] (1.15) Feitos todos estes cálculos, xa estamos en condicións para ver que efectivamente as tres funcións vistas (ademais da densidade e da distribución) caracterizan unha variable aleatoria. De aquí en diante veremos diferentes técnicas e procedementos para construír estimadores destas funcións que teñan en conta a posible censura que poden ter os datos dunha mostra. Comezaremos vendo os métodos de estimación paramétricos no próximo capítulo. 8 1. Funcións relevantes na análise de supervivencia Capítulo 2 Métodos de estimación paramétricos 2.1. Introdución Os modelos paramétricos son aqueles que supoñen unha distribución coñecida para a variable aleatoria que se está a estudar. Polo tanto, son especialmente útiles cando a distribución de procedencia dos datos é coñecida de antemán, ben sexa grazas á realización doutros estudos previos ou a aplicación de test de bondande de axuste (veremos un exemplo no capítulo 3). As funcións que caracterizan estes modelos están determinadas por un conxunto nito de parámetros, sendo estes parámetros os que se deben estimar para atopar a distribución máis axeitada aos datos. Existen diferentes métodos de estimación dos parámetros, coma o método dos momentos ou o de máxima verosimilitude (o cal veremos nas seguintes seccións). No contexto da análise de supervivencia, calquera distribución con soporte os reais non negativos pode ser considerada como un modelo para o tempo de supervivencia. Sen embargo, hai unha serie de modelos que son os máis frecuentes, como os modelos exponencial e Weibull que veremos neste capítulo. En canto á estimación dos parámetros, a presenza de datos censurados vainos obrigar a adaptar os métodos de estimación para ter en conta a censura. Nun primeiro razoamento, poderiamos plantexarnos eliminar os datos censurados da mostra ou interpretalos coma fallos, sen embargo, isto produciría un sesgo cara a valores máis pequenos da variable, o cal nos levaría a subestimar a distribución teórica da que proceden as observacións. Tendo isto en conta, precisamos construír novos estimadores que teñan en conta a censura, para o cal comezamos vendo como adaptar o método de máxima verosimilitude na seguinte sección. 9 16 2. Métodos de estimación paramétricos unha ecuación que soamente depende de α , pero que non é posible resolver analiticamente. Entón, para calcular ˆα é necesario empregar métodos numéricos coma o método de NewtonRaphson ou o de Nelder-Mead. Unha vez estimado o parámetro de forma, o de escala obtémolo substituíndo na expresión (2.25). Estimados ambos parámetros, tamén podemos construír os seus respectivos intervalos de conanza empregando as propiedades asintóticas dos EMV. Desta forma, os intervalos de conanza asintóticos para os parámetros de escala e forma dunha distribución Weibull (baixo censura) de nivel (1α ) son da forma ˆα−zα/2qd V ar(ˆα),ˆα+zα/2qd V ar(ˆα) (2.27) ˆ λ−zα/2qd V ar(ˆ λ),ˆ λ+zα/2qd V ar(ˆ λ) (2.28) onde d V ar(ˆα) e d V ar(ˆ λ) se extraen dos elementos da diagonal da inversa da matriz de información de Fisher. Na seguinte sección faremos un estudo de simulación co obxectivo de ilustrar as propiedades vistas dos estimadores. 2.5. Estudo de simulación Vexamos cun estudo de simulación realizado con o cálculo de todas estas estimacións e como afecta a porcentaxe de censura. Para elo, empregamos a función mle da librería stats4 . Esta función resolve o problema de minimizar a log-verosimilitude negativa (o cal é equivalente a maximizar a log-verosimilitude). Mediante 1000 mostras censuradas (censura aleatoria) para cada tamaño da mostra n procedentes dunha Weibull(5,1) imos analizar as propiedades dos estimadores dos parámetros, dos intervalos de conanza e o efecto que ten aumentar a censura (de 25% a 50%) sobre estas estimacións (ver código no Apéndice A.3.3): 2.5.1. Propiedades dos estimadores En base ao cadro 2.1 podemos comprobar como a estimación dos parámetros converxe aos valores reais dos mesmos ao aumentar o tamaño da mostra, xa que vemos que tanto o nesgo coma o erro típico diminúen, tendendo a 0. Este comportamento reéxase tamén no erro cadrático medio (ECM), polo que dicimos que os estimadores son consistentes en media cadrática. Por outra banda, podemos apreciar que ao aumentar a censura os estimadores teñen nesgo e erro típico máis grandes (e polo tanto tamén ECM), producindo peores estimacións debido ao maior descoñecemento dos datos (sobre todo en mostras pequenas). 2.5. Estudo de simulación 17 25% censura n Nesgo( ˆα ) qd V ar(ˆα) ECM( ˆα ) Nesgo( ˆ λ ) qd V ar(ˆ λ) ECM( ˆ λ ) 10 0.97526 1.95216 4.76209 0.00972 0.07198 0.00528 30 0.31002 0.94399 0.98724 0.00147 0.04278 0.00183 300 0.02708 0.25307 0.06478 0.00060 0.01379 0.00019 3000 0.00353 0.08274 0.00686 0.00003 0.00421 0.00002 50% censura n Nesgo( ˆα ) qd V ar(ˆα) ECM( ˆα ) Nesgo( ˆ λ ) qd V ar(ˆ λ) ECM( ˆ λ ) 10 1.44182 3.33090 13.17373 0.00844 0.09681 0.00944 30 0.36805 1.21499 1.61166 0.00120 0.05040 0.00254 300 0.03730 0.32491 0.10696 0.00140 0.01703 0.00029 3000 0.00581 0.10381 0.01081 0.00021 0.00551 0.00003 Cadro 2.1: Estimacións do nesgo, erro típico e erro cadrático medio dos estimadores dos parámetros cunha censura do 25% e 50% 2.5.2. Propiedades dos intervalos de conanza Ao longo deste capítulo temos construído intervalos de conanza para os estimadores mediante o método asintótico, o cal se basea na distribución asintótica do estimador de máxima verosimilitude. Sen embargo, existen outros métodos para construír intervalos de conanza, como pode ser o método que emprega o que se coñece como perl de verosimilitude , o cal resulta máis robusto traballando con mostras pequenas. Sexa Θ = (θ1, ..., θk) o espazo de parámetros, defínese o perl de verosimilitude (prole likelihood en inglés) dun parámetro θj como PL(θj) = m´ax θ1,...,θj−1,θj+1,...,θk L(θ1, ..., θj−1, θj+1, ..., θk) (2.29) onde L fai referencia á función de verosimilitude 2.4. Polo tanto, estamos a maximizar a verosimilitude co valor do parámetro θj xado. Mediante o perl de verosimilitude podemos construír intervalos de conanza, os cales conterán os valores máis verosímiles dos parámetros, é dicir, os que teñan un perl de verosimilitude máis grande. Polo tanto, tómase como intervalo de conanza para θj de nivel (1 −α) o seguinte conxunto: {θj: 2(L(ˆ Θ) −PL(θj)< χ2 1,α} (2.30) onde ˆ Θ é o estimador de máxima verosimilitude dos parámetros e χ2 1,α é o cuantil (1 −α) da distribución chi-cadrado cun grao de liberdade. En xeral, estes intervalos serán asimétricos, a diferenza dos construídos mediante o método asintótico. A función connt de permítenos obter este tipo de intervalos, os cales incluímos xunto aos asintóticos no noso estudo de simulación. 18 2. Métodos de estimación paramétricos Intervalos por perl de verosimilitude 25% ˆαˆ λ n Cobertura Lonx. media Cobertura Lonx. media 10 0.937 6.88890 0.935 0.32393 30 0.930 3.46309 0.940 0.17199 300 0.959 1.02313 0.939 0.05264 3000 0.947 0.32164 0.954 0.01659 50% ˆαˆ λ n Cobertura Lonx. media Cobertura Lonx. media 10 0.914 9.46189 0.929 0.44844 30 0.943 4.29943 0.964 0.22312 300 0.935 1.25773 0.945 0.06584 3000 0.941 0.39513 0.947 0.02070 Cadro 2.2: Cobertura e lonxitude media dos intervalos por perl de verosimilitude construídos para α e λ cunha censura do 25% e 50% Intervalos asintóticos 25% ˆαˆ λ n Cobertura Lonx. media Cobertura Lonx. media 10 0.959 6.86132 0.917 0.27602 30 0.949 3.46437 0.937 0.16311 300 0.962 1.02327 0.938 0.05237 3000 0.947 0.32164 0.954 0.01658 50% ˆαˆ λ n Cobertura Lonx. media Cobertura Lonx. media 10 0.953 9.36289 0.891 0.36832 30 0.956 4.30045 0.967 0.20713 300 0.938 1.25800 0.941 0.06534 3000 0.943 0.39513 0.947 0.02068 Cadro 2.3: Cobertura e lonxitude media dos intervalos asintóticos construídos para α e λ cunha censura do 25% e 50% A partir dos dous tipos de intervalos calculamos a cobertura (proporción de intervalos que conteñen o valor real do parámetro) e a lonxitude media dos mesmos para os distintos tamaños da mostra. Nos cadros 2.2 e 2.3 podemos ver os resultados obtidos. En ambos casos, vemos como 2.5. Estudo de simulación 19 os resultados son coherentes coa converxencia dos estimadores, xa que a cobertura está sempre preto do 95% e a lonxitude media redúcese cando aumenta n , apreciándose moi pouca diferencia entre os dous tipos de intervalos. Ademais, da mesma forma que nas estimacións puntuais, o aumento da censura tamén afecta aos intervalos de conanza, para os cales se pode observar unha maior lonxitude media, debida á menor información que nos aportan os datos. Se ben é certo que a distribución Weibull é moi útil e empregada na análise de supervivencia, existen outros modelos alternativos que tamén teñen aplicación neste contexto. Este é o caso da distribución Lognormal, da cal se mostrou o comportamento do seu risco na gura 1.2, ou a distribución Gompertz. Estes e outros modelos adicionais poden verse máis a fondo en [10]. No seguinte capítulo introducimos os métodos de estimación non paramétricos. 20 2. Métodos de estimación paramétricos Capítulo 3 Métodos de estimación non paramétricos Cando os tempos de supervivencia non seguen unha distribución teórica ou non é posible axustarlle unha, os métodos non paramétricos son máis ecientes. Ademais, empregar estes métodos pode ser de axuda para posteriormente asignarlle unha distribución aos datos. Como xa se mencionou no capítulo 2, a función de supervivencia empírica (1.5) non é aplicable no caso de traballar con mostras censuradas. Ademais, ao igual que no caso paramétrico, se eliminamos os datos censurados ou os interpretamos como fallos estaremos subestimando a curva de supervivencia, o cal non é adecuado. Polo tanto, ao longo deste capítulo veremos os dous métodos de estimación non paramétrica da función de supervivencia censurada máis habituais, as táboas de vida e o estimador de Kaplan-Meier. 3.1. Táboas de vida A estimación da función de supervivencia mediante táboas de vida, tamén coñecidas como táboas actuariais, é unha das técnicas máis antigas para analizar e describir a supervivencia e mortalidade dunha poboación de estudo. Estas táboas son ferramentas comúns en moitos ámbitos, como as investigacións médicas, demográcas, sociolóxicas e mesmo políticas. Cando se está a traballar con mostras cun gran número de observacións este método é o máis recomendable, xa que proporciona unha curva de supervivencia clara e sintetiza nunha táboa de fácil interpretación toda a información. O primeiro paso na construción dunha táboa de vida é dividir o tempo de supervivencia nunha serie de intervalos de tempo predeterminados [tj, tj+1) , empregando como extremos 0≤ 21 22 3. Métodos de estimación non paramétricos t0< t1... < tJ≤ ∞ . Feito isto, a probabilidade de supervivencia en cada extremo tj calcúlase mediante o produto das seguintes probabilidades condicionadas: S(tj) = P(T > tj) = P(T > t0)·P(T > t1|T > t0)·... ·P(T > tj|T > tj−1) (3.1) Polo tanto, vemos que a idea que está detrás deste estimador baséase en que a probabilidade de sobrevivir a un intervalo pode ser vista como a probabilidade de sobrevivir ao intervalo anterior e entón sobrevivir un intervalo máis. Desta forma, a posible censura estase tendo en conta á hora de estimar. As diferentes columnas das táboas de vida conteñen a seguinte información: Intervalos de tempo (prexados) Ij= [tj, tj+1) . Nestes intervalos agrúpanse todas as observacións da mostra, censuradas e non censuradas. Número de datos censurados no intervalo j-ésimo, cj . As observacións censuradas inclúen aquelas ás que se lles perde o seguemento e as que aínda non fallaron ao rematar o estudo. Número de fallos ou mortes en cada intervalo, dj . Individuos vivos no comezo do intervalo, Nj . O número de individuos no comezo do primer intervalo, N0 , coincide co tamaño da mostra. Podemos obter este valor para o resto de intervalos tendo en conta que Nj=Nj−1−cj−1−dj−1 . Individuos a risco, nj . Este dato correspóndese co número de datos susceptibles de fallar ou ser censurados no intervalo Ij . Coincide con Nj , salvo que se asume que as observacións censuradas distribúense de forma uniforme ao longo do intervalo, polo tanto considérase que soamente están expostas ao risco a metade do intervalo. Desta forma nj=Nj−cj 2 . Proporción de fallos condicionada ˆqj . É a estimación da probabilidade de fallar no intervlo Ij condicionada a estar vivo ao comezo do mesmo. Vén dada por ˆqj=dj nj Proporción de superviventes condicionada ˆpj= 1 −ˆqj . Proporción de superviventes acumulada ˆ S(tj) . Trátase dun estimador da función de supervivencia no instante tj . Tense que ˆ S(t0)=1 e ˆ S(tj) = ˆpj−1ˆ S(tj−1) . No cadro 3.1 vemos un exemplo dunha táboa de vida extraída de [4]. Nela reíctense os datos de supervivencia de 2418 varóns doentes de anxina de peito e as estimacións realizadas a partir destes datos. Os tempos de supervivencia, calculados dende o momento do diagnóstico, foron agrupados en 16 intervalos. 3.1. Táboas de vida 23 [tj, tj+1)cjdjNjnjˆqjˆpjˆ S(tj) [0,1) 0 456 2418 2418 0,1886 0,8114 1,0000 [1,2) 39 226 1962 1942,5 0,1163 0,8837 0,8114 [2,3) 22 152 1697 1686 0,0902 0,9098 0,7170 [3,4) 23 171 1523 1511,5 0,1131 0,8869 0,6524 [4,5) 24 135 1329 1317 0,1025 0,8975 0,5786 [5,6) 107 125 1170 1116,5 0,1120 0,8880 0,5193 [6,7) 133 83 938 871,5 0,0952 0,9048 0,4611 [7,8) 102 74 722 671 0,1103 0,8897 0,4172 [8,9) 68 51 546 512 0,0996 0,9004 0,3712 [9,10) 64 42 427 395 0,1063 0,8937 0,3342 [10,11) 45 43 321 298,5 0,1441 0,8559 0,2987 [11,12) 53 34 233 206,5 0,1646 0,8354 0,2557 [12,13) 33 18 146 129,5 0,1390 0,8610 0,2136 [13,14) 27 9 95 81,5 0,1104 0,8896 0,1839 [14,15) 23 6 59 47,5 0,1263 0,8737 0,1636 [15,∞) 0 30 30 30 1,0000 0,0000 0,1429 Cadro 3.1: Táboa de vida de 2418 homes con anxina de peito A expresión máis empregada para a estimación da varianza deste estimador é a dada pola fórmula de Greenwood: d V ar[ˆ S(tj)] = ˆ S(tj)2 j−1 X i=1 ˆqj njˆpj (3.2) Con isto, poderiamos construír intervalos de conanza (asumindo normalidade asintótica) de forma análoga á vista para os métodos de estimación paramétricos. Polo tanto, un intervalo de conanza asintótico de nivel (1α ) para a función de supervivencia nun instante t é da forma ˆ S(t)−zα/2qd V ar[ˆ S(t)] ,ˆ S(t) + zα/2qd V ar[ˆ S(t)] (3.3) onde ˆ S(t) é o estimador obtido a partir da táboa de vida. Cando cada un dos intervalos nos que está dividido o tempo de supervivencia contén unha única observación temos o caso límite das táboas de vida (a lonxitude dos intervalos tende a 0). Na seguinte sección estudamos este caso, que é de especial interese xa que dá lugar ao coñecido como estimador de Kaplan-Meier. 24 3. Métodos de estimación non paramétricos 3.2. Estimador de Kaplan-Meier O estimador de Kaplan-Meier, tamén coñecido como o estimador produto-límite , foi proposto en 1958 [8] polos matemáticos Edward L. Kaplan e Paul Meier e trátase do estimador non paramétrico máis habitual para estimar a función de supervivencia. Antes de introducilo formalmente, vexamos cun exemplo sinxelo a idea que hai detrás deste estimador, inspirándonos no exemplo que os propios Kaplan e Meier empregan no seu artículo. Exemplo 3.1. Supoñamos que 10 pacientes se unen a un ensaio clínico a comezos do ano 2018; durante ese ano 6 dos pacientes falecen e 4 sobreviven. A nais de ano 20 pacientes adicionais se incorporan ao estudo. En 2019, 3 dos doentes que comezaron en 2018 e 15 dos que entraron máis tarde falecen, deixando polo tanto 1 e 5 superviventes no momento de rematar o estudo. Vexamos como estimar S(2) . Unha posible estimación sería empregando a mostra reducida, S∗(2) = 1 10 = 0,1 , que ignora os 20 pacientes incorporados en 2019 que soamente foron seguidos durante un ano. Sen embargo, Kaplan e Meier pensaron que este grupo de pacientes tamén podían contribuír na estimación, facendo a seguinte interpretación: S(2) = P(T > 2) = P(T > 2|T > 1) ·P(T > 1) É dicir, a probabilidade de sobrevivir dous anos pode ser vista como a probabilidade de sobrevivir un ano e entón sobrevivir un ano máis. Polo tanto, vemos que a idea que está detrás deste estimador é a mesma que a do estimador das táboas de vida, xa que a única diferencia é que neste caso a lonxitude dos intervalos tende a 0 (de ahí a denominación produto-límite ), contendo cada un un único dato. Estimando en (3.4) as probabilidades mediante as proporcións mostrais obtense o estimador de Kaplan-Meier (KM) ˆ S(2) . Polo tanto, neste caso temos que ˆ S(2) = 1 4·4+5 10+20 = 0,25 ·0,3=0,075 . Desta forma, os datos censurados están sendo tidos en conta para estimar a curva de supervivencia. En vista do exemplo, este procedemento pódese xeneralizar da seguinte forma: en base a unha mostra de n individuos que contén k≤n tempos non censurados t(1) < t(2) < ... < t(k) , a estimación da supervivencia nun instante t é o produto das estimacións das probabilidades de supervivencia ata o tempo t , é dicir, ˆ S(t) = Y {i|t(i)≤t} ˆ P(T > t(i)|T > t(i−1)) = Y {i|t(i)≤t} ˆ P(T > t(i)) ˆ P(T > t(i−1)) (3.4) Para tempos anteriores a t(1) o estimador vale 1, mentres que en t(1) temos que ˆ S(t(1)) = ˆ P(T > t(1)) . Partindo de aquí e introducindo algo de notación chegaremos á expresión máis habitual do estimador de Kaplan-Meier de S(t) . Para cada fallo t(i) denimos: 3.2. Estimador de Kaplan-Meier 25 di : número de individuos que experimentan o fallo nese instante. Vale 1 salvo que haxa empates en t(i) . ni : número de individuos a risco nese instante, é dicir, aqueles que aínda non fallaron nin foron censurados en tempos anteriores. As observacións con tempos exactamente iguais a t(i) cóntanse coma individuos a risco en ni . Polo tanto, o estimador de Kaplan-Meier da función de supervivencia defínese como: ˆ S(t) = Y {i|t(i)≤t} ni−di ni =Y {i|t(i)≤t}1−di ni (3.5) Á vista da expresión (3.5) vemos que o estimador de Kaplan-Meier é unha función escalonada, constante entre dous tempos observados e continua pola dereita. Ademais, vemos que o valor do estimador non cambia nas observacións censuradas, as cales si inúen nos ni e polo tanto no tamaño das descontinuidades de salto. Ilustraremos cun exemplo simple o cálculo e representación deste estimador. Exemplo 3.2. Supoñamos un ensaio clínico no que se quere estudar o tempo de remisión de enfermos tratados contra o cancro, no cal participan 10 persoas. Ao longo do estudo seis pacientes recaen aos 3, 6.5, 6.5, 10, 12 e 15 meses; a outro pérdeselle o seguemento aos 8.4 meses; e os tres restantes seguen en remisión ao nal do estudo despois de 4, 5 e 10 meses. No cadro 3.1 móstrase o cálculo do estimador de Kaplan-Meier a partir da mostra censurada mencionada. Neste cadro distínguense 4 columnas: na primeira colócanse os tempos de remisión, nas dúas seguintes calcúlanse os ni e di e na última calcúlase o valor do estimador de Kaplan-Meier. Vemos que nas las correspondentes a tempos censurados (sinalados cun + no expoñente) omítense os cálculos xa que di= 0 e polo tanto o valor do estimador non varía neses tempos. Na gura 3.1 represéntase gracamente o estimador de Kaplan-Meier calculado. A partir do estimador de Kaplan-Meier poderíamos estimar outras medidas de interese de forma doada. Por exemplo, poderiamos estimar a esperanza de vida, ˆµ=R∞ 0ˆ S(t)dt , calculando a área baixo a curva de Kaplan-Meier que non é máis que a suma de áreas de rectángulos (no exemplo 3.2 obtemos ˆµ= 7,416 meses). Outro exemplo sería a estimación da mediana, ˆ tmed = ´ınf {t:ˆ S(t)≤0,5} , que no exemplo anterior é ˆ tmed = 10 meses. Como último exemplo, poderiamos estimar funcións relacionadas coa supervivencia, como é o caso do risco acumulado, ˆ H(t) = −ln ˆ S(t) , aínda que para esta función é máis habitual empregar directamente o estimador de Nelson-Aalen (ver [9, páx 86]). 32 3. Métodos de estimación non paramétricos > summary(datos.km) Call: survfit(formula = datos_s ~ 1, data = datos, conf.type = "log-log") time n.risk n.event survival std.err lower 95% CI upper 95% CI 0.5 125 12 0.904 0.0263 0.837 0.944 1.0 113 5 0.864 0.0307 0.790 0.913 2.0 108 6 0.816 0.0347 0.736 0.874 3.0 102 1 0.808 0.0352 0.727 0.867 4.0 101 3 0.784 0.0368 0.701 0.846 5.0 98 2 0.768 0.0378 0.684 0.833 6.0 96 1 0.760 0.0382 0.675 0.826 ... 170.0 38 2 0.288 0.0405 0.212 0.369 0 50 100 150 0.0 0.2 0.4 0.6 0.8 1.0 Tempo (dias) S(t) Figura 3.3: Estimador de Kaplan-Meier construído a partir dos datos pharmacoSmoking Como podemos ver, estamos a traballar cun caso no que a obervación máis grande está censurada, polo que o estimador de Kaplan-Meier non está ben denido a partir dese tempo (neste caso 182 días). Como xa se mencionou, isto pode dar problemas para estimar por exemplo a esperanza de vida µ . Unha alternativa pode ser calcular a esperanza de vida restrinxida ao máximo tempo supervivencia coa opción print.rmean=T . O resultado son 77.5 días para os 3.3. O paquete survival de . Aplicación aos datos pharmacoSmoking 33 pacientes que volven a fumar antes dos 182 días. > print(datos.km,print.rmean=T) Call: survfit(formula = datos_s ~ 1, data = datos, conf.type = "log-log") n events rmean* se(rmean) median 0.95LCL 0.95UCL [1,] 125 89 77.5 6.8 49 21 65 * restricted mean with upper limit = 182 Outra característica da función survfit() é a súa capacidade de estimar distintas curvas de supervivencia en función doutras covariables. Por exemplo, con este conxunto de datos poderiamos estimar dúas curvas de supervivencia distintas, dependendo do tratamento empregado. Vemos que a información das saídas é a mesma, só que calculada para cada caso. > datost.km=survfit(datos_s~grp,data=datos,conf.type='log-log') > summary(datost.km) Call: survfit(formula = datos_s ~ grp, data = datos, conf.type = "log-log") grp=combination time n.risk n.event survival std.err lower 95% CI upper 95% CI 0.5 61 4 0.934 0.0317 0.835 0.975 2.0 57 3 0.885 0.0408 0.774 0.944 4.0 54 1 0.869 0.0432 0.755 0.932 5.0 53 2 0.836 0.0474 0.717 0.908 8.0 51 2 0.803 0.0509 0.680 0.883 ... 170.0 25 1 0.393 0.0625 0.272 0.513 grp=patchOnly time n.risk n.event survival std.err lower 95% CI upper 95% CI 0.5 64 8 0.875 0.0413 0.766 0.935 1.0 56 5 0.797 0.0503 0.676 0.877 2.0 51 3 0.750 0.0541 0.625 0.839 3.0 48 1 0.734 0.0552 0.608 0.826 4.0 47 2 0.703 0.0571 0.575 0.799 ... 170.0 13 1 0.188 0.0488 0.103 0.291 34 3. Métodos de estimación non paramétricos > datost.km Call: survfit(formula = datos_s ~ grp, data = datos, conf.type = "log-log") n events median 0.95LCL 0.95UCL grp=combination 61 37 65 42 NA grp=patchOnly 64 52 23 14 49 0 50 100 150 0.0 0.2 0.4 0.6 0.8 1.0 Tempo (dias) S(t) Combinacion Parches Figura 3.4: Estimador de Kaplan-Meier diferenciando os tratamentos Desta forma, poderiamos comparar distintas curvas de supervivencia. Neste caso, vemos claramente como o tratamento consistente en empregar unha combinación de tres medicamentos é máis efectivo que os parches de nicotina. En efecto, a gura 3.4 indícanos que o tempo ata a recaída é menor nos pacientes que empregan os parches, que se corresponden coa curva vermella. Este mesmo procedemento poderiamolo empregar con outras covariables coma o xénero, a raza ou a condición de ter deixado de fumar con anterioridade, por exemplo. É por isto que este paquete tamén conta con máis ferramentas que axudan a realizar contrastes de hipóteses ou comparacións entre poboacións, mediante diferentes tests coma o log-rank. Ademais, tamén conta con funcións que permiten modelar o que se coñecen como modelos Cox, que permiten ver que factores ou covariables inúen máis na supervivencia dunha poboación (en [11] poden verse todas estas ferramentas con detalle). 3.4. Probas de bondade de axuste 35 3.4. Probas de bondade de axuste Como xa se mencionou, os estimadores non paramétricos poden axudar a axustarlle unha distribución paramétrica aos datos. Neste capítulo veremos un método gráco para comprobar a bondade de axuste dos datos aos modelos exponencial e Weibull a partir do estimador de Kaplan-Meier da función de supervivencia. 3.4.1. Modelo exponencial Como xa temos visto, a función de supervivencia exponencial é S(t) = e−λt . Polo tanto, tomando logaritmos chegamos á seguinte expresión, que non é máis que a relación da función de supervivencia coa función de risco acumulado: −log [S(t)] = H(t) = λt (3.9) Tendo isto en conta, a partir dun conxunto de tempos de supervivencia podemos representar gracamente −log [ ˆ S(ti)] fronte aos ti , onde ˆ S(ti) é o estimador de Kaplan-Meier. Desta forma, se os datos cos que estamos a traballar proceden efectivamente dun modelo exponencial, a nube de puntos resultante debería seguir unha liña recta (a cal pode axustarse por mínimos cadrados). Ademais, a pendente desta recta pode proporcionar un estimador do parámetro λ , ˆ λ . Vexamos este procedemento cun exemplo, empregando a base de datos pharmacoSmoking e as estimacións vistas na sección 3.3. Exemplo 3.3. A partir do estimador de Kaplan-Meier xa calculado extraemos a información que nos interesa, os diferentes tempos de supervivencia observados e o valor do estimador, para despois aplicarlles as transformacións correspondentes. > EstSup=datos.km$surv # valores do estimador Kaplan-Meier > TSurv=datos.km$time # diferentes tempos observados > logEstSup=-log(EstSup) Feito isto, temos que axustar o modelo de regresión linear relativo á expresión (3.9), o cal se realiza coa función lm . Este modelo vainos permitir representar a recta de regresión xunto á nube de puntos. Na gura 3.5 vemos como os datos non se axustan ben á recta calculada, polo que a hipótese de procedencia dun modelo exponencial non parece certa. De feito, no exemplo 3.4 veremos que os datos axústanse ben a unha Weibull. 36 3. Métodos de estimación non paramétricos ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●●● ●● ●●● ● ●●● ●● ●●●●●● ●●● ● 0 50 100 150 0.2 0.4 0.6 0.8 1.0 1.2 TSurv logEstSup Figura 3.5: Nube de puntos de −log [ ˆ S(ti)] fronte a ti xunto coa recta de regresión axustada 3.4.2. Modelo Weibull No caso Weibull, recordamos que a función de supervivencia é S(t) = e−(λt)α . De aquí, se tomamos dúas veces logaritmos e operamos no segundo membro obtemos: log [−log [S(t)]] = log [H(t)] = αlog λ+αlog t (3.10) Polo tanto, de forma análoga ao caso exponencial, podemos representar gracamente log [−log [ ˆ S(ti)]] fronte a log ti , obtendo unha nube de puntos lineal se a hipótese de procedencia dunha Weibull é certa. Ademais, ao axustar unha recta de regresión podemos obter estimadores dos parámetros da Weibull; se ˆ β0 e ˆ β1 representan os estimadores do intercepto e da pendente da recta, respectivamente, temos que ˆα=ˆ β1 e ˆ λ= exp ( ˆ β0 ˆ β1) . Vexamos este procedemento aplicándoo aos mesmos datos do exemplo 3.3. Exemplo 3.4. A partir dos tempos de supervivencia e os valores do estimador de Kaplan-Meier empregados no exemplo 3.3, calculamos as transformacións necesarias e plantexamos o modelo de regresión linear relativo á expresión (3.10), obtendo o seguinte resultado: > loglogEstSup=log(-log(EstSup)) > logTSurv=log(TSurv) 3.4. Probas de bondade de axuste 37 ● ● ●●●●●●●●● ● ● ●● ●●● ●● ● ● ● ●● ● ● ●● ● ● ●● ● ●● ● ●● 012345 −2.0 −1.5 −1.0 −0.5 0.0 logTSurv loglogEstSup Figura 3.6: Nube de puntos de log [−log [ ˆ S(ti)]] fronte a log ti xunto coa recta de regresión axustada Na gura 3.6 vemos como os datos se axustan bastante ben á recta calculada, o cal pode ser un indicio de que efectivamente a distribución dos datos é unha Weibull. Ademais, podemos extraer información como a pendente ( ˆ β1= 0,4385 ) e o intercepto da recta ( ˆ β0=−2,0032 ), que nos permiten estimar os parámetros de forma e escala da Weibull. > axuste=lm(loglogEstSup~logTSurv) > axuste Call: lm(formula = loglogEstSup ~ logTSurv) Coefficients: (Intercept) logTSurv -2.0032 0.4385 Desta forma, temos que ˆα=ˆ β1= 0,4385 e ˆ λ= exp ( ˆ β0 ˆ β1) = 0,0103 . Aínda así, ao coñecer que os datos se axustan ben a unha Weibull, poderiamos calcular os estimadores dos parámetros empregando o método de estimación por máxima verosimilitude visto no capítulo 2. De feito, se calculamos os estimadores de máxima verosimilitude obtemos que ˆαemv = 0,4899 e ˆ λemv = 0,0095 , valores similares aos obtidos mediante o método gráco. Para rematar, na gura 3.7 38 3. Métodos de estimación non paramétricos móstrase como ambas curvas de supervivencia estimadas teñen unha forma moi semellante ao estimador de Kaplan-Meier, sendo a curva estimada por máxima verosimilitude lixeiramente maior que a estimada polo método gráco. 0 50 100 150 0.0 0.2 0.4 0.6 0.8 1.0 Tempo (dias) S(t) Kaplan−Meier Weibull (grafico) Weibull (emv) Figura 3.7: Comparativa do estimador de Kaplan-Meier xunto ás curvas de supervivencia Weibull estimadas mediante o método gráco (vermella) e máxima verosimilitude (azul) Anexo A Código de empregado A.1. Figuras da introdución Para a realización das guras 1, 2 e 3 empregamos o paquete ggplot2 , que incorpora diversas ferramentas para a representación e edición dos grácos. library(ggplot2) A.1.1. Figuras 1 e 2. ## Función para crear gráficos de censuras tipo I y II ## graf_censura<-function(cen,t1,titulo,xlabel,ylabel){ m<-data.frame(t1,cen) censura<-as.factor(m$cen) ggplot(m)+ geom_point(aes(x=t1,y=1:length(t1),shape=censura),size=4,stroke=1)+ scale_shape_manual(values=c(20,4),name = "Censura",labels= c(quote(Censurado), quote(Observado)))+ geom_linerange(aes(xmin=0, xmax=t1, y=1:length(t1)))+ geom_vline(xintercept=max(t1),linetype='dashed')+ coord_cartesian(ylim = c(0,length(t1)+1),xlim=c(0,max(t1)+5))+ labs(x=xlabel, y=ylabel, title=titulo)+ theme_classic()+ theme(plot.title = element_text(family = "Bookman", face = "bold", 39 40 A. Código de empregado color = 1, size = 15, hjust = 0.5, vjust = 2.5, angle = 0, lineheight = 3, margin = margin(0, 0, 0, 0)))+ scale_y_continuous(breaks=c(1:length(t1)),labels=c('A','B','C', 'D','E','F'))+scale_x_continuous(breaks=seq(0,max(t1),5)) } graf_censura(c(1,1,0,1,0,1),c(10,15,30,25,30,19),'Tempo de aparición dun tumor en ratas expostas a canceríxenos','Tempo (semanas)','Ratas') graf_censura(c(1,1,0,1,0,1),c(10,15,25,25,25,19),'Tempo de aparición dun tumor en ratas expostas a canceríxenos','Tempo (semanas)','Ratas') A.1.2. Figura 3 ## Función para crear gráfico de censura tipo III ## graf_censura3<-function(cen,t1,t0,titulo,xlabel,ylabel){ m<-data.frame(t1,cen) censura<-as.factor(m$cen) ggplot(m)+ geom_point(aes(x=t0+t1,y=1:length(t1),shape=censura),size=4,stroke=1)+ scale_shape_manual(values=c(20,4),name = "Censura", labels=c(quote(Censurado), quote(Observado)))+ geom_linerange(aes(xmin=t0, xmax=t0+t1, y=1:length(t1)))+ geom_vline(xintercept=max(t0+t1),linetype='dashed')+ coord_cartesian(ylim = c(0,length(t1)+1),xlim=c(0,max(t0+t1)+5))+ labs(x=xlabel, y=ylabel, title=titulo)+ theme_classic()+ theme(plot.title = element_text(family = "Bookman", face = "bold", color = 1, size = 15, hjust = 0.5, vjust = 2.5, angle = 0, A.2. Figuras do Capítulo 1 41 lineheight = 3, margin = margin(0, 0, 0, 0)))+ scale_y_continuous(breaks=c(1:length(t1)),labels=c('A','B','C', 'D','E','F'))+scale_x_continuous(breaks=seq(0,max(t1+t0),1)) } graf_censura3(c(1,0,1,0,1,0),c(4,4,6,8,3,3),c(2,3,4,5,9,10),'Tempo de remisión de pacientes tratados de leucemia aguda','Tempo (meses)','Pacientes') A.2. Figuras do Capítulo 1 A.2.1. Figura 1.1 ## Gráfica da función de distribución e supervivencia (exponencial) ## par(mfrow=c(1,2)) curve(pexp(x,1),0,8,main='Distribucion',xlab='t',ylab='F(t)') curve(pexp(x,1,lower.tail=F),0,8,main='Curva de supervivencia', xlab='t',ylab='S(t)') A.2.2. Figura 1.2 ## Curvas de risco ## par(mfrow=c(2,2)) # creamos a ventá gráfica # Risco exponencial(1) curve(dexp(x)/pexp(x,lower.tail=F),0,12,main='Risco exponencial',xlab='' ,ylab='',xaxt='n',yaxt='n') title(ylab='Taxa de fallo',line=1,cex.lab=1) title(xlab='Tempo',line=1,cex.lab=1) # Risco Weibull(3,1) 3*x^2=dweibull(x,3,1)/pweibull(x,3,1,lower.tail=F) # para evitar problemas de precisión. curve(3*x^2,0,12,col='red',main='Risco Weibull',xlab='',ylab='' ,xaxt='n',yaxt='n') 48 A. Código de empregado # Resultados por perfil cob_alfa=cob_alfa/s # cobertura para alfa cob_lambda=cob_lambda/s # cobertura para lambda mean(lon_alfa_pv) # lonxitude media para alfa mean(lon_lambda_pv) # lonxitude media para lambda # Resultados asintoticos cob_alfa_as=cob_alfa_as/s cob_lambda_as=cob_lambda_as/s mean(lon_alfa_as) mean(lon_lambda_as) A.4. Capítulo 3 A.4.1. Exemplo 3.2 ## Exemplo estimador de Kaplan-Meier ## library(survival) # libraría necesaria para facer análise de supervivencia # Introducimos os datos # t=c(3,4,5.7,6.5,6.5,8.4,10,10,12,15) # tempos de supervivencia c=c(1,0,0,1,1,0,1,0,1,1) # indicador de censura datos=Surv(t,c) # Estimador KM con intervalos de confianza por defecto datos.km=survfit(datos~1,type='kaplan-meier') # Estimador KM empregando a transformación log-log nos i.c. datosl.km=survfit(datos~1,type='kaplan-meier',conf.type='log-log') summary(datos.km) summary(datosl.km) # Representación sen rexións de confianza plot(datos.km,xlab='Meses',ylab='S(t)',conf.int=F) # Representación cos intervalos de confianza calculados A.4. Capítulo 3 49 par(mfrow=c(1,2)) plot(datos.km,xlab='Meses',ylab='S(t)') plot(datosl.km,xlab='Meses',ylab='S(t)') A.4.2. Estudo dos datos pharmacoSmoking ## Aplicación a datos pharmacoSmoking ## library(survival) library(asaur) library(stats4) # Preparamos os datos data(pharmacoSmoking) datos=pharmacoSmoking attach(datos) ttr[ttr==0]=0.5 datos$ttr=ttr datos=datos[,c(2,3,4)] head(datos) # Creamos obxecto de supervivencia datos_s=Surv(datos$ttr,datos$relapse) # Estimador de Kaplan-Meier datos.km=survfit(datos_s~1,conf.type='log-log') summary(datos.km) plot(datos.km,xlab='Tempo (dias)',ylab='S(t)') print(datos.km,print.rmean=TRUE) # Estimador de Kaplan-Meier por tratamento datost.km=survfit(datos_s~grp,data=datos,conf.type='log-log') summary(datost.km) plot(datost.km,xlab='Tempo (dias)',ylab='S(t)',col=c(1,2)) legend(x = "topright", legend = c("Combinacion", "Parches"), fill = c(1, 2)) # Bondade de axuste exponencial EstSup=datos.km$surv # valores do estimador Kaplan-Meier 50 A. Código de empregado TSurv=datos.km$time # diferentes tempos observados logEstSup=-log(EstSup) axuste_exp=lm(logEstSup~0+TSurv) # modelo de regresion con intercepto=0 axuste_exp plot(logEstSup~TSurv,pch=16) abline(axuste_exp) # Bondade de axuste Weibull loglogEstSup=log(logEstSup) logTSurv=log(TSurv) axuste_weibull=lm(loglogEstSup~logTSurv) # modelo de regresion axuste_weibull plot(loglogEstSup~logTSurv,pch=16) abline(axuste_weibull) # Calculo de estimadores por maxima verosimilitude s1=subset(datos,relapse=='1') s0=subset(datos,relapse=='0') lvn=function(a,l){ -sum(dweibull(s1$ttr,a,1/l, log=TRUE))-sum(pweibull(s0$ttr,a,1/l,log=TRUE,lower=FALSE))} emv=mle(lvn,list(a=0.4385,l=0.0103),method='Nelder-Mead') # Figura 3.7 plot(datos.km,xlab='Tempo (dias)',ylab='S(t)',conf.int=FALSE) curve(pweibull(x,0.4385,1/0.0103,lower.tail=FALSE),0,182,col='red',add=TRUE) curve(pweibull(x,coef(emv)[1],1/coef(emv)[2],lower.tail=FALSE),0,182, col='blue',add=TRUE) legend(x='topright',legend=c('Kaplan-Meier','Weibull (grafico)','Weibull (emv)'), lty=c(1,1,1),col=c('black','red','blue')) Bibliografía [1] D. J. Bartholomew, A problem in life testing , Journal of the American Statistical Association 52 (1957), no. 279, 350355. [2] Byron W. Brown, Myles Hollander, and Ramesh M. Korwar, Nonparametric tests of independence for censored data with application to heart transplant studies , 1973. [3] B. Efron, The two-sample problem with censored data , Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability 4 (1967), 831852. [4] Edmund A. Gehan, Estimating survival functions from the life table , Journal of Chronic Diseases 21 (1969), no. 9, 629644. [5] Richard Gill, Censoring and stochastic integrals , vol. 125, 01 1980. [6] Myles Hollander and Frank Proschan, Testing to determine the underlying distribution using randomly censored data , Biometrics 35 (1979), no. 2, 393401. [7] Ricardo Vélez Ibarrola and Alfonso García Pérez, Principios de inferencia estadística , 1 ª ed., 2 ª reimp ed., Universidad Nacional de Educación a Distancia, Madrid, 1997. [8] E. L. Kaplan and Paul Meier, Nonparametric estimation from incomplete observations , Journal of the American Statistical Association 53 (1958), no. 282, 457481. [9] John P. Klein and Melvin L. Moeschberger, Survival analysis: techniques for censores and truncateddata , Springer, New York, 1997. [10] Elisa T. Lee, Statistical methods for survival data analysis , 2nd ed ed., Wiley, New York, 1992. [11] Dirk Foster Moore, Applied survival analysis using r , Springer, Switzerland, 2016. [12] Terry M. Therneau, Survival analysis [r package survival version 3.3-1] , -03-03 2022. [13] W. Weibull, A statistical distribution of wide applicability , Journal of Applied Mechanics 18 , 293297. 51