ANÁLISE DE SOBREVIVÊNCIA E CONFIABILIDADE

AULA 5: MODELOS DE REGRESSÃO PARAMÉTRICOS (UMA COVARIÁVEL)

\[\\[0.05in]\]

1 PACOTES NECESSÁRIOS

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

require(flexsurv)
require(reslife)
require(dplyr)
options(scipen = 999)

• 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.

• Pacote dplyr: Responsável por permitir as manipulações em bases de dados.

• options(scipen = 999) : Comando responsável por evitar a notação científica.

2 BASES DE DADOS

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

2.1 LEUCEMIA AGUDA

Considere os tempos de sobrevivência, em semanas, de 17 pacientes com leucemia aguda (Lawless, 2003) apresentados na Tabela abaixo. Para esses pacientes, suas contagens de glóbulos brancos (WBC) foram registradas na data do diagnóstico e estas, com seus correspondentes logaritmos, na base 10, encontram-se também na Tabela:

2.2 MALÁRIA

Um estudo experimental realizado com camundongos para verificar a eficácia da imunização pela malária foi conduzido no Centro de Pesquisas René Rachou, Fiocruz, MG. Nesse estudo, quarenta e quatro camundongos foram aleatoriamente divididos em três grupos e todos foram infectados pela malária (Plasmodium berguei).

Os camundongos do grupo 1 foram imunizados 30 dias antes da infecção. Além da infecção pela malária, os camundongos dos grupos 1 e 3 foram, também, infectados pela esquistossomose (Schistossoma mansoní). A resposta de interesse nesse estudo foi o tempo decorrido desde a infecção pela malária até a morte do camundongo. Este tempo foi medido em dias e o estudo foi acompanhado por 30 dias. Os tempos de sobrevivência observados para os três grupos encontram-se abaixo. 0 simbolo + indica censura.

2.3 AIDS

O estudo é constituído de pacientes portadores de HIV atendidos entre 1986 e 2000 no Instituto de Pesquisa Clínica Evandro Chagas (Ipec/Fiocruz). Dessa coorte, obteve-se uma amostra de 156 indivíduos que foram diagnosticados como portadores de AIDS (critério CDC 1993) durante o período de acompanhamento. O objetivo é mensurar o tempo, em dias, do diagnóstico até a morte dos pacientes. Para comparação da sobrevivência, vamos ver como o tipo de tratamento impacta a sobrevivência: terapia antirretroviral: A = nenhum, B = mono, C = combinada, D = potente. A base está armazenada no objeto AIDS.csv

3 INTRODUÇÃO

Os estudos na área médica muitas vezes envolvem covariáveis que podem estar relacionadas com o tempo de sobrevivência. A insuficiência renal tem alta prevalência em todo o mundo, estimando-se cerca de um milhão de pacientes utilizando alguma forma de diálise. Qual seria o efeito da idade em que o paciente inicia a diálise sobre o tempo de sobrevivência?

Existem diferenças no tempo de sobrevivência conforme a doença de base causadora da insuficiência renal? É possível estimar o efeito da idade na sobrevivência controlado pela doença de base? Todas essas questões passam pela inclusão de covariáveis na modelagem das funções de sobrevivência (confiabilidade). No caso da AIDS, por exemplo, a contagem de células CD4 e CD8 ao diagnóstico são duas covariáveis que a literatura médica mostra serem importantes fatores de prognóstico para o tempo até a ocorrência de AIDS em pacientes infectados pelo HIV.

Em confiabilidade, pode ser de interesse avaliar como o tempo de falha de um determinado dispositivo pode ser afetado de acordo com quantidades físicas, como temperatura, pressão, voltagem, etc… As técnicas não-paramétricas apresentadas anteriormente não permitem a inclusão direta de covariáveis na análise. Estas técnicas foram importantes para descrever os dados de sobrevivência (confiabilidade) pela sua simplicidade e facilidade de aplicação, pois não envolvem nenhuma estrutura paramétrica. No entanto, este fato inviabiliza uma análise mais elaborada incluindo covariáveis.

Uma forma simples de se concluir isto é dividir os dados em estratos de acordo com as covariáveis e usar as técnicas não-paramétricas. A simplicidade dos cálculos e a facilidade de entendimento são as grandes vantagens da análise estratificada. No entanto, ela apresenta sérias limitações. A mais importante é que uma análise envolvendo várias covariáveis produz um número muito grande de estratos que podem conter poucas observações, ou talvez nenhuma. Isto faz com que as comparações fiquem impossíveis de serem realizadas.

A forma mais eficiente de acomodar o efeito dessas covariáveis é utilizar um modelo de regressão apropriado para dados censurados. Nas aulas passadas, já introduzimos os modelos paramétricos, porém não tínhamos acrescentado ainda o efeito de covariáveis. Em análise de sobrevivência e confiabilidade, existem duas classes de modelos de regressão propostos na literatura: os modelos paramétricos e os modelos semiparamétricos. Os modelos paramétricos, são mais eficientes, porém menos flexíveis do que os modelos semiparamétricos.

A segunda classe de modelos, também denominada simplesmente de modelo de regressão de Cox, tem sido bastante utilizada em estudos clínicos. Além da flexibilidade, este modelo permite incorporar facilmente covariáveis dependentes do tempo, que ocorrem com frequência em várias áreas de aplicação. A aula de hoje é dedicada ao estudo dos modelos de regressão paramétricos. Os modelos de regressão semiparamétricos serão estudados na parte final da disciplina.

