ANÁLISE DE SOBREVIVÊNCIA E CONFIABILIDADE

AULA 3: MODELOS PROBABILÍSTICOS: ESTIMAÇÃO

\[\\[0.05in]\]

1 PACOTES NECESSÁRIOS

Para esta terceira aula prática de Análise de Sobrevivência e Confiabilidade, precisaremos dos seguintes pacotes:

require(flexsurv)
require(reslife)

• Pacote flexsurv: Responsável pelo ajuste de modelos paramétricos flexíveis aos dados de tempo até a ocorrência de um evento de interesse.

• Pacote reslife: Responsável por fornecer uma estimativa para o tempo médio residual de ocorrência de um evento de interesse.

2 BASES DE DADOS

Esta seção é dedicada à descrição das bases de dados utilizadas nesta aula.

2.1 CÂNCER DE BEXIGA

Essa base de dados é proveniente de assessorias estatísticas realizadas no Departamento de Estatística da UFPR. Os dados contém os tempos de reincidência, em meses, de um grupo de 20 pacientes com câncer de bexiga que foram submetidos a um procedimento cirúrgico realizado por laser. Os tempos obtidos foram:

3, 5, 6, 7, 8, 9, 10, 10+, 12, 15, 15+, 18, 19, 20, 22, 25, 28, 30, 40, 45+,

em que o simbolo + indica censura.

2.2 TRATAMENTO QUIMIOTERÁPICO

No estudo analisado neste exemplo, são apresentados os tempos, em dias até a ocorrência dos primeiros sinais de alterações indesejadas no estado geral de saúde de 45 pacientes de ambos os sexos que receberam tratamento quimioterápico após terem sido submetidos à cirurgia de intestino. Houve um acompanhamento total de 250 dias desde a entrada do primeiro paciente até o término do estudo.

7, 8, 10, 12, 13, 14+, 19, 23, 25+, 26, 27, 31, 31+, 49, 59+, 64+, 87, 89, 107, 117, 119, 230+, 233+, 130, 148, 153, 156, 159, 191, 222, 200+, 203+, 210+, 220+, 220+, 228+, 235+, 240+ , 240+, 240+, 241+, 245+, 247+, 248+, 250+,

em que o símbolo + indica a censura.

2.3 ISOLANTE ELÉTRICO

0 fabricante de um tipo de isolante elétrico quer conhecer o comportamento de seu produto funcionando a uma temperatura de 200Cº. Um teste de vida foi realizado nestas condições uasando-se 60 isoladores elétricos. O teste terminou quando 45 deles havia falhado (censura do tipo II). As 15 unidades que não haviam falhado ao final do teste foram, desta forma, censuradas no tempo t= 2729 horas:

151, 164, 336, 365, 403, 454, 455, 473, 538, 577, 592, 628, 632, 647, 675, 727, 785, 801, 811, 816, 867, 893, 930, 937, 976, 1008, 1040, 1051, 1060, 1183, 1329, 1334, 1379, 1380, 1633, 1769, 1827, 1831, 1849, 2016, 2282, 2415, 2430, 2686, 2729, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+,

em que o símbolo + indica a censura.

3 AJUSTE DOS MODELOS PARAMÉTRICOS

Nesta seção, aprenderemos como ajustar todos os modelos paramétricos vistos em sala de aula. Mas antes de fazer esses ajustes, vamos carregar a base de dados dos pacientes com câncer de bexiga:

 tempos<-c(3,5,6,7,8,9,10,10,12,15,15,18,19,20,22,25,28,30,40,45)
 cens<-c(1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0)
 dados1<-data.frame(tempos,cens)

O ajuste de todos os modelos paramétricos será feito a partir da função flexsurvreg. Essa função possui, basicamente, os seguintes argumentos:

• formula: a fórmula do modelo. No caso de modelos sem covariáveis, a sua estrutura é, por exemplo, Surv(tempos, cens) ~ 1

• dist: a distribuição utilizada. Pode ser uma das seguintes: gengamma.orig, weibull, gamma, exp, llogis e lnorm.

• data: a base de dados.

