Respostas da lista 2 da disciplina Análise de Sobrevivência e Confiabilidade - GET00147, matéria ofertada ao curso de bacharelado em Estatística pela Univesidade Federal Fluminense no semestre 2024.1. - N terminado

Questão 1

Enunciado: Dentre os modelos probabilísticos em análise de sobrevivência e confiabilidade, vimos os modelos Exponencial, Weibull, LogNormal, LogLogístico e Gama. Esses modelos são utilizados para representar o tempo T até a ocorrência de um evento de interesse. Explique quais as características da função de risco em cada um desses modelos.

No modelo exponencial, (\(T \sim Exp(\alpha), \alpha > 0\)) a funcão de risco é constante. Sendo assim, o risco da ocorrência do desfecho de interesse em uma unidade elementar dado que ela não apresentou o evento até um instante de tempo, é o mesmo para todo tempo.

No modelo weibull (\(T \sim Weibull(\alpha, \gamma), \alpha > 0, \gamma>0\)), a função de risco é alterada principalmente pelo valor de \(\gamma\). De modo que, se \(\gamma < 1\), a função de risco é estritamente decrescente, ou seja, o risco do evento da ocorrência instantânea do evento de interesse decresce ao passar do tempo. Já se \(\gamma >1\), a função de risco é estritamente crescente, implicando que o risco da ocorrência do evento de interesse aumenta conforme o tempo. Observe também, que se \(\gamma = 0\), então \(T \sim Weibull(\alpha, \gamma) = Exp(\alpha)\), então o risco será constante em \(\alpha\).

O modelo log-normal (\(T \sim lognormal(\mu, \sigma), \mu \in \mathbb{R}, \sigma > 0\)) apresenta uma função de risco não monótona. Nesse modelo, a função de risco inicialmente cresce, alcançando um ponto máximo, após o qual passa a decrescer estritamente. Isso implica que a probabilidade de ocorrência do evento de interesse começa baixa, aumentando ao longo do tempo. Entretanto, em determinado momento, a probabilidade atinge um valor máximo e, a partir daí, o risco começa a decrescer. Em outras palavras, o risco inicialmente aumenta, estabiliza-se no ponto de máximo risco, e depois diminui à medida que o tempo avança.

Para a distribuição log-logistica, quando \(\gamma \leq 1\), a função de risco é descrescente, já quando \(\gamma > 1\), a função de risco torna-se similar ao da distribuição log-normal.

No modelo gamma, quando \(\alpha < 1\), a função de risco é estritamente decrescente. Se \(\alpha > 1\), a função de risco é estritamente crescente.

Questão 2

Enunciado:

alpha = 100
gamma = 2 

sobrevivencia = function(t, gamma, alpha){exp(-(t/alpha)^gamma)}

a) Qual a probabilidade de um rato sobreviver sem tumor aos primeiros 30 dias? E aos primeiros 45 dias?

#Probabilidade de um rato sobreviver a 30 dias:
sobrevivencia(30, gamma, alpha)
## [1] 0.9139312
#Probabilidade de um rato sobreviver a 45 dias
sobrevivencia(45, gamma, alpha)
## [1] 0.8166865

b) Qual é o tempo médio até o aparecimento do tumor?

#O tempo médio até o aparecimento do tumor em dias: 
alpha*gamma(1+(1/gamma)) 
## [1] 88.62269

c) tempo mediano até o aparecimento do tumor?

alpha*log(2)^(1/gamma)
## [1] 83.25546

d)

Para a função de risco de um modelo Weibull, onde \(T \sim Weibull(\alpha, \gamma), \alpha > 0, \gamma > 0\), temos:

#Funcao de risco: 

risk = function(t,alpha,gamma){gamma*(t^(gamma-1))/(alpha^gamma)}
#risco em 30:
risk(30, alpha = 100,gamma = 2)*100
## [1] 0.6

Logo, a probabilidade de um rator desenvolver um tumor instantâneamente no dia 30 dado que ele não desenvolveu até então é de 0,6%