Para os dados de sobrevivência (confiabilidade), considere \(T_i\) a variável aleatória que responde pelo tempo até a ocorrência do evento de interesse do elemento \(i\). Além disso, a essa variável aleatória, estão associadas as seguintes funções:

• Função densidade de probabilidade \(f_i(t|\boldsymbol{x_i})\)
• Função de sobrevivência (confiabilidade) \(S_i(t|\boldsymbol{x_i})\)
• Função de risco \(\lambda_i(t|\boldsymbol{x_i})\)

O subscrito \(i\) é justificado pelo fato de que cada elemento agora tem suas covariáveis \(\boldsymbol{x_i}\), para \(i=1,2,...,n\). A função de verossimilhança é:

\[L(\boldsymbol{\theta})\propto\prod_{i=1}^n\left[f_i(t_i)\right]^{\delta_i}\left[S_i(t_i)\right]^{1-\delta_i}\]

4 AJUSTE DOS MODELOS PARAMÉTRICOS DE REGRESSÃO

Esta seção é dedicada a relembrar a notação dos modelos de regressão paramétricos para dados de sobrevivência e confiabilidade e a conexão com os comandos do R. Para o ajuste de todos os modelos, será utilizada a função flexsurvreg com a mesma estrutura vista na aula passada. A grande diferença será na substituição do argumento Surv(tempos, cens) ~ 1 por Surv(tempos, cens) ~ covariaveis.

Para avaliar a significância dos coeficientes \(\beta\) estimados, vamos utilizar a função tidy(nome_do_modelo).

4.1 MODELO DE REGRESSÃO EXPONENCIAL

No contexto com covariáveis, a variável \(T_i\) possui a distribuição \(T_i\sim Exponencial(\alpha_i)\), para todo \(i=1,2,...,n\), em que \(\alpha_i\) irá variar de acordo com os valores das covariáveis observadas para o elemento \(i\). Para introduzir covariáveis, nós temos que considerar que o parâmetro \(\alpha_i\) é sempre positivo. A solução natural é permitir que as covariáveis entrem no modelo exponencialmente:

\[\alpha_i =exp(\boldsymbol{X_i}\boldsymbol{\boldsymbol{\beta}})\]

A função flexsurvreg para o modelo exponencial solta as seguintes saídas:

• rate: É o valor de \(\exp\{\hat{\beta_0}\}\). Para obter o valor de \(\hat{\beta_0}\), basta tomar o logaritmo de rate.

• \(\hat{\beta_1}\), \(\hat{\beta_2}\), …, \(\hat{\beta_k}\): Estimativas para os demais coeficientes regressores.

4.2 MODELO DE REGRESSÃO WEIBULL

No contexto com covariáveis, a variável \(T_i\) possui a distribuição \(T_i\sim Weibul(\alpha_i,\gamma)\), para todo \(i=1,2,...,n\), em que \(\alpha_i\) irá variar de acordo com os valores das covariáveis observadas para o elemento \(i\). Para introduzir covariáveis, nós temos que considerar que o parâmetro \(\alpha_i\) é sempre positivo. A solução natural é permitir que as covariáveis entrem no modelo exponencialmente:

\[\alpha_i =exp(\boldsymbol{X_i}\boldsymbol{\boldsymbol{\beta}})\]

A função flexsurvreg para o modelo Weibull solta as seguintes saídas:

• scale: É o valor de \(\exp\{\hat{\beta_0}\}\). Para obter o valor de \(\hat{\beta_0}\), basta tomar o logaritmo de scale.

• shape: É o valor de \(\hat{\gamma}\).

• \(\hat{\beta_1}\), \(\hat{\beta_2}\), …, \(\hat{\beta_k}\): Estimativas pros demais coeficientes regressores.

4.3 MODELO DE REGRESSÃO LOGNORMAL

No contexto com covariáveis, a variável \(T_i\) possui a distribuição \(T_i\sim LogNormal(\mu_i,\gamma)\), para todo \(i=1,2,...,n\), em que \(\mu_i\) irá variar de acordo com os valores das covariáveis observadas para o elemento \(i\).

Para introduzir covariáveis, nós temos que considerar que o parâmetro \(\mu_i\in \mathbb{R}\). A solução natural é permitir que as covariáveis entrem no modelo linearmente:

\[\mu_i =\boldsymbol{X_i}\boldsymbol{\boldsymbol{\beta}}\]

A função flexsurvreg para o modelo Log-Normal solta as seguintes saídas:

• meanlog: É o valor de \(\hat{\beta_0}\).

• sdlog: É o valor de \(\hat{\sigma}\).

• \(\hat{\beta_1}\), \(\hat{\beta_2}\), …, \(\hat{\beta_k}\): Estimativas pros demais coeficientes regressores.

4.4 AIC e BIC

Para comparar os diferentes modelos ajustados, podemos utilizar o AIC e BIC, vistos na aula passada.

5 ANÁLISE DE RESÍDUOS

Os resíduos de Cox-Snell auxiliam a examinar o ajuste global do modelo. Esses resíduos são quantidades determinadas por: \[\hat{e_i}=\hat{\Lambda}(t_i|\boldsymbol{x_i})\] em que \(\hat{\Lambda}(\cdot)\) é a função de risco acumulado obtida do modelo ajustado. Os resíduos \(\hat{e_i}\) devem seguir uma distribuição exponencial padrão, isto é, exponencial com parâmetro \(\alpha=1\).