Em todos os modelos ajustados, a função flexsurvreg estima os parâmetros do modelo pelo método da máxima verossimilhança, utilizando métodos numéricos. Além disso, a função reporta um intervalo de confiança para os parâmetros e também uma estimativa para os respectivos erros padrões. Na saída, são também reportados os seguintes:

• N: tamanho da amostra.
• Events: número de ocorrências do evento.
• Censored: número de observações censuradas.
• Total time at risk: soma de todos os tempos observados (censurados ou não).
• Log-likelihood: valor da função de log-verossimilhança.
• AIC: Critério de Informação de Akaike

Para mais informações e detalhes sobre essa função, utilize o comando ?flexsurvreg()

Já as medidas de interesse com respeito à análise de sobrevivência e confiabilidade, como os tempos médio e mediano de ocorrência do evento de interesse e percentis, são obtidos a partir da função summary. Essa função recebe como argumento o modelo criado e os outros argumentos dessa função são:

• type = "mean": fornece o tempo médio de ocorrência do evento de interesse.
• type = "median": fornece o tempo mediano de ocorrência do evento de interesse.
• type = "quantile": fornece um percentil do tempo de ocorrência do evento de interesse. Para escolher o percentil, acrescente o argumento quantiles = 0,8, por exemplo.
• cl=0.95 = nível de confiança usado nos intervalos.
• se=T = reporta o erro padrão atrelado à estimativa.

Para mais informações e detalhes sobre essa função, utilize o comando ?summary.flexsurvreg()

O tempo médio residual de ocorrência do evento de interesse no tempo t é o tempo médio para a ocorrência do evento dado que o elemento não apresentou o evento até o instante t. É como se fosse a vida média residual. No R, essa quantidade é calculada com a função reslifefsr do pacote reslife. A sintaxe é simples, basta digitar reslifefsr(ajuste1,20), que retorna o tempo médio residual de ocorrência do evento de interesse atrelado ao modelo ajuste1 no tempo t=20.

Para mais informações e detalhes sobre essa função, utilize o comando ?reslifefsr()

3.1 DISTRIBUIÇÃO EXPONENCIAL

Ao assumir que o tempo \(T\) possui uma distribuição exponencial (\(T\sim Exponencial(\alpha))\), a sua função densidade de probabilidade é dada por:} \[f(t) = \alpha exp\{-\alpha t\},\] em que \(\alpha > 0\) é o parâmetro de taxa e \(t\geq0\).

Utilizando as relações entre as funções básicas estudadas anteriormente, as funções de sobrevivência (confiabilidade), de risco e de risco acumulado são

\[S(t)=exp\{-\alpha t\}\]

\[\lambda(t)=\alpha\]

\[\Lambda(t)=\alpha t\]

Para o ajuste do modelo exponencial, vamos utilizar a função flexsurvreg com dist="exp":

ajuste_exp<-flexsurvreg(Surv(tempos, cens) ~ 1, dist = "exp", data=dados1)
ajuste_exp

## Call:
## flexsurvreg(formula = Surv(tempos, cens) ~ 1, data = dados1, 
##     dist = "exp")
## 
## Estimates: 
##       est     L95%    U95%    se    
## rate  0.0490  0.0305  0.0788  0.0119
## 
## N = 20,  Events: 17,  Censored: 3
## Total time at risk: 347
## Log-likelihood = -68.27389, df = 1
## AIC = 138.5478

Assim, \(\hat{\alpha}=0,0490\) é o estimador de máxima verossimilhança para \(\alpha\).

Os comandos a seguir fazem os gráficos das funções de sobrevivência, risco e risco acumulado para este modelo, comparando com a estimação não-paramétrica:

plot(ajuste_exp,xlab="Tempos (meses)",ylab="Função de Sobrevivênca S(t)")

plot(ajuste_exp,xlab="Tempos (meses)",ylab=expression("Função de Risco " ~ lambda(t) ),type="hazard")

plot(ajuste_exp,xlab="Tempos (meses)",ylab=expression("Função de Risco Acumulado " ~ Lambda(t) ),type="cumhaz")

