Lendo a base de dados

DIR = 'data'
FILE = 'amostra_semprocesso.csv'
URL = file.path(DIR,FILE)

test.sp = read.csv(URL, row.names = 1, encoding = 'UTF-8')
knitr::kable(head(test.sp))

Vamos analisar a proporção pela classe de interesse que iremos trabalhar.

knitr::kable(table(test.sp$Êxito))

Realizando o attach das variaveis do banco:

names(test.sp)

##  [1] "Ano.de.Início"                       "Mês.de.Início"                      
##  [3] "Dia.de.Início"                       "Acidente.com.morte...lesões.físicas"
##  [5] "Ano.de.Fabricação"                   "Ano.Modelo"                         
##  [7] "Distribuição"                        "Êxito"                              
##  [9] "Fase.Atual"                          "Modelo.Veículo"                     
## [11] "Natureza"                            "Natureza.Financeira"                
## [13] "Regional"                            "UF.Foro"                            
## [15] "time"

attach(test.sp)

Definindo a base de dados

A base de dados é composta de um conjunto de informações a respeito de processos juridicos (Civis e Trabalhista) referentes a empresa Ford, contendo a informção do tempo de duração daquele processo juridico, mais outras informações que serão listadas abaixo:

• Ano.de.Início: Ano de inicio daquele processo juridico;

• Mês.de.Início: Mês de inicio daquele processo juridico;

• Dia.de.Início: Dia de inicio daquele processo juridico;

• Ano.de.Fabricação: Ano de fabricação do veículo relacionado ao processo;

• Ano.Modelo: Ano do modelo do carro relacionado ao processo;

• Acidente.com.morte…lesões.físicas: Variavel dicotomica, com sim ou não;

  • sim: Houve acidente com morte ou lesoes fisicas relacionado ao processo
  • não: Não houve acidente com morte ou lesoes fisicas relacionado ao processo

• Distribuição: Numero de distribuição daquele processo;

• Fase.Atual: Fase atual que o processo encontra-se no momento;

• Modelo.Veículo: Modelo do veículo relacionado ao processo;

• Natureza: Natureza daquele processo:

  • Trabalhista: processos de natureza trabalhista
  • Civel: processos de natureza Civel

• Natureza.Financeira: Natureza financeira relacionado aquele processo:

  • Ativo: Natureza financeira Ativa
  • Neutro: Natureza financeira Neutro
  • Passivo: Natureza financeira Passivo

• Regional: Regional do processo;

• UF.Foro: UF do foro relacionado ao processo;

• Êxito: Variável refernete ao êxito daquele processo para empresa, ou não:

  • (1) Não êxito para a empresa

• Time: tempo de duração ate o momento daquele processo.

Definindo o problema

A nossa tarefa ao decorrer desse trabalho será elaborar um modelo de Análise de Sobrevivência com o intuito de predizer a probabilidade do processo sobreviver ate um tempo t sem resultar em perda para a empresa. Essa informação é bastante importante para uma empresa, com o intuito de saber em qual tempo, em media ou em termos da mediana, que aquele processo tem de sobreviver sem resultar em perda, informação util com o intuito de monitorar esse tempo.

Será utilizado a variável Êxito, onde a censura será qualquer processo que não ocasionou em perda para empresa e a falha será o não exito para a empresa, ou seja, a perca daquele processo.

names(test.sp)

##  [1] "Ano.de.Início"                       "Mês.de.Início"                      
##  [3] "Dia.de.Início"                       "Acidente.com.morte...lesões.físicas"
##  [5] "Ano.de.Fabricação"                   "Ano.Modelo"                         
##  [7] "Distribuição"                        "Êxito"                              
##  [9] "Fase.Atual"                          "Modelo.Veículo"                     
## [11] "Natureza"                            "Natureza.Financeira"                
## [13] "Regional"                            "UF.Foro"                            
## [15] "time"

Separando as informações

min(test.sp$time)

## [1] 0

Podemos notar que o minimo do tempo é igual a 0 dias, para resolver essa situação, vamos somar 1 dia a todas as informações na base dados. Então:

tempos <- test.sp$time + 1
cens <- test.sp$`Êxito`#as.factor(cens)

Log-Rank Teste

survdiff(Surv(tempos,cens)~Natureza.Financeira,rho=0)  #Teste log-rank

## Call:
## survdiff(formula = Surv(tempos, cens) ~ Natureza.Financeira, 
##     rho = 0)
## 
##                                 N Observed Expected (O-E)^2/E (O-E)^2/V
## Natureza.Financeira=ATIVO      48       39     35.2     0.399     0.402
## Natureza.Financeira=NEUTRO   2703     2130    537.4  4719.478  5149.694
## Natureza.Financeira=PASSIVO 18617     8839  10435.3   244.198  4853.721
## 
##  Chisq= 5151  on 2 degrees of freedom, p= <2e-16