No R, esses resíduos são extraidos com a função coxsnell_flexsurvreg(nome_do_modelo)

Podemos fazer dois gráficos:

• Gráfico 1: O gráfico de \(\hat{e_i}\) versus \(\hat{\Lambda}(\hat{e_i}|\boldsymbol{x_i})\) (Obtido a partir do Kaplan-Meier) deve ser uma reta com inclinação 1.

• Gráfico 2: O gráfico das curvas de sobrevivência desses resíduos, obtidas por Kaplan-Meier e pelo modelo exponencial padrão. Quanto mais próximas elas se apresentarem, melhor é considerado o ajuste do modelo aos dados.

6 INTERPRETAÇÃO COEFICIENTES

Uma proposta razoável de interpretação é a de se fazer uso da razão de tempos medianos (Hosmer e Lemeshow, 1999)}.

6.1 VARIÁVEL QUALITATIVA

Pode-se mostrar para uma covariável binária \(x_j\) que a razão de tempos medianos é:

\[\frac{t_{0,5}(x=1,\hat{\boldsymbol{\beta}})}{t_{0,5}(x=0,\hat{\boldsymbol{\beta}})}=e^{\hat{\beta_j}}\]

Isto é, fazer o exponencial da estimativa para \(\beta_j\) já retorna essa interpretação. Essa análise se estende à variáveis qualitativas com mais de duas categorias.

IMPORTANTE: Por causa da parametrização que estamos utilizando para a distribuição exponencial, devemos fazer, para este caso, \(e^{-\hat{\beta_j}}\)

6.2 VARIÁVEL QUANTITATIVA

Para uma variável contínua \(x_j\), pode-se mostrar que a razão de tempos medianos é: \[\frac{t_{0,5}(x=c+1,\hat{\boldsymbol{\beta}})}{t_{0,5}(x=c,\hat{\boldsymbol{\beta}})}=e^{\hat{\beta_j}}\]

Isto é, fazer o exponencial da estimativa para \(\beta_j\) já retorna essa interpretação. Ela traz o impacto na razão de tempos medianos ao aumentar a variável contínua em uma unidade.

IMPORTANTE: Por causa da parametrização que estamos utilizando para a distribuição exponencial, devemos fazer, para este caso, \(e^{-\hat{\beta_j}}\)

7 UMA VARIÁVEL QUANTITATIVA

Vamos ilustrar esse cenário com a base de dados de leucemia aguda. Para isso, vamos carregar a base de dados:

temp<-c(65,156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65)
 cens<-rep(1,17)
 log10_wbc<-c(3.36,2.88,3.63,3.41,3.78,4.02,4.00,4.23,3.73,3.85,3.97,
         4.51,4.54,5.00,5.00,4.72,5.00)
 dados1<-data.frame(temp,cens,log10_wbc)

7.1 MODELO DE REGRESSÃO EXPONENCIAL

Segue o código para o ajuste desse modelo:

ajuste_exp<-flexsurvreg(data=dados1,Surv(temp,cens)~log10_wbc,dist="exponential")
ajuste_exp

## Call:
## flexsurvreg(formula = Surv(temp, cens) ~ log10_wbc, data = dados1, 
##     dist = "exponential")
## 
## Estimates: 
##            data mean   est         L95%        U95%        se        
## rate               NA  0.00020810  0.00000727  0.00595491  0.00035611
## log10_wbc  4.09588235  1.10929791  0.29870332  1.91989251  0.41357627
##            exp(est)    L95%        U95%      
## rate               NA          NA          NA
## log10_wbc  3.03222876  1.34810960  6.82022534
## 
## N = 17,  Events: 17,  Censored: 0
## Total time at risk: 1062
## Log-likelihood = -83.87705, df = 2
## AIC = 171.7541

tidy(ajuste_exp)

## # A tibble: 2 × 5
##   term      estimate std.error statistic  p.value
##   <chr>        <dbl>     <dbl>     <dbl>    <dbl>
## 1 rate      0.000208  0.000356     NA    NA      
## 2 log10_wbc 1.11      0.414         2.68  0.00731

De acordo com o modelo ajustado, temos que:

\(exp\{\hat{\beta_0}\}=0,0002\rightarrow \hat{\beta_0}=ln(0,0002)=-8,477\) e \(\hat{\beta_1}=1,109\)

Assim, para 2 pacientes (i=1,2) com log10_wbc=3 e log10_wbc=4, temos que:

\[\hat{\alpha_1}=exp\{\hat{\beta_0}+\hat{\beta_1}x_{i1}\}=exp{-8,477+1,109*3}=0,0058\]

\[\hat{\alpha_2}=exp\{\hat{\beta_0}+\hat{\beta_1}x_{i2}\}=exp{-8,477+1,109*4}=0,0175\]

Portanto, \(T_1\sim Exp(\hat{\alpha_1}=0,0058)\) e \(T_1\sim Exp(\hat{\alpha_2}=0,0175)\)

Com essas funções, podemos obter todas as demais para estes dois pacientes (função de sobrevivência, risco, …).

Nesse modelo, o coeficiente \(\hat{\beta_1}\) é significativamente diferente de 0.

7.2 MODELO DE REGRESSÃO WEIBULL

Os seguintes comandos fazem o ajuste desse modelo:

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