Os comandos a seguir calculam, na sequência, o tempo médio de ocorrência do evento de interesse, o tempo mediano de ocorrência do evento de interesse e o percentil 80 do tempo de ocorrência do evento de interesse

Mean_exp<-summary(ajuste_exp,type="mean",cl=0.95,se=T)
Mean_exp

##  
##        est     lcl      ucl       se
## 1 20.41176 12.7985 32.74889 5.132306

Median_exp<-summary(ajuste_exp,type="median",cl=0.95,se=T)
Median_exp

##  
##        est      lcl      ucl       se
## 1 14.14836 9.178828 22.41646 3.488359

quantile80_exp<-summary(ajuste_exp,type="quantile",
                        quantiles=0.8,cl=0.95,se=T)
quantile80_exp

##  
##   quantile      est      lcl      ucl       se
## 1      0.8 32.85147 20.99045 53.03236 8.166447

Já o comando a seguir fornece o tempo médio residual de ocorrência do evento de interesse no instante t=20 meses.

Resid20_exp<-reslifefsr(ajuste_exp,20)
Resid20_exp

## [1] 20.41176

3.2 DISTRIBUIÇÃO WEIBULL

Ao assumir que o tempo \(T\) possui uma distribuição Weibull (\(T\sim Weibull(\alpha,\gamma)\)), a sua função densidade de probabilidade é dada por:

\[f(t) = \frac{\gamma t^{\gamma-1}}{\alpha^\gamma} exp\left\{-\left(\frac{t}{\alpha }\right)^\gamma\right\},\]

em que \(\alpha > 0\) é o parâmetro de escala, \(\gamma > 0\) é o parâmetro de forma e \(t\geq0\).

Utilizando as relações entre as funções básicas estudadas anteriormente, as funções de sobrevivência (confiabilidade), de risco e de risco acumulado são

\[S(t)=exp\left\{-\left(\frac{t}{\alpha }\right)^\gamma\right\}\]

\[\lambda(t)=\frac{\gamma t^{\gamma-1}}{\alpha^\gamma}\]

\[\Lambda(t) = \left(\frac{t}{\alpha}\right)^{\gamma-1}\]

Para o ajuste do modelo Weibull, vamos utilizar a função flexsurvreg com dist="weibull":

ajuste_wei<-flexsurvreg(Surv(tempos, cens) ~ 1, dist = "weibull", data=dados1)
ajuste_wei

## Call:
## flexsurvreg(formula = Surv(tempos, cens) ~ 1, data = dados1, 
##     dist = "weibull")
## 
## Estimates: 
##        est     L95%    U95%    se    
## shape   1.543   1.066   2.235   0.291
## scale  21.339  15.591  29.206   3.417
## 
## N = 20,  Events: 17,  Censored: 3
## Total time at risk: 347
## Log-likelihood = -66.13336, df = 2
## AIC = 136.2667

Assim, \(\hat{\gamma}=1,543\) é o estimador de máxima verossimilhança para \(\gamma\) e \(\hat{\alpha}=1,543\) é o estimador de máxima verossimilhança para \(\alpha\).

Os comandos a seguir fazem os gráficos das funções de sobrevivência, risco e risco acumulado para este modelo, comparando com a estimação não-paramétrica:

plot(ajuste_wei,xlab="Tempos (meses)",ylab="Função de Sobrevivênca S(t)")

plot(ajuste_wei,xlab="Tempos (meses)",ylab=expression("Função de Risco " ~ lambda(t) ),type="hazard")

plot(ajuste_wei,xlab="Tempos (meses)",ylab=expression("Função de Risco Acumulado " ~ Lambda(t) ),type="cumhaz")

Os comandos a seguir calculam, na sequência, o tempo médio de ocorrência do evento de interesse, o tempo mediano de ocorrência do evento de interesse e o percentil 80 do tempo de ocorrência do evento de interesse:

Mean_wei<-summary(ajuste_wei,type="mean",cl=0.95,se=T)
Mean_wei

##  
##        est      lcl      ucl       se
## 1 19.20078 13.98625 26.24761 3.102184

