ANÁLISE DE SOBREVIVÊNCIA E CONFIABILIDADE
AULA 8: MODELO DE COX ESTRATIFICADO
\[\\[0.05in]\]1 PACOTES NECESSÁRIOS
Para esta oitava aula prática de Análise de Sobrevivência e Confiabilidade, precisaremos do seguinte pacote:
require(survival)- Pacote
survival: É o pacote básico para o estudo de Análise de Sobrevivência e Confiabilidade. Contém as principais rotinas, incluindo definição de objetosSurv, curvas Kaplan-Meier, modelos Cox e modelos paramétricos.
2 BASES DE DADOS
2.1 LEUCEMIA PEDIÁTRICA
Com o objetivo de entender melhor quais fatores afetam o tempo de sobrevivência de uma criança brasileira com leucemia linfoblástica aguda, um grupo de 128 crianças, com idade inferior a 15 anos, foi acompanhado em alguns hospitais de Belo Horizonte. A variável resposta de interesse foi ○ tempo a partir da remissão (ausência da doença) até a recidiva ou morte (a que ocorrer primeiro)}. Das 128 crianças, 120 entraram em remissão, e são estas que formam o conjunto de dados em estudo. As covariáveis consideradas foram:
A base de dados está armazenada no objeto leucemia.txt.
Os comandos a seguir fazem a leitura da base:
dados<-read.table("leucemia.txt",header=T)
dados$leuinic <- factor(dados$leuinic)
dados$idadec <- factor(dados$idadec)
dados$zpesoc <- factor(dados$zpesoc)
dados$zestc <- factor(dados$zestc)
dados$pasc <- factor(dados$pasc)
dados$vacc <- factor(dados$vacc)
dados$riskc <- factor(dados$riskc)
dados$r6c <- factor(dados$r6c)
attach(dados)Para o ajuste dos modelos de Cox e Cox estratificado, consideraremos
a seleção das variáveis leuinic, idadec,
zpeso, pasc e vacc.
3 AJUSTE DO MODELO DE COX
Esta seção é dedicada ao ajuste do modelo de Cox tradicional.
3.1 ESTIMAÇÃO DOS PARÂMETROS \(\pmb{\beta}\)
Os comandos a seguir fazem o ajuste do modelo de cox tradicional:
mod1<-coxph(Surv(tempos,cens)~leuinic + idadec + zpesoc + pasc + vacc,data=dados, x = T,method="breslow")
mod1## Call:
## coxph(formula = Surv(tempos, cens) ~ leuinic + idadec + zpesoc +
## pasc + vacc, data = dados, x = T, method = "breslow")
##
## coef exp(coef) se(coef) z p
## leuinic1 1.1091 3.0316 0.3941 2.814 0.00489
## idadec1 0.7105 2.0351 0.3706 1.917 0.05519
## zpesoc1 -2.0554 0.1280 0.4963 -4.141 3.45e-05
## pasc1 -1.2247 0.2939 0.4559 -2.686 0.00723
## vacc1 1.3239 3.7581 0.4143 3.196 0.00139
##
## Likelihood ratio test=32.03 on 5 df, p=5.852e-06
## n= 103, number of events= 393.2 QUALIDADE DO AJUSTE
Os comandos a seguir constroem os resíduos de Schoenfeld:
Schoenfeld<-resid(mod1,type="scaledsch")
par(mfrow=c(2,3))
plot(cox.zph(mod1, transform="identity"))
Como observado, há indícios de que, para as variáveis
leuinic e pasc, temos problemas na verificação
da suposição de riscos proporcionais. Para confirmar isso, podemos
conduzir um teste de hipótese cuja hipótese nula é a proporcionalidade
dos riscos para cada uma das covariáveis:
cox.zph(mod1,transform = "identity",terms=F)## chisq df p
## leuinic1 10.287 1 0.0013
## idadec1 4.678 1 0.0305
## zpesoc1 2.837 1 0.0921
## pasc1 4.305 1 0.0380
## vacc1 0.175 1 0.6760
## GLOBAL 14.686 5 0.0118De acordo com os valores-p dessa tabela, a variável
Leuinic, pasc e idadec violam a
suposição de riscos proporcionais. Como essa suposição é principalmente
violada pela variável Leuinic, ajustaremos, na próxima
seção, um modelo de cox estratificado nessa variável
4 AJUSTE DO MODELO DE COX ESTRATIFICADO
Na análise dos dados de leucemia pediátrica, apresentada e discutida anteriormente, foi possível observar, quando da análise dos resíduos de Schoenfeld, uma indicação de violação da suposição de taxas de falha proporcionais para, em especial, a covariável LEUINI (contagem de leucócitos iniciais no sangue periférico). Uma possibilidade de análise desses dados seria, desse modo, estratificá-los de acordo com a covariável LEUINI, uma vez que a suposição de taxas de falha proporcionais pode não ser válida entre as crianças com LEUINI < 75000 mm³ e aquelas com LEUINI > 73000 mm³. Esta seção é dedicada ao ajuste do modelo de Cox estratificado.
4.1 ESTIMAÇÃO DOS PARÂMETROS \(\pmb{\beta}\)
O modelo de Cox estratificado em que as crianças são separadas em dois estratos distintos, de acordo com as categorias da covariável LEUINI fica expresso, nesse caso, por:
\[\lambda(t|\pmb{x_{ij}})=\lambda_{0_j}(t)exp\{\pmb{x_{ij}\beta}\}\]
para \(j=1,2\) e \(i=1,2,...,n_j\), em que \(n_j\) é o número de observações no \(j\)-ésimo estrato.
Para ajustar esse modelo assumindo que o vetor \(\pmb{\beta}\) é comum aos dois estratos,
fez-se uso dos mesmos comandos do modelo anterior, porem colocando
strata(leuinic) ao invés de leuinic:
mod2<-coxph(Surv(tempos,cens)~idadec + zpesoc + pasc + vacc + strata(leuinic),data=dados, x = T,method="breslow")
mod2## Call:
## coxph(formula = Surv(tempos, cens) ~ idadec + zpesoc + pasc +
## vacc + strata(leuinic), data = dados, x = T, method = "breslow")
##
## coef exp(coef) se(coef) z p
## idadec1 0.7993 2.2240 0.3836 2.084 0.03716
## zpesoc1 -2.4091 0.0899 0.5210 -4.624 3.77e-06
## pasc1 -1.2490 0.2868 0.4649 -2.686 0.00723
## vacc1 1.3589 3.8920 0.4194 3.240 0.00119
##
## Likelihood ratio test=28.67 on 4 df, p=9.107e-06
## n= 103, number of events= 39Comparando este modelo com o modelo anterior, é possível observar
resultados muito similares. Observe que a variável leunic
some, pois ela foi usada na estratificação.
4.2 O VETOR \(\pmb{\beta}\) É COMUM?
Note que de base de cada estrato, . Esta suposição pode ser testada usando-se, por exemplo, o teste da razão de verossimilhança, cuja estatística de teste é dada, nesse caso, por
\[TRV=-2\left[l(\hat{\pmb{\beta}})-\sum_{j=1}^ml_j(\hat{\pmb{\beta_j}})\right]\]
sendo \(l(\hat{\pmb{\beta}})\) o logaritmo da função de verossimilhança parcial sob o modelo que assume \(\pmb{\beta}\) comuns em cada estrato e \(l_j(\hat{\pmb{\beta_j}})\), O logaritmo da função de verossimilhança parcial sob o modelo que assume \(\pmb{\beta}\)’s distintos em cada estrato. Sob a hipótese nula de que os betas são comuns e para grandes amostras, a estatística TRV segue uma distribuição qui-quadrado com (m - 1)p graus de liberdade, em que m é o número de estratos e p a dimensão do vetor \(\pmb{\beta}\).
No nosso caso a TRV terá, sob H0, distribuição qui-quadrado com (m-1)p = (2-1)4=4 graus de liberdade. Os comandos a seguir fazem este teste:
## Construindo uma base de dados só com `leuinic=0`
leucc1<-as.data.frame(cbind(tempos[leuinic==0],cens[leuinic==0],idadec[leuinic==0],zpesoc[leuinic==0],pasc[leuinic==0],vacc[leuinic==0]))
## Construindo uma base de dados só com `leuinic=1`
leucc2<-as.data.frame(cbind(tempos[leuinic==1],cens[leuinic==1],idadec[leuinic==1],zpesoc[leuinic==1],pasc[leuinic==1],vacc[leuinic==1]))
## Ajustando o modelo de cox para `leuinic=0`
mod3<-coxph(Surv(V1,V2)~V3+V4+V5+V6,data=leucc1,x = T,method="breslow")
## Ajustando o modelo de cox para `leuinic=0`
mod4<-coxph(Surv(V1,V2)~V3+V4+V5+V6,data=leucc2,x = T,method="breslow")
## Construindo a estatística de teste
trv<-2*(-mod2$loglik[2]+mod3$loglik[2]+mod4$loglik[2])
## Cálculo do p-valor:
1-pchisq(trv,4)## [1] 0.2513484O resultado do teste mostra não haver evidencias de que os \(\pmb{\beta}\) sejam distintos entre os estratos.
4.3 QUALIDADE DO AJUSTE
Os comandos a seguir constroem os resíduos de Schoenfeld:
Schoenfeld<-resid(mod2,type="scaledsch")
par(mfrow=c(2,2))
plot(cox.zph(mod2, transform="identity",terms=F))
cox.zph(mod2, transform="identity",terms=F)## chisq df p
## idadec1 0.852963 1 0.356
## zpesoc1 0.143580 1 0.705
## pasc1 2.903555 1 0.088
## vacc1 0.000494 1 0.982
## GLOBAL 4.249793 4 0.373Dos resultados obtidos, nenhuma séria violação à suposição de riscos proporcionais é sugerido para as covariáveis consideradas no modelo.
4.4 INTERPRETAÇÃO DOS COEFICIENTES
Com o modelo final, podemos fazer as interpretações dos coeficientes:
mod2## Call:
## coxph(formula = Surv(tempos, cens) ~ idadec + zpesoc + pasc +
## vacc + strata(leuinic), data = dados, x = T, method = "breslow")
##
## coef exp(coef) se(coef) z p
## idadec1 0.7993 2.2240 0.3836 2.084 0.03716
## zpesoc1 -2.4091 0.0899 0.5210 -4.624 3.77e-06
## pasc1 -1.2490 0.2868 0.4649 -2.686 0.00723
## vacc1 1.3589 3.8920 0.4194 3.240 0.00119
##
## Likelihood ratio test=28.67 on 4 df, p=9.107e-06
## n= 103, number of events= 39Considerando as demais variáveis constantes,
O paciente com idade \(>\) 96 meses possui um risco 2.22 vezes maior que o do paciente com idade \(\leq\) 96 meses.
O risco do paciente com peso padronizado $ > -2$ é 92% menor que o do paciente com peso padronizado \(\leq -2\).
O risco do paciente com percentual de linflobastos medulares \(> 5\%\) é 72 % menor que o do paciente com percentual de linflobastos medulares \(\leq 5\%\).
O paciente com percentual de vacúolos \(> 15\%\) é 3,89 vezes maior que o risco do paciente com percentual de vacúolos \(\leq 15\%\)
4.5 E A VARIÁVEL
LEUINI?
O modelo de Cox estratificado não fornece uma estimativa do efeito da covariável usada para a estratificação, no caso, a contagem de leucócitos iniciais (LEUINI). Este fato não representa, contudo. uma limitação desse modelo, pois as funções de taxa de falha acumulada de base, bem como as funções de sobrevivência de base, fornecidas por este mesmo modelo para cada uma das categorias da covariável estratificadora, permitem uma avaliação indireta desse efeito. Com esta finalidade, e para o modelo ajustado, as funções de taxa de falha acumulada de base e de sobrevivência de base, para ambos os estratos, foram obtidas no R utilizando-se os comandos apresentados a seguir:
# risco acumulado de base
H0<-basehaz(mod2,centered=F)
# risco acumulado de base do estrato 1 (leuinic=0)
H01<-as.matrix(H0[1:81,1])
# risco acumulado de base do estrato 2 (leuinic=1)
H02<-as.matrix(H0[82:101,1])
# tempos do estrato 1
tempo1<- H0$time[1:81]
# sobrevivência de base estrato 1
S01<-exp(-H01)
# tempos do estrato 2
tempo2<- H0$time[82:101]
# sobrevivência de base estrato 2
S02<-exp(-H02)
par(mfrow=c(1,2))
plot(tempo2,H02,lty=4,type="s",xlab="Tempos",xlim=range(c(0,4)),ylab=expression(Lambda[0]* (t)))
lines(tempo1,H01,type="s",lty=1)
legend(0.0,9,lty=c(1,4),c("Leuini < 75000","Leuini > 75000"),lwd=1,bty="n",cex=0.8)
plot(c(0,tempo1),c(1,S01),lty=1,type="s",xlab="Tempos",xlim=range(c(0,4)),ylab="So(t)")
lines(c(0,tempo2),c(1,S02),lty=4,type="s")
legend(2.2,0.9,lty=c(1,4),c("Leuini < 75000","Leuini > 75000"),lwd=1,bty="n",cex=0.8)
Pode-se observar que o estrato com valores mais altos de leucometria inicial apresenta risco de recidiva ou morte maior entre crianças com leucemia.