Análise de Sobrevivencia
Desgaste de funcionários usando análise de sobrevivência.
Sobre o conjunto de dados
Banco de dados sobre Atrito e desempenho de funcionários do IBM HR Analytics. Possui 35 variáveis, sendo 9 do tipo categórico e 26 do tipo numérico, com 1470 linhas. Descrição das variáveis
Age - Idade
Attrition - Atrito, variável resposta
BusinessTravel - Viagem a trabalho
DailyRate - Avaliação diária
Department - Departamento
DistanceFromHome - Distância de casa
Education - Educação
EducationField - Campo de Educação
EmployeeCount - Contagem de funcionários
EnvironmentSatisfaction - Nivel de satisfação
Gender - Gênero
HourlyRate - Taxa por hora
JobInvolvement - Envolvimento com o trabalho
JobLevel - Nível de trabalho
JobRole - Cargo de trabalho
JobSatisfaction - Satisfação com o trabalho
MaritalStatus - Estado civil
MonthlyIncome - Renda mensal
MonthlyRate - Taxa mensal
NumCompaniesWorked - Número de empresas trabalhadas
Over18 - Mais de 18 anos
OverTime - Hora extra
PercentSalaryHike - Percentual de aumento salarial
PerformanceRating - Classificação de desempenho
RelationshipSatisfaction - satisfação de relacionamento
StandardHours - horario padrão
StockOptionLevel -nivel de opção de ações
TotalWorkingYears - Total de anos trabalhados
TrainingTimesLastYear - Tempo de treinamento no ultimo ano
WorkLifeBalance - Equilibrio da vida profissional
YearsAtCompany - Anos na Empresa
YearsInCurrentRole - Anos na função atual
YearsSinceLastPromotion - Anos desde a ultima promoção
YearsWithCurrManager - Anos com o gerente atual
library(readr)
atrito <- read_csv("/Users/Clevia/Documents/Analise de Sobrevivencia/WA_Fn-UseC_-HR-Employee-Attrition.csv")Análise descritiva das variáveis
library(dplyr) #graficos
require(ggplot2)#graficos
library(summarytools)# análise descritiva
library(tidyverse)
library(survival) #analise de sobrevivencia
library(ggfortify) #gráficooptions(warn=-1) #suppress warnings
print(dfSummary(atrito,style = 'grid', valid.col = FALSE),method = 'render')Data Frame Summary
atrito
Dimensions: 1470 x 35Duplicates: 0
| No | Variable | Stats / Values | Freqs (% of Valid) | Graph | Missing | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | Age [numeric] |
|
43 distinct values |  |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 2 | Attrition [character] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 3 | BusinessTravel [character] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 4 | DailyRate [numeric] |
|
886 distinct values |  |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 5 | Department [character] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 6 | DistanceFromHome [numeric] |
|
29 distinct values |  |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 7 | Education [numeric] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 8 | EducationField [character] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 9 | EmployeeCount [numeric] | 1 distinct value |
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 10 | EmployeeNumber [numeric] |
|
1470 distinct values |  |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 11 | EnvironmentSatisfaction [numeric] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 12 | Gender [character] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 13 | HourlyRate [numeric] |
|
71 distinct values |  |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 14 | JobInvolvement [numeric] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 15 | JobLevel [numeric] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 16 | JobRole [character] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 17 | JobSatisfaction [numeric] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 18 | MaritalStatus [character] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 19 | MonthlyIncome [numeric] |
|
1349 distinct values |  |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 20 | MonthlyRate [numeric] |
|
1427 distinct values |  |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 21 | NumCompaniesWorked [numeric] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 22 | Over18 [character] | 1. Y |
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 23 | OverTime [character] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 24 | PercentSalaryHike [numeric] |
|
15 distinct values |  |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 25 | PerformanceRating [numeric] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 26 | RelationshipSatisfaction [numeric] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 27 | StandardHours [numeric] | 1 distinct value |
|
|
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 28 | StockOptionLevel [numeric] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 29 | TotalWorkingYears [numeric] |
|
40 distinct values |  |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 30 | TrainingTimesLastYear [numeric] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 31 | WorkLifeBalance [numeric] |
|
|
 |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 32 | YearsAtCompany [numeric] |