Median_wei<-summary(ajuste_wei,type="median",cl=0.95,se=T)
Median_wei

##  
##        est      lcl      ucl       se
## 1 16.82823 11.97493 23.21932 2.937372

quantile80_wei<-summary(ajuste_wei,type="quantile",
                        quantiles=0.8,cl=0.95,se=T)
quantile80_wei

##  
##   quantile      est      lcl      ucl      se
## 1      0.8 29.04556 21.71682 38.51505 4.50519

Já o comando a seguir fornece o tempo médio residual de ocorrência do evento de interesse no instante t=20 meses.

Resid20_wei<-reslifefsr(ajuste_wei,20)
Resid20_wei

## [1] 11.57553

3.3 DISTRIBUIÇÃO LOGNORMAL

Ao assumir que o tempo \(T\) possui uma distribuição lognormal (\(T\sim lognormal(\mu,\sigma)\)), a sua função densidade de probabilidade é dada por:

\[f(t)=\frac{1}{\sqrt{2\pi} t \sigma}exp\left\{-\frac{1}{2}\left(\frac{ln(t)-\mu}{\sigma}\right)^2\right\},\]

em que \(\mu \in \textbf{R}\) é a média do logaritmo do tempo de ocorrência do evento de interesse, \(\sigma > 0\) é o desvio-padrão e \(t \geq 0\).

Utilizando as relações entre as funções básicas estudadas anteriormente, as funções de sobrevivência (confiabilidade), de risco e de risco acumulado são

(t) = ()]

\[\lambda(t) = \frac{f(t)}{S(t)}\]

\[\Lambda(t) = -ln(S(t))\]

Para o ajuste do modelo lognormal, vamos utilizar a função flexsurvreg com dist="lnorm":

ajuste_lognorm<-flexsurvreg(Surv(tempos, cens) ~ 1, dist = "lnorm", data=dados1)
ajuste_lognorm

## Call:
## flexsurvreg(formula = Surv(tempos, cens) ~ 1, data = dados1, 
##     dist = "lnorm")
## 
## Estimates: 
##          est    L95%   U95%   se   
## meanlog  2.717  2.372  3.063  0.176
## sdlog    0.765  0.544  1.075  0.133
## 
## N = 20,  Events: 17,  Censored: 3
## Total time at risk: 347
## Log-likelihood = -65.7399, df = 2
## AIC = 135.4798

Assim, \(\hat{\mu}=2,717\) é o estimador de máxima verossimilhança para \(\mu\) e \(\hat{\sigma}=0,765\) é o estimador de máxima verossimilhança para \(\sigma\).

Os comandos a seguir fazem os gráficos das funções de sobrevivência, risco e risco acumulado para este modelo, comparando com a estimação não-paramétrica:

plot(ajuste_lognorm,xlab="Tempos (meses)",ylab="Função de Sobrevivênca S(t)")

plot(ajuste_lognorm,xlab="Tempos (meses)",ylab=expression("Função de Risco " ~ lambda(t) ),type="hazard")

plot(ajuste_lognorm,xlab="Tempos (meses)",ylab=expression("Função de Risco Acumulado " ~ Lambda(t) ),type="cumhaz")

Os comandos a seguir calculam, na sequência, o tempo médio de ocorrência do evento de interesse, o tempo mediano de ocorrência do evento de interesse e o percentil 80 do tempo de ocorrência do evento de interesse

Mean_lognorm<-summary(ajuste_lognorm,type="mean",cl=0.95,se=T)
Mean_lognorm

##  
##        est      lcl      ucl       se
## 1 20.28026 13.96182 32.47712 4.719672

Median_lognorm<-summary(ajuste_lognorm,type="median",cl=0.95,se=T)
Median_lognorm

##  
##        est      lcl      ucl       se
## 1 15.13751 10.51528 20.67908 2.670368

quantile80_lognorm<-summary(ajuste_lognorm,type="quantile",
                        quantiles=0.8,cl=0.95,se=T)
quantile80_lognorm

##  
##   quantile    est      lcl      ucl       se
## 1      0.8 28.814 19.57879 44.53052 6.323311

