3.2 - Estimador de Kaplan-Meier

Você está aqui

O estimador limite-produto, ou Kaplan-Meier como é usualmente chamado, é um estimador não-paramétrico para a função de confiabilidade. Ele é uma adaptação da função de confiabilidade empírica, que na ausência de censuras é definida como


$$\widehat{R}(t) = \dfrac{\mbox{número de itens em operação até o tempo~t}}{\mbox{número total de itens em teste}}~~~~~~~~~~(3.2.1)$$

Essa estimativa de R(t) é uma função escada com degraus nos tempos observados de falha. O estimador de Kaplan-Meier, na sua construção, considera tantos intervalos de tempo quantos forem o número de tempos distintos de falha. Os limites dos intervalos de tempo são os tempos distintos de falha. Vamos definir a seguir a expressão geral deste estimador, assim como foi proposto por seus autores. A forma geral é importante na medida em que permite levar em conta a censura do tipo aleatório, ou seja, aquela que ocorre antes do término do teste. Por exemplo, um item pode ser retirado do estudo por ter falhado devido a uma causa diferente da estudada. Ressaltamos, no entanto, que para os mecanismos de censura do tipo I e II que ocorrem com frequência em estudos de confiabilidade, o estimador Kaplan-Meier mantém a mesma forma da função de confiabilidade empírica dada pela expressão (3.2.1).

Suponha que existem n itens em teste e houve k tempos de falha distintos t1 < t2 < ... < tk , para 1 ≤ k ≤n. Pode ser que ocorra mais de uma falha para um mesmo tempo, o que é chamado de empate. Desta forma, vamos usar a seguinte notação:

dj : número falhas no tempo ti;
nj : número de itens em risco no tempo tj (não falhou e não foi censurado antes de tj), para j = 1, ... ,k.

O estimador de Kaplan-Meier de R(t) é definido como sendo 1, para 0 ≤ t < t1 e para t > t1 é definido pela expressão


$$\widehat{R}(t)=\left(\dfrac{n_1-d_1}{n_1}\right)\left(\dfrac{n_2-d_2}{n_2}\right)~\ldots~\left(\dfrac{n_j - d_j}{n_j}\right)~~~~~~~~~~(3.2.2)$$

sendo tj o maior tempo de falha menor que t.

Em alguns casos, podemos ter um tempo de censura coincidindo com algum tempo de falha tj, j =1, ... ,k. Neste caso, adota-se a convenção de que o tempo de censura ocorre imediatamente após o tempo de falha e portanto o item com tempo censurado deve ser considerado em risco neste instante. Esta convenção faz sentido uma vez que, se observamos uma unidade com censura em um tempo L, muito provavelmente esta unidade estará em operação por um tempo maior que L. Pode ocorrer também que o maior tempo observado na amostra seja uma censura e não um tempo de falha. Nestes casos, a estimativa (3.2.2) é definida apenas até esse tempo e não atinge o valor 0 para nenhum tempo.

A estimativa $ \widehat{R}(t) $ tem saltos somente nos tempos de falha e decresce por um fator (nj - dj)/nj imediatamente após cada tempo de falha tj, j = 1, ... ,k. Com isso, temos que


