Melakukan pengecekan missing value

melkukan pengecekan missing value pada data untuk setiap variabel.

colSums(is.na(attrition_clean))

#>               Attrition          BusinessTravel              Department 
#>                       0                       0                       0 
#> EnvironmentSatisfaction                OverTime          YearsAtCompany 
#>                       0                       0                       0

##Exploratory Data Analysis (EDA)

Summary Data

Pada bagian ini akan ditunjukkan beberapa statistik persebaran data untuk masing-masing variabel pada dataset attrition_clean. Dengan menggunakan fungsi summary() dihasilkan output sebagai berikut:

summary(attrition_clean)

#>  Attrition            BusinessTravel                  Department 
#>  No :1233   Non-Travel       : 150   Human Resources       : 63  
#>  Yes: 237   Travel_Frequently: 277   Research & Development:961  
#>             Travel_Rarely    :1043   Sales                 :446  
#>                                                                  
#>                                                                  
#>                                                                  
#>  EnvironmentSatisfaction OverTime   YearsAtCompany  
#>  1:284                   No :1054   Min.   : 0.000  
#>  2:287                   Yes: 416   1st Qu.: 3.000  
#>  3:453                              Median : 5.000  
#>  4:446                              Mean   : 7.008  
#>                                     3rd Qu.: 9.000  
#>                                     Max.   :40.000

Insight:

• Variabel kategorik berupa Attrition, BusinessTravel, Department, EnvironmentSatisfaction, dan Overtime merupakan tipe data kategorik. Sehingga summary hanya menunjukkan count dari masing-masing level kategori
  • Untuk variabel Attrition, jumlah pegawai yang mengalami event sebanyak 237 pegawai sedangkan yang tidak sebanyak 1233
  • Pada variabel BusinessTravel, tipe pegawai kebanyakan jarang melakukan perjalanan bisnis
  • Kemudian Department Research & Development sebagai department dengan jumlah pegawai terbanyak.
  • Nilai 3 merupakan nilai EnvironmentSatisfaction tertinggi yang dialami oleh pegawai
  • Sebanyak 1054 pegawai tidak bekerja lembur/Overtime.
• Variabel numerik pada data yaitu pada variabel YearsAtCompany memiliki rata-rata 7 tahun. Di mana pegawai memiliki lama bekerja di perusahaan selama 7 tahun, dengan rentang berkisar antara 0-40 tahun.

Distribusi Komponen Event pada data

Pada bagian ini akan ditunjukkan plot jumlah dan juga persentase pada variabel Attrition yang menjadi komponen event pada data kita. Hal ini bertujuan untuk melihat jumlah dan sebaran pegawai yang mengalami event berupa attrition maupun tidak.

Selanjutnya akan dibentuk tabel frekuensi dan persentase dari level kategori yang ada pada variabel Attrition.

attrition_jumlah <- attrition_clean %>% 
  group_by(Attrition) %>% 
  summarise(Count=n())

attrition_percentage <- attrition_jumlah %>%
  mutate(Percent=round(prop.table(Count),3) * 100)
attrition_percentage

#> # A tibble: 2 × 3
#>   Attrition Count Percent
#>   <fct>     <int>   <dbl>
#> 1 No         1233    83.9
#> 2 Yes         237    16.1

Berikut merupakan plot untuk menampilkan distribusi komponen Event pada data.

# Menunjukkan plot untuk Employee Attrition dalam jumlah
attrition_jumlah <- attrition_jumlah %>% ggplot(aes(x=Attrition, y=Count, fill = Attrition)) + 
  geom_bar(stat="identity", color="grey40") + 
  theme_economist() + 
  scale_fill_manual(values=c("#9FF781", "#FA5858")) +
  coord_flip() + 
  geom_text(aes(x=Attrition, y=0.01, label= Count),
            hjust=-0.8, vjust=-1, size=4, 
            colour="black", fontface="bold",
         angle=360) + labs(title="Employee Attrition (Jumlah)", x="Employee Attrition",y="Jumlah") + 
  theme(plot.title=element_text(hjust=0.5))

# Menunjukkan plot untuk Employee Attrition dalam persen
attrition_percentage <- attrition_percentage %>% 
  ggplot(aes(x = Attrition, y = Percent, fill = Attrition)) + 
  geom_bar(stat="identity", color="grey40") +
  geom_text(aes(x=Attrition, y=0.01, label= sprintf("%.2f%%", Percent)),
            hjust=0.5, vjust=-3, size=4, 
            colour="black", fontface="bold") + 
  theme_economist() +
  scale_fill_manual(values=c("#9FF781", "#FA5858")) +
  labs(x="Employee Attrition", y="Percentage") + 
  labs(title="Employee Attrition (%)") + 
  theme(plot.title=element_text(hjust=0.5))

plot_grid(attrition_jumlah, attrition_percentage, align="h", ncol=2)

Insight : Jumlah pegawai yang mengalami event jauh lebih sedikit dari yang tidak mengalami event. Di mana untuk pegawai yang mengalami event berjumlah 237 pegawai atau sebanyak 16.10% dari jumlah total pegawai.

Distribusi variabel Department terhadap Attrition level

Pada tahapan EDA berikut, akan diperlihatkan distribusi variabel departement untuk maisng-masing level kategori yang ada pada komponen event atau variabel Attrition. Berikut merupakan hasil codenya:

employee_department <- attrition_clean %>% 
  select(Department, Attrition) %>% 
  group_by(Department, Attrition) %>%
  summarize(Count = n()) %>%
ggplot(aes(x=reorder(Department, Count), y=Count, fill=Attrition)) + 
  geom_bar(stat="identity", position="dodge") + facet_wrap(~Attrition) + 
  theme_economist() + 
  theme(axis.text.x = element_text(angle = 90, size = 9), plot.title=element_text(hjust=0.5)) + 
  scale_fill_manual(values=c("#9FF781", "#FA5858")) + 
labs(y="Jumlah Pegawai", x="Department", title="Jumlah pegawai berdasarkan Attrition dan Department") + 
  geom_text(aes(x=Department, y=0.01, label= paste0(round(Count,2))),
            hjust=-0.5, vjust=0, size=3.5, 
            colour="black", fontface = "bold", angle=90)


employee_department

Insight

• Dari plot visualisasi di atas diketahui bahwa Departemen Research & Development jumlah pegawai yang mengalami attrition sebanyak 133 atau sebanyak 133/961 (14%) dari total karyawan di departemen tersebut.

• Selain itu, departemen Sales sebanyak 92 pegawai mengalami attrition atau sebanyak 92/446 (20,6%) dari total karyawan di departemen tersebut.

• Jumlah pegawai yang mengalami attrition pada Departemen HR sebanyak 12 atau 12/63 (19%)

EnvironmentSatisfaction Overview berdasarkan Department dan Attrition

Pada tahapan EDA selanjutnya akan ditunjukan dua buah visualisasi bersesuaian.

1. Pertama mengenai visualisai rata-rata kepuasan pegawai berdasarkan departemen dan lingkungan kerja atau variabel EnvironmentSatisfaction

Selanjutnya akan dibentuk tabel aggregasi untuk rata-rata nilai EnvironmentSatisfaction untuk setiap level Departemen dan level status Attrition.

env <- attrition_new %>% select(EnvironmentSatisfaction, Department, Attrition) %>% group_by(Department, Attrition) %>%
summarize(Rata_rata=mean(EnvironmentSatisfaction))
env

#> # A tibble: 6 × 3
#> # Groups:   Department [3]
#>   Department             Attrition Rata_rata
#>   <chr>                  <chr>         <dbl>
#> 1 Human Resources        No             2.76
#> 2 Human Resources        Yes            2.33
#> 3 Research & Development No             2.79
#> 4 Research & Development Yes            2.47
#> 5 Sales                  No             2.73
#> 6 Sales                  Yes            2.47

Berikut merupakan plot untuk menampilkan rata-rata kepuasan pegawai berdasarkan departemen dan rata-rata kepuasan terhadap lingkungan kerja atau variabel EnvironmentSatisfaction

ggplot(env, aes(x=Department, y=Rata_rata)) +
  geom_line(aes(group=Attrition), color="#58ACFA", linetype="dashed") + 
  geom_point(aes(color=Attrition), size=3) +  
  theme_economist() + 
  theme(plot.title=element_text(hjust=0.5), axis.text.x=element_text(angle=0)) + 