risk(45, alpha = 100, gamma = 2)*100
## [1] 0.9
risk(60, alpha = 100, gamma = 2)*100
## [1] 1.2

Do mesmo modo, segue que a probabilidade de um rato apresentar o tumor dado que não apresentou até o tempo de 45 dias é 0,9%, e para 60 dias é de 1,2%.

Questão 3

Sabemos que se \(T \sim Exp(\alpha), \alpha > 0\), então para \(t \geq 0\),

\[ f(t) = \alpha e^{-\alpha t} \\ S(t) = e^{-\alpha t} \] Logo, a função de verossimilhança \(L(\alpha)\),

\[ L(\alpha) = \prod_{i=1}^n f(t_i)^{\delta_i} \cdot S(t_i)^{1-\delta_i} \\ =\prod_{i=1}^n (\alpha e^{-\alpha t_i})^{\delta_i} \cdot(e^{-\alpha t_i})^{1-\delta_i} \\ = \prod_{i=1}^n \alpha^{\delta_i} \cdot e^{-\alpha t_i} \\ = \alpha^{\sum_{i=1}^n \delta_i} \cdot e^{-\alpha \sum_{i=1}^n t_i} \] Seja \(l(\alpha) = ln(L(\alpha))\), segue que

\[ l(\alpha) = {\sum_{i=1}^n \delta_i}\cdot log(\alpha) - \alpha \sum_{i=1}^n t_i \] Queremos achar o argumento que maxima a funcão de verossimilhança. Para isso, derivamos \(l(\alpha)\) em \(\alpha\):

\[ \frac{dl(\alpha)}{d\alpha} = \frac{{\sum_{i=1}^n \delta_i}}{\alpha} - \sum_{i=1}^n t_i = 0 \\ \iff \sum_{i=1}^n \delta_i - \alpha \sum_{i=1}^n t_i = 0 \\ \iff \hat\alpha = \frac{\sum_{i=1}^n \delta_i}{\sum_{i=1}^n t_i} = \frac{r}{\sum_{i=1}^n t_i}, \] tal que r são ocorrências dos eventos de interesse.

Questão 4

Critério de Informação de Akaike (AIC)

O Critério de Informação de Akaike (AIC) é dado pela seguinte fórmula:

\[ \text{AIC} = 2k - 2\ln(\hat{L}) \]

onde: - \(k\) é o número de parâmetros no modelo, - \(\hat{L}\) é o valor da função de verossimilhança maximizada.

Critério de Informação Bayesiano (BIC)

O Critério de Informação Bayesiano (BIC) é dado pela seguinte fórmula:

\[ \text{BIC} = k \ln(n) - 2\ln(\hat{L}) \]

onde: - \(k\) é o número de parâmetros no modelo, - \(n\) é o número de observações, - \(\hat{L}\) é o valor da função de verossimilhança maximizada.

Questão 6

O fabricante de um tipo de isolante elétrico quer conhecer o comportamento de seu produto funcionando a uma temperatura de 200°C. Um teste de vida foi realizado nessas condições, usando 60 isoladores elétricos. O teste terminou quando 45 deles haviam falhado (censura do tipo II). As 15 unidades que não haviam falhado ao final do teste foram, dessa forma, censuradas no tempo 𝑡= 2729.

#Construindo o dataset da questão
tempos <- c(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)
cens <- c(rep(1, 45), rep(0, 15))
dados <- data.frame(tempos, cens)
head(dados)
##   tempos cens
## 1    151    1
## 2    164    1
## 3    336    1
## 4    365    1
## 5    403    1
## 6    454    1

a) Ajuste o modelo probabilístico Weibull a esses dados.

