Motivação

No capitulo anterior, diferentes distribuições de probabilidade foram discutidas para o uso na modelagem de dados de confiabilidade, para quais a estimativa dos seus parâmetros estava associada a resolução do método da máxima verossimilhança. Dentre as considerações decisivas para a correta aplicação do método reside o fato dos dados coletados para a modelagem estarem associados a um mesmo fenômeno. Vejamos a razão para a supracitada condição no exemplo 1.

Exemplo 1

Dados referentes ao tempo até a falha (em meses) de duas bombas centrífugas são simulados para qual se deseja estimar os estimadores de máxima verossimilhança considerando uma distribuição exponencial. Não há dados censurados.

exp1 = rexp(n = 19, rate = 1)
exp2 = rexp(n = 19, rate = 0.2)
exp3 = c(exp1, exp2)

No caso da ausência de censura, o estimador de máxima verossimilhança (\(\lambda\)) para a distribuição exponencial pode ser estimado pela Equação 1.

\[ \lambda = \frac{n}{\sum(x_i)} \]

Nesta, \(n\) corresponde ao número de observações e \(x_i\) ao i-nesimo tempo até a falha. Em nosso exemplo \(\lambda\) assume o valor de 0.29, através do qual pode-se obter a curva parametrizada apresentada na Figura 1.

Modelo Exponencial para os tempos de falha das bombas 1 e 2

Modelo Exponencial para os tempos de falha das bombas 1 e 2

Neste ponto, pode ser curioso compreender se existem diferenças significativas entre as duas bombas capazes de tornar a modelagem conjunta do tempo até falha incorreta. Em uma breve abordagem gráfica, podemos plotar os dados de ambas as bombas em diferentes paineis e, conforme a Figura 2, verificar que o comportamento quanto ao tempo de falha das bombas se diferencia fortemente entre as mesmas, tornando infundada a análise dos tempos conjuntamente.

Curvas de Densidade de falha para as Bombas 1 e 2

Curvas de Densidade de falha para as Bombas 1 e 2

Duas abordagens partem desta conclusão: (1) a construção de modelos individuais em cada bomba, ou (2) a incorporação de variáveis/propriedades na modelagem de forma a explicar a diferença entre os tipos de bombas. Nesta última opção não estamos apenas aumentando a complexidade do modelo (i.e. sua flexibilidade) como também permitindo a realização de inferências quanto a que aspectos impactam o tempo até a falha do sistema estudado e como estes impactam.

Práticas

Prática 1 - Seleção do modelo com melhor adequação aos dados

Contextualização A seleção do modelo de regressão mais representativo para um conjunto de dados de confiablidade passa não somente pela avaliação da qualidade do ajuste das distribuições selecionadas, como também pela avaliação das considerações associadas a cada modelo. Diferentes distribuições podem ser avaliadas para tal proposito, contudo em todas o ponto comum estar em regredir apenas o parâmetro de escala como uma função de covariáveis interessantes para compreensão do sistema. Especificamente, a distribuição exponencial e a Weibull assumem a seguinte forma em seus modelos de regressão:

\[ R(t|x) = exp \left \{-\frac{t}{exp(X\beta)}\right \}\]

\[ R(t|x) = exp \left \{-(\frac{t}{exp(X\beta)})^{1/\sigma}\right \}\]

Objetivos Instrumentalizar os alunos no ferramentário necessário para seleção da distribuição de proabilidade mais adequada na construção de um modelo de regressão de análise de confiabilidade.

Banco de dados Consideremos os dados apresentados na Tabela 1 referentes ao tempo de falha de transmissores de pressão em uma indústria química associado a classe de pressão dos transmissores (SINTEF, 2010)

Proposta de Desenvolvimento