Considerando um nível de significancia de 5%, para o teste de log-rank, sob os dados observados, temos então que rejeitamos a hipotese nula, ou seja, temos evidencias suficientes para afirmar que existe uma diferença estatisticamente significante entre os grupos da variavel Natureza.Financeira.

Teste de Peto

survdiff(Surv(tempos,cens)~Natureza.Financeira,rho=1) #Teste Petoajust1

## Call:
## survdiff(formula = Surv(tempos, cens) ~ Natureza.Financeira, 
##     rho = 1)
## 
##                                 N Observed Expected (O-E)^2/E (O-E)^2/V
## Natureza.Financeira=ATIVO      48     26.6       17      5.44      8.11
## Natureza.Financeira=NEUTRO   2703   1779.5      437   4120.49   5350.07
## Natureza.Financeira=PASSIVO 18617   5004.9     6357    287.50   5273.54
## 
##  Chisq= 5364  on 2 degrees of freedom, p= <2e-16

Considerando também um nivel de significancia de 5%, agora com respeito ao teste de Peto, sob os dados observados, temos que a hipotese nula é rejeitada, ou seja, temos evidencias suficientes para afirmar que existe realmente diferença entre os grupos da variavel Natureza.Financeira.

Modelagem

ajust1<-survreg(Surv(tempos,cens)~1,dist="exponential")
ajust2<-survreg(Surv(tempos,cens)~1,dist='weibull')
ajust3<-survreg(Surv(tempos,cens)~1,dist='lognormal')

ekm<-survfit(Surv(tempos,cens)~1) #estimador de kaplan-meier
time<-ekm$time #tempos

#Estimativas do par?metro no modelo Exponencial
lambda<-exp(ajust1$coefficients[1])
lambdae<-1/lambda

#Estimativas do par?metro no modelo Weibull
lambda1<-exp(ajust2$coefficients[1])
lambdaw<-1/lambda1
gama<-1/ajust2$scale

#Estimativas do par?metro no modelo Lognormal
mu<- ajust3$coefficients
sigma<- ajust3$scale

st<-ekm$surv #sobreviv?ncia kaplan meier
ste<- exp(-time*lambdae) #sobreviv?ncia exponencial
stw<- exp(-(time*lambdaw)^gama) #sobreviv?ncia weibull
stln<- pnorm((-log(time)+ mu)/sigma) #sobreviv?ncia lognormal

Comparando modelos

par(mfrow=c(1,3))
plot(ste,st,pch=16,ylim=range(c(0.0,1)), xlim=range(c(0,1)), ylab = "S(t): Kaplan-Meier",
     xlab="S(t): exponencial")
lines(c(0,1), c(0,1), type="l", lty=1)
plot(stw,st,pch=16,ylim=range(c(0.0,1)), xlim=range(c(0,1)), ylab = "S(t): Kaplan-Meier",
     xlab="S(t): Weibull")
lines(c(0,1), c(0,1), type="l", lty=1)
plot(stln,st,pch=16,ylim=range(c(0.0,1)), xlim=range(c(0,1)), ylab = "S(t): Kaplan-Meier",
     xlab="S(t): log-normal")
lines(c(0,1), c(0,1), type="l", lty=1)

par(mfrow=c(1,3))
plot(ekm, conf.int=F, xlab="Tempos", ylab="S(t)")
lines(c(0,time),c(1,ste), lty=2)
legend("bottomleft",lty=c(1,2),c("Kaplan-Meier", "Exponencial"),bty="n",cex=0.8)
plot(ekm, conf.int=F, xlab="Tempos", ylab="S(t)")
lines(c(0,time),c(1,stw), lty=2)
legend("bottomleft",lty=c(1,2),c("Kaplan-Meier", "Weibull"),bty="n",cex=0.8)
plot(ekm, conf.int=F, xlab="Tempos", ylab="S(t)")
lines(c(0,time),c(1,stln), lty=2)
legend("bottomleft",lty=c(1,2),c("Kaplan-Meier", "Log-normal"),bty="n",cex=0.8)

#Exponencial
lste=-log(st)
#Weibull:
lstw=log(-log(st))
#Lognormal:
lstln=qnorm(1 - st)

min(lste)

## [1] 0.0197081

max(lste)

## [1] Inf

lste[lste == Inf] = sort(lste, T)[2]
lstw[lstw == Inf] = sort(lstw, T)[2]
lstln[lstln == Inf] = sort(lstln, T)[2]

par(mfrow=c(1,3))
plot(time, lste,pch=16,xlab="tempos",ylab="-log(S(t))",main="Exponencial")
abline(lm(lste ~ time))
plot(log(time),lstw,pch=16,xlab="log(tempos)",ylab="log(-log(S(t)))",main="Weibull")
abline(lm(lstw ~ log(time)))
plot(log(time),lstln, pch=16,xlab="log(tempos)", ylab=expression(Phi^-1 * (1-S(t))),main="Log-Normal")
abline(lm(lstln ~ log(time)))