#Carregando os pacotes necessários
require(flexsurv)
## Carregando pacotes exigidos: flexsurv
require(reslife)
## Carregando pacotes exigidos: reslife
## Carregando pacotes exigidos: pracma
## Carregando pacotes exigidos: gsl
## 
## Anexando pacote: 'gsl'
## Os seguintes objetos são mascarados por 'package:pracma':
## 
##     Ci, erf, erfc, eta, expint_E1, expint_Ei, fact, hypot, psi, Si,
##     sinc, zeta
ajuste_wei=flexsurvreg(Surv(tempos, cens) ~ 1, dist = "weibull", data=dados)
ajuste_wei
## Call:
## flexsurvreg(formula = Surv(tempos, cens) ~ 1, data = dados, dist = "weibull")
## 
## Estimates: 
##        est       L95%      U95%      se      
## shape     1.281     1.001     1.640     0.161
## scale  1993.215  1586.793  2503.733   231.902
## 
## N = 60,  Events: 45,  Censored: 15
## Total time at risk: 90799
## Log-likelihood = -385.697, df = 2
## AIC = 775.394

As estimativas para os parametros são \(\hat{\gamma} = 1,281\) e \(\hat{\alpha} = 1993,215\), que são obtidas pelo método da máximo verossimilhança composta. Na saida, também é possível observar o erro padrão estima \(se\) obtido utilizando resultados assintóticos da distribuição amostral dos parâmetros. Seja \(\theta\) \[ \mathbf{\hat{\theta}} \sim N \left( \mathbf{\theta}, \text{Var}(\mathbf{\hat{\theta}}) \right), \] tal que …

. Do mesmo modo, é construido o intervalo de confiança

Para o tempo mediano, sabemos que no caso da distribuição weibull, o tempo mediano \(t_{0.5}\) é dado por \(\alpha ln(2)^{\frac{1}{\gamma}}\). Logo, por invariância, \(\hat{t}_{0.5} = \hat{\alpha} \ln(2)^{\frac{1}{\hat{\gamma}}}\).

Median_wei<-summary(ajuste_wei,type="median",cl=0.95,se=T)
Median_wei
##  
##       est      lcl      ucl       se
## 1 1497.36 1170.148 1892.429 188.2431

Questao 8

dados = readr::read_csv("AIDS.csv")
## Rows: 193 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): tratam
## dbl (2): tempo, status
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
dados$tratam = factor(dados$tratam, levels = c("A", "B", 
                                                "C", "D"))
head(dados)
## # A tibble: 6 × 3
##   tempo status tratam
##   <dbl>  <dbl> <fct> 
## 1   852      1 B     
## 2   123      1 A     
## 3  1145      1 B     
## 4  2755      0 B     
## 5  2117      0 B     
## 6   329      0 A

Modelo exponencial

ajuste_exp = flexsurvreg(Surv(tempo, status) ~ tratam, dist = "exponential", data = dados)
resultado_exp = tidy(ajuste_exp)
resultado_exp
## # A tibble: 4 × 5
##   term    estimate std.error statistic   p.value
##   <chr>      <dbl>     <dbl>     <dbl>     <dbl>
## 1 rate     0.00216  0.000381     NA    NA       
## 2 tratamB -1.59     0.226        -7.07  1.58e-12
## 3 tratamC -2.68     0.445        -6.01  1.80e- 9
## 4 tratamD -3.01     1.02         -2.97  3.00e- 3

Para obter a estimavida de \(\beta_0\), temos que \(\exp(\beta_0) = 0.002156625\), então \(\hat{\beta_0} = ln(0.002156625) = -6.139\). Além disso, podemos ver que os coeficientes para o Tratamento B, C e D são estatisticamente significativos.

1 - exp(resultado_exp$estimate[2])
##   tratamB 
## 0.7968292

Logo, a mediana do tempo até óbito no tratamente monovalente é 79% menor do que sem tratamento. *Entao é melhor nao fazer tratamento?

1 - exp(resultado_exp$estimate[3])
##   tratamC 
## 0.9311423
1- exp(resultado_exp$estimate[4])
##   tratamD 
## 0.9508962

Para um paciente do tratamento B \(\hat{\alpha}_{grupoB} = \hat{\beta}_{0} + \hat{\beta}_{grupoB}\):