## Warning: package 'knitr' was built under R version 3.4.3
Dados de tempo-até-falha de bombas centrífugas
time cens pres
1.9188729 1 15
0.3518219 0 15
0.7451203 1 15
0.6853860 1 15
1.2972413 1 15
0.3802726 1 15
1.7514308 0 15
0.5983319 1 15
0.1489464 0 15
0.0619509 1 15
0.1348417 1 15
0.6985645 0 15
0.1870072 1 15
1.9137282 1 15
0.9425777 1 15
0.9392172 1 50
0.3457727 0 50
0.1748592 1 50
1.2915819 0 50
0.4882132 1 50
0.2467309 1 50
0.1637658 0 50
0.3002241 0 50
0.0541954 1 50
0.3883420 1 50
0.3669609 1 50
0.5335229 1 50
0.6270733 1 50
0.1070003 0 50
1.4513119 0 50
0.0724456 0 70
0.4359634 1 70
2.2555652 0 70
1.2354384 1 70
0.6180993 1 70
0.0213377 0 70
5.2453314 1 70
1.1406079 0 70
0.0232129 1 70
0.6566084 0 70
1.0657155 0 70
1.8496329 1 70
0.1442509 1 70
1.0934230 1 70
0.5169809 0 70
0.4062703 1 70
1.2748278 1 70
2.1892359 1 70
1.3939613 1 70
0.3647871 1 70

Iremos supor, para fins de simplicação, as distribuições exponencial e weibull em nossas análises. Inicialmente, avaliando se ambas, quando linearizadas, se comportam em torno de uma reta suporte padrão.

Uma vez que ambas as distribuições se comportam de forma linear quando sua resposta é linearizada, podemos manter ambas em nossa proposta de modelagem. A construções do modelo de regressão vem com o uso da função survreg e associação dos dados ekm_data com a nossa variável (classe de pressão).

Para a seleção da distribuição de probabilidade a ser utilizada, pode se aplicar o teste de razão de verossimilhanças, o qual consiste na razão entre a logaritmo da verossimilhança de um modelo mais simples (que se deseja testar) e outro mais generico (Equação 1). Uma vez que esta razão segue uma distribuição qui-quadrada, pode ter seu p-valor calculado por está distribuição considerando o número de graus de liberdade como a diferença do número de parâmetros.

\[TRV = 2\cdot log(\frac{L(\theta)_{simples}}{L(\theta)_{Generalizado}})\]

Para um p-valor de 0.5642792, o teste não encontra diferenças significativas entre o uso da distribuição exponencial e Weibull, de forma que seguiremos com a distribuição exponencial por conta do principio da parcimônia. O fato do parâmetro de forma da Weibull se aproximar de 1, confirma nossa análise - vide análise a seguir.

## 
## Call:
## survreg(formula = ekm_data ~ pratical_data1$pres, dist = "weibull")
##                       Value Std. Error      z     p
## (Intercept)         -0.1373    0.33876 -0.405 0.685
## pratical_data1$pres  0.0075    0.00645  1.163 0.245
## Log(scale)          -0.0779    0.13240 -0.588 0.556
## 
## Scale= 0.925 
## 
## Weibull distribution
## Loglik(model)= -39.7   Loglik(intercept only)= -40.3
##  Chisq= 1.28 on 1 degrees of freedom, p= 0.26 
## Number of Newton-Raphson Iterations: 6 
## n= 50

Para avaliarmos o impacto da pressão como variável no modelo exponencia, criaremos 3 curvas referentes à pressão nas classes 15, 50, e 70 kgf/cm².

A análise da Figura nos revela que transmissores de pressão de classes de pressão mais elevadas possuem maior tempo médio até falha.

Em uma última análise, semelhante a regressão linear clássica, avaliamos os resíduos do modelo proposto. Dentres os diferentes tipos de resíduos apresentados na literatura, se destacam os resíduos de cox-snell (GIOLO, 2006). Estes são definidos para a distribuição exponencial e Weibull como:

\[e_i = [t_i exp(-X\beta)]^{\gamma}\]