Já o comando a seguir fornece o tempo médio residual de ocorrência do evento de interesse no instante t=20 meses.

Resid20_lognorm<-reslifefsr(ajuste_lognorm,20)
Resid20_lognorm

## [1] 17.15706

3.4 DISTRIBUIÇÃO LOGLOGÍSTICA

Ao assumir que o tempo \(T\) possui uma distribuição loglogística (\(T\sim loglogística(\alpha,\gamma)\)), a sua função densidade de probabilidade é dada por:

\[f(t)=\frac{\gamma}{\alpha^\gamma}t^{\gamma-1}\left(1+(t/\alpha)^\gamma\right)^{-2}\]

em que \(\alpha > 0\) é o parâmetro de escala, \(\gamma > 0\) é o parâmetro de forma e \(t \geq 0\).

Utilizando as relações entre as funções básicas estudadas anteriormente, as funções de sobrevivência (confiabilidade), de risco e de risco acumulado são

\[S(t)=\frac{1}{1+(t/\alpha)^\gamma}\]

\[\lambda(t)=\frac{\gamma(t/\alpha)^{\gamma-1}}{\alpha\left[1+(t/\alpha)^\gamma\right]}\]

\[\Lambda(t) = -ln(S(t))\]

Para o ajuste do modelo loglogístico, vamos utilizar a função flexsurvreg com dist="llogi":

ajuste_loglogi<-flexsurvreg(Surv(tempos, cens) ~ 1, dist = "llogi", data=dados1)
ajuste_loglogi

## Call:
## flexsurvreg(formula = Surv(tempos, cens) ~ 1, data = dados1, 
##     dist = "llogi")
## 
## Estimates: 
##        est    L95%   U95%   se   
## shape   2.22   1.51   3.28   0.44
## scale  15.45  10.85  21.98   2.78
## 
## N = 20,  Events: 17,  Censored: 3
## Total time at risk: 347
## Log-likelihood = -66.03053, df = 2
## AIC = 136.0611

Assim, \(\hat{\gamma}=2,22\) é o estimador de máxima verossimilhança para \(\gamma\) e \(\hat{\alpha}=15,45\) é o estimador de máxima verossimilhança para \(\alpha\).

Os comandos a seguir fazem os gráficos das funções de sobrevivência, risco e risco acumulado para este modelo, comparando com a estimação não-paramétrica:

plot(ajuste_loglogi,xlab="Tempos (meses)",ylab="Função de Sobrevivênca S(t)")

plot(ajuste_loglogi,xlab="Tempos (meses)",ylab=expression("Função de Risco " ~ lambda(t) ),type="hazard")

plot(ajuste_loglogi,xlab="Tempos (meses)",ylab=expression("Função de Risco Acumulado " ~ Lambda(t) ),type="cumhaz")

Os comandos a seguir calculam, na sequência, o tempo médio de ocorrência do evento de interesse, o tempo mediano de ocorrência do evento de interesse e o percentil 80 do tempo de ocorrência do evento de interesse:

Mean_loglogi<-summary(ajuste_loglogi,type="mean",cl=0.95,se=T)
Mean_loglogi

##  
##        est      lcl      ucl       se
## 1 22.09722 14.85464 41.00807 6.997423

Median_loglogi<-summary(ajuste_loglogi,type="median",cl=0.95,se=T)
Median_loglogi

##  
##        est      lcl      ucl       se
## 1 15.44711 11.19229 21.26528 2.700148

quantile80_loglogi<-summary(ajuste_loglogi,type="quantile",
                        quantiles=0.8,cl=0.95,se=T)
quantile80_loglogi

##  
##   quantile      est      lcl      ucl       se
## 1      0.8 28.81187 19.20576 44.78709 6.688949

Já o comando a seguir fornece o tempo médio residual de ocorrência do evento de interesse no instante t=20 meses.

Resid20_loglogi<-reslifefsr(ajuste_loglogi,20)
Resid20_loglogi

## [1] 21.69073

3.5 DISTRIBUIÇÃO GAMMA