alpha_b = exp(-6.139 + resultado_exp$estimate[2])
alpha_b 
##      tratamB 
## 0.0004382555

logo, \(T \sim Exp(\hat{\alpha}_{grupoB} = 0.0004382555)\).

Modelo weibull

ajuste_wei = flexsurvreg(Surv(tempo, status) ~ tratam, dist = "weibull", data = dados)
resultado_wei = tidy(ajuste_wei)
resultado_wei
## # A tibble: 5 × 5
##   term    estimate std.error statistic   p.value
##   <chr>      <dbl>     <dbl>     <dbl>     <dbl>
## 1 shape      0.884    0.0766     NA    NA       
## 2 scale    450.      90.3        NA    NA       
## 3 tratamB    1.68     0.263       6.40  1.53e-10
## 4 tratamC    2.92     0.534       5.47  4.58e- 8
## 5 tratamD    3.37     1.18        2.86  4.25e- 3

Para obter a estimativaa de \(\hat\beta_0\), temos que \(\exp(\hat\beta_0) = 0.884\), então \(\hat{\beta}_0 = \ln(0.884) = -0.1232982\). Além disso, \(\hat{\beta}_B = 1.68\), \(\hat{\beta}_C = 2.92\), \(\hat{\beta}_D = 3.37\) e \(\hat{\gamma} = 450\)

Podemos ver que os coeficientes para o Tratamento B, C e D são estatisticamente significativos.

Para um paciente do tratamento C \(\hat{\alpha}_{C} = \hat{\beta}_{0} + \hat{\beta}_{C}\):, então:

alpha_c = exp(-0.1232982 + resultado_wei$estimate[4])
alpha_c 
##  tratamC 
## 16.38454

Logo, \(T \sim Weibull(\hat\alpha_c = 16.384, \hat\gamma = 450)\)

Modelo log-normal

ajuste_lnorm = flexsurvreg(Surv(tempo, status) ~ tratam, dist = "lnorm", data = dados)
resultado_lnorm  = tidy(ajuste_lnorm )
resultado_lnorm 
## # A tibble: 5 × 5
##   term    estimate std.error statistic   p.value
##   <chr>      <dbl>     <dbl>     <dbl>     <dbl>
## 1 meanlog     5.52     0.248     NA    NA       
## 2 sdlog       1.54     0.122     NA    NA       
## 3 tratamB     1.83     0.302      6.05  1.49e- 9
## 4 tratamC     3.12     0.471      6.64  3.18e-11
## 5 tratamD     2.63     0.719      3.65  2.59e- 4

Assim, temos \(\exp(\hat\beta_0) = 5.518\), \(\hat{\beta}_B = 1.825\), \(\hat{\beta}_C = 3.124\), \(\hat{\beta}_D = 2.626\) e \(\hat{\sigma} = 1.536948\)

Podemos ver que os coeficientes para o Tratamento B, C e D são estatisticamente significativos.

Para um paciente do tratamento D \(\hat{\mu}_{D} = \hat{\beta}_{0} + \hat{\beta}_{D}\):, então:

mu_D = 5.518 + resultado_wei$estimate[5]
mu_D 
##  tratamD 
## 8.885079

Logo, \(T \sim LogNormal(\hat\mu_c = 8.885, \hat\sigma = 1.537)\)

Comparação dos modelos

AIC<-c(ajuste_exp$AIC,
       ajuste_wei$AIC,
       ajuste_lnorm$AIC)
BIC<-c(ajuste_exp$BIC,
       ajuste_wei$BIC,
       ajuste_lnorm$BIC)
Modelos<-c("Modelo Exponencial", "Modelo Weibull",
           "Modelo Log-Normal")

comparacao = data.frame(Modelos, AIC,BIC)
comparacao
##              Modelos      AIC      BIC
## 1 Modelo Exponencial 1493.752 1506.803
## 2     Modelo Weibull 1493.616 1509.930
## 3  Modelo Log-Normal 1489.209 1505.522

