Análise de Sobrevivencia - Lista 3
2024-06-19
Vamos estudar e avaliar os fatores associados a sobrevivencia dos pacientes a partir do momento que eles entram na UTI com dados provenientes de 200 pacientes internados com cancer no hospital INCA.
library('survival')
library('flexsurv')
library('survminer')## Loading required package: ggplot2## Loading required package: ggpubr##
## Attaching package: 'survminer'## The following object is masked from 'package:survival':
##
## myelomadados<-read.table("UTI.txt",header=T)
head(dados)## tempo status sexo idade gptumor desnut comorbi leucopenia
## 1 182 0 M 41 Loco nao nao nao
## 2 66 1 M 24 Loco nao nao nao
## 3 182 0 F 62 Loco nao nao nao
## 4 182 0 F 72 Loco nao nao nao
## 5 182 0 M 66 Hemato nao nao nao
## 6 7 1 M 47 Mtx nao nao naoAnálise não paramétrica
Vamos iniciar então análise não-paramétrica dos dados
Temos uma variável quantitativa na base então primeiro vamos discretizar ela utilizando a mediana como corte. Essa abordagem diminui a precisão por desconsiderar a grandeza dos valores
dados = dados |> dplyr::mutate(idade_fa = ifelse(idade<=median(idade), "Menos que 61", "Mais que 61"))Iniciando a análise vamos verificar graficamente as curvas de sobrevivência kaplan-meier para cada variável independente
ekm1<- survfit(data=dados,formula = Surv(tempo,status)~sexo,
conf.type = "log-log", conf.int=0.95)
ggsurvplot(ekm1,legend.title = "Sexo",
legend.labs = c("Feminino","Masculino"),
legend = c(0.2, 0.25),
conf.int=T)+
xlab("Tempo (dias)")+
ylab("Sobrevivencia estimada")
ekm2<- survfit(data=dados,formula = Surv(tempo,status)~idade_fa,
conf.type = "log-log", conf.int=0.95)
ggsurvplot(ekm2,legend.title = "Idade",
legend.labs = c("Mais que 61","Menos que 61"),
legend = c(0.2, 0.25),
conf.int=T)+
xlab("Tempo (dias)")+
ylab("Sobrevivencia estimada")
ekm3<- survfit(data=dados,formula = Surv(tempo,status)~desnut,
conf.type = "log-log", conf.int=0.95)
ggsurvplot(ekm3,legend.title = "Desnutrido",
legend.labs = c("Não","Sim"),
legend = c(0.2, 0.25),
conf.int=T)+
xlab("Tempo (dias)")+
ylab("Sobrevivencia estimada")
ekm4<- survfit(data=dados,formula = Surv(tempo,status)~comorbi,
conf.type = "log-log", conf.int=0.95)
ggsurvplot(ekm4,legend.title = "Comorbidade",
legend.labs = c("Não","Sim"),
legend = c(0.2, 0.25),
conf.int=T)+
xlab("Tempo (dias)")+
ylab("Sobrevivencia estimada")
ekm5<- survfit(data=dados,formula = Surv(tempo,status)~leucopenia,
conf.type = "log-log", conf.int=0.95)
ggsurvplot(ekm5,legend.title = "Leucopenia",
legend.labs = c("Não","Sim"),
legend = c(0.2, 0.25),
conf.int=T)+
xlab("Tempo (dias)")+
ylab("Sobrevivencia estimada")
ekm6<- survfit(data=dados,formula = Surv(tempo,status)~gptumor,
conf.type = "log-log", conf.int=0.95)
ggsurvplot(ekm6,legend.title = "Tipo de tumor",
legend.labs = c("Hemato","Loco","Mtx"),
legend = c(0.2, 0.25),
conf.int=T)+
xlab("Tempo (dias)")+
ylab("Sobrevivencia estimada")
Agora olharemos para os testes de comparação das curvas de sobrevivência de kaplan-meier para cada variável independente afim de verificar suas associações com a variável resposta utilizando as abordagens de logrank e peto.
Como temos uma variável qualitativa com mais de 2 curvas de sobrevivência vamos fazer um teste principal para verificar se há diferença entre os grupos
Teste de logrank
lr_gptumor<-survdiff(Surv(tempo,status)~gptumor, data=dados,rho=0)
lr_gptumor## Call:
## survdiff(formula = Surv(tempo, status) ~ gptumor, data = dados,
## rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## gptumor=Hemato 33 21 15.6 1.85 2.21
## gptumor=Loco 127 57 74.9 4.27 14.37
## gptumor=Mtx 40 30 17.5 8.93 10.98
##
## Chisq= 15.6 on 2 degrees of freedom, p= 4e-04Teste de Peto
pt_gptumor<-survdiff(Surv(tempo,status)~gptumor, data=dados,rho=1)
pt_gptumor## Call:
## survdiff(formula = Surv(tempo, status) ~ gptumor, data = dados,
## rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## gptumor=Hemato 33 16.1 11.8 1.58 2.42
## gptumor=Loco 127 41.3 54.6 3.27 13.55
## gptumor=Mtx 40 22.6 13.5 6.04 9.40
##
## Chisq= 14.3 on 2 degrees of freedom, p= 8e-04Usando os dois testes, de logrank e peto, para essa comparação e ao observar o p-valor nós rejeitamos a hipótese nula a nível 5% de significância nos dando evidências de que as curvas de sobrevivências dos 3 grupos são diferentes. Portanto abaixo faremos 3 subtestes para detectar essas diferenças (Gptumor loco x hemato,Gptumor hemato x mtx,Gptumor loco x mtx) além dos teste de comparação para o resto das variáveis.
lr_sexo<-survdiff(Surv(tempo,status)~sexo, data=dados,rho=0)
pt_sexo<-survdiff(Surv(tempo,status)~sexo, data=dados,rho=1)
lr_idade<-survdiff(Surv(tempo,status)~idade_fa, data=dados,rho=0)
pt_idade<-survdiff(Surv(tempo,status)~idade_fa, data=dados,rho=1)
gpt1 = dplyr::filter(dados,gptumor!="Mtx")
gpt2 = dplyr::filter(dados,gptumor!="Loco")
gpt3 = dplyr::filter(dados,gptumor!="Hemato")
lr_gptumor1<-survdiff(Surv(tempo,status)~gptumor, data=gpt1,rho=0)
pt_gptumor1<-survdiff(Surv(tempo,status)~gptumor, data=gpt1,rho=1)
lr_gptumor2<-survdiff(Surv(tempo,status)~gptumor, data=gpt2,rho=0)
pt_gptumor2<-survdiff(Surv(tempo,status)~gptumor, data=gpt2,rho=1)
lr_gptumor3<-survdiff(Surv(tempo,status)~gptumor, data=gpt3,rho=0)
pt_gptumor3<-survdiff(Surv(tempo,status)~gptumor, data=gpt3,rho=1)
lr_desnut<-survdiff(Surv(tempo,status)~desnut, data=dados,rho=0)
pt_desnut<-survdiff(Surv(tempo,status)~desnut, data=dados,rho=1)
lr_comorbi<-survdiff(Surv(tempo,status)~comorbi, data=dados,rho=0)
pt_comorbi<-survdiff(Surv(tempo,status)~comorbi, data=dados,rho=1)
lr_leucopenia<-survdiff(Surv(tempo,status)~leucopenia, data=dados,rho=0)
pt_leucopenia<-survdiff(Surv(tempo,status)~leucopenia, data=dados,rho=1)
comparar = rbind(c(lr_sexo$chisq,lr_sexo$pvalue,pt_sexo$chisq,pt_sexo$pvalue),
c(lr_idade$chisq,lr_idade$pvalue,pt_idade$chisq,pt_idade$pvalue),
c(lr_gptumor1$chisq,lr_gptumor1$pvalue,pt_gptumor1$chisq,pt_gptumor1$pvalue), c(lr_gptumor2$chisq,lr_gptumor2$pvalue,pt_gptumor2$chisq,pt_gptumor2$pvalue), c(lr_gptumor3$chisq,lr_gptumor3$pvalue,pt_gptumor3$chisq,pt_gptumor3$pvalue),
c(lr_comorbi$chisq,lr_comorbi$pvalue,pt_comorbi$chisq,pt_comorbi$pvalue),
c(lr_leucopenia$chisq,lr_leucopenia$pvalue,pt_leucopenia$chisq,pt_leucopenia$pvalue),
c(lr_desnut$chisq,lr_desnut$pvalue,pt_desnut$chisq,pt_desnut$pvalue))
rownames(comparar) = c('Sexo','Idade','Gptumor loco x hemato','Gptumor hemato x mtx','Gptumor loco x mtx','Comorbi','Leucopenia','Desnut')
colnames(comparar) = c('logrank-est','logrank-pvalor','peto-est','peto-pvalor')
comparar## logrank-est logrank-pvalor peto-est peto-pvalor
## Sexo 1.6650144 0.1969276469 2.4218076 0.119656786
## Idade 5.2278616 0.0222278441 6.0524047 0.013887371
## Gptumor loco x hemato 5.1431232 0.0233386252 5.3367700 0.020880126
## Gptumor hemato x mtx 0.6913963 0.4056898656 0.3470285 0.555800084
## Gptumor loco x mtx 14.3971719 0.0001480245 12.9304158 0.000323286
## Comorbi 2.7750710 0.0957424085 2.4285432 0.119143575
## Leucopenia 8.8275851 0.0029671077 8.9845132 0.002722773
## Desnut 1.0091821 0.3150988611 2.0421060 0.152997696Olhando o p-valor das variáveis sexo, idade, comorbi e leucopenia e considerando o teste de logrank a nível 5% de significância nós rejeitamos a hipótese nula para Leucopenia ou seja temos evidências de que há diferença nessa curva de sobrevivência. Não foi o caso para sexo, comorbi nem idade Agora olhando o teste de peto tivemos resultados semelhantes porém desta vez rejeitamos a hipótese nula para Leucopenia e idade, diferentemente do teste de logrank. Novamente não rejeitamos a hipótese nula para sexo e comorbi
Analisando agora os 3 subtestes entre os grupos, nós iremos usar o nível de significância corrigido pelo método de bonferroni para não inflar o nível de significância na realização desses novos testes de hipóteses. Então teremos que \({\alpha}*=\frac{\alpha}{k} = \frac{0.05}{3} = 0.017\) e com isso rejeitamos hipoétese nula apenas no subteste entre os Grupos ‘Loxo x mtx’, logo temos evidências de que existe diferenças entre essas curvas de sobrevivência.Para os Grupos ‘loco x hemato’ e ‘hemato x mtx’ não rejeitamos a hipótese nula, logo temos evidências de que não há diferença nas curvas.
Modelo de regressão paramétrico
Reutilizando resultado anterior e de acordo com a literatura as variáveis com p-valor < 0.25 se mantém para a abordagem dos modelos de regressão. ## Devo olhar o p-valor do grupo tumor por si só ou olhar o p-valor da cada grupo?
Agora iremos verificar se cada variável é significativa
mod1 = flexsurvreg(data=dados,Surv(tempo,status)~sexo,dist="gengamma.orig")## Warning in .hess_to_cov(opt$hessian, hess.control$tol.solve,
## hess.control$tol.evalues): Hessian not positive definite: smallest eigenvalue
## is -2.8e+02 (threshold: -1.0e-05). This might indicate that the optimization
## did not converge to the maximum likelihood, so that the results are invalid.
## Continuing with the nearest positive definite approximation of the covariance
## matrix.mod2 = flexsurvreg(data=dados,Surv(tempo,status)~idade_fa,dist="gengamma.orig")
mod3 = flexsurvreg(data=dados,Surv(tempo,status)~gptumor,dist="gengamma.orig")## Warning in .hess_to_cov(opt$hessian, hess.control$tol.solve,
## hess.control$tol.evalues): Hessian not positive definite: smallest eigenvalue
## is -7.8e+01 (threshold: -1.0e-05). This might indicate that the optimization
## did not converge to the maximum likelihood, so that the results are invalid.
## Continuing with the nearest positive definite approximation of the covariance
## matrix.mod4 = flexsurvreg(data=dados,Surv(tempo,status)~comorbi,dist="gengamma.orig")
mod5 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia,dist="gengamma.orig")
mod6 = flexsurvreg(data=dados,Surv(tempo,status)~desnut,dist="gengamma.orig")## Warning in .hess_to_cov(opt$hessian, hess.control$tol.solve,
## hess.control$tol.evalues): Hessian not positive definite: smallest eigenvalue
## is -1.6e+02 (threshold: -1.0e-05). This might indicate that the optimization
## did not converge to the maximum likelihood, so that the results are invalid.
## Continuing with the nearest positive definite approximation of the covariance
## matrix.tidy(mod1)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 7.98e- 2 2.52e- 3 NA NA
## 2 scale 2.79e-15 6.27e-18 NA NA
## 3 k 2.21e+ 1 1.95e+ 0 NA NA
## 4 sexoM -6.71e- 1 4.32e- 1 -1.55 0.121tidy(mod2)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 7.68e- 2 1.52e- 2 NA NA
## 2 scale 7.69e-17 1.01e-15 NA NA
## 3 k 2.42e+ 1 9.50e+ 0 NA NA
## 4 idade_faMenos que 61 1.02e+ 0 4.22e- 1 2.42 0.0154tidy(mod3)## # A tibble: 5 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 8.22e- 2 2.55e- 3 NA NA
## 2 scale 2.38e-15 9.76e-18 NA NA
## 3 k 2.24e+ 1 1.99e+ 0 NA NA
## 4 gptumorLoco 1.29e+ 0 5.54e- 1 2.32 0.0202
## 5 gptumorMtx -5.11e- 1 6.44e- 1 -0.794 0.427tidy(mod4)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 7.66e- 2 9.48e- 3 NA NA
## 2 scale 1.51e-16 1.23e-15 NA NA
## 3 k 2.39e+ 1 6.16e+ 0 NA NA
## 4 comorbisim -1.27e+ 0 7.88e- 1 -1.61 0.107tidy(mod5)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 7.72e- 2 7.42e- 3 NA NA
## 2 scale 1.99e-16 1.24e-15 NA NA
## 3 k 2.42e+ 1 5.10e+ 0 NA NA
## 4 leucopeniasim -1.76e+ 0 6.70e- 1 -2.63 0.00863tidy(mod6)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 7.47e- 2 2.28e- 3 NA NA
## 2 scale 3.26e-17 7.26e-20 NA NA
## 3 k 2.47e+ 1 2.22e+ 0 NA NA
## 4 desnutsim -9.84e- 1 7.86e- 1 -1.25 0.211Olhando para o p-valor de cada termo e a nível 10% de significância rejeitamos a hipótese nula de que as variáveis idade , leucopenia e (grupo tal) tumor são iguais a zero, logo são significativas para explicar o modelo.
mod_1 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor,dist="gengamma.orig")
mod_2 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade,dist="gengamma.orig")
mod_3 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+gptumor,dist="gengamma.orig")
mod_4 = flexsurvreg(data=dados,Surv(tempo,status)~gptumor+idade,dist="gengamma.orig")
## AIC BIC
aic_bic_2 = rbind(c(mod_1$AIC,mod_1$BIC),
c(mod_2$AIC,mod_2$BIC),
c(mod_3$AIC,mod_3$BIC),
c(mod_4$AIC,mod_4$BIC))
rownames(aic_bic_2) = c("Completo","Sem gptumor","Sem idade","Sem leucopenia")
colnames(aic_bic_2) = c("AIC","BIC")
aic_bic_2## AIC BIC
## Completo 1198.508 1221.596
## Sem gptumor 1208.127 1224.618
## Sem idade 1206.196 1225.986
## Sem leucopenia 1201.553 1221.343Observando os AIC e BIC dos modelos com as variáveis significativas observamos que nesse momento o que melhor se ajusta é o completo (com idade, leucopenia e gptumor).
mod_5 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor,dist="gengamma.orig")
mod_6 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+sexo,dist="gengamma.orig")
mod_7 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+comorbi,dist="gengamma.orig")
mod_8 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut,dist="gengamma.orig")
## AIC BIC
aic_bic_3 = rbind(c(mod_5$AIC,mod_5$BIC),
c(mod_6$AIC,mod_6$BIC),
c(mod_7$AIC,mod_7$BIC),
c(mod_8$AIC,mod_8$BIC))
rownames(aic_bic_3) = c("Completo","Com sexo","Com comorbi","Com desnut")
colnames(aic_bic_3) = c("AIC","BIC")
aic_bic_3## AIC BIC
## Completo 1198.508 1221.596
## Com sexo 1198.632 1225.019
## Com comorbi 1199.250 1225.636
## Com desnut 1198.469 1224.856Ao retornar no modelo as variáveis anteriormente consideradas não significativas e olhando para os AIC e BIC de cada modelo, nós chegamos a conclusão de que o modelo com a variável desnut melhor se ajusta agora no modelo completo. Voltaremos com ela. Agora nos iremos verificar se alguma das variáveis do nosso modelo completo atual precisa sair.
mod5 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut,dist="gengamma.orig")
mod6 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+desnut,dist="gengamma.orig")## Warning in .hess_to_cov(opt$hessian, hess.control$tol.solve,
## hess.control$tol.evalues): Hessian not positive definite: smallest eigenvalue
## is -2.2e+02 (threshold: -1.0e-05). This might indicate that the optimization
## did not converge to the maximum likelihood, so that the results are invalid.
## Continuing with the nearest positive definite approximation of the covariance
## matrix.mod7 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+gptumor+desnut,dist="gengamma.orig")
mod8 = flexsurvreg(data=dados,Surv(tempo,status)~idade+gptumor+desnut,dist="gengamma.orig")
## AIC BIC
aic_bic_4 = rbind(c(mod5$AIC,mod5$BIC),
c(mod6$AIC,mod6$BIC),
c(mod7$AIC,mod7$BIC),
c(mod8$AIC,mod8$BIC))
rownames(aic_bic_4) = c("Completo","Sem gptumor","Sem idade","Sem leucopenia")
colnames(aic_bic_4) = c("AIC","BIC")
aic_bic_4## AIC BIC
## Completo 1198.469 1224.856
## Sem gptumor 1208.384 1228.174
## Sem idade 1206.076 1229.164
## Sem leucopenia 1201.672 1224.760O AIC e BIC do modelo completo continua sendo o melhor, logo não iremos tirar nenhuma variável do nosso modelo atual (leucopenia, idade gptumor e desnut). Com esse modelo então, iremos testar suas interações
mod9 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut,dist="gengamma.orig")
mod10 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut+leucopenia*idade,dist="gengamma.orig")## Warning in .hess_to_cov(opt$hessian, hess.control$tol.solve,
## hess.control$tol.evalues): Hessian not positive definite: smallest eigenvalue
## is -1.1e+02 (threshold: -1.0e-05). This might indicate that the optimization
## did not converge to the maximum likelihood, so that the results are invalid.
## Continuing with the nearest positive definite approximation of the covariance
## matrix.mod11 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut+leucopenia*gptumor,dist="gengamma.orig")## Warning in .hess_to_cov(opt$hessian, hess.control$tol.solve,
## hess.control$tol.evalues): Hessian not positive definite: smallest eigenvalue
## is -1.0e+03 (threshold: -1.0e-05). This might indicate that the optimization
## did not converge to the maximum likelihood, so that the results are invalid.
## Continuing with the nearest positive definite approximation of the covariance
## matrix.mod12 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut+leucopenia*desnut,dist="gengamma.orig")## Warning in .hess_to_cov(opt$hessian, hess.control$tol.solve,
## hess.control$tol.evalues): Hessian not positive definite: smallest eigenvalue
## is -5.0e+02 (threshold: -1.0e-05). This might indicate that the optimization
## did not converge to the maximum likelihood, so that the results are invalid.
## Continuing with the nearest positive definite approximation of the covariance
## matrix.mod13 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut+idade*gptumor,dist="gengamma.orig")## Warning in .hess_to_cov(opt$hessian, hess.control$tol.solve,
## hess.control$tol.evalues): Hessian not positive definite: smallest eigenvalue
## is -1.2e+02 (threshold: -1.0e-05). This might indicate that the optimization
## did not converge to the maximum likelihood, so that the results are invalid.
## Continuing with the nearest positive definite approximation of the covariance
## matrix.mod14 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut+idade*desnut,dist="gengamma.orig")## Warning in .hess_to_cov(opt$hessian, hess.control$tol.solve,
## hess.control$tol.evalues): Hessian not positive definite: smallest eigenvalue
## is -2.8e+02 (threshold: -1.0e-05). This might indicate that the optimization
## did not converge to the maximum likelihood, so that the results are invalid.
## Continuing with the nearest positive definite approximation of the covariance
## matrix.mod15 = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut+gptumor*desnut,dist="gengamma.orig")## Warning in .hess_to_cov(opt$hessian, hess.control$tol.solve,
## hess.control$tol.evalues): Hessian not positive definite: smallest eigenvalue
## is -1.5e+02 (threshold: -1.0e-05). This might indicate that the optimization
## did not converge to the maximum likelihood, so that the results are invalid.
## Continuing with the nearest positive definite approximation of the covariance
## matrix.## AIC BIC
aic_bic_5 = rbind(c(mod9$AIC,mod9$BIC),
c(mod10$AIC,mod10$BIC),
c(mod11$AIC,mod11$BIC),
c(mod12$AIC,mod12$BIC),
c(mod13$AIC,mod13$BIC),
c(mod14$AIC,mod14$BIC),
c(mod15$AIC,mod15$BIC))
rownames(aic_bic_5) = c("Completo","mod10","mod11","mod12","mod13","mod14","mod15")
colnames(aic_bic_5) = c("AIC","BIC")
aic_bic_5## AIC BIC
## Completo 1198.469 1224.856
## mod10 1196.397 1226.082
## mod11 1198.860 1231.843
## mod12 1199.092 1228.777
## mod13 1195.145 1228.128
## mod14 1199.855 1229.540
## mod15 1201.197 1234.180Testando as initerações pudemos observar que o modelo 13 com a interação de idade e gptumor teve um melhor ajuste. Logo iremos adiciona-lá no nosso modelo. Segue abaixo o ajuste
tidy(mod13)## # A tibble: 10 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 8.65e- 2 3.18e- 3 NA NA
## 2 scale 2.36e-14 1.05e-15 NA NA
## 3 k 2.30e+ 1 2.22e+ 0 NA NA
## 4 leucopeniasim -1.44e+ 0 6.77e- 1 -2.13 0.0331
## 5 idade -4.62e- 3 2.27e- 2 -0.204 0.838
## 6 gptumorLoco 5.14e+ 0 1.75e+ 0 2.93 0.00340
## 7 gptumorMtx -5.95e- 1 1.97e+ 0 -0.303 0.762
## 8 desnutsim -1.09e+ 0 7.06e- 1 -1.54 0.124
## 9 idade:gptumorLoco -6.87e- 2 2.95e- 2 -2.33 0.0199
## 10 idade:gptumorMtx -3.19e- 3 3.45e- 2 -0.0927 0.926Vamos agora verificar qual a melhor distribuição para o nosso modelo
mod_gam = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut+idade*gptumor,dist="gengamma.orig")## Warning in .hess_to_cov(opt$hessian, hess.control$tol.solve,
## hess.control$tol.evalues): Hessian not positive definite: smallest eigenvalue
## is -1.2e+02 (threshold: -1.0e-05). This might indicate that the optimization
## did not converge to the maximum likelihood, so that the results are invalid.
## Continuing with the nearest positive definite approximation of the covariance
## matrix.mod_exp = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut+idade*gptumor,dist="exponential")
mod_weibull = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut+idade*gptumor,dist="weibull")
mod_lnorm = flexsurvreg(data=dados,Surv(tempo,status)~leucopenia+idade+gptumor+desnut+idade*gptumor,dist="lnorm")
aic_bic_6 = rbind(c(mod_gam$AIC,mod_gam$BIC),
c(mod_exp$AIC,mod_exp$BIC),
c(mod_weibull$AIC,mod_weibull$BIC),
c(mod_lnorm$AIC,mod_lnorm$BIC))
rownames(aic_bic_6) = c("Completo","Exponencial","Weibull","Log-normal")
colnames(aic_bic_6) = c("AIC","BIC")
aic_bic_6## AIC BIC
## Completo 1195.145 1228.128
## Exponencial 1304.832 1331.218
## Weibull 1211.666 1241.351
## Log-normal 1189.059 1218.743Ao testar outros modelos vemos que o modelo log-normal é o que tem o melhor AIC-BIC, portanto será nosso modelo final. Iremos realizar então a análise de resíduos do nosso modelo final (desnut não significativo porém AIC menor, mantemos?)
CoxSnell <-coxsnell_flexsurvreg(mod_lnorm)$est
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")
De acordo com os gráficos, o modelo esta bem ajustado (ta?)
mod_lnorm## Call:
## flexsurvreg(formula = Surv(tempo, status) ~ leucopenia + idade +
## gptumor + desnut + idade * gptumor, data = dados, dist = "lnorm")
##
## Estimates:
## data mean est L95% U95% se
## meanlog NA 4.90e+00 2.35e+00 7.44e+00 1.30e+00
## sdlog NA 2.49e+00 2.15e+00 2.89e+00 1.89e-01
## leucopeniasim 9.50e-02 -1.41e+00 -2.77e+00 -4.99e-02 6.94e-01
## idade 5.83e+01 -7.74e-03 -5.28e-02 3.73e-02 2.30e-02
## gptumorLoco 6.35e-01 4.97e+00 1.58e+00 8.35e+00 1.73e+00
## gptumorMtx 2.00e-01 -1.13e+00 -4.98e+00 2.72e+00 1.96e+00
## desnutsim 7.50e-02 -1.17e+00 -2.57e+00 2.24e-01 7.12e-01
## idade:gptumorLoco 3.85e+01 -6.62e-02 -1.23e-01 -8.95e-03 2.92e-02
## idade:gptumorMtx 1.14e+01 6.06e-03 -6.13e-02 7.34e-02 3.44e-02
## exp(est) L95% U95%
## meanlog NA NA NA
## sdlog NA NA NA
## leucopeniasim 2.44e-01 6.27e-02 9.51e-01
## idade 9.92e-01 9.49e-01 1.04e+00
## gptumorLoco 1.43e+02 4.85e+00 4.24e+03
## gptumorMtx 3.24e-01 6.89e-03 1.52e+01
## desnutsim 3.10e-01 7.68e-02 1.25e+00
## idade:gptumorLoco 9.36e-01 8.84e-01 9.91e-01
## idade:gptumorMtx 1.01e+00 9.41e-01 1.08e+00
##
## N = 200, Events: 108, Censored: 92
## Total time at risk: 19800
## Log-likelihood = -585.5293, df = 9
## AIC = 1189.059exp(-1.40977)## [1] 0.2441994exp(-0.00774)## [1] 0.9922899exp(4.96546)## [1] 143.3745exp(-1.12781)## [1] 0.3237415exp(-1.17147)## [1] 0.309911exp(-0.06620)## [1] 0.9359437exp(0.00606)## [1] 1.006078Sobre a interpretação isoladamente conseguimos interpretar apenas a variável idade. Para pacientes i e j onde ambos