Ao assumir que o tempo \(T\) possui uma distribuição Gama (\(T\sim gama(\alpha,\beta)\)), a sua função densidade de probabilidade é dada por:

\[\frac{\beta^\alpha}{\Gamma(\alpha)}t^{\alpha-1}exp\{-\beta t\}\]

em que \(\alpha>0\) é o parâmetro de forma, \(\beta >0\) é o parâmetro de taxa e \(t\geq0\).

Utilizando as relações entre as funções básicas estudadas anteriormente, as funções de sobrevivência (confiabilidade), de risco e de risco acumulado são

\[S(t) = \int_{t}^\infty\frac{\beta^\alpha}{\Gamma(\alpha)}u^{\alpha-1}exp\{-\beta u\}du\]

\[\lambda(t) = \frac{f(t)}{S(t)}\]

\[\Lambda(t) = -ln(S(t))\]

Para o ajuste do modelo gama, vamos utilizar a função flexsurvreg com dist="gama":

ajuste_gama<-flexsurvreg(Surv(tempos, cens) ~ 1, dist = "gamma", data=dados1)
ajuste_gama

## Call:
## flexsurvreg(formula = Surv(tempos, cens) ~ 1, data = dados1, 
##     dist = "gamma")
## 
## Estimates: 
##        est     L95%    U95%    se    
## shape  2.1439  1.1580  3.9692  0.6737
## rate   0.1115  0.0540  0.2303  0.0413
## 
## N = 20,  Events: 17,  Censored: 3
## Total time at risk: 347
## Log-likelihood = -65.84754, df = 2
## AIC = 135.6951

Assim, \(\hat{\alpha}=2,14\) é o estimador de máxima verossimilhança para \(\alpha\) e \(\hat{\beta}=0,11\) é o estimador de máxima verossimilhança para \(\beta\).

Os comandos a seguir fazem os gráficos das funções de sobrevivência, risco e risco acumulado para este modelo, comparando com a estimação não-paramétrica:

plot(ajuste_gama,xlab="Tempos (meses)",ylab="Função de Sobrevivênca S(t)")

plot(ajuste_gama,xlab="Tempos (meses)",ylab=expression("Função de Risco " ~ lambda(t) ),type="hazard")

plot(ajuste_gama,xlab="Tempos (meses)",ylab=expression("Função de Risco Acumulado " ~ Lambda(t) ),type="cumhaz")

Os comandos a seguir calculam, na sequência, o tempo médio de ocorrência do evento de interesse, o tempo mediano de ocorrência do evento de interesse e o percentil 80 do tempo de ocorrência do evento de interesse:

Mean_gama<-summary(ajuste_gama,type="mean",cl=0.95,se=T)
Mean_gama

##  
##        est      lcl      ucl       se
## 1 19.22915 14.11851 26.05764 3.139777

Median_gama<-summary(ajuste_gama,type="median",cl=0.95,se=T)
Median_gama

##  
##        est      lcl      ucl       se
## 1 16.33618 11.67002 22.54571 2.781996

quantile80_gama<-summary(ajuste_gama,type="quantile",
                        quantiles=0.8,cl=0.95,se=T)
quantile80_gama

##  
##   quantile      est     lcl      ucl       se
## 1      0.8 28.55305 20.3679 41.10753 5.192332

Já o comando a seguir fornece o tempo médio residual de ocorrência do evento de interesse no instante t=20 meses.

Resid20_gama<-reslifefsr(ajuste_gama,20)
Resid20_gama

## [1] 12.22874

3.6 DISTRIBUIÇÃO GAMA GENERALIZADA

Ao assumir que o tempo \(T\) possui uma distribuição Gama Generalizada (\(T\sim gama(\alpha,\beta)\)), a sua função densidade de probabilidade é dada por:

\[f(t)=\frac{\gamma}{\Gamma(\kappa)\alpha^{\gamma\kappa}}t^{\gamma\kappa-1}exp\{-(\alpha t)^\gamma\}\]

Para esta distribuição, tem-se um parâmetro de escala \(\alpha\) e dois de forma \(\gamma\) e \(\kappa\), o que a torna bastante flexível

