ANÁLISE DE SOBREVIVÊNCIA E CONFIABILIDADE
AULA 6: SELEÇÃO DE VARIÁVEIS
\[\\[0.05in]\]1 PACOTES NECESSÁRIOS
Para esta sexta 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 BASE DE DADOS: ALEITAMENTO MATERNO
As Organizações Internacionais de Saúde recomendam o leite materno como a única fonte de alimentação para crianças até 6 meses de idade. Identificar fatores determinantes do aleitamento materno em diferentes populações é, portanto, fundamental para alcançar tal recomendação. Muitos trabalhos investigam os determinantes do aleitamento exclusivamente materno em populações urbanas de baixa renda. Os resultados desses trabalhos mostram que a condição sócio-econômica e o peso ao nascimento da criança estão associados com o aleitamento exclusivamente materno.
Além disso, as mulheres que têm acesso a maternidades que promovem programas de aleitamento obtêm melhores resultados. Nessa mesma linha de pesquisa, um outro estudo foi realizado pelos professores Eugênio Goulart e Cláudia Lindgren, do Departamento de Pediatria da UFMG, no Centro de Saúde São Marcos, localizado em Belo Horizonte. Este centro é um ambulatório municipal que atende essencialmente a população de baixa renda. O estudo teve como objetivos principais conhecer a prática do aleitamento materno de mães que utilizam este centro, assim como os possíveis fatores de risco ou de proteção para o desmame precoce.
Um inquérito epidemiológico composto por questões demográficas e comportamentais foi aplicado a 150 mães de crianças menores de 2 anos de idade. A variável resposta de interesse foi o tempo máximo de aleitamento materno, em meses, ou seja, o tempo contado a partir do nascimento até o desmame completo da criança. No estudo foram registradas 11 covariáveis e a variável resposta. Algumas crianças não foram acompanhadas até o desmame e, portanto, há a presença de censuras.
A variável resposta é, portanto, o par \((t_i,\delta_i)\), com \(i=1,2,...,150\), de modo que:
- \(t_i\) representa o tempo observado, em meses.
- \(\delta_i\) é a indicadora de ocorrência do desmame.
As covariáveis explicativas são:
- \(V_1\): Experiência anterior amamentação = sim e não.
- \(V_2\): Número de filhos vivos = \(\leq 2\) e >2
- \(V_3\): Conceito materno sobre o tempo ideal de amamentação = > 6 meses e \(\leq\) 6 meses
- \(V_4\): Dificuldades para amamentar nos primeiros dias pós-parto = não e sim.
- \(V_5\): Tipo de serviço em que realizou o pré-natal = público e privado/convênios.
- \(V_6\): Recebeu exclusivamente leite materno na maternidade = sim e não.
- \(V_7\): A criança teve contato com o pai = sim e não.
- \(V_8\): Renda per capita (em SM/mes) = \(\geq\) 1 SM e < 1 SM.
- \(V_9\): Peso ao nascimento = \(\geq\) 2,5 kg e < 2,5 kg.
- \(V_{10}\): Tempo de separação mãe-filho pós-parto = \(\leq\) 6 horas e >6 horas.
- \(V_{11}\): Permanência no berçário = não e sim.
A base de dados está armazenada no arquivo
dados.csv
Os comandos a seguir fazem a leitura dessa base de dados:
dados<-read.csv("dados.csv")3 FILTRO INICIAL
Como em qualquer análise estatística, um passo importante é conduzir uma boa análise descritiva dos dados. Em análise de sobrevivência e confiabilidade não é diferente. Medidas resumo podem ser utilizadas para todas as variáveis, fazendo as devidas adaptações quando necessário.
Trazendo para o contexto de análise de sobrevivência e confiabilidade, podemos usar as técnicas não paramétricas para investigar isoladamente a associação de cada variável qualitativa com a variável resposta. Para cada variável qualitativa, podemos construir as estimativas para as funções de sobrevivência (confiabilidade) segundo Kaplan-Meier. Na sequência, utilizar os testes para comparar as funções de sobrevivência (confiabilidade) nas categorias da covariável qualitativa. Esse passo é importante para verificar como a covariável qualitativa afeta a variável resposta.
Existe uma proposta na literatura de levar para a modelagem de regressão apenas as covariáveis qualitativas que tiverem um valor-p inferior a 0,25 em pelo menos um dos testes de comparação das curvas de sobrevivência (confiabilidade). Este passo acaba se tornando um filtro inicial para a seleção de variáveis qualitativas que poderão entrar no modelo de regressão final ajustado.
No caso de uma variável quantitativa, pode ser possível fazer um gráfico de dispersão entre ela e a variável resposta e demarcar as observações que são censuras ou não. Esse passo pode ajudar a investigar, de forma subjetiva, como a variável quantitativa pode afetar a variável resposta. Além disso, pode ser possível transformá-la em categórica e utilizar os testes de comparação de curvas de sobrevivência (confiabilidade). Essa abordagem pode não ser tão precisa, uma vez que é removida a grandeza dos valores.
Para a nossa base de dados de aleitamento materno, temos apenas variáveis qualitativas. Dessa forma, os comandos a seguir fornecem os testes de comparação das curvas de sobrevivência para todas elas, considerando as abordagens Logrank e Peto:
dados<-read.csv("dados.csv")
V1_LR<-survdiff(Surv(tempo,cens)~V1, data=dados,rho=0)
V1_PT<-survdiff(Surv(tempo,cens)~V1, data=dados,rho=1)
V2_LR<-survdiff(Surv(tempo,cens)~V2, data=dados,rho=0)
V2_PT<-survdiff(Surv(tempo,cens)~V2, data=dados,rho=1)
V3_LR<-survdiff(Surv(tempo,cens)~V3, data=dados,rho=0)
V3_PT<-survdiff(Surv(tempo,cens)~V3, data=dados,rho=1)
V4_LR<-survdiff(Surv(tempo,cens)~V4, data=dados,rho=0)
V4_PT<-survdiff(Surv(tempo,cens)~V4, data=dados,rho=1)
V5_LR<-survdiff(Surv(tempo,cens)~V5, data=dados,rho=0)
V5_PT<-survdiff(Surv(tempo,cens)~V5, data=dados,rho=1)
V6_LR<-survdiff(Surv(tempo,cens)~V6, data=dados,rho=0)
V6_PT<-survdiff(Surv(tempo,cens)~V6, data=dados,rho=1)
V7_LR<-survdiff(Surv(tempo,cens)~V7, data=dados,rho=0)
V7_PT<-survdiff(Surv(tempo,cens)~V7, data=dados,rho=1)
V8_LR<-survdiff(Surv(tempo,cens)~V8, data=dados,rho=0)
V8_PT<-survdiff(Surv(tempo,cens)~V8, data=dados,rho=1)
V9_LR<-survdiff(Surv(tempo,cens)~V9, data=dados,rho=0)
V9_PT<-survdiff(Surv(tempo,cens)~V9, data=dados,rho=1)
V10_LR<-survdiff(Surv(tempo,cens)~V10, data=dados,rho=0)
V10_PT<-survdiff(Surv(tempo,cens)~V10, data=dados,rho=1)
V11_LR<-survdiff(Surv(tempo,cens)~V11, data=dados,rho=0)
V11_PT<-survdiff(Surv(tempo,cens)~V11, data=dados,rho=1)
COMPARACAO_SOBREVIVENCIA <- rbind(c(V1_LR$chisq,V1_LR$pvalue,V1_PT$chisq,V1_PT$pvalue),
c(V2_LR$chisq,V2_LR$pvalue,V2_PT$chisq,V2_PT$pvalue),
c(V3_LR$chisq,V3_LR$pvalue,V3_PT$chisq,V3_PT$pvalue),
c(V4_LR$chisq,V4_LR$pvalue,V4_PT$chisq,V4_PT$pvalue),
c(V5_LR$chisq,V5_LR$pvalue,V5_PT$chisq,V5_PT$pvalue),
c(V6_LR$chisq,V6_LR$pvalue,V6_PT$chisq,V6_PT$pvalue),
c(V7_LR$chisq,V7_LR$pvalue,V7_PT$chisq,V7_PT$pvalue),
c(V8_LR$chisq,V8_LR$pvalue,V8_PT$chisq,V8_PT$pvalue),
c(V9_LR$chisq,V9_LR$pvalue,V9_PT$chisq,V9_PT$pvalue),
c(V10_LR$chisq,V10_LR$pvalue,V10_PT$chisq,V10_PT$pvalue),
c(V11_LR$chisq,V11_LR$pvalue,V11_PT$chisq,V11_PT$pvalue))
rownames(COMPARACAO_SOBREVIVENCIA)<-c("V1","V2","V3","V4","V5","V6",
"V7","V8","V9","V10","V11")
colnames(COMPARACAO_SOBREVIVENCIA)<-c("Logrank_Est","Logrank_pvalor","Peto_Est","Peto_pvalor")
COMPARACAO_SOBREVIVENCIA## Logrank_Est Logrank_pvalor Peto_Est Peto_pvalor
## V1 3.948917 0.0469015372 5.665980 0.0172970509
## V2 2.598897 0.1069381287 2.394640 0.1217517744
## V3 6.148444 0.0131527781 7.934190 0.0048509325
## V4 12.259259 0.0004629552 14.893645 0.0001137455
## V5 1.375355 0.2408939587 1.285637 0.2568536295
## V6 7.465871 0.0062879584 7.019143 0.0080642837
## V7 1.840854 0.1748508878 1.450923 0.2283799192
## V8 2.114927 0.1458689925 2.219480 0.1362789295
## V9 1.873899 0.1710292070 2.380470 0.1228606074
## V10 2.604718 0.1065461194 1.582213 0.2084422050
## V11 2.932307 0.0868233511 1.671408 0.1960700654Assim, segundo o critério da literatura, todas as covariáveis qualitativas obtiveram um p-valor < 0,25 em pelo menos um dos testes de comparação. Dessa forma, todas elas seguem para a abordagem dos modelos de regressão.
Na aula passada, vimos o ajuste de três modelos de regressão paramétricos:
- Modelo de regressão Exponencial.
- Modelo de regressão Weibull.
- Modelo de regressão LogNormal.
Temos que fazer uma seleção de covariáveis e também de um modelo probabilístico para a variável resposta. A estratégia é trabalhar com um modelo que generaliza as três distribuições de probabilidade e, posteriormente, verificar se o ajuste de alguma delas é preferível.
Lembre-se que a função densidade de probabilidade de uma variável aleatória T com distribuição gama generalizada é:
\[f(t)=\frac{\gamma\alpha^{\gamma\beta}}{\Gamma(\beta)}t^{\gamma\beta-1}exp\{-(\alpha t)^\gamma\}\]
Para esta distribuição, tem-se um parâmetro de taxa \(\alpha\) e dois de forma \(\gamma\) e \(\beta\), o que a torna bastante flexível.
- para \(\gamma=\beta=1\), tem-se que \(T\sim exponencial (\alpha)\)
- para \(\beta=1\), tem-se que \(T\sim Weibull(\gamma,\alpha)\) (reparametrizada)
- para \(\gamma=1\), tem-se que \(T\sim Gama(\alpha,\beta)\)
Além disso, pode-se mostrar, que a distribuição log-normal aparece como um caso limite da distribuição gama generalizada quando \(\beta \rightarrow \infty\)}
Assim, a estratégia de seleção de covariáveis é feita no modelo de regressão gama generalizado. Esse modelo é mais robusto e, como dito, generaliza as distribuições Exponencial, Weibull e Lognormal. Com o modelo final ajustado, testamos se o modelo gama generalizado é o mais adequado ou alguma de suas versões particulares é preferível. Temos então uma estrutura de modelos encaixados. No último tópico da aula de hoje, estudaremos como fazer essa escolha do melhor modelo probabilístico.
4 ESTRATÉGIA PARA A SELEÇÃO DE VARIÁVEIS
As covariáveis potencialmente importantes para descrever o comportamento da resposta são eleitas para a abordagem dos modelos de regressão. Dependendo do número de covariáveis, há um número gigante de possíveis modelos que podem ser formados. É certamente impraticável ajustar todos estes possíveis modelos a fim de ser selecionado o que melhor explique a variável resposta. Nessas situações, rotinas automáticas para seleção de covariáveis podem ser utilizadas,tais como os métodos forward, backward ou stepwise.
Estes métodos estão implementados e, portanto, disponíveis em pacotes estatísticos. Entretanto, tais rotinas possuem algumas desvantagens. Tipicamente, elas tendem a identificar um particular conjunto de covariáveis, em vez de possíveis conjuntos igualmente bons para explicar a resposta. Esse fato impossibilita que dois ou mais conjuntos de covariáveis igualmente bons sejam apresentados para o pesquisador, para a escolha do mais relevante em sua área de aplicação. Isto significa que esses métodos são automáticos e fazem com que o pacote estatístico escolha o modelo.
Na realidade, o que se defende aqui é que o estatístico e o pesquisador tenham uma postura pro-ativa neste processo. Isto implica, por exemplo, que covariáveis importantes em termos clínicos ou de engenharia devem ser consideradas em cada passo de inclusão ou exclusão no processo de seleção de covariáveis. Frente a estas limitações das rotinas automáticas, optou-se por utilizar métodos que envolvem a interferência mais de perto do analista. Neste estudo optou-se por utilizar uma estratégia de seleção de modelos derivada da proposta de Collett (2003a).
4.1 PASSO 1
Nesse passo, vamos ajustar todos os modelos contendo uma única covariável. Manteremos apenas as covariáveis que forem significativas ao nível de 0,10. É aconselhável utilizar o teste de Wald nesse passo
Os comandos a seguir ajusta todos esses modelos:
mod1<-flexsurvreg(data=dados,Surv(tempo,cens)~V1,dist="gengamma.orig")
mod2<-flexsurvreg(data=dados,Surv(tempo,cens)~V2,dist="gengamma.orig")
mod3<-flexsurvreg(data=dados,Surv(tempo,cens)~V3,dist="gengamma.orig")
mod4<-flexsurvreg(data=dados,Surv(tempo,cens)~V4,dist="gengamma.orig")
mod5<-flexsurvreg(data=dados,Surv(tempo,cens)~V5,dist="gengamma.orig")
mod6<-flexsurvreg(data=dados,Surv(tempo,cens)~V6,dist="gengamma.orig")
mod7<-flexsurvreg(data=dados,Surv(tempo,cens)~V7,dist="gengamma.orig")
mod8<-flexsurvreg(data=dados,Surv(tempo,cens)~V8,dist="gengamma.orig")
mod9<-flexsurvreg(data=dados,Surv(tempo,cens)~V9,dist="gengamma.orig")
mod10<-flexsurvreg(data=dados,Surv(tempo,cens)~V10,dist="gengamma.orig")
mod11<-flexsurvreg(data=dados,Surv(tempo,cens)~V11,dist="gengamma.orig")
tidy(mod1)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 0.227 0.314 NA NA
## 2 scale 0.000377 0.00927 NA NA
## 3 k 9.29 22.9 NA NA
## 4 V1sim 0.701 0.317 2.21 0.0272tidy(mod2)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 0.339 0.181 NA NA
## 2 scale 0.205 1.04 NA NA
## 3 k 4.51 4.05 NA NA
## 4 V2menor ou igual a 2 -0.609 0.382 -1.59 0.111tidy(mod3)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 0.244 0.254 NA NA
## 2 scale 0.00250 0.0417 NA NA
## 3 k 8.24 15.3 NA NA
## 4 V3menor ou igual a 6 meses -0.725 0.307 -2.36 0.0182tidy(mod4)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 0.216 0.0763 NA NA
## 2 scale 0.000266 0.00183 NA NA
## 3 k 10.8 7.00 NA NA
## 4 V4sim -1.07 0.298 -3.58 0.000342tidy(mod5)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 0.257 0.183 NA NA
## 2 scale 0.00323 0.0338 NA NA
## 3 k 7.21 9.04 NA NA
## 4 V5público 0.500 0.371 1.35 0.178tidy(mod6)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 0.361 0.330 NA NA
## 2 scale 0.117 0.946 NA NA
## 3 k 4.25 6.52 NA NA
## 4 V6sim 0.791 0.297 2.67 0.00770tidy(mod7)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 0.442 0.347 NA NA
## 2 scale 2.68 13.5 NA NA
## 3 k 2.95 3.70 NA NA
## 4 V7sim -1.07 0.778 -1.38 0.169tidy(mod8)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 0.302 0.314 NA NA
## 2 scale 0.0454 0.543 NA NA
## 3 k 5.59 9.97 NA NA
## 4 V8menor que 1 SM -0.515 0.302 -1.71 0.0882tidy(mod9)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 0.299 0.305 NA NA
## 2 scale 0.0285 0.340 NA NA
## 3 k 5.67 9.96 NA NA
## 4 V9menor que 2.5 kg 0.809 0.501 1.61 0.107tidy(mod10)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 0.475 0.364 NA NA
## 2 scale 1.22 5.15 NA NA
## 3 k 2.62 3.16 NA NA
## 4 V10menor ou igual a 6 horas 0.424 0.298 1.42 0.155tidy(mod11)## # A tibble: 4 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 shape 0.435 0.348 NA NA
## 2 scale 1.10 5.63 NA NA
## 3 k 3.00 3.86 NA NA
## 4 V11sim -0.432 0.295 -1.46 0.143De acordo com o critério do passo 1, vamos manter as variáveis V1, V3, V4, V6 e V8.
5 PASSO 2
As covariáveis significativas no passo 1 são, então, ajustadas conjuntamente. Na presença de certas covariáveis, outras podem deixar de ser significativas. Conseqüentemente, ajustam-se modelos reduzidos, excluindo-se uma única covariável de cada vez. Verificam-se as covariáveis que, quando retiradas, provocam um aumento em pelo menos um dos critérios AIC e BIC. Somente aquelas que provocarem esse aumento permanecerão no modelo.
mod1<-flexsurvreg(data=dados,Surv(tempo,cens)~V1+V3+V4+V6+V8,dist="gengamma.orig")
mod2<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8,dist="gengamma.orig")
mod3<-flexsurvreg(data=dados,Surv(tempo,cens)~V1+V4+V6+V8,dist="gengamma.orig")
mod4<-flexsurvreg(data=dados,Surv(tempo,cens)~V1+V3+V6+V8,dist="gengamma.orig")
mod5<-flexsurvreg(data=dados,Surv(tempo,cens)~V1+V3+V4+V8,dist="gengamma.orig")
mod6<-flexsurvreg(data=dados,Surv(tempo,cens)~V1+V3+V4+V6,dist="gengamma.orig")
Comparacao<-rbind(c(mod1$AIC,mod1$BIC),
c(mod2$AIC,mod2$BIC),
c(mod3$AIC,mod3$BIC),
c(mod4$AIC,mod4$BIC),
c(mod5$AIC,mod5$BIC),
c(mod6$AIC,mod6$BIC))
rownames(Comparacao)<-c("Modelo completo","Sem V1","Sem V3", "Sem V4", "Sem V6", "Sem V8")
colnames(Comparacao)<-c("AIC","BIC")
Comparacao## AIC BIC
## Modelo completo 442.5370 466.6221
## Sem V1 442.4092 463.4836
## Sem V3 444.9025 465.9769
## Sem V4 449.3428 470.4173
## Sem V6 446.0072 467.0817
## Sem V8 444.5520 465.6264Assim, de acordo com esse critério, a variável V1 é retirada nesse passo. Observe que ela é a única que diminui o AIC e BIC do modelo completo.
5.1 PASSO 3
Ajusta-se um novo modelo com as covariáveis retidas no passo 2. Neste passo, retorna-se com as covariáveis excluídas no passo 1 para confirmar que elas não são devem permanecer no modelo. Utilizamos, novamente, o AIC e BIC para essa decisão.
mod1<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8,dist="gengamma.orig")
mod2<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8+V2,dist="gengamma.orig")
mod3<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8+V5,dist="gengamma.orig")
mod4<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8+V7,dist="gengamma.orig")
mod5<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8+V9,dist="gengamma.orig")
mod6<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8+V10,dist="gengamma.orig")
mod7<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8+V11,dist="gengamma.orig")
Comparacao<-rbind(c(mod1$AIC,mod1$BIC),
c(mod2$AIC,mod2$BIC),
c(mod3$AIC,mod3$BIC),
c(mod4$AIC,mod4$BIC),
c(mod5$AIC,mod5$BIC),
c(mod6$AIC,mod6$BIC),
c(mod7$AIC,mod7$BIC))
rownames(Comparacao)<-c("Modelo completo","Com V2","Com V5", "Com V7", "Com V9", "Com V10", "Com V11")
colnames(Comparacao)<-c("AIC","BIC")
Comparacao## AIC BIC
## Modelo completo 442.4092 463.4836
## Com V2 443.3129 467.3979
## Com V5 444.4153 468.5004
## Com V7 442.6437 466.7287
## Com V9 443.3263 467.4114
## Com V10 444.1726 468.2577
## Com V11 444.2530 468.3381De acordo com esse passo, nenhuma variável retorna ao modelo.
5.2 PASSO 4
Ajusta-se um modelo incluindo-se as covariáveis retidas no passo 3. Neste passo, é testado se alguma delas pode ser retirada do modelo. Utilizamos, novamente, o AIC e BIC para essa decisão.
mod1<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8,dist="gengamma.orig")
mod2<-flexsurvreg(data=dados,Surv(tempo,cens)~V4+V6+V8,dist="gengamma.orig")
mod3<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V6+V8,dist="gengamma.orig")
mod4<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V8,dist="gengamma.orig")
mod5<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6,dist="gengamma.orig")
Comparacao<-rbind(c(mod1$AIC,mod1$BIC),
c(mod2$AIC,mod2$BIC),
c(mod3$AIC,mod3$BIC),
c(mod4$AIC,mod4$BIC),
c(mod5$AIC,mod5$BIC))
rownames(Comparacao)<-c("Modelo completo","Sem V3","Sem V4", "Sem V6", "Sem V8")
colnames(Comparacao)<-c("AIC","BIC")
Comparacao## AIC BIC
## Modelo completo 442.4092 463.4836
## Sem V3 444.2275 462.2913
## Sem V4 453.5129 471.5767
## Sem V6 445.5023 463.5661
## Sem V8 445.9702 464.0340Assim, as variáveis permanecem no modelo.
6 PASSO 5
Utilizando as covariáveis que sobreviveram ao passo 4, ajusta-se o modelo final para os efeitos principais. Para completar a modelagem, deve-se verificar a possibilidade de inclusão de termos de interação dupla entre as covariáveis inclufdas no modelo. O modelo final fica determinado pelos efeitos principais identificados no passo 4 e os termos de interação significativos identificados neste passo. Utilizamos, novamente, o AIC e BIC para essa decisão.
mod1<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8,dist="gengamma.orig")
mod2<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8+V3*V4,dist="gengamma.orig")
mod3<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8+V3*V6,dist="gengamma.orig")
mod4<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8+V3*V8,dist="gengamma.orig")
mod5<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8+V4*V6,dist="gengamma.orig")
mod6<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8+V4*V8,dist="gengamma.orig")
mod7<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8+V6*V8,dist="gengamma.orig")
Comparacao<-rbind(c(mod1$AIC,mod1$BIC),
c(mod2$AIC,mod2$BIC),
c(mod3$AIC,mod3$BIC),
c(mod4$AIC,mod4$BIC),
c(mod5$AIC,mod5$BIC),
c(mod6$AIC,mod6$BIC),
c(mod7$AIC,mod7$BIC))
rownames(Comparacao)<-c("Modelo completo","Com V3_V4","Com V3_V6",
"Com V3_V8","Com V4_V6","Com V4_V8","Com V6_V8")
colnames(Comparacao)<-c("AIC","BIC")
Comparacao## AIC BIC
## Modelo completo 442.4092 463.4836
## Com V3_V4 443.6904 467.7755
## Com V3_V6 442.6568 466.7418
## Com V3_V8 444.1216 468.2067
## Com V4_V6 443.6785 467.7636
## Com V4_V8 443.6650 467.7501
## Com V6_V8 444.1542 468.2393Nenhuma interação deve ser adicionada. Portanto, o modelo final inclui as variáveis V3, V4, V6 e v8.
7 ESCOLHA DO MODELO PROBABILÍSTICO
Após os passos mencionados anteriormente, teremos o ajuste final do modelo de regressão gama generalizado. Na sequência, trocamos apenas a distribuição gama generalizada pelas distribuições:
- Exponencial
- Weibull
- LogNormal
O melhor modelo será aquele que possuir os menores valores de AIC e BIC. Para o melhor modelo, é feita a análise de resíduos e interpretação dos coeficientes estimados, como já fizemos anteriormente para modelos mais simples.
mod_gengamma<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8,dist="gengamma.orig")
mod_exp<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8,dist="exponential")
mod_wei<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8,dist="weibull")
mod_lognorm<-flexsurvreg(data=dados,Surv(tempo,cens)~V3+V4+V6+V8,dist="lnorm")
Comparacao<-rbind(c(mod_gengamma$AIC,mod_gengamma$BIC),
c(mod_exp$AIC,mod_exp$BIC),
c(mod_wei$AIC,mod_wei$BIC),
c(mod_lognorm$AIC,mod_lognorm$BIC))
rownames(Comparacao)<-c("Gama Generalizada", "Exponencial", "Weibull", "LogNormal")
colnames(Comparacao)<-c("AIC","BIC")
Comparacao## AIC BIC
## Gama Generalizada 442.4092 463.4836
## Exponencial 442.9825 458.0357
## Weibull 444.8711 462.9349
## LogNormal 440.8734 458.9372De acordo com esse critério, escolhemos o modelo final com a distribuição LogNormal.
8 ANÁLISE DE RESÍDUOS
Seguimos com a análise de resíduos do modelo escolhido:
CoxSnell <-coxsnell_flexsurvreg(mod_lognorm)$estPodemos fazer dois gráficos:
- Gráfico 1:
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")
- 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")
De acordo com os gráficos, o modelo está bem ajustado.
9 INTERPRETAÇÃO DOS COEFICIENTES
Ao retornar no modelo lognormal, temos que:
mod_lognorm## Call:
## flexsurvreg(formula = Surv(tempo, cens) ~ V3 + V4 + V6 + V8,
## data = dados, dist = "lnorm")
##
## Estimates:
## data mean est L95% U95% se
## meanlog NA 2.7623 2.0595 3.4650 0.3586
## sdlog NA 1.4217 1.1914 1.6966 0.1282
## V3menor ou igual a 6 meses 0.4667 -0.5693 -1.1272 -0.0115 0.2846
## V4sim 0.4400 -1.0459 -1.6038 -0.4881 0.2846
## V6sim 0.6267 0.6766 0.1098 1.2434 0.2892
## V8menor que 1 SM 0.4867 -0.6596 -1.2131 -0.1061 0.2824
## exp(est) L95% U95%
## meanlog NA NA NA
## sdlog NA NA NA
## V3menor ou igual a 6 meses 0.5659 0.3239 0.9886
## V4sim 0.3514 0.2011 0.6138
## V6sim 1.9672 1.1161 3.4673
## V8menor que 1 SM 0.5170 0.2973 0.8993
##
## N = 150, Events: 65, Censored: 85
## Total time at risk: 820.8
## Log-likelihood = -214.4367, df = 6
## AIC = 440.8734- Váriavel V3: Conceito materno sobre o tempo ideal de amamentação. Temos que \(exp\{\hat{\beta}_{\leq 6}\} = 0,56\). Assim, comparado com o grupo de mães que acreditam que o tempo de amamentação é > 6 meses, o tempo mediano de desmame das maes que acreditam que o tempo de amamentação é \(\leq 6\) é 44% menor.
- Váriavel V4: Dificuldade para amamentar nos primeiros dias pós-parto. Temos que \(exp\{\hat{\beta}_{sim}\} = 0,35\). Assim, comparado com o grupo de mães que não tiverarm dificuldade de amamentação, o tempo mediano de desmame das mães que tiveram dificuldade de amamentação é 65% menor.
- Variável V6: Recebeu exclusivamente leite materno na maternidade. Temos que \(exp\{\hat{\beta}_{sim}\} = 1,96\). Assim, comparado com o grupo de mães que não receberam exclusivamente leite materno na maternidade, o tempo mediano de desmame das mães que receberam é 1,96 vezes maior.
- Variável V8: Renda per capita. Temos que \(exp\{\hat{\beta}_{< 1 }\} = 0,51\). Assim, comparado com o grupo de mães que possuem mais de 1 SM, o tempo mediano de desmame das mães que possuem menos de 1 SM é 49% menor.