Análise de
Sobrevivência
IM - UFRJ
Professor: Dani Gamerman
1. CONCEITOS ESTATÍSTICOS EM
SOBREVIVÊNCIA
1.1 Introdução
Análise de Sobrevivência é o estudo de indivíduos
(itens observados) onde um evento bem definido (falha)
ocorre depois de algum tempo (tempo de falha).
Exemplos:
(i) tempo de falha de equipamentos industriais
(engenharia)
(ii) tempo de sobrevida de um paciente (medicina)
(iii) tempo de duração do período de desemprego ou
greve (economia)
Definindo algumas características:
(i) as variáveis de resposta são não-negativas
(ii) principalmente univariados e contínuos
(iii) presença comum de censura
Primeiramente, é importante ter uma definição precisa
de tempo de falha. Isto requer especificações sobre:
 origem do tempo de falha
 unidade de medida do tempo (calendários, tempo de
operação, milhagem, número de ciclos)
 falha (mais fácil na medicina: morte, não tão fácil na
engenharia).
Exemplo: tempo de falha de um carro
Questões a serem feitas: Quando começar a contar o
tempo? Como medir o tempo de falha? O que é falha?
Nós devemos estudar melhor todas essas
especificações.
1.2 Resultados em Sobrevivência
1) Descritiva vs Inferência Estatística
Em algumas aplicações, características descritivas
simples como média simples, função de sobrevivência e
gráficos de probabilidades são suficientes. Em outras
aplicações, intervalos de confiança ou taxa de influência
de determinadas variáveis são exigidas.
2) Censura
Um sistema pode falhar mesmo antes que todos os
itens tenham falhado em um determinado tempo. Este fato
tem determinadas razões.
Os itens são normalmente censurados à direita. Mas
podem ainda ser censurados à esquerda ou podemos
determinar um intervalo de censura(mais difícil de
analisar)
3) Paramétrico ou Não-Paramétrico
Ambos serão vistos no curso. E ainda os SemiParamétricos serão apresentados.
4) Amostras Simples vs Modelos de Regressão
Se os itens pertencem à mesma população (são
similares), então uma análise de amostras simples deve
ser utilizada. Se os itens não são da mesma população e se
suas diferenças podem ser contadas (máquinas submetidas
a pressões distintas), então estas diferenças devem ser
consideradas na Análise.
Variáveis para medir tais diferenças (pressão) são
denominadas variáveis explanatórias ou covariáveis.
Ambos serão vistos no decorrer do curso.
5) Clássica vs Bayesiana
Ambas serão rapidamente revisadas e vistas neste
curso.
Modelos Bayesianos tiveram
no passado a
desvantagem de que sua Análise através de Modelos de
Regressão requeriam muito esforço computacional . Isto,
com o tempo, foi se tornando cada vez menos importante.
Métodos bayesianos são importantes em estudos de
sobrevivência porque freqüentemente temos informações
de experiências anteriores que podem usualmente serem
combinadas com os dados e incorporadas à análise.
1.3 Sistema Reparáveis e Não-Reparáveis
Outra importante definição, embora somente depois
será vista neste curso.
Considere o experimento com um sistema reparável e
os seguintes tempos de falha cumulativos:
203, 286, 481, 873, 1177, 1438, 1852, 2091, 2295, 2632.
Vejamos alguns itens interessantes sobre se a taxa de
falha aumenta ou diminui com o tempo.
Se o sistema é tomado como não reparável, então o
tempo entre falhas é considerado.
Um exemplo de análise simples na Figura 1.1 indica o
aumento da taxa de falha.
Conforme o tempo passa, é mais provável que uma
máquina deste tipo venha a falhar.
Para sistemas reparáveis o exemplo da análise simples
na Figura 1.2 indica taxa de falha constante em relação ao
tempo.
Se os tempos entre falhas são ordenados e o sistema é
reparável, a Figura 1.3 indica que a taxa de falha é
decrescente.
Sistemas reparáveis com processo de tempo de falha
ideais não serão vistos neste curso.
H
(
t
1
)
0. 0.5 1.0 1.5 2.0 2.5 3.0
0
10 0
20 0
30 0
40 0
t1
H
(
t
2
)
0. 0.2 0.4 0.6 0.8 1.0
Figura 1.1 – Função de Risco Acumulada dos Tempos entre Falhas.
0
50
1
0
0
1
0
5
0
2
0
0
0
2
0
5
0
00
t2
Figura 1.2 – Proporção de Falhas vs Tempo de Falha.
H
(
t
3
)
0. 0.2 0.4 0.6 0.8 1.0
0
50
1
0
0
1
0
5
0
2
0
0
0
2
0
5
0
00
t3
Figura 1.3 – Proporção de Falhas vs Tempo de Falha (Ordenados)
1.4 Componentes e Sistemas Reparáveis
Neste curso, trataremos de componentes reparáveis
sem referência aos sistemas que podem conter tais
componentes. Porém, é importante afirmar que Sistemas
Reparáveis também são uma importante área de estudo.
Os sistemas mais simples são:
 Sistemas em Série: o sistema só funciona se todas as
componentes funcionarem.
 Sistemas Paralelos: o sistema só falha se todas as
componentes falharem
 Sistemas k out of n : o sistema só funciona se pelo
menos k das n componentes funcionarem(se k=1
paralelo e se k=n  série).
Esses sistemas formam grande parte dos grupos de
sistemas chamados Sistemas Coerentes.
Para entender o conceito de sistemas coerentes é útil
definir o indicador de funcionamento xi para a
componente i e a função de estrutura das n componentes
do sistema:
(x1,x2,...,xn) =
1, se o sistema funciona.
0, caso contrário.
Sistemas coerentes têm função de estrutura que
satisfaz:
(i) (1,1,..., 1)=1
(ii) (0, 0,..., 0)=0
(iii)  é não-decrescente nesses argumentos.
Outros sistemas são:
 Sistemas multi-estados: onde as componentes podem
estar em vários estados (não só funcionando ou falhando)
 Sistemas load-sharing: onde a carga do sistema é
distribuída entre as componentes que funcionam.
1.5 Distribuições Binomial e Hipergeométrica
Análises estatísticas simples são obtidas se os tempos
de falha são dicotomizados:
 funcionar até certo tempo (digamos tempo de falha):
defeituoso.
 funcionar além deste tempo: não-defeituoso.
Distribuições Binomial e Geométrica são para um
número X de itens defeituosos em uma amostra de
tamanho n e probabilidade p do item ser defeituoso.
Se uma amostra com reposição (ou que não seja de
uma população muito grande), a distribuição Binomial é
obtida como:
P(X=k) = n pk(1-p)n-k onde 0 k n
k
E(X)=np e Var(X)=np(1-p)
Quando p é desconhecido, podemos estimá-lo como:
P = X/n
e
Var(X) = np(1-p)
Para n grande, np5 e n(1-p)5, a Binomial é
aproximada pela distribuição Normal com momentos
conforme os descritos anteriormente.
Para uma aproximação ao nível de significância
100(1-)%, o intervalo de confiança para p é dado por:
(X/n - z/2 (X(n-X)/n3)1/2, X/n + z/2 (X(n-X)/n3)1/2)
onde z/2 é o quartil(1-/2) da distribuição N(0,1).
Outra aproximação para n grande e =np (p pequeno)
é a distribuição de Poisson com média .
Equivalentemente testando a independência e
constância de p.
Se as amostras são sem reposição e de uma população
finita, a distribuição Hipergeométrica é obtida como:
n N-n
P(X=k)= k K-k
k=0,1,..,K
N
K
onde N é o tamanho da população e K é a população de
itens defeituosos.
E(X)=np e Var(X)=np(1-p)N-n
N-1
para p=K/N
1.6 Processos de Poisson
Usado particularmente para sistemas reparáveis.
Assume-se primeiramente que é observada uma série
de ocorrências em linha. (As ocorrências devem ser falhas
sucessivas do sistema e a linha representa o tempo real).
Assume-se que:
(i) as falhas ocorrem em intervalos disjuntos e
independentes
(ii) ocorrência = falha
(iii) a taxa de falha é uma constante .
Então se X é o número de falhas num intervalo de
tamanho s, X tem distribuição de Poisson com média s.
Também, os tempos entre falhas são independentes e
exponencialmente distribuídos com MTFB -1.
Isto indica que a exponencial com linha base (solo) é a
distribuição para o tempo de falha.
Isto pode ser generalizado para permitir taxas de falha
não constantes.
 Processo de Poisson não-homogêneo.
2. DISTRIBUIÇÕES DE PROBABILIDADE
2.1 Introdução
Em muitas áreas de aplicação da estatística, o ponto
inicial para avaliação da variável de interesse é a
distribuição Normal. Isto pode resultar de uma
consideração pragmática pura ou da argumentação
baseada no Teoria do Limite central, que diz que se uma
variável aleatória é a soma de um grande número de
pequenos efeitos , então a distribuição é aproximadamente
Normal.
No contexto de sobrevivência, o caso da normalidade é
muito menos usado. Para que possamos entender, tempos
de vida e resistência são quantidades positivas. Do ponto
de vista do modelo, é natural começar a pensar no
processo de Poisson, idéias já discutidas na sessão 1.6,
baseado na distribuição Exponencial. Contudo esta
distribuição tem uma limitada aplicabilidade na prática,
generalizações da Exponencial como a Gamma e a
Weibull já provarão ter maior valor prático em modelos
de sobrevivência. Estas e outras distribuições de
probabilidade comumente encontradas em estudos de
sobrevivência são discutidas nas sessões 2.3 à 2.7.
Outros aspectos centrais da discussão sobre análise de
sobrevivência são as funções de sobrevivência e de risco,
e a natural ocorrência de dados censurados. Estes assuntos
são discutidos nas sessões 2.2 e 2.8.
Finalmente neste capítulo colocamos os resultados
probabilísticos em contexto de análise de dados. Contudo
métodos gerais para ajustar distribuições de
probabilidades são desenvolvidos no Capítulo 3, algumas
técnicas básicas são apresentadas na sessão 2.9 à 2.11.
2.2 Conceitos Iniciais para a Distribuição de
Sobrevivência
Chamaremos de T a variável aleatória que representará
o tempo de falha dentro do nosso estudo. Aqui a noção de
tempo é usada de maneira genérica. Ele pode ser
realmente tempo ou qualquer outra variável não negativa,
desde que haja um número qualquer de falhas ou quebras
associado a variável.
Denotamos : F (t )  Pr(T  t )
sendo a distribuição de T e denotamos :
S (t )  Pr(T  t )  1  F (t )
sendo a Função de Sobrevivência de T. Note que alguns
autores definem F(t) e S(t) por Pr(T<=t) e Pr(T>t)
respectivamente. Na prática isto não faz diferença para os
resultados que se seguem quando T é uma variável
continua, este caso será considerado a partir de agora. Nós
iremos assumir que T tem a função de densidade :
dF (t )
dS (t )

dt
dt
tal que a probabilidade da unidade falhar em um curto
espaço de tempo [t , t+t) é :
f (t ) 
Pr(t  T  t )  f (t )t
Considere a probalidade condicional do item falhar
naquele instante [t , t+t) e não ter falhado até o tempo t :
Pr(t  T  t | t ) 
f (t )t
S (t )
Podemos pensar como a probabilidade do item
iminentemente falhar em t. A função h(t) é dada por :
h(t ) 
f (t )
S (t )
esta é a função de risco, de taxa de falha, ou hazard e é
um indicador natural da propensão a falha após uma
unidade de tempo ter transcorrido. A função de taxa de
falha acumulada é dada por :
t
H (t )   h(u )du
e concluímos que :
0
S (t )  exp(H (t ))
2.1
Note que f, F, S, h e H são funções tais que o
conhecimento de uma delas nos permite o cálculo de
todas as outras.
Alguns casos típicos são discutidos aqui :
• Se h(t)= é constante então H(t)= t e S(t)=exp(-t ) é a
distribuição de sobrevivência exponencial com parâmetro
. A densidade correspondente é f(t)= exp(-t).
• Se h(t) é uma função crescente de t , então T é dito ter
uma taxa de falha crescente (IFR). Isto á apropriado
quando a unidade medida tem relação com fadiga ou
danos cumulativos.
• Se h(t) é uma função decrescente de t, então T é dito ter
um taxa de falha decrescente (DFR). Isto pode ocorrer,
por exemplo, quando o processo diminui a quantidade
produzida ao longo do tempo diminuindo o risco de falha.
Isto é comum em alguns ambientes de produção de
componentes eletrônicos.
• Outro caso comum mencionado é o “bat-tub harzard”
onde a função de taxa de falha é decrescente inicialmente
e depois torna-se crescente. Isto costuma acontecer em
linha de produção onde os componentes iniciais tem uma
qualidade melhor que os finais provocando este tipo de
oscilação na taxa de falha.
2.3 A Distribuição Exponencial
Como mencionado na sessão 1.6, a distribuição
exponencial é o ponto natural de início para uma
distribuição de sobrevivência. Relembrando temos que a
Distribuição de Sobrevivência, hazard e função de
densidade tem a seguinte forma:
S (t )  exp( t )
h(t )  
2.2
f (t )   exp( t ),
onde  é um parâmetro positivo, freqüentemente chamado
de taxa, e onde t>0. Note também que a distribuição
exponencial tem média 1/  e variância 1/ 2. A forma da
densidade é a mesma para todos os , e 1/  age como um
parâmetro de escala. Então, por exemplo, se o tempo de
sobrevivência, T, de um certo tipo de componente é
medido em minutos e ele é distribuído exponencialmente
com taxa igual , então T*=T/60 medido em horas é
distribuído exponencialmente com taxa 60 . Uma outra
formulação comum é termos a parametrização = 1/ no
lugar de . A Figura 2.1 mostra duas densidades da
distribuição exponencial com diferentes taxas.
As funções hazard correspondentes são apresentadas
na Figura 2.2. Nós iremos mostrar que a Distribuição
Exponencial é um caso especial das Famílias de
Distribuições Weibull e Gamma.
2.4 Distribuições Weibull e Gumbel
Uma variável aleatória Weibull (W. Weibull
(1939,1951)) possui a seguinte função de sobrevivência:

t
S (t )  exp(   )
 
2.3
para t>0 e onde  e  são parâmetros positivos, sendo 
um parâmetro de escala e  um parâmetro de forma. Note
que quando =1, obtemos uma Distribuição Exponencial
com parâmetro =1/ .
A função de falha (hazard) da Weibull é :
h(t )    t 1
Então temos DFR para <1, constante para =1 e IFR
para >1. Em particular, para 1<<2 a função de falha se
aproxima de uma função linear e para =2 a função é
linear; para >2 a função cresce rapidamente acima de
uma função linear. A função de taxa de falha (hazard) da
Weibull para diferentes valores dos parâmetros é mostrada
na Figura 2.3.
A função de densidade da Weibull é
f (t )  
  1
t

t
exp(   )
 
2.4
para t > 0.
A média e a variância são dadas por :
media  ( 1  1)
var iancia  (2 1  1)
onde,

( x)   u x 1 exp( u )du
0
2.5
veja , por exemplo, Abramowitz and Stegun (1972,
CAPITULO 6). Um programa em fortran para calcular a
equação 2.5 é dado em Griffiths and Hill (1985, pp. 2436), que é baseado em um programa anterior de Pike e Hill
(1966). Quando  é grande (>5), a média e a variância são
aproximadamente  e 1.642/  respectivamente. A forma
da densidade depende de . Na Figura 2.4 são mostradas
algumas funções de densidade da Weibull para diferentes
valores de .
A Distribuição Weibull é provavelmente a mais utilizada
das distribuições em análise de sobrevivência. Uma
possível explicação para isto se deve ao seu
comportamento nos extremos da distribuição, à
possibilidade de variarmos o seu formato e em particular
a possibilidade de utilizá-la como uma generalização da
Exponencial.
A Distribuição de Gumbel tem a seguinte função de
sobrevivência :
S ( x)  exp( exp[( x   ) /  ])
2.6
para    x   , onde  é o parâmetro de locação e  é o
parâmetro de escala. Esta distribuição também começa
com limite de distribuição mínimo, veja Galambos
(1978), e tem uma taxa de falha exponencial crescente.
Em alguns casos permite valores negativos com
probabilidades positivas. Mais comumente a distribuição
de Gumbel é gerada através de Log(t) quando T tem uma
distribuição Weibull. A relação entre os parâmetros da
Gumbel e da Weibull é a seguinte :
  log( )
  1/ 