Como queremos o modelo com menor AIC e BIC, escolhemos o modelo Log-Normal. Podemos, então, prosseguir com a análise de resíduos de Cox-Snell.

CoxSnell <-coxsnell_flexsurvreg(ajuste_lnorm)$est
#Gráficos
par(mfrow = c(2,1))
ekm<- survfit (Surv(CoxSnell,dados$status)~1,type="kaplan-meier")
plot(ekm, fun="cumhaz",xlab="Resíduos de Cox-Snell",ylab="Função de risco acumulado")
abline(0, 1, col="red")
plot(ekm, xlab="Resíduos de Cox-Snell",ylab="Função de sobrevivência")
CoxSnell<-sort(CoxSnell)
exp1<-exp(-CoxSnell)
points(CoxSnell,exp1, col="red",type="o")

Interpretação do modelo

Uma possível forma de interpretar os coeficientes estimados pelo modelo, é olhando a razão dos tempos medianos. De forma que, dado \(x_j\) uma covariável binária (0 se não possui, 1 se possui) do modelo, então seja \(t_{0,5}\) o tempo mediano, \[ \frac{t_{0,5}(x = 1, \hat{\beta})}{t_{0,5}(x = 0, \hat{\beta})} = \exp(\hat{\beta}) \]

ajuste_lnorm
## Call:
## flexsurvreg(formula = Surv(tempo, status) ~ tratam, data = dados, 
##     dist = "lnorm")
## 
## Estimates: 
##          data mean  est      L95%     U95%     se       exp(est)  L95%   
## meanlog       NA     5.5178   5.0313   6.0043   0.2482       NA        NA
## sdlog         NA     1.5369   1.3152   1.7960   0.1222       NA        NA
## tratamB   0.5181     1.8254   1.2336   2.4172   0.3019   6.2050    3.4335
## tratamC   0.1813     3.1245   2.2019   4.0470   0.4707  22.7477    9.0426
## tratamD   0.0725     2.6264   1.2174   4.0354   0.7189  13.8241    3.3785
##          U95%   
## meanlog       NA
## sdlog         NA
## tratamB  11.2139
## tratamC  57.2246
## tratamD  56.5656
## 
## N = 193,  Events: 90,  Censored: 103
## Total time at risk: 181080
## Log-likelihood = -739.6043, df = 5
## AIC = 1489.209

Olhando para o tratamento B, temos que \(\hat\beta = 6.2050\). Logo, o tempo mediano de morte de pessoas com HIV tratadas com terapia antirretroviral monovalente é 6 vezes maior do que o tempo mediano para os que não receberam nenhum tratamento.
Já no tratamento C, vemos que \(\hat\beta = 22.7477\). Sendo assim, os pacientes com HIV que receberam o tratamento antirretroviral combinada, apresentaram um tempo mediano até a morte aproximadamente 23 vezes maior do que os que não receberam. Vendo o tratamento D, segue que \(\hat\beta = 13.8241\). Logo, é esperado os pacientes que recebem o tratamento antirretroviral potente apresentam tempo mediano até a morte de aproximadamente 14 vezes maior do que os que não receberam nenhum tratamento.

Gráficos

plot(ajuste_lnorm, newdata = list(tratam = "A"), type = "survival",
     col = "red", col.obs = "white", ci = FALSE, xlab = "Dias", 
     ylab = "Função de Sobrevivência Estimada")
lines(ajuste_lnorm, newdata = list(tratam = "B"), type = "survival",
     col = "green", col.obs = "white", ci = FALSE)
lines(ajuste_lnorm, newdata = list(tratam = "C"), type = "survival",
     col = "blue", col.obs = "white", ci = FALSE)
lines(ajuste_lnorm, newdata = list(tratam = "D"), type = "survival",
     col = "purple", col.obs = "white", ci = FALSE)

legend("topright", legend = c("Nenhum", "Mono", "Combinada", "Potente"), 
       col = c("red", "green", "blue", "purple"), lty = 1)