## Call:
## flexsurvreg(formula = Surv(temp, cens) ~ log10_wbc, data = dados1, 
##     dist = "weibull")
## 
## Estimates: 
##            data mean   est         L95%        U95%        se        
## shape              NA       1.022       0.688       1.517       0.206
## scale              NA    4632.134     162.697  131881.406    7914.679
## log10_wbc       4.096      -1.098      -1.917      -0.280       0.418
##            exp(est)    L95%        U95%      
## shape              NA          NA          NA
## scale              NA          NA          NA
## log10_wbc       0.333       0.147       0.756
## 
## N = 17,  Events: 17,  Censored: 0
## Total time at risk: 1062
## Log-likelihood = -83.87136, df = 3
## AIC = 173.7427

tidy(ajuste_wei)

## # A tibble: 3 × 5
##   term      estimate std.error statistic  p.value
##   <chr>        <dbl>     <dbl>     <dbl>    <dbl>
## 1 shape         1.02     0.206     NA    NA      
## 2 scale      4632.    7915.        NA    NA      
## 3 log10_wbc    -1.10     0.418     -2.63  0.00853

De acordo com o modelo ajustado, temos que:

\(exp\{\hat{\beta_0}\}=4632.13\rightarrow \hat{\beta_0}=ln(4632.134)=8,441\), \(\hat{\beta_1}=-1.109\) e \(\hat{\gamma}=1,021\)

Assim, para 2 pacientes (i=1, 2) com log10_wbc=3 e log10_wbc=4, temos que:

\[\hat{\alpha_1}=exp\{\hat{\beta_0}+\hat{\beta_1}x_{i1}\}=exp\{8,441-1.109*3\}=165,892\]

\[\hat{\alpha_2}=exp\{\hat{\beta_0}+\hat{\beta_1}x_{i2}\}=exp\{8,441-1.109*4\}=54,681) \]

Portanto, \(T_1\sim Weibull(\hat{\alpha_1}=165,892 ; \hat{\gamma}=1,021)\) e \(T_2\sim Weibull(\hat{\alpha_2}=54,681 ; \hat{\gamma}=1,021)\)

Com essas funções, podemos obter todas as demais para estes dois pacientes (função de sobrevivência, risco, …).

Nesse modelo, o coeficiente \(\hat{\beta_1}\) é significativamente diferente de 0.

7.3 MODELO DE REGRESSÃO LOGNORMAL

O ajuste desse modelo aos dados é:

ajuste_lognorm<-flexsurvreg(data=dados1,Surv(temp,cens)~log10_wbc,dist="lognormal")
ajuste_lognorm

## Call:
## flexsurvreg(formula = Surv(temp, cens) ~ log10_wbc, data = dados1, 
##     dist = "lognormal")
## 
## Estimates: 
##            data mean  est      L95%     U95%     se       exp(est)  L95%   
## meanlog         NA    11.0738   7.3193  14.8284   1.9156       NA        NA
## sdlog           NA     1.1576   0.8272   1.6201   0.1985       NA        NA
## log10_wbc   4.0959    -1.8829  -2.7897  -0.9761   0.4626   0.1521    0.0614
##            U95%   
## meanlog         NA
## sdlog           NA
## log10_wbc   0.3768
## 
## N = 17,  Events: 17,  Censored: 0
## Total time at risk: 1062
## Log-likelihood = -83.7589, df = 3
## AIC = 173.5178

tidy(ajuste_lognorm)

## # A tibble: 3 × 5
##   term      estimate std.error statistic    p.value
##   <chr>        <dbl>     <dbl>     <dbl>      <dbl>
## 1 meanlog      11.1      1.92      NA    NA        
## 2 sdlog         1.16     0.199     NA    NA        
## 3 log10_wbc    -1.88     0.463     -4.07  0.0000470

De acordo com o modelo ajustado, temos que:

\(\hat{\beta_0}=11,073\), \(\hat{\beta_1}=-1.882\) e \(\hat{\sigma}=1,157\)

Assim, para 2 pacientes (i=1 e 2) com log10_wbc=3 e log10_wbc=4, temos que:

\[\hat{\mu_1}=\hat{\beta_0}+\hat{\beta_1}x_{i1}= 11,073-1,882*3=5,425\]

\[\hat{\mu_2}=\hat{\beta_0}+\hat{\beta_1}x_{i1}= 11,073-1,882*4=3.542\]

Portanto, \(T_1\sim LogNormal(\hat{\mu_1}=5,425,\hat{\sigma}=1,157)\) e \(T_2\sim LogNormal(\hat{\mu_1}=3.542,\hat{\sigma}=1,157)\)

Com essas funções, podemos obter todas as demais para estes dois pacientes (função de sobrevivência, risco, …).

Nesse modelo, o coeficiente \(\hat{\beta_1}\) é significativamente diferente de 0.

7.4 COMPARAÇÃO DOS MODELOS

Os comandos a seguir calculam o AIC e o BIC para estes modelos:

AIC<-c(ajuste_exp$AIC,
       ajuste_wei$AIC,
       ajuste_lognorm$AIC)
BIC<-c(ajuste_exp$BIC,
       ajuste_wei$BIC,
       ajuste_lognorm$BIC)
Modelos<-c("Modelo Exponencial", "Modelo Weibull",
           "Modelo Log-Normal")
Comparacao<-data.frame(Modelos,AIC,BIC)
Comparacao

##              Modelos      AIC      BIC
## 1 Modelo Exponencial 171.7541 173.4205
## 2     Modelo Weibull 173.7427 176.2424
## 3  Modelo Log-Normal 173.5178 176.0174