labs(title="Environment Satisfaction berdasarkan Department", y="Rata-rata Environment Satisfaction", x="Departemen") + scale_color_manual(values=c("#58FA58", "#FA5858"))

Insight :

• Terlihat bahwa untuk masing-masing level Attrition, pegawai yang mengalami attrition dalam hal ini keluar dari perusahaan memiliki nilai EnvironmentSatisfaction yang rendah dibandingkan dengan pegawai yang tidak keluar dari perusahaan

• Departemen HR memiliki rata-rata nilai terendah dalam hal kepuasan lingkungan kerja pada pegawai yang memiliki level attrition = yes dari perusahaan.

• Departemen Sales memiliki rata-rata nilai terendah dalam hal kepuasan lingkungan kerja pada pegawai yang memiliki level attrition = no dari perusahaan.

2. Kedua mengenai visualisai nilai masing-masing kepuasan pegawai berdasarkan departemen terhadap lingkungan kerja atau variabel EnvironmentSatisfaction

attrition_yes <- attrition_clean %>% 
  filter(Attrition == "Yes")

department_env <- attrition_yes %>% 
  select(Department, EnvironmentSatisfaction) %>% 
  group_by(Department, EnvironmentSatisfaction) %>%
  summarize(count=n()) %>% 
  ggplot(aes(x=EnvironmentSatisfaction, y=count, fill=Department))+
  geom_bar(stat='identity') + 
  facet_wrap(~Department) + 
  theme_economist() + 
  theme(legend.position="bottom", plot.title=element_text(hjust=0.5)) + 
  scale_fill_manual(values=c("#84b08b", "#848bb0", "#bd95b6")) + 
  geom_label(aes(label=count, fill = Department), colour = "white", fontface = "italic") + 
  labs(title="Environment Satisfaction by Attrition Status and Department", x="Enivironment Satisfaction", y="Jumlah pegawai")

department_env

Insight :

• Pada keseluruhan departemen baik Human Resources, Research & Developmen, serta Sales, tingkat EnvironmentSatisfaction 1 terbanyak.

• Departemen Human Resorces dan Sales memiliki distribusi tingkat EnvironmentSatisfaction yang cenderung rata pada setiap levelnya.

Tingkat Attrition berdasarkan Overtime status

overtime_attrition <- attrition_yes %>%
  select(OverTime, Attrition) %>% 
  group_by(Attrition, OverTime) %>% 
  summarize(Jumlah = n()) %>%
  mutate(Persentase = round(prop.table(Jumlah),2) * 100)
overtime_attrition

#> # A tibble: 2 × 4
#> # Groups:   Attrition [1]
#>   Attrition OverTime Jumlah Persentase
#>   <fct>     <fct>     <int>      <dbl>
#> 1 Yes       No          110         46
#> 2 Yes       Yes         127         54

overtime_attrition %>% 
  ggplot(aes(x = OverTime, y = Jumlah, fill = OverTime)) + 
  geom_bar(stat="identity") + 
  scale_fill_manual(values=c("#75ab67", "#cca06a")) + 
  geom_label(aes(label=paste0(Jumlah)), fill="#FFF9F5", colour = "black", fontface = "italic")+
  labs(title="Tingkat Attrition berdasarkan Overtime status", 
       subtitle="Dalam Jumlah", 
       x="Overtime Status", 
       y="Number of Employees") + 
  theme_economist() + 
  theme(legend.position="bottom", 
        plot.title=element_text(hjust=0.5))

Melihat Data Waktu Survival dengan Indikator Teramati dan Tersensor

attrition_clean$event <- with(attrition_clean, ifelse(Attrition=="Yes", 1, 0)) #mendefinisikan kolom attrition pada data sebagai kejadian yang diamati
survival <- Surv(attrition_clean$YearsAtCompany, attrition_clean$event)
head(survival, 52)

#>  [1]  6  10+  0   8+  2+  7+  1+  1+  9+  7+  5+  9+  5+  2+  4  10+  6+  1+ 25+
#> [20]  3+  4+  5  12+  0+  4  14+ 10   9+ 22+  2+  1+  4+ 10+  1   2   5+  3   2+
#> [39]  1+  5+  1+  1+  1   9+ 12+ 22   9+  1+  9+  1+  1   2

Insight:

• Dari output di atas diketahui bahwa waktu survival dari pegawai pertama adalah 6 dan event teramati.
  • Artinya pegawai pertama keluar dari perusahaan setelah 6 tahun bekerja di perusahaan tersebut.
• Sementara pegawai ke-14 memiliki waktu survival 2 dan event tidak teramati atau tesensor artinya pegawai ke-14 keluar dari perusahaan namun karena sebab lain seperti meninggal atau dipindah tugaskan ke cabang perusahaan lain. Dengan demikian, pegawai ke-14 dikatakan tersensor di T=C=2.

Fungsi Survival tanpa Mempertimbangkan Faktor Lain

Pada bagian ini kita akan melakukan fit untuk fungsi survival dengan menggunakan fungsi survit(). Berikut kita akan membuat sebuah fungsi survival sederhana yang tidak mempertimbangkan perbedaan faktor lain. Sehingga kita hanya akan menggunakan interceptnya ~1.

model_none <- survfit(Surv(YearsAtCompany, event)~1, data = attrition_clean)
model_none

#> Call: survfit(formula = Surv(YearsAtCompany, event) ~ 1, data = attrition_clean)
#> 
#>         n events median 0.95LCL 0.95UCL
#> [1,] 1470    237     40      32      NA

Dari hasil analisis pada output dapat dilihat bahwa banyaknya observasi ada 1470 dan hanya 237 pegawai yang mengalami event. Akan dilihat ringkasan yang lebih lengkap dari model, sebagai berikut:

summary(model_none)

#> Call: survfit(formula = Surv(YearsAtCompany, event) ~ 1, data = attrition_clean)
#> 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0   1470      16    0.989 0.00271        0.984        0.994
#>     1   1426      59    0.948 0.00583        0.937        0.960
#>     2   1255      27    0.928 0.00690        0.914        0.941
#>     3   1128      20    0.911 0.00769        0.896        0.927
#>     4   1000      19    0.894 0.00851        0.877        0.911
#>     5    890      21    0.873 0.00947        0.855        0.892
#>     6    694       9    0.862 0.01007        0.842        0.882
#>     7    618      11    0.846 0.01091        0.825        0.868
#>     8    528       9    0.832 0.01173        0.809        0.855
#>     9    448       8    0.817 0.01264        0.793        0.842
#>    10    366      18    0.777 0.01516        0.748        0.807
#>    11    246       2    0.770 0.01568        0.740        0.802
#>    13    200       2    0.763 0.01644        0.731        0.796
#>    14    176       2    0.754 0.01736        0.721        0.789
#>    15    158       1    0.749 0.01789        0.715        0.785
#>    16    138       1    0.744 0.01857        0.708        0.781
#>    17    126       1    0.738 0.01934        0.701        0.777
#>    18    117       1    0.732 0.02018        0.693        0.772
#>    19    104       1    0.725 0.02117        0.684        0.767
#>    20     93       1    0.717 0.02233        0.674        0.762
#>    21     66       1    0.706 0.02449        0.660        0.756
#>    22     52       1    0.692 0.02753        0.641        0.749
#>    23     37       1    0.674 0.03253        0.613        0.741
#>    24     35       1    0.654 0.03686        0.586        0.731
#>    31     16       1    0.614 0.05256        0.519        0.726
#>    32     13       1    0.566 0.06641        0.450        0.713
#>    33     10       1    0.510 0.08037        0.374        0.694
#>    40      1       1    0.000     NaN           NA           NA

Insight dari output di atas:

• Pada awal pengamatan yaitu waktu 0 (time=0), terdapat 1470 pegawai (n.risk=1470) yang beresiko mengalami event. Enam belas pegawai di antaranya mengalami event pada saat itu (n.event=16). Sehingga diperoleh 98.9% pegawai (survival=0.989) yang berhasil survive, yaitu melewati waktu 0 tahun tanpa mengalami event.

• Pada waktu 16 (time=16), terdapat 138 pegawai (n.risk=138) yang beresiko mengalami event. Satu pegawai di antaranya mengalami event pada saat itu (n.event=1). Sehingga diperoleh 74.4% pegawai (survival=0.744) yang berhasil survive, yaitu melewati waktu 16 tahun tanpa mengalami event.