Para o ajuste do modelo gama generalizado, vamos utilizar a função flexsurvreg com dist="gengamma.orig":

ajuste_gengama<-flexsurvreg(Surv(tempos, cens) ~ 1, dist = "gengamma.orig", data=dados1)
ajuste_gengama

## Call:
## flexsurvreg(formula = Surv(tempos, cens) ~ 1, data = dados1, 
##     dist = "gengamma.orig")
## 
## Estimates: 
##        est       L95%      U95%      se      
## shape  4.14e-01  2.15e-03  7.99e+01  1.11e+00
## scale  5.52e-02  8.03e-26  3.80e+22  1.55e+00
## k      1.07e+01  4.57e-04  2.51e+05  5.50e+01
## 
## N = 20,  Events: 17,  Censored: 3
## Total time at risk: 347
## Log-likelihood = -65.69348, df = 3
## AIC = 137.387

Assim, \(\hat{\alpha}=0,052\) é o estimador de máxima verossimilhança para \(\alpha\), \(\hat{\gamma}=0,414\) é o estimador de máxima verossimilhança para \(\gamma\) e \(\hat{\kappa}=1,07\) é o estimador de máxima verossimilhança para \(\kappa\).

Os comandos a seguir fazem os gráficos das funções de sobrevivência, risco e risco acumulado para este modelo, comparando com a estimação não-paramétrica:

plot(ajuste_gengama,xlab="Tempos (meses)",ylab="Função de Sobrevivênca S(t)")

plot(ajuste_gengama,xlab="Tempos (meses)",ylab=expression("Função de Risco " ~ lambda(t) ),type="hazard")

plot(ajuste_gengama,xlab="Tempos (meses)",ylab=expression("Função de Risco Acumulado " ~ Lambda(t) ),type="cumhaz")

Os comandos a seguir calculam, na sequência, o tempo médio de ocorrência do evento de interesse, o tempo mediano de ocorrência do evento de interesse e o percentil 80 do tempo de ocorrência do evento de interesse:

#Mean_gengama<-summary(ajuste_gengama,type="mean",cl=0.95,se=T)
Mean_gengama<-NA
Median_gengama<-summary(ajuste_gengama,type="median",cl=0.95,se=T)
Median_gengama

##  
##       est lcl ucl  se
## 1 15.6301   0 Inf NaN

quantile80_gengama<-summary(ajuste_gengama,type="quantile",
                        quantiles=0.8,cl=0.95,se=T)
quantile80_gengama

##  
##   quantile      est     lcl ucl  se
## 1      0.8 28.55103 26.1183 Inf NaN

Já o comando a seguir fornece o tempo médio residual de ocorrência do evento de interesse no instante t=20 meses.

Resid20_gengama<-reslifefsr(ajuste_gengama,20)
Resid20_gengama

## [1] 14.38189

3.7 COMPARAÇÃO ENTRE OS MODELOS

Para encerrar a aula de hoje, vamos comparar, graficamente, as funções de sobrevivência, risco e risco acumulado, dos 6 modelos. lwd.ci=0, lty.ci=0)

plot(ajuste_exp,xlab="Tempos (meses)",ylab="Função de Sobrevivênca S(t)",lwd.ci=0, lty.ci=0, col="red")
lines(ajuste_wei,lwd.ci=0, lty.ci=0, col="blue")
lines(ajuste_lognorm,lwd.ci=0, lty.ci=0, col="green")
lines(ajuste_loglogi,lwd.ci=0, lty.ci=0, col="orange")
lines(ajuste_gama,lwd.ci=0, lty.ci=0, col="purple")
lines(ajuste_gengama,lwd.ci=0, lty.ci=0, col="gray")
legend(30,1,legend=c("Exponencial","Weibull","Log-Normal","Log-Logísitca","Gama","Gama-Generalizada"),fill=c("red","blue","green","orange","purple","gray"))