|
37 distinct values |  |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 33 | YearsInCurrentRole [numeric] |
|
19 distinct values |  |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 34 | YearsSinceLastPromotion [numeric] |
|
16 distinct values |  |
0 (0.0%) | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| 35 | YearsWithCurrManager [numeric] |
|
18 distinct values |  |
0 (0.0%) |
Generated by summarytools 1.0.1 (R version 4.1.3)
2022-12-06
Proporção da variável resposta no banco de dados
Aproximadamente 84% não há atrito e aproximadamente 16% há atrito.
atrito$YearsAtCompany[atrito$YearsAtCompany %in% c("0")]<- "1"
#table(atrito$YearsAtCompany)
#levels(atrito$YearsAtCompany)number_of_attirttions <- atrito %>% group_by(Attrition) %>% summarise(Count=n())
number_of_attirttions## # A tibble: 2 x 2
## Attrition Count
## <chr> <int>
## 1 No 1233
## 2 Yes 237#visualization
ggplot(data=number_of_attirttions, aes(x=Attrition, y=Count)) +
geom_histogram(stat = "identity",fill="blue")
Atrito por hora extra
attrition_OverT= atrito %>% dplyr::select(Attrition, OverTime) %>% filter ( Attrition =="Yes" ) %>%group_by(OverTime)%>%summarise(Count=n())
attrition_OverT## # A tibble: 2 x 2
## OverTime Count
## <chr> <int>
## 1 No 110
## 2 Yes 127ggplot(data=attrition_OverT, aes(x=OverTime, y=Count ,fill=OverTime))+
geom_col()
Vemos que funcionário que fazem hora extra tem mais atrito do que funcionários que não fazem.
Atrito por departamento
Aqui observamos o atrito entre funcionários por departamento. Vemos que no departamento de Pesquisa de Desenvolvimento há um maior número de atrito.
attach(atrito)
dapartment_attrition <- atrito %>%
dplyr::select(Department, Attrition) %>%
filter(Attrition== "Yes") %>%
group_by(Department) %>% summarise(Count=n())
#jobe_role <-attrition_employee %>% select(JobRole,Attrition) %>% filter(Attrition=='Yes') %>% group_by(JobRole)%>% summarise(Count=n())
dapartment_attrition## # A tibble: 3 x 2
## Department Count
## <chr> <int>
## 1 Human Resources 12
## 2 Research & Development 133
## 3 Sales 92#visualisation of the data
ggplot(data=dapartment_attrition, aes(x=Department, y=Count , fill=Department))+
geom_col()
Atrito por gênero
Observamos que há um maior número de atrito entre os funcionários do gênero masculino.
attrition_gender= atrito %>% dplyr::select(Attrition, Gender) %>% filter ( Attrition =="Yes" ) %>%group_by(Gender)%>%summarise(Count=n())
attrition_gender## # A tibble: 2 x 2
## Gender Count
## <chr> <int>
## 1 Female 87
## 2 Male 150ggplot(data=attrition_gender, aes(x=Gender, y=Count ,fill=Gender))+
geom_col()
Atrito por viagem a trabalho
attrition_businessT = atrito %>% dplyr::select(Attrition, BusinessTravel) %>% filter ( Attrition =="Yes" ) %>%group_by(BusinessTravel)%>%summarise(Count=n())
attrition_businessT## # A tibble: 3 x 2
## BusinessTravel Count
## <chr> <int>
## 1 Non-Travel 12
## 2 Travel_Frequently 69
## 3 Travel_Rarely 156ggplot(data=attrition_businessT, aes(x=BusinessTravel, y=Count ,fill=BusinessTravel))+
geom_col()
Vemos que há um número maior de atrito entre funcionários que raramente viajam a trabalho e um número menor de atrito entre funcionários que não viajam.
Atrito por função
jobrole_attrition= atrito %>% dplyr::select( Attrition,Department ,MonthlyIncome , JobRole) %>% filter (Attrition=="Yes") %>% group_by(JobRole) %>% summarise(Count=n())
jobrole_attrition## # A tibble: 9 x 2
## JobRole Count
## <chr> <int>
## 1 Healthcare Representative 9
## 2 Human Resources 12
## 3 Laboratory Technician 62
## 4 Manager 5
## 5 Manufacturing Director 10
## 6 Research Director 2
## 7 Research Scientist 47
## 8 Sales Executive 57
## 9 Sales Representative 33ggplot(data=jobrole_attrition, aes(x=JobRole, y=Count , fill= JobRole ))+
geom_col() +
facet_wrap(~ JobRole)
Vemos que técnico de laboratório, executivo de vendas, cientista de pesquisa e representante de vendas têm a maior taxa de atrito entre as diferentes funções.
#str(atrito)
atrito$YearsAtCompany <- as.numeric(atrito$YearsAtCompany)
#str(atrito)atrito$event <- with(atrito,ifelse(Attrition=="Yes", 1, 0))
time<-atrito$YearsAtCompany #define tempo
event<-atrito$event #define evento
#group<-atrito$OverTime #define groups to compare multiple survival curves
survival<-Surv(time,event) #criação de survival object
head(survival, n = 15)## [1] 6 10+ 1 8+ 2+ 7+ 1+ 1+ 9+ 7+ 5+ 9+ 5+ 2+ 4#survivalCriando um modelo de sobrevivência
model<- survfit(Surv(time,event) ~ 1)
print(model, print.rmean=TRUE)## Call: survfit(formula = Surv(time, event) ~ 1)
##
## n events rmean* se(rmean) median 0.95LCL 0.95UCL
## [1,] 1470 237 28.6 0.94 40 32 NA
## * restricted mean with upper limit = 40summary(model)## Call: survfit(formula = Surv(time, event) ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1470 75 0.949 0.00574 0.938 0.960
## 2 1255 27 0.929 0.00683 0.915 0.942
## 3 1128 20 0.912 0.00764 0.897 0.927
## 4 1000 19 0.895 0.00846 0.878 0.912
## 5 890 21 0.874 0.00943 0.855 0.892
## 6 694 9 0.862 0.01004 0.843 0.882
## 7 618 11 0.847 0.01088 0.826 0.869
## 8 528 9 0.833 0.01171 0.810 0.856
## 9 448 8 0.818 0.01262 0.793 0.843
## 10 366 18 0.777 0.01515 0.748 0.808
## 11 246 2 0.771 0.01567 0.741 0.802
## 13 200 2 0.763 0.01643 0.732 0.796
## 14 176 2 0.755 0.01736 0.721 0.790
## 15 158 1 0.750 0.01789 0.716 0.786
## 16 138 1 0.745 0.01857 0.709 0.782
## 17 126 1 0.739 0.01934 0.702 0.778
## 18 117 1 0.732 0.02018 0.694 0.773
## 19 104 1 0.725 0.02118 0.685 0.768
## 20 93 1 0.717 0.02234 0.675 0.763
## 21 66 1 0.707 0.02450 0.660 0.756
## 22 52 1 0.693 0.02754 0.641 0.749
## 23 37 1 0.674 0.03255 0.613 0.741
## 24 35 1 0.655 0.03688 0.587 0.731
## 31 16 1 0.614 0.05260 0.519 0.726
## 32 13 1 0.567 0.06646 0.450 0.713
## 33 10 1 0.510 0.08044 0.375 0.695
## 40 1 1 0.000 NaN NA NASaída do modelo,
tempo = os pontos de tempo em que a curva tem um passo
n.risco = número de funcionários em risco de desligamento no tempo
t n.event = o número de eventos que ocorrem no tempo t
survival = estimativa de Kaplan-Meier (probabilidade de sobrevivência) no tempo t
std.err = erro padrão da estimativa K-M no tempo t
lower IC 95% = limite de confiança inferior para a curva de sobrevivência ou probabilidade
IC 95% upper = limite de confiança superior para a curva ou probabilidade de sobrevivência
O método padrão para a transformação para construção do intervalo de confiança é “log”
Plot do modelo
autoplot(model)+ labs(x="Tempo em Anos",y="Sobrevivência",title="Gráfico de Sobrevivência")
Estimador do risco acumulados de Kaplan-Meier
autoplot(model, conf.int = T, surv.size= 1, fun = "cumhaz")+ labs(x="Tempo de Sobrevivência",y="Risco",title="Estimador de Risco Acumulado")
Comparação entre curvas de sobrevivência
Hora Extra
#comparação entre curvas
model_comp<- survfit(Surv(time,event) ~ atrito$OverTime, data = atrito)
summary(model_comp)## Call: survfit(formula = Surv(time, event) ~ atrito$OverTime, data = atrito)
##
## atrito$OverTime=No
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1054 37 0.965 0.00567 0.954 0.976
## 2 908 12 0.952 0.00668 0.939 0.965
## 3 816 6 0.945 0.00722 0.931 0.959
## 4 728 11 0.931 0.00830 0.915 0.947
## 5 648 7 0.921 0.00903 0.903 0.939
## 6 509 4 0.914 0.00966 0.895 0.933
## 7 457 3 0.908 0.01020 0.888 0.928
## 8 388 3 0.901 0.01090 0.879 0.922
## 9 332 6 0.884 0.01256 0.860 0.909
## 10 269 8 0.858 0.01525 0.829 0.888
## 11 181 1 0.853 0.01588 0.823 0.885
## 13 146 2 0.842 0.01768 0.808 0.877
## 14 126 2 0.828 0.01977 0.790 0.868
## 17 89 1 0.819 0.02162 0.778 0.862
## 18 82 1 0.809 0.02355 0.764 0.856
## 20 65 1 0.796 0.02627 0.747 0.850
## 22 35 1 0.774 0.03398 0.710 0.843
## 23 24 1 0.741 0.04535 0.658 0.836
## 32 8 1 0.649 0.09535 0.486 0.865
## 33 6 1 0.541 0.12671 0.342 0.856
## 40 1 1 0.000 NaN NA NA
##
## atrito$OverTime=Yes
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 416 38 0.909 0.0141 0.881 0.937
## 2 347 15 0.869 0.0168 0.837 0.903
## 3 312 14 0.830 0.0190 0.794 0.868
## 4 272 8 0.806 0.0203 0.767 0.847
## 5 242 14 0.759 0.0226 0.716 0.805
## 6 185 5 0.739 0.0238 0.694 0.787
## 7 161 8 0.702 0.0259 0.653 0.755
## 8 140 6 0.672 0.0276 0.620 0.728
## 9 116 2 0.660 0.0283 0.607 0.718
## 10 97 10 0.592 0.0325 0.532 0.660
## 11 65 1 0.583 0.0333 0.521 0.652
## 15 50 1 0.572 0.0346 0.508 0.644
## 16 43 1 0.558 0.0363 0.492 0.634
## 19 32 1 0.541 0.0391 0.469 0.623
## 21 22 1 0.516 0.0444 0.436 0.611
## 24 13 1 0.477 0.0560 0.379 0.600
## 31 7 1 0.408 0.0792 0.279 0.597autoplot(model_comp, conf.int = F, surv.size= 1)+ labs(x="Tempo de Sobrevivência em anos",y="Probabilidade de Sobrevivência",title="Hora Extra",colour="Hora Extra")
Teste de Peto
survdiff(survival ~ atrito$OverTime, rho=1)## Call:
## survdiff(formula = survival ~ atrito$OverTime, rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## atrito$OverTime=No 1054 99.6 156.1 20.5 83
## atrito$OverTime=Yes 416 116.3 59.8 53.4 83
##
## Chisq= 83 on 1 degrees of freedom, p= <2e-16Como o p-valor é menor que 0,05 podemos concluir pelo teste de Peto que há diferença entre os extratos
Viagem á trabalho
#comparação entre curvas
model_comp<- survfit(Surv(time,event) ~ atrito$BusinessTravel, data = atrito)
summary(model_comp)## Call: survfit(formula = Surv(time, event) ~ atrito$BusinessTravel,
## data = atrito)
##
## atrito$BusinessTravel=Non-Travel
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 150 5 0.967 0.0147 0.938 0.996
## 2 125 1 0.959 0.0165 0.927 0.992
## 3 119 1 0.951 0.0182 0.916 0.987
## 4 103 1 0.942 0.0202 0.903 0.982
## 5 94 2 0.922 0.0242 0.875 0.970
## 10 38 2 0.873 0.0405 0.797 0.956
##
## atrito$BusinessTravel=Travel_Frequently
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 277 24 0.913 0.0169 0.881 0.947
## 2 241 7 0.887 0.0192 0.850 0.925
## 3 211 6 0.862 0.0212 0.821 0.904
## 4 191 6 0.835 0.0232 0.790 0.881
## 5 168 4 0.815 0.0247 0.768 0.865
## 6 142 4 0.792 0.0265 0.741 0.846
## 7 123 3 0.772 0.0281 0.719 0.830
## 8 107 3 0.751 0.0300 0.694 0.812
## 9 90 2 0.734 0.0316 0.675 0.799
## 10 71 5 0.682 0.0369 0.614 0.759
## 11 46 1 0.668 0.0389 0.595 0.748
## 13 39 1 0.650 0.0415 0.574 0.737
## 18 26 1 0.625 0.0469 0.540 0.724
## 21 15 1 0.584 0.0595 0.478 0.713
## 23 8 1 0.511 0.0858 0.367 0.710
##
## atrito$BusinessTravel=Travel_Rarely
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1043 46 0.956 0.00636 0.944 0.968
## 2 889 19 0.935 0.00776 0.920 0.951
## 3 798 13 0.920 0.00871 0.903 0.937
## 4 706 12 0.905 0.00966 0.886 0.924
## 5 628 15 0.883 0.01092 0.862 0.905
## 6 479 5 0.874 0.01156 0.851 0.897
## 7 429 8 0.857 0.01270 0.833 0.883
## 8 361 6 0.843 0.01376 0.817 0.871
## 9 306 6 0.827 0.01505 0.798 0.857
## 10 257 11 0.791 0.01779 0.757 0.827
## 11 171 1 0.787 0.01828 0.752 0.823
## 13 139 1 0.781 0.01900 0.745 0.819
## 14 126 2 0.769 0.02062 0.729 0.810
## 15 112 1 0.762 0.02155 0.721 0.805
## 16 100 1 0.754 0.02264 0.711 0.800
## 17 89 1 0.746 0.02392 0.700 0.794
## 19 73 1 0.735 0.02568 0.687 0.788
## 20 66 1 0.724 0.02760 0.672 0.780
## 22 37 1 0.705 0.03308 0.643 0.773
## 24 25 1 0.677 0.04209 0.599 0.764
## 31 11 1 0.615 0.07002 0.492 0.769
## 32 9 1 0.547 0.08958 0.397 0.754
## 33 8 1 0.478 0.10114 0.316 0.724
## 40 1 1 0.000 NaN NA NAautoplot(model_comp, conf.int = F, surv.size= 1)+ labs(x="Tempo de Sobrevivência em anos",y="Probabilidade de Sobrevivência",title="Viagem a trabalho",colour="Viagem a Trabalho")
Teste de Peto
survdiff(survival ~ atrito$BusinessTravel, rho=1)## Call:
## survdiff(formula = survival ~ atrito$BusinessTravel, rho = 1)
##
## N Observed Expected (O-E)^2/E
## atrito$BusinessTravel=Non-Travel 150 11.2 22.4 5.599
## atrito$BusinessTravel=Travel_Frequently 277 63.3 41.2 11.934
## atrito$BusinessTravel=Travel_Rarely 1043 141.3 152.3 0.789
## (O-E)^2/V
## atrito$BusinessTravel=Non-Travel 7.03
## atrito$BusinessTravel=Travel_Frequently 16.56
## atrito$BusinessTravel=Travel_Rarely 3.01
##
## Chisq= 20.6 on 2 degrees of freedom, p= 3e-05Como o p-valor é menor que 0,05 podemos concluir pelo teste de Peto que há diferença entre os extratos
Gênero
model_comp<- survfit(Surv(time,event) ~ atrito$Gender, data= atrito)
summary(model_comp)## Call: survfit(formula = Surv(time, event) ~ atrito$Gender, data = atrito)
##
## atrito$Gender=Female
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 588 25 0.957 0.00832 0.941 0.974
## 2 509 10 0.939 0.01006 0.919 0.959
## 3 465 7 0.925 0.01124 0.903 0.947
## 4 409 5 0.913 0.01218 0.890 0.937
## 5 374 9 0.891 0.01392 0.864 0.919
## 6 287 3 0.882 0.01478 0.853 0.911
## 7 258 5 0.865 0.01635 0.833 0.897
## 8 222 4 0.849 0.01781 0.815 0.885
## 9 185 4 0.831 0.01965 0.793 0.870
## 10 149 6 0.797 0.02313 0.753 0.844
## 11 113 1 0.790 0.02397 0.745 0.839
## 13 92 1 0.782 0.02521 0.734 0.833
## 15 72 1 0.771 0.02709 0.720 0.826
## 17 56 1 0.757 0.02990 0.701 0.818
## 22 18 1 0.715 0.04969 0.624 0.819
## 24 11 1 0.650 0.07670 0.516 0.819
## 32 3 1 0.433 0.18417 0.188 0.997
## 33 2 1 0.217 0.17877 0.043 1.000
## 40 1 1 0.000 NaN NA NA
##
## atrito$Gender=Male
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 882 50 0.943 0.00779 0.928 0.959
## 2 746 17 0.922 0.00919 0.904 0.940
## 3 663 13 0.904 0.01029 0.884 0.924
## 4 591 14 0.882 0.01153 0.860 0.905
## 5 516 12 0.862 0.01269 0.837 0.887
## 6 407 6 0.849 0.01352 0.823 0.876
## 7 360 6 0.835 0.01448 0.807 0.864
## 8 306 5 0.821 0.01547 0.792 0.852
## 9 263 4 0.809 0.01645 0.777 0.842
## 10 217 12 0.764 0.01997 0.726 0.804
## 11 133 1 0.758 0.02063 0.719 0.800
## 13 108 1 0.751 0.02160 0.710 0.795
## 14 94 2 0.735 0.02392 0.690 0.784
## 16 76 1 0.726 0.02549 0.677 0.777
## 18 66 1 0.715 0.02737 0.663 0.770
## 19 60 1 0.703 0.02939 0.647 0.763
## 20 55 1 0.690 0.03151 0.631 0.755
## 21 42 1 0.674 0.03478 0.609 0.745
## 23 26 1 0.648 0.04200 0.570 0.735
## 31 13 1 0.598 0.06159 0.489 0.732autoplot(model_comp, conf.int = F, surv.size= 1)+ labs(x="Tempo de Sobrevivência",y="Probabilidade de Sobrevivência",title="Gênero",colour="Genero")
Teste de Peto
survdiff(survival ~ atrito$Gender, rho=1)## Call:
## survdiff(formula = survival ~ atrito$Gender, rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## atrito$Gender=Female 588 78.3 88.1 1.090 2.08
## atrito$Gender=Male 882 137.6 127.8 0.751 2.08
##
## Chisq= 2.1 on 1 degrees of freedom, p= 0.1Como o p-valor é maior que 0,05 podemos concluir pelo teste de Peto que não há diferença significativa entre os extratos
Departamento
model_comp<- survfit(Surv(time,event) ~ Department, data= atrito)
summary(model_comp)## Call: survfit(formula = Surv(time, event) ~ Department, data = atrito)
##
## Department=Human Resources
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 63 4 0.937 0.0307 0.878 0.999
## 2 58 2 0.904 0.0372 0.834 0.980
## 3 50 2 0.868 0.0436 0.787 0.958
## 4 42 1 0.847 0.0472 0.760 0.945
## 5 38 1 0.825 0.0510 0.731 0.931
## 7 24 1 0.791 0.0593 0.683 0.916
## 20 7 1 0.678 0.1163 0.484 0.949
##
## Department=Research & Development
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 961 47 0.951 0.00696 0.938 0.965
## 2 812 13 0.936 0.00803 0.920 0.952
## 3 735 8 0.926 0.00871 0.909 0.943
## 4 650 11 0.910 0.00976 0.891 0.929
## 5 571 13 0.889 0.01110 0.868 0.911
## 6 441 4 0.881 0.01171 0.859 0.904
## 7 389 8 0.863 0.01310 0.838 0.889
## 8 332 4 0.853 0.01394 0.826 0.880
## 9 280 4 0.841 0.01501 0.812 0.870
## 10 230 14 0.789 0.01935 0.752 0.828
## 13 127 1 0.783 0.02017 0.745 0.824
## 17 81 1 0.773 0.02212 0.731 0.818
## 18 74 1 0.763 0.02416 0.717 0.812
## 22 32 1 0.739 0.03315 0.677 0.807
## 31 9 1 0.657 0.08285 0.513 0.841
## 33 5 1 0.526 0.13494 0.318 0.869
## 40 1 1 0.000 NaN NA NA
##
## Department=Sales
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 446 24 0.946 0.0107 0.925 0.967
## 2 385 12 0.917 0.0133 0.891 0.943
## 3 343 10 0.890 0.0154 0.860 0.921
## 4 308 7 0.870 0.0168 0.837 0.903
## 5 281 7 0.848 0.0183 0.813 0.885
## 6 226 5 0.829 0.0197 0.792 0.869
## 7 205 2 0.821 0.0203 0.782 0.862
## 8 176 5 0.798 0.0223 0.755 0.843
## 9 150 4 0.777 0.0241 0.731 0.825
## 10 121 4 0.751 0.0265 0.701 0.805
## 11 82 2 0.733 0.0288 0.678 0.791
## 13 66 1 0.722 0.0305 0.664 0.784
## 14 57 2 0.696 0.0343 0.632 0.767
## 15 51 1 0.683 0.0362 0.615 0.757
## 16 41 1 0.666 0.0390 0.594 0.747
## 19 31 1 0.644 0.0432 0.565 0.735
## 21 21 1 0.614 0.0509 0.522 0.722
## 23 13 1 0.567 0.0653 0.452 0.710
## 24 12 1 0.519 0.0750 0.391 0.689
## 32 5 1 0.415 0.1106 0.247 0.700autoplot(model_comp, conf.int = F, surv.size= 1)+ labs(x="Tempo de Sobrevivência",y="Probabilidade de Sobrevivência",title="Departamentos",colour="Departamentos")
Teste de Peto
survdiff(survival ~ atrito$Department, rho=1)## Call:
## survdiff(formula = survival ~ atrito$Department, rho = 1)
##
## N Observed Expected (O-E)^2/E
## atrito$Department=Human Resources 63 11.1 9.42 0.317
## atrito$Department=Research & Development 961 121.7 139.10 2.181
## atrito$Department=Sales 446 83.0 67.35 3.656
## (O-E)^2/V
## atrito$Department=Human Resources 0.374
## atrito$Department=Research & Development 6.905
## atrito$Department=Sales 5.980
##
## Chisq= 6.9 on 2 degrees of freedom, p= 0.03Como o p-valor é menor que 0,05 podemos concluir pelo teste de Peto que há diferença significativa entre os extratos
Cargo de Trabalho
model_comp<- survfit(Surv(time,event) ~ JobRole, data= atrito)
summary(model_comp)## Call: survfit(formula = Surv(time, event) ~ JobRole, data = atrito)
##
## JobRole=Healthcare Representative
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 131 3 0.977 0.0131 0.952 1.000
## 8 60 1 0.961 0.0206 0.921 1.000
## 9 48 1 0.941 0.0283 0.887 0.998
## 10 41 2 0.895 0.0415 0.817 0.980
## 18 13 1 0.826 0.0765 0.689 0.990
## 40 1 1 0.000 NaN NA NA
##
## JobRole=Human Resources
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 52 4 0.923 0.0370 0.853 0.998
## 2 47 2 0.884 0.0446 0.801 0.976
## 3 39 2 0.838 0.0526 0.741 0.948
## 4 32 1 0.812 0.0571 0.708 0.932
## 5 28 1 0.783 0.0620 0.671 0.915
## 7 16 1 0.734 0.0750 0.601 0.897
## 20 1 1 0.000 NaN NA NA
##
## JobRole=Laboratory Technician
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 259 28 0.892 0.0193 0.855 0.931
## 2 198 7 0.860 0.0220 0.818 0.905
## 3 171 4 0.840 0.0237 0.795 0.888
## 4 145 7 0.800 0.0270 0.748 0.854
## 5 123 5 0.767 0.0296 0.711 0.827
## 6 86 1 0.758 0.0306 0.701 0.821
## 7 77 5 0.709 0.0356 0.642 0.782
## 9 50 2 0.681 0.0395 0.608 0.763
## 10 37 2 0.644 0.0451 0.561 0.739
## 17 7 1 0.552 0.0935 0.396 0.769
##
## JobRole=Manager
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 5 85 1 0.988 0.0117 0.966 1
## 7 72 1 0.975 0.0179 0.940 1
## 10 60 1 0.958 0.0238 0.913 1
## 24 21 1 0.913 0.0500 0.820 1
## 32 8 1 0.799 0.1153 0.602 1
##
## JobRole=Manufacturing Director
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 145 2 0.986 0.00969 0.967 1.000
## 3 127 1 0.978 0.01234 0.955 1.000
## 4 114 1 0.970 0.01492 0.941 1.000
## 5 103 1 0.960 0.01749 0.927 0.995
## 10 42 4 0.869 0.04629 0.783 0.965
## 33 1 1 0.000 NaN NA NA
##
## JobRole=Research Director
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 22 9 1 0.889 0.105 0.706 1
## 31 3 1 0.593 0.252 0.258 1
##
## JobRole=Research Scientist
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 292 14 0.952 0.0125 0.928 0.977
## 2 238 6 0.928 0.0156 0.898 0.959
## 3 207 3 0.915 0.0172 0.882 0.949
## 4 176 3 0.899 0.0191 0.862 0.937
## 5 146 6 0.862 0.0235 0.817 0.909
## 6 100 3 0.836 0.0271 0.785 0.891
## 7 78 2 0.815 0.0304 0.757 0.877
## 8 65 3 0.777 0.0359 0.710 0.851
## 9 51 1 0.762 0.0383 0.690 0.841
## 10 42 5 0.671 0.0509 0.579 0.779
## 13 16 1 0.629 0.0626 0.518 0.765
##
## JobRole=Sales Executive
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 326 7 0.979 0.00803 0.963 0.994
## 2 301 8 0.953 0.01197 0.929 0.976
## 3 278 3 0.942 0.01323 0.917 0.969
## 4 258 5 0.924 0.01529 0.894 0.954
## 5 237 5 0.904 0.01728 0.871 0.939
## 6 187 5 0.880 0.01991 0.842 0.920
## 7 167 2 0.870 0.02102 0.830 0.912
## 8 143 5 0.839 0.02429 0.793 0.888
## 9 119 3 0.818 0.02658 0.768 0.872
## 10 95 4 0.784 0.03053 0.726 0.846
## 11 58 2 0.757 0.03495 0.691 0.828
## 13 45 1 0.740 0.03801 0.669 0.818
## 14 38 2 0.701 0.04489 0.618 0.795
## 15 32 1 0.679 0.04853 0.590 0.781
## 16 22 1 0.648 0.05528 0.548 0.766
## 19 16 1 0.608 0.06499 0.493 0.749
## 21 8 1 0.532 0.09101 0.380 0.744
## 23 3 1 0.354 0.15692 0.149 0.844
##
## JobRole=Sales Representative
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 83 17 0.795 0.0443 0.713 0.887
## 2 53 4 0.735 0.0501 0.643 0.840
## 3 35 7 0.588 0.0639 0.475 0.728
## 4 21 2 0.532 0.0690 0.413 0.686
## 5 15 2 0.461 0.0759 0.334 0.637
## 9 6 1 0.384 0.0944 0.237 0.622autoplot(model_comp, conf.int = F, surv.size= 1)+ labs(x="Tempo de Sobrevivência",y="Probabilidade de Sobrevivência",title="Cargo de Trabalho",colour="Cargo")
Teste de Peto
survdiff(survival ~ atrito$JobRole, rho=1)## Call:
## survdiff(formula = survival ~ atrito$JobRole, rho = 1)
##
## N Observed Expected (O-E)^2/E
## atrito$JobRole=Healthcare Representative 131 7.56 22.14 9.5991
## atrito$JobRole=Human Resources 52 11.15 6.61 3.1235
## atrito$JobRole=Laboratory Technician 259 58.45 31.09 24.0642
## atrito$JobRole=Manager 102 3.86 23.52 16.4309
## atrito$JobRole=Manufacturing Director 145 8.57 22.99 9.0391
## atrito$JobRole=Research Director 80 1.36 15.40 12.7948
## atrito$JobRole=Research Scientist 292 43.16 35.62 1.5968
## atrito$JobRole=Sales Executive 326 50.01 51.33 0.0338
## atrito$JobRole=Sales Representative 83 31.74 7.17 84.1631
## (O-E)^2/V
## atrito$JobRole=Healthcare Representative 12.35
## atrito$JobRole=Human Resources 3.56
## atrito$JobRole=Laboratory Technician 31.35
## atrito$JobRole=Manager 23.19
## atrito$JobRole=Manufacturing Director 11.37
## atrito$JobRole=Research Director 15.97
## atrito$JobRole=Research Scientist 2.14
## atrito$JobRole=Sales Executive 0.05
## atrito$JobRole=Sales Representative 95.06
##
## Chisq= 187 on 8 degrees of freedom, p= <2e-16Como o p-valor é menor que 0,05 podemos concluir pelo teste de Peto que há diferença significativa entre os extratos
ajust1<-survreg(survival~1,dist="exponential")
alpha1<-exp(ajust1$coefficients[1])ajust2<-survreg(survival~1,dist="weibull")
alpha2<-exp(ajust2$coefficients[1]);
gama<-1/ajust2$scale
ajust3<-survreg(survival~1,dist="lognorm")
#S(T) KAPLA VS ESTIMATIVAS
ekm<- survfit(survival~1, data= atrito)
time<-ekm$time
st<-ekm$surv
ste<- exp(-time/alpha1)
stw<- exp(-(time/alpha2)^gama)
stln<- pnorm((-log(time)+ 4.778)/2.0347)
cbind(time,st,ste,stw,stln)## time st ste stw stln
## [1,] 1 0.9489796 0.9773530 0.9736467 0.9905693
## [2,] 2 0.9285633 0.9552188 0.9503263 0.9776568
## [3,] 3 0.9120994 0.9335860 0.9283537 0.9647216
## [4,] 4 0.8947695 0.9124430 0.9073750 0.9522360
## [5,] 5 0.8736570 0.8917789 0.8872166 0.9402959
## [6,] 6 0.8623271 0.8715828 0.8677717 0.9289012
## [7,] 7 0.8469783 0.8518440 0.8489663 0.9180230
## [8,] 8 0.8325411 0.8325523 0.8307454 0.9076246
## [9,] 9 0.8176743 0.8136975 0.8130655 0.8976691
## [10,] 10 0.7774608 0.7952697 0.7958912 0.8881217
## [11,] 11 0.7711400 0.7772592 0.7791929 0.8789506
## [12,] 12 0.7711400 0.7596566 0.7629449 0.8701274
## [13,] 13 0.7634286 0.7424526 0.7471250 0.8616262
## [14,] 14 0.7547533 0.7256383 0.7317133 0.8534241
## [15,] 15 0.7499764 0.7092047 0.7166921 0.8455004
## [16,] 16 0.7445418 0.6931433 0.7020453 0.8378366
## [17,] 17 0.7386327 0.6774457 0.6877582 0.8304159
## [18,] 18 0.7323196 0.6621036 0.6738173 0.8232231
## [19,] 19 0.7252781 0.6471089 0.6602101 0.8162446
## [20,] 20 0.7174794 0.6324538 0.6469250 0.8094677
## [21,] 21 0.7066085 0.6181306 0.6339510 0.8028811
## [22,] 22 0.6930199 0.6041318 0.6212780 0.7964744
## [23,] 23 0.6742896 0.5904500 0.6088963 0.7902379
## [24,] 24 0.6550242 0.5770781 0.5967970 0.7841630
## [25,] 25 0.6550242 0.5640090 0.5849714 0.7782414
## [26,] 26 0.6550242 0.5512359 0.5734113 0.7724656
## [27,] 27 0.6550242 0.5387520 0.5621091 0.7668288
## [28,] 29 0.6550242 0.5146261 0.5402488 0.7559465
## [29,] 30 0.6550242 0.5029713 0.5296770 0.7506895
## [30,] 31 0.6140852 0.4915805 0.5193354 0.7455482
## [31,] 32 0.5668478 0.4804477 0.5092177 0.7405177
## [32,] 33 0.5101631 0.4695670 0.4993180 0.7355934
## [33,] 34 0.5101631 0.4589327 0.4896306 0.7307712
## [34,] 36 0.5101631 0.4383812 0.4708708 0.7214167
## [35,] 37 0.5101631 0.4284531 0.4617880 0.7168771
## [36,] 40 0.0000000 0.3999978 0.4356694 0.7037696Ajuste de modelos
ajuste_ex <-survreg(survival~ OverTime + BusinessTravel + Department + JobRole, data = atrito, dist="exponential")
summary(ajuste_ex)##
## Call:
## survreg(formula = survival ~ OverTime + BusinessTravel + Department +
## JobRole, data = atrito, dist = "exponential")
## Value Std. Error z p
## (Intercept) 21.036 2086.276 0.01 0.99196
## OverTimeYes -1.099 0.131 -8.37 < 2e-16
## BusinessTravelTravel_Frequently -0.904 0.314 -2.88 0.00401
## BusinessTravelTravel_Rarely -0.553 0.301 -1.84 0.06591
## DepartmentResearch & Development -15.155 2086.276 -0.01 0.99420
## DepartmentSales -14.991 2086.276 -0.01 0.99427
## JobRoleHuman Resources -16.852 2086.276 -0.01 0.99355
## JobRoleLaboratory Technician -1.877 0.357 -5.25 1.5e-07
## JobRoleManager 0.580 0.667 0.87 0.38419
## JobRoleManufacturing Director -0.142 0.460 -0.31 0.75667
## JobRoleResearch Director 1.277 0.782 1.63 0.10259
## JobRoleResearch Scientist -1.293 0.364 -3.55 0.00038
## JobRoleSales Executive -1.268 0.981 -1.29 0.19616
## JobRoleSales Representative -2.888 0.988 -2.92 0.00345
##
## Scale fixed at 1
##
## Exponential distribution
## Loglik(model)= -996.6 Loglik(intercept only)= -1132
## Chisq= 270.78 on 13 degrees of freedom, p= 3e-50
## Number of Newton-Raphson Iterations: 17
## n= 1470ajuste_w <-survreg(survival~ OverTime + BusinessTravel + Department + JobRole, data = atrito, dist="weibull")
summary(ajuste_w)##
## Call:
## survreg(formula = survival ~ OverTime + BusinessTravel + Department +
## JobRole, data = atrito, dist = "weibull")
## Value Std. Error z p
## (Intercept) 18.457 1461.371 0.01 0.98992
## OverTimeYes -0.969 0.126 -7.67 1.7e-14
## BusinessTravelTravel_Frequently -0.794 0.280 -2.83 0.00460
## BusinessTravelTravel_Rarely -0.493 0.266 -1.85 0.06401
## DepartmentResearch & Development -12.980 1461.371 -0.01 0.99291
## DepartmentSales -12.818 1461.371 -0.01 0.99300
## JobRoleHuman Resources -14.542 1461.371 -0.01 0.99206
## JobRoleLaboratory Technician -1.713 0.322 -5.33 9.9e-08
## JobRoleManager 0.567 0.588 0.96 0.33524
## JobRoleManufacturing Director -0.141 0.406 -0.35 0.72889
## JobRoleResearch Director 1.163 0.691 1.68 0.09240
## JobRoleResearch Scientist -1.199 0.323 -3.71 0.00021
## JobRoleSales Executive -1.159 0.866 -1.34 0.18086
## JobRoleSales Representative -2.680 0.875 -3.06 0.00218
## Log(scale) -0.126 0.052 -2.42 0.01545
##
## Scale= 0.882
##
## Weibull distribution
## Loglik(model)= -993.8 Loglik(intercept only)= -1131.1
## Chisq= 274.46 on 13 degrees of freedom, p= 5.2e-51
## Number of Newton-Raphson Iterations: 17
## n= 1470ajuste_ln <-survreg(survival~ OverTime + BusinessTravel + Department + JobRole, data = atrito, dist="lognormal")
summary(ajuste_ln)##
## Call:
## survreg(formula = survival ~ OverTime + BusinessTravel + Department +
## JobRole, data = atrito, dist = "lognormal")
## Value Std. Error z p
## (Intercept) 11.4313 547.6868 0.02 0.98335
## OverTimeYes -1.0104 0.1290 -7.83 4.8e-15
## BusinessTravelTravel_Frequently -0.6602 0.2562 -2.58 0.00997
## BusinessTravelTravel_Rarely -0.3532 0.2348 -1.50 0.13257
## DepartmentResearch & Development -6.4138 547.6866 -0.01 0.99066
## DepartmentSales -6.0718 547.6867 -0.01 0.99115
## JobRoleHuman Resources -7.9262 547.6867 -0.01 0.98845
## JobRoleLaboratory Technician -1.6696 0.2821 -5.92 3.3e-09
## JobRoleManager 0.6899 0.4999 1.38 0.16758
## JobRoleManufacturing Director -0.0607 0.3481 -0.17 0.86164
## JobRoleResearch Director 1.3047 0.5772 2.26 0.02381
## JobRoleResearch Scientist -1.0480 0.2846 -3.68 0.00023
## JobRoleSales Executive -1.1503 0.7717 -1.49 0.13606
## JobRoleSales Representative -2.8245 0.7874 -3.59 0.00033
## Log(scale) 0.3589 0.0490 7.33 2.3e-13
##
## Scale= 1.43
##
## Log Normal distribution
## Loglik(model)= -985.6 Loglik(intercept only)= -1119.1
## Chisq= 266.95 on 13 degrees of freedom, p= 1.9e-49
## Number of Newton-Raphson Iterations: 13
## n= 1470Exponencial, Weibull e Lognormal
library(survival)
library(rms)
ajusteE <- psm( survival~OverTime + BusinessTravel + Department + JobRole , data = atrito, dist="exponential")
residE <- residuals(ajusteE)
ajusteW <- psm(survival~OverTime + BusinessTravel + Department + JobRole , data = atrito, dist="weibull")
residW <- residuals(ajusteW)
ajusteLN <- psm(survival~OverTime + BusinessTravel + Department + JobRole , data = atrito, dist="lognormal")
residLN <- residuals(ajusteLN)Escolha da melhor distribuição
Vemos que a Weibull se ajusta melhor aos dados
### Escolha da melhor distribuição
par(mfrow=c(1,3))
survplot(residE,main="Exponential",ylab="Complement of residual CDF")
survplot(residW,main="Weibull",ylab="Complement of residual CDF")
survplot(residLN,main="Lognormal",ylab="Complement of residual CDF")
residuos deviance
Nos residuos vemos que a distribuição dos pontos são semelhantes e ambos estão no intervalo de -3 a 3 então podemos considerar que não ha pontos de influência
### residuos deviance
residE_dev = residuals(ajusteE, type = 'deviance')
residW_dev = residuals(ajusteW, type = 'deviance')
residLN_dev = residuals(ajusteLN, type = 'deviance')par(mfrow=c(1,3))
plot( residE_dev, ylim = c(-3,3) )
plot( residW_dev, ylim = c(-3,3) )
plot( residLN_dev, ylim = c(-3,3) )
pontos influentes sobre os valores preditos
Vemos que a log normal possui mais pontos dispersos
### residuos ldresp: pontos influentes sobre os valores preditos
residE_ldresp = residuals(ajusteE, type = 'ldresp')
residW_ldresp = residuals(ajusteW, type = 'ldresp')
residLN_ldresp = residuals(ajusteLN, type = 'ldresp')
par(mfrow=c(1,3))
plot( residE_ldresp )
plot( residW_ldresp )
plot( residLN_ldresp )
# pontos influentes sobre os parâmetros
### residuos ldcase: pontos influentes sobre os parâmetros
residE_ldcase = residuals(ajusteE, type = 'ldcase')
residW_ldcase = residuals(ajusteW, type = 'ldcase')
residLN_ldcase = residuals(ajusteLN, type = 'ldcase')
par(mfrow=c(1,3))
plot( residE_ldcase )
plot( residW_ldcase )
plot( residLN_ldcase )
Modelo proporcional de cox
#### Modelo proporcional de cox
library(ggfortify)
library(survminer)
mod_cox = coxph( survival~OverTime + BusinessTravel + Department + JobRole , data = atrito )
summary(mod_cox)## Call:
## coxph(formula = survival ~ OverTime + BusinessTravel + Department +
## JobRole, data = atrito)
##
## n= 1470, number of events= 237
##
## coef exp(coef) se(coef) z
## OverTimeYes 1.152e+00 3.164e+00 1.318e-01 8.738
## BusinessTravelTravel_Frequently 9.558e-01 2.601e+00 3.157e-01 3.027
## BusinessTravelTravel_Rarely 5.659e-01 1.761e+00 3.017e-01 1.876
## DepartmentResearch & Development 1.449e+01 1.961e+06 1.274e+03 0.011
## DepartmentSales 1.422e+01 1.500e+06 1.274e+03 0.011
## JobRoleHuman Resources 1.634e+01 1.247e+07 1.274e+03 0.013
## JobRoleLaboratory Technician 2.042e+00 7.706e+00 3.818e-01 5.348
## JobRoleManager -5.881e-01 5.554e-01 6.837e-01 -0.860
## JobRoleManufacturing Director 3.049e-01 1.357e+00 4.768e-01 0.640
## JobRoleResearch Director -1.263e+00 2.827e-01 7.976e-01 -1.584
## JobRoleResearch Scientist 1.456e+00 4.290e+00 3.885e-01 3.748
## JobRoleSales Executive 1.547e+00 4.697e+00 9.935e-01 1.557
## JobRoleSales Representative 3.137e+00 2.304e+01 1.003e+00 3.128
## Pr(>|z|)
## OverTimeYes < 2e-16 ***
## BusinessTravelTravel_Frequently 0.002469 **
## BusinessTravelTravel_Rarely 0.060677 .
## DepartmentResearch & Development 0.990929
## DepartmentSales 0.991097
## JobRoleHuman Resources 0.989771
## JobRoleLaboratory Technician 8.87e-08 ***
## JobRoleManager 0.389742
## JobRoleManufacturing Director 0.522478
## JobRoleResearch Director 0.113247
## JobRoleResearch Scientist 0.000178 ***
## JobRoleSales Executive 0.119443
## JobRoleSales Representative 0.001758 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## OverTimeYes 3.164e+00 3.161e-01 2.44343 4.096
## BusinessTravelTravel_Frequently 2.601e+00 3.845e-01 1.40066 4.829
## BusinessTravelTravel_Rarely 1.761e+00 5.678e-01 0.97494 3.181
## DepartmentResearch & Development 1.961e+06 5.099e-07 0.00000 Inf
## DepartmentSales 1.500e+06 6.668e-07 0.00000 Inf
## JobRoleHuman Resources 1.247e+07 8.021e-08 0.00000 Inf
## JobRoleLaboratory Technician 7.706e+00 1.298e-01 3.64614 16.285
## JobRoleManager 5.554e-01 1.801e+00 0.14541 2.121
## JobRoleManufacturing Director 1.357e+00 7.372e-01 0.53280 3.454
## JobRoleResearch Director 2.827e-01 3.537e+00 0.05921 1.350
## JobRoleResearch Scientist 4.290e+00 2.331e-01 2.00341 9.188
## JobRoleSales Executive 4.697e+00 2.129e-01 0.67018 32.919
## JobRoleSales Representative 2.304e+01 4.341e-02 3.22725 164.436
##
## Concordance= 0.798 (se = 0.015 )
## Likelihood ratio test= 262 on 13 df, p=<2e-16
## Wald test = 223.4 on 13 df, p=<2e-16
## Score (logrank) test = 291.5 on 13 df, p=<2e-16test.ph = cox.zph(mod_cox)
test.ph## chisq df p
## OverTime 0.0888 1 0.766
## BusinessTravel 3.0550 2 0.217
## Department 2.7986 2 0.247
## JobRole 14.6302 8 0.067
## GLOBAL 21.1156 13 0.071testando proporcionalidade
### testando proporcionalidade
ggcoxzph(test.ph)
ggcoxdiagnostics( mod_cox, type = 'deviance', linear.predictions = F )
Conclusão
A partir da análise acima, podemos concluir que ‘Hora Extra’, ‘Departamento’, ‘Viagem a trabalho’ e ‘Cargo de trabalho’ afetam a taxa de atrito de funcionários. As curvas de sobrevivência mostram como os funcionários tem probabilidade de atrito em cada ponto do tempo. A mesma conclusão pode ser tirada do teste de Peto que resulta em um valor p significativo (<0,05). A distribuição que melhor de ajustou aos dados foi a Weibull.