$$\widehat{R}(t_j^+)=\left\{\begin{array}{lcl}\dfrac{n_j-d_j}{n_j},~~~~~~~~~~~~~~~~~\mbox{se}~j= 1,\\\\ \dfrac{n_j - d_j}{n_j}~\widehat{R}(t_{j-1}^+),~~~~~~\mbox{se}~j = 2, \ldots, k, \end{array} \right.~~~~~~~~~~(3.2.3)$$

em que $ \widehat{R}(t_j^+) $ representa a estimativa da função de confiabilidade imediatamente após o tempo tj, j = 1,... ,k.

Exemplo 3.2.1:

Considerando os dados do Exemplo 2.1, que diz respeito ao número de ciclos de válvulas, temos para o primeiro tempo de falha: t1 = 5.625, número de itens em risco n1 = 30 e número de falhas d1 = 1.

clique aqui para efetuar o download dos dados utilizados nesse exemplo

Portanto, pela equação (3.2.3) temos


$$\widehat{R}(5.625^{+}) = \dfrac{30-1}{30} = \dfrac{29}{30} = 0,967$$

Para o segundo tempo de falha t2 = 11223, temos n2 = 29 e d2 = 1. Logo,


$$\widehat{R}(11.223^{+}) = \dfrac{29-1}{29}\widehat{R}(5.625^{+}) = 0,965 \ast 0,967 = 0,933$$

Procedendo da mesma maneira para os demais tempos de falha, obtemos as estimativas de confiabilidade mostradas na Tabela 3.2.1.

Tabela 3.2.1: Estimativa de Kaplan-Meier para os dados das válvulas.

tj dj nj $ \widehat{R}(t_j^{+}) $
5.625 1 30 0,967
11.223 1 29 0,933
12.128 1 28 0,9
13.566 1 27 0,867
14.921 1 26 0,833
16.513 1 25 0,8
22.138 1 24 0,767
26.791 1 23 0,733
27.144 1 22 0,7
27.847 1 21 0,667
28.613 1 20 0,633
31.224 1 19 0,6
36.229 1 18 0,567
38.590 1 17 0,533
39.580 1 16 0,5
40.278 1 15 0,467
41.324 1 14 0,433
44.540 1 13 0,4

O gráfico da função de confiabilidade estimada $ \widehat{R}(t) $ obtida a partir do estimador de Kaplan-Meier é apresentado na Figura 3.2.1. Note que este gráfico não atinge o valor 0, isto sempre ocorrerá quando o maior tempo observado for censurado.

Figura 3.2.1: Estimador de Kaplan-Meier para a função de confiabilidade.

Assim como a tabela de vida, o estimador de Kaplan-Meier está sujeito a variação amostral, o que torna desejável ter uma idéia da precisão destes estimadores.

Uma estimativa para a variância do estimador de Kaplan-Meier é dada pela fórmula de Greenwood


$$\widehat{Var}\left(\widehat{R}(t)\right) = \widehat{R}(t)^2\left[\dfrac{d_1}{n_1(n_1 - d_1)} + \dfrac{d_2}{n_2(n_2 - d_2)} + \ldots + \dfrac{d_j}{n_j(n_j - d_j)}\right]~~~~~~~~~~(3.2.4)$$

sendo tj o maior tempo de falha menor ou igual que t.

A partir do cálculo deste valor, um intervalo de confiança aproximado para R(t) com 95% de confiança é dado por


$$\widehat{R}(t) \pm 1,96 \sqrt{\widehat{Var}\left(\widehat{R}(t)\right)}~~~~~~~~~~(3.2.5)$$

Por exemplo, a estimativa da variância de $ \widehat{R}(13.000) $ é obtida a partir de (3.2.4) como


$$\widehat{Var}\left(\widehat{R}(13.000)\right) = (0,9)^2 \ast \left[\dfrac{1}{30(29-1)} + \dfrac{1}{29(29-1)} + \dfrac{1}{28(28-1)}\right]$$


$$= 0,81 \ast \left[\dfrac{1}{870} + \dfrac{1}{812} + \dfrac{1}{756}\right]$$


$$= 0,81 \ast 0,0037 = 0,003$$

Extraindo a raiz quadrada, temos que o desvio padrão de $ \widehat{R}(13.000) $ é dado por


$$\widehat{\sigma}\left(\widehat{R}(13.000)\right) = \sqrt{\widehat{Var}\left(\widehat{R}(13.000)\right)} = \sqrt{0,003}= 0,0548$$

Portanto, os limites do intervalo de confiança 95% para $ \widehat{R}(13.000) $ são dados por


$$LI = \widehat{Var}\left(\widehat{R}(13.000)\right) - 1,96\sqrt{\widehat{Var}\left(\widehat{R}(t)\right)} = 0,9 - 1,96 \ast 0,0548 = 0,7926$$


$$LS = \widehat{Var}\left(\widehat{R}(13.000)\right) + 1,96\sqrt{\widehat{Var}\left(\widehat{R}(t)\right)} = 0,9 + 1,96 \ast 0,0548 = 1,0074$$

Como o valor obtido para o limite superior é maior que 1, adotamos neste caso o valor 1. Portanto, o intervalo de confiança obtido para $ \widehat{R}(13.000) $ é dado por


$$IC_{95\%}\left(\widehat{R}(13.000)\right) = [0,7926;~1]$$

Procedendo da mesma maneira para todos os tempos de falha, obtemos os resultados da Tabela 3.2.2.

Tabela 3.2.2: Tabela de confiabilidades.

Tempo de falha Confiabilidade Desvio padrão Limite Inferior Limite Superior
5.625+ 0,967 0,033 0,902 1
11.223+ 0,933 0,045 0,844 1
12.128+ 0,9 0,055 0,793 1
13.566+ 0,867 0,062 0,745 0,988
14.921+ 0,833 0,068 0,7 0,967
16.513+ 0,8 0,073 0,657 0,943
22.138+ 0,767 0,077 0,615 0,918
26.791+ 0,733 0,081 0,575 0,892
27.144+ 0,7 0,084 0,536 0,864
27.847+ 0,667 0,086 0,498 0,835
28.613+ 0,633 0,088 0,461 0,806
31.224+ 0,6 0,089 0,425 0,775
36.229+ 0,567 0,09 0,389 0,744
38.590+ 0,533 0,091 0,355 0,712
39.580+ 0,5 0,091 0,321 0,679
40.278+ 0,467 0,091 0,288 0,645
41.324+ 0,433 0,09 0,256 0,611
44.540+ 0,4 0,089 0,225 0,575

Entretanto, para valores extremos de $ \widehat{R}(t) $ (próximo de 0 ou 1), este intervalo de confiança pode apresentar um limite inferior negativo ou um limite superior maior do que 1. Este problema pode ser resolvido utilizando uma transformação para R(t) dada por


$$U(t) = \ln[-\ln(R(t))]~~~~~~~~~~(3.2.6)$$

Portanto, uma estimativa para U(t) é dada por


$$\widehat{U}(t) = \ln[-\ln(\widehat{R}(t))]~~~~~~~~~~(3.2.7)$$

cuja variância é dada por


$$\widehat{Var}\left(\widehat{U}(t)\right)=\dfrac{\dfrac{d_1}{n_1(n_1-d_1)}+\cdots+\dfrac{d_j}{n_j(n_j -d_j)}}{\left[\ln\left(\dfrac{n_1-d_1}{n_1}\right)+\cdots+\ln\left(\dfrac{n_j-d_j}{n_j}\right)\right]^2} ~~~~~~~~(3.2.8)$$

Portanto, um intervalo aproximado com 95% de confiança para U(t) é obtido como


$$\widehat{U}(t) \pm 1,96~\sqrt{\widehat{Var}\left(\widehat{U}(t)\right)}~~~~~~~~~~(3.2.9)$$

Consequentemente, o correspondente intervalo de confiança aproximado com 95% para R(t) é dado por


$$\widehat{R}(t)^{\exp\left(\pm1,96~\sqrt{\widehat{Var}(\widehat{U}(t))}\right)}~~~~~~~~~~(3.2.10)$$

cujos limites estão sempre no intervalo [0,1].

Para o cálculo do intervalo 95% para R(13.000) do nosso exemplo, temos que a variância de $ \widehat{U}(13.000) $ é obtida a partir de (3.2.8) como


$$\widehat{Var}\left(\widehat{U}(13.000)\right)=\dfrac{\dfrac{1}{30(30-1)}+\dfrac{1}{29(29-1)}+\dfrac{1}{28(28-1)}}{\left[\ln\left(\dfrac{30-1}{30} \right)+\ln\left(\dfrac{29-1}{29}\right)+\ln\left(\dfrac{28-1}{28}\right)\right]^2}$$


$$= \dfrac{0,0037}{0,0111} = 0,3333$$

Logo, pela equação (3.2.10) temos que os limites do intervalo de confiança 95% para R(13.000) são dados por


$$LI = 0,9^{\exp\left(1,96 \times \sqrt{0,3333} \right)} = 0,9^{3,1004} = 0,7213$$


$$LS = 0,9^{\exp\left(-1,96 \times \sqrt{0,3333} \right)} = 0,9^{0,3225} = 0,9666$$

Resultando portanto no intervalo [0,7213; 0,9666] para R(13.000) com confiança de 95%, contrastando com o intervalo [0,7926; 1,0074] obtido pela fórmula (3.2.4).

Para entender como executar essa função do Software Action, você pode consultar o manual do usuário.

 

Justificativa da Fórmula de Greenwood

 

A ideia principal consiste na aplicação do método delta. Para uma função suave $ f $, temos que


$$ f(X) \approx f(c) + f^{\prime} (c) (X - c)$$

no qual $ f^{\prime} $ é a derivada de $ f $ e $ c $ uma constante "próxima" ao valor esperado de $ X $. Com isso, obtemos que 


$$\mathbb{E} \left[ f(X) \right] \approx f(c) + f^{\prime}(c) \left[ \mathbb{E}(X) - c \right] $$

e


$$ Var \left[f(X) \right] \approx \left[ f^{\prime} (c) \right]^2 Var(X). $$

 

Para derivarmos a fórmula tradicional de Grennwood utilizaremos $ f(x) = \ln (x) $. Ao aplicarmos o logaritmo natural no estimador de Kaplan-Meier, obtemos


$$ \ln \widehat{R}(t) = \sum_{t_j \leq t} \ln \left(1- \frac{d_j}{n_j} \right).$$

A distribuição condicional de $ d_j $ dado que temos $ n_j $ sobre risco em $ t_j $ é binomial com parâmetros $ p_j $ e $ n_j $. Desta forma, temos que $ \mathbb{E} (d_j \mid n_j) = n_j p_j $ e $ Var (d_j \mid n_j) = n_j p_j (1-p_j) $. Através do método delta com $ f(x) = \ln (x) $ e $ c = p_j $, concluímos que 


$$ \ln \widehat{R}(t) = \sum_{t_j \leq t } \left( \ln(1-p_j) - \frac{1}{1-p_j} \left(\frac{d_j}{n_j} -p_j \right) \right)$$

Confiabilidade

Sobre o Portal Action

O Portal Action é mantido pela Estatcamp - Consultoria Estatística e Qualidade, com o objetivo de disponibilizar uma ferramenta estatística em conjunto com uma fonte de informação útil aos profissionais interessados.

Facebook

CONTATO

  •  Maestro Joao Seppe, 900, São Carlos - SP | CEP 13561-180
  • Telefone: (16) 3376-2047
  • E-Mail: [email protected]