O modelo Exponencial se ajusta melhor aos dados segundo ambos os critérios.

7.5 ANÁLISE DE RESÍDUOS

O seguinte comando extrai os resíduos de Cox-Snell:

CoxSnell <-coxsnell_flexsurvreg(ajuste_exp)$est

Podemos fazer dois gráficos:

• Gráfico 1:

ekm<- survfit (Surv(CoxSnell,dados1$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")

• Gráfico 2:

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")

Ambos os gráficos apontam para uma boa qualidade de ajuste.

7.6 INTERPRETAÇÃO DOS COEFICIENTES

Ao retornar no modelo exponencial, temos que:

\[e^{-\hat{\beta_1}}=e^{-1,10929} = 0,3297\]

Assim, temos a seguinte interpretação:

Ao aumentar em uma unidade a variável log10_wbc o tempo mediano de sobrevivência é diminuído em aproximadamente 67%(100%-33%). Logo, a variável log10_wbc é um fator de risco para a sobrevivência em pacientes com leucemia agura.

Os comandos a seguir fazem os gráficos de sobrevivência para duas pessoas: uma com log10_wbc=3 e a outra com log10_wbc=4.

plot(ajuste_exp,newdata = list(log10_wbc=3),type="survival",
     col="red",col.obs="white", ci=F, xlab="Semanas",ylab="Função de sobrevivência estimada")
lines(ajuste_exp,newdata = list(log10_wbc=4),type="survival",
     col="green",col.obs="white", ci=F)
legend(5,0.3,col=c("red","green"),legend=c("log10_wbc=3","log10_wbc=4"),lty=1)

Os comandos a seguir fazem os gráficos da função de risco para duas pessoas: uma com log10_wbc=3 e a outra com log10_wbc=4.

plot(ajuste_exp,newdata = list(log10_wbc=3),type="hazard",
     col="red",col.obs="white", ci=F, xlab="Semanas",ylab="Função de risco estimada")
lines(ajuste_exp,newdata = list(log10_wbc=4),type="hazard",
     col="green",col.obs="white", ci=F)
legend(5,0.04,col=c("red","green"),legend=c("log10_wbc=3","log10_wbc=4"),lty=1)

Os comandos a seguir calculam os tempos médios e medianos de sobrevivência dos pacientes com leucemia:

##### Media
media1<-summary(ajuste_exp,newdata = list(log10_wbc=3),type="mean")
media1

## log10_wbc=3 
##        est      lcl      ucl
## 1 172.3632 57.92481 497.3538

media2<-summary(ajuste_exp,newdata = list(log10_wbc=4),type="mean")
media2

## log10_wbc=4 
##        est      lcl      ucl
## 1 56.84373 35.12009 92.05879

##### Mediana
median1<-summary(ajuste_exp,newdata = list(log10_wbc=3),type="median")
median1

## log10_wbc=3 
##        est      lcl      ucl
## 1 119.4731 42.30748 318.5014

median2<-summary(ajuste_exp,newdata = list(log10_wbc=4),type="median")
median2

## log10_wbc=4 
##        est      lcl      ucl
## 1 39.40107 25.11928 63.91947

8 UMA VARIÁVEL QUALITATIVA

Para uma única variável qualitativa, vamos considerar a base de dados de malária com a variável explicativa grupos. Para isso, vamos carregar a base de dados:

tempos<-c(7,8,8,8,8,12,12,17,18,22,30,30,30,30,30,30,8,8,9,10,10,14,
          15,15,18,19,21,22,22,23,25,8,8,8,8,8,8,9,10,10,10,11,17,19)
cens<-c(rep(1,10), rep(0,6),rep(1,15),rep(1,13))
grupos<-c(rep("Grupo 1",16), rep("Grupo 2",15), rep("Grupo 3",13))
dados2<-data.frame(tempos,cens,grupos)

Como a variável é qualitativa, é bem interessante começar a análise fazendo um teste de comparação de curvas de sobrevivência, que já aprendemos. Esse passo é feito para verificar se a variável tem potencial para participar do modelo final. Seguem os comandos:

#### Logrank
survdiff(Surv(tempos,cens)~grupos, data=dados2)

## Call:
## survdiff(formula = Surv(tempos, cens) ~ grupos, data = dados2)
## 
##                 N Observed Expected (O-E)^2/E (O-E)^2/V
## grupos=Grupo 1 16       10    17.00    2.8816    6.4111
## grupos=Grupo 2 15       15    14.51    0.0167    0.0317
## grupos=Grupo 3 13       13     6.49    6.5190   10.4447
## 
##  Chisq= 12.6  on 2 degrees of freedom, p= 0.002

Assim, há evidências de que alguma curva de sobrevivência se difere das demais em algum instante do tempo. As comparações múltiplas podem ser utilizadas para detectar essa diferença, mas pularemos esta parte (já foi visto).

8.1 MODELO DE REGRESSÃO EXPONENCIAL

Segue o ajuste para estes dados:

ajuste_exp<-flexsurvreg(data=dados2,Surv(tempos,cens)~grupos,dist="exponential")
ajuste_exp

## Call:
## flexsurvreg(formula = Surv(tempos, cens) ~ grupos, data = dados2, 
##     dist = "exponential")
## 
## Estimates: 
##                data mean  est      L95%     U95%     se       exp(est)  L95%   
## rate                NA     0.0333   0.0179   0.0620   0.0105       NA        NA
## gruposGrupo 2   0.3409     0.6328  -0.1674   1.4329   0.4082   1.8828    0.8459
## gruposGrupo 3   0.2955     1.0683   0.2439   1.8927   0.4206   2.9104    1.2762
##                U95%   
## rate                NA
## gruposGrupo 2   4.1910
## gruposGrupo 3   6.6373
## 
## N = 44,  Events: 38,  Censored: 6
## Total time at risk: 673
## Log-likelihood = -143.8657, df = 3
## AIC = 293.7315

tidy(ajuste_exp)

## # A tibble: 3 × 5
##   term          estimate std.error statistic p.value
##   <chr>            <dbl>     <dbl>     <dbl>   <dbl>
## 1 rate            0.0333    0.0105     NA    NA     
## 2 gruposGrupo 2   0.633     0.408       1.55  0.121 
## 3 gruposGrupo 3   1.07      0.421       2.54  0.0111

De acordo com o modelo ajustado, temos que:

\(exp\{\hat{\beta_0}\}=0.033\rightarrow \hat{\beta_0}=ln(0.033)=-3,411\), \(\hat{\beta_{grupo2}}=0,633\) e \(\hat{\beta_{grupo3}}=1,068\)

Assim, para 3 ratos (i=1,2 e 3) que pertencem, respectivamente, aos grupos 1, 2 e 3, tem-se que

\[\hat{\alpha_1}=exp\{\hat{\beta_0}+\hat{\beta_{grupo2}}x_{i1}+\hat{\beta_{grupo3}}x_{i2}\}=exp\{-3,411\}=0,033\]

\[\hat{\alpha_2}=exp\{\hat{\beta_0}+\hat{\beta_{grupo2}}x_{i1}+\hat{\beta_{grupo3}}x_{i2}\}=exp\{-3,411+0,633\}=0,061\]

\[\hat{\alpha_3}=exp\{\hat{\beta_0}+\hat{\beta_{grupo2}}x_{i1}+\hat{\beta_{grupo3}}x_{i2}\}=exp\{-3,411+1,068\}=0,096\]

Portanto, \(T_1\sim Exp(\hat{\alpha_1}=0,033)\), \(T_2\sim Exp(\hat{\alpha_2}=0,061)\) e \(T_3\sim Exp(\hat{\alpha_3}=0,096)\)

Com essas funções, podemos obter todas as demais para estes dois pacientes (função de sobrevivência, risco, …).

Conforme o modelo ajustado, notase que apenas o efeito associado ao grupo 3 (\(\hat{\beta_{grupo3}}\)) é significamente diferente de 0.

8.2 MODELO DE REGRESSÃO WEIBULL

Segue o ajuste para estes dados:

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

## Call:
## flexsurvreg(formula = Surv(tempos, cens) ~ grupos, data = dados2, 
##     dist = "weibull")
## 
## Estimates: 
##                data mean  est      L95%     U95%     se       exp(est)  L95%   
## shape               NA     2.3609   1.8250   3.0542   0.3101       NA        NA
## scale               NA    26.5290  20.4037  34.4931   3.5533       NA        NA
## gruposGrupo 2   0.3409    -0.4299  -0.7702  -0.0895   0.1736   0.6506    0.4629
## gruposGrupo 3   0.2955    -0.8714  -1.2225  -0.5204   0.1791   0.4183    0.2945
##                U95%   
## shape               NA
## scale               NA
## gruposGrupo 2   0.9144
## gruposGrupo 3   0.5943
## 
## N = 44,  Events: 38,  Censored: 6
## Total time at risk: 673
## Log-likelihood = -129.2443, df = 4
## AIC = 266.4886

tidy(ajuste_wei)

## # A tibble: 4 × 5
##   term          estimate std.error statistic     p.value
##   <chr>            <dbl>     <dbl>     <dbl>       <dbl>
## 1 shape            2.36      0.310     NA    NA         
## 2 scale           26.5       3.55      NA    NA         
## 3 gruposGrupo 2   -0.430     0.174     -2.48  0.0133    
## 4 gruposGrupo 3   -0.871     0.179     -4.87  0.00000114

De acordo com o modelo ajustado, temos que:

\(exp\{\hat{\beta_0}\}=26.529\rightarrow \hat{\beta_0}=ln(26.529)=3,278\), \(\hat{\beta_{grupo2}}=-0,429\), \(\hat{\beta_{grupo3}}=-0,871\) e \(\hat{\gamma}=2,361\)

Assim, para 3 ratos (i=1,2 e 3) que pertencem, respectivamente, aos grupos 1, 2 e 3, tem-se que

\[\hat{\alpha_1}=exp\{\hat{\beta_0}+\hat{\beta_{grupo2}}x_{i1}+\hat{\beta_{grupo3}}x_{i2}\}=exp\{3,278\}=26,523\]

\[\hat{\alpha_2}=exp\{\hat{\beta_0}+\hat{\beta_{grupo2}}x_{i1}+\hat{\beta_{grupo3}}x_{i2}\}=exp\{3,278-0,429\}=17,270\]

\[\hat{\alpha_3}=exp\{\hat{\beta_0}+\hat{\beta_{grupo2}}x_{i1}+\hat{\beta_{grupo3}}x_{i2}\}=exp\{3,278-0,871\}=11,100\]

Portanto, \(T_1\sim Weibull(\hat{\alpha_1}=26,523 ; \hat{\gamma}=2,361)\), \(T_2\sim Weibull(\hat{\alpha_2}=17,270 ; \hat{\gamma}=2,361)\) e \(T_3\sim Weibull(\hat{\alpha_3}=11,100 ; \hat{\gamma}=2,361)\)

Com essas funções, podemos obter todas as demais para estes dois pacientes (função de sobrevivência, risco, …).

Conforme o modelo ajustado, nota-se que os efeitos associados ao grupo 2 (\(\hat{\beta_{grupo2}}\)) e ao grupo 3(\(\hat{\beta_{grupo3}}\)) são significamente diferentes de 0.

8.3 MODELO DE REGRESSÃO LOGNORMAL

Segue o ajuste para estes dados:

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

## Call:
## flexsurvreg(formula = Surv(tempos, cens) ~ grupos, data = dados2, 
##     dist = "lnorm")
## 
## Estimates: 
##                data mean  est      L95%     U95%     se       exp(est)  L95%   
## meanlog             NA     2.8754   2.6180   3.1329   0.1314       NA        NA
## sdlog               NA     0.5015   0.3969   0.6337   0.0598       NA        NA
## gruposGrupo 2   0.3409    -0.1800  -0.5415   0.1816   0.1845   0.8353    0.5819
## gruposGrupo 3   0.2955    -0.5864  -0.9614  -0.2114   0.1913   0.5563    0.3824
##                U95%   
## meanlog             NA
## sdlog               NA
## gruposGrupo 2   1.1991
## gruposGrupo 3   0.8094
## 
## N = 44,  Events: 38,  Censored: 6
## Total time at risk: 673
## Log-likelihood = -128.5089, df = 4
## AIC = 265.0179

tidy(ajuste_lognorm)

## # A tibble: 4 × 5
##   term          estimate std.error statistic  p.value
##   <chr>            <dbl>     <dbl>     <dbl>    <dbl>
## 1 meanlog          2.88     0.131     NA     NA      
## 2 sdlog            0.502    0.0598    NA     NA      
## 3 gruposGrupo 2   -0.180    0.184     -0.976  0.329  
## 4 gruposGrupo 3   -0.586    0.191     -3.07   0.00218

De acordo com o modelo ajustado, temos que:

\(\hat{\beta_0}=2,875\), \(\hat{\beta_{grupo2}}=-0.179\), \(\hat{\beta_{grupo3}}=-0.586\) e \(\hat{\sigma}=0,501\)

Assim, para 3 ratos (i=1,2 e 3) que pertencem, respectivamente, aos grupos 1, 2 e 3, tem-se que

\[\hat{\mu_1}=\hat{\beta_0}+\hat{\beta_{grupo2}}x_{i1}+\hat{\beta_{grupo3}}x_{i2}= 2,875\]

\[\hat{\mu_2}=\hat{\beta_0}+\hat{\beta_{grupo2}}x_{i1}+\hat{\beta_{grupo3}}x_{i2}= 2,875-0,179=2,696\]

\[\hat{\mu_3}=\hat{\beta_0}+\hat{\beta_{grupo2}}x_{i1}+\hat{\beta_{grupo3}}x_{i2}= 2,875-0,586=2,289\]

Portanto, \(T_1\sim LogNormal(\hat{\mu_1}=2,875,\hat{\sigma}=0,501)\), \(T_2\sim LogNormal(\hat{\mu_2}=2,696,\hat{\sigma}=0,501)\) e \(T_3\sim LogNormal(\hat{\mu_3}=2,289,\hat{\sigma}=0,501)\)

Com essas funções, podemos obter todas as demais para estes dois pacientes (função de sobrevivência, risco, …).

Conforme o modelo ajustado, nota-se que apenas o efeito associado ao grupo 3 (\(\hat{\beta_{grupo3}}\)) é significamente diferente de 0.

8.4 COMPARAÇÃO DOS MODELOS

Os comandos a seguir calculam o AIC e o BIC para este modelo:

AIC<-c(ajuste_exp$AIC,
       ajuste_wei$AIC,
       ajuste_lognorm$AIC)
BIC<-c(ajuste_exp$BIC,
       ajuste_wei$BIC,
       ajuste_lognorm$BIC)
Modelos<-c("Modelo Exponencial", "Modelo Weibull",
           "Modelo Log-Normal")
Comparacao<-data.frame(Modelos,AIC,BIC)
Comparacao

##              Modelos      AIC      BIC
## 1 Modelo Exponencial 293.7315 299.0841
## 2     Modelo Weibull 266.4886 273.6254
## 3  Modelo Log-Normal 265.0179 272.1546

Logo, o modelo Log-Normal desponta como o mais adequado a esses dados.

8.5 ANÁLISE DE RESÍDUOS

O seguinte comando extrai os resíduos de Cox-Snell:

CoxSnell <-coxsnell_flexsurvreg(ajuste_lognorm)$est

Podemos fazer dois gráficos*:

• Gráfico 1:

ekm<- survfit (Surv(CoxSnell,dados2$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")

• Gráfico 2:

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")

Ambos os gráficos apontam para uma boa qualidade de ajuste.

8.6 INTERPRETAÇÃO DOS COEFICIENTES

Ao retornar no modelo lognormal, temos que:

ajuste_lognorm

## Call:
## flexsurvreg(formula = Surv(tempos, cens) ~ grupos, data = dados2, 
##     dist = "lnorm")
## 
## Estimates: 
##                data mean  est      L95%     U95%     se       exp(est)  L95%   
## meanlog             NA     2.8754   2.6180   3.1329   0.1314       NA        NA
## sdlog               NA     0.5015   0.3969   0.6337   0.0598       NA        NA
## gruposGrupo 2   0.3409    -0.1800  -0.5415   0.1816   0.1845   0.8353    0.5819
## gruposGrupo 3   0.2955    -0.5864  -0.9614  -0.2114   0.1913   0.5563    0.3824
##                U95%   
## meanlog             NA
## sdlog               NA
## gruposGrupo 2   1.1991
## gruposGrupo 3   0.8094
## 
## N = 44,  Events: 38,  Censored: 6
## Total time at risk: 673
## Log-likelihood = -128.5089, df = 4
## AIC = 265.0179

Como vimos, apenas o efeito associado ao grupo 3 foi considerado significativo pelo teste de Wald. Assim, temos a seguinte interpretação:

Com respeito ao grupo 1, o tempo mediano de sobrevivência do grupo 3 é diminuído em aproximadamente 45%.

Os comandos a seguir fazem os gráficos de sobrevivência para para 3 ratos que pertencem, respectivamente, aos grupos 1, 2 e 3,

plot(ajuste_lognorm,newdata = list(grupos="Grupo 1"),type="survival",
     col="red",col.obs="white", ci=F, xlab="Semanas",ylab="Função de sobrevivência estimada",t=c(0:30))
lines(ajuste_lognorm,newdata = list(grupos="Grupo 2"),type="survival",
     col="green",col.obs="white", ci=F, t=c(0:30))
lines(ajuste_lognorm,newdata = list(grupos="Grupo 3"),type="survival",
     col="blue",col.obs="white", ci=F, t=c(0:30))

legend(2,0.3,col=c("red","green","blue"),legend=c("Grupo 1","Grupo 2", "Grupo 3"),lty=1)

Os comandos a seguir fazem os gráficos da função de risco para para 3 ratos que pertencem, respectivamente, aos grupos 1, 2 e 3,

plot(ajuste_lognorm,type="hazard",newdata = list(grupos="Grupo 1"), col="red",col.obs="white", ci=F, xlab="Semanas",ylab="Função de risco estimada", xlim=c(0,30), t=c(0:30))
lines(ajuste_lognorm,newdata = list(grupos="Grupo 2"),type="hazard",
     col="green",col.obs="white", ci=F, t=c(0:30))
lines(ajuste_lognorm,newdata = list(grupos="Grupo 3"),type="hazard",
     col="blue",col.obs="white", ci=F, t=c(0:30))
legend(15,0.04,col=c("red","green","blue"),legend=c("Grupo 1","Grupo 2", "Grupo 3"),lty=1)

Os comandos a seguir calculam os tempos médios e medianos de sobrevivência dos pacientes com leucemia:

##### Media
media1<-summary(ajuste_lognorm,newdata = list(grupos="Grupo 1"),type="mean")
media1

## grupos=Grupo 1 
##        est      lcl      ucl
## 1 20.10971 15.48984 26.27365

media2<-summary(ajuste_lognorm,newdata = list(grupos="Grupo 2"),type="mean")
media2

## grupos=Grupo 2 
##        est      lcl     ucl
## 1 16.79779 13.03515 21.7447

media3<-summary(ajuste_lognorm,newdata = list(grupos="Grupo 3"),type="mean")
media3

## grupos=Grupo 3 
##       est      lcl      ucl
## 1 11.1873 8.343218 15.29365

##### Mediana
median1<-summary(ajuste_lognorm,newdata = list(grupos="Grupo 1"),type="median")
median1

## grupos=Grupo 1 
##        est      lcl      ucl
## 1 17.73317 13.88561 23.07366

median2<-summary(ajuste_lognorm,newdata = list(grupos="Grupo 2"),type="median")
median2

## grupos=Grupo 2 
##        est      lcl      ucl
## 1 14.81265 11.55929 19.07105

median3<-summary(ajuste_lognorm,newdata = list(grupos="Grupo 3"),type="median")
median3

## grupos=Grupo 3 
##        est      lcl      ucl
## 1 9.865202 7.499062 13.06826

9 EXERCÍCIOS

1- Base de dados AIDS

1. Ajuste o modelo de regressão exponencial aos dados considerando a variável tratam como explicativa.

• Encontre a estimativa de \(\beta_0\)
• Investigue os coeficientes que foram significativos.
• Qual a distribuição estimada de um paciente que pertence ao tratamento B?

2. Ajuste o modelo de regressão weibull aos dados considerando a variável tratam como explicativa.

• Encontre a estimativa de \(\beta_0\) e \(\gamma\)
• Investigue os coeficientes que foram significativos.
• Qual a distribuição estimada de um paciente que pertence ao tratamento C?

3. Ajuste o modelo de regressão Log-Normal aos dados considerando a variável tratam como explicativa.

• Encontre a estimativa de \(\beta_0\) e \(\sigma\)
• Investigue os coeficientes que foram significativos.
• Qual a distribuição estimada de um paciente que pertence ao tratamento D?

4. Qual dos três modelos melhor se ajusta aos dados? Faça a análise de resíduos nesse modelo.

5. Qual a interpretação dos coeficientes estimados no modelo final?

6. Faça os gráficos da função de sobrevivência estimada e da função de risco para cada um dos tratamentos.