• Sebanyak 67.4% pegawai masih belum keluar dari perusahaan 23 tahun setelah studi dilakukan atau dapat dinyatakan bahwa peluang untuk survive (tidak keluar dari perusahaan) hingga tahun ke-23 adalah sekitar 0.67.

• Interpretasi untuk titik waktu yang lain serupa.

• Di akhir masa studi, yaitu saat time=40, terdapat satu pegawai yang beresiko untuk mengalami event, dan satu orang tersebut akan mengalami event. Sehingga tidak menyisakan pegawai yang masih survive

• Dapat dikatakan bahwa tidak ada pegawai yang tersensor dikarenakan masa studi telah selesai. pegawai dapat tersensor dikarenakan sebab lain yang membuat pegawai keluar dari pengamatan pada masa studi sedang berlangsung.

Kaplan Meier Curve untuk Fungsi Survival tanpa Mempertimbangkan Faktor Lain

Fungsi survival di atas dapat disajikan secara grafis, sebagai berikut:

1. Menggunakan fungsi plot() pada base untuk melakukan visualisasi terhadap fungsi survival yang sudah dibuat sebelumnya. Karena merupakan plot dari base, hasil visualiasi masih sangat sederhana.

plot(model_none, xlab="Waktu survival", ylab="S(t)",main="Survival Plot")

Insight:

• Garis utuh menyatakan taksiran fungsi survival dan garis putus-putus menyatakan batas atas dan batas bawah (interval kepercayaan 95%) untuk fungsi survival pada data employee attrition tanpa mempertimbangkan faktor lain.

2. Cara lain untuk bisa melakukan visualisasi terhadap fungsi survival dapat menggunakan fungsi ggsurvplot() yang berasal dari library survminer. Hal ini memungkinkan kita melakukan customisasi plot yang kita inginkan.

library(survminer)
ggsurvplot(model_none, conf.int=TRUE, pval=TRUE, risk.table=TRUE,
           palette= "#6aa329", 
           title="Kaplan-Meier Curve untuk Employee Attrition tanpa melibatkan Faktor Lain",
           break.time.by = 4,
           ggtheme = theme_minimal(),
           font.tickslab = c(12, "bold", "#941241"),
           risk.table.fontsize = 3,
           tables.col = "darkred")

Fungsi Survival berdasarkan Overtime Status

Sebelumnya sudah diestimasi fungsi survival tanpa memperhitungkan faktor lain. Namun mungkin saja faktor-faktor lain dapat menjelaskan cepat atau lambatnya seorang pegawai untuk keluar dari perusahaan atau mengalami event. Misalkan ingin diketahui fungsi survival dengan mempertimbangkan efek Overtime.

model_overtime <- survfit(Surv(YearsAtCompany, event) ~ OverTime, data = attrition_clean)
model_overtime

#> Call: survfit(formula = Surv(YearsAtCompany, event) ~ OverTime, data = attrition_clean)
#> 
#>                 n events median 0.95LCL 0.95UCL
#> OverTime=No  1054    110     40      32      NA
#> OverTime=Yes  416    127     24      16      NA

Insight Dari hasil analisis pada output dapat ditarik beberapa poin penting di antaranya:

• Dilihat bahwa banyaknya pegawai yang bekerja tanpa Overtime terdapat sebanyak 1054 pegawai

• Banyaknya pegawai yang bekerja Overtime ada sebanyak 416 pegawai.

• Dapat dilihat juga bahwa pada kelompok pegawai yang Overtime terdapat 127 pegawai yang keluar dari perusahaan atau mengalami event. Jumlah tersebut lebih banyak dari pada banyaknya event yang terjadi pada kelompok pegawai yang tidak bekerja Overtime.

Sehingga dapat diduga bahwa pegawai yang bekerja Overtime dapat memberikan efek pada cepatnya seorang pegawai keluar dari perusahaan. Berikut adalah ringkasan lebih lengkap dari fungsi survival berdasarkan Overtime:

summary(model_overtime)

#> Call: survfit(formula = Surv(YearsAtCompany, event) ~ OverTime, data = attrition_clean)
#> 
#>                 OverTime=No 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0   1054       6    0.994 0.00232        0.990        0.999
#>     1   1024      31    0.964 0.00578        0.953        0.976
#>     2    908      12    0.951 0.00677        0.938        0.965
#>     3    816       6    0.944 0.00730        0.930        0.959
#>     4    728      11    0.930 0.00836        0.914        0.947
#>     5    648       7    0.920 0.00909        0.902        0.938
#>     6    509       4    0.913 0.00971        0.894        0.932
#>     7    457       3    0.907 0.01025        0.887        0.927
#>     8    388       3    0.900 0.01094        0.879        0.922
#>     9    332       6    0.884 0.01260        0.859        0.909
#>    10    269       8    0.857 0.01527        0.828        0.888
#>    11    181       1    0.853 0.01590        0.822        0.884
#>    13    146       2    0.841 0.01770        0.807        0.876
#>    14    126       2    0.828 0.01978        0.790        0.867
#>    17     89       1    0.818 0.02163        0.777        0.862
#>    18     82       1    0.808 0.02356        0.763        0.856
#>    20     65       1    0.796 0.02627        0.746        0.849
#>    22     35       1    0.773 0.03397        0.709        0.843
#>    23     24       1    0.741 0.04532        0.657        0.835
#>    32      8       1    0.648 0.09528        0.486        0.865
#>    33      6       1    0.540 0.12663        0.341        0.855
#>    40      1       1    0.000     NaN           NA           NA
#> 
#>                 OverTime=Yes 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0    416      10    0.976 0.00751        0.961        0.991
#>     1    402      28    0.908 0.01423        0.881        0.936
#>     2    347      15    0.869 0.01684        0.836        0.902
#>     3    312      14    0.830 0.01903        0.793        0.868
#>     4    272       8    0.805 0.02034        0.766        0.846
#>     5    242      14    0.759 0.02265        0.716        0.804
#>     6    185       5    0.738 0.02383        0.693        0.786
#>     7    161       8    0.702 0.02593        0.653        0.754
#>     8    140       6    0.672 0.02757        0.620        0.728
#>     9    116       2    0.660 0.02829        0.607        0.718
#>    10     97      10    0.592 0.03254        0.531        0.659
#>    11     65       1    0.583 0.03329        0.521        0.652
#>    15     50       1    0.571 0.03460        0.507        0.643
#>    16     43       1    0.558 0.03626        0.491        0.634
#>    19     32       1    0.540 0.03909        0.469        0.623
#>    21     22       1    0.516 0.04437        0.436        0.611
#>    24     13       1    0.476 0.05595        0.378        0.599
#>    31      7       1    0.408 0.07916        0.279        0.597

Insight:

• Pada awal pengamatan kelompok pegawai yang bekerja dengan Overtime = "No", yaitu waktu 0 (time=0), terdapat 1054 pegawai (n.risk=1054) yang beresiko mengalami event. Enam pegawai di antaranya mengalami event pada saat itu (n.event=6). Sehingga diperoleh 99.4% pegawai (survival=0.994) yang berhasil survive, yaitu melewati waktu 0 tahun tanpa mengalami event.

• Pada awal pengamatan kelompok pegawai yang bekerja dengan Overtime = "Yes", yaitu waktu 0 (time=0), terdapat 416 pegawai (n.risk=416) yang beresiko mengalami event. Sepuluh pegawai di antaranya mengalami event pada saat itu (n.event=10). Sehingga diperoleh 97.6% pegawai (survival=0.976) yang berhasil survive, yaitu melewati waktu 0 tahun tanpa mengalami event.

• Interpretasi untuk titik waktu yang lain serupa.

• Di akhir masa studi kelompok pegawai yang bekerja dengan Overtime="No", yaitu saat (time=40), terdapat 1 pegawai yang beresiko untuk mengalami event, dan 1 orang tersebut akan mengalami event. Sehingga tidak menyisakan pegawai yang masih survive.

• Di akhir masa studi kelompok pegawai yang bekerja dengan Overtime="Yes", yaitu saat time=31, terdapat 7 pegawai yang beresiko untuk mengalami event, dan 1 orang pegawai yang akan mengalami event.

• Dapat dilihat bahwa kelompok pegawai yang bekerja Overtime hanya dapat survive di perusahaan sampai tahun ke-31 sedangkan kelompok pegawai yang bekerja tidak Overtime dapat survive di perusahaan sampai tahun ke-40.

Kaplan Meier Curve untuk Fungsi Survival berdasarkan Overtime Status