Onde \(\gamma\) corresponde ao parâmetro de forma (unitário para a exponencial) e \(\beta\) corresponde ao vetor de parâmetros do modelo de regressão. Estes residuos quando plotados contra a exponencial da resposta do estimador Kaplan-Meier para os resíduos devem seguir uma reta suporte padrão caso sigam a distribuição exponencial esperada para estes resíduos, o que pode ser verificado pelo gráfico resultante a seguir.

Prática 2 - Seleção de variáveis em um modelo de regressão para confiablidade

Contextualização A regressão de dados como entendemos é um problema de natureza multivariada, em sua forma mais simples a estimativa de parâmetros passa pelo seu termo constante e uma covariável. Ao elevarmos o número de parâmetros de nosso modelo, surge o questionamento para compreender quais destes são úteis para compreensão do fenômeno e que permitem manter a parcimônia na construção do modelo.

Na análise de confiabilidade está questão se mantém, com o agravante que a interpretação do impacto de cada variável passa a ser não-linear, pela transformação \(y = log(t)\).

A parcimônia guia a seleção de variáveis (GIOLO, 2006). Para após obtermos o logaritmo da varossimilhança (loglik) do modelo sem covariáveis - modelo mais simples - comparamos o loglik de modelos com 1 variável individualmente e os comparamos através do teste de razão de verossimilhanças (sesguindo uma distribuição qui-quadrado com graus de liberadade igual a diferença no número de parâmetros dos modelos). Para as variáveis mantidas na primeira etapa, novos modelos são criados e comparados com aquelas mais simples pelo mesmo teste.

Objetivos A presente pratica busca fornecer aos alunos ferramentário necessário para a tomada de decisão na seleção de modelos de regressão para análise de confiabilidade.

Banco de dados Dados provenientes de Dantas et al. (2009) foram adaptados de forma a criar um banco de dados dinâmico contendo dados de tempo-até-falha para poços de petróleo considerando como covariáveis para modelagem, a idade do poço (variável contínua) e o método de elevação (variação discreta - 1 para Bombeio Mecânica e 0 para Bombeio por cavidade progressiva). Os dados apresentam censuras a direita.

Dados de tempo-até-falha de poços de petróleo
time cens idade metodo
1.1858030 1 5.457808 1
1.0862391 0 2.536018 1
2.0431615 1 4.388090 1
1.0817025 1 4.570102 1
2.0023768 1 5.692060 1
1.8529707 0 4.409240 1
1.6663323 1 3.713536 1
1.0448526 0 3.982563 1
1.2232016 1 5.051435 1
2.4908416 0 4.012460 1
1.9967279 1 2.177059 1
1.4972317 0 5.093868 1
1.6008150 1 5.111796 1
1.4582668 1 2.971650 1
2.0487817 1 4.210862 1
2.6018241 1 5.656920 1
1.3913145 1 3.522883 1
1.7826698 1 3.440202 1
2.8705286 1 3.930378 1
2.2028964 0 3.208559 1
1.5320724 1 4.041846 1
1.4571594 1 4.952663 1
1.7643201 0 4.873417 1
2.0787548 1 3.705733 1
1.9274031 1 5.206371 1
1.8433588 0 5.292382 1
2.4267677 0 4.161374 1
1.1653153 1 3.748309 1
1.6580680 1 3.795499 1
2.2520952 0 5.907351 1
1.2114073 1 2.094101 1
1.1001786 1 4.487515 1
0.6709064 1 3.533209 1
1.9052090 0 3.096921 1
2.1137369 0 2.630434 1
1.0011750 0 2.520831 1
1.2564062 1 3.336696 1
2.1388844 1 4.525506 1
1.4995945 1 5.565852 1
2.0938607 1 2.838182 1
0.3142331 1 7.360565 0
1.6074203 1 11.676149 0
1.2483420 0 17.860893 0
0.5759666 0 12.963800 0
0.9307933 1 10.491123 0
1.9258344 1 8.433101 0
3.4772590 1 10.229112 0
0.2327047 0 14.152591 0
0.0388480 1 6.424880 0
0.1013466 0 6.668703 0
0.2955996 1 6.077177 0
0.7102152 1 10.387602 0
1.8412575 1 8.155387 0
0.4939213 1 7.369601 0
0.3062962 0 13.402704 0

