QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois...

12
ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA ACELERADO PARA DADOS DE SOBREVIVÊNCIACORRELACIONADOS: UMA APLICAÇÃO PARA POÇOS DE PETRÓLEO Patrícia Borchardt Santos Universidade Federal do Rio Grande do Norte Campus Universitário - Lagoa Nova CEP: 59072-970 - Natal, RN - Brasil [email protected] Dione Maria Valença Universidade Federal do Rio Grande do Norte Campus Universitário - Lagoa Nova CEP: 59072-970 - Natal, RN - Brasil [email protected] RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório para tratar de dados de sobrevivência correlacionados. O primeiro método, que está implementado no software SAS, utiliza a quadratura Gauss-Hermite adaptada para obter a verossimilhança marginalizada. O segundo método, implementado no software livre R, está baseado no método da verossimilhança penalizada para estimar os parâmetros do modelo. Faremos uma aplicação dos modelos usando dados reais sobre o tempo de funcionamento de poços petrolíferos da Bacia Potiguar (RN/CE). PALAVRAS CHAVE. Modelos de tempo de falha acelerado. Dados correlacionados. Petróleo. Aplicações à Indústria. ABSTRACT We presented in this work two methods of estimation for accelerated failure time models with random effects to process grouped survival data. The first method, which is implemented in software SAS uses an adapted Gauss-Hermite quadrature to determine marginalized likelihood. The second method, implemented in the free software R, is based on the method of penalized likelihood to estimate the parameters of the model. We will implement the models using actual data on the time of operation of oil wells from the Potiguar Basin (RN / CE). KEYWORDS. Accelerated failure time models. Grouped data. Oil. Applications to industry. 1. Introdução Em estudos de análise de sobrevivência, estamos interessados no tempo até que determinado evento ocorra, geralmente chamado de tempo até a falha, tempo de sobrevivência ou mesmo tempo de vida. Como exemplos, podemos ter o tempo até a cura de uma doença ou o tempo até a morte de um paciente. Para trabalharmos com dados de sobrevivência, consideramos em geral a suposição de que os tempos de vida são todos independentes. Entretanto, existem situações em que tal condição não pode ser assumida. A dependência dos tempos pode ser XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2331

Transcript of QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois...

Page 1: QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório

ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA ACELERADO PARA DADOS DE SOBREVIVÊNCIACORRELACIONADOS:

UMA APLICAÇÃO PARA POÇOS DE PETRÓLEO

Patrícia Borchardt SantosUniversidade Federal do Rio Grande do Norte

Campus Universitário - Lagoa NovaCEP: 59072-970 - Natal, RN - Brasil

[email protected]

Dione Maria ValençaUniversidade Federal do Rio Grande do Norte

Campus Universitário - Lagoa NovaCEP: 59072-970 - Natal, RN - Brasil

[email protected]

RESUMO

Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório para tratar de dados de sobrevivência correlacionados. O primeiro método, que está implementado no software SAS, utiliza a quadratura Gauss-Hermite adaptada para obter a verossimilhança marginalizada. O segundo método, implementado no software livre R, está baseado no método da verossimilhança penalizada para estimar os parâmetros do modelo. Faremos uma aplicação dos modelos usando dados reais sobre o tempo de funcionamento de poços petrolíferos da Bacia Potiguar (RN/CE).

PALAVRAS CHAVE. Modelos de tempo de falha acelerado. Dados correlacionados. Petróleo. Aplicações à Indústria.

ABSTRACT

We presented in this work two methods of estimation for accelerated failure time models with random effects to process grouped survival data. The first method, which is implemented in software SAS uses an adapted Gauss-Hermite quadrature to determine marginalized likelihood. The second method, implemented in the free software R, is based on the method of penalized likelihood to estimate the parameters of the model. We will implement the models using actual data on the time of operation of oil wells from the Potiguar Basin (RN / CE).

KEYWORDS. Accelerated failure time models. Grouped data. Oil. Applications to industry.

1. IntroduçãoEm estudos de análise de sobrevivência, estamos interessados no tempo até que

determinado evento ocorra, geralmente chamado de tempo até a falha, tempo de sobrevivência ou mesmo tempo de vida. Como exemplos, podemos ter o tempo até a cura de uma doença ou o tempo até a morte de um paciente. Para trabalharmos com dados de sobrevivência, consideramos em geral a suposição de que os tempos de vida são todos independentes. Entretanto, existem situações em que tal condição não pode ser assumida. A dependência dos tempos pode ser

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2331

Page 2: QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório

causada quando observamos tempos de sobrevivência em grupos, podendo ser, por exemplo, indivíduos de uma mesma família. Outra situação que sugere dependência entre os tempos é quando observamos mais de um tempo para um único indivíduo.

