Análise de Sobrevivência - Projeto 02
1 Modelos de regressão
1.1 Descrição do estudo e das variáveis
library(readr)
library(survival)
library(flexsurv)
library(knitr)
library(ggplot2)
library(RColorBrewer)
library(pacman)
library(tidyverse)
library(scales)
pacman::p_load(knitr, captioner, bundesligR, stringr)
table_nums <- captioner::captioner(prefix = "Tabela")
tab.1_cap <- table_nums(name = "tab_1",
caption = "Descrição das variáveis da base de dados")
Drogas_SP <- read_csv("Drogas-SP.csv")[,-5]Os dados são provenientes da Pesquisa Nacional de Saúde do Escolar (PeNSE) de 2015, conduzida pelo IBGE, e se referem ao tempo até a experimentação de drogas ilícitas em uma amostra de estudantes do 9º ano do ensino fundamental de escolas privadas da cidade de São Paulo.
Objetivo: modelar o tempo até a experimentação de drogas ilícitas utilizando técnicas de análise de sobrevivência.
Tabela 1: Descrição das variáveis da base de dados
desc_variaveis <- data.frame(Variável = colnames(Drogas_SP),
Renomeando = c("tempo",
"censura",
"sexo",
"cor_raca",
"mora_mae",
"tempo_livre_responsaveis",
"bebida_alcoolica",
"cigarro",
"responsavel_fuma",
"sentido_sozinho",
"amigos_drogas"),
Descrição = c("Tempo até o primeiro uso de drogas ilícitas, como: maconha, cocaína, crack, loló, lança-perfume, ecstasy, oxy, etc (variável resposta)",
"Indicador de censura ",
"Qual é o seu sexo?",
"Qual é a sua cor ou raça? ",
"Você mora com sua mãe? ",
"NOS ÚLTIMOS 30 DIAS, com que frequência seus pais ou responsáveis sabiam realmente o que você estava fazendo em seu tempo livre? ",
"Alguma vez na vida você tomou uma dose de bebida alcoólica? (Uma dose equivale a uma lata de cerveja ou uma taça de vinho ou uma dose de cachaça ou uísque etc) ",
"Alguma vez na vida, você já fumou cigarro, mesmo uma ou duas tragadas? ",
"Algum de seus pais ou responsáveis fuma? ",
"NOS ÚLTIMOS 12 MESES com que frequência tem se sentido sozinho(a)? ",
"Quantos amigos seus usam drogas? "),
Categorias =c("-","0=censura e 1=falha","1=Masculino e 2=Feminino","1=Branca, 2=Preta, 3=Amarela, 4=Parda, 5=Indígena e 99=Não informado","1=Sim, 2=Não e 99=Não informado","1=Nunca, 2=Raramente, 3=Às vezes, 4=Na maior parte do tempo, 5=Sempre e 99=Não informado","1=Sim, 2=Não e 99=Não informado", "1=Sim, 2=Não e 99=Não informado","1=Nenhum deles, 2=Só meu pai ou responsável do sexo masculino, 3=Só minha mãe ou responsável do sexo feminino, 4=Meu pai e minha mãe ou responsáveis, 5=Não sei e 99=Não informado","1=Nunca, 2=Raramente, 3=Às vezes, 4=Na maioria das vezes, 5=Sempre e 99=Não informado","1=Nenhum, 2=Poucos, 3=Alguns, 4=A maioria, 5=Todos, 6=Não sei e 99=Não informado"))
knitr::kable(desc_variaveis)| Variável | Renomeando | Descrição | Categorias |
|---|---|---|---|
| Tempo | tempo | Tempo até o primeiro uso de drogas ilícitas, como: maconha, cocaína, crack, loló, lança-perfume, ecstasy, oxy, etc (variável resposta) | - |
| Censura | censura | Indicador de censura | 0=censura e 1=falha |
| VB01001 | sexo | Qual é o seu sexo? | 1=Masculino e 2=Feminino |
| VB01002 | cor_raca | Qual é a sua cor ou raça? | 1=Branca, 2=Preta, 3=Amarela, 4=Parda, 5=Indígena e 99=Não informado |
| VB01006 | mora_mae | Você mora com sua mãe? | 1=Sim, 2=Não e 99=Não informado |
| VB07002 | tempo_livre_responsaveis | NOS ÚLTIMOS 30 DIAS, com que frequência seus pais ou responsáveis sabiam realmente o que você estava fazendo em seu tempo livre? | 1=Nunca, 2=Raramente, 3=Às vezes, 4=Na maior parte do tempo, 5=Sempre e 99=Não informado |
| VB05002 | bebida_alcoolica | Alguma vez na vida você tomou uma dose de bebida alcoólica? (Uma dose equivale a uma lata de cerveja ou uma taça de vinho ou uma dose de cachaça ou uísque etc) | 1=Sim, 2=Não e 99=Não informado |
| VB04001 | cigarro | Alguma vez na vida, você já fumou cigarro, mesmo uma ou duas tragadas? | 1=Sim, 2=Não e 99=Não informado |
| VB04006A | responsavel_fuma | Algum de seus pais ou responsáveis fuma? | 1=Nenhum deles, 2=Só meu pai ou responsável do sexo masculino, 3=Só minha mãe ou responsável do sexo feminino, 4=Meu pai e minha mãe ou responsáveis, 5=Não sei e 99=Não informado |
| VB12001 | sentido_sozinho | NOS ÚLTIMOS 12 MESES com que frequência tem se sentido sozinho(a)? | 1=Nunca, 2=Raramente, 3=Às vezes, 4=Na maioria das vezes, 5=Sempre e 99=Não informado |
| VB06006 | amigos_drogas | Quantos amigos seus usam drogas? | 1=Nenhum, 2=Poucos, 3=Alguns, 4=A maioria, 5=Todos, 6=Não sei e 99=Não informado |
colnames(Drogas_SP)<- c("tempo",
"censura",
"sexo",
"cor_raca",
"mora_mae",
"tempo_livre_responsaveis",
"bebida_alcoolica",
"cigarro",
"responsavel_fuma",
"sentido_sozinho",
"amigos_drogas")A base de dados é composta por 445 observações e 13 variáveis: as 9 covariáveis e a variável resposta que é representada pelo tempo até o primeiro uso de drogas ilícitas e uma variável indicadora de ocorrência de censura.
Observação: a variável V0007 foi removida da base por não conter descrição e conter somente uma categoria (2).
1.2 Estimação de S(t) por Kaplan-Meier para cada uma das covariáveis
Drogas_SP <- Drogas_SP %>% mutate(sexo_cat = ifelse(sexo==1,"Masculino","Feminino"),
cor_raca_cat = ifelse(cor_raca ==1, "Branca",ifelse(cor_raca==2,"Preta",ifelse(cor_raca==3,"Amarela",ifelse(cor_raca==4,"Parda","Indígena")))),
mora_mae_cat = ifelse(mora_mae == 1,"Sim","Não"),
tempo_livre_responsaveis_cat = ifelse(tempo_livre_responsaveis==1,"Nunca",ifelse(tempo_livre_responsaveis==2,"Raramente",ifelse(tempo_livre_responsaveis==3,"Às vezes",ifelse(tempo_livre_responsaveis==4,"Na maior parte do tempo","Sempre")))),
bebida_alcoolica_cat = ifelse(bebida_alcoolica==1,"Sim","Não"),
cigarro_cat = ifelse(cigarro==1,"Sim","Não"),
responsavel_fuma_cat = ifelse(responsavel_fuma==1,"Nenhum deles",ifelse(responsavel_fuma==2,"Só meu pai ou responsável do sexo masculino",ifelse(responsavel_fuma==3,"Só minha mãe ou responsável do sexo feminino",ifelse(responsavel_fuma==4,"Meu pai e minha mãe ou responsáveis",ifelse(responsavel_fuma==5,"Não sei","Não informado"))))),
sentido_sozinho_cat = ifelse(sentido_sozinho==1,"Nunca",ifelse(sentido_sozinho==2,"Raramente",ifelse(sentido_sozinho==3,"Às vezes",ifelse(sentido_sozinho==4,"Na maioria das vezes","Sempre")))),
amigos_drogas_cat = ifelse(amigos_drogas == 1,"Nenhum",ifelse(amigos_drogas==2,"Poucos",ifelse(amigos_drogas==3,"Alguns",ifelse(amigos_drogas==4,"A maioria",ifelse(amigos_drogas==5,"Todos","Não sei"))))))
tempos<- Drogas_SP$tempo
cens<-Drogas_SP$censura
ekmV1<-survfit(Surv(tempos,cens)~Drogas_SP$sexo_cat)
ekmV2<-survfit(Surv(tempos,cens)~Drogas_SP$cor_raca_cat)
ekmV3<-survfit(Surv(tempos,cens)~Drogas_SP$mora_mae_cat)
ekmV4<-survfit(Surv(tempos,cens)~Drogas_SP$tempo_livre_responsaveis_cat)
ekmV5<-survfit(Surv(tempos,cens)~Drogas_SP$bebida_alcoolica_cat)
ekmV6<-survfit(Surv(tempos,cens)~Drogas_SP$cigarro_cat)
ekmV7<-survfit(Surv(tempos,cens)~Drogas_SP$responsavel_fuma_cat)
ekmV8<-survfit(Surv(tempos,cens)~Drogas_SP$sentido_sozinho_cat)
ekmV9<-survfit(Surv(tempos,cens)~Drogas_SP$amigos_drogas_cat )
plot(ekmV1,lty=c(1,2),xlab="Tempo (anos)",ylab="S(t) estimada",main="sexo")
legend("bottomleft",lty=c(1,2),c("Feminino","Masculino"),lwd=1,bty="n")
- A partir dos 15 anos, a curva de sobrevivência do sexo feminino possui
queda abrupta, ficando abaixo da curva do sexo masculino. Isso nos diz
que a partir dessa idade, as probabilidades de sobrevivência para o sexo
feminino são inferiores às probabilidades de sobrevivência do sexo
masculino
plot(ekmV2,lty=c(1,2,3,4,5),xlab="Tempo (anos)",ylab="S(t) estimada",main="cor_raca", col=c(1,2,3,4,5))
legend("bottomleft",lty=c(1,2,3,4,5),c("Amarela","Branca","Indígena","Parda","Preta"),lwd=1,bty="n", col=c(1,2,3,4,5))
- A curva de sobrevivência por Kaplan-Meier da cor branca é a primeira a
apresentar queda na probabilidade de sobrevivência.
plot(ekmV3,lty=c(1,2),xlab="Tempo (anos)",ylab="S(t) estimada",main="mora_mae")
legend("bottomleft",lty=c(1,2),c("Não","Sim"),lwd=1,bty="n")- A diferença entre as probabilidades de sobrevivência dos estudantes que não moram com a mãe e dos que moram com a mãe é maior ao se aproximar dos 15 anos.
plot(ekmV4,lty=c(1,2,3,4,5),xlab="Tempo (anos)",ylab="S(t) estimada",main="tempo_livre_responsaveis",col=c(1,2,3,4,5))
legend("bottomleft",lty=c(1,2,3,4,5),c("Às vezes","Na maior parte do tempo","Nunca","Raramente","Sempre"),lwd=1,bty="n",col=c(1,2,3,4,5))- As curvas de sobrevivência das categorias Às vezes, Na maior parte do tempo e Sempre ficam acima das curvas das categorias Nunca e Raramente em grande parte dos anos.
plot(ekmV5,lty=c(1,2),xlab="Tempo (anos)",ylab="S(t) estimada",main="bebida_alcoolica")
legend("bottomleft",lty=c(1,2),c("Não","Sim"),lwd=1,bty="n")
- A curva de sobrevivência dos estudantes que alguma vez na vida você
tomaram uma dose de bebida alcoólica possui mais falhas, bem como decai
mais cedo e rapidamente que a curva dos estudantes que não tomaram.
- As probabilidades de sobrevivência para os estudantes que tomaram bebida alcóolica é menor.
plot(ekmV6,lty=c(1,2),xlab="Tempo (anos)",ylab="S(t) estimada",main="cigarro")
legend("bottomleft",lty=c(1,2),c("Não","Sim"),lwd=1,bty="n")- A curva de sobrevivência dos estudantes que alguma vez na vida você fumaram cigarro possui mais falhas decai mais rapidamente que a curva dos estudantes que não fumaram. As probabilidades de sobrevivência para os estudantes que tomaram bebida alcóolica é menor.
plot(ekmV7,lty=c(1,2,3,4,5,6),xlab="Tempo (anos)",ylab="S(t) estimada",main="responsavel_fuma", col=c(1,2,3,4,5,6))
legend("bottomleft",lty=c(1,2,3,4,5,6),c("Meu pai e minha mãe ou responsáveis","Não informado","Não sei","Nenhum deles","Só meu pai ou responsável do sexo masculino","Só minha mãe ou responsável do sexo feminino"),lwd=1,bty="n", col=c(1,2,3,4,5,6))
- As menores probabilidades de sobrevivência são observadas para os
estudantes que não sabem se o responsável fuma, seguido dos estudantes
que declararam que só a mãe ou responsável do sexo feminino fuma.
plot(ekmV8,lty=c(1,2,3,4,5),xlab="Tempo (anos)",ylab="S(t) estimada",main="sentido_sozinho", col=c(1,2,3,4,5))
legend("bottomleft",lty=c(1,2,3,4,5),c("Às vezes","Na maioria das vezes","Nunca","Raramente","Sempre"),lwd=1,bty="n", col=c(1,2,3,4,5))- As curvas dos estudantes que declararam que nos últimos 12 meses se sentiram sozinhos sempre ou na maioria das vezes estão abaixo das demais, o que indica menores probabilidades de sobrevivência.
plot(ekmV9,lty=c(1,2,3,4,5),xlab="Tempo (anos)",ylab="S(t) estimada",main="amigos_drogas", col=c(1,2,3,4,5))
legend("bottomleft",lty=c(1,2,3,4,5),c("A maioria","Alguns","Não sei","Nenhum","Poucos"),lwd=1,bty="n", col=c(1,2,3,4,5))
- A curva de sobrevivência dos estudantes que declararam que a maioria
dos amigos usam drogas destoa das demais, possuindo queda abrupta.
1.3 Teste Logrank para cada uma das covariáveis
logrank1 <- survdiff(Surv(tempos,cens)~Drogas_SP$sexo_cat,rho=0)
logrank2 <- survdiff(Surv(tempos,cens)~Drogas_SP$cor_raca_cat,rho=0)
logrank3 <- survdiff(Surv(tempos,cens)~Drogas_SP$mora_mae_cat,rho=0)
logrank4 <- survdiff(Surv(tempos,cens)~Drogas_SP$tempo_livre_responsaveis_cat,rho=0)
logrank5 <- survdiff(Surv(tempos,cens)~Drogas_SP$bebida_alcoolica_cat,rho=0)
logrank6 <- survdiff(Surv(tempos,cens)~Drogas_SP$cigarro_cat,rho=0)
logrank7 <- survdiff(Surv(tempos,cens)~Drogas_SP$responsavel_fuma_cat,rho=0)
logrank8 <- survdiff(Surv(tempos,cens)~Drogas_SP$sentido_sozinho_cat,rho=0)
logrank9 <- survdiff(Surv(tempos,cens)~Drogas_SP$amigos_drogas_cat,rho=0)
lr_pv = data.frame(
Covariaveis = c("sexo",
"cor_raca",
"mora_mae",
"tempo_livre_respon",
"bebida_alcoolica",
"cigarro",
"responsavel_fuma",
"sentido_sozinho" ,
"amigos_drogas"),
"P-Valores"= round(c(1-pchisq(logrank1$chisq,1),
1-pchisq(logrank2$chisq,4),
1-pchisq(logrank3$chisq,1),
1-pchisq(logrank4$chisq,4),
1-pchisq(logrank5$chisq,1),
1-pchisq(logrank6$chisq,1),
1-pchisq(logrank7$chisq,5),
1-pchisq(logrank8$chisq,4),
1-pchisq(logrank9$chisq,4)),4),
Categorias=c("1=Masculino e 2=Feminino",
"1=Branca, 2=Preta, 3=Amarela, 4=Parda, 5=Indígena",
"1=Sim, 2=Não",
"1=Nunca, 2=Raramente, 3=Às vezes, 4=Na maior parte do tempo, 5=Sempre",
"1=Sim, 2=Não",
"1=Sim, 2=Não",
"1=Nenhum deles, 2=Só meu pai ou responsável do sexo masculino, 3=Só minha mãe ou responsável do sexo feminino, 4=Meu pai e minha mãe ou responsáveis, 5=Não sei e 99=Não informado",
"1=Nunca, 2=Raramente, 3=Às vezes, 4=Na maioria das vezes,
5=Sempre",
"1=Nenhum, 2=Poucos, 3=Alguns, 4=A maioria, 5=Todos, 6=Não sei"),
Resultado = ifelse(c(1-pchisq(logrank1$chisq,1),
1-pchisq(logrank2$chisq,4),
1-pchisq(logrank3$chisq,1),
1-pchisq(logrank4$chisq,4),
1-pchisq(logrank5$chisq,1),
1-pchisq(logrank6$chisq,1),
1-pchisq(logrank7$chisq,5),
1-pchisq(logrank8$chisq,4),
1-pchisq(logrank9$chisq,4))<=0.1,"Ao nível de significancia de 0.1, rejeitamos a hipótese nula
de igualdade entre as curvas de sobrevivência das categorias da variável","Ao nível de significancia de 0.1, não rejeitamos a hipótese nula
de igualdade entre as curvas de sobrevivência das categorias da variável"))
knitr::kable(lr_pv,col.names =c("Covariáveis","P-valores", "Categorias","Resultados"))| Covariáveis | P-valores | Categorias | Resultados |
|---|---|---|---|
| sexo | 0.5581 | 1=Masculino e 2=Feminino | Ao nível de significancia de 0.1, não rejeitamos a hipótese nula |
| de igualdade entre as curvas de sobrevivência das categorias da variável | |||
| cor_raca | 0.9660 | 1=Branca, 2=Preta, 3=Amarela, 4=Parda, 5=Indígena | Ao nível de significancia de 0.1, não rejeitamos a hipótese nula |
| de igualdade entre as curvas de sobrevivência das categorias da variável | |||
| mora_mae | 0.6251 | 1=Sim, 2=Não | Ao nível de significancia de 0.1, não rejeitamos a hipótese nula |
| de igualdade entre as curvas de sobrevivência das categorias da variável | |||
| tempo_livre_respon | 0.0036 | 1=Nunca, 2=Raramente, 3=Às vezes, 4=Na maior parte do tempo, 5=Sempre | Ao nível de significancia de 0.1, rejeitamos a hipótese nula |
| de igualdade entre as curvas de sobrevivência das categorias da variável | |||
| bebida_alcoolica | 0.0000 | 1=Sim, 2=Não | Ao nível de significancia de 0.1, rejeitamos a hipótese nula |
| de igualdade entre as curvas de sobrevivência das categorias da variável | |||
| cigarro | 0.0000 | 1=Sim, 2=Não | Ao nível de significancia de 0.1, rejeitamos a hipótese nula |
| de igualdade entre as curvas de sobrevivência das categorias da variável | |||
| responsavel_fuma | 0.3852 | 1=Nenhum deles, 2=Só meu pai ou responsável do sexo masculino, 3=Só minha mãe ou responsável do sexo feminino, 4=Meu pai e minha mãe ou responsáveis, 5=Não sei e 99=Não informado | Ao nível de significancia de 0.1, não rejeitamos a hipótese nula |
| de igualdade entre as curvas de sobrevivência das categorias da variável | |||
| sentido_sozinho | 0.0090 | 1=Nunca, 2=Raramente, 3=Às vezes, 4=Na maioria das vezes, | |
| 5=Sempre | Ao nível de significancia de 0.1, rejeitamos a hipótese nula | ||
| de igualdade entre as curvas de sobrevivência das categorias da variável | |||
| amigos_drogas | 0.0000 | 1=Nenhum, 2=Poucos, 3=Alguns, 4=A maioria, 5=Todos, 6=Não sei | Ao nível de significancia de 0.1, rejeitamos a hipótese nula |
| de igualdade entre as curvas de sobrevivência das categorias da variável |
1.4 Seleção de covariáveis para o modelo
1.4.1 Passo 1
- Passo 1: Ajustar modelos de regressão univariados usando a distribuição Gama Generalizada
Drogas_SP_reg <- Drogas_SP %>% mutate (V1 = as.factor(sexo),
V2 = as.factor(cor_raca),
V3 = as.factor(mora_mae),
V4 = as.factor(tempo_livre_responsaveis),
V5 = as.factor(bebida_alcoolica),
V6 = as.factor(cigarro),
V7 = as.factor(responsavel_fuma),
V8 = as.factor(sentido_sozinho),
V9 = as.factor(amigos_drogas)) %>%
select(tempo, censura, V1,V2,V3,V4,V5,V6,V7,V8,V9)
m0<-flexsurvreg(formula=Surv(tempo,censura)~1,dist='gengamma',data=Drogas_SP) #modelo nulo
m1<-flexsurvreg(formula=Surv(tempo,censura)~V1,dist='gengamma',data=Drogas_SP_reg) #modelo com V1
m2<-flexsurvreg(formula=Surv(tempo,censura)~V2,dist='gengamma',data=Drogas_SP_reg) #modelo com V2
m3<-flexsurvreg(formula=Surv(tempo,censura)~V3,dist='gengamma',data=Drogas_SP_reg) #modelo com V3
m4<-flexsurvreg(formula=Surv(tempo,censura)~V4,dist='gengamma',data=Drogas_SP_reg) #modelo com V4
m5<-flexsurvreg(formula=Surv(tempo,censura)~V5,dist='gengamma',data=Drogas_SP_reg) #modelo com V5
m6<-flexsurvreg(formula=Surv(tempo,censura)~V6,dist='gengamma',data=Drogas_SP_reg) #modelo com V6
m7<-flexsurvreg(formula=Surv(tempo,censura)~V7,dist='gengamma',data=Drogas_SP_reg) #modelo com V7
m8<-flexsurvreg(formula=Surv(tempo,censura)~V8,dist='gengamma',data=Drogas_SP_reg) #modelo com V8
m9<-flexsurvreg(formula=Surv(tempo,censura)~V9,dist='gengamma',data=Drogas_SP_reg) #modelo com V9
TRV1=2*(m1$loglik-m0$loglik)
TRV2=2*(m2$loglik-m0$loglik)
TRV3=2*(m3$loglik-m0$loglik)
TRV4=2*(m4$loglik-m0$loglik)
TRV5=2*(m5$loglik-m0$loglik)
TRV6=2*(m6$loglik-m0$loglik)
TRV7=2*(m7$loglik-m0$loglik)
TRV8=2*(m8$loglik-m0$loglik)
TRV9=2*(m9$loglik-m0$loglik)
selecao1 <- data.frame(Modelo=c("Nulo","V1: sexo","V2: cor_raca","V3: mora_mae","V4: tempo_livre_responsaveis","V5: bebida_alcoolica","V6: cigarro","V7: responsavel_fuma","V8: sentido_sozinho","V9: amigos_drogas"),
TRV = c("-",round(c(TRV1,TRV2,TRV3,TRV4,TRV5,TRV6,TRV7,TRV8,TRV9),4)),
pvalor = c("-", round(c(1-pchisq(TRV1,df=1), 1-pchisq(TRV2,df=1), 1-pchisq(TRV3,df=1), 1-pchisq(TRV4,df=1), 1-pchisq(TRV5,df=1), 1-pchisq(TRV6,df=1), 1-pchisq(TRV7,df=1), 1-pchisq(TRV8,df=1), 1-pchisq(TRV9,df=1)),4)),
Significancia = c("-",ifelse(round(c(1-pchisq(TRV1,df=1), 1-pchisq(TRV2,df=1), 1-pchisq(TRV3,df=1), 1-pchisq(TRV4,df=1), 1-pchisq(TRV5,df=1), 1-pchisq(TRV6,df=1), 1-pchisq(TRV7,df=1), 1-pchisq(TRV8,df=1), 1-pchisq(TRV9,df=1)),4)
<=0.1,"Ao nível de 10%, É significativa","Ao nível de 10%, NÃO É significativa")))
knitr::kable(selecao1)| Modelo | TRV | pvalor | Significancia |
|---|---|---|---|
| Nulo | - | - | - |
| V1: sexo | 0.5913 | 0.4419 | Ao nível de 10%, NÃO É significativa |
| V2: cor_raca | 0.6609 | 0.4162 | Ao nível de 10%, NÃO É significativa |
| V3: mora_mae | 0.0393 | 0.8429 | Ao nível de 10%, NÃO É significativa |
| V4: tempo_livre_responsaveis | 15.8959 | 1e-04 | Ao nível de 10%, É significativa |
| V5: bebida_alcoolica | 41.9055 | 0 | Ao nível de 10%, É significativa |
| V6: cigarro | 71.8734 | 0 | Ao nível de 10%, É significativa |
| V7: responsavel_fuma | 3.7302 | 0.0534 | Ao nível de 10%, É significativa |
| V8: sentido_sozinho | 10.8572 | 0.001 | Ao nível de 10%, É significativa |
| V9: amigos_drogas | 89.3055 | 0 | Ao nível de 10%, É significativa |
1.4.2 Passo 2
Passo 2: Ajustar conjuntamente todas as cováriaveis significativas ao nível de significância de \(0.1\) no passo anterior e, em seguida, ajustar modelos reduzidos eliminando uma covariável de cada vez.
Ao nível de significância de \(0.1\), as covariáveis significativas no passo 1 foram tempo_livre_responsaveis (V4), bebida_alcoolica (V5), cigarro (V6), responsavel_fuma(V7), sentido_sozinho (V8) e amigos_drogas (V9). Portanto, estas covariáveis compõem o modelo completo.
m.comp <-flexsurvreg(formula=Surv(tempo,censura)~V4+V5+V6+V7+V8+V9,dist='gengamma',data=Drogas_SP_reg) #modelo completo
m.1 <- flexsurvreg(formula=Surv(tempo,censura)~V5+V6+V7+V8+V9,dist='gengamma',data=Drogas_SP_reg) #modelo sem V4
m.2 <- flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V7+V8+V9,dist='gengamma',data=Drogas_SP_reg) #modelo sem V5
m.3 <- flexsurvreg(formula=Surv(tempo,censura)~V4+V5+V7+V8+V9,dist='gengamma',data=Drogas_SP_reg) #modelo sem V6
m.4 <- flexsurvreg(formula=Surv(tempo,censura)~V4+V5+V6+V8+V9,dist='gengamma',data=Drogas_SP_reg) #modelo sem V7
m.5 <- flexsurvreg(formula=Surv(tempo,censura)~V4+V5+V6+V7+V9,dist='gengamma',data=Drogas_SP_reg) #modelo sem V8
m.6 <- flexsurvreg(formula=Surv(tempo,censura)~V4+V5+V6+V7+V8,dist='gengamma',data=Drogas_SP_reg) #modelo sem V9
TRV.1=2*(m.comp$loglik-m.1$loglik) #TRV comparando o modelo sem V4 e o modelo completo (H0:beta_V4=0)
TRV.2=2*(m.comp$loglik-m.2$loglik) #TRV comparando o modelo sem V5 e o modelo completo (H0:beta_V5=0)
TRV.3=2*(m.comp$loglik-m.3$loglik) #TRV comparando o modelo sem V6 e o modelo completo (H0:beta_V6=0)
TRV.4=2*(m.comp$loglik-m.4$loglik) #TRV comparando o modelo sem V7 e o modelo completo (H0:beta_V7=0)
TRV.5=2*(m.comp$loglik-m.5$loglik) #TRV comparando o modelo sem V8 e o modelo completo (H0:beta_V8=0)
TRV.6=2*(m.comp$loglik-m.6$loglik) #TRV comparando o modelo sem V9 e o modelo completo (H0:beta_V9=0)
selecao2 <- data.frame(Modelo=c("Completo: V4+V5+V6+V8+V9","Tirando tempo_livre_responsaveis (V4): V5+V6+V8+V9 "," Tirando bebida_alcoolica (V5): V4+V6+V7+V8+V9"," Tirando cigarro(V6): V4+V5+V8+V9 ","Tirando responsavel_fuma (V7) : V4+V5+V6+V8+V9","Tirando sentido_sozinho (V8): ","Tirando amigos_drogas (V9): V4+V5+V6+V8"),
TRV = c("-",round(c(TRV.1,TRV.2,TRV.3,TRV.4,TRV.5,TRV.6),4)),
pvalor = c("-", round(c(1-pchisq(TRV.1,df=1),1-pchisq(TRV.2,df=1),1-pchisq(TRV.3,df=1),1-pchisq(TRV.4,df=1),1-pchisq(TRV.5,df=1),1-pchisq(TRV.6,df=1)),4)),
Significancia = c("-",ifelse(round(c(1-pchisq(TRV.1,df=1), 1-pchisq(TRV.2,df=1), 1-pchisq(TRV.3,df=1), 1-pchisq(TRV.4,df=1), 1-pchisq(TRV.5,df=1),1-pchisq(TRV.6,df=1)),4)
<=0.1,"Ao nível de 10%, É significativa","Ao nível de 10%, NÃO É significativa")))
knitr::kable(selecao2)| Modelo | TRV | pvalor | Significancia |
|---|---|---|---|
| Completo: V4+V5+V6+V8+V9 | - | - | - |
| Tirando tempo_livre_responsaveis (V4): V5+V6+V8+V9 | 4.6341 | 0.0313 | Ao nível de 10%, É significativa |
| Tirando bebida_alcoolica (V5): V4+V6+V7+V8+V9 | 0.9119 | 0.3396 | Ao nível de 10%, NÃO É significativa |
| Tirando cigarro(V6): V4+V5+V8+V9 | 4.5017 | 0.0339 | Ao nível de 10%, É significativa |
| Tirando responsavel_fuma (V7) : V4+V5+V6+V8+V9 | 1.3784 | 0.2404 | Ao nível de 10%, NÃO É significativa |
| Tirando sentido_sozinho (V8): | 2.8897 | 0.0891 | Ao nível de 10%, É significativa |
| Tirando amigos_drogas (V9): V4+V5+V6+V8 | 29.4385 | 0 | Ao nível de 10%, É significativa |
- Ao nível de significância de \(0.1\), a retirada das variáveis tempo_livre_responsaveis (V4), cigarro (V6), sentido_sozinho (V8) e amigos_drogas (V9) impacta o modelo, portanto, estas covariáveis são mantidas no modelo completo.
1.4.3 Passo 3
Passo 3: Ajustar um modelo com as covariáveis retidas no passo 2 (ao nível de significância de \(0.1\)). Retornar as covariáveis excluídas no passo 2 para verificar se não são significativas.
Ao nível de significância de \(0.1\), as covariáveis retidas no passo 2 tempo_livre_responsaveis (V4), cigarro (V6), sentido_sozinho (V8) e amigos_drogas (V9), portanto, estas covariáveis compõem o modelo completo.
As covariáveis excluídas no passo 2: bebida_alcoolica (V5) e responsavel_fuma(V7).
m.comp1 <- flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9,dist='gengamma',data=Drogas_SP_reg) #modelo com variáveis retidas no passo 2
m.compV5 <- flexsurvreg(formula=Surv(tempo,censura)~V4+V5+V6+V8+V9,dist='gengamma',data=Drogas_SP_reg) #voltar com V5
m.compV7 <- flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V7+V8+V9,dist='gengamma',data=Drogas_SP_reg) #voltar com V7
TRV.V5=2*(m.compV5$loglik-m.comp1$loglik) #TRV comparando o modelo sem V5 e o modelo voltando com V5 (H0:beta_V5=0)
TRV.V7=2*(m.compV7$loglik-m.comp1$loglik) #TRV comparando o modelo sem V7 e o modelo voltando com V7 (H0:beta_V7=0)
selecao3 <- data.frame(Modelo=c("Completo: V4+V6+V8+V9","Voltando com bebida_alcoolica (V5): V4+V5+V6+V8+V9","Voltando com responsavel_fuma (V7): V4+V6+V7+V8+V9"),
TRV = c("-",round(c(TRV.V5,TRV.V7),4)),
pvalor = c("-",round(c(1-pchisq(TRV.V5,df=1),1-pchisq(TRV.V7,df=1)),4)),
Significancia = c("-",ifelse(round(c(1-pchisq(TRV.V5,df=1), 1-pchisq(TRV.V7,df=1)),4)
<=0.1,"Ao nível de 10%, É significativa","Ao nível de 10%, NÃO É significativa")))
knitr::kable(selecao3)| Modelo | TRV | pvalor | Significancia |
|---|---|---|---|
| Completo: V4+V6+V8+V9 | - | - | - |
| Voltando com bebida_alcoolica (V5): V4+V5+V6+V8+V9 | 0.9931 | 0.319 | Ao nível de 10%, NÃO É significativa |
| Voltando com responsavel_fuma (V7): V4+V6+V7+V8+V9 | 1.4596 | 0.227 | Ao nível de 10%, NÃO É significativa |
- Ao nível de \(0.1\), nenhuma covariável excluída no passo 2 retorna ao modelo.
1.4.4 Passo 4
- Passo 4: Retornar as covariáveis excluídas no passo 1 (V1,V2,V3) para verificar se não são significativas.
m.comp2<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9,dist='gengamma',data=Drogas_SP_reg) #modelo obtido
m.comp2.V1<-flexsurvreg(formula=Surv(tempo,censura)~V1+V4+V6+V8+V9,dist='gengamma',data=Drogas_SP_reg)
m.comp2.V2<-flexsurvreg(formula=Surv(tempo,censura)~V2+V4+V6+V8+V9,dist='gengamma',data=Drogas_SP_reg)
m.comp2.V3<-flexsurvreg(formula=Surv(tempo,censura)~V3+V4+V6+V8+V9,dist='gengamma',data=Drogas_SP_reg)
TRV.V1=2*(m.comp2.V1$loglik-m.comp2$loglik) #TRV comparando o modelo sem V1 e o modelo voltando com V1 (H0:beta_V1=0)
TRV.V2=2*(m.comp2.V2$loglik-m.comp2$loglik) #TRV comparando o modelo sem V2 e o modelo voltando com V1 (H0:beta_V2=0)
TRV.V3=2*(m.comp2.V3$loglik-m.comp2$loglik) #TRV comparando o modelo sem V3 e o modelo voltando com V3 (H0:beta_V3=0)
selecao4 <- data.frame(Modelo=c("Completo: V4+V6+V8+V9","Voltando com sexo (V1): V1+V4+V6+V8+V9 ","Voltando com cor_raca (V2): V2+V4+V6+V8+V9 ","Voltando com mora_mae (V3): V3+V4+V6+V8+V9"),
TRV = c("-",round(c(TRV.V1,TRV.V2,TRV.V3),4)),
pvalor = c("-",round(c(1-pchisq(TRV.V1,df=1),1-pchisq(TRV.V2,df=1),1-pchisq(TRV.V3,df=1)),4)),
Significancia = c("-",ifelse(round(c(1-pchisq(TRV.V1,df=1), 1-pchisq(TRV.V2,df=1), 1-pchisq(TRV.V3,df=1)),4)
<=0.1,"Ao nível de 10%, É significativa","Ao nível de 10%, NÃO É significativa")))
knitr::kable(selecao4)| Modelo | TRV | pvalor | Significancia |
|---|---|---|---|
| Completo: V4+V6+V8+V9 | - | - | - |
| Voltando com sexo (V1): V1+V4+V6+V8+V9 | 0.0111 | 0.9159 | Ao nível de 10%, NÃO É significativa |
| Voltando com cor_raca (V2): V2+V4+V6+V8+V9 | 1.3141 | 0.2517 | Ao nível de 10%, NÃO É significativa |
| Voltando com mora_mae (V3): V3+V4+V6+V8+V9 | 0.413 | 0.5204 | Ao nível de 10%, NÃO É significativa |
- Ao nível de \(0.1\), nenhuma covariável excluída no passo 1 retorna ao modelo.
1.4.5 Passo 5
- Passo 5: Testar se alguma covariável pode sair do modelo.
m.comp2<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9,dist='gengamma',data=Drogas_SP_reg) #modelo obtido
m.semV4<-flexsurvreg(formula=Surv(tempo,censura)~V6+V8+V9,dist='gengamma',data=Drogas_SP_reg) #modelo sem V4
m.semV6<-flexsurvreg(formula=Surv(tempo,censura)~V4+V8+V9,dist='gengamma',data=Drogas_SP_reg) #modelo sem V6
m.semV8<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V9,dist='gengamma',data=Drogas_SP_reg) #modelo sem V8
m.semV9<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8,dist='gengamma',data=Drogas_SP_reg) #modelo sem V9
TRV.semV4=2*(m.comp2$loglik-m.semV4$loglik) #TRV comparando o modelo completo e o sem V4 (H0:beta_V4=0)
TRV.semV6=2*(m.comp2$loglik-m.semV6$loglik) #TRV comparando o modelo completo e o sem V6 (H0:beta_V6=0)
TRV.semV8=2*(m.comp2$loglik-m.semV8$loglik) #TRV comparando o modelo completo e o sem V8 (H0:beta_V8=0)
TRV.semV9=2*(m.comp2$loglik-m.semV9$loglik) #TRV comparando o modelo completo e o sem V9 (H0:beta_V9=0)
selecao5 <- data.frame(Modelo=c("Completo: V4+V6+V8+V9","Tirando tempo_livre_responsaveis (V4): V6+V8+V9 ","Tirando cigarro (V6): V4+V8+V9","Tirando sentido_sozinho (V8): V4+V6+V9","Tirando amigos_drogas (V9): V4+V6+V8"),
TRV = c("-",round(c(TRV.semV4,TRV.semV6,TRV.semV8,TRV.semV9),4)),
pvalor = c("-",round(c(1-pchisq(TRV.semV4,df=1),1-pchisq(TRV.semV6,df=1),1-pchisq(TRV.semV8,df=1),1-pchisq(TRV.semV9,df=1)),4)),
Significancia = c("-",ifelse(round(c(1-pchisq(TRV.semV4,df=1), 1-pchisq(TRV.semV6,df=1),1-pchisq(TRV.semV8,df=1),1-pchisq(TRV.semV9,df=1)),4)
<=0.1,"Ao nível de 10%, É significativa","Ao nível de 10%, NÃO É significativa")))
knitr::kable(selecao5)| Modelo | TRV | pvalor | Significancia |
|---|---|---|---|
| Completo: V4+V6+V8+V9 | - | - | - |
| Tirando tempo_livre_responsaveis (V4): V6+V8+V9 | 6.7723 | 0.0093 | Ao nível de 10%, É significativa |
| Tirando cigarro (V6): V4+V8+V9 | 17.9297 | 0 | Ao nível de 10%, É significativa |
| Tirando sentido_sozinho (V8): V4+V6+V9 | 3.5034 | 0.0612 | Ao nível de 10%, É significativa |
| Tirando amigos_drogas (V9): V4+V6+V8 | 33.676 | 0 | Ao nível de 10%, É significativa |
- Ao nível de significância de \(0.1\), conclui-se que a retirada de qualquer covariável (V4, V6, V8 ou V9) é significativa, logo, o modelo completo mantém as covariáveis tempo_livre_responsaveis (V4), cigarro (V6), sentido_sozinho (V8) e amigos_drogas (V9).
1.4.6 Passo 6
- Passo 6: Verificar a possibilidade de inclusão de interação dupla entre as variáveis “sobreviventes”.
m.comp2<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9,dist='gengamma',data=Drogas_SP_reg) #modelo obtido
m.V4V6<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V6,dist='gengamma',data=Drogas_SP_reg) #modelo com interação entre V4 e V6
m.V4V8<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8,dist='gengamma',data=Drogas_SP_reg) #modelo com interação entre V4 e V8
m.V4V9<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9,dist='gengamma',data=Drogas_SP_reg) #modelo com interação entre V4 e V9
m.V6V8<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V8,dist='gengamma',data=Drogas_SP_reg) #modelo com interação entre V6 e V8
m.V6V9<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V9,dist='gengamma',data=Drogas_SP_reg) #modelo com interação entre V6 e V9
m.V8V9<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V8*V9,dist='gengamma',data=Drogas_SP_reg) #modelo com interação entre V8 e V9
TRV.V4V6=2*(m.V4V6$loglik-m.comp2$loglik) #TRV comparando o modelo sem interacao e o modelo com V4*V6 (H0:beta_V4*V6=0)
TRV.V4V8=2*(m.V4V8$loglik-m.comp2$loglik) #TRV comparando o modelo sem interacao e o modelo com V4*V8 (H0:beta_V4*V8=0)
TRV.V4V9=2*(m.V4V9$loglik-m.comp2$loglik) #TRV comparando o modelo sem interacao e o modelo com V4*V9 (H0:beta_V4*V9=0)
TRV.V6V8=2*(m.V6V8$loglik-m.comp2$loglik) #TRV comparando o modelo sem interacao e o modelo com V6*V8 (H0:beta_V6*V8=0)
TRV.V6V9=2*(m.V6V9$loglik-m.comp2$loglik) #TRV comparando o modelo sem interacao e o modelo com V6*V9 (H0:beta_V6*V9=0)
TRV.V8V9= 2*(m.V8V9$loglik-m.comp2$loglik) #TRV comparando o modelo sem interacao e o modelo com V8*V9 (H0:beta_V8*V9=0)
selecao6 <- data.frame(Modelo=c("Completo: V4+V6+V8+V9",
"Com interação entre tempo_livre_responsaveis e cigarro (V4 e V6): V4+V6+V8+V9+V4*V6 ",
"Com interação entre tempo_livre_responsaveis e sentido_sozinho (V4 e V8): V4+V6+V8+V9+V4*V8 ",
"Com interação entre tempo_livre_responsaveis e amigos_drogas (V4 e V9): V4+V6+V8+V9+V4*V9 ",
"Com interação entre cigarro e sentido_sozinho (V6 e V8): V4+V6+V8+V9+V6*V8",
"Com interação entre cigarro e amigos_drogas (V6 e V9): V4+V6+V8+V9+V6*V9",
"Com interação entre sentido_sozinho e amigos_drogas (V8 e V9): V4+V6+V8+V9+V8*V9"),
TRV = c("-",round(c(TRV.V4V6,TRV.V4V8,TRV.V4V9,TRV.V6V8,TRV.V6V9,TRV.V8V9),4)),
pvalor = c("-",round(c(1-pchisq(TRV.V4V6,df=1),1-pchisq(TRV.V4V8,df=1),1-pchisq(TRV.V4V9,df=1),1-pchisq(TRV.V6V8,df=1),1-pchisq(TRV.V6V9,df=1),1-pchisq(TRV.V8V9,df=1)),4)),
Significancia = c("-",ifelse(round(c(1-pchisq(TRV.V4V6,df=1),1-pchisq(TRV.V4V8,df=1),1-pchisq(TRV.V4V9,df=1),1-pchisq(TRV.V6V8,df=1),1-pchisq(TRV.V6V9,df=1),1-pchisq(TRV.V8V9,df=1)),4)
<=0.1,"Ao nível de 10%, É significativa","Ao nível de 10%, NÃO É significativa")))
knitr::kable(selecao6)| Modelo | TRV | pvalor | Significancia |
|---|---|---|---|
| Completo: V4+V6+V8+V9 | - | - | - |
| Com interação entre tempo_livre_responsaveis e cigarro (V4 e V6): V4+V6+V8+V9+V4*V6 | 2.4047 | 0.121 | Ao nível de 10%, NÃO É significativa |
| Com interação entre tempo_livre_responsaveis e sentido_sozinho (V4 e V8): V4+V6+V8+V9+V4*V8 | 26.8446 | 0 | Ao nível de 10%, É significativa |
| Com interação entre tempo_livre_responsaveis e amigos_drogas (V4 e V9): V4+V6+V8+V9+V4*V9 | 7.8189 | 0.0052 | Ao nível de 10%, É significativa |
| Com interação entre cigarro e sentido_sozinho (V6 e V8): V4+V6+V8+V9+V6*V8 | 4.1242 | 0.0423 | Ao nível de 10%, É significativa |
| Com interação entre cigarro e amigos_drogas (V6 e V9): V4+V6+V8+V9+V6*V9 | 0.2517 | 0.6159 | Ao nível de 10%, NÃO É significativa |
| Com interação entre sentido_sozinho e amigos_drogas (V8 e V9): V4+V6+V8+V9+V8*V9 | 22.3374 | 0 | Ao nível de 10%, É significativa |
Ao nível de significância de \(0.1\), as interações:
tempo_livre_responsaveis e sentido_sozinho (V4 e V8)
tempo_livre_responsaveis e amigos_drogas (V4 e V9)
cigarro e sentido_sozinho (V6 e V8)
sentido_sozinho e amigos_drogas (V8 e V9)
são significativas, portanto, serão adicionadas ao modelo.
1.4.7 Passo 7
Passo 7: Modelo final escolhido.
\(Modelo \, Final = V4+V6+V8+V9+V4*V8+V4*V9+V6*V8+V8*V9\)
O modelo final ficou composto pelas covariáveis: tempo_livre_responsaveis (V4), cigarro (V6), sentido_sozinho (V8) e amigos_drogas (V9).
Os termos de interação \(V4*V8\), \(V4*V9\), \(V6 e V8\) e \(V6*V9\) foram significativos ao nível de \(0.1\), portanto, foram incluídos no modelo.
modeloGama<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9+V6*V8+V8*V9,dist='gengamma',data=Drogas_SP_reg) #modelo final1.5 Investigando a utilização de modelos mais simples
Uma vez escolhido o conjunto de covariáveis, o interesse se concentra em investigar a utilização de modelos mais simples, casos especiais da gama generalizada, mas não menos adequados aos dados.
Ao rodar a função flexsurvreg para os modelos Exponencial, Weibull e Lognormal utilizando o modelo \(V4+V6+V8+V9+V4*V8+V4*V9+V6*V8+V8*V9\), obtivemos erros de convergência e singularidade:
Warning: Optimisation has probably not converged to the maximum likelihood - Hessian is not positive definite.
Error in solve.default(hessian, tol = tol.solve) : Lapack routine dgesv: system is exactly singular: U[41,41] = 0
A fim de possibilitar a continuidade do trabalho, adotamos a seguinte metodologia: como os últimos termos a serem inseridos no modelo final foram os de interação, removemos cada termo de interação individualmente e ajustamos os modelos Exponencial, Weibull e Lognormal, verificando a presença de erros.
Ao remover cada termo de interação, individualmente, obtivemos os seguintes resultados:
#------------------------------ Modelo final com todas as interações
#Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9+V6*V8+V8*V9,dist='exponential',data=Drogas_SP_reg)
#Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9+V6*V8+V8*V9,dist='weibull',data=Drogas_SP_reg)
#Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9+V6*V8+V8*V9,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
#------------------------------- Modelo final sem interação V8*V9
#Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9+V6*V8,dist='exponential',data=Drogas_SP_reg)
#Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9+V6*V8,dist='weibull',data=Drogas_SP_reg)
#Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9+V6*V8,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
#------------------------------- Modelo final sem interação V6*V8
#Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9+V8*V9,dist='exponential',data=Drogas_SP_reg)
#Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9+V8*V9,dist='weibull',data=Drogas_SP_reg)
#Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9+V8*V9,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
#------------------------------- Modelo final sem interação V4*V9
#Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V6*V8+V8*V9,dist='exponential',data=Drogas_SP_reg)
#Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V6*V8+V8*V9,dist='weibull',data=Drogas_SP_reg)
#Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V6*V8+V8*V9,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
#------------------------------- Modelo final sem interação V4*V8
#Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9+V6*V8+V8*V9,dist='exponential',data=Drogas_SP_reg)
#Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9+V6*V8+V8*V9,dist='weibull',data=Drogas_SP_reg)
#Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9+V6*V8+V8*V9,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
selecao_modelo1 <- data.frame(Modelo=c("Retirando do modelo final a interação $V8*V9$",
"Retirando do modelo final a interação $V6*V8$",
"Retirando do modelo final a interação $V4*V9$",
"Retirando do modelo final a interação $V4*V8$"),
Resultado=c("Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal",
"Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal",
"Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal",
"Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal"))
knitr::kable(selecao_modelo1)| Modelo | Resultado |
|---|---|
| Retirando do modelo final a interação \(V8*V9\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final a interação \(V6*V8\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final a interação \(V4*V9\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final a interação \(V4*V8\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
- Dado que a remoção dos termos de interação, individualmente, não acarretou em melhoras na execução dos modelos Exponencial, Weibull e Lognormal, o próximo passo consistiu na remoção conjunta de termos de interação. Para isso, as combinações dos termos de interação tomados 2 a 2, 3 a 3 e 4 a 4 foram retiradas do modelo e foi verificada a presença de erros.
# combn(4,2)
# combn(4,3)
# -------------------------- Modelo final sem V8*V9 e V6*V8
# Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V4*V9,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
# -------------------------- Modelo final sem V8*V9 e V4*V9
# Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V6*V8,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V6*V8,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V6*V8,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
# -------------------------- Modelo final sem V8*V9 e V4*V8
# Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9+V6*V8,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9+V6*V8,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9+V6*V8,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
#------------------------------ Modelo final sem V6*V8 e V4*V9
# Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V8*V9,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V8*V9,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8+V8*V9,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
#------------------------------ Modelo final sem V6*V8 e V4*V8
# Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9+V8*V9,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9+V8*V9,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9+V8*V9,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
#------------------------------ Modelo final sem V4*V8 e V4*V9
# Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V8+V8*V9,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V8+V8*V9,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V8+V8*V9,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
#------------------------------ Modelo final sem V8*V9, V6*V8 e V4*V9
# Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V8,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
# --------------------------- Modelo final sem V8*V9, V6*V8 e V4*V8
# Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V4*V9,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
#---------------------------- Modelo final sem V8*V9, V4*V9 e V4*V8
# Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V8,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V8,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V8,dist='lnorm',data=Drogas_SP_reg)
# SEM ERROS
#---------------------------- Modelo final sem V6*V8, V4*V8 e V4*V9
# Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V8*V9,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V8*V9,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V8*V9,dist='lnorm',data=Drogas_SP_reg)
# Erro de convergência e singularidade da Matriz Hessiana
#---------------------------- Modelo final sem V8*V9, V6*V8, V4*V9 e V4*V8
# Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9,dist='lnorm',data=Drogas_SP_reg)
# SEM ERROS
selecao_modelo2 <- data.frame(Modelo=c("Retirando do modelo final as interações $V8*V9$ e $V6*V8$",
"Retirando do modelo final as interações $V8*V9$ e $V4*V9$",
"Retirando do modelo final as interações $V8*V9$ e $V4*V8$",
"Retirando do modelo final as interações $V6*V8$ e $V4*V9$",
"Retirando do modelo final as interações $V6*V8$ e $V4*V9$",
"Retirando do modelo final as interações $V4*V9$ e $V4*V8$",
"Retirando do modelo final as interações $V8*V9$, $V6*V8$ e $V4*V9$",
"Retirando do modelo final as interações $V8*V9$, $V6*V8$ e $V4*V8$",
"Retirando do modelo final as interações $V8*V9$, $V4*V8$ e $V4*V9$",
"Retirando do modelo final as interações $V6*V8$, $V4*V8$ e $V4*V9$",
"Retirando do modelo final as interações $V8*V9$, $V6*V8$, $V4*V9$ e $V4*V8$"),
Resultado=c("Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal",
"Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal",
"Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal",
"Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal",
"Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal",
"Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal",
"Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal",
"Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal",
"Não apresentou erros para os modelos Exponencial, Weibull e Lognormal",
"Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal",
"Não apresentou erros para os modelos Exponencial, Weibull e Lognormal"))
knitr::kable(selecao_modelo2)| Modelo | Resultado |
|---|---|
| Retirando do modelo final as interações \(V8*V9\) e \(V6*V8\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final as interações \(V8*V9\) e \(V4*V9\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final as interações \(V8*V9\) e \(V4*V8\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final as interações \(V6*V8\) e \(V4*V9\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final as interações \(V6*V8\) e \(V4*V9\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final as interações \(V4*V9\) e \(V4*V8\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final as interações \(V8*V9\), \(V6*V8\) e \(V4*V9\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final as interações \(V8*V9\), \(V6*V8\) e \(V4*V8\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final as interações \(V8*V9\), \(V4*V8\) e \(V4*V9\) | Não apresentou erros para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final as interações \(V6*V8\), \(V4*V8\) e \(V4*V9\) | Apresentou erros de convergência e singularidade para os modelos Exponencial, Weibull e Lognormal |
| Retirando do modelo final as interações \(V8*V9\), \(V6*V8\), \(V4*V9\) e \(V4*V8\) | Não apresentou erros para os modelos Exponencial, Weibull e Lognormal |
- De acordo com os resultados obtidos nesses testes, temos que o modelo apenas com o termo de interação \(V6*V8\) e o modelo sem interações que não apresentaram erros. Portanto, optamos por adotar o \(Modelo \, Final = V4+V6+V8+V9+V6*V8\), composto pelas covariáveis tempo_livre_responsaveis (V4), cigarro (V6), sentido_sozinho (V8) e amigos_drogas (V9) e pelo termo de interação \(V6 e V8\).
#---------------------------- Modelo final sem V8*V9, V4*V9 e V4*V8
# modelo final Gama Generalizado
modeloGama<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V8,dist='gengamma',data=Drogas_SP_reg)
# Modelo exponencial
modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V8,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V8,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V8,dist='lnorm',data=Drogas_SP_reg)
TRV.exp=2*(modeloGama$loglik-modeloexp$loglik) #TRV comparando o modelo exponencial e o modelo Gama (H0: O modelo exponencial é adequado)
TRV.Wei=2*(modeloGama$loglik-modeloWei$loglik)#TRV comparando o modelo Weibull e o modelo Gama (H0: O modelo Weibull é adequado)
TRV.Logn=2*(modeloGama$loglik-modeloLogn$loglik)#TRV comparando o modelo Weibull e o modelo Gama (H0: O modelo Weibull é adequado)
selecao_modelo3 <- data.frame(Modelos = c("Exponencial","Weibull","Lognormal"),
TRV = round(c(TRV.exp,TRV.Wei,TRV.Logn),4),
pvalor = round(c(1-pchisq(TRV.exp,df=2),1-pchisq(TRV.Wei,df=1),1-pchisq(TRV.Logn,df=1)),4),
Significancia = ifelse (round(c(1-pchisq(TRV.exp,df=2),1-pchisq(TRV.Wei,df=1),1-pchisq(TRV.Logn,df=1)),4)<=0.1,"Ao nível de 10%, NÃO É adequado","Ao nível de 10%, É adequado"))
knitr::kable(selecao_modelo3)| Modelos | TRV | pvalor | Significancia |
|---|---|---|---|
| Exponencial | 176.8423 | 0.0000 | Ao nível de 10%, NÃO É adequado |
| Weibull | 1.2317 | 0.2671 | Ao nível de 10%, É adequado |
| Lognormal | 0.3215 | 0.5707 | Ao nível de 10%, É adequado |
#---------------------------- Modelo final sem interações
# modelo final Gama Generalizado
# modeloGama<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9,dist='gengamma',data=Drogas_SP_reg)
# Modelo exponencial
# modeloexp<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9,dist='exponential',data=Drogas_SP_reg)
# Modelo Weibull
# modeloWei<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9,dist='weibull',data=Drogas_SP_reg)
# Modelo Log-normal
# modeloLogn<-flexsurvreg(formula=Surv(tempo,censura)~V4+V6+V8+V9,dist='lnorm',data=Drogas_SP_reg)
# TRV.exp=2*(modeloGama$loglik-modeloexp$loglik) #TRV comparando o modelo exponencial e o modelo Gama (H0: O modelo exponencial é adequado)
# TRV.Wei=2*(modeloGama$loglik-modeloWei$loglik)#TRV comparando o modelo Weibull e o modelo Gama (H0: O modelo Weibull é adequado)
# TRV.Logn=2*(modeloGama$loglik-modeloLogn$loglik)#TRV comparando o modelo Weibull e o modelo Gama (H0: O modelo Weibull é adequado)
#
# selecao_modelo3 <- data.frame(Modelos = c("Exponencial","Weibull","Lognormal"),
# TRV = round(c(TRV.exp,TRV.Wei,TRV.Logn),4),
# pvalor = round(c(1-pchisq(TRV.exp,df=2),1-pchisq(TRV.Wei,df=1),1-pchisq(TRV.Logn,df=1)),4),
# Significancia = ifelse (round(c(1-pchisq(TRV.exp,df=2),1-pchisq(TRV.Wei,df=1),1-pchisq(TRV.Logn,df=1)),4)<=0.1,"Ao nível de 10%, NÃO É adequado","Ao nível de 10%, É adequado"))
#
# knitr::kable(selecao_modelo3)Ao nível de significância de \(0.1\), rejeita-se a hipótese nula de adequação do modelo Exponencial. Os modelos Weibull e Lognormal se adequam aos dados.
Escolhemos então o modelo Lognormal para modelar os tempos até a experimentação de drogas ilícitas pois foi o que apresentou maior p-valor (\(0.5707\).
Desse modo, todas as análises seguintes são baseadas nesse modelo.
1.6 Estimando os parâmetros do modelo Lognormal
modeloLogn<-survreg(formula=Surv(tempo,censura)~V4+V6+V8+V9+V6*V8,dist='lognormal',data=Drogas_SP_reg)
summary(modeloLogn)##
## Call:
## survreg(formula = Surv(tempo, censura) ~ V4 + V6 + V8 + V9 +
## V6 * V8, data = Drogas_SP_reg, dist = "lognormal")
## Value Std. Error z p
## (Intercept) 3.24792 284.92012 0.01 0.991
## V42 -0.03061 0.04815 -0.64 0.525
## V43 -0.04666 0.04781 -0.98 0.329
## V44 -0.01312 0.04012 -0.33 0.744
## V45 0.03221 0.04634 0.70 0.487
## V62 0.12562 0.05713 2.20 0.028
## V82 0.05505 0.04564 1.21 0.228
## V83 -0.00849 0.04225 -0.20 0.841
## V84 0.03972 0.04402 0.90 0.367
## V85 -0.04165 0.05029 -0.83 0.408
## V92 -0.55839 284.92012 0.00 0.998
## V93 -0.58017 284.92012 0.00 0.998
## V94 -0.65450 284.92012 0.00 0.998
## V96 0.00487 527.77047 0.00 1.000
## V62:V82 -0.08653 0.07070 -1.22 0.221
## V62:V83 0.03219 0.07249 0.44 0.657
## V62:V84 -0.04878 0.07804 -0.63 0.532
## V62:V85 -0.02531 0.08123 -0.31 0.755
## Log(scale) -2.35354 0.11434 -20.58 <2e-16
##
## Scale= 0.095
##
## Log Normal distribution
## Loglik(model)= -103.1 Loglik(intercept only)= -164.8
## Chisq= 123.44 on 17 degrees of freedom, p= 3.4e-18
## Number of Newton-Raphson Iterations: 20
## n= 445
- Ao nível de significância de \(0.1\), \(V9\) é retirada do modelo \(V4+V6+V8+V9+V6*V8\) por não ser significativa (ou seja, não se rejeita a hipótese nula de que o coeficiente é igual a zero)
modeloLogn<-survreg(formula=Surv(tempo,censura)~V6+V8+V6*V8,dist='lognormal',data=Drogas_SP_reg)
summary(modeloLogn)##
## Call:
## survreg(formula = Surv(tempo, censura) ~ V6 + V8 + V6 * V8, data = Drogas_SP_reg,
## dist = "lognormal")
## Value Std. Error z p
## (Intercept) 2.64948 0.03285 80.65 <2e-16
## V62 0.22528 0.05612 4.01 6e-05
## V82 0.07200 0.04702 1.53 0.126
## V83 0.00652 0.04285 0.15 0.879
## V84 0.05539 0.04498 1.23 0.218
## V85 -0.04956 0.05046 -0.98 0.326
## V62:V82 -0.12910 0.06683 -1.93 0.053
## V62:V83 -0.00111 0.06965 -0.02 0.987
## V62:V84 -0.09851 0.07398 -1.33 0.183
## V62:V85 -0.06128 0.07811 -0.78 0.433
## Log(scale) -2.24840 0.11509 -19.54 <2e-16
##
## Scale= 0.106
##
## Log Normal distribution
## Loglik(model)= -122.4 Loglik(intercept only)= -164.8
## Chisq= 84.79 on 9 degrees of freedom, p= 1.8e-14
## Number of Newton-Raphson Iterations: 6
## n= 445
- Ao nível de significância de \(0.1\), \(V4\) é retirada do modelo \(V4+V6+V8+V6*V8\) por não ser significativa (ou seja, não se rejeita a hipótese nula de que o coeficiente é igual a zero).
modeloLogn<-survreg(formula=Surv(tempo,censura)~V6+V8+V6*V8,dist='lognormal',data=Drogas_SP_reg)
summary(modeloLogn)##
## Call:
## survreg(formula = Surv(tempo, censura) ~ V6 + V8 + V6 * V8, data = Drogas_SP_reg,
## dist = "lognormal")
## Value Std. Error z p
## (Intercept) 2.64948 0.03285 80.65 <2e-16
## V62 0.22528 0.05612 4.01 6e-05
## V82 0.07200 0.04702 1.53 0.126
## V83 0.00652 0.04285 0.15 0.879
## V84 0.05539 0.04498 1.23 0.218
## V85 -0.04956 0.05046 -0.98 0.326
## V62:V82 -0.12910 0.06683 -1.93 0.053
## V62:V83 -0.00111 0.06965 -0.02 0.987
## V62:V84 -0.09851 0.07398 -1.33 0.183
## V62:V85 -0.06128 0.07811 -0.78 0.433
## Log(scale) -2.24840 0.11509 -19.54 <2e-16
##
## Scale= 0.106
##
## Log Normal distribution
## Loglik(model)= -122.4 Loglik(intercept only)= -164.8
## Chisq= 84.79 on 9 degrees of freedom, p= 1.8e-14
## Number of Newton-Raphson Iterations: 6
## n= 445
Ao nível de significância de \(0.1\), a interação \(V62:V82\) é significativa. De acordo com a literatura, recomenda-se manter as variáveis utilizadas em uma interação significativa ainda que elas sozinhas não sejam significativas. Portanto, utilizaremos o modelo \(V6+V8+V6*V8\).
Observações: Dentre os modelos considerados nessa etapa, o modelo com menor Bayesian Information Criterion (BIC) foi o modelo \(V6+V8+V6*V8\), o que corrobora sua escolha. Além disso, ao repetir o teste TRV para seleção de modelos mais simples utilizando como modelo final \(V6+V8+V6*V8\), a distribuição Lognormal continuou apresentando o maior p-valor e, portanto, sendo escolhida como a adequada para o processo de modelagem.
1.7 Adequação do modelo ajustado
1.7.1 Resíduos de Cox-Snell
- Os resíduos \(\hat{e_i}\) vêm de uma população homogênea e devem seguir uma distribuição exponencial padrão se o modelo for adequado.
xb<-modeloLogn$linear.predictors
sigma<-modeloLogn$scale
res<-(log(Drogas_SP_reg$tempo)-(xb))/sigma
ei<- -log(1-pnorm(res)) #residuos de cox-Snell = funcao de taxa de falha acumulada
# Para o modelo de regressão lognormal:
ekm1<-survfit(Surv(ei,Drogas_SP_reg$censura)~1)
t<-ekm1$time
st<-ekm1$surv
sexp<-exp(-t)
par(mfrow=c(1,2))
plot(st,sexp,xlab="S(ei): Kaplan-Meier",ylab="S(ei): Exponencial padrao",pch=16)
plot(ekm1,conf.int=F,mark.time=F, xlab="Residuos de Cox-Snell", ylab="Sobrevivencia estimada")
lines(t,sexp,lty=4)
legend("bottomleft",lty=c(1,4),c("Kaplan-Meier","Exponencial padrao"),cex=0.8,bty="n")
- Acredita-se que o modelo de regressão Lognormal se encontra
razoavelmente ajustado aos dados.
1.7.2 Resíduos padronizados
- Se o modelo de regressão lognormal for apropriado, \(\hat{\upsilon_i}\) é normal padrão.
xb<-modeloLogn$linear.predictors
sigma<-modeloLogn$scale
nu<-(log(Drogas_SP_reg$tempo)-(xb))/sigma
resid<-exp(nu)
ekm<- survfit(Surv(resid,Drogas_SP_reg$censura)~1)
resid<-ekm$time
sln<-pnorm(-log(resid)) #funcao de sobrevivencia da lognormal padrao
par(mfrow=c(1,2))
plot(ekm$surv,sln, xlab="S(ei): Kaplan-Meier",ylab="S(ei): Lognormal padrão",pch=16)
plot(ekm,conf.int=F,mark.time=F,xlab="Residuos (ei)",ylab="Sobrevivencia estimada",pch=16)
lines(resid,sln,lty=2)
legend("bottomleft",lty=c(1,2),c("Kaplan-Meier","Log-normal padrao"),cex=0.8,bty="n")1.7.3 Resíduos deviance
- Se o modelo de regressão lognormal for apropriado, os gráficos dos resíduos deviance vs. tempo devem apresentar comportamento aleatório em torno de zero.
res<-(log(Drogas_SP_reg$tempo)-(xb))/sigma
ei<- -log(1-pnorm(res)) #residuos de cox-Snell
mi<-Drogas_SP_reg$censura-ei
di<-sign(mi)*(-2*(mi+Drogas_SP_reg$censura*log(Drogas_SP_reg$censura-mi)))^0.5
par(mfrow=c(1,2))
plot(Drogas_SP$cigarro,di,xlab="cigarro (V6)")
plot(Drogas_SP$sentido_sozinho,di,xlab="sentido_sozinho (V8)")
- A partir da análise dos gráficos acima, nota-se um aglomerado na
categoria 2 de V6, o que não corrobora a hipótese de comportamento
aleatório em torno de zero. Entretanto, na categoria 1 de V6 e nas
demais categorias de V8, podemos observar que os pontos estão
distribuídos de forma mais apropriada.
1.8 Interpretação dos resultados
- Tomando-se o exponencial dos coeficientes estimados, obtém-se a razão dos tempos medianos de sobrevivência.
interpretacao1 <- data.frame(Exp_beta = round(c(exp(modeloLogn$coefficients[2:6])),2))
knitr::kable(interpretacao1)| Exp_beta | |
|---|---|
| V62 | 1.25 |
| V82 | 1.07 |
| V83 | 1.01 |
| V84 | 1.06 |
| V85 | 0.95 |
V62: o tempo mediano até a experimentação de drogas ilícitas para os estudantes que não fumaram cigarro, é 1.25 vezes o tempo mediano dos estudantes que alguma vez na vida já fumaram, mesmo uma ou duas tragadas.
V82: o tempo mediano até a experimentação de drogas ilícitas para os estudantes que nos últimos 12 meses RARAMENTE se sentiram sozinhos é \(1.07\) vezes o tempo mediano dos estudantes que NUNCA se sentiram sozinhos.
V83: o tempo mediano até a experimentação de drogas ilícitas para os estudantes que nos últimos 12 meses ÀS VEZES se sentiram sozinhos é \(1.01\) vezes o tempo mediano dos estudantes que NUNCA se sentiram sozinhos.
V84: o tempo mediano até a experimentação de drogas ilícitas para os estudantes que nos últimos 12 meses NA MAIORIA DAS VEZES se sentiram sozinhos é \(1.06\) vezes o tempo mediano dos que NUNCA se sentiram sozinhos.
V85: o tempo mediano até a experimentação de drogas ilícitas para os estudantes que nos últimos 12 meses SEMPRE se sentiram sozinhos é \(0.95\) vezes o tempo mediano dos que NUNCA se sentiram sozinhos.
A experimentação de cigarro e o fato do estudante sempre se sentir sozinho com maior frequência acarretam em menores tempos medianos até a experimentação de drogas ilícitas.
interpretacao2 <- data.frame(Interacao=c("Fixado sentido_sozinho=5 (Sempre)","Fixado sentido_sozinho=2 (Raramente) "),
Exp_beta = round(c(exp(modeloLogn$coefficients[2]+modeloLogn$coefficients[10]),
exp(modeloLogn$coefficients[2]+modeloLogn$coefficients[7])),2))
knitr::kable(interpretacao2)| Interacao | Exp_beta |
|---|---|
| Fixado sentido_sozinho=5 (Sempre) | 1.18 |
| Fixado sentido_sozinho=2 (Raramente) | 1.10 |
Fixando a variável sentido_sozinho (V8) na categoria 5 (Sempre), não ter fumado alguma vez aumenta em \(18\%\) o tempo mediano da experimentação de drogas ilícitas em comparação a ter fumado alguma vez.
Fixando a variável sentido_sozinho (V8) na categoria 2 (Raramente), não ter fumado alguma vez aumenta em \(10\%\) o tempo mediano da experimentação de drogas ilícitas em comparação a ter fumado alguma vez.
1.9 Curvas de sobrevivencia estimadas pelo modelo final
1.9.1 Comparando V6=1 com V6=2:
- Fixou-se sentido_sozinho (V8) na categoria de referência (1- Nunca se sentiu sozinho nos últimos 12 meses).
b0<-modeloLogn$coefficients[1]
beta.V62<-modeloLogn$coefficients[2]
beta.V82<-modeloLogn$coefficients[3]
beta.V83<-modeloLogn$coefficients[4]
beta.V84<-modeloLogn$coefficients[5]
beta.V85<-modeloLogn$coefficients[6]
beta.V62.V82<-modeloLogn$coefficients[7]
beta.V62.V83<-modeloLogn$coefficients[8]
beta.V62.V84<-modeloLogn$coefficients[9]
beta.V62.V85<-modeloLogn$coefficients[10]
sigma<-modeloLogn$scale
V62<-0 # 1 se V6=2 e 0 c.c
V82<-0 #1, se V8=2, 0 caso contrário.
V83<-0 #1, se V8=3, 0 caso contrário.
V84<-0 #1, se V8=4, 0 caso contrário.
V85<-0 #1, se V8=5, 0 caso contrário.
V62.V82<- 0
V62.V83<- 0
V62.V84<- 0
V62.V85<- 0
tempo<-0:25
sobrev.1<-pnorm((-log(tempo)+b0+beta.V62*V62+
beta.V82*V82+
beta.V83*V83+
beta.V84*V84+
beta.V85*V85+
beta.V62.V82*V62.V82+
beta.V62.V83*V62.V83+
beta.V62.V84*V62.V84+
beta.V62.V85*V62.V85)/sigma)
V62<-1
sobrev.2<-pnorm((-log(tempo)+b0+beta.V62*V62+
beta.V82*V82+
beta.V83*V83+
beta.V84*V84+
beta.V85*V85+
beta.V62.V82*V62.V82+
beta.V62.V83*V62.V83+
beta.V62.V84*V62.V84+
beta.V62.V85*V62.V85)/sigma)
par(mfrow=c(1,1))
plot(sobrev.1,tempo*0,pch=" ",ylim=range(c(0,1)), xlim=range(c(0,25)),
xlab="Tempo (anos)",ylab="S(t) estimada",bty="n",main="Cigarro")
lines(tempo,sobrev.1,lty=1)
lines(tempo,sobrev.2,lty=2)
legend("topright",lty=c(1,2),c("Sim","Não"),lwd=1, bty="n",cex=0.8)A função de sobrevivência nos indica que o tempo até a experimentação de drogas ilícitas é consideravelmente menor para os estudantes que já fumaram cigarro alguma vez. A distância entre as curvas fornece uma dimensão dessa diferença.
O risco de experimentar drogas ilícitas começa a aumentar mais cedo para os estudantes que já fumaram cigarro alguma vez.
1.9.2 Comparando V8=5 com V8=1:
- Fixou-se cigarro (V6) na categoria (2- Não).
V62<-1 # 1 se V6=2 e 0 c.c
V82<-0 #1, se V8=2, 0 caso contrário.
V83<-0 #1, se V8=3, 0 caso contrário.
V84<-0 #1, se V8=4, 0 caso contrário.
V85<-1 #1, se V8=5, 0 caso contrário.
V62.V82<- 0
V62.V83<- 0
V62.V84<- 0
V62.V85<- 1
tempo<-0:25
sobrev.1<-pnorm((-log(tempo)+b0+beta.V62*V62+
beta.V82*V82+
beta.V83*V83+
beta.V84*V84+
beta.V85*V85+
beta.V62.V82*V62.V82+
beta.V62.V83*V62.V83+
beta.V62.V84*V62.V84+
beta.V62.V85*V62.V85)/sigma)
V85<-0
sobrev.2<-pnorm((-log(tempo)+b0+beta.V62*V62+
beta.V82*V82+
beta.V83*V83+
beta.V84*V84+
beta.V85*V85+
beta.V62.V82*V62.V82+
beta.V62.V83*V62.V83+
beta.V62.V84*V62.V84+
beta.V62.V85*V62.V85)/sigma)
par(mfrow=c(1,1))
plot(sobrev.1,tempo*0,pch=" ",ylim=range(c(0,1)), xlim=range(c(0,25)),
xlab="Tempo (anos)",ylab="S(t) estimada",bty="n",main="NOS ÚLTIMOS 12 MESES com que frequência tem se sentido sozinho(a)?")
lines(tempo,sobrev.1,lty=1)
lines(tempo,sobrev.2,lty=2)
legend("topright",lty=c(1,2),c("Sempre","Nunca"),lwd=1, bty="n",cex=0.8)- A função de sobrevivência nos indica que o tempo até a experimentação de drogas ilícitas é menor para os estudantes que nos últimos 12 meses declararam que sempre se sentem sozinhos. Nota-se a proximidade das curvas.
1.10 Estimação para um indivíduo
1.10.1 Perfil 1: indivíduo com valores V6=2 (Não), V8=1 (Nunca)
V62<-1 # 1 se V6=2 e 0 c.c
V82<-0 #1, se V8=2, 0 caso contrário.
V83<-0 #1, se V8=3, 0 caso contrário.
V84<-0 #1, se V8=4, 0 caso contrário.
V85<-0 #1, se V8=5, 0 caso contrário.
V62.V82<- 0
V62.V83<- 0
V62.V84<- 0
V62.V85<- 0
xb<-b0+beta.V62*V62+
beta.V82*V82+
beta.V83*V83+
beta.V84*V84+
beta.V85*V85+
beta.V62.V82*V62.V82+
beta.V62.V83*V62.V83+
beta.V62.V84*V62.V84+
beta.V62.V85*V62.V85
par(mfrow=c(3,1))
# Função de sobrevivência
tempo<-0:30
sobreviv<-pnorm((-log(tempo)+xb)/sigma)
plot(tempo,sobreviv,type="l",ylim=range(c(0,1)), xlim=range(c(0,30)),
xlab="Tempo (anos)",ylab="S(t) estimada",bty="n",main="Função de sobrevivência")
# Função Taxa de falha
fdensidade<-(1/(sqrt(2*pi)*tempo*sigma))*exp(-0.5*((log(tempo)-xb)/sigma)^2)
plot(tempo,fdensidade,type="l",xlab="Tempo (anos)",ylab="f(t) estimada",bty="n",
main="Função densidade")
h_t<-fdensidade/sobreviv
plot(tempo,h_t,type="l",xlab="Tempo (anos)",ylab="h(t) estimada",bty="n",main=
"Função taxa de falha")# Percentis da distribuição do tempo de vida:
t_p<-function(p){
exp(qnorm(p)*sigma+xb)
}
knitr::kable(data.frame(Percentis=c("Percentil 20%","Mediana","Percentil 90%"),Valor=c(t_p(0.2),t_p(0.5),t_p(0.9))))| Percentis | Valor |
|---|---|
| Percentil 20% | 16.21459 |
| Mediana | 17.72117 |
| Percentil 90% | 20.28843 |
Os indivíduos do perfil 1 (Nunca fumaram; Nos últimos 12 meses, nunca se sentem sozinhos), começam a apresentar um decaimento da função de sobrevivência por volta de 15 anos. Aos 20 anos, a função de sobrevivência já se encontra próxima a zero.
A partir dos 15 anos, o risco instantâneo de experimentar drogas ilícitas começa a aumentar.
Aos 20 anos, o risco de experimentar drogas ilícitas é de 90%.
1.10.2 Perfil 2: indivíduo com valores V6=1 (Sim), V8=5 (Sempre)
V62<-0 # 1 se V6=2 e 0 c.c
V82<-0 #1, se V8=2, 0 caso contrário.
V83<-0 #1, se V8=3, 0 caso contrário.
V84<-0 #1, se V8=4, 0 caso contrário.
V85<-1 #1, se V8=5, 0 caso contrário.
V62.V82<- 0
V62.V83<- 0
V62.V84<- 0
V62.V85<- 0
xb<-b0+beta.V62*V62+
beta.V82*V82+
beta.V83*V83+
beta.V84*V84+
beta.V85*V85+
beta.V62.V82*V62.V82+
beta.V62.V83*V62.V83+
beta.V62.V84*V62.V84+
beta.V62.V85*V62.V85
par(mfrow=c(3,1))
# Função de sobrevivência
tempo<-0:30
sobreviv<-pnorm((-log(tempo)+xb)/sigma)
plot(tempo,sobreviv,type="l",ylim=range(c(0,1)), xlim=range(c(0,30)),
xlab="Tempo (anos)",ylab="S(t) estimada",bty="n",main="Função de sobrevivência")
# Função Taxa de falha
fdensidade<-(1/(sqrt(2*pi)*tempo*sigma))*exp(-0.5*((log(tempo)-xb)/sigma)^2)
plot(tempo,fdensidade,type="l",xlab="Tempo (anos)",ylab="f(t) estimada",bty="n",
main="Função densidade")
h_t<-fdensidade/sobreviv
plot(tempo,h_t,type="l",xlab="Tempo (anos)",ylab="h(t) estimada",bty="n",main=
"Função taxa de falha")# Percentis da distribuição do tempo de vida:
t_p<-function(p){
exp(qnorm(p)*sigma+xb)
}
knitr::kable(data.frame(Percentis=c("Percentil 20%","Mediana","Percentil 90%"),Valor=c(t_p(0.2),t_p(0.5),t_p(0.9))))| Percentis | Valor |
|---|---|
| Percentil 20% | 12.31815 |
| Mediana | 13.46269 |
| Percentil 90% | 15.41302 |
Os indivíduos do perfil 2 (Já fumaram alguma vez; Nos últimos 12 meses, sempre se sentem sozinhos), apresentam decaimento da função de sobrevivência por volta de 13 anos, mais cedo que os indivíduos do perfil 1.
O risco instantâneo de experimentar drogas ilícitas começa a apresentar aumento entre 10 e 15 anos, também mais cedo que os indivíduos do perfil 1.
Aos 15 anos, o risco de experimentar drogas ilícitas é de 90%.