Proposta de desenvolvimento Para o presente banco de dados, por questão de simplicidade, trabalharemos com a distribuição weibull. Para avaliarmos se a suposição é válida para o conjunto de dados, a linearização do modelo de Weibull contra o log dos tempos nos permite visualizar o ajuste, bastando para tal que os dados se apresentem na forma de uma reta de inclinação padrão.

library(survival)

ekm_data <- Surv(pratical_data2$time, pratical_data2$cens)

ekm <- survfit(ekm_data ~ 1)

ekm_surv <- ekm$surv
ekm_time <- ekm$time

plot(log(ekm_time), log(-log(ekm_surv)))

Bem como explicado anteriormente, a definição das variáveis que irão compor o modelo passam primeiramente pela construção de modelos com cada variável individualmente e sua comparação com o modelo sem covariáveis.

model_weib_X1 <- survreg(ekm_data ~ pratical_data2$metodo, dist = "weibull")

model_weib_X2 <- survreg(ekm_data ~ pratical_data2$idade, dist = "weibull")

model_weib_X1_X2 <- survreg(ekm_data ~ pratical_data2$metodo + 
    pratical_data2$idade, dist = "weibull")

loglik_X1 <- model_weib_X1$loglik
loglik_X2 <- model_weib_X2$loglik
loglik_X1_X2 <- model_weib_X1_X2$loglik

TRV_X1 = 2 * (-loglik_X1[1] + loglik_X1[2])
TRV_X2 = 2 * (-loglik_X2[1] + loglik_X2[2])

Como o teste de razão de verossimilhanças segue a distribuição qui-quadrada o cálculo do p-valor do teste pode ser realizado por meio das seguintes expressões para cada variável.

teste_x1 = 1 - pchisq(TRV_X1, 1)
teste_x2 = 1 - pchisq(TRV_X2, 1)

Encontramos, por meio do teste que ambas as variáveis são significativas tendo o método de elevação do poço um p-valor de 0.1377709 e a idade do poço um p-valor de 0.7634866. Ao observarmos o loglik de ambas as variáveis veremos que a idade (X2) o menor valor entre ambas as variáveis, se tornando então o parâmetro de comparação para o teste de razão de verossimilhança do modelo com ambas as variáveis combinadas.

TRV_X1_X2 = 2 * (-loglik_X2[2] + loglik_X1_X2[2])

teste_x1_x2 = 1 - pchisq(TRV_X1_X2, 1)

O teste da razão de verossimilhanças para o modelo com ambas as variáveis obtem então o seguinte p-valor: 0.0080007, indicando que apenas a variável X2 (idade) deve ser mantido no modelo.

Podemos compreender a cerca das variáveis através da visualização gráfica das respostas das curvas paramétrizadas. Supondo que ambas as variáveis fossem significativas, podemos propor as seguintes visualizações (1) curva de confiabilidade para um poço por método BCP com idades de 2 e 20 anos, e (2) um poço de 6 anos com método BCP e BM.

model_exp_X1_new <- exp(-((ekm_time)/(exp(model_weib_X1_X2$coefficients[1] + 
    model_weib_X1_X2$coefficients[2] * 0 + model_weib_X1_X2$coefficients[3] * 
    2)))^(1/model_weib_X1_X2$scale))

model_exp_x1_old <- exp(-((ekm_time)/(exp(model_weib_X1_X2$coefficients[1] + 
    model_weib_X1_X2$coefficients[2] * 0 + model_weib_X1_X2$coefficients[3] * 
    20)))^(1/model_weib_X1_X2$scale))