Consideramos como motivação um estudo retrospectivo sobre o tempo de funcionamento de poços de petróleo, no qual foi coletada uma amostra composta por 616 poços-coluna. No período de janeiro de 2000 à dezembro de 2006 foi observado o tempo de funcionamento de cada um dos poços-coluna dentro da sua normalidade até apresentarem falha relacionada a equipamentos de sub-superfície, que cause parada total do funcionamento do poço-coluna. Após detectada a falha, os equipamentos eram consertados e o tempo até uma nova falha era observado. Este procedimento foi feito para toda a amostra durante o período observado. Sendo assim, o evento de interesse pôde ser medido mais de uma vez em cada poço-coluna, caracterizando eventos recorrentes.

Quando as observações pertencentes a um mesmo grupo apresentam-se correlacionadas, os modelos tradicionais para análise de sobrevivência como o modelo de riscos proporcionais de Cox (Cox, 1972) ou o modelo de tempo de falha acelerado (Kalbfleisch, D. e Prentice, 2002) não podem ser diretamente aplicados. Neste enfoque, Klein (1992) propôs um modelo de riscos proporcionais de Cox com fragilidade Gama, em que as estimativas dos parâmetros foram obtidas usando um algoritmo EM. Clayton (1985) estenderam o modelo de riscos proporcionais pela adição de um termo de fragilidade multiplicativo, com distribuição Gama. Therneau (2000) fornecem uma ampla discussão sobre modelos de Cox aplicados a dados de sobrevivência correlacionados em que utilizam o método da verossimilhança penalizada para estimar os parâmetros do modelo.

No contexto paramétrico, Valença (2003) apresenta um modelo de tempo de falha acelerado (MTFA) para tratar de dados correlacionados e considerando um distribuição normal para o efeito aleatório, propõe um procedimento para obtenção de estimadores de máxima verossimilhança dos parâmetros da verossimilhança marginalizada. Considera uma aproximação da integral por uma quadratura adaptada, e a maximização da verossimilhança através de um método iterativo. O procedimento NLMIXED do SAS e o algoritmo quase-Newton foram usados para obter as estimativas. Lambert et al. (2004) em um trabalho com o mesmo enfoque, utiliza também um MTFA com um efeito aleatório. Os autores avaliam diferentes combinações da distribuição assumida para o efeito aleatório e para a função risco base, para ajustar dados de sobrevivência de pacientes com transplante renal agrupados em diferentes centros de transplante. Para o efeito aleatório os autores consideram as distribuições Gama, Normal Inversa e Log-Normal e usam também o proc NLMIXED do SAS.

Neste trabalho usamos MTFA com efeito aleatório para tratar a dependência entre tempos de sobrevivência e descrevemos os principais aspectos teóricos da abordagem usada por Valença (2003) e Lambert et al. (2004), implementado através do procedimento NLMIXED do SAS. Apresentamos brevemente a abordagem de estimação implementada no R, que se baseia na verossimilhança penalizada e realizamos um estudo de simulação para investigar a performance do método proposto.

Este trabalho está organizado da seguinte forma: Na Seção 2 introduzimos um MTFA com efeito aleatório para dados de sobrevivência correlacionados. Na Seção 3 apresentamos a abordagem para a estimação do modelo utilizando o software SAS. Uma idéia da estimação implementada no R e um estudo de simulação são dados na Seção 4. Na Seção 5 realizamos uma aplicação dos modelos usando dados reais sobre poços de petróleo. Finalizamos este trabalho com algumas conclusões na Seção 6.

2. Modelo de Tempo de Falha Acelerado com Efeito Aleatório

