Survival Analysis
Apr 3, 2022
1 Import Data
library(magrittr)
library(tidyverse)
rm(list = ls())
ChurnStudy <- read.csv("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 : chr "Master's Degree" "Bachelor's Degree" "Bachelor's Degree" "Bachelor's Degree" ...
## $ Xgender : chr "F" "F" "M" "M" ...
## $ 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 102 Exploratory Data Analysis - EDA
# Danh gia muc do thieu cua du lieu
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# Remove na
ChurnStudy %<>% na.omit()
ChurnStudy %>% summary()## Churn Xeducation Xgender Xsatisfaction
## Min. :0.0 Length:1984 Length:1984 Min. :1.000
## 1st Qu.:0.0 Class :character Class :character 1st Qu.:2.000
## Median :0.5 Mode :character Mode :character 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# So luong nguoi mua theo hinh thuc tra cham
sum(ChurnStudy$Xmonthly.charges < 0)## [1] 7777*100/nrow(ChurnStudy)## [1] 3.881048# Xmonthly.charges
# Chuyen nguoi mua tra cham thanh tra nhanh
convert_to_pos <- function(x) {case_when(x < 0 ~ -x,
x > 0 ~ x)}
ChurnStudy %<>% mutate(Xmonthly.charges = convert_to_pos(Xmonthly.charges))# Hoc Van
hoc_van <- ChurnStudy$Xeducation %>% unique()
hoc_van## [1] "Master's Degree" "Bachelor's Degree" "Doctorate Degree"convert_hoc_van <- function(x) {case_when(x == hoc_van[1] ~ "Master",
x == hoc_van[2] ~ "Bachelor",
x == hoc_van[3] ~ "Doctorate")}
ChurnStudy %<>% mutate(Xeducation = convert_hoc_van(Xeducation))library(skimr)
skim(ChurnStudy)| Name | ChurnStudy |
| Number of rows | 1984 |
| Number of columns | 10 |
| _______________________ | |
| Column type frequency: | |
| character | 2 |
| numeric | 8 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| Xeducation | 0 | 1 | 6 | 9 | 0 | 3 | 0 |
| Xgender | 0 | 1 | 1 | 1 | 0 | 2 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Churn | 0 | 1 | 0.50 | 0.50 | 0.00 | 0.00 | 0.50 | 1.00 | 1.00 | ▇▁▁▁▇ |
| Xsatisfaction | 0 | 1 | 2.94 | 1.42 | 1.00 | 2.00 | 3.00 | 4.00 | 5.00 | ▇▇▅▇▆ |
| Xsatisfaction2 | 0 | 1 | 3.51 | 1.41 | 1.00 | 2.00 | 3.00 | 5.00 | 6.00 | ▇▆▆▆▂ |
| Xservice.calls | 0 | 1 | 0.47 | 0.99 | 0.00 | 0.00 | 0.00 | 1.00 | 5.00 | ▇▁▁▁▁ |
| Xtenure2 | 0 | 1 | 8.45 | 2.96 | 1.00 | 7.00 | 9.00 | 11.00 | 12.00 | ▂▁▃▅▇ |
| Xpurch.last.month | 0 | 1 | 1.48 | 1.46 | 0.00 | 1.00 | 1.00 | 1.00 | 10.00 | ▇▂▁▁▁ |
| Xmonthly.charges | 0 | 1 | 145.48 | 92.07 | 0.17 | 68.74 | 131.48 | 221.42 | 381.82 | ▇▇▆▅▁ |
| time | 0 | 1 | 8.45 | 2.96 | 1.00 | 7.00 | 9.00 | 11.00 | 12.00 | ▂▁▃▅▇ |
library(hrbrthemes)
theme_set(theme_ipsum())
# Education ~ Gender
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 = "Edu level") +
theme(legend.position = "top")
library(hrbrthemes)
theme_set(theme_ipsum())
# Churn ~ Edu
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 = "Edu level") +
theme(legend.position = "top")
# Churn ~ Edu
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 = "Edu level") +
theme(legend.position = "top")
# Churn ~ Gender ~ Edu
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 = "Edu level") +
theme(legend.position = "top") +
facet_wrap(~ Xgender)
# Churn ~ Gender ~ Edu ~ Xmonthly.charges
# Số tiền mua hàng trung bình theo từng nhóm giới tính và học vấn
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 = "Edu Level") ## Warning: `summarise_each_()` was deprecated in dplyr 0.7.0.
## Please use `across()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
# Churn ~ Gender ~ Edu ~ Xpurch.last.month
# 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 = "Edu Level") 
3 SA Models
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)
# Churn ~ Gender ~ Edu ~ Xpurch.last.month
# Số lần mua hàng trung bình mỗi tháng thì nam giới cũng cao hơn:
# Lieu 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?.
# Du lieu huan luyen
set.seed(1020)
train <- ChurnStudy %>% sample_n(round(0.75*nrow(.)))
# Du lieu kiem dinh
test <- dplyr::setdiff(ChurnStudy, train)3.1.1 Evaluation the Churn/ Stay Rate: Whole data
3.1.1.1 Survival test
library(survival)
fit <- survfit(Surv(time, Churn) ~ 1, data = train) # time, Churn <- 2 variables
print(fit)## Call: survfit(formula = Surv(time, Churn) ~ 1, data = train)
##
## n events median 0.95LCL 0.95UCL
## [1,] 1488 752 11 11 11u <- 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 11 0.993 0.00222 0.988 0.997
## 2 1458 52 0.957 0.00528 0.947 0.968
## 3 1397 30 0.937 0.00636 0.924 0.949
## 4 1352 25 0.919 0.00712 0.905 0.933
## 5 1322 43 0.889 0.00822 0.873 0.906
## 6 1241 45 0.857 0.00922 0.839 0.875
## 7 1144 55 0.816 0.01032 0.796 0.836
## 8 1015 72 0.758 0.01162 0.736 0.781
## 9 852 84 0.683 0.01303 0.658 0.709
## 10 679 122 0.561 0.01468 0.533 0.590
## 11 444 99 0.436 0.01590 0.405 0.468
## 12 235 114 0.224 0.01639 0.194 0.259Có 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.
3.1.1.2 Plot
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")
Note: tốc độ rời bỏ (độ dốc) bắt đầu tăng mạnh từ tháng thứ 6 trở đi
# Plot + khoảng tin cậy 95% cho survival rate
library(survminer)
ggsurvplot(fit,
pval = TRUE,
conf.int = TRUE,
ggtheme = theme_minimal(),
palette = c("#E7B800"))## Warning in .pvalue(fit, data = data, method = method, pval = pval, pval.coord = pval.coord, : There are no survival curves to be compared.
## This is a null model.
### Evaluation the Churn/ Stay Rate btw Gender #### Survival test
fit <- survfit(Surv(time, Churn) ~ Xgender, data = train) # replace 1 by Xgender
print(fit)## Call: survfit(formula = Surv(time, Churn) ~ Xgender, data = train)
##
## n events median 0.95LCL 0.95UCL
## Xgender=F 337 133 12 11 12
## Xgender=M 1151 619 11 11 11summary(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 337 1 0.997 0.00296 0.991 1.000
## 2 328 11 0.964 0.01032 0.944 0.984
## 3 316 4 0.951 0.01185 0.928 0.975
## 4 309 2 0.945 0.01255 0.921 0.970
## 5 306 9 0.917 0.01522 0.888 0.948
## 6 294 7 0.896 0.01695 0.863 0.929
## 7 274 11 0.860 0.01943 0.822 0.899
## 8 241 8 0.831 0.02124 0.790 0.874
## 9 199 22 0.739 0.02643 0.689 0.793
## 10 148 28 0.599 0.03202 0.540 0.666
## 11 87 12 0.517 0.03540 0.452 0.591
## 12 47 18 0.319 0.04265 0.245 0.414
##
## Xgender=M
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1151 10 0.991 0.00274 0.986 0.997
## 2 1130 41 0.955 0.00611 0.943 0.967
## 3 1081 26 0.932 0.00744 0.918 0.947
## 4 1043 23 0.912 0.00842 0.895 0.928
## 5 1016 34 0.881 0.00963 0.863 0.900
## 6 947 38 0.846 0.01082 0.825 0.867
## 7 870 44 0.803 0.01204 0.780 0.827
## 8 774 64 0.737 0.01361 0.711 0.764
## 9 653 62 0.667 0.01494 0.638 0.697
## 10 531 94 0.549 0.01653 0.517 0.582
## 11 357 87 0.415 0.01765 0.382 0.451
## 12 188 96 0.203 0.01742 0.172 0.240summary(fit)$table ## records n.max n.start events rmean se(rmean) median 0.95LCL
## Xgender=F 337 337 337 133 10.216340 0.15134105 12 11
## Xgender=M 1151 1151 1151 619 9.688488 0.08983954 11 11
## 0.95UCL
## Xgender=F 12
## Xgender=M 113.1.1.3 Plot
# Tốc độ và tỉ lệ rời bỏ cho hai nhóm giới tính
ggsurvplot(fit,
pval = TRUE,
conf.int = TRUE,
risk.table = TRUE,
ggtheme = theme_minimal())
Note: Ty lệ sống còn (trong tình huống này được hiểu là ở lại) survival rate (SR). KH nữ có tỉ lệ rời bỏ thấp hơn so với nhóm khách hàng nam giới
3.1.1.4 Log-rank test (statistical significant)
# Su dung kiem dinh thong ke de kiem tra 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 337 133 166 6.55 9.95
## Xgender=M 1151 619 586 1.86 9.95
##
## Chisq= 10 on 1 degrees of freedom, p= 0.002print(survdiff)## function (formula, data, subset, na.action, rho = 0, timefix = TRUE)
## {
## call <- match.call()
## m <- match.call(expand.dots = FALSE)
## m$rho <- NULL
## if (!inherits(formula, "formula"))
## stop("The 'formula' argument is not a formula")
## Terms <- if (missing(data))
## terms(formula, "strata")
## else terms(formula, "strata", data = data)
## m$formula <- Terms
## m[[1L]] <- quote(stats::model.frame)
## m <- eval(m, parent.frame())
## y <- model.extract(m, "response")
## if (!inherits(y, "Surv"))
## stop("Response must be a survival object")
## if (attr(y, "type") != "right")
## stop("Right censored data only")
## ny <- ncol(y)
## n <- nrow(y)
## if (!is.logical(timefix) || length(timefix) > 1)
## stop("invalid value for timefix option")
## if (timefix)
## y <- aeqSurv(y)
## offset <- attr(Terms, "offset")
## if (!is.null(offset)) {
## offset <- as.numeric(m[[offset]])
## if (length(attr(Terms, "factors")) > 0)
## stop("Cannot have both an offset and groups")
## if (any(offset < 0 | offset > 1))
## stop("The offset must be a survival probability")
## expected <- sum(-log(offset))
## observed <- sum(y[, ny])
## if (rho != 0) {
## num <- sum(1/rho - ((1/rho + y[, ny]) * offset^rho))
## var <- sum(1 - offset^(2 * rho))/(2 * rho)
## }
## else {
## var <- sum(-log(offset))
## num <- var - observed
## }
## chi <- num * num/var
## rval <- list(n = n, obs = observed, exp = expected, var = var,
## chisq = chi)
## }
## else {
## strats <- attr(Terms, "specials")$strata
## if (length(strats)) {
## temp <- untangle.specials(Terms, "strata", 1)
## dropx <- temp$terms
## if (length(temp$vars) == 1)
## strata.keep <- m[[temp$vars]]
## else strata.keep <- strata(m[, temp$vars], shortlabel = TRUE)
## }
## else strata.keep <- rep(1, nrow(m))
## if (length(strats))
## ll <- attr(Terms[-dropx], "term.labels")
## else ll <- attr(Terms, "term.labels")
## if (length(ll) == 0)
## stop("No groups to test")
## else groups <- strata(m[ll])
## fit <- survdiff.fit(y, groups, strata.keep, rho)
## if (is.matrix(fit$observed)) {
## otmp <- apply(fit$observed, 1, sum)
## etmp <- apply(fit$expected, 1, sum)
## }
## else {
## otmp <- fit$observed
## etmp <- fit$expected
## }
## df <- (etmp > 0)
## if (sum(df) < 2)
## chi <- 0
## else {
## temp2 <- ((otmp - etmp)[df])[-1]
## vv <- (fit$var[df, df])[-1, -1, drop = FALSE]
## chi <- sum(solve(vv, temp2) * temp2)
## }
## rval <- list(n = table(groups), obs = fit$observed, exp = fit$expected,
## var = fit$var, chisq = chi)
## if (length(strats))
## rval$strata <- table(strata.keep)
## }
## na.action <- attr(m, "na.action")
## if (length(na.action))
## rval$na.action <- na.action
## rval$call <- call
## class(rval) <- "survdiff"
## rval
## }
## <bytecode: 0x7fad2304c6d8>
## <environment: namespace:survival>p-vale rat be, p = 0.000962 < 0.05 . bac bo HO: khong co su khac biet giua nam-nu. Vay, 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.2 Evaluation the Churn/ Stay Rate btw Edu level
3.1.2.1 Survival test
fit_educ <- survfit(Surv(time, Churn) ~ Xeducation, data = train) # replace 1 by Xgender
fit_educ %>% summary()## Call: survfit(formula = Surv(time, Churn) ~ Xeducation, data = train)
##
## Xeducation=Bachelor
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 1186 11 0.991 0.00278 0.985 0.996
## 2 1158 38 0.958 0.00584 0.947 0.970
## 3 1113 25 0.937 0.00712 0.923 0.951
## 4 1076 25 0.915 0.00818 0.899 0.931
## 5 1046 38 0.882 0.00950 0.863 0.900
## 6 975 37 0.848 0.01061 0.828 0.869
## 7 898 44 0.807 0.01180 0.784 0.830
## 8 797 60 0.746 0.01326 0.720 0.772
## 9 665 59 0.680 0.01462 0.652 0.709
## 10 536 94 0.561 0.01643 0.529 0.594
## 11 359 78 0.439 0.01773 0.405 0.475
## 12 189 91 0.228 0.01841 0.194 0.267
##
## Xeducation=Doctorate
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2 86 3 0.965 0.0198 0.927 1.000
## 3 82 2 0.942 0.0254 0.893 0.993
## 6 77 4 0.893 0.0338 0.829 0.962
## 7 68 2 0.866 0.0376 0.796 0.943
## 8 63 2 0.839 0.0411 0.762 0.924
## 9 53 5 0.760 0.0502 0.667 0.865
## 10 40 5 0.665 0.0592 0.558 0.792
## 11 24 9 0.415 0.0754 0.291 0.593
## 12 11 7 0.151 0.0662 0.064 0.357
##
## Xeducation=Master
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2 214 11 0.949 0.0151 0.919 0.979
## 3 202 3 0.935 0.0169 0.902 0.968
## 5 198 5 0.911 0.0195 0.873 0.950
## 6 189 4 0.892 0.0213 0.851 0.934
## 7 178 9 0.847 0.0250 0.799 0.897
## 8 155 10 0.792 0.0287 0.738 0.850
## 9 134 20 0.674 0.0345 0.609 0.745
## 10 103 23 0.523 0.0385 0.453 0.604
## 11 61 12 0.420 0.0408 0.347 0.508
## 12 35 16 0.228 0.0418 0.159 0.3273.1.2.2 Plot
# option 1
ggsurvplot(fit_educ,
pval = TRUE,
conf.int = TRUE,
ggtheme = theme_minimal()) -> p
p
# option 2
p$plot + facet_wrap(~ Xeducation)
##### Log-rank test (statistical significant)
# 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 1186 600 596.6 0.01991 0.11429
## Xeducation=Doctorate 88 39 43.1 0.38790 0.48351
## Xeducation=Master 214 113 112.4 0.00367 0.00513
##
## Chisq= 0.5 on 2 degrees of freedom, p= 0.83.1.3 Evaluation the Churn/ Stay Rate btw Gender-Edu
3.1.3.1 Survival test
fit2 <- survfit(Surv(time, Churn) ~ Xgender + Xeducation, data = train)3.1.3.2 Plot
# option 1
p1 <- ggsurvplot(fit2,
pval = TRUE,
conf.int = TRUE,
risk.table = TRUE,
ggtheme = theme_minimal())
p1
# option 2
p1$plot + theme(legend.position = "right")
# option 3
p2 <- ggsurvplot(fit2,
pval = TRUE,
conf.int = TRUE,
ggtheme = theme_minimal())# option 4
p2$plot +
theme_minimal() +
theme(legend.position = "top") +
facet_grid(Xgender ~ Xeducation)
3.1.3.3 Log-rank test (statistical significant)
# 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 271 104 136.77 7.852341 1.14e+01
## Xgender=F, Xeducation=Doctorate 21 7 7.96 0.115979 1.35e-01
## Xgender=F, Xeducation=Master 45 22 21.25 0.026496 3.17e-02
## Xgender=M, Xeducation=Bachelor 915 496 459.78 2.852909 8.73e+00
## Xgender=M, Xeducation=Doctorate 67 32 35.13 0.278428 3.45e-01
## Xgender=M, Xeducation=Master 169 91 91.11 0.000128 1.75e-04
##
## Chisq= 13.2 on 5 degrees of freedom, p= 0.02p3 <- ggsurvplot(fit2,
fun = "event",
conf.int = TRUE,
ggtheme = theme_bw())
p3$plot +
theme_bw() +
theme(legend.position = "right")+
facet_grid(Xeducation ~ Xgender)
### Evaluation the Churn/ Stay Rate btw group Xsatisfaction
#### Survival test
# Churn / stay rate giua cac 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?
fit3 <- survfit(Surv(time, Churn) ~ Xsatisfaction, data = train)3.1.3.4 Plot
ggsurvplot(fit3,
pval = TRUE,
conf.int = TRUE,
risk.table = TRUE,
ggtheme = theme_minimal())
3.1.3.5 Log-rank test (statistically significant variable)
# 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 322 191 166 3.926 5.989
## Xsatisfaction=2 325 204 160 12.338 18.659
## Xsatisfaction=3 214 101 108 0.496 0.695
## Xsatisfaction=4 352 144 179 6.964 10.866
## Xsatisfaction=5 275 112 139 5.314 7.766
##
## Chisq= 34.6 on 4 degrees of freedom, p= 6e-073.1.4 Evaluation the Churn/ Stay Rate btw Xservice.calls
3.1.4.1 Descriptive
ChurnStudy %>%
group_by(Xservice.calls) %>%
count() %>%
ggplot(aes(Xservice.calls, n)) +
geom_col() +
geom_text(aes(label = n), vjust = -0.5)
3.1.4.2 Survival test
# Churn / stay rate giua cac nhóm có mức độ thỏa mãn về số cuộc gọi phản hồi về dịch vụ và chăm sóc
# Dán lại nhãn:
train %<>% mutate(Calls = case_when(Xservice.calls == 0 ~ "Yes",
Xservice.calls != 0 ~ "No"))
fit4 <- survfit(Surv(time, Churn) ~ Calls , data = train)3.1.4.3 Plot
ggsurvplot(fit4,
pval = TRUE,
conf.int = TRUE,
risk.table = TRUE,
ggtheme = theme_minimal())
3.1.4.4 Log-rank test (statistically significant variable)
# 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) ~ 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 373 236 195 8.53 13.8
## Calls=Yes 1115 516 557 2.99 13.8
##
## Chisq= 13.8 on 1 degrees of freedom, p= 2e-043.2 Cox Regression
3.2.1 Cox proportional-hazards model
# Hồi quy Cox với hai biến định lượng là Xmonthly.charges và Xpurch.last.month
res.cox1 <- coxph(Surv(time, Churn) ~ Xgender + Xmonthly.charges +
Xsatisfaction + Xpurch.last.month, data = train)
res.cox1 %>% summary()## Call:
## coxph(formula = Surv(time, Churn) ~ Xgender + Xmonthly.charges +
## Xsatisfaction + Xpurch.last.month, data = train)
##
## n= 1488, number of events= 752
##
## coef exp(coef) se(coef) z Pr(>|z|)
## XgenderM 0.1679329 1.1828573 0.0957866 1.753 0.07957 .
## Xmonthly.charges 0.0088192 1.0088582 0.0004082 21.607 < 2e-16 ***
## Xsatisfaction -0.0773196 0.9255940 0.0260482 -2.968 0.00299 **
## Xpurch.last.month -0.0013978 0.9986031 0.0226881 -0.062 0.95087
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## XgenderM 1.1829 0.8454 0.9804 1.4271
## Xmonthly.charges 1.0089 0.9912 1.0081 1.0097
## Xsatisfaction 0.9256 1.0804 0.8795 0.9741
## Xpurch.last.month 0.9986 1.0014 0.9552 1.0440
##
## Concordance= 0.746 (se = 0.01 )
## Likelihood ratio test= 524.5 on 4 df, p=<2e-16
## Wald test = 493 on 4 df, p=<2e-16
## Score (logrank) test = 555.3 on 4 df, p=<2e-16Note: 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.
3.2.2 Testing proportional Hazards assumption
test.ph <- cox.zph(res.cox1)
par(mfrow = c(2,2))
plot(test.ph)
ggcoxzph(test.ph)
p = 0.1483 > 0.05 -> insignificant, the Schoenfeld residuals are independent of time
3.2.3 Survival rate
# uoc luong ty le o lai/ roi bo theo thoi gian
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 11 0.995 0.00143 0.993 0.998
## 2 1458 52 0.973 0.00400 0.965 0.981
## 3 1397 30 0.960 0.00523 0.950 0.970
## 4 1352 25 0.949 0.00622 0.937 0.961
## 5 1322 43 0.929 0.00788 0.913 0.944
## 6 1241 45 0.907 0.00963 0.888 0.926
## 7 1144 55 0.878 0.01181 0.856 0.902
## 8 1015 72 0.838 0.01477 0.809 0.867
## 9 852 84 0.784 0.01847 0.749 0.821
## 10 679 122 0.689 0.02431 0.643 0.738
## 11 444 99 0.583 0.02971 0.528 0.644
## 12 235 114 0.335 0.03623 0.271 0.4143.2.4 Plot
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")
ggsurvplot(survfit(res.cox1),
data = train,
color = "#2E9FDF",
ggtheme = theme_minimal())## Warning: Now, to change color palette, use the argument palette= '#2E9FDF'
## instead of color = '#2E9FDF'
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,
data = train,
legend.labs = c("Sex = M", "Sex = F"),
ggtheme = theme_minimal())## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
4 References
https://rpubs.com/chidungkt/339028 https://bvag.com.vn/wp-content/uploads/2013/01/k2_attachments_PHAN-TICH-SONG-SOT.pdf http://rstudio-pubs-static.s3.amazonaws.com/357614_0475597715cc46df9ae7aeb369fc6731.html#1_survival_analysis:_time_-_event_analysis https://epirhandbook.com/vn/survival-analysis.html