2. Xử lí số liệu và thực hiện các phân tích khám phá sơ bộ EDA
# Load một số gói cho xử lí và phân tích sơ bộ:
library(magrittr)
library(tidyverse)
# Đọc dữ liệu và đánh giá sơ bộ về dữ liệu:
rm(list = ls())
ChurnStudy <- read.csv("F:/R_project/w1/churn_data.csv")
ChurnStudy %>% str()
## 'data.frame': 2000 obs. of 10 variables:
## $ Churn : int 1 0 0 1 1 1 1 0 0 0 ...
## $ Xeducation : Factor w/ 3 levels "Bachelor's Degree",..: 3 1 1 1 1 1 1 1 1 1 ...
## $ Xgender : Factor w/ 2 levels "F","M": 1 1 2 2 2 2 2 2 2 1 ...
## $ Xsatisfaction : int 5 2 5 3 3 5 1 1 5 5 ...
## $ Xsatisfaction2 : int 2 4 5 2 5 4 3 5 5 3 ...
## $ Xservice.calls : int 0 0 0 0 3 1 0 0 2 0 ...
## $ Xtenure2 : int 9 10 11 1 3 7 5 5 12 9 ...
## $ Xpurch.last.month: int 1 3 0 1 1 1 1 1 0 3 ...
## $ Xmonthly.charges : num 260.9 -18.2 243.8 45.6 269 ...
## $ time : int 9 10 11 1 3 7 5 5 12 9 ...
ChurnStudy %>% dim()
## [1] 2000 10
Dữ liệu thiếu là một vấn đề thường xuyên gặp trong thực tế. Do vậy chúng ta cần đánh giá mức độ thiếu của dữ liệu:
sapply(ChurnStudy, function(x) {x %>% is.na() %>% sum()})
## Churn Xeducation Xgender Xsatisfaction
## 0 0 0 7
## Xsatisfaction2 Xservice.calls Xtenure2 Xpurch.last.month
## 9 0 0 0
## Xmonthly.charges time
## 0 0
Do tỉ lệ NA là thấp nên chúng ta có thể chấp nhận phương án là bỏ các quan sát có NA:
ChurnStudy %<>% na.omit()
Xem qua số liệu:
ChurnStudy %>% summary()
## Churn Xeducation Xgender Xsatisfaction
## Min. :0.0 Bachelor's Degree:1593 F: 452 Min. :1.000
## 1st Qu.:0.0 Doctorate Degree : 122 M:1532 1st Qu.:2.000
## Median :0.5 Master's Degree : 269 Median :3.000
## Mean :0.5 Mean :2.942
## 3rd Qu.:1.0 3rd Qu.:4.000
## Max. :1.0 Max. :5.000
## Xsatisfaction2 Xservice.calls Xtenure2 Xpurch.last.month
## Min. :1.000 Min. :0.0000 Min. : 1.000 Min. : 0.000
## 1st Qu.:2.000 1st Qu.:0.0000 1st Qu.: 7.000 1st Qu.: 1.000
## Median :3.000 Median :0.0000 Median : 9.000 Median : 1.000
## Mean :3.512 Mean :0.4718 Mean : 8.449 Mean : 1.479
## 3rd Qu.:5.000 3rd Qu.:1.0000 3rd Qu.:11.000 3rd Qu.: 1.000
## Max. :6.000 Max. :5.0000 Max. :12.000 Max. :10.000
## Xmonthly.charges time
## Min. :-97.53 Min. : 1.000
## 1st Qu.: 67.80 1st Qu.: 7.000
## Median :131.48 Median : 9.000
## Mean :143.74 Mean : 8.449
## 3rd Qu.:221.42 3rd Qu.:11.000
## Max. :381.82 Max. :12.000
Biến Xmonthly.charges là số tiền thanh toán ở tháng trước nhưng lại có số âm. Đây là những tình huống mà khách hàng mua theo hình thức trả chậm (mua chịu). Số những cases này là 77:
# Số người mua theo hình thức trả chậm:
sum(ChurnStudy$Xmonthly.charges < 0)
## [1] 77
# Chiếm tỉ lệ 3.88%:
77*100 / nrow(ChurnStudy)
## [1] 3.881048
Để đơn giản cho phân tích chúng ta chuyển các thương vụ mua theo hình thức trả chậm này (số tiền thanh toán có giá trị âm) thành dương:
convert_to_pos <- function(x) {case_when(x < 0 ~ -x,
x > 0 ~ x)}
# Sử dụng hàm vừa có:
ChurnStudy %<>% mutate(Xmonthly.charges = convert_to_pos(Xmonthly.charges))
Học vấn của khách hàng:
hoc_van <- ChurnStudy$Xeducation %>% unique()
hoc_van
## [1] Master's Degree Bachelor's Degree Doctorate Degree
## Levels: Bachelor's Degree Doctorate Degree Master's Degree
# Viết hàm dán lại nhãn học vấn cho gọn:
convert_hoc_van <- function(x) {case_when(x == hoc_van[1] ~ "Master",
x == hoc_van[2] ~ "Bachelor",
x == hoc_van[3] ~ "Doctorate")}
# Dùng hàm này dán lại nhãn cho biến học vấn:
ChurnStudy %<>% mutate(Xeducation = convert_hoc_van(Xeducation))
Các phân tích khám phá (EDA - Exploratory Data Analysis):
library(hrbrthemes)
theme_set(theme_ipsum()) # Cố định theme.
# Cơ cấu học vấn theo giới tính không có nhiều khác biệt:
ChurnStudy %>%
group_by(Xgender, Xeducation) %>%
count() %>%
ggplot(aes(Xgender, n, fill = Xeducation)) +
geom_col(position = "fill") +
coord_flip() +
scale_y_percent() +
scale_fill_ipsum(name = "Education Level") +
theme(legend.position = "top")
# Rời bỏ cũng không khác biệt nhiều giữa các nhóm học vấn khác nhau:
ChurnStudy %>%
mutate(Churn = as.factor(Churn)) %>%
group_by(Churn, Xeducation) %>%
count() %>%
ggplot(aes(Churn, n, fill = Xeducation)) +
geom_col(position = "fill") +
coord_flip() +
scale_y_percent() +
scale_fill_ipsum(name = "Education Level") +
theme(legend.position = "top")
# Cũng không có khác biệt dáng kể về giới tính và học vấn:
ChurnStudy %>%
mutate(Churn = as.factor(Churn)) %>%
group_by(Churn, Xeducation, Xgender) %>%
count() %>%
ggplot(aes(Churn, n, fill = Xeducation)) +
geom_col(position = "fill") +
coord_flip() +
scale_y_percent() +
scale_fill_ipsum(name = "Education Level") +
theme(legend.position = "top") +
facet_wrap(~ Xgender)
Số tiền mua hàng trung bình theo từng nhóm giới tính và học vấn. Nhìn chung nam giới ở tất cả các nhóm học vấn dành nhiều tiền hơn để mua hàng:
ChurnStudy %>%
group_by(Xgender, Xeducation) %>%
summarise_each(funs(mean), Xmonthly.charges) %>%
ggplot(aes(Xgender, Xmonthly.charges, fill = Xeducation)) +
geom_col(position = "dodge") +
scale_fill_ipsum(name = "Education Level")
Tương tự là số lần mua hàng trung bình mỗi tháng thì nam giới cũng cao hơn:
ChurnStudy %>%
group_by(Xgender, Xeducation) %>%
summarise_each(funs(mean), Xpurch.last.month) %>%
ggplot(aes(Xgender, Xpurch.last.month, fill = Xeducation)) +
geom_col(position = "dodge") +
scale_fill_ipsum(name = "Education Level")
Số tiền và số lần mua hàng trung bình mỗi tháng có khác biệt bất kể trình độ học vấn cho ta nhận định ban đầu rằng rất có thể tỉ lệ rời bỏ / ở lại giữa hai nhóm khách hàng là khác biệt. Nhận định này cần được củng cố vững chắc hơn bằng cách bằng chứng thống kê từ phân tích SA dưới đây.
3. Các mô hình thống kê cho phân tích sống còn SA
Mục này trình bày hai cách tiếp cận thường sử dụng cho phân tích SA là: (1) mô hình Kaplan - Meir, và (2) hồi quy Cox (Cox Regression).
Cách tiếp cận của Kaplan - Meir là một phương pháp phi tham số áp dụng trong trường hợp mô hình chỉ có các biến định tính còn hồi quy Cox áp dụng cho trường hợp mà mô hình ngoài biến định tính còn có biến định lượng.
3.1 Mô hình KM (Kaplan - Meir Model)
Mục dưới đây chúng ta sẽ áp dụng một số mô hình thống kê thường được áp dụng cho phân tích SA. Cách tiếp cận của Kaplan và Meir thường được sử dụng trong tình huống chúng ta muốn trả lời câu hỏi: liệu tỉ lệ rời bỏ / ở lại có khác biệt giữa các nhóm khách hàng khi căn cứ vào giới tính, trình độ học vấn của họ hay không?.
Trước hết thực hiện phân chia dữ liệu:
# Dữ liệu huấn luyện:
set.seed(1020)
train <- ChurnStudy %>% sample_n(round(0.75*nrow(.)))
# Dữ liệu kiểm định:
test <- dplyr::setdiff(ChurnStudy, train)
3.1.1 Đánh giá chung về tỉ lệ rời bỏ / ở lại
Dưới đây là sử dụng cách tiếp cận của KM để nghiên cứu tỉ lệ ỏ lại / rời bỏ chung của khách hàng:
#-------------------------------------------------------------------
# Thực hiện SA cho toàn bộ không phân biệt giới tính hay học vấn
#-------------------------------------------------------------------
library(survival)
fit <- survfit(Surv(time, Churn) ~ 1, data = train)
print(fit)
## Call: survfit(formula = Surv(time, Churn) ~ 1, data = train)
##
## n events median 0.95LCL 0.95UCL
## 1488 733 11 11 11
# Kết quả của phân tích SA:
u <- summary(fit)
u
## Call: survfit(formula = Surv(time, Churn) ~ 1, data = train)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1488 14 0.991 0.00250 0.986 0.996
## 2 1457 42 0.962 0.00498 0.952 0.972
## 3 1405 35 0.938 0.00629 0.926 0.950
## 4 1352 23 0.922 0.00701 0.908 0.936
## 5 1322 41 0.894 0.00809 0.878 0.910
## 6 1245 43 0.863 0.00908 0.845 0.881
## 7 1151 56 0.821 0.01022 0.801 0.841
## 8 1021 72 0.763 0.01155 0.740 0.786
## 9 855 85 0.687 0.01301 0.662 0.713
## 10 667 123 0.560 0.01480 0.532 0.590
## 11 437 100 0.432 0.01603 0.402 0.465
## 12 224 99 0.241 0.01690 0.210 0.277
Có thể thấy tỉ lệ sống còn (trong tình huống này được hiểu là ở lại) survival rate (SR) giảm khi thời gian trôi qua: từ 0.991 ở tháng thứ nhất và đến tháng 12 chỉ còn 0.241.
Để tìm ra mức độ và diễn biến của tỉ lệ ở lại chúng ta nên hình ảnh hóa bằng đường cong sống còn (survival curve):
df1 <- data.frame(Rate = u$surv, Tenure = u$time)
df1 %>%
ggplot(aes(Tenure, Rate)) +
geom_line() +
geom_point(color = "red") +
scale_x_continuous(breaks = seq(1, 12, by = 1)) +
theme_ipsum(grid = "XY")
Có thể thấy tốc độ rời bỏ (độ dốc) bắt đầu tăng mạnh từ tháng thứ 6 trở đi. Điều này gợi ý rằng nhà cung cấp dịch vụ nên có các chính sách chăm sóc khách hàng cho những khách hàng còn ở lại từ thời điểm này.
Có thể cải tiến survival curve trên với gói survminer. Điểm khác là graph kiểu này nó bổ sung thêm khoảng tin cậy 95% cho survival rate:
library(survminer)
ggsurvplot(fit,
pval = TRUE,
conf.int = TRUE,
ggtheme = theme_minimal(),
palette = c("#E7B800"))
3.1.2 Tỉ lệ rời bỏ / ở lại có khác biệt giữa các nhóm giới tính không?
Kế tiếp chúng ta phân tích KM theo hai nhóm giới tính. Lúc này cần lưu ý rằng kết quả của survival rate theo thời gian được tính cho cả hai nhóm giới tính:
#-------------------------------------
# Thực hiện SA theo nhóm giới tính
#-------------------------------------
fit <- survfit(Surv(time, Churn) ~ Xgender, data = train)
print(fit)
## Call: survfit(formula = Surv(time, Churn) ~ Xgender, data = train)
##
## n events median 0.95LCL 0.95UCL
## Xgender=F 331 128 12 11 12
## Xgender=M 1157 605 11 11 11
summary(fit)
## Call: survfit(formula = Surv(time, Churn) ~ Xgender, data = train)
##
## Xgender=F
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 331 2 0.994 0.00426 0.986 1.000
## 2 324 7 0.972 0.00905 0.955 0.990
## 3 316 5 0.957 0.01122 0.935 0.979
## 5 307 8 0.932 0.01397 0.905 0.960
## 6 295 7 0.910 0.01594 0.879 0.942
## 7 273 13 0.867 0.01919 0.830 0.905
## 8 236 6 0.845 0.02070 0.805 0.886
## 9 200 19 0.764 0.02565 0.716 0.816
## 10 150 28 0.622 0.03204 0.562 0.688
## 11 86 14 0.521 0.03650 0.454 0.597
## 12 45 19 0.301 0.04374 0.226 0.400
##
## Xgender=M
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1157 12 0.990 0.00298 0.984 0.995
## 2 1133 35 0.959 0.00585 0.948 0.971
## 3 1089 30 0.933 0.00741 0.918 0.947
## 4 1045 23 0.912 0.00840 0.896 0.929
## 5 1015 33 0.882 0.00958 0.864 0.901
## 6 950 36 0.849 0.01072 0.828 0.870
## 7 878 43 0.807 0.01192 0.784 0.831
## 8 785 66 0.740 0.01353 0.713 0.767
## 9 655 66 0.665 0.01496 0.636 0.695
## 10 517 95 0.543 0.01666 0.511 0.576
## 11 351 86 0.410 0.01770 0.377 0.446
## 12 179 80 0.227 0.01811 0.194 0.265
summary(fit)$table
## records n.max n.start events *rmean *se(rmean) median 0.95LCL
## Xgender=F 331 331 331 128 10.34087 0.14545286 12 11
## Xgender=M 1157 1157 1157 605 9.68957 0.08920485 11 11
## 0.95UCL
## Xgender=F 12
## Xgender=M 11
Kế tiếp có thể hình ảnh hóa tốc độ và tỉ lệ rời bỏ cho hai nhóm giới tính. Có thể thấy ph nữ có tỉ lệ rời bỏ thấp hơn so với nhóm khách hàng nam giới:
ggsurvplot(fit,
pval = TRUE,
conf.int = TRUE,
# Bổ sung thêm Risk Table:
risk.table = TRUE,
# Thay bằng theme_minimal() của gói ggplot2:
ggtheme = theme_minimal())
Đồ thị trên cho thấy SR có khác biệt và chúng ta có thể sử dụng một kiểm định thống kê chính thức để trả lời xem sự khác biệt đó có ý nghĩa thống kê hay không bằng log-rank test:
survdiff(Surv(time, Churn) ~ Xgender, data = train)
## Call:
## survdiff(formula = Surv(time, Churn) ~ Xgender, data = train)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Xgender=F 331 128 162 7.25 10.9
## Xgender=M 1157 605 571 2.06 10.9
##
## Chisq= 10.9 on 1 degrees of freedom, p= 0.000962
Giá trị P-Value này rất bé (cũng được hiển thị ngay trên đồ thị mà chúng ta vừa có) cho thấy có bằng chứng thống kê đủ mạnh để kết luận rằng sự khác biệt về tỉ lệ rời bỏ giữa hai nhóm giới tính là có ý nghĩa.
3.1.3 Tỉ lệ rời bỏ / ở lại có khác biệt giữa các nhóm học vấn hay không?
Chúng ta cũng có thể thực hiện KM tương ứng với nhóm học vấn và kiểm định xem có sự khác biệt về tỉ lệ ở lại / rời bỏ giữa các nhóm học vấn:
#-------------------------------------
# Thực hiện SA theo nhóm học vấn
#-------------------------------------
# Thực hiện SA:
fit_educ <- survfit(Surv(time, Churn) ~ Xeducation + Xeducation, data = train)
fit_educ %>% summary()
## Call: survfit(formula = Surv(time, Churn) ~ Xeducation + Xeducation,
## data = train)
##
## Xeducation=Bachelor
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1204 14 0.988 0.00309 0.982 0.994
## 2 1175 32 0.961 0.00557 0.951 0.972
## 3 1135 29 0.937 0.00705 0.923 0.951
## 4 1090 22 0.918 0.00798 0.902 0.934
## 5 1061 37 0.886 0.00928 0.868 0.904
## 6 992 34 0.856 0.01032 0.836 0.876
## 7 918 42 0.816 0.01148 0.794 0.839
## 8 819 57 0.760 0.01291 0.735 0.785
## 9 682 63 0.689 0.01443 0.662 0.718
## 10 532 90 0.573 0.01641 0.542 0.606
## 11 362 81 0.445 0.01788 0.411 0.481
## 12 184 84 0.242 0.01900 0.207 0.282
##
## Xeducation=Doctorate
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 3 83 3 0.964 0.0205 0.925 1.000
## 6 77 5 0.901 0.0332 0.839 0.969
## 7 66 4 0.847 0.0409 0.770 0.931
## 8 59 3 0.804 0.0457 0.719 0.898
## 9 50 7 0.691 0.0557 0.590 0.809
## 10 38 7 0.564 0.0629 0.453 0.702
## 11 23 9 0.343 0.0690 0.231 0.509
## 12 11 4 0.218 0.0664 0.120 0.396
##
## Xeducation=Master
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2 198 10 0.949 0.0156 0.919 0.980
## 3 187 3 0.934 0.0176 0.900 0.969
## 4 184 1 0.929 0.0182 0.894 0.966
## 5 183 4 0.909 0.0205 0.870 0.950
## 6 176 4 0.888 0.0225 0.845 0.933
## 7 167 10 0.835 0.0267 0.784 0.889
## 8 143 12 0.765 0.0312 0.706 0.829
## 9 123 15 0.672 0.0355 0.606 0.745
## 10 97 26 0.492 0.0398 0.419 0.576
## 11 52 10 0.397 0.0419 0.323 0.488
## 12 29 11 0.246 0.0442 0.173 0.350
# Hình ảnh hóa:
ggsurvplot(fit_educ,
pval = TRUE,
conf.int = TRUE,
ggtheme = theme_minimal()) -> p
# Cách 1:
p
# Cách 2:
p$plot + facet_wrap(~ Xeducation)
# Giá trị P-value = 0.78 cho thấy không có sự khác biệt về tỉ lệ rời
# bỏ / ở lại giữa các nhóm học vấn bằng kiểm định log - rank test:
survdiff(Surv(time, Churn) ~ Xeducation, data = train)
## Call:
## survdiff(formula = Surv(time, Churn) ~ Xeducation, data = train)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Xeducation=Bachelor 1204 585 591.7 0.0751 0.4567
## Xeducation=Doctorate 86 42 41.2 0.0164 0.0203
## Xeducation=Master 198 106 100.2 0.3412 0.4634
##
## Chisq= 0.5 on 2 degrees of freedom, p= 0.776
3.1.4 Tỉ lệ rời bỏ / ở lại có khác biệt khi tính đến cả hai yếu tố là giới tính và học vấn hay không?
Tương tự chúng ta phân tích xa hơn cho cả: (1) giới tính, và (2) nhóm học vấn:
#----------------------------------------------
# Thực hiện SA theo giới tính và học vấn
#----------------------------------------------
fit2 <- survfit(Surv(time, Churn) ~ Xgender + Xeducation, data = train)
# Plot kiểu 1:
p1 <- ggsurvplot(fit2,
pval = TRUE,
conf.int = TRUE,
risk.table = TRUE,
ggtheme = theme_minimal())
p1
# Kiểu 2:
p1$plot + theme(legend.position = "right")
# Kiểu 3:
p2 <- ggsurvplot(fit2,
pval = TRUE,
conf.int = TRUE,
ggtheme = theme_minimal())
# Kiểu 4:
p2$plot +
theme_minimal() +
theme(legend.position = "top") +
facet_grid(Xgender ~ Xeducation)
# Sự khác biệt giữa các nhóm theo giới tính + học vấn là có ý nghĩa thống kê:
survdiff(Surv(time, Churn) ~ Xgender + Xeducation, data = train)
## Call:
## survdiff(formula = Surv(time, Churn) ~ Xgender + Xeducation,
## data = train)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Xgender=F, Xeducation=Bachelor 265 100 133.58 8.44202 12.1179
## Xgender=F, Xeducation=Doctorate 23 9 8.72 0.00876 0.0101
## Xgender=F, Xeducation=Master 43 19 19.99 0.04910 0.0579
## Xgender=M, Xeducation=Bachelor 939 485 458.09 1.58129 4.9575
## Xgender=M, Xeducation=Doctorate 63 33 32.45 0.00916 0.0113
## Xgender=M, Xeducation=Master 155 87 80.16 0.58299 0.7711
##
## Chisq= 12.5 on 5 degrees of freedom, p= 0.0282
# Hoặc hình ảnh hóa tỉ lệ rời bỏ tích lũy theo thời gian:
p3 <- ggsurvplot(fit2,
fun = "event",
conf.int = TRUE,
ggtheme = theme_bw())
p3$plot +
theme_bw() +
theme(legend.position = "right")+
facet_grid(Xeducation ~ Xgender)
3.1.5 Tỉ lệ rời bỏ / ở lại có khác biệt giữa các nhóm có mức độ thỏa mãn về dịch vụ - chăm sóc khách hàng mà họ được cung cấp hay không?
Mức độ thỏa mãn của khách hành được đo lường bằng biến Xsatisfaction theo thang đo Likert 5. Hiển nhiên là mức độ thỏa mãn của họ sẽ là nhân tố ảnh hưởng nhiều đến tỉ lệ ở lại / rời bỏ. Chúng ta có thể thực hiện phân tích KM để đánh giá tác động của mức độ hài lòng đến tỉ lệ rời bỏ / ở lại của khách hàng:
# Thực hiện mô hình KM và hình ảnh hóa:
fit3 <- survfit(Surv(time, Churn) ~ Xsatisfaction, data = train)
ggsurvplot(fit3,
pval = TRUE,
conf.int = TRUE,
risk.table = TRUE,
ggtheme = theme_minimal())
# Mức độ hài lòng là nhân tố ảnh hưởng đến tỉ lệ rời bỏ của khách hàng:
survdiff(Surv(time, Churn) ~ Xsatisfaction, data = train)
## Call:
## survdiff(formula = Surv(time, Churn) ~ Xsatisfaction, data = train)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Xsatisfaction=1 310 180 159 2.79 4.21
## Xsatisfaction=2 348 225 160 26.72 40.06
## Xsatisfaction=3 229 102 116 1.76 2.49
## Xsatisfaction=4 347 130 172 10.35 15.87
## Xsatisfaction=5 254 96 126 7.09 10.06
##
## Chisq= 57.1 on 4 degrees of freedom, p= 1.17e-11
3.1.6 Mức độ không hài lòng của khách hàng thể hiện qua số cuộc gọi phản hồi về dịch vụ và chăm sóc khách hàng có ảnh hưởng đến tỉ lệ rời bỏ / ở lại không?
Biến Xservice.calls đo lường số cuộc gọi của khách hàng (có thể là cuộc gọi phàn nàn về sản phẩm - dịch vụ được cung cấp) có thể là một nhân tố có ảnh hưởng đến tỉ lệ rời bỏ. Trước hết chúng ta tính toán phân bố của biến này:
ChurnStudy %>%
group_by(Xservice.calls) %>%
count() %>%
ggplot(aes(Xservice.calls, n)) +
geom_col() +
geom_text(aes(label = n), vjust = -0.5)
Có 1478 khách hàng không phản hồi. Có thể coi nhóm này là hài lòng với dịch vụ được cung cấp. Để phân tích có hay không sự khác biệt về những người có và không có cuộc gọi đến tỉ lệ rời bỏ có lẽ chúng ta nên dán lại nhãn cho biến này trước khi thực hiện KM:
# Dán lại nhãn:
train %<>% mutate(Calls = case_when(Xservice.calls == 0 ~ "Yes",
Xservice.calls != 0 ~ "No"))
# Thực hiện mô hình KM và hình ảnh hóa:
fit4 <- survfit(Surv(time, Churn) ~ Calls , data = train)
ggsurvplot(fit4,
pval = TRUE,
conf.int = TRUE,
risk.table = TRUE,
ggtheme = theme_minimal())
# Có bằng chứng thống kê cho thấy sự khác biệt về tỉ lệ rời bỏ là có ý nghĩa:
survdiff(Surv(time, Churn) ~ Calls, data = train)
## Call:
## survdiff(formula = Surv(time, Churn) ~ Calls, data = train)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Calls=No 354 222 180 9.88 15.5
## Calls=Yes 1134 511 553 3.21 15.5
##
## Chisq= 15.5 on 1 degrees of freedom, p= 8.34e-05
3.2 Hồi quy COx (Cox Regression)
Cách tiếp cận mà chúng ta vừa sử dụng ở trên (KM Model) có một nhược điểm là chúng ta không thể đưa vào mô hình các biến định lượng, ví dụ, biến Xmonthly.charges vì mô hình là một kiểu tiếp cận phi tham số (non-parametric). Để giải quyết nhược điểm của các tiếp cận này chúng ta có thể sử dụng hồi quy Cox (1972) hay còn gọi là Cox proportional-hazards model (chưa biết nên dịch là gì).
Dưới đây là hồi quy Cox với hai biến định lượng là Xmonthly.charges và Xpurch.last.month:
# Thực hiện hồi quy Cox:
res.cox1 <- coxph(Surv(time, Churn) ~ Xgender + Xmonthly.charges +
Xsatisfaction + Xpurch.last.month, data = train)
# Xem kết quả của hồi quy Cox:
res.cox1 %>% summary()
## Call:
## coxph(formula = Surv(time, Churn) ~ Xgender + Xmonthly.charges +
## Xsatisfaction + Xpurch.last.month, data = train)
##
## n= 1488, number of events= 733
##
## coef exp(coef) se(coef) z Pr(>|z|)
## XgenderM 0.2864258 1.3316594 0.0974086 2.940 0.00328 **
## Xmonthly.charges 0.0090558 1.0090969 0.0004061 22.300 < 2e-16 ***
## Xsatisfaction -0.0795540 0.9235282 0.0273389 -2.910 0.00362 **
## Xpurch.last.month -0.0104377 0.9896166 0.0236770 -0.441 0.65933
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## XgenderM 1.3317 0.7509 1.1002 1.6118
## Xmonthly.charges 1.0091 0.9910 1.0083 1.0099
## Xsatisfaction 0.9235 1.0828 0.8753 0.9744
## Xpurch.last.month 0.9896 1.0105 0.9447 1.0366
##
## Concordance= 0.761 (se = 0.013 )
## Rsquare= 0.317 (max possible= 0.998 )
## Likelihood ratio test= 567.4 on 4 df, p=0
## Wald test = 531.5 on 4 df, p=0
## Score (logrank) test = 608.7 on 4 df, p=0
Kết quả của hồi quy COx chỉ ra rằng:
Ngoại trừ Xpurch.last.month thì các biến khác là có ý nghĩa thống kê trong việc giải thích tỉ lệ ở lại - rời bỏ. Chú ý rằng Pr(>|z|) = coef / se(coef) đóng vai trò như giá trị của P-value trong OLS.
Căn cứ vào dấu có thể thấy Xsatisfaction (tức mức độ thỏa mãn của khách hàng) càng cao thì tỉ lệ rời bỏ càng thấp (tức ở lại càng cao).
Cuối cùng mức độ phù hợp của mô hình (còn gọi là Global statistical significance of the model) được đánh giá bằng cả ba thống kê là The likelihood-ratio test, Wald test, và score logrank có vai trò tương tự như thống kê F trong hồi quy OLS. Với số quan sát lớn thì các thống kê này sẽ không khác biệt nhiều (như trong tình huống của chúng ta) và tương ứng với P-value rất bé. Chú ý rằng với mẫu là bé thì Likelihood ratio test nên được sử dụng.
Chúng ta có thể ước lượng tỉ lệ ở lại / rời bỏ như theo thời gian như sau:
res.cox1 %>%
survfit() %>%
summary() -> u
u
## Call: survfit(formula = .)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1488 14 0.994 0.00172 0.990 0.997
## 2 1457 42 0.974 0.00354 0.967 0.981
## 3 1405 35 0.957 0.00461 0.948 0.966
## 4 1352 23 0.946 0.00524 0.936 0.956
## 5 1322 41 0.925 0.00625 0.913 0.937
## 6 1245 43 0.902 0.00726 0.888 0.917
## 7 1151 56 0.871 0.00854 0.854 0.888
## 8 1021 72 0.826 0.01018 0.806 0.846
## 9 855 85 0.766 0.01217 0.743 0.790
## 10 667 123 0.661 0.01520 0.632 0.692
## 11 437 100 0.544 0.01807 0.509 0.580
## 12 224 99 0.307 0.02269 0.265 0.355
# Và hình ảnh hóa tỉ lệ ở lại:
df3 <- data.frame(Rate = u$surv, Tenure = u$time)
df3 %>%
ggplot(aes(Tenure, Rate)) +
geom_line() +
geom_point(color = "red") +
scale_x_continuous(breaks = seq(1, 12, by = 1)) +
theme_ipsum(grid = "XY")
Căn cứ độ dốc có thể thấy từ tháng thứ 8 trở đi tỉ lệ rời bỏ bắt đầu tăng tốc bất kể nhóm giới tính của khách hàng.
Hình ảnh trên có thể được thể hiện theo một cách khác như sau:
ggsurvplot(survfit(res.cox1),
color = "#2E9FDF",
ggtheme = theme_minimal())
Nếu muốn hình ảnh hóa tỉ lệ rời bỏ này theo nhóm giới tính:
df_new <- data.frame(Xgender = c("M", "F"),
Xmonthly.charges = rep(mean(train$Xmonthly.charges), 2),
Xsatisfaction = rep(mean(train$Xsatisfaction), 2),
Xpurch.last.month = rep(mean(train$Xpurch.last.month), 2))
fit <- survfit(res.cox1, newdata = df_new)
ggsurvplot(fit,
legend.labs = c("Sex = M", "Sex = F"),
ggtheme = theme_minimal())
(Còn nữa ….)