Untuk memudahkan membandingkan fungsi survival dari kedua kelompok pegawai berdasarkan level Overtime tersebut, disajikan secara grafis visualisasi berikut:

ggsurvplot(model_overtime, conf.int=TRUE, pval=TRUE, risk.table=TRUE,
           palette= c("#6aa329", "#bda237"), 
           title="Kaplan-Meier Curve untuk Employee Attrition berdasarkan Overtime Status",
           break.time.by = 4,
           ggtheme = theme_minimal(),
           font.tickslab = c(12, "bold", "#941241"),
           risk.table.fontsize = 3,
           tables.theme = theme_bw(),
           tables.col = "darkred")

- Dapat dilihat dari plot yang terbentuk bahwa kelompok pegawai yang bekerja tidak Overtime (Overtime = "No") dapat bertahan lebih lama tidak terkena event ditandai dengan grafik fungsi survival yang lebih tinggi dari pada kelompok pegawai yang bekerja Overtime (Overtime = "Yes").

• Hal ini berarti peluang pegawai yang bekerja tidak Overtime untuk melewati suatu titik tanpa terkena event lebih tinggi dibandingkan dengan pegawai yang bekerja Overtime.

Fungsi Survival berdasarkan Environment_Satisfaction Status

Sebelumnya sudah diestimasi fungsi survival tanpa memperhitungkan faktor lain dan dengan memperhitungkan faktor Overtime. Namun mungkin saja ada faktor-faktor lain yang dapat menjelaskan cepat atau lambatnya seorang pegawai untuk keluar dari perusahaan atau mengalami event. Misalkan ingin diketahui fungsi survival dengan mempertimbangkan efek Environment Satisfaction.

model_envs <- survfit(Surv(YearsAtCompany, event) ~ EnvironmentSatisfaction, data = attrition_clean)
model_envs

#> Call: survfit(formula = Surv(YearsAtCompany, event) ~ EnvironmentSatisfaction, 
#>     data = attrition_clean)
#> 
#>                             n events median 0.95LCL 0.95UCL
#> EnvironmentSatisfaction=1 284     72     32      19      NA
#> EnvironmentSatisfaction=2 287     43     NA      22      NA
#> EnvironmentSatisfaction=3 453     62     NA      33      NA
#> EnvironmentSatisfaction=4 446     60     40      31      NA

Insight - Dari hasil analisis pada output dapat dilihat bahwa banyaknya pegawai yang dengan Low Environment Satisfaction terdapat sebanyak 284, Medium Environment Satisfaction sebanyak 287, High Environment Satisfaction sebanyak 453, dan Very High Environment Satisfaction sebanyak 446 pegawai.

• Dapat dilihat juga bahwa pada kelompok pegawai yang menghasilkan Low Environment Satisfaction terdapat 72 pegawai yang mengalami event dari total 284 pegawai yang berada di kelompok tersebut.

• Dapat dilihat juga bahwa pada kelompok pegawai yang menghasilkan Very High Environment Satisfaction terdapat 60 pegawai yang mengalami event dari total 446 pegawai yang berada di kelompok tersebut.

• Sehingga dapat diduga bahwa pegawai yang menghasilkan Low Environment Satisfaction dapat memberikan efek paling cepat pada seorang pegawai untuk keluar dari perusahaan. Dan pegawai yang menghasilkan Very High Environment Satisfaction dapat memberikan efek paling lambat pada seorang pegawai keluar dari perusahaan.

Berikut adalah ringkasan lebih lengkap dari fungsi survival berdasarkan EnvironmentSatisfaction:

summary(model_envs)

#> Call: survfit(formula = Surv(YearsAtCompany, event) ~ EnvironmentSatisfaction, 
#>     data = attrition_clean)
#> 
#>                 EnvironmentSatisfaction=1 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0    284       3    0.989 0.00607        0.978        1.000
#>     1    278      19    0.922 0.01601        0.891        0.954
#>     2    242       9    0.888 0.01906        0.851        0.926
#>     3    212       9    0.850 0.02200        0.808        0.894
#>     4    178       7    0.816 0.02450        0.770        0.866
#>     5    154       6    0.785 0.02676        0.734        0.839
#>     6    115       1    0.778 0.02739        0.726        0.833
#>     7    107       2    0.763 0.02874        0.709        0.822
#>     8     98       1    0.755 0.02948        0.700        0.816
#>     9     86       3    0.729 0.03214        0.669        0.795
#>    10     71       5    0.678 0.03719        0.609        0.755
#>    13     41       1    0.661 0.03978        0.588        0.744
#>    14     37       1    0.643 0.04253        0.565        0.732
#>    15     34       1    0.624 0.04530        0.542        0.720
#>    18     26       1    0.600 0.04951        0.511        0.706
#>    19     23       1    0.574 0.05381        0.478        0.690
#>    24     10       1    0.517 0.07289        0.392        0.681
#>    32      3       1    0.345 0.14884        0.148        0.803
#> 
#>                 EnvironmentSatisfaction=2 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0    287       3    0.990  0.0060        0.978        1.000
#>     1    275       9    0.957  0.0121        0.934        0.981
#>     2    243       6    0.934  0.0152        0.904        0.964
#>     3    221       5    0.912  0.0175        0.879        0.947
#>     4    201       3    0.899  0.0189        0.862        0.937
#>     5    180       2    0.889  0.0200        0.850        0.929
#>     6    149       2    0.877  0.0214        0.836        0.920
#>     7    128       6    0.836  0.0262        0.786        0.889
#>     8    108       1    0.828  0.0271        0.777        0.883
#>    10     78       4    0.786  0.0330        0.724        0.853
#>    20     18       1    0.742  0.0526        0.646        0.853
#>    22      8       1    0.649  0.0982        0.483        0.873
#> 
#>                 EnvironmentSatisfaction=3 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0    453       6    0.987 0.00537        0.976        0.997
#>     1    442      12    0.960 0.00925        0.942        0.978
#>     2    390       8    0.940 0.01138        0.918        0.963
#>     3    352       3    0.932 0.01219        0.909        0.956
#>     4    318       6    0.915 0.01391        0.888        0.942
#>     5    290       6    0.896 0.01562        0.866        0.927
#>     6    224       6    0.872 0.01802        0.837        0.908
#>     8    170       4    0.851 0.02030        0.812        0.892
#>     9    137       3    0.833 0.02253        0.790        0.878
#>    10    114       3    0.811 0.02524        0.763        0.862
#>    16     44       1    0.792 0.03066        0.734        0.855
#>    17     38       1    0.771 0.03626        0.704        0.846
#>    21     22       1    0.736 0.04870        0.647        0.838
#>    23     14       1    0.684 0.06792        0.563        0.831
#>    33      5       1    0.547 0.13384        0.339        0.884
#> 
#>                 EnvironmentSatisfaction=4 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0    446       4    0.991 0.00446        0.982        1.000
#>     1    431      19    0.947 0.01069        0.927        0.969
#>     2    380       4    0.937 0.01168        0.915        0.961
#>     3    343       3    0.929 0.01250        0.905        0.954
#>     4    303       3    0.920 0.01346        0.894        0.947
#>     5    266       7    0.896 0.01591        0.865        0.928
#>     7    188       3    0.881 0.01767        0.848        0.917
#>     8    152       3    0.864 0.01997        0.826        0.904
#>     9    132       2    0.851 0.02171        0.809        0.895
#>    10    103       6    0.801 0.02835        0.748        0.859
#>    11     71       2    0.779 0.03173        0.719        0.844
#>    13     57       1    0.765 0.03399        0.701        0.835
#>    14     50       1    0.750 0.03659        0.681        0.825
#>    31      6       1    0.625 0.11809        0.431        0.905
#>    40      1       1    0.000     NaN           NA           NA

Kaplan Meier Curve untuk Fungsi Survival berdasarkan Environment Satisfaction

Untuk memudahkan membandingkan fungsi survival dari keempat kelompok pegawai tersebut, disajikan secara grafis visualisasi berikut:

ggsurvplot(model_envs, pval=TRUE, risk.table=TRUE,
           palette= c("#6aa329", "#ed9715", "#269091", "#91267a"), 
           title="Kaplan-Meier Curve untuk Employee Attrition
           berdasarkan Overtime Status",
           break.time.by = 4,
           ggtheme = theme_minimal(),
           font.tickslab = c(12, "bold", "#941241"),
           risk.table.fontsize = 3,
           risk.table.height = 0.4,
           tables.col = "darkred",
           tables.theme = theme_bw(),
           legend = "none")

