Modelos Não Paramétricos

O objetivo em análise de sobrevivência é estimar a função de sobrevida da variável em estudo. Uma metodologia interessante para tratar esse problema é a estimação não paramétrica, pois neste caso não impomos um modelo teórico para as falhas e/ou censuras observadas. A função de sobrevivência possui algumas propriedades, tais como, continuidade à esquerda, $S(0) = 1$ e $lim_{t\rightarrow \infty} S(t) = 0$. A figura abaixo nos mostra a curva de sobrevivência de uma variável aleatória com distribuição exponencial.

Figura: Curva de sobrevivência para o modelo exponencial

Em (a) é apresentado a curva teórica, em (b) dividimos o tempo em partições iguais com amplitude de 0,5, em (c) dividimos em amplitudes de 0,1. Percebemos que se a amplitude tender a zero teremos a verdadeira curva de sobrevivência. Então poderemos definir a partição $\{0,\xi_1,\xi_2,\ldots,\xi_k\}$ um estimador para a curva de sobrevida seria uma função escada dada por

$$\hat{S}(\xi_j) = 1 -\dfrac{\sum_{i=1}^j D_i}{n}$$

sendo $D_i$ as observações que falharam no intervalo $(\xi_{i-1},\xi_{i}]$. No caso em que não temos censuras este seria o estimador da função de sobrevida, porém em análise de sobrevivência é necessário considerar as censuras no processo de estimação da curva de sobrevivência. Para isso definiremos algumas quantidades.

• $I_j = (\xi_{j-1},\xi_j]$;

• $N_j:$ O número de u.e. em risco (que não falhou nem foi censurada) no instante $\xi_{j-1}$;

• $D_j:$ O número de falhas em $I_j$;

• $W_j:$ O número de censuras em $I_j$.

A função de sobrevivência pode ser escrita como

$$S(\xi_j) = P(T^0 > \xi_j)$$

$$= P(T^0 > \xi_j \vert T^0 > \xi_{j-1})P(T^0 > \xi_{j-1} \vert T^0 > \xi_{j-2})\ldots P(T^0 > \xi_1 \vert T^0 > \xi_{0})$$

$$= \prod_{l=1}^{j} P(T^0 > \xi_l \vert T^0 > \xi_{l-1})$$

$$= \prod_{l=1}^{j} \left[1 - P(T^0 \leq \xi_l \vert T^0 > \xi_{l-1})\right]$$

$$= \prod_{l=1}^{j} [1 - q_l]$$

Para o caso em que não existem censuras $q_l$ pode ser estimado por $\dfrac{D_l}{N_l}$, ou seja o número de falhas dividido pelo número de unidades em risco em $I_l$. Neste mesmo caso verifica-se que $N_1 = n , N_2 = N_1 - D_1, \ldots, N_j = N_{j-1} - D_{j-1}$, então teremos

$$\widehat{S}(\xi_j) = \prod_{l=1}^{j} [1 - \hat{q}_l]$$

$$= \prod_{l=1}^{j} \dfrac{N_{j-1} - D_{j-1} - D_j}{N_{j-1} - D_{j-1}}$$

$$= \prod_{l=1}^{j} \dfrac{N_{j} - D_j}{N_{j-1} - D_{j-1}}$$

$$= \left(\dfrac{N_{1} - D_1}{N_{1}}\right)\left(\dfrac{N_{2} - D_2}{N_{1} - D_{1}}\right)\ldots\left(\dfrac{N_{j} - D_j}{N_{j-1} - D_{j-1}}\right)$$

$$= \dfrac{N_{j+1}}{n}$$

como já visto anteriormente, pois $N_{j+1} = n - \sum_{l=1}^{j} D_l$. Para o caso de censuras, precisamos corrigir o denominador do estimador de $q_l$, pois o número de u.e. em risco está sendo sobre estimado. Se assumirmos que as unidades foram censuradas uniformemente ao longo do intervalo $I_l$, um estimador razoável para $q_l$ seria

$$\hat{q}_l = \dfrac{D_l}{N_l - W_l/2}$$    assim, $$\widehat{S}(t) = \sum_{l:\xi_l < t}[1 - \hat{q}_l]$$