plot(ajuste_exp,xlab="Tempos (meses)",ylab="Função de Sobrevivênca S(t)",lwd.ci=0, lty.ci=0, col="red",type="hazard",ylim=c(0,0.2))
lines(ajuste_wei,lwd.ci=0, lty.ci=0, col="blue",type="hazard")
lines(ajuste_lognorm,lwd.ci=0, lty.ci=0, col="green",type="hazard")
lines(ajuste_loglogi,lwd.ci=0, lty.ci=0, col="orange",type="hazard")
lines(ajuste_gama,lwd.ci=0, lty.ci=0, col="purple",type="hazard")
lines(ajuste_gengama,lwd.ci=0, lty.ci=0, col="gray",type="hazard")
legend(0,0.2,legend=c("Exponencial","Weibull","Log-Normal","Log-Logísitca","Gama","Gama-Generalizada"),fill=c("red","blue","green","orange","purple","gray"))

plot(ajuste_exp,xlab="Tempos (meses)",ylab="Função de Sobrevivênca S(t)",lwd.ci=0, lty.ci=0, col="red",type="cumhaz")
lines(ajuste_wei,lwd.ci=0, lty.ci=0, col="blue",type="cumhaz")
lines(ajuste_lognorm,lwd.ci=0, lty.ci=0, col="green",type="cumhaz")
lines(ajuste_loglogi,lwd.ci=0, lty.ci=0, col="orange",type="cumhaz")
lines(ajuste_gama,lwd.ci=0, lty.ci=0, col="purple",type="cumhaz")
lines(ajuste_gengama,lwd.ci=0, lty.ci=0, col="gray",type="cumhaz")
legend(0,4,legend=c("Exponencial","Weibull","Log-Normal","Log-Logísitca","Gama","Gama-Generalizada"),fill=c("red","blue","green","orange","purple","gray"))

Na sequência, os comandos a seguir trazem uma comparação de resumos numéricos atrelados à algumas quantidades de interesse:

Resumo<-cbind(c(Mean_exp[[1]][1],Mean_wei[[1]][1],Mean_lognorm[[1]][1],Mean_loglogi[[1]][1],Mean_gama[[1]][1],Mean_gengama[[1]][1]),
              c(Median_exp[[1]][1],Median_wei[[1]][1],Median_lognorm[[1]][1],Median_loglogi[[1]][1],Median_gama[[1]][1],Median_gengama[[1]][1]),
              c(quantile80_exp[[1]][2],quantile80_wei[[1]][2],quantile80_lognorm[[1]][2],quantile80_loglogi[[1]][2],quantile80_gama[[1]][2],quantile80_gengama[[1]][2]),
              c(Resid20_exp,Resid20_wei,Resid20_lognorm,Resid20_loglogi,Resid20_gama,Resid20_gengama))
colnames(Resumo)<-c("Média","Mediana","Percentil 80","Residual em 20")
rownames(Resumo)<-c("Exponencial","Weibull","Log-Normal","Log-Logística","Gama","Gama Generalizada")

Resumo

##                   Média    Mediana  Percentil 80 Residual em 20
## Exponencial       20.41176 14.14836 32.85147     20.41176      
## Weibull           19.20078 16.82823 29.04556     11.57553      
## Log-Normal        20.28026 15.13751 28.814       17.15706      
## Log-Logística     22.09722 15.44711 28.81187     21.69073      
## Gama              19.22915 16.33618 28.55305     12.22874      
## Gama Generalizada NA       15.6301  28.55103     14.38189

3.8 EXERCÍCIOS

1- Para a base de dados do tratamento quimioterápico, ajuste todos os 6 modelos vistos hoje.

• Compare os gráficos da função de sobrevivência, risco e risco acumulado.
• Compare os tempos médio e mediano de ocorrência do evento de interesse.
• Compare as estimativas da função de sobrevivência no tempo 270 dias. Dica: use o comando summary(ajuste,t=270)

2- Para a base de dados de isolante elétrico, ajuste todos os 6 modelos vistos hoje.

• Compare os gráficos da função de confiabilidade, risco e risco acumulado.
• Compare os tempos médio e mediano de ocorrência do evento de interesse.
• Se a garantia é de 500 horas, compare as estimativas das probabilidades de falha até este instante.