Seja Tij o tempo de vida do indivíduo j no grupo i, com j=1,...,ni e i=1,...,k. Seja xij um vetor de covariáveis. Denotemos os tempos de censura por Cij, e considere as respostas dadas por Yij= min(log Tij,log Cij) . O indicador de falhas ijδ é definido como )CT(Iδ ijijij ≤= , sendo I(.) a

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2332

Page 3: QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório

função indicadora. Os dados são então representados pelos pares de variáveis aleatórias ( )ijij δ,Y e as covariáveis xij. Consideramos aqui que os tempos de sobrevivência estão sujeitos a censuras à direita, assumimos que o mecanismo de censura é aleatório e que a censura é não informativa, ou seja, a sua distribuição não depende de parâmetros desconhecidos.

Podemos estender o modelo de tempo de falha acelerado pela adição de um efeito aleatório associado ao grupo. Consideramos, então, a seguinte representação para o modelo log-linear com efeito aleatório

ijijT

iij σεβUTlog ++= x (1)sendo s'ε ij erros aleatórios independentes e identicamente distribuídos, com média e variância conhecidas. O vetor ( ) T

pβ,..,β,ββ 21= e σ são parâmetros desconhecidos, e para cada grupo temos um efeito aleatório Ui, representado por variáveis aleatórias independentes e identicamente distribuídas com densidade g, médiaα e variânciaθ . Assumimos que 0=)ε,U(Cov iji , e que, condicionado ao efeito aleatório Ui, as respostas dentro do grupo i são independentes. Assumimos também que os efeitos aleatórios Ui são independentes dos tempos de censura. Este modelo se reduz ao modelo de tempo de falha acelerado usual quando 0=θ .

Assim, dado Ui, e conforme a distribuição assumida para o erro ijε , temos um particular modelo log-linear. Com uma formulação semelhante Klein, et al. (1999) estudam um modelo para o caso em que ijε e Ui possuem distribuição Normal. Outros autores, usando modelos equivalentes, discutem o ajuste do modelo de tempo de falha acelerado Weibull, em geral assumindo distribuição gama para o efeito aleatório (ver por exemplo Keiding et al., 1997 e Morris e Christiansen, 1995).

3. Abordagem implementada no software SAS

Nos modelos com efeito aleatório, os métodos de máxima verossimilhança baseiam-se, em geral, na verossimilhança marginal. Descrevemos a seguir esta verossimilhança e uma proposta para o ajuste do modelo.

3.1 Verossimilhança Marginal

Considere o modelo descrito em (1) e seja ),,,( θσβαλ T= o vetor de parâmetros desconhecidos que desejamos estimar. A verossimilhança condicional ao efeito do grupo, para o indivíduo j no grupo i é dada por

ijδijiij

ijδijiij )x,u|y(S)x,u|y(f −1

com f e S denotando, respectivamente, as funções de densidade e sobrevivência condicionais delog Tij dado o efeito do i-ésimo grupo Ui.

Denotamos o vetor de tempos observados no grupo i como Tiniiii )Y,...,Y,Y(Y 21= e

consideramos ijδ e xi definidos analogamente. Pela suposição de independência condicional ao efeito do grupo, temos que a verossimilhança condicional para os indivíduos do grupo i é da forma

,)x,u|y(S)x,u|y(f)u|σ,β(L ijδijiij

ni

j

ijδijiijii

=∏= 1

1

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2333

Page 4: QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório

para i=1,..., k. Assim, a verossimilhança relativa à distribuição marginal de ),( iiY δ , denotada por )(λiL , é representada por

.du)θ,α;u(g)u|σ,β(L)λ(L iiiii ∫=

Observe que, apesar da possível correlação existente dentro dos grupos, assumimos aqui a independência entre os vetores ),(),...,,( 11 kkYY δδ . Desta forma, o logaritmo da verossimilhança marginal que desejamos obter para toda a amostra, é dado por

.du)θ,α;u(g)u|σ,β(Llog)λ(l iiii

k

i∫∑

==

1(2)

A solução da integral em (2) e a maximização de )(λl com respeito a λ , dependem das distribuições postuladas para o efeito aleatório e para os tempos de vida.

Este problema pode ser visto como um problema geral que ocorre em modelos mistos. Em algumas situações, as distribuições são convenientemente estabelecidas de forma conjugada, como na mistura Weibull e gama, para possibilitar uma solução analítica da integral. Contudo, em situações mais gerais não é possível obter uma forma fechada para a verossimilhança. Se supomos a distribuição Normal, por exemplo, para o efeito aleatório, não obtemos uma integral analiticamente tratável.

Nossa proposta considera a distribuição Normal para o efeito aleatório e utiliza o procedimento NLMIXED para a obtenção da verossimilhança marginal, cuja abordagem básica adotada é a aproximação da integral por uma quadratura Gauss-Hermite adaptada (Liu e Pierce, 1994). A maximização de )(λl é encontrada através do algoritmo quase-Newton.

4. Abordagem implementada no software R

O segundo procedimento de estimação está implementado no R e utiliza o comando survreg com fragilidade baseado no método da verossimilhança penalizada (Therneau et. al, 2003). Como opção para o efeito aleatório temos as distribuições Gama, Gaussiana e t. Zhang (2007) comparou esse método com dois outros através de estudos de simulação, onde assumiu para a fragilidade as distribuições Gama e Log-Normal.

Embora o procedimento da verossimilhança penalizada seja bem explorado para o modelo de regressão de Cox (Therneau, 2000), não encontramos na literatura documentação abrangente para o MTFA.

Neste trabalho, como uma alternativa de compensar essa lacuna, faremos um estudo de simulação para averiguar a performance desse método.

4.1 Simulação

Apresentamos um estudo de simulação considerando um MTFA para avaliar o desempenho do método discutido acima, considerando a distribuição Normal para o efeito aleatório, e comparamos com um método que supõe que os tempos de vida são todos independentes, ambos implementados no software R (versão 2.7.1).

Para este estudo assumimos grupos com três tamanhos distintos: 10, 100 e 500, com 5 indivíduos em cada grupo. Consideramos que temos uma única covariável xij associada a cada indivíduo j no grupo i, com distribuição Normal padrão. Os tempos de censura foram gerados a partir de uma distribuição uniforme (0, q), onde q é escolhido de forma a fornecer 0% e 30% de

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2334

Page 5: QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório

censura na amostra. Os erros aleatórios ijε foram gerados a partir de variáveis aleatórias independentes e identicamente distribuídas com distribuição valor extremo padrão.

Realizamos 1000 réplicas para cada amostra e os verdadeiros valores para os coeficientes do modelo foram 5=α , 5,0=β e 1=σ . Além disso, consideramos quatro valores para a variância do efeito aleatório: θ =0 (sem dependência entre os tempos), θ =0,25, θ=0,5 e θ =1.

As estimativas da média, erro padrão (EP) e raiz do erro quadrático médico (REQM) para os dados simulados estão nas Tabelas 1 e 2.

Tabela 1: Estimativas da média, erro padrão e raiz do erro quadrático médio dos parâmetros com ,5=α ,5,0=β ,1=σ 1000 réplicas, amostras divididas em 10, 100 e 500 grupos, sem censura e

variando o valor de θ .

θ Par Método k =10 k =100 k =500Média EP REQM Média EP REQM Média EP REQM

0

θ Ind - - - - - - - - -Pen 0,0640 0,0941 0,1138 0,0189 0,0254 0,0317 0,0079 0,0099 0,0126

α Ind 4,9897 0,1543 0,1547 4,9992 0,0496 0,0496 4,9997 0,0213 0,0213

Pen 4,9717 0,1560 0,1586 4,9887 0,0489 0,0502 4,9952 0,0209 0,0214

β Ind 0,5056 0,1527 0,1528 0,4991 0,0446 0,0446 0,5008 0,0198 0,0199

Pen 0,4960 0,1493 0,1493 0,4986 0,0446 0,0446 0,4996 0,0200 0,0200

σ Ind 0,9728 0,1084 0,1118 0,9955 0,0339 0,0342 0,9999 0,0161 0,0161

Pen 0,9251 0,1233 0,1443 0,9779 0,0422 0,0476 0,9904 0,0186 0,0210

0,25

θ Ind - - - - - - - - -

Pen 0,2944 0,2204 0,2248 0,2956 0,0747 0,0875 0,2933 0,0327 0,0542

α Ind 5,0335 0,2107 0,2134 5,0562 0,0695 0,0894 5,0556 0,0322 0,0643

Pen 4,9579 0,2166 0,2206 4,9610 0,0725 0,0823 4,9611 0,0304 0,0493

β Ind 0,5081 0,1809 0,1810 0,5043 0,0538 0,0540 0,4995 0,0242 0,0242

Pen 0,5039 0,1724 0,1724 0,5027 0,0507 0,0508 0,4991 0,0223 0,0223

σ Ind 1,0875 0,1317 0,1581 1,1350 0,0436 0,1418 1,1379 0,0198 0,1393

Pen 0,8954 0,1303 0,1671 0,9081 0,0415 0,1008 0,9102 0,0182 0,0917

0,5

θ Ind - - - - - - - - -

Pen 0,5560 0,3716 0,3757 0,5496 0,1013 0,1125 0,5546 0,0477 0,0724

α Ind 5,1017 0,2801 0,2980 5,1050 0,0873 0,1366 5,1077 0,0384 0,1143

Pen 4,9392 0,2770 0,2836 4,9519 0,0816 0,0946 4,9569 0,0388 0,0580

β Ind 0,4910 0,1963 0,1966 0,4998 0,0621 0,0621 0,5002 0,0284 0,0284

Pen 0,5065 0,1771 0,1772 0,5036 0,0499 0,0499 0,4999 0,0232 0,0232

σ Ind 1,1848 0,1597 0,2441 1,2473 0,0594 0,2543 1,2565 0,0251 0,2578

Pen 0,8877 0,1243 0,1675 0,8883 0,0419 0,1192 0,8902 0,0175 0,1112

1

θ Ind - - - - - - - - -

Pen 1,1361 0,6010 0,6107 1,0688 0,1727 0,1859 1,0653 0,0787 0,1022

α Ind 5,1511 0,3612 0,3915 5,1924 0,1184 0,2259 5,2004 0,0515 0,2069

Pen 4,8595 0,2826 0,3133 4,9538 0,1077 0,1172 4,9524 0,0511 0,0698

β Ind 0,5064 0,2329 0,2330 0,4979 0,0727 0,0727 0,4994 0,0337 0,0337

Pen 0,5052 0,1760 0,1744 0,4946 0,0499 0,0502 0,4992 0,0238 0,0238

σ Ind 1,3450 0,2036 0,4006 1,4434 0,0750 0,4497 1,4571 0,0355 0,4585

Pen 0,8756 0,1221 0,1735 0,8766 0,0372 0,1289 0,8783 0,0169 0,1229

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2335

Page 6: QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório

Tabela 2: Estimativas da média, erro padrão e raiz do erro quadrático médio dos parâmetros com ,5=α ,5,0=β ,1=σ 1000 réplicas, amostras divididas em 10, 100 e 500 grupos, com 30% de

censura e variando o valor de θ .

θ Par Método k =10 k =100 k =500Média EP REQM Média EP REQM Média EP REQM

0

θ Ind - - - - - - - - -

Pen 0,0878 0,1310 0,1577 0,0240 0,0355 0,0429 0,0109 0,0146 0,0183

α Ind 5,0029 0,1728 0,1729 4,9989 0,0533 0,0533 4,9996 0,0252 0,0252

Pen 4,9572 0,1783 0,1833 4,9865 0,0557 0,0573 4,9938 0,0254 0,0261

β Ind 0,5246 0,1879 0,1895 0,5041 0,0567 0,0569 0,5012 0,0252 0,0252

Pen 0,5164 0,1913 0,1920 0,5007 0,0581 0,0581 0,5000 0,0260 0,0260

σ Ind 0,9745 0,1360 0,1384 0,9980 0,0432 0,0432 0,9988 0,0185 0,0185

Pen 0,9222 0,1544 0,1729 0,9759 0,0505 0,0559 0,9877 0,0240 0,0267

0,25

θ Ind - - - - - - - - -

Pen 0,3004 0,2699 0,2745 0,2711 0,0832 0,0858 0,2715 0,0362 0,0421

α Ind 5,0282 0,2429 0,2446 5,0283 0,0787 0,0836 5,0298 0,0342 0,0453

Pen 4,9266 0,2334 0,2447 4,9266 0,0731 0,1036 4,9228 0,0334 0,0841

β Ind 0,5070 0,2096 0,2097 0,4886 0,0620 0,0630 0,4878 0,0283 0,0308

Pen 0,4967 0,2019 0,2019 0,4826 0,0607 0,0631 0,4861 0,0263 0,0298

σ Ind 1,0617 0,1436 0,1563 1,0823 0,0486 0,0956 1,0865 0,0209 0,0890

Pen 0,8863 0,1501 0,1833 0,8994 0,0535 0,1139 0,8986 0,0233 0,1041

0,5

θ Ind - - - - - - - - -

Pen 0,5451 0,3867 0,3893 0,5024 0,1095 0,1096 0,4977 0,0476 0,0477

α Ind 5,0482 0,2927 0,2967 5,0539 0,0903 0,1051 5,0561 0,0400 0,0689

Pen 4,8932 0,2824 0,3019 4,8996 0,0860 0,1322 4,9016 0,0388 0,1058

β Ind 0,5020 0,2320 0,2320 0,4816 0,0696 0,0720 0,4781 0,0313 0,0381

Pen 0,5014 0,2031 0,2031 0,4843 0,0600 0,0620 0,4811 0,0266 0,0326

σ Ind 1,1348 0,1616 0,2105 1,1603 0,0516 0,1684 1,1603 0,0230 0,1620

Pen 0,8601 0,1490 0,2044 0,8732 0,0456 0,1347 0,8722 0,0201 0,1294

1

θ Ind - - - - - - - - -

Pen 0,9729 0,5759 0,5765 0,9346 0,1701 0,1823 0,9280 0,0779 0,1061

α Ind 5,0945 0,3623 0,3744 5,1126 0,1148 0,1608 5,1118 0,0488 0,1220

Pen 4,8879 0,3476 0,3652 4,8720 0,1111 0,1695 4,8735 0,0485 0,1355

β Ind 0,4939 0,2601 0,2602 0,4697 0,0773 0,0830 0,4690 0,0338 0,0459

Pen 0,4837 0,2174 0,2181 0,4763 0,0612 0,0656 0,4809 0,0278 0,0338

σ Ind 1,2505 0,1975 0,3192 1,2966 0,0597 0,3025 1,3006 0,0270 0,3018

Pen 0,8398 0,1350 0,2095 0,8515 0,0429 0,1546 0,8503 0,0194 0,1509

4.1.2 Resultados

As Tabelas 1 e 2 mostram que, em geral, as estimativas fornecidas pelo método que supõe independência são maiores que as dadas pelo método que trata a dependência entre os tempos. A maior diferença entre os dois métodos foi observada quanto às estimativas de σ . Com o aumento da variância do efeito aleatório, as estimativas obtidas pelo modelo sob independência tendem a crescer, tornando-se mais distantes do verdadeiro valor. O contrário acontece com o método que utiliza a penalização da verossimilhança, cujas estimativas de σ diminuem com o aumento de θ , mas também afastam-se do valor desejado. Essa regularidade foi notada para ambos os percentuais de censura.

Quanto ao crescimento dos grupos e, consequentemete, com o aumento da amostra, evidenciamos que os valores de σ crescem e se afastam do valor almejado se ignorarmos a

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2336

Page 7: QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório

dependência entre os tempos (exceto quando temos θ =0) . As estimativas fornecidas supondo dependência também crescem com o aumento dos grupos, mas, ao contrário das obtidas sob independência, tornam-se cada vez mais próximas de 1. Essa tendência é melhor evidenciada na ausência de censura.

Aumentando a proporção de censura, as estimativas dadas pelos dois métodos tendem a diminuir. Para o modelo com efeito aleatório elas tornam-se mais distantes do verdadeiro valor. O contrário ocorre para o método usual, uma vez que esse método superestima as estimativas de σ . De forma que podemos acreditar que, com exceção do caso em que θ =0, o aumento da censura aperfeiçoa a estimava do σ para tal método. É importante ressaltar, porém, que embora as estimativas sob este método estejam se aproximando dos valores desejados, o método que trata a dependência ainda fornece estimativas menos viesadas para σ .

Nas Figuras 1, 2 e 3 apresentamos alguns gráficos para ilustrar as interpretações obtidas.

0.0 0.2 0.4 0.6 0.8 1.0

0.8

0.9

1.0

1.1

1.2

1.3

1.4

k=10, sem censuraθ

Est

imat

iva

IndependênciaPenalizada

0.0 0.2 0.4 0.6 0.8 1.0

0.8

0.9

1.0

1.1

1.2

1.3

1.4

k=10, 30% de censuraθ

Est

imat

iva

IndependênciaPenalizada

Figura 1: Estimativas de σ para o método que assume independência e para o método baseado na verossimilhança penalizada para 0% e 30% de censura, com 10 grupos.

0.0 0.2 0.4 0.6 0.8 1.0

0.8

0.9

1.0

1.1

1.2

1.3

1.4

k=100, sem censuraθ

Est

imat

iva

IndependênciaPenalizada

0.0 0.2 0.4 0.6 0.8 1.0

0.8

0.9

1.0

1.1

1.2

1.3

1.4

k=100, 30% de censuraθ

Est

imat

iva

IndependênciaPenalizada

Figura 2: Estimativas de σ para o método que assume independência e para o método baseado na verossimilhança penalizada para 0% e 30% de censura, com 100 grupos.

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2337

Page 8: QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório

0.0 0.2 0.4 0.6 0.8 1.0

0.8

0.9

1.0

1.1

1.2

1.3

1.4

k=500, sem censuraθ

Est

imat

iva

IndependênciaPenalizada

0.0 0.2 0.4 0.6 0.8 1.0

0.8

0.9

1.0

1.1

1.2

1.3

1.4

k=500, 30% de censuraθ

Est

imat

iva

IndependênciaPenalizada

Figura 3: Estimativas de σ para o método que assume independência e para o método baseado na verossimilhança penalizada para 0% e 30% de censura, com 500 grupos.

5. Aplicação

Para ilustrar o método, faremos uma aplicação usando o conjunto de dados que mencionamos na introdução sobre o tempo de funcionamento de poços de petróleo. Como vimos, após detectada a falha os equipamentos eram consertados e o tempo até uma nova falha era observado, de modo que dos 616 poços-coluna obtivemos 2374 observações, sendo 563 censuradas (23,7%). O objetivo do estudo era verificar quais covariáveis influenciam no tempo de funcionamento dos poços de petróleo. As covariáveis selecionadas foram:

• Produção base do poço, medida em 3m /dia;• Método de elevação: método artificial utilizado para elevar o fluido, contido no fundo do

poço, até a superfície. Neste estudo, os dois métodos considerados foram BM - bombeio mecânico e BCP - bombeio por cavidade progressiva, (com variável indicadora BM assumindo o valor 1 se o método é Bombeio Mecânico e 0 caso contrário);

• Idade do poço medida no momento da falha, em anos;• Unidade administrativa: unidade que administra os poços de acordo com a sua

localização geográfica. São quatro unidades: OP-ARG (Unidade Operacional Alto do Rodrigues) OP-CAM (Unidade Operacional Canto do Amaro), OP-ET (Unidade Operacional Campo de Estreito) e OP-RFQ (Unidade Operacional Fazenda Riacho da Forquilha).

• Profundidade da bomba: profundidade onde se encontra instalada a bomba de produção do poço, em metros.

Para explorar um pouco mais os dados, estimamos a função de sobrevivência usando o estimador não paramétrico de Kaplan-Meier para os tempos de falha dos poços por unidade administrativa e por método de elevação (Figura 4), observamos que o Bombeio Mecânico (BM) como método de elevação parece proporcionar maior tempo de funcionamento dos poços do que o Bombeio por Cavidade Progressiva (BCP). Quanto às unidades administrativas, a região do Canto do Amaro (CAM) apresenta maior sobrevivência do que as regiões Alto do Rodrigues (ARG), Estreito (ET) e Fazenda Riacho da Forquilha (RFQ).

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2338

Page 9: QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório

0 10000 30000 50000

0.0

0.2

0.4

0.6

0.8

1.0

Curva Estimada KM

Tempo Entre Falhas

Sob

revi

vênc

ia

BCPBM

0 10000 30000 50000

0.0

0.2

0.4

0.6

0.8

1.0

Curva Estimada KM

Tempo Entre falhas

Sob

revi

vênc

ia

OP-ARGOP-CAMOP-ETOP-RFQ

Figura 4: Estimador de Kaplan-Meier para os tempos de falha dos poços por Unidade Administrativa e por Método de Elevação

Seja ijy o logaritmo do tempo decorrido até o poço i parar de funcionar ou ser censurado, na j-ésima recorrência.

Considerando que observações de um mesmo poço-coluna podem estar relacionadas, assumimos um modelo Weibull com efeito aleatório para os dados. Após uma análise prévia optamos pelo seguinte modelo:

+++++++= irfqieticamijidiBMijprodiij RFQβETβCAMβIDADEβBMβPRODβUTlog++++ iijrfq*prodiijet*prodiijcam*prodiprofb RFQ*PRODβET*PRODβCAM*PRODβPROFBβ

ijiiprofb*rfqiiprofb*etiiprofb*cam εσPROFB*RFQβPROFB*ETβPROFB*CAMβ +++

com )θ,α(N~U i , representando o efeito aleatório do poço i e com ijε representando o erro aleatório do modelo, com distribuição valor extremo padrão.

Os resultados para o ajuste sob homogeneidade e para os dois ajustes do modelo com efeito aleatório estão na Tabela 3.

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2339

Page 10: QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório

Tabela 3: Estimação de máxima verossimilhança dos parâmetros para dados sobre poços de petróleo utilizando a distribuição Weibull.

Sob independência Com efeito aleatório (SAS) Com efeito aleatório (R)Parâmetros Estim. EP P-valor Estim. EP P-valor Estim. EP P-valor

α 7,4688 0,2012 <0,0001 7,4689 0,2588 <0,0001 7,2224 0,2465 <0,0001

prodβ -0,0502 0,0084 <0,0001 -0,0500 0,0091 <0,0001 -0,0457 0,0084 <0,0001

BMβ 0,5197 0,1002 <0,0001 0,5291 0,1291 <0,0001 0,6710 0,1218 <0,0001

idβ 0,0587 0,0059 <0,0001 0,0793 0,0078 <0,0001 0,0852 0,0069 <0,0001

camβ 1,3421 0,2472 <0,0001 1,3430 0,3090 <0,0001 1,4302 0,2954 <0,0001

etβ 0,9612 0,1685 <0,0001 0,9575 0,2275 <0,0001 0,8075 0,2189 <0,0001

rfqβ 1,8245 0,2572 <0,0001 1,8259 0,3304 <0,0001 1,7816 0,3172 <0,0001

profbβ 0,0022 0,0003 <0,0001 0,0021 0,0004 <0,0001 0,0021 0,0004 <0,0001

camprod*β 0,0320 0,0143 0,0252 0,0392 0,0157 0,0128 0,0347 0,0145 0,0162

etprod*β 0,0186 0,0119 0,1190 0,0176 0,0138 0,2035 0,0158 0,0128 0,2160

rfqprod*β -0,0777 0,0160 <0,0001 -0,0458 0,0187 0,0144 -0,0380 0,0171 0,0263

profbcam*β -0,0021 0,0004 <0,0001 -0,0020 0,0005 0,0002 -0,0021 0,0005 0,0001

profbet*β -0,0027 0,0004 <0,0001 -0,0028 0,0006 <0,0001 -0,0026 0,0005 <0,0001

profbrfq*β -0,0027 0,0005 <0,0001 -0,0027 0,0006 <0,0001 -0,0028 0,0006 <0,0001σ 1,34 <0,0001 1,1951 0,0243 <0,0001 1,0800θ - - - 0,4529 0,0597 <0,0001 0,4740 ? ?

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2340

Page 11: QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório

A estimativa da variância do efeito aleatório é significativa ( θ̂ = 0,4529 para a proposta implementada no SAS e θ̂ = 0,4740 para a proposta via R).

Com respeito aos objetivos do estudo, notamos que os poços-coluna nos quais o método de elevação utilizado foi o Bombeio Mecânico apresentaram maior tempo de funcionamento do que aqueles cujo método de elevação utilizado foi o Bombeio por Cavidade Progressiva BMβ̂ = 0,5291 e BMβ̂ = 0,6710, proposta SAS e R, respectivamente, ambas com p-valor menor que 0,0001. Também observamos que o tempo de funcionamento dos poços-coluna localizados nas unidades operativas ET, CAM e RFQ apresentaram maior tempo de funcionamento do que na unidade ARG.

Verificamos também que à medida que a produção aumenta, o funcionamento do poço-coluna diminui, pois a estimativa da produção foi negativa e significativa. De acordo com o aumento da idade o poço-coluna tende a funcionar por mais tempo. O mesmo acontece com a profundidade da bomba, ou seja, quanto mais profunda estiver instalada a bomba, maior é o tempo livre de falha dos poços-coluna.

Além disso, consideramos relevante a interação entre a produção do poço e a unidade que o administra. É possível concluir que os poços com produção elevada têm menor tempo de funcionamento se estiverem localizados nas unidades CAM, RFQ e ARG. No entanto, não é possível perceber que a produção afeta significativamente o tempo de funcionamento dos poços localizados na unidade ET.

Por fim, observamos a relação entre a profundidade da bomba e a unidade administrativa. Constatamos que poços mais profundos tendem a ter maior tempo de sobrevivência se estiverem nas unidades CAM e ARG. Nas unidades ET e RFQ o contrário acontece, de modo que o aumento da profundidade da bomba diminui o tempo de vida dos poços.

6. Conclusões Finais

Neste trabalho estudamos dois procedimentos para ajustar um MTFA com efeito aleatório a dados de sobrevivência correlacionados, um implementado através do procedimento NLMIXED do SAS, e o outro implementada no R. No primeiro caso descrevemos os aspectos teóricos do procedimento, e no segundo caso R realizamos um estudo de simulação para investigar a performance do método. Aplicamos ambos os métodos a um conjunto de dados relativos a tempos entre falhas de equipamento de sub-superficie de poços de petróleo e observamos que apesar de não haver descrição detalhada na literatura a implementação usada no software livre R apresenta resultados próximos daqueles obtidos pelo SAS. Notamos também que embora nessa aplicação as estimativas dos parâmetros sejam próximas tanto para o modelo que assume independência quanto para os modelos com efeito aleatório, lembramos que através da simulação foi observado que o modelo usual pode trazer grandes distorções principalmente na estimativa de σ e, consequentemente, prejuízo na interpretação.

AgradecimentosAs autoras agradecem à PETROBRAS por ter financiado parcialmente este trabalho através do projeto SAFES – Sistema de Análises de Falha de Equipamentos de sub-superfícies.

Referências

Clayton, D. e J. Cuzick, (1985), Multivariate generalizations of the proportional hazards model. Journal of the Royal Statistical Society, Series A 148, 82-117.Cox, D. R. (1972). Regression models and life tables (with discussion). Journal of the Royal Sta-tistical Society, Series B, 34, 187-220.

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2341

Page 12: QUALI:ESTIMAÇÃO EM MODELOS DE TEMPO DE FALHA … · RESUMO Apresentamos neste trabalho dois métodos de estimação para modelos de tempo de falha acelerado com efeito aleatório

Dantas, M. A. (2008), Modelos de Dados de Falhas de Equipamentos de Sub-superfície em Poços de Petróleo da Bacia Potiguar, Dissertação de mestrado, PEP UFRN.Keiding, N., Andersen, P. K. e Klein, J. P. (1997). The role of frailty models and accelerated failure time models in describing heterogeneity due to omitted covariates. Statistics in Medicine, 16, 215-224.Klein, J. (1992), Semiparametric estimation of random effects using the Cox model based on the EM algorithm. Biometrics, 795-806.Klein, J. P., Pelz, C. e Zhang, M. (1999). Modelling random effects for censored data by a mul-tivariate normal regression model. Biometrics, 55, 497-506.Lambert, P., D. Collett, A. Kimber, e R. Johnson (2004), Parametric accelerated failure time models with random effects and an application to kidney transplant survival. Statistics in Medicine. 23 (20).Liu, Q. e Pierce, D.A. (1994), A note on Gauss-Hermite quadrature. Biometrika, volume 81, número 3, 624-629.Morris, C. e Christiansen, C. (1995). FittingWeibull duration models with random effects. Life-time Data Analysis, 1, 347-359.Therneau, T., P. Grambsch, Modeling survival data: extending the Cox model. Springer, New York, 2000.Therneau, T. M. e Grambsch, P. M., Pankratz, V. S., (2003), Penalized survival models and frailty. J. Comput. Graph. Statist. 12 (1), 156 – 175.Valença, D. M. (2003), Teste de Homogeneidade e Estimação para Dados de Sobrevivência Agrupados e com Erros de Medida, Tese de doutorado, IME USP.Zhang, J. e Peng, Y. (2007), An alternative estimation method for the accelerated failure time frailty model. Computational Statistics and Data Analysis, Elsevier, vol. 51, nº 9, 4413-4423.

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 2342