plot(ajuste_lnorm,type="hazard",nwdata = list(tratam="A"), col="red",col.obs="white", ci=F, xlab="Dias",ylab="Função de risco estimada", ylim = c(0,0.004))
## Warning in plot.window(...): "nwdata" não é um parâmetro gráfico
## Warning in plot.xy(xy, type, ...): "nwdata" não é um parâmetro gráfico
## Warning in axis(side = side, at = at, labels = labels, ...): "nwdata" não é um
## parâmetro gráfico
## Warning in axis(side = side, at = at, labels = labels, ...): "nwdata" não é um
## parâmetro gráfico
## Warning in box(...): "nwdata" não é um parâmetro gráfico
## Warning in title(...): "nwdata" não é um parâmetro gráfico
lines(ajuste_lnorm,newdata = list(tratam="B"),type="hazard",
     col="green",col.obs="white", ci=F)
lines(ajuste_lnorm,newdata = list(tratam="C"),type="hazard",
     col="blue",col.obs="white", ci=F)
lines(ajuste_lnorm,newdata = list(tratam="D"),type="hazard",
     col="purple",col.obs="white", ci=F)

legend("topright", legend = c("Nenhum", "Mono", "Combinada", "Potente"), 
       col = c("red", "green", "blue", "purple"), lty = 1)

Questão 9

#Base de dados 
temp<-c(65,156,100,134,16,108,121,4,39,143,
56,26,22,1,1,5,65,56,65,17,7,16,22,
3,4,2,3,8,4,3,30,4,43)
cens<-c(rep(1,17),rep(1,16))
log10wbc<-c(3.36,2.88,3.63, 3.41,3.78,4.02,4,4.23,
3.73,3.85,3.97,4.51,4.54,5,5,4.72,5,3.64,
3.48,3.6,3.18,3.95,3.72,4,4.28,4.43,4.45,
4.49,4.41,4.32,4.90,5,5)
AG<-c(rep(0,17),rep(1,16))
dados<-data.frame(temp,cens,log10wbc,AG)
head(dados)
##   temp cens log10wbc AG
## 1   65    1     3.36  0
## 2  156    1     2.88  0
## 3  100    1     3.63  0
## 4  134    1     3.41  0
## 5   16    1     3.78  0
## 6  108    1     4.02  0

Escolha das covariáveis

Primeiro, podemos testar separadamente se as variáveis são significativas para o desfecho de estudo.

ajuste_wbc = survdiff(Surv(temp,cens)~log10wbc, data=dados,rho=0)
ajuste_ag = survdiff(Surv(temp,cens)~AG, data=dados,rho=0)

ajuste_wbc$pvalue
## [1] 0.01896658
ajuste_ag$pvalue
## [1] 0.003653063
mod1<-flexsurvreg(data=dados,Surv(temp,cens)~log10wbc,dist="gengamma.orig")
mod2<-flexsurvreg(data=dados,Surv(temp,cens)~AG,dist="gengamma.orig")
mod3<-flexsurvreg(data=dados,Surv(temp,cens)~AG+log10wbc,dist="gengamma.orig")
covar = c("log10wbc", "AG", "Ambas")

#Visualizando os ajustes
flexsurv::tidy(mod3)
## # A tibble: 5 × 5
##   term     estimate std.error statistic  p.value
##   <chr>       <dbl>     <dbl>     <dbl>    <dbl>
## 1 shape      0.234      0.121     NA    NA      
## 2 scale      0.0570     0.559     NA    NA      
## 3 k         13.1       13.4       NA    NA      
## 4 AG        -1.02       0.416     -2.45  0.0142 
## 5 log10wbc  -1.12       0.368     -3.05  0.00228
tibble::tibble(variaveis = covar, aic = c(mod1$AIC,mod2$AIC,mod3$AIC)  )
## # A tibble: 3 × 2
##   variaveis   aic
##   <chr>     <dbl>
## 1 log10wbc   306.
## 2 AG         299.
## 3 Ambas      302.