Insight:

• Dapat dilihat dari plot yang terbentuk bahwa kelompok pegawai yang menghasilkan Low Environment Satisfaction tidak dapat bertahan lebih lama untuk tidak terkena event ditandai dengan grafik fungsi survival yang paling rendah di antara kelompok pegawai yang lain.

• Hal ini berarti peluang pegawai yang menghasilkan Medium Environment Satisfaction, High Environment Satisfaction, Very High Environment Satisfaction untuk melewati suatu titik tanpa terkena event lebih tinggi dibandingkan dengan pegawai yang bekerja menghasilkan Law Environment Satisfaction.

Fungsi Survival berdasarkan Department

Sebelumnya sudah diestimasi fungsi survival tanpa memperhitungkan faktor lain, dengan memperhitungkan faktor Overtime, dan dengan memperhitungkan faktor Environment Satisfaction. Namun mungkin saja ada faktor-faktor lain yang dapat menjelaskan cepat atau lambatnya seorang pegawai untuk keluar dari perusahaan atau mengalami event. Misalkan ingin diketahui fungsi survival dengan mempertimbangkan efek Department. Berikut merupakan code yang dapat dijalankan:

model_dept <- survfit(Surv(YearsAtCompany, event) ~ Department, data = attrition_clean)
model_dept

#> Call: survfit(formula = Surv(YearsAtCompany, event) ~ Department, data = attrition_clean)
#> 
#>                                     n events median 0.95LCL 0.95UCL
#> Department=Human Resources         63     12     NA      20      NA
#> Department=Research & Development 961    133     40      33      NA
#> Department=Sales                  446     92     32      23      NA

Insight:

• Dari hasil analisis pada output dapat dilihat bahwa banyaknya pegawai yang berasal dari HR Department ada sebanyak 63, Research & Development Department ada sebanyak 961, dan Sales Department ada sebanyak 446.

• Dapat dilihat juga bahwa terdapat 92 pegawai yang mengalami event dari total 446 pegawai yang berada di kelompok Sales Department.

• Dapat dilihat juga bahwa terdapat 133 pegawai yang mengalami event dari total 961 pegawai yang berada di kelompok Research & Development Department.

• Dapat diduga bahwa pegawai yang berasal dari Research & Development Department dapat memberikan efek paling lambat pada seorang pegawai keluar dari perusahaan. Dan pegawai yang berasal dari Sales Department dapat memberikan efek paling cepat pada seorang pegawai keluar dari perusahaan.

Berikut adalah ringkasan lebih lengkap dari fungsi survival berdasarkan Department:

summary(model_dept)

#> Call: survfit(formula = Surv(YearsAtCompany, event) ~ Department, data = attrition_clean)
#> 
#>                 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
#>     0    961       9    0.991 0.00311        0.985        0.997
#>     1    934      38    0.950 0.00706        0.937        0.964
#>     2    812      13    0.935 0.00811        0.919        0.951
#>     3    735       8    0.925 0.00879        0.908        0.942
#>     4    650      11    0.909 0.00982        0.890        0.929
#>     5    571      13    0.889 0.01115        0.867        0.911
#>     6    441       4    0.881 0.01176        0.858        0.904
#>     7    389       8    0.862 0.01314        0.837        0.889
#>     8    332       4    0.852 0.01397        0.825        0.880
#>     9    280       4    0.840 0.01504        0.811        0.870
#>    10    230      14    0.789 0.01936        0.752        0.828
#>    13    127       1    0.783 0.02018        0.744        0.823
#>    17     81       1    0.773 0.02212        0.731        0.817
#>    18     74       1    0.762 0.02416        0.716        0.811
#>    22     32       1    0.739 0.03313        0.676        0.806
#>    31      9       1    0.657 0.08279        0.513        0.841
#>    33      5       1    0.525 0.13483        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
#>     0    446       7    0.984 0.00589        0.973        0.996
#>     1    429      17    0.945 0.01086        0.924        0.967
#>     2    385      12    0.916 0.01344        0.890        0.943
#>     3    343      10    0.889 0.01548        0.859        0.920
#>     4    308       7    0.869 0.01691        0.836        0.903
#>     5    281       7    0.847 0.01836        0.812        0.884
#>     6    226       5    0.829 0.01977        0.791        0.868
#>     7    205       2    0.820 0.02039        0.781        0.861
#>     8    176       5    0.797 0.02232        0.755        0.842
#>     9    150       4    0.776 0.02412        0.730        0.825
#>    10    121       4    0.750 0.02651        0.700        0.804
#>    11     82       2    0.732 0.02885        0.678        0.791
#>    13     66       1    0.721 0.03047        0.664        0.783
#>    14     57       2    0.696 0.03425        0.632        0.766
#>    15     51       1    0.682 0.03619        0.615        0.757
#>    16     41       1    0.665 0.03895        0.593        0.746
#>    19     31       1    0.644 0.04320        0.564        0.734
#>    21     21       1    0.613 0.05087        0.521        0.721
#>    23     13       1    0.566 0.06526        0.452        0.710
#>    24     12       1    0.519 0.07495        0.391        0.689
#>    32      5       1    0.415 0.11050        0.246        0.699

Untuk memudahkan membandingkan fungsi survival dari ketiga kelompok pegawai tersebut, disajikan secara grafis visualisasi berikut:

ggsurvplot(model_dept, pval=TRUE, risk.table="percentage",
           palette= c("#e0ba31", "#2cc774", "#c72c84"), 
           title="Kaplan-Meier Curve untuk Employee Attrition
           berdasarkan Departemen",
           break.time.by = 5,
           ggtheme = theme_minimal(),
           font.tickslab = c(12, "bold", "#941241"),
           risk.table.fontsize = 3,
           risk.table.height = 0.3,
           tables.col = "black",
           tables.theme = theme_bw(),
           legend = "none")

Insight: - Dapat dilihat dari plot yang terbentuk bahwa kelompok pegawai yang berasal dari Sales Department tidak dapat bertahan lebih lama untuk tidak terkena event ditandai dengan grafik fungsi survival yang paling rendah di antara kelompok pegawai yang lain.

• Dapat dilihat dari plot yang terbentuk bahwa kelompok pegawai yang berasal dari Research & Development Department dapat bertahan lebih lama untuk tidak terkena event ditandai dengan grafik fungsi survival yang paling tinggi di antara kelompok pegawai yang lain.

Fungsi Survival berdasarkan Business Travel

Sebelumnya sudah diestimasi fungsi survival tanpa memperhitungkan faktor lain, dengan memperhitungkan faktor Overtime, dengan memperhitungkan faktor Environment Satisfaction, dan dengan memperhitungkan faktor department. Namun mungkin saja ada faktor lain yang dapat menjelaskan cepat atau lambatnya seorang pegawai untuk keluar dari perusahaan atau mengalami event. Misalkan ingin diketahui fungsi survival dengan mempertimbangkan efek Business Travel. Berikut code untuk membentuk fungsi survivalnya:

model_travel<- survfit(Surv(YearsAtCompany, event) ~ BusinessTravel, data = attrition_clean)
model_travel

#> Call: survfit(formula = Surv(YearsAtCompany, event) ~ BusinessTravel, 
#>     data = attrition_clean)
#> 
#>                                     n events median 0.95LCL 0.95UCL
#> BusinessTravel=Non-Travel         150     12     NA      NA      NA
#> BusinessTravel=Travel_Frequently  277     69     NA      21      NA
#> BusinessTravel=Travel_Rarely     1043    156     33      31      NA

Insight:

• Dari hasil analisis pada output dapat dilihat bahwa banyaknya pegawai yang berasal dari kelompok Non-Travel ada sebanyak 150, Travel Frequently ada sebanyak 277, dan Travel Rarely ada sebanyak 1043.

Dapat dilihat juga bahwa terdapat 12 pegawai yang mengalami event dari total 150 pegawai yang berada di kelompok Non-Travel.

• Dapat dilihat juga bahwa terdapat 69 pegawai yang mengalami event dari total 277 pegawai yang berada di kelompok Travel Frequently.

• Sehingga dapat diduga bahwa pegawai yang berasal dari kelompok Travel Frequently dapat memberikan efek paling cepat pada seorang pegawai keluar dari perusahaan. Dan pegawai yang berasal dari kelompok Non-Travel dapat memberikan efek paling lambat pada seorang pegawai keluar dari perusahaan.