model_exp_X2_BM <- exp(-((ekm_time)/(exp(model_weib_X1_X2$coefficients[1] + 
    model_weib_X1_X2$coefficients[2] * 1 + model_weib_X1_X2$coefficients[3] * 
    6)))^(1/model_weib_X1_X2$scale))

model_exp_x2_BCM <- exp(-((ekm_time)/(exp(model_weib_X1_X2$coefficients[1] + 
    model_weib_X1_X2$coefficients[2] * 0 + model_weib_X1_X2$coefficients[3] * 
    6)))^(1/model_weib_X1_X2$scale))

par(mfrow = c(2, 1))

plot(ekm_time, model_exp_X1_new, type = "l", col = "blue", ylim = c(0, 
    1), xlab = "Tempo / anos", ylab = "R(t)")
par(new = TRUE)
lines(ekm_time, model_exp_x1_old, type = "l", col = "red", ylim = c(0, 
    1))
legend(x = 1.25, y = 0.8, col = c("blue", "red"), lty = c(1, 
    1), legend = c("Poço jovem", "Poço Antigo"))

plot(ekm_time, model_exp_X2_BM, type = "l", col = "blue", ylim = c(0, 
    1), xlab = "Tempo / anos", ylab = "R(t)")
par(new = TRUE)
lines(ekm_time, model_exp_x2_BCM, type = "l", col = "red", ylim = c(0, 
    1))
legend(x = 1.25, y = 1, col = c("blue", "red"), lty = c(1, 1), 
    legend = c("BCP", "BM"))

Prática 3 - Impacto de dados censurados

Contextualização A presença de observações de tempo-até-falha com censuras, nas quais se destacam as censuras à direita, são absorvidas pela teoria da confiabilidade (a.k.a análise de sobrevivência) através de ponderações impostas na função de verossimilhança. Conforme a Equação 9, a função de verossimilhança pode ser dividida em dois termos - relacionados (1) a contribução dos dados não censurados (função de densidade de probabilidade) e (2) a contribuição dos dados que não sofreram o evento, a falha (função de confiabilidade).

\[L(\theta) = \prod_{i = 1}^n[f(y_i|x_i)]^{(\delta_i)}[S(y_i|x_i)]^{(1-\delta_i)}\]

A variável \(y\) na função de verossimilhança corresponde a \(y = log(t)\), sendo \(y\) a variável tempo observado nos dados coletados de confiabilidade. Esta é uma transformação necessária para converter o tradicional modelo de regressão gaussiano (\(y = \beta_0 + \beta_1\cdot x + \varepsilon\)) em um modelo de natureza assimétrica, tal como são os modelos paramétricos envolvidos na análise de confiabilidade. Supondo um modelo de natureza exponencial, por exemplo, chegamos a seguinte expressão a ser maximizada pelo método da máxima verossimilhança (Equação 10).

\[logL(\theta) = \sum_{i = 1}^{n}[\delta_i(y_i-\beta_0-\beta_1\cdot x)-exp(y_i-\beta_0-\beta_1\cdot x)]\]

Pela Equação 10, podemos verificar não somente uma parcela não linear - exponencial - como também um termo linear presente apenas para observações de tempo-até-falha não censuradas.

Objetivos A presente prática busca sensibilizar os estudantes quanto ao impacto do percentual de censuras presentes nos dados de confiabilidade com foco na construção de um modelo de regressão.

Banco de dados Para possibilitarmos diferentes estudos de caso, o banco de dados possuirá natureza randômica. O mesmo corresponde ao tempo-até-falha de bombas centrifugas, adaptado de Fulano (2008), operando em diferentes alturas manometricas (15 metros e 50 metros) com diferentes simulações para o número de censuras.