Por meio dos resultados obtidos, é interessante prosseguir tanto com o segundo modelo quanto com o primeiro. Escolhendo então o modelo que melhor se ajusta aos dados, ficamos com o modelo 2, que tem apenas a covariável AG, pois tem o menor AIC.

Escolha da distribuição

  ajustes <- list(
    expon = flexsurvreg(data=dados, Surv(temp,cens)~AG, dist="exp"),
    weibull = flexsurvreg(data=dados, Surv(temp,cens)~AG, dist="weibull"),
    lnorm = flexsurvreg(data=dados, Surv(temp,cens)~AG, dist="lnorm"),
    llogis = flexsurvreg(data=dados, Surv(temp,cens)~AG, dist="llogis"),
    gamma = flexsurvreg(data=dados, Surv(temp,cens)~AG, dist="gamma"),
    gengamma = flexsurvreg(data=dados,Surv(temp,cens)~AG,dist="gengamma.orig"),
    expon2 = flexsurvreg(data=dados,Surv(temp,cens)~AG+log10wbc, dist="exp")
  )
resultados <- tibble::tibble(
  nome_distribuicao = c("Exponencial", "Weibull", "Lognormal", "Log-logística", "Gamma", "Gengamma","exp2"),
  AIC = sapply(ajustes, function(x) x$AIC),
  BIC = sapply(ajustes, function(x) x$BIC)
)
resultados
## # A tibble: 7 × 3
##   nome_distribuicao   AIC   BIC
##   <chr>             <dbl> <dbl>
## 1 Exponencial        303.  306.
## 2 Weibull            304.  309.
## 3 Lognormal          309.  314.
## 4 Log-logística      310.  314.
## 5 Gamma              304.  309.
## 6 Gengamma           299.  305.
## 7 exp2               299.  304.

Por apresentar o menor AIC e BIC, será escolhido o modelo exponencial com duas variaveis .

Análise dos resíduos

CoxSnell <-coxsnell_flexsurvreg(ajustes$expon)$est
#
ekm<- survfit (Surv(CoxSnell,dados$cens)~1,type="kaplan-meier")
plot(ekm, fun="cumhaz",xlab="Resíduos de Cox-Snell",ylab="Função de risco acumulado")
abline(0, 1, col="red")

#
plot(ekm, xlab="Resíduos de Cox-Snell",ylab="Função de sobrevivência")
CoxSnell<-sort(CoxSnell)
exp1<-exp(-CoxSnell)
points(CoxSnell,exp1, col="red",type="o")

CoxSnell2 <-coxsnell_flexsurvreg(ajustes$expon2)$est
plot(ekm, xlab="Resíduos de Cox-Snell",ylab="Função de sobrevivência")
CoxSnell<-sort(CoxSnell)
exp1<-exp(-CoxSnell)
points(CoxSnell,exp1, col="red",type="o")

coeficientes

ajustes$expon2
## Call:
## flexsurvreg(formula = Surv(temp, cens) ~ AG + log10wbc, data = dados, 
##     dist = "exp")
## 
## Estimates: 
##           data mean  est       L95%      U95%      se        exp(est)  L95%    
## rate            NA   0.001081  0.000112  0.010471  0.001253        NA        NA
## AG        0.484848   1.018063  0.305328  1.730799  0.363647  2.767829  1.357070
## log10wbc  4.135758   0.700033  0.139238  1.260828  0.286125  2.013818  1.149397
##           U95%    
## rate            NA
## AG        5.645162
## log10wbc  3.528341
## 
## N = 33,  Events: 33,  Censored: 0
## Total time at risk: 1349
## Log-likelihood = -146.545, df = 3
## AIC = 299.0901

Como \(e^{\hat\beta_{AG}} = 2,7678\), temos que o tempo mediano da morte de pacientes com leucemia aguda que aprentam o antígeno Calla na superfice dos bastos é de aproximadamente 2,76 vezes maior do que os que não apresentam. Sendo \(e^{\hat\beta_{logWBC}} = 2,013\), então, a cada uma