Berikut adalah ringkasan lebih lengkap dari fungsi survival berdasarkan Business Travel:

summary(model_travel)

#> Call: survfit(formula = Surv(YearsAtCompany, event) ~ BusinessTravel, 
#>     data = attrition_clean)
#> 
#>                 BusinessTravel=Non-Travel 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0    150       1    0.993 0.00664        0.980        1.000
#>     1    143       4    0.966 0.01514        0.936        0.996
#>     2    125       1    0.958 0.01688        0.925        0.991
#>     3    119       1    0.950 0.01856        0.914        0.987
#>     4    103       1    0.941 0.02054        0.901        0.982
#>     5     94       2    0.921 0.02450        0.874        0.970
#>    10     38       2    0.872 0.04063        0.796        0.955
#> 
#>                 BusinessTravel=Travel_Frequently 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0    277       7    0.975 0.00943        0.956        0.993
#>     1    266      17    0.912 0.01708        0.880        0.947
#>     2    241       7    0.886 0.01930        0.849        0.925
#>     3    211       6    0.861 0.02131        0.820        0.904
#>     4    191       6    0.834 0.02333        0.789        0.881
#>     5    168       4    0.814 0.02479        0.767        0.864
#>     6    142       4    0.791 0.02661        0.740        0.845
#>     7    123       3    0.772 0.02820        0.718        0.829
#>     8    107       3    0.750 0.03005        0.693        0.811
#>     9     90       2    0.733 0.03161        0.674        0.798
#>    10     71       5    0.682 0.03687        0.613        0.758
#>    11     46       1    0.667 0.03893        0.595        0.748
#>    13     39       1    0.650 0.04152        0.573        0.736
#>    18     26       1    0.625 0.04684        0.539        0.724
#>    21     15       1    0.583 0.05942        0.478        0.712
#>    23      8       1    0.510 0.08574        0.367        0.709
#> 
#>                 BusinessTravel=Travel_Rarely 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0   1043       8    0.992 0.00270        0.987        0.998
#>     1   1017      38    0.955 0.00645        0.943        0.968
#>     2    889      19    0.935 0.00783        0.920        0.950
#>     3    798      13    0.920 0.00877        0.903        0.937
#>     4    706      12    0.904 0.00971        0.885        0.923
#>     5    628      15    0.882 0.01096        0.861        0.904
#>     6    479       5    0.873 0.01160        0.851        0.896
#>     7    429       8    0.857 0.01273        0.832        0.882
#>     8    361       6    0.843 0.01378        0.816        0.870
#>     9    306       6    0.826 0.01507        0.797        0.856
#>    10    257      11    0.791 0.01780        0.757        0.826
#>    11    171       1    0.786 0.01829        0.751        0.823
#>    13    139       1    0.780 0.01901        0.744        0.819
#>    14    126       2    0.768 0.02063        0.729        0.810
#>    15    112       1    0.761 0.02156        0.720        0.805
#>    16    100       1    0.754 0.02264        0.711        0.799
#>    17     89       1    0.745 0.02392        0.700        0.794
#>    19     73       1    0.735 0.02568        0.686        0.787
#>    20     66       1    0.724 0.02760        0.672        0.780
#>    22     37       1    0.704 0.03307        0.642        0.772
#>    24     25       1    0.676 0.04207        0.598        0.764
#>    31     11       1    0.615 0.06998        0.492        0.768
#>    32      9       1    0.546 0.08952        0.396        0.753
#>    33      8       1    0.478 0.10108        0.316        0.724
#>    40      1       1    0.000     NaN           NA           NA

Untuk memudahkan membandingkan fungsi survival dari ketiga kelompok pegawai tersebut, disajikan secara grafis visualisasi berikut:

ggsurvplot(model_travel, pval=TRUE, risk.table="percentage",
           palette= c("#e04c31", "#31e09d", "#5731e0"), 
           title="Kaplan-Meier Curve untuk Employee Attrition
           berdasarkan Business Travel",
           break.time.by = 5,
           ggtheme = theme_minimal(),
           font.tickslab = c(12, "bold", "#941241"),
           risk.table.fontsize = 3,
           risk.table.height = 0.3,
           tables.col = "black",
           tables.theme = theme_bw(),
           legend = "none")

Insight:

• Dapat dilihat dari plot yang terbentuk bahwa kelompok pegawai yang berasal dari Travel Frequently tidak dapat bertahan lebih lama untuk tidak terkena event ditandai dengan grafik fungsi survival yang paling rendah di antara kelompok pegawai yang lain.

• Dapat dilihat dari plot yang terbentuk bahwa kelompok pegawai yang berasal dari Non-Travel dapat bertahan lebih lama untuk tidak terkena event ditandai dengan grafik fungsi survival yang paling tinggi di antara kelompok pegawai yang lain.

Visualisasi fungsi kumulatif hazard berdasarkan Overtime

plot(model_overtime, 
     lty = 1:2,
     col = 1:2,
     lwd = 2,
     main = expression(paste("Estimasi Kaplan Meier ", hat (H)(t), " berdasarkan Variabel Overtime")),
     xlab = "Waktu Survival",
     ylab = "Cumulative Hazard",
     fun = "cumhaz",
     legend = T)
legend('topleft', c("No", "Yes"), lty = 1:2, col = 1:2, lwd = 2)

Visualisasi fungsi kumulatif hazard berdasarkan EnvironmentSatisfaction

plot(model_envs,
     lty = 1:4,
     col = 1:4,
     lwd = 2,
     main = expression(paste("Estimasi Kaplan-Meier ",hat (H)(t), " berdasarkan Variabel EnvironmentSatisfaction")),
     xlab = "Waktu survival",
     ylab = "Cumulative Hazard",
     fun = "cumhaz")
legend('topleft', 
       c("Envsatisfaction1", "Envsatisfaction2","Envsatisfaction3","Envsatisfaction4"),
       lty = 1:4,
       col=1:4,
       lwd = 2,
       bty="n")

Visualisasi fungsi kumulatif hazard berdasarkan Department

plot(model_dept,
     lty = 1:3,
     col = 1:3,
     lwd = 2,
     main = expression(paste("Estimasi Kaplan-Meier ",hat (H)(t), " berdasarkan Variabel Department")),
     xlab = "Waktu survival",
     ylab = "Cumulative Hazard",
     fun = "cumhaz")
legend('topleft', 
       c("Human Resources ","Research & Development ","Sales"),
       lty = 1:3,
       col=1:3,
       lwd = 2,
       bty="n")

Visualisasi fungsi kumulatif hazard berdasarkan BusinessTravel

plot(model_travel,
     lty = 1:3,
     col = 1:3,
     lwd = 2,
     main = expression(paste("Estimasi Kaplan-Meier ",hat (H)(t), " berdasarkan Variabel BusinessTravel")),
     xlab = "Waktu survival",
     ylab = "Cumulative Hazard",
     fun = "cumhaz")
legend('topleft', 
       c("Non-Travel ","Travel_Frequently ","Travel_Rarely "),
       lty = 1:3,
       col = 1:3,
       lwd = 2,
       bty="n")

Log-Rank test untuk faktor Overtime

Tujuan : untuk mengetahui apakah pegawai yang mengalami kelebihan jam kerja dari jam kerja rata-rata dengan pegawai yang bekerja sesuai dengan jam kerja rata-rata memiliki survival experience yang berbeda atau tidak. Dalam hal ini, akan dibandingkan 2 kelompok pegawai yang berbeda.

Berikut adalah R code untuk pengujian ini :

h_overtime <- survdiff(Surv(YearsAtCompany, event)~ OverTime, data = attrition_clean)
h_overtime

#> Call:
#> survdiff(formula = Surv(YearsAtCompany, event) ~ OverTime, data = attrition_clean)
#> 
#>                 N Observed Expected (O-E)^2/E (O-E)^2/V
#> OverTime=No  1054      110    171.4      22.0      81.7
#> OverTime=Yes  416      127     65.6      57.5      81.7
#> 
#>  Chisq= 81.7  on 1 degrees of freedom, p= <0.0000000000000002

Insight: - Dengan signifikansi \(\alpha\) = 0.05 didapatkan nilai p-value < 0.05 sehingga \(H_0\) ditolak.

• Dapat disimpulkan bahwa terdapat perbedaan survival experience di antara grup pegawai yang mengalami kelebihan jam kerja dari jam kerja rata-rata dengan pegawai yang bekerja sesuai dengan jam kerja rata-rata.