Dados de tempo-até-falha de bombas centrífugas
time pres cens1 cens2 cens3
0.4251181 15 1 1 0
1.1726788 15 1 1 1
0.1086399 15 1 0 1
0.1032818 15 1 0 0
1.4067380 15 1 0 1
0.8290527 15 1 1 0
0.0970119 15 1 0 1
0.3124697 15 1 0 0
0.6223221 15 1 1 1
0.1759310 15 1 1 0
0.3245958 15 1 1 0
0.8696158 15 1 0 0
0.1165252 15 1 1 0
0.5105044 15 1 1 0
0.0255073 15 1 0 0
0.5444694 50 1 1 1
1.2598571 50 1 0 1
0.9089320 50 1 0 0
0.4502028 50 1 1 1
0.4689589 50 1 1 1
2.7432119 50 1 0 0
0.7529394 50 1 0 0
1.3531882 50 1 1 1
1.3390934 50 1 1 1
2.2957073 50 1 1 0
0.2958243 50 1 0 1
1.0258274 50 1 1 1
2.5148300 50 1 0 0
0.3671352 50 1 0 1
0.8623877 50 1 1 0
2.9491330 50 1 1 1
1.1496040 50 1 0 0
0.0823096 50 1 1 0
0.9572549 50 1 0 0
0.0589603 50 1 0 0
0.3947582 50 1 1 0
2.2795474 50 1 1 0
1.7933264 50 1 1 0
2.3512324 50 1 0 0
0.1026425 50 1 1 0
1.7369990 50 1 1 0
1.9141504 50 1 1 0
0.1297978 50 1 1 1
0.2780801 50 1 1 0
0.1564974 50 1 1 0
1.8644049 50 1 1 0
0.7512054 50 1 0 1
0.2936622 50 1 1 1
0.6679693 50 1 1 0
1.3229677 50 1 0 0

Proposta de desenvolvimento É possível observar pelos dados da Tabela 3, que os quatro casos de censura aumentam progressivamente a quantidade de observações censuradas. Por simplicidade, podemos assumir um modelo de regressão exponencial. A construção dos quatro modelos de regressão, passa, então, pela definição da curva de Kaplan-Meier para cada um dos conjuntos de dados (diferenciados pela coluna relativa às censuras nos dados) e a construção dos quatro modelos.

 

# Carregar pacote de análise de confiabilidade

library(survival)

# Construção dos 'objetos de sobrevivência'

ekm_data1 <- Surv(pratical_data1$time, pratical_data1$cens1)
ekm_data2 <- Surv(pratical_data1$time, pratical_data1$cens2)
ekm_data3 <- Surv(pratical_data1$time, pratical_data1$cens3)

# Construção das curvas de sobrevivência não paramétricas

ekm1 <- survfit(ekm_data1 ~ 1)
ekm2 <- survfit(ekm_data2 ~ 1)
ekm3 <- survfit(ekm_data3 ~ 1)

# Avaliação das curvas

par(mfrow = c(1, 3))
plot(ekm1, main = "Cens 1")
plot(ekm2, main = "Cens 2")
plot(ekm3, main = "Cens 3")

# Construção dos modelos de regressão

model_exp1 <- survreg(ekm_data1 ~ pratical_data1$pres, dist = "exponential")
model_exp2 <- survreg(ekm_data2 ~ pratical_data1$pres, dist = "exponential")
model_exp3 <- survreg(ekm_data3 ~ pratical_data1$pres, dist = "exponential")
Coeficientes para os modelos de regressão.
bias Coef..Pressão
Cens 1 -1,1084323 0,0240317
Cens 2 -0,4094080 0,0193373
Cens 3 0,0365573 0,0209399

Proposta de Discussão A estimativa dos coeficientes do modelo de regressão para confiabilidade, como pode ser observado pelos resultados da Tabela 3, estão diretamente relacionados com a quantidade de censuras e com a proporção de censuras associadas a cada nível da variável modelada. O uso de intervalos de confiança pode ser valioso para a compreensão das diferenças entre os valores dos coeficientes. Pode-se afirmar que quão maior a proporção de censuras, em termos, gerais, se obtem uma estimativa mais conservadora do impacto das variáveis no confiabilidade do objeto de estudo.