Análise de Sobrevivência - Projeto 01
1 Análise descritiva
1.1 Descrição da base
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.
A base de dados contém 445 observações e 11 variáveis.
A variável V0007 foi removida por não conter descrição.
1.2 Descrição das variáveis
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=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=Sim, 2=Não 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=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 |
| VB04006A | responsavel_fuma | Algum de seus pais ou responsáveis fuma? | 1=Sim, 2=Não 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")1.3 Falhas e censuras
falhas_censuras <- data.frame(Categoria = c("Censura","Falha"),Percentual=round(as.matrix((table(Drogas_SP$censura)/445)*100),2))
ggplot(falhas_censuras, aes(x="", y=Percentual, fill=Categoria))+
geom_bar(width = 1, stat = "identity")+
coord_polar("y", start=0)+
scale_fill_brewer(palette="Set1")+
theme_minimal()+
ggtitle("Percentual de falhas e censuras")+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
plot.title = element_text(hjust = 0.5)) +
xlab(element_blank())+
ylab(element_blank())+
geom_text(aes(x =c(1,1.2) ,y = c(50, 5),
label = paste(`Percentual`,"%")), size=5)De acordo com a Figura 1,
- 90.34% dos estudantes da amostra ainda não experimentaram drogas ilícitas.
- 9.66% dos estudantes da amostra já experimentaram drogas ilícitas.
1.4 Variável resposta
A variável resposta ‘tempo’, que corresponde ao tempo até o primeiro uso de drogas no contexto de análise de sobrevivência, nada mais é que a idade dos estudantes no momento que a coleta dos dados pela PeNSE foi realizada.
par(mfrow=c(4,1))
ggplot(data = Drogas_SP, aes (y = tempo))+
geom_boxplot(fill="#E69F00") +
theme_classic() +
scale_fill_brewer(palette="Set1")+
ggtitle("Boxplot do tempo até o primeiro uso de drogas (Idade dos estudantes na PeNSE)\n")+
theme(plot.title = element_text(hjust = 0.5))+
ylab("Idade\n")De acordo com o boxplot, podemos observar que:
- As idades mínima e máxima observadas nos estudantes da amostra são 11 e 16 anos, respectivamente, e ambas são consideradas outliers ;
- De acordo com o primeiro quartil, 25% dos estudantes da amostra possuem idade igual ou inferior a 13 anos;
- De acordo com a mediana, 50% dos estudantes da amostra possuem idade igual ou inferior a 14 anos;
- De acordo com o terceiro quartil, 75% dos estudantes da amostra possuem idade igual ou inferior a 14 anos;
- De acordo com o intervalo interquartil, 50% dos estudantes da amostra possuem idade entre 13 e 14 anos.
1.5 Covariáveis
As demais variáveis se referem às seguintes perguntas feitas na pesquisa:
1.5.1 Qual é o seu sexo?
sexo <- data.frame(Sexo = c("Masculino","Feminino"),Percentual=round(as.matrix((table(Drogas_SP$sexo)/445)*100),2))
ggplot(sexo, aes(x="", y=Percentual, fill=Sexo))+
geom_bar(width = 1, stat = "identity")+
coord_polar("y", start=0)+
scale_fill_brewer(palette="Set1")+
theme_minimal()+
ggtitle("Percentual de obs. por categorias de 'sexo'" )+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
plot.title = element_text(hjust = 0.5)) +
xlab(element_blank())+
ylab(element_blank())+
geom_text(aes(y = `Percentual`/2 + c(0, cumsum(`Percentual`)[-length(`Percentual`)]),
label = paste(`Percentual`,"%")), size=5)1.5.2 Qual é a sua cor ou raça?
cor_raca <- data.frame(Cor=c("Branca","Preta","Amarela","Parda","Indigena"),Percentual = round(as.matrix((table(Drogas_SP$cor_raca)/445)*100),2))
ggplot(cor_raca, aes(x=Cor, y=Percentual,fill=Cor)) +
geom_bar(stat="identity",color="black", width=0.8)+
ylim(0,100)+
geom_text(aes(label= paste0(Percentual,"%")), vjust=-1.5, color="black", size=3.5)+
scale_fill_brewer(palette="Set1")+
xlab("\nCor")+
ylab("Percentual (%)\n")+
ggtitle("Percentual de obs. por categorias de 'cor_raca'\n" )+
theme_classic()+
theme(plot.title = element_text(hjust = 0.5))1.5.3 Você mora com sua mãe?
mora_mae <- data.frame(Categoria = c("Sim","Não"),Percentual=round(as.matrix((table(Drogas_SP$mora_mae)/445)*100),2))
ggplot(mora_mae, aes(x="", y=Percentual, fill=Categoria))+
geom_bar(width = 1, stat = "identity")+
coord_polar("y", start=0)+
scale_fill_brewer(palette="Set1")+
theme_minimal()+
ggtitle("Percentual de obs. por categorias de 'mora_mae'" )+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
plot.title = element_text(hjust = 0.5)) +
xlab(element_blank())+
ylab(element_blank())+
geom_text(aes(x=c(1,1.3),y = c(45.17,100),
label = paste(`Percentual`,"%")), size=5)1.5.4 NOS ÚLTIMOS 30 DIAS, com que frequência seus pais ou responsáveis sabiam realmente o que você estava fazendo em seu tempo livre?
tempo_livre_responsaveis <- data.frame(Categoria=c("Nunca", "Raramente", "Às vezes", "Na maior parte do tempo", "Sempre"),Percentual = round(as.matrix((table(Drogas_SP$tempo_livre_responsaveis)/445)*100),2))
ggplot(tempo_livre_responsaveis, aes(x=Categoria, y=Percentual,fill=Categoria)) +
geom_bar(stat="identity",color="black", width=0.8)+
ylim(0,100)+
geom_text(aes(label= paste0(Percentual,"%")), vjust=-1.5, color="black", size=3.5)+
scale_fill_brewer(palette="Set1")+
xlab("\nCor")+
ylab("Percentual (%)\n")+
ggtitle("Percentual de obs. por categorias de 'tempo_livre_responsaveis'\n" )+
theme_classic()+
theme(plot.title = element_text(hjust = 0.5))1.5.5 Alguma vez na vida você tomou uma dose de bebida alcoólica?
bebida_alcoolica <- data.frame(Categoria = c("Sim","Não"),Percentual = round(as.matrix((table(Drogas_SP$bebida_alcoolica)/445)*100),2))
ggplot(bebida_alcoolica, aes(x="", y=Percentual, fill=Categoria))+
geom_bar(width = 1, stat = "identity")+
coord_polar("y", start=0)+
scale_fill_brewer(palette="Set1")+
theme_minimal()+
ggtitle("Percentual de obs. por categorias de 'bebida_alcoolica'" )+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
plot.title = element_text(hjust = 0.5)) +
xlab(element_blank())+
ylab(element_blank())+
geom_text(aes(x=c(0.9,0.9),y = c(25,75),
label = paste(`Percentual`,"%")), size=5)1.5.6 Alguma vez na vida, você já fumou cigarro, mesmo uma ou duas tragadas?
cigarro <- data.frame(Categoria = c("Sim","Não"),Percentual= round(as.matrix((table(Drogas_SP$cigarro)/445)*100),2))
ggplot(cigarro, aes(x="", y=Percentual, fill=Categoria))+
geom_bar(width = 1, stat = "identity")+
coord_polar("y", start=0)+
scale_fill_brewer(palette="Set1")+
theme_minimal()+
ggtitle("Percentual de obs. por categorias de 'cigarro'" )+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
plot.title = element_text(hjust = 0.5)) +
xlab(element_blank())+
ylab(element_blank())+
geom_text(aes(x=c(1.2,0.9),y = c(10,55),
label = paste(`Percentual`,"%")), size=5)1.5.7 Algum de seus pais ou responsáveis fuma?
responsavel_fuma <- data.frame(Categoria = c("Nenhum deles", "Só meu pai ou responsável do sexo masculino", "Só minha mãe ou responsável do sexo feminino", "Meu pai e minha mãe ou responsáveis", "Não sei", "Não informado"),Percentual = round(as.matrix((table(Drogas_SP$responsavel_fuma)/445)*100),2))
ggplot(responsavel_fuma, aes(x=Categoria, y=Percentual,fill=Categoria)) +
geom_bar(stat="identity",color="black", width=0.8)+
ylim(0,100)+
geom_text(aes(label= paste0(Percentual,"%")), vjust=-1.5, color="black", size=3.5)+
scale_fill_brewer(palette="Set1")+
xlab("\nCor")+
ylab("Percentual (%)\n")+
ggtitle("Percentual de obs. por categorias de 'responsavel_fuma'\n" )+
theme_classic()+
theme(plot.title = element_text(hjust = 0.5),axis.text.x = element_blank())1.5.8 NOS ÚLTIMOS 12 MESES com que frequência tem se sentido sozinho(a)?
sentido_sozinho <- data.frame(Categoria = c("Nunca", "Raramente", "Às vezes", "Na maioria das vezes", "Sempre"),Percentual = round(as.matrix((table(Drogas_SP$sentido_sozinho)/445)*100),2))
ggplot(sentido_sozinho, aes(x=Categoria, y=Percentual,fill=Categoria)) +
geom_bar(stat="identity",color="black", width=0.8)+
ylim(0,100)+
geom_text(aes(label= paste0(Percentual,"%")), vjust=-1.5, color="black", size=3.5)+
scale_fill_brewer(palette="Set1")+
xlab("\nCor")+
ylab("Percentual (%)\n")+
ggtitle("Percentual de obs. por categorias de 'sentido_sozinho'\n" )+
theme_classic()+
theme(plot.title = element_text(hjust = 0.5))1.5.9 Quantos amigos seus usam drogas?
amigos_drogas <- data.frame(Categoria= c("Nenhum", "Poucos", "Alguns", "A maioria","Não sei"),Percentual = round(as.matrix((table(Drogas_SP$amigos_drogas)/445)*100),2))
ggplot(amigos_drogas, aes(x=Categoria, y=Percentual,fill=Categoria)) +
geom_bar(stat="identity",color="black", width=0.8)+
ylim(0,100)+
geom_text(aes(label= paste0(Percentual,"%")), vjust=-1.5, color="black", size=3.5)+
scale_fill_brewer(palette="Set1")+
xlab("\nCor")+
ylab("Percentual (%)\n")+
ggtitle("Percentual de obs. por categorias de 'amigos_drogas'\n" )+
theme_classic()+
theme(plot.title = element_text(hjust = 0.5))2 Estimação não-paramétrica
tempos<- Drogas_SP$tempo
cens<-Drogas_SP$censura2.1 Estimação da função de sobrevivência S(t) por Kaplan-Meier com transformação log[-log(S(t))]
\(S(t) = Pr(T > t)\)
ekm<-survfit(Surv(tempos,cens)~1,conf.type="log-log")
plot(ekm,xlab="Tempo (anos)",ylab="S(t) Kaplan-Meier",conf.int=T, main="S(t) estimada por Kaplan-Meier com transf. log[-log(S(t))] + IC\n")A probabilidade de sobrevivência nos dá a probabilidade de um um indivíduo do estudo sobreviver após um tempo \(t\) especificado.
O gráfico acima contém as probabilidades de sobrevivência correspondentes a cada tempo de falha ordenado.
Por ser obtida de forma empírica, a função de sobrevivência é uma função de degraus: começa com uma linha horizontal na probabilidade de sobrevivência 1 e e então vai descendo para probabilidades de sobrevivência menores a medida que falhas vão sendo observadas.
A probabilidade de sobrevivência estimada após o primeiro tempo de falha ordenado (11 anos) continua bem próxima de \(1\).
Aos 15 anos ocorre o quinto (e último) tempo de falha ordenado e a probabilidade de sobrevivência estimada fica acima de \(70\%\), ou seja, estima-se que a probabilidade de um estudante não experimentar drogas ilícitas após os 15 anos é superior a \(70\%\).
Conforme os anos passam, percebemos que as quedas nas probabilidades de sobrevivência foram aumentando.
Como o maior tempo (16 anos) é uma censura, o valor de \(S(t)\) a partir desse ponto é indeterminado.
2.1.1 ICs normais com transformacao log[-log(S(t))] para a estimativa de Kaplan-Meier de S(t)
A saída abaixo contém as estimativas pontuais e intervalares (com transformação log-log) de Kaplan-Meier. As colunas são relativas a/ao:
- tempo de falha;
- número de indivíduos em risco de sofrer o evento;
- número de indivíduos que sofreram o evento;
- estimativa da probabilidade de sobrevivência;
- erro padrão da estimativa;
- limite inferior do iC de 95% para a probabilidade de sobrevivência;
- limite superior do IC de 95% para a probabilidade de sobrevivência.
summary(ekm)## Call: survfit(formula = Surv(tempos, cens) ~ 1, conf.type = "log-log")
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 11 445 1 0.998 0.00224 0.984 1.000
## 12 444 5 0.987 0.00547 0.970 0.994
## 13 438 18 0.946 0.01073 0.920 0.963
## 14 285 14 0.900 0.01583 0.864 0.926
## 15 30 5 0.750 0.06261 0.601 0.849
Podemos observar que,
- Aos 11 anos ocorre a primeira falha, o que ocasiona a queda da probabilidade de sobrevivência de \(1\) para \(0.998\).
- Aos 12 anos são observadas mais 5 falhas e \(S(t)\) passa de \(0.998\) para \(0.987\).
- Aos 13 anos é observado o maior número de falhas, 18! \(S(t)\) passa de \(0.987\) para \(0.946\).
- Aos 14 anos são observadas 14 falhas e \(S(t)\) passa de \(0.946\) para \(0.900\).
- Por fim, 5 falhas são observadas aos 15 anos e \(S(t)\) sofre uma queda mais acentuada, indo \(0.900\) para \(0.750\).
- Utilizando a estimação de Kaplan-Meier, podemos afirmar com confiança de 95% que a real probabilidade de não experimentar drogas ilícitas (sobreviver) após os 15 anos está entre \(0.601\) e \(0.849\).
2.2 Estimação da função de sobrevivência S(t) por Nelson-Aalen
Abaixo, temos o gráfico de \(S(t)\) estimada por Nelson-Aalen:
ss<-survfit(coxph(Surv(tempos,cens)~1,method="breslow"))
plot(ss, main="S(t) estimada por Nelson-Aalen\n")- Podemos concluir o mesmo observado para o gráfico de \(S(t)\) estimada por Kaplan-Meier, visto que não se encontram diferenças perceptíveis entre os dois gráficos.
2.2.1 ICs para a estimativa de Nelson-Aalen de S(t)
A saída abaixo contém as estimativas pontuais e intervalares de Nelson-Aalen. As colunas são relativas a/ao:
- tempo de falha;
- número de indivíduos em risco de sofrer o evento;
- número de indivíduos que sofreram o evento;
- estimativa da probabilidade de sobrevivência;
- erro padrão da estimativa;
- limite inferior do iC de 95% para a probabilidade de sobrevivência;
- limite superior do IC de 95% para a probabilidade de sobrevivência.
summary(ss)## Call: survfit(formula = coxph(Surv(tempos, cens) ~ 1, method = "breslow"))
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 11 445 1 0.998 0.00224 0.993 1.000
## 12 444 5 0.987 0.00544 0.976 0.997
## 13 438 18 0.947 0.01055 0.926 0.968
## 14 285 14 0.901 0.01553 0.872 0.932
## 15 30 5 0.763 0.05838 0.657 0.887
Podemos observar que,
- As probabilidades de sobrevivência estimadas por Nelson-Aalen foram iguais ou superiores às estimadas por Kaplan-Meier.
- Os erros-padrão das estimativas por Nelson-Aalen são iguais ou menores aos obtidos por Kaplan-Meier.
- Os intervalos de 95% de confiança por Nelson-Aalen tem menor ampplitude que os intervalos de confiança por Kaplan-Meier.
- A maior diferença entre as probabilidades de sobrevivência estimadas por Nelson-Aalen e Kaplan-Meier foi para \(S(15)\).
- Utilizando a estimação de Nelson-Aalen, podemos afirmar com confiança de 95% que a real probabilidade de não experimentar drogas ilícitas (sobreviver) após os 15 anos está entre \(0.657\) e \(0.887\).
2.2.2 Comparando S(t) por Kaplan-Meier e Nelson-Aalen
plot(ekm,ylab="S(t) estimada",xlab="Tempo (anos)",conf.int=F,col=4)
lines(ss,conf.int=F,col=6)
Legenda<-c("Kaplan-Meier","Nelson-Aalen")
legend("bottomright", Legenda,lwd=c(1,1),cex=1.2,inset=0.00,col=c(4,6),bty="n")- Como comentado anteriormente, as estimativas são bem próximas, com a maior diferença sendo observada a partir dos 15 anos.
2.3 Comparando S(t) estimada por Kaplan-Meier com transformação log[-log(S(t))] entre os sexos masculino e feminino
Existe diferença entre as curvas de sobrevivência de homens e mulheres?
grupos<-Drogas_SP$sexo
grupos <- ifelse(grupos==1,"Masculino","Feminino")
ekm<-survfit(Surv(tempos,cens)~grupos,conf.type="log-log")
plot(ekm,lty=c(2,1),xlab="Tempo (anos)",ylab="S(t) estimada")
legend(1,0.3,lty=c(2,1),c("Feminino","Masculino"),lwd=1,bty="n")Podemos observar que,
A curva do sexo feminino possui menos degraus: algum tempo de falha observado para o sexo masculino não ocorreu para o sexo feminino.
Até os 15 anos, a curva de sobrevivência por Kaplan-Meier do sexo feminino fica acima da curva de sovrevivência do sexo masculino, o que significa que as probabilidades de sobrevivência do sexo feminino são maiores que as do sexo masculino.
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.
2.3.1 ICs com transformacao log[-log(S(t))] para a estimativa de Kaplan-Meier de S(t) para os sexos feminino e masculino
summary(ekm)## Call: survfit(formula = Surv(tempos, cens) ~ grupos, conf.type = "log-log")
##
## grupos=Feminino
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 12 223 1 0.996 0.00447 0.969 0.999
## 13 221 8 0.959 0.01323 0.924 0.979
## 14 134 7 0.909 0.02230 0.854 0.944
## 15 7 2 0.650 0.15608 0.277 0.865
##
## grupos=Masculino
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 11 222 1 0.995 0.00449 0.968 0.999
## 12 221 4 0.977 0.00996 0.947 0.991
## 13 217 10 0.932 0.01685 0.890 0.959
## 14 151 7 0.889 0.02264 0.836 0.926
## 15 23 3 0.773 0.06547 0.612 0.874
- Como foi observado no gráfico, as probabilidades de sobrevivência para os dois grupos não diferem muito nos tempos de falha em comum, com exceção do tempo de falha \(t=15\).
- Para o tempo de falha \(t=15\), a probabilidade de sobrevivência estimada para o sexo feminino foi \(0.650\), enquanto a do sexo masculino foi \(0.773\).
- Essa diferença pode estar relacionada ao número de indivíduos sob risco de sofrer o evento: para ambos os sexos o número de eventos em \(t=15\) é similar (2 e 3). Entretanto, o número de indivíduos do sexo masculino sob risco de sofrer o evento é 3 vezes maior que o do sexo feminino, afetando a probabilidade estimada.
\(\hat{S}(t) = \prod_{t_i<t}{1- \frac{d_i}{n_i}}\)
onde \(d_i\) é o número de eventos em \(t_i\) e \(n_i\) o número de indivíduos sob risco em \(t_i\).
2.3.2 Teste logrank para comparar as curvas de sobrevivência nos dois grupos
survdiff(Surv(tempos,cens)~grupos,rho=0)## Call:
## survdiff(formula = Surv(tempos, cens) ~ grupos, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## grupos=Feminino 223 18 19.8 0.171 0.343
## grupos=Masculino 222 25 23.2 0.147 0.343
##
## Chisq= 0.3 on 1 degrees of freedom, p= 0.6
O teste de logrank tem estatística de teste igual a \(0.3\), e p-valor \(0.6\). Logo, ao nível de 5%, não rejeitamos a hipótese nula de igualdade entre as curvas de sobrevivência dos sexos masculino e feminino.
3 Estimação paramétrica
3.1 Estimação pela distribuição Exponencial
- A forma da Exponencial utilizada é dada por:
\(\hat{f(t)}= \frac{1}{\hat{\alpha}} e^{\frac{-t}{\hat{\alpha}}}\)
3.1.1 EMV e IC do parâmetro de escala da Exponencial
Utilização do Método Delta: O R ajusta uma distribuição Exponencial com parâmetro \(\theta\) tal que, na nossa parametrização, o parâmetro de escala \(\alpha\) é obtido por \(\alpha = e^{\theta}\).
A estimativa de máxima verossimilhança para \(\alpha\), o parâmetro de escala da distribuição Exponencial, sua variância e o IC de \(95\%\) são dados por:
tempos<- Drogas_SP$tempo
cens<-Drogas_SP$censura
# Estimacao usando a distribuicao Exponencial:
ajuste_exp<-survreg(Surv(tempos,cens)~1,dist='exponential')
# como a funçãoo usa outra parametrização da exponencial:
alfa_exp<-exp(ajuste_exp$coefficients[1]) # EMV do parametro de escala da Exponencial | estimativa pontual do parâmetro alfa
varalfa_exp<-ajuste_exp$var[1, 1]*(exp(ajuste_exp$coefficient[1]))^2 # EMV da variancia do parametro de escala
tab_exp <- data.frame(EMV = exp(ajuste_exp$coefficients[1]), Variancia=ajuste_exp$var[1, 1]*(exp(ajuste_exp$coefficient[1]))^2 , LI= alfa_exp-qnorm(0.975)*sqrt(varalfa_exp),LS=alfa_exp+qnorm(0.975)*sqrt(varalfa_exp) )
tab_exp <- round(tab_exp,4)
knitr::kable(tab_exp,row.names = FALSE)| EMV | Variancia | LI | LS |
|---|---|---|---|
| 141.814 | 467.7023 | 99.4269 | 184.201 |
- Com base nas estimativas obtidas, podemos estimar quantidades de interesse para os dados.
3.1.2 Tempo médio até a experimentação de drogas ilícitas
- Uma vez estimados os parâmetros, o tempo médio até a experimentação de drogas ilícitas pode ser estimado:
#Estimando o tempo medio ate a reincidencia:
E_T= round(data.frame(E_T = alfa_exp),2)
knitr::kable(E_T, row.names = FALSE)| E_T |
|---|
| 141.81 |
- Estima-se o tempo médio até a experimenação de drogas ilícitas seja de \(141.81\) anos, visto que na distribuição Exponencial \(\hat{E}(T)=\hat{\alpha}\).
3.1.3 Estimação da função de taxa de falha
# Estimando a função de taxa de falha:
lambda_exp<-function(t){1/alfa_exp*(t^0)}
# Estimando a funcao de taxa de falha:
plot(lambda_exp,ylab=expression(lambda(t)),xlab="t (em anos)",xlim=c(0,16),
main="Função da taxa de falha - Exponencial") #plotando a funcao de taxa de falha- A função de taxa de falha, que nesse caso é constante (propriedade de falta de memória da distribuição Exponencial), é estimada como sendo \(0.007\).
- Esse é o risco instantâneo de uso de drogas ilícitas, independente do tempo (em anos).
- Por exemplo, um estudante com 5 anos tem o mesmo risco instantâneo estimado de usar drogas ilícitas que um estudante de 15 anos.
3.1.4 Estimação dos percentis 100p%
# Estimando o percentil 100p%:
perc_exp<-function(p){-alfa_exp*log(1-p)}#criando uma funcao no R para o percentil
perc_exp_est <- data.frame(cbind(perc_exp(p=0.1),perc_exp(p=0.2),perc_exp(p=0.3)))
colnames(perc_exp_est)<- c("10%","20%","30%")
perc_exp_est <- round(perc_exp_est,2)
knitr::kable(perc_exp_est, row.names = FALSE)| 10% | 20% | 30% |
|---|---|---|
| 14.94 | 31.64 | 50.58 |
- Aos \(14.94\) anos, estima-se que 10% terão utilizado drogas ilícitas;
- Aos \(31.64\) anos, estima-se que 20% terão utilizado drogas ilícitas;
- Aos \(50.58\) anos, estima-se que 30% terão utilizado drogas ilícitas.
O gráfico abaixo mostra como os percentis variam em função do valor de p.
plot(perc_exp,ylab="tp (em anos)",xlab="p",
main="Percentil 100% - Exponencial") #plotando o percentil versus p3.1.5 Estimação da função de sobrevivência
#Estimando a funcao de sobrevivencia:
S_exp<-function(t){exp(-t/alfa_exp)}#criando uma funcao no R para a sobrevivencia
var_S_exp<-function(t){varalfa_exp*((t/alfa_exp^2)*exp(-t/alfa_exp))^2}#criando uma funcao no R para a variancia da funcao de sobrêvivencia
f_sobrev <- data.frame(est_pontual = c(S_exp(t=15),S_exp(t=18),S_exp(t=22)),
LI = c(S_exp(t=15)-qnorm(0.975)*sqrt(var_S_exp(t=15)),S_exp(t=18)-qnorm(0.975)*sqrt(var_S_exp(t=18)),S_exp(t=22)-qnorm(0.975)*sqrt(var_S_exp(t=22))),
LS = c(S_exp(t=15)+qnorm(0.975)*sqrt(var_S_exp(t=15)), S_exp(t=18)+qnorm(0.975)*sqrt(var_S_exp(t=18)), S_exp(t=22)+qnorm(0.975)*sqrt(var_S_exp(t=22))))
f_sobrev <- round(f_sobrev,4)
rownames(f_sobrev)<-c("t=15","t=18","t=22")
knitr::kable(f_sobrev,row.names = TRUE)| est_pontual | LI | LS | |
|---|---|---|---|
| t=15 | 0.8996 | 0.8712 | 0.9281 |
| t=18 | 0.8808 | 0.8474 | 0.9142 |
| t=22 | 0.8563 | 0.8166 | 0.8960 |
Estimativas pontuais e intervalares (95% de confiança) foram obtidas para os \(t\) igual a 10, 12 e 15 anos, sendo:
\(\hat{S}(15)=0.8996 \quad(0.8712; 0.9281):\) estima-se que a chance de um indivíduo não ter experimentado drogas ilícitas após os 15 anos é de \(89.96\%\).
\(\hat{S}(18)=0.8808 \quad(0.8474; 0.9142):\) estima-se que a chance de um indivíduo não ter experimentado drogas ilícitas após os 18 anos é de \(88.08\%\).
\(\hat{S}(22)=0.8563 \quad(0.8166; 0.8960):\) estima-se que a chance de um indivíduo não ter experimentado drogas ilícitas após os 22 anos é de \(85.63\%\).
É interessante observar a relação entre e \(\hat{S}(15)=0.8996\) e \(\hat{t}_{0.1004} \approx 15\) , estimado anteriormente: Se para 15 anos estima-se que \(89.96\%\) dos pacientes não terão experimentado drogas ilícitas, esse é o tempo em que, aproximadamente, 10% deles já terão experimentado esse evento.
O gráfico mostra como a função de sobrevivência \(S(t)\) varia ao longo de \(t\) (em anos):
plot(S_exp,ylab="S(t)",xlab="t (em anos)",xlim=c(0,100),
main="Função de sobrevivência - Exponencial") #plotando a funcao de sobrêvivencia3.2 Estimação usando a distribuição Weibull
- A forma da Weibull utilizada é dada por:
\(\hat{f}(t)=\frac{\hat{\gamma}}{\hat{\alpha}\hat{\gamma}}t^{\hat{\gamma}-1}e^{-(\frac{t}{\hat{\alpha}})^\hat{\gamma}}\)
3.2.1 EMVs e ICs para os parâmetros da Weibull
- Assim como na distribuição Exponencial, a parametrização utilizada pelo R para a distribuição Weibull é diferente da que utilizamos.
- O R ajusta uma distribuição Weibull com parâmetros \(\theta_1\) e \(\theta_2\) tais que, na nossa parametrização, \(\alpha=exp(\theta_1)\) e \(\\gamma=\frac{1}{\theta_2}\).
- As estimativas de máxima verossimilhança para \(\alpha\) e \(\gamma\), suas variâncias e respectivos intervalos de confiança são:
ajuste_wei<-survreg(Surv(tempos,cens)~1,dist='weibull')
alfa_wei<-exp(ajuste_wei$coefficients[1]) # EMV do parametro de escala da Weibull
gama_wei<-1/ajuste_wei$scale# EMV do parametro de forma da Weibull
varalfa_wei<-ajuste_wei$var[1,1]*exp(ajuste_wei$coefficients[1])^2
vargama_wei<-ajuste_wei$var[2,2]*(-1/ajuste_wei$scale^2)^2
tab_weib <- data.frame(EMV = c(gama_wei,alfa_wei),
Variancia = c(varalfa_wei,vargama_wei),
LI = c(alfa_wei-qnorm(0.975)*sqrt(varalfa_wei),gama_wei-qnorm(0.975)*sqrt(vargama_wei)),
LS = c(alfa_wei+qnorm(0.975)*sqrt(varalfa_wei),gama_wei+qnorm(0.975)*sqrt(vargama_wei)))
rownames(tab_weib) <-c("alfa","gama")
knitr::kable(tab_weib)| EMV | Variancia | LI | LS | |
|---|---|---|---|---|
| alfa | 15.53066 | 0.0833058 | 15.65154 | 16.78294 |
| gama | 16.21724 | 768.3947748 | -38.79941 | 69.86073 |
3.2.2 Tempo médio até a experimentação de drogas ilícitas
- Uma vez estimados os parâmetros, o tempo médio até a experimentação de drogas ilícitas pode ser estimado:
#Estimando o tempo medio ate a reincidencia:
E_T= round(data.frame(E_T = alfa_wei*gamma(1+1/gama_wei)),2)
knitr::kable(E_T, row.names = FALSE)| E_T |
|---|
| 15.68 |
- Portanto, estima-se o tempo médio até a experimenação de drogas ilícitas seja de \(15.68\) anos.
3.2.3 Estimação da função de taxa de falha
# Estimando a funcao de taxa de falha:
lambda_wei<-function(t){gama_wei/alfa_wei^gama_wei*t^(gama_wei-1)}#criando uma funcao no R para a taxa de falha
tab_lambda_wei <- round(cbind.data.frame(lambda_wei = c(lambda_wei(t=15), lambda_wei(t=18),lambda_wei(t=22))),2)
rownames(tab_lambda_wei) <- c("15","18","22")
knitr::kable(tab_lambda_wei)| lambda_wei | |
|---|---|
| 15 | 0.31 |
| 18 | 4.36 |
| 22 | 80.49 |
- O risco instantâneo de experimentação de drogas ilícitas é estimado em \(0.31\), \(4.36\) e \(80.49\) em 15, 18 e 22 anos, respectivamente.
plot(lambda_wei,ylab=expression(lambda(t)),xlab="t (em anos)",xlim=c(0,40),
main = "Função da taxa de falha - Weibull") #plotando a funcao de taxa de falha- A função de taxa de falha estimada pelo ajuste da distribuição Weibull é crescente, diferente do observado quando se utilizou a distribuição Exponencial.
- Logo, os indivíduos tornam-se mais suscetíveis a uma falha instantânea à medida em que o tempo (em anos) aumenta.
3.2.4 Estimação dos percentis 100p%
# Estimando o percentil 100p%:
perc_wei<-function(p){alfa_wei*(-log(1-p))^(1/gama_wei)}#criando uma funcao no R para o percentil
perc_wei_est <- data.frame(cbind(perc_wei(p=0.1),perc_wei(p=0.2),perc_wei(p=0.3)))
colnames(perc_wei_est)<- c("10%","20%","30%")
perc_wei_est <- round(perc_wei_est,2)
knitr::kable(perc_wei_est, row.names = FALSE)| 10% | 20% | 30% |
|---|---|---|
| 14.03 | 14.72 | 15.18 |
- Aos \(14.03\) anos, estima-se que 10% terão utilizado drogas ilícitas;
- Aos \(14.72\) anos, estima-se que 20% terão utilizado drogas ilícitas;
- Aos \(15.18\) anos, estima-se que 30% terão utilizado drogas ilícitas.
O gráfico abaixo mostra como os percentis variam em função do valor de p.
plot(perc_wei,ylab="tp (em anos)",xlab="p",
main="Percentil 100% - Weibull") #plotando o percentil versus p3.2.5 Estimação da função de sobrevivência
#Estimando a funcao de sobrevivencia:
S_wei<-function(t){exp(-(t/alfa_wei)^gama_wei)}#criando uma funcao no R para a sobrevivencia
f_sobrev <- data.frame(est_pontual = c(S_wei(t=15),S_wei(t=18),S_wei(t=22)))
f_sobrev <- round(f_sobrev,4)
rownames(f_sobrev)<-c("t=15","t=18","t=22")
knitr::kable(f_sobrev,row.names = TRUE)| est_pontual | |
|---|---|
| t=15 | 0.7425 |
| t=18 | 0.0064 |
| t=22 | 0.0000 |
Estimativas pontuais foram obtidas para os \(t\) igual a 15, 18 e 22 anos, sendo:
- \(\hat{S}(15)=0.7425\) estima-se que a chance de um indivíduo não ter experimentado drogas ilícitas após os 15 anos é de \(74.25\%\).
- \(\hat{S}(18)=0.0064\) estima-se que a chance de um indivíduo não ter experimentado drogas ilícitas após os 18 anos é de \(0.64\%\).
- \(\hat{S}(22)=0.000\) estima-se que a chance de um indivíduo não ter experimentado drogas ilícitas após os 22 anos é de \(0\%\).
O gráfico mostra como a função de sobrevivência \(S(t)\) varia ao longo de \(t\) (em anos):
plot(S_wei,ylab="S(t)",xlab="t (em anos)",xlim=c(0,40),
main="Função de sobrevivência - Weibull") #plotando a funcao de sobrevivencia3.3 Estimação usando a distribuição Log-normal
- A forma da Log-normal utilizada é dada por:
\(\hat{f}(t)=\frac{1}{\sqrt{2\pi}t\hat{\sigma}}exp{-\frac{1}{2}( \frac{log(t)-\hat{\mu}}{\hat{\sigma}})^2 }\)
3.3.1 EMVs e Ics para parâmetros da Log-normal
- As estimativas de máxima verossimilhança para \(\mu\) e \(\sigma\), suas variâncias e respectivos intervalos de confiança são:
ajuste_logn<-survreg(Surv(tempos,cens)~1,dist='lognorm')
mu_logn<-ajuste_logn$coefficients[1] # EMV do parametro de locacao da Log-normal
sigma_logn<-ajuste_logn$scale# EMV do parametro de escala da Log-normal
var_mu<-ajuste_logn$var[1,1]# variancia da estimativa do parametro de locacao
var_sigma<-ajuste_logn$var[2,2]
tab_logn <- data.frame(EMV = c(mu_logn,sigma_logn),
Variancia = c(var_mu,var_sigma),
LI = c(mu_logn-qnorm(0.975)*sqrt(var_mu),mu_logn+qnorm(0.975)*sqrt(var_mu)),
LS = c(sigma_logn-qnorm(0.975)*sqrt(var_sigma),sigma_logn+qnorm(0.975)*sqrt(var_sigma)))
rownames(tab_logn) <-c("mu","sigma")
knitr::kable(tab_logn)| EMV | Variancia | LI | LS | |
|---|---|---|---|---|
| mu | 2.7944794 | 0.0005022 | 2.750557 | -0.1043686 |
| sigma | 0.1234007 | 0.0135050 | 2.838402 | 0.3511701 |
- Essas estimativas podem ser utilizadas para se obter as demais medidas de interesse.
3.3.2 Tempo médio até a experimentação de drogas ilícitas
- Uma vez estimados os parâmetros, o tempo médio até a experimentação de drogas ilícitas pode ser estimado:
#Estimando o tempo medio ate a reincidencia:
E_T=round(data.frame(E_T=exp(mu_logn+(sigma_logn^2)/2)),2)
knitr::kable(E_T,row.names = FALSE)| E_T |
|---|
| 16.48 |
- Portanto, estima-se o tempo médio até a experimenação de drogas ilícitas seja de \(16.48\) anos.
3.3.3 Estimação da função de taxa de falha
f_logn<-function(t){(1/(sqrt(2*pi)*t*sigma_logn))*exp(-0.5*((log(t)-mu_logn)/sigma_logn)^2)}#criando uma funcao no R para a densidade
S_logn<-function(t){pnorm((-log(t)+mu_logn)/sigma_logn)}
lambda_logn<-function(t){f_logn(t)/S_logn(t)}#criando uma funcao no R para a taxa de falha
tab_lambda_logn <- round(cbind.data.frame(lambda_logn = c(lambda_logn(t=15), lambda_logn(t=18),lambda_logn(t=22))),2)
rownames(tab_lambda_logn)<- c("15","18","22")
knitr::kable(tab_lambda_logn)| lambda_logn | |
|---|---|
| 15 | 0.22 |
| 18 | 0.61 |
| 22 | 1.01 |
O risco instantâneo de experimentação de drogas ilícitas é estimado em \(0.22\), \(0.61\) e \(1.01\) em 15, 18 e 22 anos, respectivamente.
O gráfico da função da taxa de falha ao longo do tempo t(em anos) pode ser visto abaixo:
plot(lambda_logn,ylab=expression(lambda(t)),xlab="t (em anos)",xlim=c(0,45),
main = " Função da taxa de falha - Log-normal") #plotando a funcao de taxa de falha3.3.4 Estimação dos percentis 100p%
# Estimando o percentil 100p%:
perc_logn<-function(p){exp(qnorm(p)*sigma_logn+mu_logn)}#criando uma funcao no R para o percentil
perc_logn_est <- data.frame(cbind(perc_logn(p=0.1),perc_logn(p=0.2),perc_logn(p=0.3)))
colnames(perc_logn_est)<- c("10%","20%","30%")
perc_logn_est <- round(perc_logn_est,2)
knitr::kable(perc_logn_est, row.names = FALSE)| 10% | 20% | 30% |
|---|---|---|
| 13.96 | 14.74 | 15.33 |
- Aos \(13.96\) anos, estima-se que 10% terão utilizado drogas ilícitas;
- Aos \(14.74\) anos, estima-se que 20% terão utilizado drogas ilícitas;
- Aos \(15.33\) anos, estima-se que 30% terão utilizado drogas ilícitas.
O gráfico abaixo mostra como os percentis variam em função do valor de p.
plot(perc_logn,ylab="t (em anos)",xlab="p",
main = " Percentil 100% - Log-normal") #plotando o percentil versus p3.3.5 Estimação da função de sobrevivência
#Estimando a funcao de sobrevivencia:
S_logn<-function(t){pnorm((-log(t)+mu_logn)/sigma_logn)}#criando uma funcao no R para a sobrevivencia
f_sobrev <- data.frame(est_pontual = c(S_logn(t=15),S_logn(t=18),S_logn(t=22)))
f_sobrev <- round(f_sobrev,4)
rownames(f_sobrev)<-c("t=15","t=18","t=22")
knitr::kable(f_sobrev,row.names = TRUE)| est_pontual | |
|---|---|
| t=15 | 0.7582 |
| t=18 | 0.2186 |
| t=22 | 0.0081 |
Estimativas pontuais foram obtidas para os \(t\) igual a 15, 18 e 22 anos, sendo:
- \(\hat{S}(15)=0.7582\) estima-se que a chance de um indivíduo não ter experimentado drogas ilícitas após os 15 anos é de \(75.82\%\).
- \(\hat{S}(18)=0.2186\) estima-se que a chance de um indivíduo não ter experimentado drogas ilícitas após os 18 anos é de \(21.86\%\).
- \(\hat{S}(22)=0.0081\) estima-se que a chance de um indivíduo não ter experimentado drogas ilícitas após os 22 anos é de \(0.81\%\).
O gráfico mostra como a função de sobrevivência \(S(t)\) varia ao longo de \(t\) (em anos):
plot(S_logn,ylab="S(t)",xlab="t (em anos)",xlim=c(0,100),
main = " Função de sobrevivência - Log-normal") #plotando a funcao de sobrevivencia3.4 Comparando os resultados das três distribuições
Comparação entre os resultados da estimação paramétrica para as distribuições Exponencial, Weibull e Log-normal
3.5 Plotando as funções das três distribuições conjuntamente
#Funcao densidade:
f_exp<-function(t){(1/alfa_exp)*exp(-t/alfa_exp)} #funcao densidade da Exponencial
f_wei<-function(t){(gama_wei/alfa_wei^gama_wei)*t^(gama_wei-1)*exp(-(t/alfa_wei)^gama_wei)} #funcao densidade da Weibull
plot(f_exp,type="l",col=4,ylab="f(t)",xlab="t",xlim=c(0,40), ylim=c(0,0.4), main= "Comparando funções de densidade")#plotando as funcoes
t=seq(0,40,0.1)
lines(t,f_wei(t),type="l",col=6)
lines(t,f_logn(t),type="l",col=9)
Legenda<-c(expression(Exponencial),expression(Weibull),expression(Log-normal))
legend("topright", Legenda,lwd=c(1,1,1),cex=1.2,inset=0.00,col=c(4,6,9),bty="n")#Funcao de taxa de falha:
plot(lambda_exp,type="l",col=4,ylab=expression(lambda(t)),xlab="t",xlim=c(0,40),ylim=c(0,1), main= "Comparando funções de taxa de falha")#plotando as funcoes
t=seq(0,40,0.1)
lines(t,lambda_wei(t),type="l",col=6)
lines(t,lambda_logn(t),type="l",col=9)
Legenda<-c(expression(Exponencial),expression(Weibull),expression(Log-normal))
legend("topright", Legenda,lwd=c(1,1,1),cex=1.2,inset=0.00,col=c(4,6,9),bty="n")#Funcao de sobrevivencia:
plot(S_exp,type="l",col=4,ylab="S(t)",xlab="t",xlim=c(0,40), main= "Comparando funções de sobrevivência")#plotando as funcoes
t=seq(0,8040,0.1)
lines(t,S_wei(t),type="l",col=6)
lines(t,S_logn(t),type="l",col=9)
Legenda<-c(expression(Exponencial),expression(Weibull),expression(Log-normal))
legend("topright", Legenda,lwd=c(1,1,1),cex=1.2,inset=0.00,col=c(4,6,9),bty="n")3.6 Teste de Razão de Verossimilhança
- Hipóteses:
\(H0: \text{O modelo de interesse é adequado}\)
\(H1: \text{O modelo de interesse não é adequado}\)
- Estimou-se os parâmetros por MV para cada uma das três distribuições (Exponencial, Weibull e Lognormal), e para a distribuição de referência (Gama generalizada).
- E então, as estatisticas TRV e seus p-valores foram calculados:
ajuste_exp<-survreg(Surv(tempos,cens)~1,dist='exponential')
ajuste_wei<-survreg(Surv(tempos,cens)~1,dist='weibull')
ajuste_logn<-survreg(Surv(tempos,cens)~1,dist='lognorm')
ajuste_gamg<-flexsurvreg(Surv(tempos,cens)~1,dist="gengamma")
TRV_exp=2*(ajuste_gamg$loglik-ajuste_exp$loglik[2])
TRV_wei=2*(ajuste_gamg$loglik-ajuste_wei$loglik[2])
TRV_logn=2*(ajuste_gamg$loglik-ajuste_logn$loglik[2])
TRV <- data.frame(pvalor= c(1-pchisq(TRV_exp,df=2),1-pchisq(TRV_wei,df=1),1-pchisq(TRV_logn,df=1)))
rownames(TRV) <- c("Exponencial","Weibull","Log-normal")
knitr::kable(TRV)| pvalor | |
|---|---|
| Exponencial | 0.0000000 |
| Weibull | 0.1511699 |
| Log-normal | 0.8446717 |
- Analisando os p-valores concluímos que, ao nível de \(5\%\), as distribuições weibull e log-normal são adequadas.