A função de densidade da Gumbel é a seguinte :
f ( x)   1 exp(( x   ) /  ) S ( x)
para    x   , e tem a mesma forma para todos os
parâmetros.
Note que a média e a variância da Gumbel são  - e
(2/6)2, respectivamente, onde =0.5772 é a constante de
Euler, e a distribuição é negativamente inclinada. A
densidade e a taxa de falha (harzard) para a distribuição
de Gumbel com =0 e =1 é mostrada nas Figuras 2.5 e
2.6 respectivamente.
2.5 Distribuições Normal e Lognormal
A distribuição Normal é a distribuição mais comumente
utilizada em Estatística.
Em confiabilidade é geralmente usada como um
modelo para log T. A função de densidade da
distribuição lognormal é descrita pela equação abaixo:
f(t) 
1
  (log ²
exp

2²
2²t


E(T )  exp ² / 2
Var(T)  exp2²exp(²)  1
As funções de Sobrevivência e de Risco podem ser
escritas somente em termos de integrais.
Algumas densidades e funções de risco são plotadas na
figura 2.9 e 2.10. A função de Risco é crescente para
valores de t próximos de zero
e eventualmente
decrescente quando t   .
(a)
(b )
(d )
(c )
Figura 2.9 Funções de densidade para 4 distribuições log normal com média 1 e
(a)  2 (b)   0 (c)   1(d )   1.25
(a)
(c )
(b )
(d )
figura 2.10 Funções de risco para 4 distribuições log normal com média 1 e
(a ) ρ = 0.25 (b) ρ = 0.5 (c) ρ = 1 (d) ρ = 1.25
2.6
Distribuições
Gama
Generalizada
e
Gama
A distribuição Gama é descrita pela equação
abaixo:
λρ ρ 1
f (t) =
t exp( λt )
Γ( p)
E (T ) =
ρ
λ
e Var (T) =
ρ
σ2
Densidades são positivas (ver figura 2.11) mas
tendem para normal quando é grande. As funções de
Sobrevivência e risco podem ser escritas somente em
termos de integrais.
A função de risco é decrescente para ρ < 1 ,
constante para ρ = 1 (exponencial) e crescente para   
(ver figura 2.12).
A Gama é obtida como a distribuição do -ésimo
tempo de falha em um processo de Poisson .
(a)
(c )
(b )
(d )
(b )
figura 2.11 Funções de densidade para 4 distribuições gama com média 1 e
(a ) ρ = 0.5 (b) ρ = 1.5 (c) ρ = 2.5 (d) ρ = 5.0
(a)
(b )
(c )
(d )
figura 2.12 Funções de risco para 4 distribuições gama com média 1 e
(a ) ρ = 0.5 (b) ρ = 1.5 (c) ρ = 2.5 (d) ρ = 5.0
A distribuição gama generalizada é descrita pela equação
abaixo:
 t 
f(t) 
exp ( t 

A distribuição gama generalizada inclui os seguintes
casos especiais:
(i) Gama : 
(ii) Weibull:   
(iii) Exponencia
l:     
2.7 Distribuição Exponencial por Partes
Uma generalização da distribuição exponencial
Baseado na partição de (0, ) em I1 , I 2 ,
Temos abaixo a função de Risco:
 , t  I1  [0, t1 ]
 , t  I  [0, t ]

2
2
h(t )   2

m , t  I m  [0, t m ]
Vantagem: Pode-se aproximar qualquer função de
Risco desejada.
Desvantagem: Grande número de parâmetros (“nãoparamétrica”)
S (t )  ki exp- i (t - ti -1 ), t
f (t )  i ki exp- i (t - ti -1 )

onde ki  exp

E (T )  1 
m
1
i -1
i, i  1,...,m


  j (t j - t j -1 ), i  2,...,m e k1  1
j 1
ki (1 1  11 )
i2
P r(T  t / T  ti 1 )  exp i (t ti 1 ), t
Ii
2.8 Censura
Observações incompletas freqüentemente ocorrem nos
estudos de sobrevivência e confiabilidade.
Nos testes de confiabilidade é comum aguardar até
todos os itens falharem.
Nos estudos de sobrevivência, pacientes abandonam o
tratamento ou continuam vivos depois do final dos
estudos. Isso resulta em algumas observações
incompletas, ditas censuradas.
Tipos comuns de censura a direita:
Tipo I: Observações são acompanhadas até um tempo c
fixado inicialmente.
Tipo II: Observações são acompanhadas até obter-se um
número pré-determinado de falhas.
Tipo III: Aleatória à direita: Associado aos tempos
i = min( Tde
i , Ci )
Y
C
'
s
falha Ti existem i onde observa-se apenas
e
1 , Ti  Ci (não censurado)
Xi  
0 , Ti  Ci (censurado)
Yi
onde
é o tempo de falha observado, Ti e Ci ' s
independentes
2.9 Métodos dos Momentos para Dados
Simples: Sem Censura
Métodos informais (métodos formais serão
apresentados no próximo capítulo)
Baseado nos momentos e estimativas simples
Suponha que t1,...,tn sejam tempos de falha observados.
Exemplo 2.1: T – número de milhões de revoluções de
rolimã até a falha.
Dados: 17.88, 28.92, 33.00, 41.52, 42.12, 45.60, 48.40,
51.84, 51.96, 54.12, 55.56, 67.80, 68.64, 68.64, 68.88,
84.12, 93.12, 98.64, 105.12, 105.84, 127.92, 128.04,
173.40 (ordenados por conveniência)
Média amostral: t = 72.22 e desvio padrão amostral (s.d.):
st = 37.49
Recai segundo §2.3 que a média e o s.d. coincidem no
modelo exponencial.
Neste caso, t /st se aproxima de 2, logo o modelo
exponencial não é apropriado.
Para ajustar Weibull e lognormal, é mais fácil trabalhar
com xi = log ti.
Novamente, média amostral: x = 4.150 e s.d. amostral:
sx = 0.534.
Estes cálculos valem para a média μ –γσ e s.d. πσ/√6
~
da Gumbel, trazendo assim os momentos estimados σ =
0.416 e μ~ = 4.390.
~
Em termos dos parâmetros da Weibull, temos: α =
exp( μ~ ) = 80.64
σ =2.40, diferente do 1.
e ~η = 1/ ~
De forma similar, os parâmetros estimados da
σ = 0.534.
lognormal são μ~ = 4.150 e ~
Outra aproximação é baseada na função de
sobrevivência empírica Sˆ (t ) dada por
número _ de _ ti ' s  Yt
Sˆ t  
n
Sˆ (t ) é um estimador não paramétrico de S(t)
n Sˆ (t ) possui distribuição binomial com média S(t) e ,
segundo §1.15, um IC a 100%(1-α) para S(t) dado por



1

ˆ t 1  Sˆ t   2

 Sˆ t 1  Sˆ t  
S
ˆ
ˆ
S t   z 
 , S t   z 

2
2
n
n




1
2
De forma similar, a função de taxa de falha acumulada
empírica é dada por
Hˆ t    log Sˆ t 
1  Sˆ t 
Onde o s.d. é dado por 

ˆ
 nS t  
1
2
Os gráficos destas funções empíricas podem ser usados
para checar a adequação das hipóteses dos parâmetros.
Assuma os tempos ordenados t(1) < ... < t(n) .
i -1
Sˆ ti    1 n
• Salto de 1 / n tempo t(i).
• Realocado por 1- (i – 0.5) / n . (Outra forma possível)
Modelo Weibull: S(t) = exp {- (t / α)}
Log S(t) = - (t / α)  log{-log S(t)}=  log t –  log α
Se o modelo Weibull é apropriado, então o gráfico de
é aproximadamente uma linha reta.
logt , log  log Sˆ t

i 

   
i
Inicialmente os parâmetros estimados serão obtidos a
partir de :
- log α = intercepto
 = coeficiente angular
 log t   
 , – função
Modelo lognormal: S t   1  



de distribuição de N(0,1).
 11  S t  
log t




Se o modelo lognormal for apropriado, então o gráfico
de (log t(i), -1{(i-0.5)/n}) será uma linha reta.
Inicialmente os parâmetros estimados serão dados por :
-μ / σ= intercepto
1 / σ = coeficiente angular
Estes gráficos são dados nas figuras 2.13 e 2.14 do
exemplo 2.1.
Inicialmente os parâmetros estimados são (entre
parênteses por momentos):
Weibull:  =2.3(2.4) e α = 77.3(80.6)
Lognormal: μ = 4.2(4.15) e σ = 0.56(0.53)
Diferentemente do baseado pelo método dos
momentos, o gráfico das probabilidades pode ser usado
com censura.
Eles são definidos por t < t (r ), para r tempos de falha
(não censurados).
2.10 Estimador do Produto-Limite
O estimador do produto-limite (PL) ou de KaplanMeyer é um estimador não paramétrico da função de
sobrevivência.
Ele coincide com a função empírica de sobrevivência
quando não há censura.
a1 < ... < ak – k tempos de falha distintos (a0 = 0)
d1, ...dk – número de falhas em cada tempo de falha (d0=0)
n1 < ... < nk – número de itens em risco em cada tempo de
falha (nk = 0)
 dj 
ˆ
1  
O estimador do PL é: S t  
 j:a j t  n j 

Esta é uma função escada começando do 1 para t = 0 e
alterando-se a cada ak.
É como se a distribuição de falhas se concentrasse nos
pontos a1, ... , ak.
De acordo com a teoria assintótica , média e variância
de Sˆ t  são dados por
S(t) e
S t 
dj

  n n
j:a j t
j
j
dj
H(t) pode ser estimado de forma similar por
Hˆ t    log Sˆ t 
De forma mais simples e intuitiva, podemos estimar
~
dj
H(t) usando H
t  


j:a j t n j
~
dj
que é relacionado ao estimador h (a j )  n j .
Pode-se utilizar análise gráfica do estimador do PL
para avaliação da adequação de modelos Weibul e lognormal.
Exemplo 2.3: Resistência de corda a uma certa tensão
(em unidades codificadas).
Principais interesses:
• Qual a confiabilidade de uma corda após 53 unidades
de tensão ?
• O modelo de distribuição Weibull é apropriado ?
Da tabela 2.2 , Sˆ (53)  0.6849
2
2
e Var(Sˆ (53))  0.6849 0.0112 0.0725
Um IC de 95% para S(53) é dado por
(0.6849-1.69x0.0725, 0.6849+1.69x0.0725)=(0.54,0.83)
Fora 3 pontos isolados a figura 2.17 parece com uma
linha reta.
Investigação similar com modelo lognormal apresenta
os mesmos resultados.
Tabela 2.1 Resistência de 48 cordas
Observações não censuradas
36,3
41,7
43,9
52,4
52,6
52,7
54,8
54,8
55,1
57,1
57,3
57,7
60,7
49,9
53,1
55,4
57,8
50,1
53,6
55,9
58,1
50,8
53,6
56,0
58,9
Observações censuradas pela direita
29,6
33,4
35,0 40,0
41,9
42,5
51,9
53,9
56,1
59,0
52,1
53,9
56,5
59,1
52,3
54,1
56,9
59,6
52,3
54,6
57,1
60,4
Tabela 2.2 Cálculo da amostral da
Exemplo 2.3
j
aj
nj
dj
0
48
0
1
36,3
44
1
2
41,7
42
1
3
43,9
39
1
4
49,9
38
1
5
50,1
37
1
6
50,8
36
1
7
51,9
35
1
8
52,1
34
1
9
52,3
33
2
10
52,4
31
1
11
52,6
30
1
12
52,7
29
1
13
53,1
28
1
(nj-dj)/nj
1,0000
0,9773
0,9762
0,9744
0,9737
0,9730
0,9722
0,9714
0,9706
0,9394
0,9677
0,9667
0,9655
0,9643
S(aj+0)
1,0000
0,9773
0,9540
0,9295
0,9051
0,8806
0,8562
0,8317
0,8072
0,7583
0,7338
0,7094
0,6849
0,6605
dj/(nj(nj-dj))
0,0000
0,0005
0,0006
0,0007
0,0007
0,0008
0,0008
0,0008
0,0009
0,0020
0,0011
0,0011
0,0012
0,0013
A aproximação do erro padrão de S(53) , vem da equação (2.21) e tabela 2.2
1
-1
-3
log(-log(l1))
Weibull
3.0
3.5
4.0
4.5
5.0
4.5
5.0
log(t)
0 1 2
-2
l2
Lognormal
3.0
3.5
4.0
log(t)
3. MÉTODOS ESTATÍSTICOS
AMOSTRAS SIMPLES
PARA
3.1 Introdução
Final do último capítulo: métodos estatísticos simples.
Este capítulo: métodos mais formais, máxima
verossimilhança, inferência bayesiana dinâmica.
3.2 Estimação por Máxima Verossimilhança:
Generalidades
Suponha uma amostra de tempos de vidas t1 ,...,t n
de uma certa população. Todos os ti `s possuem densidade
f (t \ θ) , onde    ,..., m 
1
Caso as observações não sejam censuradas então, a função de
verossimilhança é
n
L    f ti \  
i 1
para observações censuradas (a direita):
A contribuição para a verossimilhança é a
probabilidade de sobrevivência após o tempo de censura.
Separando os dados em conjuntos disjuntos: C - itens
censurados e U - para itens não censurados.



L    f ti \   S ti \  
 iU
 iC

Para outras formas de censura existem outras expressões.
É mais conveniente trabalhar com l    logL 
Estimativa da máxima verossimilhança (EMV) de
θ é θˆ maximizando L(θ), ou l(θ).
Normalmente são obtidos resolvendo
l  
0
 j
, j  1,...,m
Se φ = g(θ) for uma transformação uma a uma, φˆ = g(θˆ ) é o EMV de φ.
Assumindo que q(p;θ) como o quantil 1-p de T, ou
seja, Pr(T≥ q(p;θ) = S(q(p;θ)) = p
 
O EMV de qp;  é q p;ˆ .
P ela teoria assintótica : θˆ ~ N (θ, V ), quando n → ∞
V=J
-1
com J =
∂l 2 (θ)
∂θ jθ k
, 1 < j, k < m
jk
Notação: o elemento(j, k) da matrizV e denotadocomovjk .
P ela teoria assintótica : ˆ ~ N , V , quando n  
V  J -1
 l 2   
 , 1  j, k  m
com J   
   
j k  jk

Notação: o elemento(j, k) da matrizV e denotadocomovjk .
J é a informação observada da matriz
Em particular,
Para uma transformação φ = g (θ) , o resultado assintótico é :
∂g (θˆ )
φˆ ~ N φ, ∑∑
∂j
j=1 k =1
m
m
∂g (θˆ )
v jk
∂θ k
Em particular, para um m = 1, a variância assintótica de φˆ reduz - se à
∂g (θˆ )
∂θ
2
v11
O EMV possui muitas vantagens sobre todos os outros
métodos clássicos de estimação:
• Ele é universal;
• Ele é invariante;
• Ele possui boas propriedades assintóticas:
Consistência, normalidade assintótica e eficiência;
• Distribuição assintótica é facilmente encontrada.
3.3 Máxima Verossimilhança (MV) estimação
: ilustrações
Suponha que temposde vida t1 ,...,t n são observados, sendo xi  logt i  ,
Onde r itens não são censurados e n - r são censurados.
Distribuição Exponencial :
n
l(λ ) = rlog(λ ) λ ∑t i
i =1
∂l(λ ) r
=
∂λ
λ
n
∑t
i
i =1
∂l 2 (λ )
r
=
<0
∂λ2
λ2
A solução para l′
= 0 é um ponto de máximo dado por λ =
r
n
∑t
i =1
λ2
com variância assintótica
r
Distribuição Weibull
1
αη
l(η, α ) = r log η rη log α + (η 1)∑log t i
i∈U
ou no formato da Gumbel :
rμ
σ
(l μ, σ) = r log σ + ∑x i
i∈U σ
n
∑exp
i =1
μ
xi
σ
n
xi μ
∂l
= r + ∑exp
∂μ
σ
i =1
∂l
= rσ
∂σ
∑x
i∈U|
n
i
+ rμ + ∑(x i
i =1
μ )exp
μ
xi
σ
n
∑t
i =1
η
i
i
P ort ant o,o EMV sat isfaz
x
1 n
ˆμ = σˆ log ∑exp i
r i =1
σˆ
n
xi
∑xi exp σˆ
i =1
n
xi
∑exp σˆ
i =1
1 n
∑x +σˆ
r i =1 i
= 0 ( não envolve μˆ )
Mét odos numéricos são usados para det erminar σˆ e depois μˆ
A segunda derivada do EMV ( usado para a variância assintótica ) são :
∂2 l
r
=
∂μ 2
σˆ 2
∂2 l
=
∂μ∂σ
∂2 l
= r
∂σ 2
n
∑
i =1
n
∑
i =1
μˆ
xi
σˆ
exp
μˆ
xi
σˆ
μˆ
xi
2
exp
σˆ
μˆ
xi
σˆ
Exemplo3.1: Os dados do exemplo 2.2 possuem n = 13 r = 10 e
∑t
i
= 23.05
()
Modelo exponencial : λˆ = 0.434 e s.e λˆ = 0.137.
Figura 3.1 mostra um plot de l(λ ) .
Exemplo3.2 : Os dados no exemplo 2.1 acompanhados das est imativas de ML
(momentos& est imativas gráficas)
Modelo Gumbel : μˆ = 4.405 (4.390 & 4.35) e σˆ = 0.476 (0.416 & 0.435)
Modelo Log - Normal: μˆ = 4.150 (4.150 & 4.2 ) e σˆ = 0.522 (0.534 & 0.56)
Os quantis podem ser est imados através de uma distribuição ajustada,
visto que log N (μˆ , σˆ ) .
L
o
g
V
e
r
s
i
m
l
h
a
n
ç
-26 -24 -2 -20
0 .
0
2 .0
4 .
0
6 .
1
8 .1
0 .2
L a mb d a
Figura 3.1 Log-Verossimilhança para o tempo de vida de
um componente de um avião com distribuição
Exponencial. A reta no gráfico foi feita para mostrar o
intervalo de confiança 95% para lambda, baseado em W.
O EMV do quantil 0.1 de uma distribuição log - Normal é
exp(μˆ - 1.282σˆ ) com Φ(- 1.282) = 0.1
A tabela 3.2 abaixo reúne as estimativas dos quantis (s.e estão
entre parenteses)
Quantil
Weibull Log-Normal
mediana 68.7(8.0) 63.4(6.9)
10% inferior 28.1(6.3) 32.5(4.8)
Note a grande diferença na calda
3.4 Intervalos de Confianças e Testes
  ,,n está divida dentro (A, B das
dimensões ma e mb.
Nós
testar a Hipótese H: A A
^ interessa
^
( θ ( A ) , θ ( B )  é um EMV de (A,B
A é um EMV de B de H.
Estão disponíveis dois procedimentos:
1) Onde temos
^ H , que:
^
^ para
(
B
)
(
(
A
)
W(A) = 2{i  θ , θ   iA, θ B ) A ~ c2 ma),
é aproximadamente.
Grandes valores de W  grandes diferenças em
comum com log – máxima verossimilhança  grande
suporte contra H.
O teste da relação da MV rejeita H se W(A  c2(
(ma)
onde c2( (ma) é 1 -  quantil de c^2ma
(A)
são dadas por:
^As regiões de confiança para θ
{θ ( A ) :W(A,c2 ma)}
^
(A)
2) Seja VA a variância
^
^ assintótica de θ ^ 
Então:W*( θ ( A )  θ ( A )  ATVA-1 θ ( A )  A,
distribuição aproximada c2 ma) de H.
As regiões de confianças e os testes são obtidos a partir
das informações acima com W* substituindo W.
Em particular, para θ(A) escalar, a região de confiança
torna-se:
^
^
½
(
A
)
(
[
z /2 VA ,  A) + z /2 VA ½] ( um intervalo)
onde é usado o fato de c2( 1) = [N(0,1)].
Embora W e W* sejam assintoticamente equivalentes e
geralmente similares, preferimos W pela reparametrização acima e por ser invariante.
Exemplo 3.1 (cont.): Temos o modelo exponencial,
onde θ(A) = , ma = 1 é W*( ) = 2{- 18.35 - 10log  +
23.05}
O intervalo de confiança (IC = 95%) baseado em W* é
dado por: {:W()£ 3.84}=[.22,.76] e como 3.84 =
c0.52(1). (Figura 3.1)
Então o intervalo de confiança (IC = 95%) baseado em
W*é:
[.434 – 1.96 x .137, .434+1.96 x .137] =
[0.17,0.70](simétrico)
O intervalo de confiança (IC=95%) é dado por:
[0.21,0.74] baseado em 2r/ ~ c22r).
Exemplo 3.2 (cont.): Desejamos testar a hipótese H de
exponencialidade.
Usando o modelo Gumbel temos: A  , B  , ma
 e H é   
W*(1) = 49.56 >> 3.84  H0 rejeitado.
Logo, W(1) = 15.50 confirma a rejeição de H como
esperado.
O Exemplo 3.2 considera a hipótese de recursividade:
onde a Exponencial é um caso especial da Gumbel.
São mais difíceis de considerar a hipótese de não –
recursividade para o tratamento clássico estatístico.
3.5 Bondade do Ajuste
Enfoque formal: Encaixar o modelo dentro de uma
classe de modelos ( Exemplo 3.2) ou usar a forma
excelente.
Técnicas Gráficas
Plote o gráfico de QQ: Seja  e , parâmetros de
locação e escala, (onde é estimador do PL) e F é a
função de distribuição para  =  e  = 
O plote dos pontos [aj , F pj)] deveria ser linear.
^
Plote PP: junte os pontos (pj , F aj , θ ). (esta linha
deveria ser y = x)
Pode ser usado fora do modelo de locação de escala.
Plote SP: Para estabilizar a variabilidade de PP, plote a
transformação y = (2/)sin-1x em ambos os eixos.
As figuras 3.2 e 3.3 mostra os dados do exemplo 2.3
plotado em PP e SP.
Goodness-of fit
Valores Esperados
1.00
0.80
0.60
0.40
0.20
0.00
0.00
0.20
0.40
0.60
0.80
1.00
Valores Observados
Figura 3.2 Plote PP da Weibull para os dados da resistência
da madeira
Valores Esperados
1.00
0.80
0.60
0.40
0.20
0.00
0.00
0.20
0.40
0.60
0.80
1.00
Valores Observados
Figura 3.3 Plote SP da Weibull para os dados da resistência
da madeira
3.6 Elementos de Estatística Bayesiana
Incorpora informação subjetiva sobre o problema
(experiência anterior).
É feita através da especificação de uma distribuição a
priori P().
Informação a priori vaga: a análise é guiada pela
informação dos dados.
Assuma, como antes, uma amostra t = ( t1, , tn ) com
densidade f (t ;).
Isto é combinado com a priori e leva a
n
p( ) f (ti ; )
Fórmula
i 1
P( | t ) 
p(t )
de Bayes
Válido para  e t discreto e contínuo.
P( | t) é a densidade a posteriori (dado os dados t).
Como t é constante, A fórmula de Bayes pode ser
simplificada em
P( | t )  p( ) L( )
A constante removida P(t) pode ser recuperada por
P (t )   p ( ) L( )d
Estimativas a posteriori para  são obtidas através de
medidas de locação de P( | t).
Exemplo: Considere os dados do exemplo 2.2 com
modelo exponencial
É conveniente atribuir a priori (distribuição gama).
L( )  10e-23.05
Combinada com a verossimilhança de forma simples
(priori conjugada).
p( )  a1e-b
Para especificar os valores de a e b assuma que
acredita-se que  está próximo de 0.5 e que é pouco
provável que ele seja menor que 0.2.
Então, tome a moda da priori igual a 0.5 e P( < 0.2) 
0.05
a = 3 e b = 4.
A posteriori é  | t ~ Gamma (13, 27.05) com moda
0.444 e média 0.481 (figura 6.1).
p(λ | t)∝( λ3-1e -4λ )(λ10 e -23.05λ ) ∝λ12 e 27.05λ
O desvio padrão a posteriori (priori) é 0.133 (0.175).
Regiões de confiança são facilmente obtidas da
posteriori.
Particularmente úteis são as regiões de maior
densidade a posteriori (HPD).
Por exemplo, o intervalo HPD de 95% para todo  é
[0.231, 0.758].
Interpretação da região HPD é simples diferentemente
das regiões de confiança clássicas.
Inferência sobre funções paramétricas são obtidas de
maneira similar.
Assuma interesse na confiabilidade em um certo tempo
t0 .
Para o modelo exponencial isto é S ( t0 ; ) = e-lt0, uma
função de .
A posteriori completa de S ( t0 ; ) pode ser obtida.
Como exemplo, 0.01 = Pr ( < 0.225 | t) = Pr (S ( t0 ;
) > e-0.225 t0 | t)  a probabilidade a posteriori de que S
(t0; ) exceda e-0.225 t0 é 0.01.
Predição: assuma que se está interessado no tempo de
vida S de um novo item.
Inferência deve ser baseada na distribuição de S | ( t1,
... ,tn ).
p(s | t)   p(s |  , t) p(  | t) d   p(s |  ) p(  | t) d
(S independente de t dado )
Por exemplo, a densidade do tempo de vida de um
novo item é

p(s |  )   (e
0
 s
12  27.05 
)( e

)d   13e ( 27.05 s )  d  (1  0.037s) 14
0
(a constante de proporcionalidade é 0.481).
Informação a priori vaga é usualmente representada
por adequados valores pequenos dos parâmetros da priori
conjugada.
No exemplo, pequenos valores de a e b, como 0.5.
O intervalo de 95% de confiança a priori é [0.001,
5.024]. (muito grande)
A posteriori é Gama (10.5, 23.55) com média 0.446 e
desvio padrão 0.138.
O limite de uma distribuição a priori vaga é uma priori
não informativa.
No exemplo, isso é obtido fazendo a, b  0 => p(l) l1
(priori imprópria).
A posteriori é Gama (10, 23.05) com média 0.434 e
desvio padrão 0.137. (similar aos resultados da inferência
por máxima verossimilhança).
A priori não informativa é um meio para obtenção da
posteriori na ausência de informação a priori.
3.7 Outros Tópicos em Inferência Bayesiana
(§6.4 do livro texto)
Especificação de prioris
Não é necessário que a especificação seja muito
precisa.
Próxima seção: análise Bayesiana com especificação a
priori parcial.
Verificação de inconsistências: análise pré-posteriori.
Análise conjugada é conveniente mas nem sempre
apropriada.
Prioris não informativas devem ser usadas com
cuidado: pode levar a absurdos.
Parâmetros de distúrbio
Suponha que os parâmetros são divididos em (y, )
onde y é de interesse
 é somente necessário (distúrbio) e pode ser eliminado
facilmente via
p ( | t )   p ( ,  | t )d
Densidade marginal
-t
Exemplo: tempos de falha Weibull com S(t) = e
A verossimilhança é
n
r

r
 s ( )
s
(

)

t
com

i
L( ,  )  ( )
t e

i 1
i
i 1
Assuma que a priori é
a 1 b
p( ,  )  p( ) e
 p( ,  | t )  p( ) 
r
Integrando com relação a  temos
a  r 1
r
t e 
 [ b  s ( )]
i
i 1

r

  i1ti 

p( | t )  p( ) [b  s ( )]ar 




r
Entre chaves: p(t | ) - verossimilhança marginal (ou
integrada) de .
Fator de Bayes
Considere um teste Bayesiano de uma dada hipótese
H.
Se H é uma região então Pr(H | t) reflete a crença em
H a posteriori.
Para o exemplo em §3.6 e H :  < 0.25, temos Pr(H | t)
= Pr(  < 0.25 | t) = 0.2 => H é rejeitada
Similarmente, pode-se usar a razão de chances a
posteriori dada por
Pr(H | t ) Pr(H ) p(t | H )

Pr(H | t ) Pr(H ) p(t | H )
A primeira razão do lado direito é a razão de chances a
priori e o segunda é o fator de Bayes.
Isso representa a razão relativa de verossimilhanças
entre as duas hipóteses H e H .
Maiores valores do fator de Bayes  maior apoio dos
dados em H.
O fator de Bayes é útil quando deseja-se testar
hipóteses nulas bilaterais H :  = 0
Note que Pr(H | t) será sempre 0.
Há alguma controvérsia a respeito de quão adequado é
testar hipóteses nulas bilaterais.
No exemplo, o fator de Bayes para H : l = 0.25 é 0.81
 a crença no valor 0.25 é reduzida pelos dados.
3.8 Modelos Bayesianos Dinâmicos
Baseado em uma distribuição exponencial por partes
para os tempos de falência com suposição explícita de
conexão entre intervalos.
Distribuição E. P. (Exponencial por Partes): o risco é
constante nos intervalos.
É geralmente razoável assumir funções de risco
contínuas  algumas conexão entre valores de risco em
intervalos sucessivos.
Forma matemática adotada para simplicidade: passeio
aleatório na escala log
log λ i = log λ i 1 + w i onde E(w i ) = 0 e Var (w i ) = Wi
wi - termo de perturbação permitindo o aumento da
incerteza como um movimento direto do intervalo.
Conseqüências deste modelo:
(1) Preserva a posição:
(2) Aumento da incerteza :
Um artifício útil para determinar valores para Wi´s:
fatores de desconto controlam a quantidade de informação
(medida de precisão) passando direto dos intervalos.
Fator de desconto é um número entre 0 e 1(geralmente
fechado para 1).
1) se o desconto é fechado para 0 – nenhuma
informação passa direto do intervalo. Estimação M.V.
(máxima verossimilhança) com distribuição E.P.
(exponencial por partes), estimador P.L. (produto-limite
ou Kaplan Meyer).
2) se o desconto é 1 – toda informação passa direto
do intervalo
Parâmetros são os mesmos  tempo de falência
exponencial.
Função de verossimilhança
De §2.7, a verossimilhança para um dado indivíduo é
∞
Πλ
xi
i
exp{ b i λ i}
i =1
onde Xi é o indicador da falência no intervalo Ii
bi é o tempo observado em Ii, para este indivíduo.
Para uma amostra de tamanho n, a verossimilhança é
∞
L = ΠLi
i =1
L i = λ idi exp{ λ i a i }
onde
onde di é o número de indivíduos observados até a falha
em Ii
ai é o tempo total observado em Ii para a amostra
Li é a verossimilhança para λi baseado nos eventos
observados em Ii dado Di-1, a informação do intervalo
anterior (ver final de §2.7).
De fato, o produto das verossimilhanças vai de i = 1 a i =
N onde N é o indexador do último intervalo com tempo
de falência observado (censurado ou não censurado).
Análise Seqüencial e Distribuições a Priori
Assuma que λ i Di 1 ~ G(α i 1, γ i 1)
A Priori assumida para
λ i Di
1
é
G ( c i α i 1 ,c i γ i 1 )
onde
0 < c i ≤1
Conseqüentemente:
1) E(λ i D i 1 ) = E(λ i 1 D i 1 )
2) Var (λ i Di 1 ) = c i 1Var (λ i 1 Di aumento
da incerteza
1)
Atualização da distribuição de λi feita direta da fórmula
de Bayes.
p(λ i Di ) ∝ E(λ i Di 1 )L i
∝ λ iciαi 1 1exp{ ci γ i 1λ i}λ idi exp{ a i λ i}
∝ λ iciαi 1 +di 1exp{ (ci γ i 1 + a i )λ i}
e λ i D i ~ G(c i α i
1
+ d i , c i γ i 1 + a i ) (análise conjugada)
As análises procedem do aumento de i para i+1 como
antes dito.
Inicia em i = 0 indo para i = N, o último intervalo com
informação de dados.
{ci} controla a suavidade da função risco:
ci → 0: sem passagem de informação
ci = 0: passagem total de informação
Especificação dos ci´s
Dado que λ i 1 D i 1 ~ G(α i 1 , γ i 1 )
Usando o método Delta
E(logλ i 1 D i 1 ) = log(
αi 1
1
) e Var (logλ i 1 D i 1 ) =
γi 1
αi 1
Da relação entre sucessivos λ´s:
αi 1
1
E(logλ i D i 1 ) = log(
) e Var (logλ i 1 D i 1 ) =
+ Wi
γi 1
αi 1
Usando novamente o método Delta dado
E( λ i D i 1 ) = (
αi 1
α
1
) e Var (logλ i D i 1 ) = ( i 1 ) 2 (
+ Wi )
γi 1
γi 1 αi 1
De onde um obtém ci = (1 + α i 1 Wi ) 1
Geralmente, Wi é selecionado proporcional
ao
tamanho de Ii.
Quanto maior o intervalo, mais informação é perdida.
A proporcionalidade constante é W, a variância da
perturbação supera uma unidade de tempo.
Fatores de desconto são associados diretamente com
W: tem que ser especificado só a primeira vez para um
modelo dado.
Inferência
Inferência é baseada na distribuição predita de um
novo tempo de falência S baseado em DN, o (total) de
informação do dado.
Interesse particular: sobrevivência predita S(s D N ) e
risco predito h(s D N )
Estes são obtidos após a integração fora dos λ´s com
respeito a sua distribuição suavizada (ou filtrada) p(λ i D N )
Estas distribuições são obtidas via um algoritmo
recursivo.
Modelo de Seleção
Um modelo é especificado pela escolha de :
1) Priori para λ 1 D 0
2) fator de desconto
3) grade de intervalos
Modelos M1 e M2 podem ser comparados via seus
fatores de Bayes
p( t M 1 )
p( t M 2 )
Cada verossimilhança marginal é obtida após a
integração fora dos parâmetros como segue
N
p(t M) = ∏ p({eventos observados em I i } D i 1 , M)
i =1
N
=∏∫
p({eventos observados em I i } λ i , D i 1 , M) p (λ i D i 1 , M) dλ i
i =1
N
=∏ ∫
L i p (λ i D i 1 , M) dλ i
i =1
facilmente obtido.
Tempo de Falência do Sistema de Telecomunicação
Período de Observação de 20 de Maio de 1985 a 31de Outubro de 1985
( Z = dias da instalação para a falência, cancelamento ou data de encerramento)
c = censurado e u = não censurado
Z
Z
Z
Z
Z
Z
164c
3u
164c
164c
163c
2u
155c
155c
155c
139u
45u
150c
150c
150c
149c
147c
101c
146c
1u
143c
139c
139c
139c
139c
138c
135c
135c
135c
1u
134c
163c
163c
163c
163c
163c
152c
152c
152c
152c
94u
149c
149c
149c
149c
149c
143c
143c
142c
10u
141c
40u
138c
138c
138c
138c
13u
134c
134c
134c
134c
77u
162c
162c
73c
63u
151c
151c
151c
151c
151c
149c
149c
149c
115u
148c
141c
141c
34u
140c
140c
138c
138c
138c
138c
137c
133c
133c
133c
133c
133c
161c
160c
160c
67u
141c
151c
151c
151c
90c
151c
148c
147c
147c
147c
147c
140c
140c
140c
140c
140c
137c
137c
137c
137c
137c
64c
133c
133c
133c
133c
156c
151c
147c
54u
137c
3.9 O Estimador Atuarial
Utilizado em tabelas de mortalidade onde muitas vezes
dados estão agrupados. Assume-se que a distribuição é
contínua e divide-se o tempo em intervalos geralmente
iguais onde a taxa de falha é constante.
Ex: Em tabela de mortalidade, divide-se tempo (tempo
de vida de população) em intervalos 0 – 1, 1 – 5 (ou 0 –
5), 5 – 10, 10 – 15, ... anos.
Raramente faz-se divisão ano a ano.
Suponha que a população tem n indivíduos morrendo
em um ano com idades y1, ..., yn. (Não há censura)
Se o indivíduo morre no intervalo i , sua contribuição à
verossimilhança é:
i 1
i  e
 j t j t j 1 
 e  y t
i
i 1

j 1
A verossimilhança total é dada pelo produto das
contribuições individuais

d

 e
n
 xik
 i 
k 1
i
i
i 1
 verossimilhança fatora em i
onde di = # de mortes no intervalo i

, y k  t i 1
 0

xik   y k  t i 1 , y k  t i 1 ; t i 
t t
,
yk  t
 i i 1
i
x
ik
indivíduo k morre antes do intervalo i)
(indivíduo k morre no intervalo i)
(indivíduo k morre depois do intervalo i)
é o tempo total em risco no intervalo i
Por analogia, na estimação do modelo exponencial, o
EMV de i é:

i 
di
 xik
Se todos os indivíduos morrem no final dos intervalos,
x
ik
k
 t i  t i 1   ri
onde ri = # de indivíduos observados no intervalo i.
Se também há censura, fórmulas não se alteram (apenas
ri será diferente).
Vamos supor agora que a censura também está sujeita
a um mecanismo:
 Aleatório cuja taxa é constante ao longo dos mesmos
intervalos
 Independente do mecanismo de falha (mortalidade)
Inferência acima não é alterada pela independência.
Normalmente, em tabelas de mortalidade, dados são
fornecidos em forma grupada, isto é, só são fornecidos:
di = # de mortes no intervalo i
mi = # de censuras no intervalo i
Temos 3 grupos de indivíduos em cada intervalo:
i) ri – di – mi - sobrevivem ao intervalo i
ii) di
- morrem no intervalo i
iii) mi
- são censurados no intervalo i
Vamos supor que taxa de censura é i no intervalo i
Já
vimos
que
verossimilhança
fatora
verossimilhanças condicionais à história passada.
em
A contribuição dada à verossimilhança do intervalo i
de cada um dos 3 grupos acima dada sobrevivência até o
início do intervalo é dada por:
(i) Pr (Y > ti , C > ti | Y > ti-1 , C > ti-1)
(ii) Pr (Y  ti , C > Y | Y > ti-1 , C > ti-1)
(iii) Pr (Y > C , C  ti | Y > ti-1 , C > ti-1)
falha
(iii)
(i)
ti
(ii)
ti-1
censura
ti
De fundo bi = ti – ti-1 temos:
(i) Pr(Y>ti | Y>ti-1).Pr(C>ti | C>ti-1) = exp{-bii}exp{-bii}
= exp{-bi(i+i)}

f Y  y | Y  t i 1    f c c | C  t i 1 dcdy
ti

(ii)
ti 1
y
ti

ti 1
y
  i exp i  y  ti1   i exp i c  ti1 dcdy
  i exp i  y  ti 1 exp i c  ti 1 dy
ti
ti 1

(iii)
i
i   i
i
i   i
1  exp bi i   i 
1  exp bi i   i  , por analogia a (ii)
Logo, a verossimilhança do intervalo i é dada por
(ri
di
mi )bi (λ i + θi ) + d i log
λi
θi
+ mi log
+ (d i + mi )log{1 exp[ bi (λ i + θi )]
λ i + θi
λ i + θi
Nenhum outro fator de verossimilhança depende de i e i
EMV de i e i podem ser obtidos a partir da verossimilhança
acima dando
 r  d i  mi
di
i  
log i
bi d i  mi  
ri





e i  
 r  d  mi 
mi

log i i
bi d i  mi  
ri

Normalmente,
d i  mi
ri
é pequeno
x2
Fazendo aproximação log1  x   x 
2

di
d
i 
 i|
bi ri  di  mi  2 bi ri
temos
Equivale a assumir que mortes e censuras se distribuem
uniformemente nos intervalos.
A probabilidade de sobrevivência a um intervalo é
S t i | Y  t i 1   e ibi

e
i bi
 e di
e pode ser estimada por
ri|
No caso específico de uma tabela de mortalidade, n é
bastante grande e a única informação é d1, d2, ..., ou seja,
número de mortos em cada intervalo.

Podemos calcular ri pois ri   d j
j i
são fornecidos.
, mas xik não
Os problemáticos são os indivíduos que morrem (ou
são censurados) no intervalo.
Suposição: Dada a massa de indivíduos, é razoável
supor que indivíduos morrem (ou são censurados)
aleatoriamente no intervalo.
Logo,
para esses indivíduos é tomada como
 # de indivíduos t  t
t  ti 1

 que morremno   i
i 1
 ti 1   d i  i

  2
2

 int ervalo

Como ri = ri+1 + di temos  x
ik
k
se não há censura
 t  t 

  xik   xik  di  i i 1  ti 1   ri 1  ti  ti 1 
 2

Di
Si
onde Di = { indivíduos que morrem no intervalo i}
Si = {indivíduos que morrem após intervalo i}
 di 
ti  ti 1   r t  t 
i 1 i
i 1
2
d

 ti  ti 1   ri 1  i 
2

= (t i
ti
)
1
ri
di
2
 ti  ti 1  ri|
A suposição acima é razoável: da é falha apenas no
intervalo 0-1 ano onde a tendência à falha é nitidamente
maior perto de 0.
dL
t i  t i 1   ri|
Logo, i é usualmente estimada por
A probabilidade de sobrevivência a um intervalo é
S ti | Y  ti 1   e i ti ti 1 

e pode ser estimada por e
Muitas vezes
i ti ti 1 
 e di
ri|
di
é pequeno. Fazendo a expansão de
ri|
Taylor em torno de 0, obtemos
e di
ri|
 1
di
ri|
~
~
que é denotado por pi ou 1  qi
e substituindo esses valores em S t i    S t j | Y  t j 1 
i

j 1
i
temos o estimador atuarial S ti    ~p j
j 1
Observe que o estimador atuarial difere do estimador
PL pela troca de ri por quanto t = ti , i = 1, 2, ...
Nos outros pontos ele é contínuo e o PL é tipo escada
1
0
3.10 O Estimador Bayesiano
Quando não se assume nenhuma distribuição específica
para os tempos de falha a própria Fd F ou a f.s. S torna-se
o parâmetro da distribuição.
Convenciona-se dizer que o problema é não paramétrico
pois a dimensão do parâmetro é infinita.
Temos que construir priori sobre F(t) ( ou S(t) ),  t
[0;)
Ferguson (1973): Seja  uma medida finita R+, isto é,
 ( [a,b] ) = c > 0,  ( [a,b]A )  c e  ( [0,) ) < 
Ex:   A  eu du   ( [0,) ) = 1 < 

A
A distribuição
P é um processo de Dirichlet se qualquer
partição B1, ..., Bk de R+, Pr(B1), ..., Pr(Bk) tem
distribuição Dirichlet com parâmetro ((B1), ..., (Bk) )
  1 , , k  ~ D1 , ,  k 
~
onde 1 = (A1), ..., k = (Ak)
e, portanto, p( ) 
~
k


i
i 1
i 1
X   X1,, X n  é uma amostra aleatória de um
~
processo de Dirichlet se:
Pr(X1A1, ..., XnAn | P(A1), ..., P(An) ) =
Pode se mostrar que P r X 1  A1  
n
 P A 
i
i 1
  A1 
  0;  

O processo de Dirichlet funciona como a priori para a
distribuição amostral. Assim, por exemplo,
F t   Pr 0; t 
Logo tem densidade
e
EF (t ) 
 ~ D1; 2 
 1
 1
f F (t )  F (t ) 1  F (t )
1
1
  0; t  

1   2   0;   
 corresponde ao valor esperado a priori para F.
2
Se não há censura na amostra, a distribuição é
conjugada e a posteriori de P também é um processo de
Dirichlet com parâmetro  onde:
n
 t    t    I  yi  t  , se  é contínuo,  não é mais.
i 1
Se há censura, a distribuição não é conjugada e a forma da
posteriori complica. Pode-se obter a esperança a
posteriori de F(t) ou S(t),  t > 0

Usando essa esperança como estimador S (t ) temos
  t;    r t  0 l    yi ;    r  yi 
S (t ) 

  0;    n
j  k 1   y j ;     r  y j   c j

yl < t  yl+1
c = k, ..., m onde
y1, ..., yk são as observações não-censuradas
yk+1, ..., ym são os diferentes valores das observações
censuradas
ck+1, ..., cm são o número de observações censuradas =
yk+1, ..., ym
r(t) = # de observações à vista em t

Se não há censura, S (t ) 
yi-1  t  yi , i = 1, ..., n
  t ;     n  j  1
  0;     n

No caso geral, se 0, S (t )  estimador PL

Quando se  > 0, S (t ) comporta-se da seguinte forma
A(n)
Abordagem completamente não paramétrica sugerida por
Hill (1968) baseada na hipótese A(n):
Sejam X1, ..., Xn+1 permutáveis com distribuição P.
Suponha que observa-se X1 = x1, ..., Xn = xn onde x1 < x2
< ... < xn (possível pela permutabilidade)
Sejam I(0), ..., I(n) intervalos dados por I(i) = (xi, xi+1)
onde x0 = 0 e xn+1 = 
Então a distribuição preditiva de Xn+1 dado X1 = x1, ...,
Xn = xn dá Pr  X n 1  I (i )  
1
n 1
Eventuais empates nos xi’s podem ser separados
acrescentando  pequeno a um deles
No caso de censura a situação complica pois não
conhecemos todos os xi’s. Resolve-se também de uma
forma não paramétrica.
Suponha
x4 e x5 são tempos de censura.
O objetivo é calcular as probabilidades preditivas para
I(0), I(1), I(2), I(3). Isso é feito considerando todos os
possíveis valores (não censurados) de x4 e x5.
Supondo que x4  I(1) e x5  I(2) temos sob a hipótese
A(5) que a probabilidade
1
de I(0) é 6 , de I(1) é 1 6  1 6  13 , de I(2) é 1 6  1 6  13
e de I(3) é 1 6
Se supomos, por exemplo,
Implícito nos cálculos acima:
Falha dos censurados poderia ocorrer em qualquer
ponto no intervalo de censura, ou seja, o ponto de censura
é trazido de volta até o início do intervalo.
Ex.: x4 poderia ocorrer em qualquer ponto de I(1)
O cálculo é portanto uma aproximação (que fornece cota
superior).
Defina Z1, ..., ZN – tempos de falha a ser observados
Z – tempos de falha a ser previsto
n < N – número de falhas observadas
x   x1 ,  , x n  – valores observados de falha
~
yn+1, ..., yN – valores observados de censura
IEC = {Zj > yj, j = n+1, ..., N} Info. Exata de Censura
IPC = { Zj > Uj, j = n+1, ..., N} Info. Parcial de Censura
onde Ui é o maior xi | xi < yi
Ex.: No exemplo, temos x1 < y4 < x2 < y5 < x3  U4 = x1
e U5 = x2
Ao invés de calcularmos Qi = Pr( Z  I(i) | IEC,
i = 0, 1, ..., n calcularemos
),
Pi = Pr( Z  I(i) | IPC, x
~
Como IEC  IPC, Qi = Pr( Z  I(i) | IEC, IPC,

)=
Pr( Z  I(i),IEC| IPC, x )
~
Pr(IEC| IPC, x )
~

Pr( IEC| IPC, x, Z  I(i),)
~
Pr(IEC| IPC, x )
 Pi
~
Como, dado IPC, pouca informação adicional sobre
IEC é fonecida por [Z  I(i)] principalmente em amostras
moderadas ou grandes, Qi  Pi.
Para o cálculo dos Pi, define-se ci - # de observações
censuradas em I(i)
i
C i   c k , i = 0, 1, ..., n
k 0
i = 1 / (N - (i - 1) - Ci)
Pi1  i1   1   j  , i = 0, 1, ...,
i
Então P0 = 0 e
n-1
j 1
Prova: (i = 0): Se c0 = 0, C0 = 0 e A(N)  P0 = 1 / (N+1)
Se c0  0, sob IPC, essas c0 observações são colocadas
no início de I(0) e portanto não trazem nenhuma info.
Podemos considerar apenas N-c0 = N-C0 observações
restantes.
A(N-C0)  P0 = 1 / (N-C0+1)
Sob IPC, c1 observações são colocadas no início de
I(1) e para o cálculo da parcela acima (condicional a Z >
x1) não trazem nenhuma info. assim como x1.
Podemos então sob A(N – 1 – c0 – c1) = A(N – 1 – C1)
obter Pr( Z  I(1) | Z > x1 , IPC, x
~ ) = 1
Logo, P    1   
1
1
0
O mesmo raciocínio pode ser seguido e usando
indução mostra-se o resultado.
A substituição de IEC por IPC faz com que
i = Pr( Z  I(i) | Z > xi , IPC, x
) produzam cotas
~
superiores para Pr( Z  I(i) | Z > xi , IEC, x
)
~
Observe que
Cotas inferiores podem ser obtidas se censuras são
trazidas para cima.
IPCI = {Zj > Uj , j = n+1, ..., N} Info. Parcial de
Censura Inferior onde Uj é o menor xi | xi > yj
Nesse caso Q0  Pr( Z  I(0) | IPCI, x
)
~
por A(N)

1
N 1
,
i
Assim como Pi +1 = λ i +1
(1 λ ) = λ
j
j=1
i
s
i +1
(1 λ )
s
j
j=1
Podemos definir iI  Pr( Z  I(i) | Z > xi , IPCI, x
~ )
s
1
Por analogia com i temos que , iI 
N  (i  1)  Ci 1
i = 0, 1, ..., n
com C-1 = 0
Definindo então Pi = Pi s
temos
, Pi I = Pr( Z  I(i) | IPCI, x
)
~
Além disso,
1I = Pr( Z  I(1) | Z > x , IPCI, x
~ )  Pr( Z  I(1) | Z >
1
s
x1 , IEC, x
~ )  1
Multiplicando as duas inequações tem-se
Em geral
4. MODELOS DE REGRESSÃO PARA
DADOS DE CONFIABILIDADE
4.1 Introdução
Até agora, itens pertencem a mesma população.
Às vezes, outras variáveis afetam os tempos de falha.
Exemplos: tensão, pressão, temperatura (confiabilidade), idade, tratamento, sexo (análise de sobrevivência).
Estas variáveis são chamadas covariáveis, e devem ser
incorporadas ao modelo.
Elas podem ser contínuas (tensão, pressão,
temperatura, idade) ou discretas (tratamento, sexo).
De agora em diante, analisaremos dados de falha com
covariáveis (regressão).
Em confiabilidade, a maioria dos modelos são
baseados na distribuição Weibull.
Outra opção é a lognormal: após a transformação log
nos tempos de falha, podemos fazer uso da teoria normal
e da regressão padrão.
Em geral, prefere-se a distribuição Weibull devido a
facilidade do seu uso com dados censurados e devido a
forma da função taxa de falha.
Propósito do estudo: determinar o quanto T é afetado
por x (covariáveis).
Exemplo: x – tensão; T pode decrescer com x.
§ 4.2 a 4.6 descrevem diferentes modelos.
§ 4.7 a 4.9 lidamos com modelos baseados na distribuição
Weibull.
Capítulo 5 lida com modelos de taxa de falha
proporcional.
Capítulo 6 lida com modelos Bayesianos dinâmicos.
4.2 Modelos de Tempo de Vida Acelerado
(ALM)
Suponhamos tempos de falha sujeitos a uma carga.
ALM: tempo de falha é o produto de uma função da
carga e do tempo de falha padrão.
Tempo de falha padrão: tempo de falha para um nível
padrão de carga
Pr (T  t | x) = S (t; x) = S0 (t yx )
aonde S0 é a sobrevivência básica (padrão) e yx é uma
função positiva de x.
Quando x está no nível básico: yx = 1. Para outros
valores de x, tempo de falha é acelerado (multiplicado)
por yx.
Exemplo: S0 é exponencial unidade (S0 (t) = e  t ) e
yx = x 
 S (t; x) = S0 (t x  ) = exp {t x  
No ALM: Tyx tem distribuição básica ou log T =  log
yx + log W, aonde W ~ P (não depende de x).
Uma especificação comum é log yx = xTb (similar aos
modelos lineares normais).
A função taxa de falha h (t; x) é
 log S 0 t y x 
 log S t; x 

 yx h0 (t yx )
t
t
A adequação do ALM pode ser preliminarmente
avaliada pelos seguintes gráficos:
1) Sendo log T =  log yx + “termo de erro”, o gráfico
de log T versus x deve fornecer uma indicação da forma
de y.
2) Se os dados podem ser agrupados em k grupos
homogêneos, os gráficos das estimativas
das funções
de sobrevivência versus log t devem ser cópias
horizontais deslocadas de S0 quando
h t; x   
Sj (t) = S0 (t yx) = S0 {exp (log t + log yj )} = s0 (log t + log yj )
3) Também no ALM com k grupos, os quantis são
proporcionais. Isto é, se qj (p) é o p-ésimo quantil do
grupo j então qj (p) yj = qk (p) yk porque S0 (qj (p) yj) =
Sj (qj (p)) = p, j.
Os quantis 0.1, 0.2, ..., 0.9 para todos os grupos podem
ser estimados de
.
Podemos fazer o gráfico dos quantis estimados para o
grupo j, j = 2, ..., k versus os do grupo 1. Os gráficos
devem ser aproximadamente lineares com inclinação y1
/yj .
Equivalentemente, os gráficos dos log-quantis devem
ser linhas paralelas.
Os gráficos também poderão ser feitos versus a média
dos quantis (ao invés de versus os quantis do grupo 1).
4.3 Modelos de Taxas de Falha Proporcionais
(PHM)
PHM: a taxa de falha é acelerada, isto é, h (t; x) = h0
(t) yx aonde h0 é a função taxa de falha básica (padrão) e
yx é uma função positiva de x.
Quando x está no nível básico: yx = 1.
y
Este nome vem do fato h t ; x1   x
h t ; x 2 
t).
yx
(não depende de
1
2
Uma especificação comum é log yx = xTb. Assim, h (t;
x) = h0 (t) exp{xTb}.
A função de sobrevivência é




S (t; x) = exp   h t; x  dt = exp  y x  h0 t  dt =S 0 t 
t
0
t
0
yx
Se S (t; x) é um ALM e um PHM então
S (t; x) = Sa 0 (t y ax) = S p 0 t y
e a única solução possível é a sobrevivência da Weibull
Sa 0 (t) = Sp 0 (t) = exp {t h } com y px  y ax .
Gráficos preliminares podem novamente serem feitos
depois da separação dos dados em k grupos.
Assim, temos que
log Sj (t) = y j log S0 (t) ou log {log S j (t)} = log y j +
log {log S0 (t)}
e os gráficos de Sˆ j devem ser múltiplos mutuamente na
escala log ou devem ser cópias verticais deslocadas de
cada um na escala log-log.
px
4.4 Modelos de Razão
Proporcionais (POM)
POM: a função taxa de falha satisfaz
de
Chances
1  S0 t 
1  S t; x 
y x
S t; x 
S0 t 
Após diferenciar em relação a t, h t; x  y x S t; x   1  S t; x 
h0 t 
S0 t 
1  S0 t 
A razão taxa de falha é yx quando t = 0 e tende a 1
quando t  .
 ocorre diminuição do efeito de x na taxa de falha a
medida que o tempo cresce.
Avaliações preliminares podem ser feitas baseadas nos
gráficos 1 Sˆ j Sˆ j .


4.5 Generalizações
A relação PHM
log {log S j (t)} = log y j + log {log S0 (t)}
pode ser generalizada como
 1 {log S j (t)} = g 1 (yj ) +  1 {log S0 (t)}.
Isto inclui o POM como o caso especial:  1 (u) = log
{(1  u) / u}.
A relação ALM log q j = log q 0  log (y j / y 0) pode
ser generalizada como  2 (q j ) =  2 (q 0) g 2 (y j / y 0).
Outras generalizações são:
1) ALM generalizado aonde y x também depende de t.
2) PHM generalizado aonde aonde s (l ) é a família
Box-Cox de transformações em s.
3) Modelo de deslocamento no tempo: h (t; x) = h0 (t +
yx) – taxa de falha de ação atrasada.
4) Modelo de taxa de falha polinomial: h (t; x) = y 0x +
y 1x t + ... + y qx t q.
5) Modelo de taxa de falha por partes: h (t; x) = h i (t;
x), para t  I i, i = 1, ..., N.
6) Covariáveis dependentes do tempo: x pode depender
do tempo em alguns casos.
4.6 Modelos Bayesianos dinâmicos (DMB) (cf
Gamerman, 1991)
DMB combina modelos de taxas de falha
proporcionais (PHM) com distribuição exponencial por
partes (similar ao modelo de taxa de falha por partes).
Do PHM, temos que h (t; x) = h0 (t) y x.
Assuma que h0 (t) = exp {b i, 0}, t  I i.
Também, generalize y x para y i, x, para que dependa de
t.
Uma especificação comum é y i, x = exp {xTb i}.
O modelo pode ser escrito como h (t; x) = exp {xTb i},
aonde b i agora inclui uma primeira componente b i, 0 e xT
agora inclui uma primeira componente 1.
O modelo é completado com a relação entre os
parâmetros em intervalos sucessivos, como em § 3.8.
b i = G i b i1 + w i aonde E(w i) = 0 e Var (w i) = Wi.
A matriz evolução G i fixa a parte determinística da
evolução e o termo de erro w i controla o aumento na
incerteza a medida que o tempo passa.
Exemplos:
1) G i = I, é a matriz identidade (passeio aleatório
simples);
 1 tamanho( I i ) 

0
1


2) G i = diag (1, G i 1) aonde Gi1  
(modelo de crescimento generalizado).
A variância pode novamente ser especificada através
dos fatores de desconto.
Pode ter uma para cada parâmetro mas o mais comum
é ter uma para o parâmetro de base b i, 0 e uma para os
coeficientes de regressão.
Geralmente, esta última é próxima a 1 (se 1, o modelo
se torna PHM).
Se não há presença de covariáveis, o modelo se torna o
mesmo de § 3.8.
4.7 Modelos
Weibull
Baseados
na
Distribuição
O Modelo de Regressão Weibull pode ser escrito como:




t 


S t ; x   exp   


  x 

onde
log x  xT 
S0 (t )  exp{t }

S t; x  S0 (ty ax )
 S0 (t )
y px
onde
onde
y ax 
1
ALM
x
y px 
1
PHM

x
1
log
T

log


W
Pode, também, ser escrito como:
x

w


  T  w
 





onde: Pr(W  w)  Pr log     Pr T   x e   exp  e w
 x   




Sabemos que W ~ Gumbel e ainda que  é somente um fator
escala, mas também pode ser dependente de x.
Estimação e Testes
Regressão Linear Simples:
log x  0  1 x
Estimador de Mínimos Quadrados pode ser usado.
Procedimento similar pode ser usado para regressões múltiplas.
Mais formalmente, o Estimador de Máxima Verossimilhança
pode ser calculado como no capítulo 3: parâmetros são
estimados, variâncias assintóticas obtidas e testes de hipóteses
calculados via Teste da Razão de Máxima Verossimilhança.
Esses cálculos podem ser manuseados por GLIM.
Plot de Resíduos
PPplot pode ser usado como em 3.5

ˆ

Os dados transformados ui  S ti ; xi ,  convergem para
uma amostra de variáveis aleatórias U(0,1).
Os ui’s são chamados resíduos generalizados.
Também, zi = log(-log ui) tem distribuição Gumbel.
Um plot de probabilidade pode ser construído como segue:
• Ordene os zi’s em zi:n, ..., zn:n
• Faça o gráfico zi:n contra log{- log(1-pi:n)} onde pi:n é:
1
n
ou
i  0.5
n
Dados não censurados
n  0.5
n  j  0.5
1

n j i , jJ n  j  1.5
Dados censurados
• Checar se o gráfico é da forma y=x
Resíduos generalizados são também proveitosos quando
plotados contra covariáveis (ver 4.8)
4.8 Um Exemplo: Resistências de Fibras de
Carbono e Pacotes
Dados das tabelas 4.1 e 4.2 são analisados.
Tabela 4.1
2,247
3,726
4,118
4,632
5,571
2,640
3,727
4,141
4,634
5,684
2,842
3,728
4,216
4,636
5,721
2,908
3,783
4,251
4,678
5,998
3,099
3,785
4,262
4,698
6,060
1,901
2,518
2,856
3,235
3,554
2,132
2,522
2,917
3,243
3,562
2,203
2,525
2,928
3,264
3,628
2,228
2,532
2,937
3,272
3,852
2,257
2,575
2,937
3,294
3,871
1,312
2,021
2,301
2,535
2,773
3,128
1,314
2,027
2,301
2,554
2,800
3,233
1,479
2,055
2,339
2,566
2,809
3,433
1,552
2,063
2,359
2,570
2,818
3,585
1,700
2,098
2,382
2,586
2,821
3,585
1,339
1,852
2,171
2,386
2,604
3,174
1,434
1,862
2,172
2,39
2,62
1,549
1,864
2,18
2,41
2,633
1,574
1,931
2,194
2,43
2,67
1,589
1,952
2,211
2,431
2,682
(a) Length 1mm
3,126 3,245 3,328
3,786 3,896 3,912
4,326 4,402 4,457
4,738 4,832 4,924
(b) Length 10mm
2,350 2,361 2,396
2,614 2,616 2,618
2,977 2,996 3,030
3,332 3,346 3,377
3,886 3,971 4,024
(c) Length 20mm
1,803 1,861 1,865
2,140 2,179 2,224
2,382 2,426 2,434
2,629 2,633 2,642
2,848 2,880 2,954
(d) Length 50mm
1,613 1,746 1,753
1,974 2,019 2,051
2,27 2,272 2,28
2,458 2,471 2,497
2,699 2,705 2,735
3,355
3,964
4,466
5,043
3,383
4,050
4,519
5,099
3,572
4,063
4,542
5,134
3,581
4,082
4,555
5,359
3,681
4,311
4,684
5,473
2,397
2,624
3,125
3,408
4,027
2,445
2,659
3,139
3,435
4,225
2,454
2,675
3,145
3,493
4,395
2,454
2,738
3,220
3,501
5,020
2,484
2,740
3,223
3,537
1,944
2,240
2,435
2,648
3,012
1,958
2,253
2,478
2,684
3,067
1,966
2,270
2,490
2,697
3,084
1,997
2,272
2,511
2,726
3,090
2,086
2,274
2,514
2,770
3,096
1,764
2,055
2,299
2,514
2,785
1,807
2,058
2,308
2,558
2,785
1,812
2,088
2,335
2,577
3,02
1,84
2,125
2,349
2,593
3,042
1,852
2,162
2,356
2,601
3,116
Tabela 4.2
2,526
2,834
3,060
2,485
2,793
3,039
2,110
2,710
2,870
1,889
2,485
2,854
(a) Length 20mm
2,546 2,628 2,628 2,669 2,669 2,710 2,731 2,731 2,731 2,752 2,752
2,834 2,854 2,875 2,875 2,895 2,916 2,916 2,957 2,977 2,998 3,060
3,080
(b) Length 50mm
2,526 2,546 2,546 2,567 2,628 2,649 2,669 2,710 2,731 2,752 2,772
2,813 2,813 2,854 2,854 2,854 2,895 2,916 2,936 2,936 2,957 2,957
3,039 3,039 3,080
(c) Length 150mm
2,260 2,340 2,440 2,510 2,510 2,570 2,570 2,610 2,610 2,610 2,650
2,710 2,710 2,750 2,750 2,750 2,750 2,770 2,770 2,790 2,830 2,830
2,870 2,900 2,900 2,920 2,940
(d) Length 300mm
2,115 2,177 2,259 2,279 2,320 2,341 2,341 2,382 2,382 2,402 2,443
2,505 2,505 2,526 2,587 2,608 2,649 2,669 2,690 2,690 2,710 2,751
2,854 2,875
2,793
3,060
2,793
3,018
2,670
2,830
2,464
2,751
Sob ALM, a distribuição de logT:
1) Variam em locação de um pro outro.
2) Tem a mesma variância.
Variância amostral dos dados da tabela 4.1são:
0.042, 0.040, 0.045 e 0.037 => estaticamente equivalentes
figura 4.1 é o gráfico 2 da seção 4.2: variação horizontal.
Figuras 4.2e 4.3: linear (particularmente na figura 4.3).
Gráficos de log( log(Sˆ j )) contra log t deve ser:
•Cópias verticais da mesma curva, sob PHM;
•Linhas paralelas, sob modelo de regressão Weibull (ver
seção 4.7)
log( log(Sˆ j )) contra log t
figura 4.4 mostra gráfico de
para dados da tabela 4.1com linhas
Weibull
(ajuste
por
Máxima
Verossimilhança
separadamente em cada amostra).
Figura 4.6 mostra o mesmo gráfico para os dados da
tabela 4.2
Sob POM, devem ser cópias verticais uma da outra (ver
figura 4.5)
Agora, nesta Seção, vamos nos concentrar em Weakest
Link Hypotesis:
WLH : SL (t )  [S1 (t )]L
onde Sr é a função de sobrevivência para fibra de
comprimento r.
WLH pode ser implementada nos seguintes modelos
Weibull:
1



 t
M 0  S L (t )  exp L 
t

  1



 t 
  exp  
t



  L 


 onde t
L



 t1L 

t L  t1L  ,0    1
M1 – como em M0 mas com log  log1   log L
M2 – como em M1 mas com
(modelo log-linear no parâmetro escala da distribuição
Gumbel);
M3 – distribuição
comprimento.
Weibull
separada
para
cada
A função de verossimilhança pode ser obtida para cada
modelo e na maioria dos casos, maximizá-la não é difícil.
Os valores maximizados da função log-verossimilhaça
são:
l0 = -229.1 (2 parâmetros)
l1 = -227.7 (3 parâmetros)
l2 = -227.6 (4 parâmetros)
l3 = -220.1 (8 parâmetros)
Teste de hipótese H1: =1 é calculado via estatística de
razão de verossimilhança:
W1  2(l1  l0 )  2x1.4  2.8  c02.05 (3  2)  3.84
 H1
é aceita.
(O Estimador de Máxima Verossimilhança sob M1 é 0.90)
Dados da tabela 4.2: valores maximizados da logverossimilhança são:
l0 = 21.8 (2 parâmetros)
l1 = 29.6 (3 parâmetros)
l2 = 31.5 (4 parâmetros)
l3 = 35.8 (8 parâmetros)
Teste de H1 agora tem W1 = 15.6 > 3.84 => M0 rejeitado
O Estimador de Máxima Verossimilhança de  sob M1 é 0.58
WLH não necessariamente implica modelo Weibull
PHM também é possível.
4.9 Análise Bayesiana de Dados de Confiança
Como mencionado anteriormente, as análises de grande
amostra de máximo verossimilhança oferecem resultados
similares à análise bayesiana com vaga informação a priori.
Exemplo: A estimativa de máximo verossimilhança 0,58 de ξ
para os dados da tabela 4.2 é uma aproximação da média da
posteriori de ξ.
Os testes de níveis
complemento bayesiano.
de
significância
não
possuem
O teste de hipótese bayesiano pode ser feito da seguinte
maneira:
1) Constrói-se a região de máxima densidade a posteriori
100(1-α)% do parâmetro de interesse (ξ, no exemplo
anteriormente citado);
2) Observa-se se esta região contém o valor do parâmetro
especificado por H (1, no exemplo anterior);
3) Aceita-se H se este for o caso ou rejeita-se H, caso contrário.
Também, para grandes amostras, a razão de máxima
posteriori pode ser aproximada pela razão de máxima
verossimilhança, se a priori for vaga.
Exemplo: Considere os dados da tabela 6.1 do modelo
Weibull.
log t     (temperatura  273.2)1   1 / 2W onde W ~ Gumbel
5. TAXA DE FALHA PROPORCIONAL
5.1 Introdução
Baseada no documento de Cox (1972).
Concentrada em análise de sobrevivência em vez de
confiança.
O modelo linear de taxa de falha proporcional assume
que
h (t ; x )  h (t ) exp{ x  }
T
0
Cada variável x afeta a taxa de falha, oscilando para cima
e para baixo.
A análise depende das pressuposições dos parâmetros
feitos em h0.
Se a taxa de falha de base acumulada , H 0 (t ; )   [G (t )]
então:
S (t; x)  exp{H0 (t ) exp( xT  )}
 exp{[G(t )] exp( xT  )}
 exp{ exp(y  log   xT  )}
onde y=logG(t).
Segue que z = ηy + log α + xT β ~ Gumbel e
y = -η-1 log α – xT(-η-1 β) + η-1z
sendo que esta é a fórmula do modelo log-linear de
regressão da Weibull.
Pode-se usar outras fórmulas para h0 e a que possuir o
melhor ajuste, pode então, ser utilizada.
5.2 Análise do Modelo Semiparamétrico de
Taxa de Falha Proporcional
De agora em diante, h0 não será especificado.
Assume-se n itens, r tempos de falha distintos e Ri é o
conjunto de risco, ou seja, o grupo de itens que tem
chances de falhar um pouco antes de t(i).
β é estimado pela função de verossimilhança:
r
L(  )  
i 1
exp( x(Ti )  )
 exp( x
l Ri
T
(l )
)
Há muitas justificativas para a verossimilhança acima:
1) Cox (1972) originalmente considerou análises condicionais
dos tempos de falha . A probabilidade condicional do item (i)
falhar dado que conhecemos as falhas anteriores é:
h(t(i ) ; x(i ) )
 h(t(l ) ; x(l ) )
l  Ri

exp( x(Ti )  )
 exp( x
l  Ri
T
(l )
)
A verossimilhança é formada pelo produto das probabilidades
condicionais.
r
(Basicamente, usa-se Pr( A1 ,..., Ar )   Pr( Ai Ai 1 ,..., A1 ) )
i 1
2) Também, é obtida por Kalbfleisch & Prentice (1973) como a
verossimilhança
marginal de β baseada nos postos das observações.
Se L(β ) é tratado como a verossimilhança, então a teoria da
máximo verossimilhança pode ser aplicada:
i) O estimador de máximo verossimilhança
maximizando L(β );
é obtido
ii) Sua variância assintótica é obtida pela 2ª derivada de l(β).
Computacionalmente:
Primeiro, calcula-se l  exp( xT  );
Então, acumula-se o decréscimo dos tempos de falha
pra obter ci  l , i = 1,...r.
lRi
L(β ) é formado pelo produto dos termos λ(i)/ci.
O log de verossimilhança e suas derivadas são:
r
l (  )   ( x(Ti )   log ci )
i 1
r
dl (  )
  ( x( i )  vi )
d
i 1
r
d 2l (  )
T


(
A

v
v
)

i
i
i
2
d
i 1
Onde
vi 
1
ci
 l xl
l  Ri
e
Ai 
1
l xl xlT

ci lRi
^
A maximização numérica é usada para obter  .
Há pacotes computacionais que podem nos auxiliar, por
exemplo, o R.
Observações Empatadas:
Dados empatados complicam os cálculos da
verossimilhança.
Uma aproximação razoável neste caso é:
exp(siT  )
L(  )  
di
T
i 1 [  exp( x( l )  )]
r
l Ri
Onde di é o número de tempos que falham em t(i) e s(i) é
a soma das covariâncias que dependem destes tempos.
A aproximação é boa se di /n(i) for pequeno.
As covariâncias que dependem dos tempos são
substituídas de x(l) por x(l) (t(i)).
5.3 Estimação de Sobrevivência e Função de
Taxa de Falha
A taxa de sobrevivência de base é obtida por um
método não paramétrico junto com o estimador do
produto limite, como:
^

j :t ( j )  t
(1 
j
^
^

l  Ri
)
onde
^
l  exp( xl  )
l
^
Também, o estimador de H0 é dado por
^
H 0 (t )   log S0 (t )
Os métodos do capítulo 3 podem ser aplicados aqui
para verificar as formas paramétricas para S0 .
.
5.4 Métodos de Verificação
As suposições da taxa de falha proporcional podem ser
verificadas pelos dados e modelos de ajuste:
hj (t; x)  h0 j (t ) exp( xT  )
Depois do ajuste, o gráfico do log vs. t deve ter paralelas
verticais.
O resíduo pode ser definido por um item não
censurado como:
^
^
e  H 0 (t ) exp( x  )
T
e por um item censurado, e deve ser somado por 1.
Se os resíduos se comportam como amostras de uma
distribuição exponencial de média 1, o ajuste é adequado.
Deve-se ter cautela, pois os resíduos podem estar
adequados mesmo quando a taxa de falha proporcional
não estiver.
Também, o gráfico dos resíduos vs. as variáveis pode
não se encaixar no modelo.
5.5 Exemplos Numéricos
Tabela 4.3 contem dados de falha sob stress
Figura 4.7 indica a influencia da taxa de stress na falha
h(t;x) = h0(t)em(x) para alguma função m de x
caso m(x)= -β log x ,
temos que β^ =0,2069 com erro padrão assintótico =
0,0457 => forte efeito da taxa de stress
Figura 5.1 mostra um gráfico de log^H0(t) vs log t
indicando Weibull
A checagem da hipótese de modelo de taxas
proporcionais é feita segundo estratificação sugerida em
§5.4
Isto leva ao gráfico na Figura 5.2 (paralelismo vertical=
OK)
Checagem da linearidade em log x é feita através de
um estimador não paramétrico dos efeitos de x
5
m x     j z j  x 
i 1
onde zj(x) é um indicador de que x está no nível j.
Este modelo não é estimável (uma constante pode sempre
ser adicionada ao αj)
Mas se torna estimável ao se impor uma restrição (Σαj
=0 ou αk =0, para algum k)
O valor maximizado de l é –177,00 em oposição a178,19 sob o modelo linear.
W = 2 (-177,00+178,19) = 2,38 < 7,81 = χ20,05(4-1) =>
modelo linear aceito
Tabela 4.4 contem dados de tempo de falha de fibras
Figura 4.8 contem PH plots para diferentes níveis de
stress.
Dando indicação de que o efeito dos lotes nos tempos
de falha do modelo proposto é
7


ht; x   h0 t exp  log S   j z j 
i 1


onde S é o nível de stress e zj é o indicador do lote j , j =
1,..,7
Tabela 4.3 Teste de força em bilhas de liga de Tungstênio
Coluna 1: Taxa de Stress em MNm-2seg-1 ; Colunas 2-13: Tempo de falha
0,1
1676 2213 2283 2297 2320 2412 2491 2527 2599 2693 2804
1
1895 1908 2178 2299 2381 2422 2441 2458 2476 2528 2560
10
2271 2357 2458 2536 2705 2783 2790 2827 2837 2875 2887
100
1997 2068 2076 2325 2384 2752 2799 2845 2899 2922 3098
1000
2540 2544 2606 2690 2863 3007 3024 3068 3126 3156 3176
2861
2970
2899
3162
3685
1
- x
- x
-x
o
-x
o
xo
o
o
0
*
-1
+
-2
x
x
*
-
x
*
*
* x
x
*
*
-
*
*
*
-
*
o
o
o
-3
Log(log(probabilidade de sobrevivência))
Figura 4.7 Plot da Weibull para dados de tungestênio,
com linhas ajustadas pela Weibull
7.4
7.5
7.6
7.7
7.8
7.9
8.0
Log do tempo de falha sob as taxas de stress 0.1(*),1(+),10(-),100(x),1000(o)
2
1
0
-1
+
-3
-2
x
x
*
x
* --x xoo
* -x
o
-x
o
* *
x-x
oo
* *
x
o
*
*x - - o o
*
o
*
o
-4
log(Função de taxa de falha acumulada)
Figura 5.1
7.4
7.6
7.8
8.0
8.2
Log do tempo de falha sob as taxas de stress 0.1(*),1(+),10(-),100(x),1000(o)
3
2
1
0
-4 -3 -2 -1
log(Função de taxa de falha acumulada)
Figura 5.2
x
x
x
xx
x
x
x
xxxxx
x
x
x
x
x
xxx
xx xx x
x
x
xx
xxxx
x
xxx
x
x
x
x x
x
x
7.6
7.8
8.0
8.2
log(tempo de falha)
Procedimentos iterativos funcionarão melhor caso
iniciem com bons valores iniciais.
Estes podem ser obtidos através da MV com Weibull
ou por mínimos quadrados.
γ^ =32,586 com erro padrão assintótico igual a 3,097 =>
forte efeito da taxa de stress.
O teste para efeito do lote é baseado na Razão de MV
entre M0: não há efeito do lote e M1: existe efeito do lote
A maximização das log verossimilhanças são –320,119
e –261,775 que nos trazem
W =2(-261,775+320119)=116,69 >> χ2α(7)=> o efeito do
lote é significante.
A hipótese de modelo de taxas proporcionais é testada
pela estratificação dos dados através dos lotes.
A figura 5.3 mostra o gráfico do log^Hj aproximado por
PH
A figura 5.4 mostra o gráfico do log^H0 vs log t , onde
linearidade indica Weibull
Também é possível investigar a interação entre o stress
e os lotes.
Isto é feito através da inclusão de 7 covariáveis zj log
S no modelo
A maximização da log-verossimilhança é –259,241 com
W=5,07 (não significante)
Ajustando uma linha reta na figura 5.4 teremos
log^H0(t)= -11,19 + 1,533logt
E
H^j(T)=15,296T1,533S32,586eα^j , j = 1,...,8 (α8 = 0)
Exemplo : o tempo médio de falha previsto para S= 23,4
e j = 4 é T = 7,4 anos.
Tabela 4.4 Dados de Falha de fibra de Kevlar
Stress em Mpa; Numero do lote; Tempo de falha em horas(* caso dado censurado)
Stress Lote Tempo
Stress Lote Tempo
Stress Lote Tempo
Stress Lote Tempo
29,7 2
2,2 29,7 5
243,9 27,6 2
694,1 25,5 1 11.487,3
29,7 7
4,0 29,7 4
254,1 27,6 4
876,7 25,5 5 11.727,1
29,7 7
4,0 29,7 1
444,4 27,6 1
930,4 25,5 4 13.501,3
29,7 7
4,6 29,7 8
590,4 27,6 6
1.254,9 25,5 1 14.032,0
29,7 7
6,1 29,7 8
638,2 27,6 4
1.275,6 25,5 4 29.808,0
29,7 6
6,7 29,7 1
755,2 27,6 4
1.536,8 25,5 1 31.008,0
29,7 7
7,9 29,7 1
952,2 27,6 1
1.755,5 23,4 7 4.000,0
29,7 5
8,3 29,7 1 1.108,2 27,6 8
2.046,2 23,4 7 5.376,0
29,7 2
8,5 29,7 4 1.148,5 27,6 4
6.177,5 23,4 6 7.320,0
29,7 2
9,1 29,7 4 1.569,3 25,5 6
225,2 23,4 3 8.616,0
29,7 2
10,2 29,7 4 1.750,6 25,5 7
503,6 23,4 5 9.120,0
29,7 3
12,5 29,7 4 1.802,1 25,5 3
1.087,7 23,4 2 14.400,0
29,7 5
13,3 27,6 3
19,1 25,5 2
1.134,3 23,4 6 16.104,0
29,7 7
14,0 27,6 3
24,3 25,5 2
1.824,3 23,4 5 20.231,0
29,7 3
14,6 27,6 3
69,8 25,5 2
1.920,1 23,4 6 20.233,0
29,7 6
15,0 27,6 2
71,2 25,5 2
2.383,0 23,4 5 35.880,0
29,7 3
18,7 27,6 3
136,0 25,5 3
2.442,5 23,4 1 41.000,0
29,7 2
22,1 27,6 2
199,1 25,5 8
2.974,6 23,4 1 41.000,0
29,7 7
45,9 27,6 2
403,7 25,5 2
3.708,9 23,4 1 41.000,0
29,7 2
55,4 27,6 2
432,2 25,5 8
4.908,9 23,4 1 41.000,0
29,7 7
61,2 27,6 1
453,4 25,5 2
5.556,0 23,4 4 41.000,0
29,7 5
87,5 27,6 2
514,1 25,5 6
6.271,1 23,4 4 41.000,0
29,7 8
98,2 27,6 6
514,2 25,5 8
7.332,0 23,4 4 41.000,0
29,7 3
101,0 27,6 6
541,6 25,5 8
7.918,7 23,4 4 41.000,0
29,7 2
111,4 27,6 2
544,9 25,5 6
7.996,0 23,4 8 41.000,0
29,7 6
144,0 27,6 8
554,2 25,5 8
9.240,3 23,4 8 41.000,0
29,7 2
158,0 27,6 1
664,5 25,5 8
9.973,0 23,4 8 41.000,0
0
1
Log(-log(probabilidade de sobrevivência))
Figura 4.8 plot da Weibull para dados de fibra de Kevlar
sob pressão , com linhas ajustadas pela Weibull
+
-5
-4 -3 -2
-1
+++
+
+
+
+
+
+
+
+
++
+
+
-
+++
x x
x
+
+
x
+
+
x
+ + - -+
xxx
x
++++
-x
+
x
+
x
++
--o o
xx
-xx
oo
x
o
x
o
o
- x x
o
x
o
x
o
x
2
4
6
8
10
12
Log do tempo de falha nos níveis de stress 29.7(+),27.6(-),25.5(x),23.4(o)
4
2
0
-2
7
7
7 3
3
5 32
2 3
6
52
62 3
7
7
3 27 2
6
6 85 5
7
3 22 2
41
2
8 4
2
66
1
85
2
2
4
1
6
8
2
1
4
26
6
8
11 4
448 8
1
8
811 44
18 4
4 1
7 73 3 3
5
5 2 2 26
22
6
8
-4
2
-6
log(Função de taxa de falha acumulada)
Figura 5.3 Log da função de taxa de falha acumulada
para dados da fibra de Kevlar: modelo estratificado .
Lotes 1-8
0
2
4
6
8
log(tempo de falha)
10
12
5
0
-5
o
o
o
o
o
oo
oo
oo
ooo
o
oo o
oo
o o oo o o
o oo o
o o oo
o oooo o
oo oo o
o o
oooooo o
ooooooo
o
oo
o oooo
o oooo oo
o ooo o
o
oo o
o
-10
log(Função de taxa de falha acumulada)
Figura 5.4 Log da Log da função de taxa de falha
acumulada para dados da fibra de Kevlar
0
2
4
6
8
log(tempo de falha)
Em geral:
i = 0, 1, ..., n-2
10
12
6. MODELOS BAYESIANOS DINÂMICOS
6.1 Introdução
Os modelos PH (Proportional Harzard Model) são razoáveis
em alguns casos, porem muitas as vezes modelos mais
genéricos são requeridos.
Os enfoques utilizados são :
1 - Parametrização das mudanças entre as faixas de
   0ekt
estimação. Exemplo :
2 - Estimação independente por partes em cada
intervalo.
O primeiro enfoque exige uma suposição forte em torno
dos dados enquanto o segundo não exige nada ,porem as
informações contidas nos intervalos anteriores são
desconsideradas.
Modelos dinâmicos permite mapearmos as mudanças entre
as diferentes faixas sem a necessidade de fortes suposições e
permitindo que os dados orientem a direção do modelo.
Modelos dinâmicos provem uma abordagem para os
modelos de PHM.
6.2 Elementos do Modelo
Distribuição Exponencial por Partes
A Distribuição Exponencial por Partes é usada para
tempo de falha baseada na função de densidade Harzard :
h(t; x)  i  e x i, t  Ii , i  1,2.. onde o primeiro componente
t
Bi  Bi,0 (baseline) e o primeiro componente x  1 .
T
Fatorização Temporal da verossimilhança
Os dados podem ser divididos ao longo do mesmo
intervalo Ii . Cada intervalo prove a verossimilhança
baseado na condicional das informações obtidos do
intervalo anterior.
N
TFL : L   Li
I 1
onde
N
N
I 1
I 1
Li   p(t w | , Di 1 ) S (tiq | , Di 1 )
Esta é a verossimilhança do intervalo Ii condicionado
em Di-1,
F(Ii) é o conjunto de itens que falharam no intervalo Ii ,
R(Ii) é o conjunto de itens que estavam em risco antes Ii
mais não falharam em Ii e
L é o conjunto de parâmetros
6.3 Desenvolvendo a Análise
Evolução Paramétrica
Como descrito em 4.6, parâmetros sucessivos podem
ser calculados via
 Gi (  )   w
i
i
i 1
onde,
E ( wi )  0
Var ( wi )  Wi
bi  tamanho ( I i )
e
Wi  biW
i
O valor de W indica o potencial de mudança dos b´s
por unidade de tempo.
W  0 conduz ao um modelo estático.
W-1  0 grande potencial de mudança conduzindo a
modelos independentes.
Assume-se que :  i | Di ~ [mi 1 , Ci 1 ] (especificação parcial)
A priori para o intervalo Ii é  i | Di 1 ~ [ i , Pi ] onde
 i  Gi (bi )mi 1
Pi  G (bi )Ci 1GiT (bi )  Wi
Novamente, o valor de Wi (ou W) podem ser definidos
como fatores de desconto.
Verossimilhança e Atualização
A verossimilhança Li pode ser escrita como :
Li 
 i
wF ( I i )
( w) (  i ( w) (t w ti 1 ))
e
e
(  i ( q ) (tiq ti 1 ))
qR ( I i )
ri
  Lij
j 1
Onde
Li 
 i
( w) (  i ( w) (t w ti 1 ))
e
wF ( I i )
Lij  [i
e
(  i ( q ) (tiq ti 1 ))
qR ( I i )
ri
  Lij
j 1
( j)
( j )  ij (  i (tij ti 1 ))
] e
xs (T ) 
Onde i  e
é o valor constante da hazard para o
item s em Ii , ri é o numero do item em Ii
dij é o indicador de falha para o item em j ,
tij é o tempo de sobrevivência do indivíduo j em Ii .
(s)
Uma analise baysiana completa envolve a especificação
da prior para 1 | Do até a integração numérica de cada
intervalo Ii , alternativamente, especificações parciais
podem ser utilizadas.
A atualização é efetuada em ciclos através dos intervalos
e dentro de cada intervalo através de seus itens. O calculo
envolvendo esta analise esta especificado abaixo:
 i | Di 1, j 1
~
G.R. : log (i j )  Z j  i
LINEA
R
BAYES
~
(i j ) | Di 1, j 1
lij
(i j ) | Di 1, j
 i | (i j ) , Di 1, j
~
i
MEDIA
SAIDA
 i | Di 1, j
~
j  j 1
Passo 1 : Dado
 i | Di 1; j 1 ~ [ ij , Cij ]
obtemos :
log (i j ) | Di 1, j 1 ~ [ f ij , qij ]
onde,
f ij  x j T aij
qij  x j T Cij x j ;
(i j ) | Di 1, j 1
Passo 2 : Aproximar
por uma Gamma
utilizando os dados do Passo 1. Parâmetros
e ij  (q )1 e( fij )
;
 ij  (qij ) 1
ij
Passo 3 : Atualizar, usando a informação do Lij para
;
(i j ) | Di 1, j ~ G[ ij   ij ,  ij  tij  ti 1 ]
Passo 4 : Use os métodos lineares de bayes e a relação
( j)
entre log i e  i para obter
 i | Di 1, j ~ [mij , Cij ];
~
Passo 5 : Defina ai,j+1 = mij e Pi,j+1=Cij e repita os
passo anteriores;
Passo 6 : Quando j=rij , todos os itens Li tem sido
processado. Evolui para o próximo intervalo.
O procedimento acima tem como
especificação parcial da distribuição :
 i | Di , i  1,2,...;
resultado
a
~
Este procedimento continua ate In , o ultimo intervalo
com informações.
Smoothing
É mais conveniente trabalhar com a distribuição de
 i | DN , i  1,.., N .
~
Ela gera resultados mais precisos e é utilizada na
predição.Existem distribuições Smoothing obtidas via
algoritmos recursivos.
 i | Dn ~ [mi N , Ci N ]
~
então
miN  mi  Ci GiT1 (bi ) Pi 11 (miN1  ai 1 )
CiN  Ci  Ci GiT1 (bi ) Pi 11 ( Pi 1  CiN1 ) Pi 11GiT1 (bi )Ci
N
N
C
 CN
N
m

m
Com valores iniciais N
N
6.4 Predição de Ocorrências Futuras
Dois casos devem ser considerados :
1) Predição dos itens que sobraram e ainda não falharam
2) Predição de novos itens.
Para o caso 2, a distribuição de sobrevivência predita é
S(W\DN) onde ,
i 1
S (W | DN )  S ( w | W  ti 1 , DN ) S (t j | W  t j 1 , DN )
j 1
com
w  Ii
O fato acima pode ser obtido da seguinte forma :

S ( w | W  ti 1 , DN )   S ( w | k ,W  tk 1 , DN ) p(k | W  tk 1 , DN )dk
0

  e( k ( wtk 1 )) p(k | W  tk 1 , DN )dk
0
p (k | W  t k 1 , DN ) Obtido de p(  k | W  t k 1 , DN )
Aproximados pela distribuição de  k | DN
.
Se
k | DN ~ G ( k ,  k ) então,
S ( w | DN )  (1 
w  ti
i
)
 i
i 1
 (1 
j 1
bi
i
)  i
w  Ii
Cálculos parecidos são utilizados para a distribuição
preditora a posteriori S ( w | D0 ) .
Devemos verificar se a especificação da priori e razoável.
6.5 Ilustrações
Exemplo 1 : Dados de 90 pacientes com câncer gástrico
divididos em dois grupos de tratamento : Quimioterapia e
Radioterapia / somente quimioterapia.
PHM ajustado com 2 covariáveis : x1 - Indicador de
tratamento e x2=x1*t. Isto é equivalente a assumir um
modelo de regressão com os coeficientes  1  1  2t
^
^
Os estimadores MLE são : 1  1.2711  2  0.0794
.
^
^
Os estimadores MLE são : 1  1.2711  2  0.0794
Analises dinâmicas foram feitas com :
3
1 - Priori Vaga 1 | D0 ~ [0,10 I 2 ];
2 - dados obtidos pelos tempos observados de falha.
3 - evolução do erro da variância , dado por :
0 
0
Wi  bi

 0 0.0016 
Para a variação
com em 4.6 .
1
foi utilizado um crescimento linear
Figura 1 mostra as estimativas do coeficiente de
regressão.
Figura 2 mostra a curva de sobrevivência predita.
Figura 3 mostra estimativas para as taxas de mudança de
1 .
Exemplo 2 : Dados de 300 desempregados que foram
acompanhados durante um ano na Inglaterra. Esta é uma
amostra de um grande estudo feito por Germerman e
West (1987). As covariáveis selecionadas são x1 ~ idade
e x2 ~ faixa de renda.
Analises dinâmicas foram feitas com :
T
2
2
1 - Priori 1 | D0 ~ [(1,0.025,0.25) , diag (1, (0.0125) , (0.025) )];
2 - Intervalos Semanais;
3 - Evolução do erro da variância Wi=7*diag(10-2,105,10-3).
Figura 4
indivíduos.
mostra preposteriori par 4 hipotéticos
Figura 5 mostra estimativas paramétricas(priori estava
muito forte).
Figura 6 mostra predição para os mesmos 4 indivíduos.