Log-Rank test untuk faktor EnvironmentSatisfaction

Tujuan : untuk mengetahui apakah terdapat perbedaan survival experience antara kelompok pegawai dengan efek environment satisfaction (low, medium, high, very high).

Berikut adalah R code untuk pengujian ini

h_satisfaction <- survdiff(Surv(YearsAtCompany, event)~ EnvironmentSatisfaction, data = attrition_clean)
h_satisfaction

#> Call:
#> survdiff(formula = Surv(YearsAtCompany, event) ~ EnvironmentSatisfaction, 
#>     data = attrition_clean)
#> 
#>                             N Observed Expected (O-E)^2/E (O-E)^2/V
#> EnvironmentSatisfaction=1 284       72     44.8    16.556    20.971
#> EnvironmentSatisfaction=2 287       43     46.1     0.207     0.265
#> EnvironmentSatisfaction=3 453       62     74.6     2.135     3.204
#> EnvironmentSatisfaction=4 446       60     71.5     1.853     2.748
#> 
#>  Chisq= 21.3  on 3 degrees of freedom, p= 0.00009

Insight: - Dengan signifikansi 0.05 didapatkan nilai p-value <0.05 sehingga \(H_0\) ditolak.

• Dapat disimpulkan bahwa terdapat perbedaan survival experience di antara minimal 1 grup pegawai dengan efek environment satisfaction yang berbeda.

• Walaupun hasil pengujian ini tidak menyatakan tingkat environment satisfaction mana yang berbeda, tapi secara implisit dari perbandingan Observed dan Expected, dapat dihasilkan beberapa hal, yaitu :

  • Dapat diduga bahwa pada Tingkat 1 (Observed lebih tinggi daripada Expected, sehingga efek dari rendahnya tingkat kepuasan sekitar terhadap keinginan untuk keluar dari perusahaan pada pegawai lebih buruk, karena jumlah employee attrition lebih tinggi dari pada yang seharusnya).

  • Sementara pada Tingkat 4 (Observed lebih rendah dibandingkan Expected, yang mana mengimplikasikan, efek dari tingginya tingkat kepuasan sekitar terhadap keinginan untuk keluar dari perusahaan pada pegawai lebih baik, karena jumlah employee attrition lebih rendah dari pada yang seharusnya).

  • Sehingga dapat disimpulkan bahwa paling tidak Environment Satisfaction tingkat 1 dan tingkat 4 memiliki survival experience yang berbeda.

Log-Rank test untuk faktor Department

Tujuan : untuk mengetahui apakah terdapat perbedaan survival experience antara kelompok pegawai yang berasal dari Departemen (HR, Research and Development, Sales).

Berikut adalah R code untuk pengujian ini :

h_department <- survdiff(Surv(YearsAtCompany, event)~ Department, data = attrition_clean)
h_department

#> Call:
#> survdiff(formula = Surv(YearsAtCompany, event) ~ Department, 
#>     data = attrition_clean)
#> 
#>                                     N Observed Expected (O-E)^2/E (O-E)^2/V
#> Department=Human Resources         63       12     10.5     0.215     0.231
#> Department=Research & Development 961      133    152.6     2.521     7.283
#> Department=Sales                  446       92     73.9     4.440     6.629
#> 
#>  Chisq= 7.4  on 2 degrees of freedom, p= 0.02

Insight: - Dengan signifikansi 0.05 didapatkan nilai p-value <0.05 sehingga \(H_0\) ditolak.

• Dapat disimpulkan bahwa terdapat perbedaan survival experience di antara minimal 1 Departemen yang berbeda.

• Walaupun hasil pengujian ini tidak menyatakan department mana yang mempunyai survival experience yang berbeda, tapi secara implisit dari perbandingan Observed dan Expected, dapat dihasilkan beberapa hal, yaitu :

  • Dapat diduga bahwa pada Departemen Sales (Observed lebih tinggi daripada Expected, sehingga efek dari pegawai yang berasal dari Departemen Sales terhadap keinginan untuk keluar dari perusahaan pada pegawai lebih buruk, karena jumlah employee attrition lebih tinggi dari pada yang seharusnya).

  • Sementara pada Departemen Resource and Development (Observed lebih rendah dibandingkan Expected, yang mana mengimplikasikan, efek dari pegawai yang berasal dari Departemen Resource and Development terhadap keinginan untuk keluar dari perusahaan lebih baik, karena jumlah employee attrition lebih rendah dari pada yang seharusnya).

  • Sehingga dapat disimpulkan bahwa paling tidak pegawai yang berasal dari Departemen Sales dan Resource and Development memiliki survival experience yang berbeda.

Log-Rank test untuk faktor BusinessTravel

Tujuan : untuk mengetahui apakah terdapat perbedaan survival experience antara kelompok pegawai dengan frekuensi melakukan perjalanan bisnis dengan 3 kategori (non_travel, Travel_Frequently, Travel_Rarely). Berikut adalah R code untuk pengujian ini :

h_travel <- survdiff(Surv(YearsAtCompany, event)~ BusinessTravel, data = attrition_clean)
h_travel

#> Call:
#> survdiff(formula = Surv(YearsAtCompany, event) ~ BusinessTravel, 
#>     data = attrition_clean)
#> 
#>                                     N Observed Expected (O-E)^2/E (O-E)^2/V
#> BusinessTravel=Non-Travel         150       12     24.5     6.413      7.34
#> BusinessTravel=Travel_Frequently  277       69     45.0    12.820     16.25
#> BusinessTravel=Travel_Rarely     1043      156    167.5     0.785      2.75
#> 
#>  Chisq= 20.6  on 2 degrees of freedom, p= 0.00003

Insight: - Dengan signifikansi 0.05 didapatkan nilai p-value <0.05 sehingga \(H_0\) ditolak.

• Dapat disimpulkan bahwa terdapat perbedaan survival experience di antara minimal 1 kategori frekuensi perjalanan bisnis yang berbeda.

• Walaupun hasil pengujian ini tidak menyatakan kategori frekuensi perjalanan mana yang mempunyai survival experience yang berbeda, tapi secara implisit dari perbandingan Observed dan Expected, dapat dihasilkan beberapa hal, yaitu :

  • Dapat diduga bahwa pada pegawai yang sering menjalankan perjalanan bisnis (Observed lebih tinggi daripada Expected, sehingga efek dari pegawai yang sering menjalankan perjalanan bisnis terhadap keinginan untuk keluar dari perusahaan lebih buruk, karena jumlah employee attrition lebih tinggi dari pada yang seharusnya).

  • Sementara pada pegawai yang jarang melakukan perjalanan bisnis (Observed lebih rendah dibandingkan Expected, yang mana mengimplikasikan, efek dari pegawai yang jarang menjalani perjalanan bisnis terhadap keinginan untuk keluar dari perusahaan lebih baik, karena jumlah employee attrition lebih rendah daripada yang seharusnya).

  • Sehingga dapat disimpulkan bahwa paling tidak pegawai yang sering dan jarang melakukan perjalanan bisnis dari memiliki survival experience yang berbeda.

Konsep Model Cox-PH

CoxPH menggambarkan hubungan antara distribusi survival function dan kovariat. Variabel predictor diekspresikan oleh fungsi hazard sebagai berikut:

\[ h(t, x) = h_0(t)e^{\beta_1x_1+...+\beta_px_p}\] Metode CoxPH merupakan model yang semi-parametrik karena terdiri dari 2 komponen:

• Komponen non-parametrik \[h_0(t)\]

yang disebut sebagai baseline hazard, yaitu hazard bila semua kovariat bernilai nol.

• Komponen parametrik

\[e^{\beta_1x_1+...+\beta_px_p}\]

disebut sebagai partial hazard yang tidak bergantung dengan waktu.

• Secara umum, model Cox membuat estimasi fungsi log-risk

\[h(t, x)\]

sebagai kombinasi linier dari kovariat statis dan baseline hazard.

Konsep Model Regresi Cox-PH

Model regresi Cox Proportional Hazard merupakan metode semiparametrik yang dapat digunakan untuk memodelkan hubungan antara variabel respon yang berupa waktu survival dengan satu atau lebih variabel prediktor. Dalam penerapannya, analisis regresi cox-PH memiliki 2 pendekatan dalam melakukan estimasi parameter, yaitu:

1. Jika baseline hazard h0(t) diasumsikan berasal dari distribusi tertentu, maka ini disebut metode parametrik. Atau,