fitgg <- flexsurvreg(formula = Surv(tempos, cens) ~ 1,dist="gengamma")
fitw <- flexsurvreg(formula = Surv(tempos, cens) ~ 1, dist="weibull")
fite <- flexsurvreg(formula = Surv(tempos, cens) ~ 1, dist="exp")
fitl <- flexsurvreg(formula = Surv(tempos, cens) ~ 1, dist="lnorm")

#Obten??o das log-verossimilhan?as
loglikm1<-fite$loglik[1]  #log-verossimilhan?a modelo exponencial   
loglikm2<-fitw$loglik[1]  #log-verossimilhan?a modelo weibull
loglikm3<-fitl$loglik[1]  #log-verossimilhan?a modelo lognormal
loglikG<-fitgg$loglik[1]   #log-verossimilhan?a modelo generalizado 

#Obten??o das estat?sticas de teste

RVe<-2*(loglikG-loglikm1) #Estat?stica da raz?o de verossimilhan?as para o modelo exponencial
RVw<-2*(loglikG-loglikm2) #Estat?stica da raz?o de verossimilhan?as para o modelo weibull
RVln<-2*(loglikG-loglikm3) #Estat?stica da raz?o de verossimilhan?as para o modelo lognormal

#Obten??o dos p-valores dos testes            
pRV1<-1-pchisq(RVe,2) #Exponencial
pRV2<-1-pchisq(RVw,1) #Weibull
pRV3<-1-pchisq(RVln,1) #Lognormal

c('pvalue_exp'=pRV1,'pvalue_weibull'=pRV2,'pvalue_lognormal'=pRV3)

##       pvalue_exp   pvalue_weibull pvalue_lognormal 
##      0.000000000      0.003475434      0.000000000

Resultado

Considerando o método 01 para comparação de desempenho das curvas de sobrevivência de cada modelo com o kaplan-meyer, podemos notar que o modelo que performou melhor foi o modelo Weibull, apresentando quase toda completude de sua curva sobre posta a curva do kaplan-meyer. O segundo modelo que melhor desempenhou foi o modelo Exponencial, porem para tempos maiores ele apresentou um vies maior. Ao analisarmos o modelo Log-Normal, podemos notar que apartir de um ponto as probabilidades estimadas pelo modelo divergiram de forma consideravel do kaplan-meyer,

Com respeito ao método 2, podemos notar que em sua maioria, os pontos do modelo Wibull se aproximaram mais da reta apartir de um certo ponto, porém para tempos menores ele apresentou uma consideravel quantidade de pontos distantes da reta y=x. Além disso, temos que o modelo exponencial apresentou grande parte dos seus pontos também sobrepostos sobre a reta y=x, exceto para tempos maiores, mas ainda sim proximos da reta, e por fim temos o modelo Log-Normal, aonde o comportamento não aparentou seguir uma linearidade, evidenciando que o modelo não ajustou-se bem a esse conjunto de dados.

Para o teste TRV, mesmo sendo impreciso para esta situação, devido a grande presença de censuras, obtivemos que para nenhum dos 3 modelos testados apresentados tiveram um pvalor superior ao nivel de significancia de 5%, ou seja, temos que pelo TRV nenhum dos 3 modelos pode ser considerado adequados.

Por fim, iremos considerar o modelo weibull por apresentar a sua curva estimada da sobrevivencia mais proximo da do kaplan-meyer e apresentar grande parte dos seus pontos, no metodo 02, dispostos sobre a reta y=x, para continuarmos as nossas analises.

Tempo medio

tmw<-as.vector((gamma(1+(1/gama)))/lambdaw)
c('dias' = tmw,'anos' = tmw/365)

##        dias        anos 
## 1103.525044    3.023356

Temos que o tempo médio para a sobrevivência de um processo juridico resultar em perda para a empresa é de aproximadamente 1104 dias, aproximadamente 3 anos.

Tempo mediano

p = 0.5

tpw<-as.vector((-log(1-p))^(1/gama)/lambdaw)
c('dias' = tpw,'anos' = tpw/365)

##       dias       anos 
## 666.085316   1.824891

Temos que o tempo mediano para a sobrevivência de um processo juridico resultar em perda para a empresa é de 666 dias, aproximadamente 1 ano e 10 meses.

Percentil 75%

p = 0.75

tpw<-as.vector((-log(1-p))^(1/gama)/lambdaw)
c('dias' = tpw,'anos' = tpw/365)

##        dias        anos 
## 1493.130187    4.090768

Temos que, com uma probabilidade de 75%, o tempo necessário de sobrevivência do processo juridico ate a sua perda para esta probabilidade é de aproximadamente 1493 dias, equivalente a 4 anos e 1 mes, aproximadamente.