também conhecido como estimador tábua de vida. E sua variância é dada por (ver Lawless 1982)

$$\widehat{\mbox{Var}}(\widehat{S}(t)) = \widehat{S}(t)^2 \sum_{i:\xi_i<t}\dfrac{\hat{q}_i}{(N_i - W_i/2)(1-\hat{q}_i)}$$

O estimador produto-limite, também conhecido como estimador Kaplan-Meier é obtido quando a partição $\{0,\xi_1,\xi_2,\ldots,\xi_k\}$ é tal que a primeira falha ocorre em $(0,\xi_1]$ a segunda falha ocorre em $(\xi_1,\xi_2]$, ... e a $k$ésima falha ocorre em $(\xi_{k-1},\xi_k]$, assim o estimador é dado por

$$\widehat{S}(t)_{KM} = \prod_{j: t_j <t} \dfrac{n_j - d_j}{n_j}$$

sendo $n_j$ o número de u.e. em risco em $t_j$ e $d_j$ o número de u.e. que falharam em $t_j$. Neste caso, estamos observando as falhas pontualmente, ou seja, no instante que ocorreu. No estimador tábua de vida, observamos as falhas em intervalos, sendo difícil conhecer o exato instante de cada falha. A variância para o estimador Kaplan-Meier da função de sobrevivência é dado por

$$\widehat{\mbox{Var}}(\widehat{S}(t)_{KM}) = \widehat{S}(t)^2 \sum_{i:\xi_i<t}\dfrac{d_i}{n_i(n_i - d_i)}$$

Na seção mostramos algumas propriedades de martingais e uma justificativa mais formal de que $\widehat{S}(t) \stackrel{a}{\rightarrow} N(S(t),\sigma_t^2)$ para $t$ fixo.

Outro estimador conhecido para $S(t)$ é o estimador de Breslow, que utiliza o estimador de Nelson-Aalen. Definiremos algumas quantidades antes de apresentar os estimadores de Nelson-Aalen e Breslow. Definimos $A(t) = \int_{0}^{t} \alpha(x)dx$, sendo $\alpha(t)=\dfrac{f(t)}{S(t)}$ a função taxa de falhas. Assim

$$\alpha(t) = \lim_{h\rightarrow 0} \dfrac{A(t+h) - A(t)}{h} \stackrel{h\rightarrow 0}{\Longrightarrow} A(t+h) - A(t) \approx h \alpha(t) .$$

Desta forma $h\alpha(t) \approx P($ocorrer falha em$(t,t+h]\vert$   está em risco$)$ e no instante de falha $t_i$ essa quantidade pode ser estimada por $d_i/n_i$. O estimador de Nelson-Aalen para $A(t)$ é dado por

$$\widehat{A}(t) = \sum_{i:t_i<t} \dfrac{d_i}{n_i}$$

O estimador de Breslow para a função de sobrevivência é

$$\widehat{S}(t)_{B} = \exp( - \widehat{A}(t))$$

Fleming-Harrington fizeram um estudo teórico do estimador de Breslow. Usando a aproximação $\exp(-x) \approx 1 - x$ quando $x$ é pequeno, verificamos que o estimador de Breslow é assintoticamente equivalente ao estimador Kaplan-Meier.

$$\widehat{S}(t)_{B} = \exp( - \widehat{A}(t)) = \prod_{i:t_i<t} \exp(-d_i/n_i) \approx \prod_{i:t_i<t} (1 - d_i/n_i) = \widehat{S}(t)_{KM},$$

verifica-se que $d_i/n_i \rightarrow 0$ pois estamos supondo que $n_i \rightarrow \infty \forall i$.

A estimação não paramétrica de $S(t)$ é muito útil mesmo quando queremos ajustar um modelo paramétrico, uma das forma de verificar se um modelo é adequado é fazer um gráfico de $\log(\widehat{A}(t))$ versus $\log(t)$. Assim, se o modelo exponencial for adequado esse gráfico deverá ser próximo da relação $\log(A(t)) = \log(\lambda) + \log(t)$