2. Distribusi dari h0(t) diabaikan. Inilah yang dilakukan pada model Cox-PH.

Dalam estimasi parameter, perlu diperhatikan dua hal berikut:

1. Apakah event teramati atau tersensor; dan

2. Apakah ada data ties atau tidak.

Event teramati berkontribusi secara langsung dalam konstruksi likelihood, sementara data tersensor berkontribusi secara tidak langsung dalam pembentukan himpunan resiko di tiap titik waktu dimana terjadi event. Adapun yang dimaksud dengan data ties adalah jika ada beberapa event terjadi bersamaan di satu titik waktu. Untuk data seperti ini, perlu ditentukan urutan kejadian yang mungkin dalam proses konstruksi likelihood.

Metode Breslow digunakan apabila terdapat data ties.

Konsep metode breslow dalam konstruksi likelihood model Cox-PH \(h(t,x)\) dimana terjadi D events pada rentang waktu pengamatan adalah sebagai berikut

\[ L_1(\beta) = \prod_{i=1}^{D} \frac{exp(\beta'Z_i)}{[\sum_{l \in R_i} exp(\beta'x_l)]^{d_i}}\]

Aplikasi pada data Attrition

Dalam pengaplikasian pada data attrition akan digunakan metode Breslow dalam melakukan analisis Regresi Cox PH. Hal ini disebabkan pada waktu \(t_i\) dimungkinkan terjadi \(d_i\) event. Kejadian-kejadian tersebut diasumsikan terjadi satu per-satu. Dengan metode coxph() menggunakan model Surv() dan memasukkan semua variabel-variabel dependen dan independen dihasilkan output sebagai berikut:

cox <- coxph(Surv(YearsAtCompany, event) ~ OverTime + EnvironmentSatisfaction + Department + BusinessTravel, data = attrition_clean, method="breslow")
summary(cox)

#> Call:
#> coxph(formula = Surv(YearsAtCompany, event) ~ OverTime + EnvironmentSatisfaction + 
#>     Department + BusinessTravel, data = attrition_clean, method = "breslow")
#> 
#>   n= 1470, number of events= 237 
#> 
#>                                     coef exp(coef) se(coef)      z
#> OverTimeYes                       1.1389    3.1233   0.1315  8.660
#> EnvironmentSatisfaction2         -0.5952    0.5514   0.1939 -3.071
#> EnvironmentSatisfaction3         -0.7577    0.4687   0.1747 -4.337
#> EnvironmentSatisfaction4         -0.7604    0.4675   0.1768 -4.300
#> DepartmentResearch & Development -0.2349    0.7907   0.3041 -0.773
#> DepartmentSales                   0.1316    1.1406   0.3095  0.425
#> BusinessTravelTravel_Frequently   1.1375    3.1188   0.3151  3.610
#> BusinessTravelTravel_Rarely       0.7251    2.0649   0.3019  2.402
#>                                              Pr(>|z|)    
#> OverTimeYes                      < 0.0000000000000002 ***
#> EnvironmentSatisfaction2                     0.002136 ** 
#> EnvironmentSatisfaction3                    0.0000145 ***
#> EnvironmentSatisfaction4                    0.0000171 ***
#> DepartmentResearch & Development             0.439803    
#> DepartmentSales                              0.670708    
#> BusinessTravelTravel_Frequently              0.000306 ***
#> BusinessTravelTravel_Rarely                  0.016311 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#>                                  exp(coef) exp(-coef) lower .95 upper .95
#> OverTimeYes                         3.1233     0.3202    2.4137    4.0417
#> EnvironmentSatisfaction2            0.5514     1.8135    0.3771    0.8063
#> EnvironmentSatisfaction3            0.4687     2.1335    0.3328    0.6602
#> EnvironmentSatisfaction4            0.4675     2.1391    0.3305    0.6611
#> DepartmentResearch & Development    0.7907     1.2648    0.4357    1.4348
#> DepartmentSales                     1.1406     0.8767    0.6219    2.0920
#> BusinessTravelTravel_Frequently     3.1188     0.3206    1.6818    5.7837
#> BusinessTravelTravel_Rarely         2.0649     0.4843    1.1427    3.7314
#> 
#> Concordance= 0.711  (se = 0.018 )
#> Likelihood ratio test= 118.9  on 8 df,   p=<0.0000000000000002
#> Wald test            = 121.5  on 8 df,   p=<0.0000000000000002
#> Score (logrank) test = 130.6  on 8 df,   p=<0.0000000000000002

Insight:

• Terlihat dari output yang dihasilkan, hazard pegawai dengan (overtime = Yes) untuk mengalami employee attrition adalah 3.1233 kali dibandingkan dengan hazard untuk pegawai tanpa overtime. Hal ini mengimplikasikan bahwa pegawai tanpa overtime memiliki risiko yang lebih rendah untuk mengalami employee attrition.

• Di antara 4 kondisi kepuasan terhadap lingkungan (EnvironmentSatisfaction), pegawai yang memberi nilai 1 untuk EnvironmentSatisfaction memiliki risiko yang paling besar untuk mengalami employee attrition. Bila dibandingkan dengan pegawai yang memberi nilai 1, hazard pada pegawai yang memberi nilai 2 adalah 0.5514 kalinya, hazard pada pegawai yang memberi nilai 3 adalah 0.4687 kalinya, dan hazard pada pegawai yang memberi nilai 4 adalah 0.4675 kalinya.

• Di antara 3 jenis departemen, pegawai dari department Sales memiliki risiko paling besar untuk mengalami employee attrition. Bila dibandingkan dengan departemen Human Resources, hazard pegawai dari departemen Research & Development adalah 0.7907 kalinya, dan hazard pegawai dari departemen Sales adalah 1.1406 kalinya.

• Di antara 3 kondisi frekuensi menjalani BusinessTravel, pegawai yang tidak menjalani business travel memiliki risiko paling kecil untuk mengalami employee attrition. Bila dibandingkan dengan hazard pegawai tanpa business travel, hazard pegawai dengan kondisi Travel Frequently adalah 3.1188 kalinya dan hazard pegawai dengan kondisi Travel Rarely adalah 2.0649 kalinya.

• Uji Simultan Analisis Regresi Cox PH :

  • Hipotesis :
    • \(H_0\) = variabel independen tidak berpengaruh terhadap variabel dependen
    • \(H_1\) = Minimal terdapat satu variabel independen yang berpengaruh terhadap variabel dependen
  • Kesimpulan:
    • Model ini cukup baik dalam menjelaskan kondisi pegawai.
    • Terlihat dari uji model p-value untuk uji rasio likelihood adalah 2e-16 yang lebih kecil dari taraf signifikansi/alpha 0.05. Hal ini mengimplikasikan bahwa minimal salah satu dari jenis variabel, dapat menjelaskan lama waktu hingga pegawai mengalami employee attrition.

• Uji Parsial(Satu-satu variabel) Analisis Regresi Cox PH :

  • Berdasarkan uji parsial koefisien regresi, terlihat bahwa masing-masing pengukuran tersebut cukup signifikan dalam menjelaskan cepat-lambatnya employee attrition, kecuali untuk departemen, yaitu departemen R&D dan Sales tidak berbeda signifikan hazardnya dengan departemen HR.

Pembentukan Model Regresi

Berdasarkan output yang dihasilkan pada sub-bab sebelumnya, dihasilkan model regresi Cox-PH sebagai berikut:

\[h(t,x) = h_0(t)exp^{(\beta_1x_1+\beta_2x_{2,2}+\beta_3x_{2,3}+\beta_4x_{2,4}+\beta_5x_{3,2}+\beta_6x_{3,3}+\beta_7x_{4,2}+\beta_8x_{4,3})}\] \[h(t,x) = h_0(t)exp^{(1.138x_1-0.595x_{2,2}-0.757x_{2,3}-0.76_4x_{2,4}-0.235x_{3,2}+0.131x_{3,3}+1.137x_{4,2}+0.725x_{4,3})}\]

di mana:

• \(x_1\) = Variabel Overtime (1 = Yes, 0 = No)
• \(x_{2,j}\) = Variabel EnvironmentSatisfaction (2 = Medium, 3 = High, 4 = Very High, 0 = lainnya)
• \(x_{3,j}\) = Variabel Department (2 = Research and Development, 3 = Sales, 0 = lainnya)
• \(x_{4,j}\) = Variabel TravelBusiness (2 = Travel Frequently, 3 = Travel Rarely, 0 = lainnya)