Ngôn ngữ lập trình trong phân tích dữ liệu

CHƯƠNG 1 - PHÂN TÍCH VỀ BỘ DỮ LIỆU LC LOAN

Bộ dữ liệu LC Loan đóng vai trò quan trọng trong việc phân tích rủi ro và hiệu quả hoạt động tín dụng. Thông qua việc khai thác thông tin chi tiết về đặc điểm khách hàng, cấu trúc khoản vay và tình trạng trả nợ, dữ liệu này giúp đánh giá chất lượng tín dụng cũng như khả năng trả nợ của người vay. Việc phân tích LC Loan không chỉ hỗ trợ xây dựng các mô hình dự báo rủi ro tín dụng mà còn cung cấp cơ sở thực tiễn cho việc hoạch định chính sách tín dụng và quản trị danh mục cho vay hiệu quả.

Giới thiệu bộ dữ liệu

Phần này thiết lập môi trường làm việc, tải dữ liệu thô và xây dựng từ điển dữ liệu ban đầu, đảm bảo tính tái lập và minh bạch của quy trình phân tích.

Tải hàm chức năng

log_info <- function(..., .sep = "") {
msg <- paste(..., sep = .sep)
cat(sprintf("[INFO] %s | %s\n", format(Sys.time(), "%Y-%m-%d %H:%M:%S"), msg))
}

safe_parse_num <- function(x) {
if (is.numeric(x)) as.numeric(x) else readr::parse_number(
as.character(x),
locale = readr::locale(grouping_mark = ",", decimal_mark = ".")
)
}

winsorize_vec <- function(x, probs = c(.01,.99)) {
if(!is.numeric(x)) return(x)
qs <- stats::quantile(x, probs, na.rm=TRUE)
pmin(pmax(x, qs[1]), qs[2])
}

mode_impute <- function(x) {
if (is.numeric(x)) return(x)
ux <- unique(x[!is.na(x)])
if (length(ux)==0) return(x)
tab <- sort(table(x), decreasing = TRUE)
x[is.na(x)] <- names(tab)[1]
x
}

emp_to_years <- function(s) {
s <- as.character(s)
dplyr::case_when(
is.na(s) ~ NA_real_,
s %in% c("10+ years","10+ Years","10+ năm","10+ năm ","10+yrs") ~ 10,
s %in% c("< 1 year","<1 year","< 1 năm","<1y") ~ 0.5,
TRUE ~ suppressWarnings(readr::parse_number(s))
)
}

need_cols <- function(df, cols) all(cols %in% names(df))

parse_issue_date <- function(x) {
if (inherits(x, "Date")) return(x)
as.Date(suppressWarnings(lubridate::parse_date_time(
as.character(x),
orders = c("my","bY","Y-m","Y-m-d","d/m/Y","m/d/Y"),
quiet = TRUE
)))
}

coalesce_issue_date <- function(df) {
candidates <- c("issue_d","issue_date","earliest_cr_line","last_pymnt_d","last_credit_pull_d")
out <- rep(as.Date(NA), nrow(df))
for (nm in candidates) if (nm %in% names(df)) {
d <- parse_issue_date(df[[nm]])
fill_idx <- is.na(out) & !is.na(d)
if (any(fill_idx)) out[fill_idx] <- d[fill_idx]
}
out
}

sample_rows <- function(df, max_n = 200000, seed = 3698) {
total <- nrow(df)
size  <- min(total, max_n)
if (!is.finite(size) || size <= 0) return(df[0, , drop = FALSE])
set.seed(seed)
dplyr::slice_sample(df, n = size)
}

Đoạn mã trên định nghĩa một tập hợp các hàm tiện ích (utility functions) phục vụ cho xử lý và làm sạch dữ liệu trước khi phân tích thống kê. Cụ thể:
log_info() dùng để ghi nhật ký thời gian và thông báo tiến trình, giúp theo dõi quy trình phân tích.
safe_parse_num() chuẩn hóa dữ liệu số có dấu phân tách hàng nghìn hoặc thập phân, đảm bảo tính nhất quán của biến định lượng.
winsorize_vec() thực hiện winsorization nhằm loại bỏ ảnh hưởng của các giá trị ngoại lai, tăng độ ổn định thống kê của mô hình.
mode_impute() thay thế giá trị thiếu của biến định tính bằng giá trị mode (tần suất cao nhất), hạn chế mất mát thông tin.
emp_to_years() quy đổi dữ liệu kinh nghiệm làm việc từ chuỗi ký tự sang đơn vị năm, giúp định lượng hóa biến định tính.
need_cols() kiểm tra sự tồn tại của các cột bắt buộc trong bộ dữ liệu, đảm bảo tính đầy đủ của biến.
parse_issue_date() và coalesce_issue_date() chuẩn hóa và gộp các cột ngày tháng khác nhau về định dạng ngày thống nhất, phục vụ cho phân tích chuỗi thời gian.
sample_rows() lấy mẫu ngẫu nhiên có giới hạn, giúp giảm tải tính toán mà vẫn giữ được tính đại diện thống kê.
Về ý nghĩa thống kê, các hàm này giúp làm sạch, chuẩn hóa và kiểm soát dữ liệu đầu vào, từ đó giảm sai số, đảm bảo tính chính xác, tin cậy và hợp lệ của kết quả phân tích hồi quy và mô hình hóa sau này. ### Tải dữ liệu thô và đánh giá sơ bộ

library(dplyr)
data_path <- if (exists("params") && !is.null(params$data_path)) params$data_path else "LC_loan.csv"
if (!file.exists(data_path)) stop("Không tìm thấy file: ", data_path)
log_info("Đọc dữ liệu từ: ", data_path)

## [INFO] 2025-11-10 11:43:28 | Đọc dữ liệu từ: LC_loan.csv

loans_raw <- readr::read_csv(data_path, show_col_types = FALSE) |>
janitor::clean_names()

n_obs  <- nrow(loans_raw)
n_vars <- ncol(loans_raw)
dup_id_n <- if ("id" %in% names(loans_raw)) {
loans_raw |> dplyr::count(id, name="n") |> dplyr::summarise(dups=sum(pmax(n-1,0))) |> dplyr::pull(dups)
} else NA_integer_

missing_by_col <- sapply(loans_raw, function(x) sum(is.na(x)))
missing_anyrow <- sum(!stats::complete.cases(loans_raw))

quality_tbl <- tibble::tibble(
column = names(loans_raw),
missing_n = as.integer(missing_by_col),
missing_pct = round(100 * missing_by_col / max(n_obs,1), 2)
) |> arrange(desc(missing_n))

Đoạn mã trên thực hiện bước đầu tiên trong quy trình phân tích dữ liệu: tải và kiểm tra chất lượng dữ liệu gốc. Về kỹ thuật, code thực hiện các thao tác sau:
Xác định đường dẫn file dữ liệu (data_path) và dừng chương trình nếu không tìm thấy tệp CSV đầu vào.
Ghi nhật ký tiến trình đọc dữ liệu bằng hàm log_info().
Đọc file CSV bằng read_csv() và chuẩn hóa tên biến (về dạng chữ thường, không dấu, không khoảng trắng) bằng clean_names() – giúp dễ thao tác và nhất quán khi gọi biến sau này.
Tính kích thước dữ liệu gồm số quan sát (n_obs) và số biến (n_vars).
Kiểm tra trùng ID (nếu có cột id) để phát hiện dữ liệu lặp lại, từ đó đảm bảo tính duy nhất của bản ghi.
Đếm số lượng và tỷ lệ giá trị thiếu theo từng cột, đồng thời tạo bảng tổng hợp quality_tbl giúp đánh giá chất lượng dữ liệu.
Về ý nghĩa thống kê, bước này có vai trò xác định tính toàn vẹn và độ tin cậy của dữ liệu đầu vào. Việc phát hiện thiếu dữ liệu, trùng lặp hoặc định dạng không nhất quán giúp người phân tích có căn cứ xử lý trước khi tiến hành các bước mô hình hóa, tránh sai lệch kết quả thống kê và đảm bảo tính hợp lệ của các suy luận sau này.

Định nghĩa siêu dữ liệu và công cụ khám phá dữ liệu

Chunk này là bước khởi tạo, định nghĩa siêu dữ liệu (var_meta) mô tả ý nghĩa kinh tế lượng của từng biến và các hàm tiện ích cần thiết để tự động tóm tắt dữ liệu.

library(tidyverse)
library(DT)

#--- 1) Từ điển mô tả/ý nghĩa ----
var_meta <- tibble::tribble(
  ~var,              ~label,                     ~unit,         ~definition,                                                                 ~construction,                                                                 ~economic_role,                                                                 ~expected_sign_on_default, ~caveats,
  "id",              "Loan ID",                  NA,            "Định danh khoản vay duy nhất.",                                            "Từ nguồn gốc dữ liệu, không biến đổi.",                                       "Khoá kỹ thuật, không dùng làm biến giải thích.",                               "0",                     "Không đưa vào mô hình.",
  "issue_d",         "Issue date",               NA,            "Tháng/năm giải ngân khoản vay.",                                            "Chuẩn hoá dạng date hoặc YM.",                                                "Kiểm soát thời kỳ/chu kỳ kinh tế.",                                           "?",                     "Hiệu ứng mùa vụ/chu kỳ, cần tạo biến năm/quý.",
  "revenue",         "Annual income",            "USD/year",    "Thu nhập năm của người vay.",                                               "Từ `annual_inc` hoặc trường tương ứng; có thể winsor để giảm ngoại lệ.",       "Khả năng trả nợ; thu nhập ↑ → rủi ro ↓.",                                      "-",                     "Thường lệch phải mạnh; cân nhắc log-transform.",
  "dti_n",           "Debt-to-income (DTI)",     "%",           "Tỷ lệ nợ/thu nhập tại thời điểm vay.",                                       "Chuẩn hoá về % hoặc [0–1]; kiểm tra outlier (DTI>60%).",                       "Gánh nặng nợ; DTI ↑ → rủi ro ↑.",                                             "+",                     "Đảm bảo đơn vị thống nhất toàn bộ pipeline.",
  "loan_amnt",       "Loan amount",              "USD",         "Số tiền vay phê duyệt.",                                                     "Đọc trực tiếp; có thể winsor theo phân vị.",                                    "Quy mô khoản vay; lớn hơn có thể ↑ rủi ro nếu năng lực trả nợ không đổi.",    "+/0",                   "Tác động phụ thuộc thu nhập/tài sản; tương tác với revenue.",
  "fico_n",          "FICO score",               "score",       "Điểm tín dụng tổng hợp (300–850).",                                         "Tách/chuẩn hoá từ trường FICO, lấy trung vị nếu là khoảng.",                   "Chất lượng tín dụng; FICO ↑ → rủi ro ↓.",                                      "-",                     "Có thể phi tuyến; dùng spline/bin hoặc bậc hai.",
  "experience_c",    "Has employment info",      "0/1",         "Cờ người vay có thông tin kinh nghiệm làm việc.",                           "Từ trường kinh nghiệm; 1 nếu có thông tin hợp lệ.",                             "Proxy về minh bạch hồ sơ.",                                                    "-",                     "Dễ bị lệch do thiếu dữ liệu.",
  "emp_length",      "Employment length",        "years",       "Số năm làm việc (0–10+, đã chuẩn hoá).",                                     "Chuẩn hoá text ('10+ years', '3 years', '<1 year') → số năm.",                 "Ổn định thu nhập; lâu năm → rủi ro ↓.",                                        "-",                     "Không tuyến tính; nên bin theo dải năm.",
  "purpose",         "Loan purpose",             "category",    "Mục đích vay do khách hàng khai báo.",                                       "Chuẩn hoá factor; gom nhóm hiếm vào 'other'.",                                  "Khác biệt rủi ro theo mục đích (small_business ↑, credit_card/refi ↓).",       "mixed",                 "Nguy cơ tự khai sai mục đích.",
  "home_ownership_n","Home ownership",           "category",    "Tình trạng sở hữu nhà: RENT/MORTGAGE/OWN…",                                  "Chuẩn hoá factor.",                                                            "Tài sản thế chấp & ổn định cư trú; OWN/MORTGAGE thường rủi ro ↓ hơn RENT.",    "-",                     "MORTGAGE kèm đòn bẩy; bối cảnh thị trường quan trọng.",
  "addr_state",      "Borrower state",           "category",    "Bang/khu vực cư trú.",                                                       "Chuẩn hoá hai chữ cái.",                                                       "Hiệu ứng địa lý/lao động/quy định.",                                           "?",                     "Dùng fixed-effects hoặc pooling theo vùng.",
  "zip_code",        "Zip (3 digits)",           "category",    "ZIP rút gọn (ẩn danh).",                                                     "Giữ 3 ký tự đầu; ẩn danh.",                                                    "Địa lý vi mô; thay `addr_state` trong vài mô hình.",                           "?",                     "Dễ trùng thông tin với state; cân nhắc một trong hai.",
  "default",         "Default flag",             "0/1",         "Biến đích: vỡ nợ/không vỡ nợ.",                                              "Chuẩn hoá nhị phân 0/1; mapping từ nhãn văn bản nếu cần.",                      "Đích mô hình.",                                                                  NA,                      "Cân bằng lớp (class imbalance) & định nghĩa thời điểm default.",
  "title",           "User loan title",          "text",        "Tiêu đề do người vay nhập.",                                                 "Chuẩn hoá text; có thể dùng sentiment/topic.",                                  "Tín hiệu định tính (noisy).",                                                   "mixed",                 "Nhiễu cao; nên regularize mạnh.",
  "desc",            "User loan description",    "text",        "Mô tả tự do của người vay.",                                                 "Tiền xử lý NLP, lexicon/sentiment/embeddings.",                                 "Tín hiệu định tính; cảm xúc/độ tin cậy.",                                       "mixed",                 "Thiếu nhiều/thiên lệch; cẩn trọng khi suy diễn."
)

#--- 2) Hàm tiện ích an toàn ----
pNA <- function(x) mean(is.na(x))*100
nuniq <- function(x) dplyr::n_distinct(x, na.rm = TRUE)
pcts  <- function(x) stats::quantile(x, probs = c(.01,.05,.25,.5,.75,.95,.99), na.rm = TRUE, names = TRUE)

# Tóm tắt top levels cho biến phân loại
top_levels <- function(x, k = 5) {
  if (!is.character(x) && !is.factor(x)) return(NA_character_)
  tb <- sort(table(x), decreasing = TRUE)
  tb <- head(tb, k)
  paste0(names(tb), " (", sprintf("%.1f", 100*as.numeric(tb)/sum(tb)), "%)", collapse = "; ")
}

#--- 3) Tạo bảng thống kê tổng hợp theo kiểu biến ----
build_dict <- function(df, meta_tbl = NULL) {
  stopifnot(is.data.frame(df))
  cols <- names(df)

  # Stats cho từng cột
  stats_tbl <- purrr::map_dfr(cols, function(v) {
    x <- df[[v]]
    type <- paste(class(x), collapse = "/")
    miss <- pNA(x)
    uniq <- nuniq(x)

    if (is.numeric(x)) {
      qs  <- pcts(x)
      tibble::tibble(
        var = v, type = "numeric", missing_pct = round(miss,2), n_unique = uniq,
        min = suppressWarnings(min(x, na.rm=TRUE)),
        p01 = unname(qs[1]), p05 = unname(qs[2]), p25 = unname(qs[3]),
        median = unname(qs[4]), mean = suppressWarnings(mean(x, na.rm=TRUE)),
        p75 = unname(qs[5]), p95 = unname(qs[6]), p99 = unname(qs[7]),
        max = suppressWarnings(max(x, na.rm=TRUE)),
        sd  = suppressWarnings(stats::sd(x, na.rm=TRUE)),
        example = NA_character_, top_levels = NA_character_
      )
    } else if (is.logical(x) || (is.integer(x) && uniq<=3)) {
      lv <- top_levels(as.character(x), k = 3)
      tibble::tibble(
        var = v, type = "binary/ordinal", missing_pct = round(miss,2), n_unique = uniq,
        min = NA, p01=NA, p05=NA, p25=NA, median=NA, mean=NA, p75=NA, p95=NA, p99=NA, max=NA, sd=NA,
        example = paste(head(unique(x), 3), collapse=", "),
        top_levels = lv
      )
    } else {
      lv <- top_levels(x, k = 5)
      tibble::tibble(
        var = v, type = "categorical/text", missing_pct = round(miss,2), n_unique = uniq,
        min = NA, p01=NA, p05=NA, p25=NA, median=NA, mean=NA, p75=NA, p95=NA, p99=NA, max=NA, sd=NA,
        example = paste(head(unique(x), 3), collapse=", "),
        top_levels = lv
      )
    }
  })

  # Gộp với meta (ý nghĩa kinh tế lượng)
  out <- stats_tbl %>%
    dplyr::left_join(meta_tbl, by = "var") %>%
    dplyr::relocate(var, label, type, unit, missing_pct, n_unique)

  out
}

Đoạn code này được thiết kế để xây dựng từ điển dữ liệu (data dictionary) — một bước quan trọng trong quy trình tiền xử lý và phân tích thống kê. Về kỹ thuật, code gồm ba phần chính: * (1) Bảng mô tả var_meta liệt kê thông tin chi tiết cho từng biến, bao gồm tên, nhãn, đơn vị đo, định nghĩa, cách tính toán, vai trò kinh tế và kỳ vọng dấu trong mô hình vỡ nợ. * (2) Các hàm phụ trợ như pNA, nuniq, pcts, top_levels được dùng để tính tỷ lệ thiếu dữ liệu, số giá trị duy nhất, các phân vị thống kê và tóm tắt nhóm giá trị phổ biến của biến phân loại. * (3) Hàm tổng hợp build_dict() thực hiện quét toàn bộ các cột của bộ dữ liệu, xác định loại biến (số, nhị phân, phân loại), tính các chỉ tiêu thống kê mô tả (min, max, trung vị, độ lệch chuẩn, v.v.), và ghép kết quả với bảng var_meta để tạo ra bảng từ điển hoàn chỉnh. Về ý nghĩa thống kê, phần này giúp hiểu rõ bản chất và chất lượng của dữ liệu, phát hiện các vấn đề như thiếu giá trị, phân phối lệch, hay biến nhiễu, đồng thời cung cấp nền tảng cho việc lựa chọn và biến đổi biến phù hợp trong mô hình kinh tế lượng dự báo rủi ro tín dụng. ### Tạo bảng tóm tắt của bộ dữ liệu

Đoạn mã tạo một bảng tóm tắt cấu trúc và loại dữ liệu của các biến trong bộ dữ liệu loans_raw và sắp xếp bảng này theo loại biến (numeric, binary/ordinal, categorical/text) để dễ xem xét.

dict_tbl <- build_dict(loans_raw, var_meta)

dict_tbl <- dict_tbl %>%
  dplyr::arrange(factor(type, levels=c("numeric","binary/ordinal","categorical/text")),
                 dplyr::desc(!is.na(mean)), var)
head(dict_tbl)

Về kỹ thuật, lệnh build_dict(loans_raw, var_meta) gọi hàm đã định nghĩa để kết hợp dữ liệu gốc (loans_raw) với bảng mô tả biến (var_meta), tạo ra một bảng tổng hợp (dict_tbl) chứa cả thông tin thống kê (phân vị, trung bình, độ lệch chuẩn, v.v.) và ý nghĩa kinh tế lượng của từng biến. Sau đó, dữ liệu được sắp xếp lại theo loại biến — gồm ba nhóm chính:numeric (biến định lượng), binary/ordinal (biến nhị phân hoặc thứ bậc), categorical/text (biến định tính hoặc văn bản). Cuối cùng, lệnh head(dict_tbl) giúp xem trước vài dòng đầu của bảng này để kiểm tra tính hợp lệ và cấu trúc hiển thị. Về ý nghĩa thống kê, bước này giúp đảm bảo tính nhất quán và chất lượng của dữ liệu trước khi mô hình hóa, đồng thời cho phép người phân tích hiểu rõ đặc điểm của từng biến — từ loại dữ liệu, mức độ thiếu, đến vai trò tiềm năng trong mô hình dự báo rủi ro tín dụng.

Chuẩn bị và định dạng bảng tóm tắt dữ liệu để xuất báo cáo (PDF)

stopifnot(exists("dict_tbl"))

dict_clean <- dict_tbl %>%
dplyr::mutate(
unit        = dplyr::if_else(is.na(.data$unit), "", .data$unit),
missing_pct = round(missing_pct, 2),
n_unique    = as.integer(round(n_unique)),
mean        = round(mean,   4),
sd          = round(sd,     4),
p25         = round(p25,    4),
median      = round(median, 4),
p75         = round(p75,    4)
)

tbl1 <- dict_clean %>%
dplyr::select(tidyselect::any_of(c("var","label","type","unit"))) %>%
kableExtra::kbl(
booktabs  = knitr::is_latex_output(),
longtable = FALSE,
linesep   = ""
)

if (knitr::is_latex_output()) {
tbl1 <- kableExtra::kable_styling(
tbl1,
latex_options = c("HOLD_position","striped"),
full_width = FALSE,
font_size = 9
)
}

tbl1

var	label	type	unit
default	Default flag	numeric	0/1
dti_n	Debt-to-income (DTI)	numeric	%
experience_c	Has employment info	numeric	0/1
fico_n	FICO score	numeric	score
id	Loan ID	numeric
loan_amnt	Loan amount	numeric	USD
revenue	Annual income	numeric	USD/year
addr_state	Borrower state	categorical/text	category
desc	User loan description	categorical/text	text
emp_length	Employment length	categorical/text	years
home_ownership_n	Home ownership	categorical/text	category
issue_d	Issue date	categorical/text
purpose	Loan purpose	categorical/text	category
title	User loan title	categorical/text	text
zip_code	Zip (3 digits)	categorical/text	category

Đoạn mã này đảm bảo tính toàn vẹn và dễ đọc của Từ điển Dữ liệu trước khi trình bày. Về mặt kỹ thuật, lệnh stopifnot(exists(“dict_tbl”)) được sử dụng để kiểm tra điều kiện tiên quyết, đảm bảo biến chứa từ điển (dict_tbl) phải tồn tại. Tiếp theo, khối lệnh dplyr::mutate thực hiện chuẩn hóa dữ liệu bằng cách chuyển đổi các cột số liệu thống kê mô tả (như missing_pct, mean, sd, p25, median, p75) về số thập phân cố định (làm tròn 2 hoặc 4 chữ số), và xử lý các giá trị NA trong cột unit thành chuỗi rỗng (““). Mục đích thống kê của việc làm tròn là để tăng cường khả năng đọc của các số liệu thống kê mô tả trong báo cáo. Cuối cùng, khối lệnh kableExtra::kbl chỉ chọn các cột cơ bản (var, label, type, unit) và sử dụng các tùy chọn kable_styling để định dạng bảng (ví dụ: font_size = 9, striped) một cách chuyên nghiệp, đảm bảo bảng từ điển được trình bày gọn gàng và thẩm mỹ, đặc biệt trong đầu ra LaTeX. ### Thống kê mô tả

stopifnot(exists("dict_clean"))

tbl2 <- dict_clean %>%
dplyr::select(tidyselect::any_of(c(
"var","missing_pct","n_unique",
"mean","sd","p25","median","p75"
))) %>%
kableExtra::kbl(
booktabs  = knitr::is_latex_output(),
longtable = TRUE,
linesep   = ""
)

if (knitr::is_latex_output()) {
tbl2 <- kableExtra::kable_styling(
tbl2,
latex_options = c("HOLD_position","striped","scale_down"),
full_width = FALSE,
font_size = 8
)
}

tbl2

var	missing_pct	n_unique	mean	sd	p25	median	p75
default	0.00	2	1.995000e-01	3.996000e-01	0.0	0.00	0.00
dti_n	0.00	6625	1.845430e+01	1.107340e+01	11.9	17.76	24.28
experience_c	0.00	2	1.000000e+00	4.900000e-03	1.0	1.00	1.00
fico_n	0.00	48	6.978592e+02	3.189370e+01	672.0	692.00	712.00
id	0.00	1048575	5.758232e+07	3.512422e+07	30555560.5	57315190.00	75243770.50
loan_amnt	0.00	1545	1.446931e+04	8.730590e+03	8000.0	12000.00	20000.00
revenue	0.00	55236	7.759825e+04	7.099969e+04	46696.0	65000.00	93000.00
addr_state	0.00	51	NA	NA	NA	NA	NA
desc	95.13	47636	NA	NA	NA	NA	NA
emp_length	0.00	12	NA	NA	NA	NA	NA
home_ownership_n	0.00	4	NA	NA	NA	NA	NA
issue_d	0.00	109	NA	NA	NA	NA	NA
purpose	0.00	14	NA	NA	NA	NA	NA
title	1.27	31030	NA	NA	NA	NA	NA
zip_code	0.00	937	NA	NA	NA	NA	NA

Bảng Thống kê Mô tả cho thấy chất lượng dữ liệu tốt với hầu hết các biến rủi ro cốt lõi có tỷ lệ thiếu (missing_pct) là 0.00. Biến mục tiêu default có tỷ lệ vỡ nợ trung bình (mean) là 19.95% (cho thấy mất cân bằng nhẹ). Các biến định lượng chính như revenue có độ lệch lớn (SD 70,999 vs. Median 65,000), gợi ý về phân phối lệch và sự hiện diện của giá trị ngoại lệ cần được xử lý. Chỉ có hai biến (desc và title) có tỷ lệ thiếu đáng kể, trong đó desc có thể bị loại bỏ do tỷ lệ thiếu quá cao (95.13%).

Trình bày Bảng tóm tắt dữ liệu

Chunk này sử dụng thư viện kableExtra để tạo ra và định dạng bảng từ điển dữ liệu đã được chuẩn bị (pdf_dict) thành một bảng chất lượng cao (thường là LaTeX) để đưa vào báo cáo chính thức.

loans0 <- loans_raw
if (!("Default" %in% names(loans0)) && ("default" %in% names(loans0))) {
  loans0 <- loans0 %>%
    mutate(
      Default = readr::parse_number(
        as.character(.data[["default"]])
      )
    )
}

num_candidates <- c(
  "loan_amnt",  
  "revenue",     
  "dti_n",       
  "fico_n")

Khối code này thực hiện các bước tiền xử lý dữ liệu ban đầu. Kỹ thuật cốt lõi là tạo bản sao dữ liệu gốc (loans0 <- loans_raw) và sử dụng lệnh if để chuẩn hóa biến mục tiêu vỡ nợ (default). Nếu biến mục tiêu chưa ở dạng chuẩn hóa (Default không tồn tại), mã sẽ chuyển đổi biến gốc thành dạng số nhị phân (0/1) bằng readr::parse_number. Ý nghĩa thống kê của bước này là đảm bảo biến mục tiêu đã sẵn sàng cho các mô hình hồi quy phân loại (như Hồi quy Logistic). Cuối cùng, mã định nghĩa vector num_candidates chứa danh sách các biến định lượng cốt lõi (loan_amnt, revenue, dti_n, fico_n) để sử dụng trong phân tích mô hình sau này.

Xử lý dữ liệu thô và mã hóa dữ liệu

Chuyển đổi loại biến và tạo biến sơ bộ

Chunk này tập trung vào việc chuẩn hóa loại dữ liệu cho các biến số, đảm bảo chúng phù hợp với yêu cầu của các thuật toán thống kê và học máy tiếp theo. Đồng thời, một số biến mới cũng được tạo ra hoặc chuyển đổi sang định dạng số hóa dễ sử dụng hơn.

loans1 <- loans0 %>%
mutate(
issue_d   = coalesce_issue_date(.),
across(all_of(num_candidates), safe_parse_num),
emp_years = if ("emp_length" %in% names(.)) emp_to_years(emp_length) else NA_real_,
Default   = if ("Default" %in% names(.)) as.integer(readr::parse_number(as.character(Default))) else NA_integer_,
purpose   = if ("purpose" %in% names(.)) as.factor(purpose) else factor(NA),
home_ownership_n = if ("home_ownership_n" %in% names(.)) as.factor(home_ownership_n) else factor(NA),
addr_state= if ("addr_state" %in% names(.)) as.factor(addr_state) else factor(NA),
title     = if ("title" %in% names(.)) as.character(title) else NA_character_,
desc      = if ("desc"  %in% names(.))  as.character(desc)  else NA_character_
)

Khối code này thực hiện chuẩn bị dữ liệu bằng cách chuyển đổi kiểu dữ liệu của các biến chính. Kỹ thuật bao gồm việc chuẩn hóa biến ngày (issue_d) và đảm bảo các biến định lượng cốt lõi (num_candidates) được chuyển đổi an toàn sang kiểu số. Ý nghĩa quan trọng là chuyển đổi các biến phân loại (purpose, home_ownership_n, addr_state) thành kiểu Factor, đây là định dạng bắt buộc cho mô hình thống kê phân loại. Đồng thời, code thực hiện Feature Engineering ban đầu bằng cách chuyển đổi độ dài kinh nghiệm (emp_length) thành số năm (emp_years). Tóm lại, khối này hoàn thành việc xác định kiểu dữ liệu để sẵn sàng cho bước mô hình hóa tiếp theo.

Tạo biến nâng cao, xử lý ngoại lai và điền khuyết cơ bản

Chunk này tiến hành các kỹ thuật xử lý dữ liệu chuyên sâu hơn như tạo ra các chỉ số tài chính quan trọng, xử lý các giá trị ngoại lai, và điền khuyết (imputation) cho các biến phân loại và ký tự để chuẩn bị dữ liệu cho mô hình hóa.

winsor_probs <- if (exists("params") && !is.null(params$winsor_probs)) params$winsor_probs else c(.01,.99)

loans2 <- loans1 %>%
mutate(
lti       = dplyr::if_else(is.finite(revenue) & revenue > 0, loan_amnt/revenue, NA_real_),
revenue_w = if ("revenue" %in% names(.)) winsorize_vec(revenue, c(0.01, 0.99)) else NA_real_,
) %>%
mutate(across(where(is.character), ~na_if(., ""))) %>%
mutate(across(where(is.factor), ~forcats::fct_explicit_na(., na_level = "missing"))) %>%
mutate(across(where(is.character), mode_impute))

Khối code này thực hiện Feature Engineering và Imputation cơ bản. Kỹ thuật cốt lõi là xác định ngưỡng Winsorization (mặc định 1% và 99%) và tạo ra hai biến mới: Tỷ lệ Khoản vay trên Thu nhập (lti) và phiên bản xử lý ngoại lệ (revenue_w) của biến thu nhập bằng phương pháp Winsorization. Ý nghĩa thống kê của bước này là tạo chỉ số rủi ro mới (lti) và ổn định dữ liệu bằng cách giảm thiểu ảnh hưởng của các giá trị ngoại lệ. Cuối cùng, code xử lý giá trị thiếu bằng cách chuyển chuỗi rỗng thành NA, sau đó điền khuyết rõ ràng các biến Factor thành “missing” và điền khuyết các biến Character bằng Mode Imputation.

Điền khuyết nâng cao cho biến số định lượng

Chunk này thực hiện kỹ thuật điền khuyết đa biến (MICE) cho các biến số định lượng cốt lõi. Đây là một phương pháp điền khuyết hiện đại nhằm duy trì mối quan hệ thống kê giữa các biến số sau khi điền khuyết giá trị thiếu.

use_mice_flag <- isTRUE(if (exists("params")) params$use_mice else FALSE)
if (use_mice_flag) {
log_info("Chạy MICE cho cụm biến numeric chính...")
if (!requireNamespace("mice", quietly = TRUE)) stop("Thiếu gói 'mice' 
                                                    (tắt params$use_mice hoặc cài đặt).")
library(mice)
mice_df <- loans2 %>% dplyr::select(all_of(c("loan_amnt","revenue","revenue_w",
                                             "dti_n","fico_n","emp_years","lti")))
imp <- mice::mice(mice_df, m=1, maxit=5, method='pmm', seed=3698, printFlag=FALSE)
mice_completed <- mice::complete(imp, 1)
loans2 <- loans2 %>% mutate(across(names(mice_completed), ~ mice_completed[[cur_column()]]))
}

Khối code này thực hiện Điền khuyết đa biến (MICE), một kỹ thuật xử lý giá trị thiếu nâng cao, nếu cờ use_mice_flag được bật. Kỹ thuật cốt lõi là sử dụng thuật toán pmm (Predictive Mean Matching) để ước tính và điền khuyết các giá trị thiếu trong các biến định lượng chính (loan_amnt, revenue, fico_n, v.v.). Ý nghĩa Thống kê: Việc này giúp tạo ra một tập dữ liệu đầy đủ và chính xác hơn, giảm thiểu sai lệch ước lượng (bias) so với các phương pháp điền khuyết đơn giản, đảm bảo dữ liệu sẵn sàng cho mô hình hóa. Mã sử dụng seed=3698 để đảm bảo tính tái lập của quá trình điền khuyết.

Lọc dữ liệu thiếu/không hợp lệ và phân nhóm biến số

Chunk này tập trung vào việc đảm bảo tính toàn vẹn của dữ liệu bằng cách loại bỏ các quan sát không hợp lệ đối với các biến tài chính cốt lõi, đồng thời tạo ra biến phân loại mới từ biến định lượng để chuẩn bị cho các phân tích thống kê và mô hình hóa.

df_cleaned <- loans2 %>%
filter(!is.na(loan_amnt), loan_amnt > 0, !is.na(revenue), revenue > 0) %>%
mutate(
fico_group = if ("fico_n" %in% names(.)) cut(
fico_n, breaks=c(0,580,670,740,800,Inf),
labels=c("Poor","Fair","Good","Very Good","Exceptional"),
right=FALSE, include.lowest=TRUE
) else factor(NA)
)

Khối code này thực hiện bước lọc dữ liệu cuối cùng và phân nhóm biến (Variable Binning). Về mặt Kỹ thuật, lệnh filter loại bỏ các hồ sơ có giá trị thiếu (NA) hoặc không hợp lệ (<= 0) ở hai biến số tiền vay (loan_amnt) và thu nhập (revenue), tạo ra tập dữ liệu sạch cuối cùng (df_cleaned). Ý nghĩa Thống kê quan trọng là việc tạo ra biến phân loại fico_group từ biến điểm FICO định lượng (fico_n) bằng hàm cut(). Biến FICO được chia thành 5 nhóm rủi ro tín dụng tiêu chuẩn (“Poor” đến “Exceptional”). Việc này giúp đơn giản hóa mô hình và tăng cường khả năng giải thích kết quả bằng cách chuyển một biến liên tục thành các nhóm rủi ro dễ hiểu cho mục đích báo cáo.

Mã hóa biến định lượng bằng kỹ thuật phân nhóm

Chunk này thực hiện việc phân nhóm (binning) nhiều biến số tài chính liên tục quan trọng thành các biến phân loại theo thứ tự (Ordinal Factor). Việc này nhằm mục đích đơn giản hóa mối quan hệ, làm cho dữ liệu trở nên dễ giải thích hơn và chuẩn bị cho các phân tích tỷ lệ vỡ nợ theo nhóm.

df_binned <- df_cleaned %>%
mutate(
dti_band = if ("dti_n" %in% names(.)) cut(
dti_n, breaks=c(-Inf,5,10,15,20,30,40,Inf),
labels=c("≤5","5–10","10–15","15–20","20–30","30–40",">40")
) else factor(NA),
emp_band = cut(emp_years, breaks=c(-Inf,1,3,5,10,Inf),
labels=c("<1y","1–3y","3–5y","5–10y","10y+")),
rev_band = cut(revenue_w,
breaks=c(0,40e3,60e3,80e3,120e3,200e3,Inf),
labels=c("<40k","40–60k","60–80k","80–120k","120–200k",">200k"),
right=FALSE),
lti_band = cut(lti, breaks=c(-Inf,.05,.10,.20,.30,Inf),
labels=c("≤5%","5–10%","10–20%","20–30%",">30%"))
)

Khối code này tạo ra tập dữ liệu cuối cùng (df_binned) bằng cách áp dụng phân nhóm (binning) cho bốn biến định lượng rủi ro cốt lõi. Kỹ thuật được sử dụng là hàm cut() để chia các biến liên tục thành các cấp độ phân loại có ý nghĩa: dti_n thành dti_band (Tỷ lệ Nợ trên Thu nhập), emp_years thành emp_band (Thời gian làm việc), revenue_w thành rev_band (Thu nhập hàng năm đã được Winsorization), và lti thành lti_band (Tỷ lệ Khoản vay trên Thu nhập). Ý nghĩa thống kê của việc phân nhóm là chuyển các biến rủi ro liên tục thành các nhóm rủi ro rời rạc dễ hiểu hơn, giúp tăng khả năng giải thích của các mô hình và cho phép tính toán các số liệu rủi ro (như Default Rate, Lift) cho từng nhóm một cách trực quan trong báo cáo. ## Thực hiện các thống kê cơ bản

Bảng tóm tắt chất lượng dữ liệu ban đầu

library(dplyr)
library(scales)
library(kableExtra)

is_tex <- knitr::is_latex_output()
is_html <- knitr::is_html_output()

n_obs  <- nrow(loans_raw)
n_vars <- ncol(loans_raw)
dup_id <- if ("id" %in% names(loans_raw)) sum(duplicated(loans_raw$id)) else NA
miss_n <- sum(!complete.cases(loans_raw))
miss_pct <- miss_n / n_obs

tbl <- tibble(
  Observations = comma(n_obs),
  Variables = comma(n_vars),
  `Duplicate IDs` = comma(dup_id),
  `Rows with missing` = paste0(comma(miss_n), " (", percent(miss_pct, 0.1), ")")
)

kbl(tbl, booktabs = is_tex, caption = "Dataset Quality Snapshot") %>%
  kable_styling(full_width = FALSE,
                latex_options = if (is_tex) "hold_position" else NULL,
                bootstrap_options = if (is_html) "striped" else NULL)

Dataset Quality Snapshot

Observations	Variables	Duplicate IDs	Rows with missing
1,048,575	15	0	997,515 (95.1%)

Khối code này thực hiện việc tính toán và trình bày Dataset Quality Snapshot (ảnh chụp chất lượng dữ liệu). Kỹ thuật là sử dụng các hàm cơ bản (nrow, sum(duplicated)) để tính toán số lượng quan sát, biến số, và số ID trùng lặp. Ý nghĩa Thống kê cốt lõi là xác định tỷ lệ hàng có giá trị thiếu (miss_pct) trên dữ liệu thô (loans_raw). Kết quả cho thấy tập dữ liệu có quy mô lớn (1,048,575 Quan sát) và không có hồ sơ trùng lặp (0 Duplicate IDs). Tuy nhiên, có tới 95.1% số hàng chứa giá trị thiếu. Tỷ lệ thiếu cao này là một vấn đề nghiêm trọng, khẳng định bước điền khuyết (Imputation) là bước quan trọng và phức tạp nhất trong toàn bộ quy trình tiền xử lý dữ liệu.

Phân tích thành phần chính và lựa chọn chiều dữ liệu tối ưu

Chunk này thực hiện phân tích thành phần chính trên các biến số định lượng cốt lõi để giảm chiều dữ liệu. Mục đích là tìm ra một tập hợp các thành phần chính (Principal Components - PC) có thể giải thích phần lớn phương sai của dữ liệu gốc, giúp đơn giản hóa việc mô hình hóa.

library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(kableExtra)

pca_df <- df_cleaned %>%
  dplyr::select(loan_amnt, revenue_w, dti_n, fico_n, emp_years, lti) %>%
  dplyr::mutate(dplyr::across(everything(), as.numeric)) %>%
  tidyr::drop_na()

keep_pca <- names(pca_df)[sapply(pca_df, function(x){
  s <- stats::sd(x, na.rm = TRUE); is.finite(s) && s > 0
})]
pca_df <- dplyr::select(pca_df, tidyselect::any_of(keep_pca))

pca_res <- NULL
if (nrow(pca_df) > 5 && ncol(pca_df) >= 2) {
  pca_res <- stats::prcomp(pca_df, center = TRUE, scale. = TRUE)
  imp <- summary(pca_res)$importance
  var_expl <- as.numeric(round(100 * imp["Proportion of Variance", ], 1))
  cum_expl <- as.numeric(round(100 * imp["Cumulative Proportion", ], 1))
  eig_vals <- as.numeric(round(pca_res$sdev^2, 3))  # eigenvalues

  pcs <- paste0("PC", seq_along(var_expl))
  df_imp <- tibble::tibble(
    PC = pcs,
    `Eigenvalue` = eig_vals,
    `Variance (%)` = var_expl,
    `Cumulative (%)` = cum_expl
  )
  target <- 80
  k_sel <- dplyr::coalesce(which(cum_expl >= target)[1], length(cum_expl))
  df_imp <- df_imp %>%
    mutate(Selected = if_else(row_number() <= k_sel, "✓", ""))

  kbl <- kbl(
    df_imp,
    booktabs = TRUE,
    caption = paste0("PCA-Tỷ lệ phương sai giải thích (k = ", k_sel,
                     ", Cumulative = ", cum_expl[k_sel], "%)")
  ) %>%
    kable_styling(latex_options = c("striped","hold_position"), font_size = 10)

  if (k_sel >= 1) {
    kbl <- kbl %>% row_spec(1:k_sel, bold = TRUE, color = "white", background = "#6E79C0")
  }
  print(kbl)

  assign("k_sel", k_sel, envir = .GlobalEnv)
  assign("cum_k", cum_expl[k_sel], envir = .GlobalEnv)
}

PCA-Tỷ lệ phương sai giải thích (k = 4, Cumulative = 85%)

PC	Eigenvalue	Variance (%)	Cumulative (%)	Selected
PC1	1.647	27.4	27.4	✓
PC2	1.503	25.0	52.5	✓
PC3	1.007	16.8	69.3	✓
PC4	0.941	15.7	85.0	✓
PC5	0.826	13.8	98.7
PC6	0.076	1.3	100.0

Khối code này thực hiện Phân tích Thành phần chính (PCA) trên sáu biến rủi ro định lượng cốt lõi (loan_amnt, fico_n, v.v.) nhằm giảm chiều dữ liệu. Về mặt kỹ thuật và thống kê, code sử dụng hàm stats::prcomp với tùy chọn scale. = TRUE (chuẩn hóa dữ liệu) để đảm bảo các biến đóng góp công bằng, tránh sai lệch. Mã tính toán Eigenvalue và phương sai tích lũy. Kết quả cho thấy cần chọn 4 Thành phần Chính (PC1 đến PC4) để giải thích 85.0% tổng phương sai. Việc này đạt được mục tiêu giảm chiều dữ liệu thành \(k=4\) biến mới trong khi vẫn giữ lại phần lớn thông tin (85%), tạo ra các biến mới ít tương quan hơn và giúp mô hình hóa hiệu quả hơn.

Bảng phân tích rủi ro tín dụng: Tỷ lệ vỡ nợ theo nhóm LTI

Chunk này tạo ra bảng tổng hợp các chỉ số rủi ro quan trọng cho từng phân khúc Tỷ lệ Vay/Thu nhập (LTI). Phân tích này là cốt lõi trong việc xây dựng các bảng điểm (scorecard) và chính sách quản lý rủi ro tín dụng.

library(dplyr)
library(kableExtra)
library(scales)
library(forcats)

tab_lti <- df_binned %>%
  mutate(lti_band = fct_explicit_na(lti_band, "Missing")) %>%
  group_by(lti_band) %>%
  summarise(n = n(),
            defaults = sum(Default == 1, na.rm = TRUE),
            .groups = "drop") %>%
  mutate(
    exposure = n / sum(n),
    rate     = defaults / n,
    se       = sqrt(rate * (1 - rate) / n),
    ci_lo    = pmax(rate - 1.96 * se, 0),
    ci_hi    = pmin(rate + 1.96 * se, 1)
  )

overall <- with(tab_lti, sum(defaults) / sum(n))
tab_lti <- tab_lti %>%
  mutate(lift = rate / overall) %>%
  bind_rows(summarise(., lti_band="Total", n=sum(n), defaults=sum(defaults),
                      exposure=1, rate=defaults/n, se=sqrt(rate*(1-rate)/n),
                      ci_lo=pmax(rate-1.96*se,0), ci_hi=pmin(rate+1.96*se,1), lift=1))

fmt <- function(x,d=2) formatC(x,format="f",digits=d,big.mark=",")
tab_fmt <- tab_lti %>%
  transmute(
    `LTI band` = lti_band,
    Obs. = formatC(n,format="d",big.mark=","),
    Exposure = percent(exposure,0.1),
    Defaults = formatC(defaults,format="d",big.mark=","),
    `Default rate` = percent(rate,0.01),
    `95% CI` = paste0(percent(ci_lo,0.01),"–",percent(ci_hi,0.01)),
    `Lift vs Overall` = fmt(lift)
  )

kbl(tab_fmt, caption="Bảng 4. Tỷ lệ vỡ nợ theo dải LTI (95% CI & Lift).", booktabs=TRUE,
    align=c("l","r","r","r","r","r","r")) %>%
  kable_styling(latex_options=c("striped","hold_position","scale_down"), font_size=9) %>%
  add_header_above(c(" " = 1, "Volume" = 3, "Risk Metrics" = 3)) %>%
  row_spec(nrow(tab_fmt), bold=TRUE) %>%
  footnote(general="Default rate = Defaults / Obs.; CI: p ± 1.96·SE; Lift = rate nhóm / toàn tập.",
           threeparttable=TRUE, general_title="")

Bảng 4. Tỷ lệ vỡ nợ theo dải LTI (95% CI & Lift).

	Volume			Risk Metrics
LTI band	Obs.	Exposure	Defaults	Default rate	95% CI	Lift vs Overall
≤5%	48,944	4.7%	5,712	11.67%	11.39%–11.95%	0.59
5–10%	139,915	13.3%	19,196	13.72%	13.54%–13.90%	0.69
10–20%	360,820	34.4%	59,701	16.55%	16.42%–16.67%	0.83
20–30%	271,946	25.9%	59,559	21.90%	21.75%–22.06%	1.10
>30%	226,950	21.6%	64,992	28.64%	28.45%–28.82%	1.44
Total	1,048,575	100.0%	209,160	19.95%	19.87%–20.02%	1.00
Default rate = Defaults / Obs.; CI: p ± 1.96·SE; Lift = rate nhóm / toàn tập.

Khối code này thực hiện quy trình tính toán rủi ro chuyên sâu trên biến LTI (Loan-to-Income). Về mặt kỹ thuật và thống kê, mã tính toán tỷ lệ vỡ nợ, khoảng tin cậy 95% (CI) và chỉ số Lift, sử dụng phương pháp thống kê chuẩn cho tỷ lệ để xác định rủi ro. Kết quả xác nhận mối quan hệ tỷ lệ thuận mạnh mẽ giữa LTI và rủi ro: nhóm LTI \(\le 5\%\) là nhóm an toàn nhất (Lift 0.59, Tỷ lệ Vỡ nợ 11.67%), trong khi nhóm LTI \(>30\%\) là nhóm rủi ro cao nhất (Lift 1.44, Tỷ lệ Vỡ nợ 28.64%). Sự khác biệt này là có ý nghĩa thống kê và khẳng định LTI band là một yếu tố phân loại rủi ro cực kỳ hiệu quả cần được đưa vào các chiến lược tín dụng.

Bảng phân tích rủi ro tín dụng: Tỷ lệ vỡ nợ theo tình trạng sở hữu nhà

Chunk này tạo ra bảng tổng hợp các chỉ số rủi ro quan trọng cho từng phân khúc Tình trạng Sở hữu nhà. Đây là phân tích rủi ro then chốt, vì sở hữu nhà được coi là thước đo mức độ ổn định tài chính và tài sản thế chấp tiềm năng của người vay.

library(dplyr)
library(kableExtra)
library(scales)
library(forcats)

tab_home <- df_binned %>%
  mutate(home_group = fct_explicit_na(home_ownership_n, "Missing")) %>%
  group_by(home_group) %>%
  summarise(n = n(),
            defaults = sum(Default == 1, na.rm = TRUE),
            .groups = "drop") %>%
  mutate(rate = defaults / n,
         se = sqrt(rate * (1 - rate) / n),
         ci_lo = pmax(rate - 1.96 * se, 0),
         ci_hi = pmin(rate + 1.96 * se, 1),
         exposure = n / sum(n))

overall <- with(tab_home, sum(defaults) / sum(n))
tab_home <- tab_home %>%
  mutate(lift = rate / overall) %>%
  bind_rows(summarise(., home_group = "Total", n = sum(n),
                      defaults = sum(defaults), exposure = 1,
                      rate = defaults / n, se = sqrt(rate * (1 - rate) / n),
                      ci_lo = pmax(rate - 1.96 * se, 0),
                      ci_hi = pmin(rate + 1.96 * se, 1), lift = 1))

tab_fmt <- tab_home %>%
  transmute(
    `Home Ownership` = home_group,
    Obs. = formatC(n, format = "d", big.mark = ","),
    Exposure = percent(exposure, 0.1),
    Defaults = formatC(defaults, format = "d", big.mark = ","),
    `Default rate` = percent(rate, 0.01),
    `95% CI` = paste0(percent(ci_lo, 0.01), "–", percent(ci_hi, 0.01)),
    `Lift vs Overall` = formatC(lift, format = "f", digits = 2)
  )

kbl(tab_fmt, caption = "Bảng 5. Tỷ lệ vỡ nợ theo tình trạng sở hữu nhà (95% CI & Lift).",
    booktabs = TRUE, align = c("l", "r", "r", "r", "r", "r", "r")) %>%
  kable_styling(latex_options = c("striped", "hold_position", "scale_down"),
                font_size = 9) %>%
  add_header_above(c(" " = 1, "Volume" = 3, "Risk Metrics" = 3)) %>%
  row_spec(nrow(tab_fmt), bold = TRUE) %>%
  footnote(general = "Default rate = Defaults / Obs.; CI = p ± 1.96·SE; Lift = rate nhóm / toàn tập.",
           threeparttable = TRUE, general_title = "")

Bảng 5. Tỷ lệ vỡ nợ theo tình trạng sở hữu nhà (95% CI & Lift).

	Volume			Risk Metrics
Home Ownership	Obs.	Exposure	Defaults	Default rate	95% CI	Lift vs Overall
MORTGAGE	517,417	49.3%	88,860	17.17%	17.07%–17.28%	0.86
OTHER	366	0.0%	70	19.13%	15.10%–23.15%	0.96
OWN	114,563	10.9%	23,418	20.44%	20.21%–20.67%	1.02
RENT	416,229	39.7%	96,812	23.26%	23.13%–23.39%	1.17
Total	1,048,575	100.0%	209,160	19.95%	19.87%–20.02%	1.00
Default rate = Defaults / Obs.; CI = p ± 1.96·SE; Lift = rate nhóm / toàn tập.

Khối code này thực hiện quy trình tính toán rủi ro chi tiết trên biến Tình trạng Sở hữu nhà (home_ownership_n). Về mặt kỹ thuật và thống kê, mã tính toán Tỷ lệ vỡ nợ (rate), khoảng tin cậy 95% (CI), và chỉ số Lift. Sau đó, mã định dạng chuyên nghiệp kết quả bằng kableExtra. Kết quả cho thấy tình trạng sở hữu nhà là yếu tố phân loại rủi ro rất hiệu quả so với tỷ lệ vỡ nợ tổng thể là 19.95%. Nhóm RENT (Thuê nhà) là phân khúc rủi ro cao nhất (23.26%, Lift 1.17). Ngược lại, nhóm MORTGAGE (Vay thế chấp) là nhóm an toàn nhất (17.17%, Lift 0.86). Sự khác biệt rõ ràng và không chồng lấn của 95% CI giữa hai nhóm này khẳng định sự khác biệt về rủi ro là có ý nghĩa thống kê cao, hỗ trợ cho việc đưa ra quyết định tín dụng phân biệt. ## Trực quan hóa bộ dữ liệu

Biểu đồ phân bổ số tiền vay

Chunk này tạo ra Biểu đồ Histogram và Mật độ (Density Plot) chuyên nghiệp để khám phá phân bố của biến số định lượng cốt lõi loan_amnt (Số tiền vay), đồng thời trực quan hóa các thống kê trung tâm (Mean, Median) và độ lệch.

library(dplyr)
library(ggplot2)
library(scales)

if (need_cols(df_cleaned, "loan_amnt")) {
  x <- df_cleaned$loan_amnt |> as.numeric() |> na.omit()
  if (length(x) > 1) {
    n <- length(x)
    q <- quantile(x, c(.25,.5,.75,.95), na.rm=TRUE)
    mu <- mean(x)
    binw <- max(1000, round(2*IQR(x)/(n^(1/3)), -2))
    tail95 <- mean(x >= q[["95%"]])

    df <- data.frame(loan_amnt = x)
    ggplot(df, aes(x = loan_amnt)) +
      annotate("rect", xmin=q[1], xmax=q[3], ymin=0, ymax=Inf, alpha=.15, fill="#DDF5EA") +
      geom_histogram(binwidth=binw, fill="#A5C9A1", color="#6E8F6C", linewidth=.25) +
      geom_density(aes(y = ..scaled.. * max(..count..)), color="#146C43", alpha=.25) +
      geom_vline(xintercept=c(q[1],q[2],q[3],mu),
                 linetype=c(3,2,3,1), color=c("#8E8E8E","#0B3B2E","#8E8E8E","#1E8449"), linewidth=c(.4,.8,.4,.7)) +
      labs(
        title="Phân bổ Số tiền vay",
        subtitle=paste0("n=",format(n,big.mark=",")," | Bin≈",dollar(binw)," | Tail≥P95: ",percent(tail95,0.1)),
        x="Số tiền vay (USD)", y="Số khách hàng",
        caption="Vùng IQR; Đứt = Median; Liền = Mean; đường trơn = mật độ"
      ) +
      scale_x_continuous(labels=dollar_format(), breaks=pretty(x,6)) +
      scale_y_continuous(labels=comma, expand=expansion(mult=c(0,.1))) +
      theme_minimal(base_size=30) +
      theme(plot.title=element_text(face="bold", size=30),
            panel.grid.minor=element_blank(),
            plot.caption=element_text(size=30,color="grey35"))
  }
}

Khối code này sử dụng ggplot2 để tạo một biểu đồ Histogram kèm đường mật độ (Density) cho biến loan_amnt, nhằm kiểm tra phân phối. Về mặt kkỹ thuật và thống kê, mã xác định số lượng quan sát (n), tính toán các tứ phân vị (q), trung bình (mu), và đặc biệt là tính toán chiều rộng cột (binwidth) bằng quy tắc tối ưu hóa Freedman-Diaconis (2*IQR(x)/(n^(1/3))) để đảm bảo chất lượng biểu đồ.
Nó cũng tính toán phần đuôi (Tail \(\ge\) P95) để đo lường tỷ lệ các khoản vay lớn. Kết quả trực quan từ biểu đồ cho thấy phân bố số tiền vay bị lệch phải rõ rệt (Right-skewed). Mặc dù Trung vị (Median, đường đứt nét) nằm quanh $12,000, nhưng Trung bình (Mean, đường liền) lại cao hơn, dịch về phía $14,469, khẳng định sự lệch phải do sự hiện diện của các khoản vay có giá trị lớn. Ý nghĩa thống kê: Phân bố lệch phải cho thấy cần phải chuẩn hóa biến này (ví dụ: logarit hóa) nếu sử dụng trong các mô hình yêu cầu tính chuẩn tắc. Chỉ số Tail \(\ge\) P95 là 5.1% cho biết các khoản vay trên $20,000 chiếm một tỷ lệ nhỏ nhưng có ảnh hưởng đáng kể đến giá trị trung bình.

Biểu đồ Pareto: Phân bổ mục đích vay (Purpose)

Chunk này tạo ra Biểu đồ Pareto để phân tích sự phân bổ tần suất của biến phân loại purpose (Mục đích vay). Biểu đồ này giúp xác định các nhóm mục đích vay chiếm ưu thế, tuân theo Quy tắc 80/20, hỗ trợ cho việc mã hóa biến số và phân tích rủi ro theo nhóm.

library(dplyr); library(ggplot2); library(scales); library(stringr); library(forcats)

top_n <- 12; other_label <- "other"; missing_label <- "missing"
pareto_threshold <- 0.80

if (need_cols(df_cleaned, "purpose")) {
  df_cnt <- df_cleaned %>%
    mutate(purpose = purpose %>%
             as.character() %>% str_squish() %>% str_to_lower() %>%
             fct_explicit_na(missing_label) %>%
             fct_lump_n(n = top_n, other_level = other_label)) %>%
    count(purpose, sort = TRUE, name = "n") %>%
    mutate(share = n / sum(n), cum_share = cumsum(share))

  if (nrow(df_cnt) > 0) {
    cut_idx <- which(df_cnt$cum_share >= pareto_threshold)[1]
    cut_x <- df_cnt$n[cut_idx]; total_n <- sum(df_cnt$n)

    ggplot(df_cnt, aes(y = fct_rev(purpose))) +
      annotate("rect", xmin = -Inf, xmax = cut_x, ymin = -Inf, ymax = Inf,
               fill = alpha("#9AC9C1", .35)) +
      geom_col(aes(x = n), fill = "#77B5AD", color = "#65ABA3", width = 0.7) +
      geom_text(aes(x = n, label = paste0(comma(n), " (", percent(share, 0.1), ")")),
                hjust = -0.1, size = 10) +
      geom_line(aes(x = cum_share * max(n) * 1.2, y = as.numeric(fct_rev(purpose))),
                color = "#FF9F40", linewidth = 1) +
      geom_point(aes(x = cum_share * max(n) * 1.2, y = as.numeric(fct_rev(purpose))),
                 color = "#FF9F40", size = 1.8) +
      geom_vline(xintercept = cut_x, linetype = 3, color = "grey40") +
      scale_x_continuous(
        name = "Số khoản vay", labels = comma, expand = expansion(mult = c(0, 0.1)),
        sec.axis = dup_axis(
          name = "Tích lũy (%)",
          labels = function(z) percent(z / max(df_cnt$n) * 1.2, accuracy = 10)
        )
      ) +
      labs(
        title = paste0("Mục đích vay – Pareto (Top ", top_n, ")"),
        subtitle = paste0("Tổng n = ", format(total_n, big.mark = ","), 
                          " · Ngưỡng Pareto = ", percent(pareto_threshold)),
        caption = paste0("Thanh = số lượng; đường = tích lũy %. Vùng mờ: phần đóng góp đầu đến 80%.")
      ) +
      theme_minimal(base_size = 30) +
      theme(plot.title = element_text(face = "bold", size=30),
            panel.grid.major.y = element_blank(),
            plot.caption = element_text(size = 30, color = "grey35"))
  }
}

Khối code này thực hiện phân tích Pareto (Pareto Analysis) trên biến Mục đích vay (purpose) để xác định các phân khúc chiếm phần lớn khối lượng. Về mặt kỹ thuật và thống kê, mã tiền xử lý và chuẩn hóa các giá trị, sau đó sử dụng hàm fct_lump_n để gom các nhóm nhỏ thành “other”, đồng thời tính toán tỷ lệ tích lũy (cum_share) để tìm ra ngưỡng 80% khối lượng (ngưỡng Pareto). Kết quả trực quan từ biểu đồ Pareto cho thấy sự tập trung rủi ro cao.
Cụ thể, hai mục đích vay chính là debt_consolidation (hợp nhất nợ, 57.8%) và credit_card (thẻ tín dụng, 22.4%) đã đóng góp hơn 80% tổng khối lượng khoản vay. Ý nghĩa thống kê của kết quả này là việc phân tích và mô hình hóa rủi ro cần ưu tiên tập trung vào hai phân khúc thống trị này, vì chúng đại diện cho phần lớn danh mục, từ đó tối ưu hóa hiệu quả mô hình và quản lý rủi ro.

Biểu đồ Pareto: Phân bổ tình trạng sở hữu nhà (Home Ownership)

Chunk này tạo ra Biểu đồ Pareto để khám phá sự phân bổ tần suất của biến home_ownership_n. Phân tích này là cần thiết vì tình trạng sở hữu nhà là một trong những yếu tố rủi ro quan trọng nhất trong việc đánh giá khả năng tài chính của người vay.

library(dplyr)
library(forcats)
library(ggplot2)
library(scales)

pareto_threshold <- 0.80

if (need_cols(df_cleaned, "home_ownership_n")) {
  df_cnt <- df_cleaned %>%
    mutate(home_top = fct_explicit_na(as.factor(home_ownership_n), "Missing")) %>%
    count(home_top, name = "n", sort = TRUE) %>%
    mutate(
      share = n / sum(n),
      cum_share = cumsum(share),
      label = factor(home_top, levels = rev(unique(home_top)))
    )

  cut_idx <- which(df_cnt$cum_share >= pareto_threshold)[1]
  cut_x   <- df_cnt$n[cut_idx]
  total_n <- sum(df_cnt$n)

  ggplot(df_cnt, aes(y = label)) +
    # Vùng mờ 80%
    annotate("rect", xmin = -Inf, xmax = cut_x, ymin = -Inf, ymax = Inf,
             fill = "#E8F6F3", alpha = 0.5) +
    # Cột chính
    geom_col(aes(x = n), fill = "#48C9B0", color = "#17A589", width = 0.7) +
    # Nhãn số lượng
    geom_text(aes(x = n, label = paste0(comma(n), " (", percent(share, 0.1), ")")),
              hjust = -0.1, size = 10, color = "#145A32") +
    # Đường Pareto (tích lũy)
    geom_line(aes(x = cum_share * max(n) * 1.2, y = as.numeric(label)),
              color = "#2E86C1", linewidth = 1.1) +
    geom_point(aes(x = cum_share * max(n) * 1.2, y = as.numeric(label)),
               color = "#2E86C1", size = 2) +
    geom_vline(xintercept = cut_x, linetype = 3, color = "#566573") +
    scale_x_continuous(
      name = "Số khách hàng",
      labels = comma,
      sec.axis = dup_axis(
        name = "Tích lũy (%)",
        labels = function(z) percent(z / max(n) * 1.2, accuracy = 10)
      )
    ) +
    labs(
      title = "Phân bổ Tình trạng Sở hữu Nhà – Biểu đồ Pareto",
      subtitle = paste0("Tổng số khách hàng: ", format(total_n, big.mark = ","), 
                        " | Ngưỡng Pareto 80%"),
      y = NULL,
      caption = "Thanh = số lượng; Đường = tích lũy %. Vùng mờ biểu thị nhóm đóng góp 80% đầu."
    ) +
    theme_minimal(base_size = 30) +
    theme(
      plot.title = element_text(face = "bold", color = "#1B4F72", size=30),
      plot.subtitle = element_text(color = "#34495E", size =30),
      panel.grid.major.y = element_blank(),
      panel.grid.major.x = element_line(linewidth = 0.25, color = "grey85"),
      plot.caption = element_text(size = 30, color = "grey40")
    )
}

Khối code này sử dụng ggplot2 để tạo Biểu đồ Pareto cho biến Sở hữu nhà nhằm đánh giá sự phân bổ khối lượng khách hàng. Về mặt kỹ thuật và thống kê, mã xử lý các giá trị thiếu thành nhóm “Missing” (fct_explicit_na), sau đó tính toán Tỷ lệ tích lũy (cum_share) để xác định ngưỡng Pareto 80% (khối lượng cần tập trung).
Kết quả trực quan cho thấy sự phân bổ khối lượng không đồng đều: hai nhóm MORTGAGE (49.3%) và RENT (Thuê nhà) (39.7%) đã chiếm 89.0% tổng số khách hàng, vượt qua ngưỡng Pareto 80%. Nhóm OWN (Sở hữu nhà) chỉ chiếm 10.9%. Ý nghĩa thống kê của kết quả này là việc phân tích và mô hình hóa nên được ưu tiên tập trung vào hai nhóm khách hàng lớn nhất (MORTGAGE và RENT), bởi vì chúng đại diện cho gần 90% rủi ro và phơi nhiễm của danh mục. ### Phân tích mối quan hệ: Khoản vay (Loan Amount) và Thu nhập (Revenue)

Chunk này tạo ra một biểu đồ tán xạ (Scatter Plot) chuyên sâu, sử dụng kỹ thuật Heatmap 2D để trực quan hóa mối quan hệ giữa Số tiền vay (loan_amnt) và Thu nhập (revenue_w) của người vay, đồng thời ước tính độ mạnh của mối quan hệ này.

library(dplyr)
library(ggplot2)
library(scales)

# Tuỳ chọn
pal <- c("#E8F6F3","#9AC9C1","#4D8F86","#FFDC9E","#FF9F40","#E6B800","#D98800","#B35A00")

# Lấy mẫu dữ liệu và chuẩn bị biến
df <- sample_rows(df_cleaned, 200000) %>%
  transmute(
    revenue = as.numeric(revenue_w),
    loan = as.numeric(loan_amnt)
  ) %>%
  filter(is.finite(revenue), is.finite(loan))

# Tính toán thống kê
r_p <- cor(df$revenue, df$loan, method="pearson", use="complete.obs")
r_s <- cor(df$revenue, df$loan, method="spearman", use="complete.obs")
beta <- coef(lm(loan ~ revenue, df))[2]

# Đường median theo nhóm thu nhập
med_line <- df %>%
  mutate(g = ntile(revenue, 20)) %>%
  group_by(g) %>%
  summarise(x = median(revenue), y = median(loan), .groups="drop")

# Vẽ biểu đồ
ggplot(df, aes(revenue, loan)) +
  stat_bin2d(bins=60, aes(fill=after_stat(count))) +
  scale_fill_gradientn(colours = pal, name="Mật độ") +
  stat_density_2d(color="black", alpha=.2, linewidth=.3) +
  geom_path(data=med_line, aes(x, y), color=pal[8], linewidth=1) +
  geom_smooth(method="lm", se=FALSE, color=pal[5], linetype=2) +
  labs(
    title="Mối quan hệ giữa Thu nhập và Khoản vay",
    subtitle=paste0("r(Pearson)=", round(r_p,2),
                    " · ρ(Spearman)=", round(r_s,2),
                    " · β̂=", dollar(beta, .01)),
    x="Thu nhập (USD)", y="Khoản vay (USD)",
    caption="Heatmap: mật độ dữ liệu; vàng: median; nét đứt: hồi quy OLS."
  ) +
  theme_minimal(base_size=26)

Khối code này tạo biểu đồ Heatmap mật độ 2D để phân tích mối quan hệ giữa Thu nhập và Số tiền vay. Kỹ thuật thống kê bao gồm lấy mẫu dữ liệu và tính toán Hệ số tương quan Pearson (\(r = 0.49\)) và Spearman (\(\rho = 0.50\)). Kết quả xác nhận mối quan hệ tương quan dương vừa phải giữa hai biến. Mật độ tập trung cao ở vùng thu nhập và khoản vay thấp. Ý nghĩa thống kê quan trọng là: phân bố dữ liệu bị lệch phải, thể hiện qua việc đường trung vị (median) nằm cao hơn đường hồi quy OLS (\(\hat{\beta} \approx 0.10\) USD). Điều này ngụ ý rằng mô hình tuyến tính có thể đánh giá thấp xu hướng thực tế do ảnh hưởng của các giá trị ngoại lệ. ### Biểu đồ Ma trận Tương quan Spearman (Spearman Correlation Heatmap)

Chunk này tạo ra Biểu đồ Ma trận Tương quan (Correlation Matrix Plot) cho tất cả các biến số định lượng cốt lõi. Biểu đồ này giúp khám phá mối quan hệ tuyến tính/đơn điệu đa biến, đặc biệt quan trọng để đánh giá vấn đề đa cộng tuyến (multicollinearity) giữa các biến độc lập.

library(dplyr); library(tidyr); library(corrplot)

# ---- Chọn gradient ----
grad_name <- "Purple-Orange Strong"   # "Blue-Red Strong" | "Purple-Orange Strong"
pal <- switch(grad_name,
  "Blue-Red Strong"     = colorRampPalette(c("#2C7BB6","#92C5DE","#F7F7F7","#F4A582","#D6604D"))(200),
  "Purple-Orange Strong"= colorRampPalette(c("#5E3C99","#B2ABD2","#F7F7F7","#FDB863","#E66101"))(200)
)

# --- Chuẩn bị dữ liệu như trước ---
num_df <- df_cleaned |> dplyr::select(dplyr::where(is.numeric))
drop_id_like <- c("^id$", "^zip", "customer_", "client_")
num_df2 <- num_df |> dplyr::select(-tidyselect::matches(paste(drop_id_like, collapse="|"), ignore.case = TRUE))
qc <- purrr::map_dfr(num_df2, ~tibble::tibble(n_non_na=sum(!is.na(.x)), n_uniq=dplyr::n_distinct(.x, na.rm=TRUE), var_x=stats::var(.x, na.rm=TRUE)), .id="col")
good_cols <- qc |> dplyr::filter(n_uniq>=2, is.finite(var_x), var_x>0) |> dplyr::pull(col)
X <- dplyr::select(num_df2, tidyselect::any_of(good_cols))

if (ncol(X) >= 2) {
  cor_mat <- suppressWarnings(stats::cor(X, use="pairwise.complete.obs", method="spearman"))
  # p-value mask để ẩn ô không ý nghĩa
  p_mat <- matrix(NA_real_, ncol=ncol(X), nrow=ncol(X), dimnames=list(colnames(X), colnames(X)))
  for (i in seq_len(ncol(X))) for (j in i:ncol(X)) {
    xi <- X[[i]]; xj <- X[[j]]; ok <- is.finite(xi)&is.finite(xj)
    if (sum(ok)>=5) { pv <- suppressWarnings(cor.test(xi[ok], xj[ok], method="spearman")$p.value)
      p_mat[i,j] <- p_mat[j,i] <- pv }
  }

  corrplot::corrplot(
    cor_mat,
    method    = "color",
    type      = "lower",
    order     = "hclust",
    col       = pal,          # gradient diverging rõ ràng
    tl.cex    = 1.55, tl.srt = 45, tl.col = "grey20",
    diag      = FALSE,
    p.mat     = p_mat, sig.level = 0.05, insig = "blank", # ẩn ô không ý nghĩa
    na.label  = " ",
    cl.lim    = c(-1, 1)      # cố định thang -1..1
  )
}

Đoạn code thực hiện tính toán và trực quan hóa ma trận tương quan giữa các biến định lượng bằng hệ số tương quan Spearman, giúp phát hiện mối quan hệ tuyến tính và phi tuyến giữa các biến. Hàm cor() được dùng để tính hệ số tương quan, sau đó kết quả được trực quan bằng corrplot() nhằm thể hiện cả mức độ (màu đậm/nhạt) và chiều hướng tương quan (màu cam dương, tím âm). Các giá trị không có ý nghĩa thống kê (p-value > 0.05) được ẩn để đảm bảo độ tin cậy. Kết quả cho thấy các biến như revenue, revenue_w và loan_amnt có tương quan dương mạnh, phản ánh thu nhập cao thường gắn với quy mô vay lớn. Ngược lại, biến lti tương quan âm với default, cho thấy người có tỷ lệ thu nhập trên nợ cao hơn thường ít có khả năng vỡ nợ hơn. ### Biểu đồ Biplot PCA: Giải thích PC1 và PC2

Chunk này tạo ra biểu đồ Biplot (Biểu đồ tải trọng và điểm số) cho hai thành phần chính đầu tiên (PC1 và PC2), nhằm trực quan hóa cách các biến gốc đóng góp vào từng thành phần và cách dữ liệu được phân bố trên không gian 2 chiều mới này.

library(ggplot2)
library(scales)

pal <- list(bg = c("#EEF0FA","#DDE1F3","#BFC4E3","#9FA6D4","#6E79C0"),
            line = "#D9B53A", arrow = "#2C3E50", text = "#2C3E50")

if (exists("pca_res")) {
  scores <- as.data.frame(pca_res$x[, 1:2])
  loads  <- as.data.frame(pca_res$rotation[, 1:2]) |> tibble::rownames_to_column("var")
  imp <- summary(pca_res)$importance
  vx <- round(100 * imp["Proportion of Variance", 1:2], 1)

  mult <- 0.8 * min(
    max(abs(scores$PC1)) / max(abs(loads$PC1)),
    max(abs(scores$PC2)) / max(abs(loads$PC2))
  )
  loads <- loads |> mutate(PC1s = PC1 * mult, PC2s = PC2 * mult)

  ggplot(scores, aes(PC1, PC2)) +
    stat_bin2d(bins = 60, aes(fill = after_stat(count))) +
    scale_fill_gradientn(colours = pal$bg, name = "Count") +
    geom_segment(data = loads, aes(x = 0, y = 0, xend = PC1s, yend = PC2s),
                 arrow = arrow(length = unit(0.018, "npc")), color = pal$arrow) +
    ggrepel::geom_text_repel(data = loads, aes(PC1s, PC2s, label = var),
                             size = 15, color = pal$text, seed = 123) +
    labs(title = "PCA Biplot: PC1 vs PC2",
         x = paste0("PC1 (", vx[1], "%)"),
         y = paste0("PC2 (", vx[2], "%)")) +
    theme_minimal(base_size = 28)
}

Đoạn code trên thực hiện trực quan hóa kết quả phân tích thành phần chính (PCA) bằng biểu đồ biplot. Cụ thể, nó lấy điểm tọa độ của các quan sát (scores) và vector tải của biến gốc (loadings) từ đối tượng pca_res. Hai thành phần đầu tiên (PC1 và PC2) được chọn vì giải thích phần lớn phương sai dữ liệu. Để hiển thị cân đối, các vector tải được nhân với hệ số mult nhằm thu nhỏ độ dài cho phù hợp khung đồ thị. Biểu đồ được vẽ bằng ggplot2, trong đó các điểm được hiển thị bằng mật độ màu (heatmap) qua stat_bin2d, còn hướng và độ mạnh ảnh hưởng của từng biến được thể hiện qua các mũi tên và nhãn biến (geom_segment, geom_text_repel). Cách hiển thị này giúp người xem dễ dàng nhận diện nhóm biến cùng hướng (tương quan dương) hoặc ngược hướng (tương quan âm) trong không gian PCA. ### Biểu đồ Lollipop: Tỷ lệ vỡ nợ theo nhóm LTI (LTI vs Default Rate)

Chunk này tạo ra Biểu đồ Lollipop (Lollipop Plot) để trực quan hóa mối quan hệ đơn điệu giữa Tỷ lệ Vay/Thu nhập (lti_band) và Tỷ lệ Vỡ nợ (Default rate). Biểu đồ này tích hợp Khoảng Tin cậy 95% và Quy mô nhóm (Exposure) để đánh giá rủi ro toàn diện.

library(dplyr)
library(forcats)
library(stringr)
library(ggplot2)
library(scales)
library(viridis)

# Hàm Wilson CI (được định nghĩa lại bên trong chunk để tránh lỗi)
wilson_ci <- function(x, n, conf = 0.95) {
  z  <- qnorm(1 - (1 - conf)/2)
  p  <- ifelse(n > 0, x / n, NA_real_)
  den <- 1 + z^2 / n
  center <- (p + z^2/(2*n)) / den
  half   <- z * sqrt(p*(1 - p)/n + z^2/(4*n^2)) / den
  cbind(lo = pmax(0, center - half), hi = pmin(1, center + half))
}

if (need_cols(df_binned, c("lti_band","Default"))) { 

  # 1) Chuẩn hoá nhóm + sắp xếp
  dat <- df_binned %>%
    mutate(lti_band = fct_explicit_na(lti_band, na_level = "Missing")) 
  
  if (!is.ordered(dat$lti_band)) { 
    lv <- levels(dat$lti_band)
    ord <- suppressWarnings(as.numeric(stringr::str_extract(lv, "^-?\\d+")))
    if (any(!is.na(ord))) dat$lti_band <- factor(dat$lti_band, levels = lv[order(ord)])
  }

  # 3) Tóm tắt theo dải LTI
  sum_lti <- dat %>% 
    group_by(lti_band) %>% 
    summarise(
      n        = n(),
      defaults = sum(Default == 1, na.rm = TRUE),
      .groups  = "drop"
    ) %>%
    mutate(
      rate = ifelse(n > 0, defaults / n, NA_real_),
      exp  = n / sum(n)
    )

  # 4) CI + baseline + p-trend
  ci_mat <- wilson_ci(sum_lti$defaults, pmax(sum_lti$n, 1))
  sum_lti$ci_lo <- ci_mat[, "lo"]; sum_lti$ci_hi <- ci_mat[, "hi"]
  overall_rate  <- with(sum_lti, sum(defaults)/sum(n))
  p_trend <- try({
    ok <- !is.na(sum_lti$rate) & sum_lti$n > 0
    if (sum(ok) >= 3) suppressWarnings(prop.trend.test(
      x = sum_lti$defaults[ok], n = sum_lti$n[ok]
    )$p.value) else NA_real_
  }, silent = TRUE)
  p_trend_txt <- if (!inherits(p_trend, "try-error") && is.finite(p_trend)) {
    paste0("Cochran–Armitage p = ", formatC(p_trend, format = "e", digits = 2)) } else NULL

  # 5) Nhãn y gộp (bin + n + %exposure)
  sum_lti <- sum_lti %>% 
    mutate(
      y_lab = paste0(
        as.character(lti_band), 
        "  (n=", formatC(n, format = "d", big.mark = ","), "; ",
        percent(exp, accuracy = 0.1), ")"
      )
    )

  # 6) Vẽ-màu theo rate (viridis), size theo exposure
  xmax <- max(sum_lti$ci_hi, na.rm = TRUE)

  gp <- ggplot(sum_lti, aes(x = rate, y = fct_rev(y_lab))) + 
    # “Que kẹo”
    geom_segment(aes(x = 0, xend = rate, yend = fct_rev(y_lab)),
                 linewidth = 1, color = "grey70") +
    # CI 95%
    geom_errorbarh(aes(xmin = ci_lo, xmax = ci_hi),
                   height = 0.18, linewidth = 0.7, color = "grey40") +
    # Điểm trung tâm-shape 21 (viền đậm, fill theo rate), size theo exposure
    geom_point(aes(size = exp, fill = rate),
               shape = 21, color = "grey20", stroke = 0.5) +
    # Baseline toàn tập
    geom_vline(xintercept = overall_rate, linetype = "dashed",
               linewidth = 0.7, color = "grey35") +
    # Nhãn % tại điểm
    geom_text(aes(label = percent(rate, accuracy = 0.01),
                  x = pmin(rate + 0.02 * xmax, xmax)),
              hjust = 0, size = 10, color = "grey15") +
    scale_x_continuous(
      labels = percent_format(accuracy = 0.5),
      limits = c(0, xmax * 1.1),
      expand = expansion(mult = c(0, 0.03))
    ) +
    scale_fill_viridis(name = "Default rate", option = "C", direction = -1) +
    scale_size(name = "Exposure", range = c(2.5, 7),
               breaks = c(0.02, 0.05, 0.10, 0.20, 0.40),
               labels = percent_format(accuracy = 1)) +
    guides(fill = guide_colorbar(barheight = unit(40, "pt"), barwidth = unit(8, "pt")),
           size = guide_legend(override.aes = list(shape = 21, fill = "grey70"))) +
    labs(
      title = "LTI vs Default-Lollipop với 95% CI (Wilson)",
      subtitle = paste0(
        "Đường gạch: mức toàn tập = ", percent(overall_rate, accuracy = 0.01),
        if (!is.null(p_trend_txt)) paste0(" · ", p_trend_txt) else ""
      ),
      x = "Default rate",
      y = NULL,
      caption = "Màu tô = Default rate; kích thước = Exposure. CI 95% theo Wilson."
    ) +
    theme_classic(base_size = 30) +
    theme(
      plot.title.position = "plot",
      plot.title = element_text(face = "bold", size=30),
      axis.title.x = element_text(margin = margin(t = 6)),
      axis.text.y  = element_text(size = 26),
      panel.grid.major.x = element_line(color = "grey90", linewidth = 0.4),
      panel.grid.minor.x = element_blank(),
      legend.position = "right",
      legend.key.height = unit(12, "pt"),
      legend.key.width  = unit(12, "pt")
    )

  print(gp) # Thêm lệnh print(gp) để đảm bảo biểu đồ hiển thị
}

Đoạn code trên thực hiện việc phân tích mối quan hệ giữa tỷ lệ nợ trên thu nhập (LTI) và tỷ lệ vỡ nợ (Default rate) thông qua biểu đồ Lollipop kèm khoảng tin cậy 95% theo Wilson.
Cụ thể, code chia dữ liệu thành các nhóm LTI (lti_band), sau đó tính số lượng quan sát (n), số trường hợp vỡ nợ (defaults) và tỷ lệ vỡ nợ (rate) cho từng nhóm. Hàm wilson_ci() được sử dụng để tính khoảng tin cậy 95% cho tỷ lệ vỡ nợ — giúp biểu đồ thể hiện rõ độ tin cậy của từng điểm ước lượng. Biểu đồ minh họa tỷ lệ vỡ nợ bằng màu sắc (fill), quy mô nhóm vay bằng kích thước điểm (size), và đường gạch dọc biểu diễn mức trung bình toàn tập.
Ngoài ra, code còn thực hiện kiểm định xu hướng Cochran–Armitage để đánh giá xem tỷ lệ vỡ nợ có tăng dần theo các nhóm LTI hay không.
Về kết quả, biểu đồ cho thấy tỷ lệ vỡ nợ tăng dần khi LTI tăng, thể hiện rõ rủi ro tín dụng cao hơn ở các khoản vay có tỷ lệ nợ trên thu nhập lớn. Nhóm LTI > 30% có tỷ lệ vỡ nợ cao nhất (~28.64%), trong khi nhóm ≤5% thấp nhất (~11.67%). Điều này gợi ý rằng khả năng trả nợ của khách hàng giảm khi gánh nặng nợ tăng, phù hợp với kỳ vọng lý thuyết về rủi ro tín dụng. ### Biểu đồ Lollipop: Tỷ lệ vỡ nợ theo mục đích vay (Purpose vs Default Rate)

Chunk này tạo ra Biểu đồ Lollipop để trực quan hóa mức độ rủi ro vỡ nợ của từng nhóm Mục đích vay đã được nhóm nhỏ (lumped) và chuẩn hóa. Biểu đồ tích hợp Khoảng Tin cậy 95% và hiển thị cả Tỷ lệ vỡ nợ lẫn Quy mô nhóm (Exposure).

suppressPackageStartupMessages({
  library(dplyr); library(forcats); library(stringr)
  library(ggplot2); library(scales); library(viridis)
})

if (need_cols(df_cleaned, c("purpose","Default"))) {

  # 1) Chuẩn hoá purpose, gom top-8 + Other + Missing
  dat <- df_cleaned %>%
    mutate(
      purpose = tolower(trimws(as.character(purpose))),
      purpose = ifelse(purpose %in% c("", "na", "null"), NA, purpose),
      purpose_top = fct_lump_n(purpose, n = 8, other_level = "other"),
      purpose_top = fct_explicit_na(purpose_top, na_level = "missing")
    )

  # 2) Tổng hợp
  sum_pur <- dat %>%
    group_by(purpose_top) %>%
    summarise(
      n        = n(),
      defaults = sum(Default == 1, na.rm = TRUE),
      .groups  = "drop"
    ) %>%
    mutate(
      rate = defaults / pmax(n, 1),
      exp  = n / sum(n)
    )

  # 3) CI Wilson + baseline
  wilson_ci <- function(x, n, conf = .95){
    z <- qnorm(1 - (1 - conf)/2)
    p <- ifelse(n > 0, x/n, NA_real_)
    den <- 1 + z^2/n
    center <- (p + z^2/(2*n))/den
    half   <- z*sqrt(p*(1-p)/n + z^2/(4*n^2))/den
    cbind(
      lo = pmax(0, center - half),
      hi = pmin(1, center + half)
    )
  }

  ci <- wilson_ci(sum_pur$defaults, pmax(sum_pur$n,1))
  sum_pur$ci_lo <- ci[,"lo"]
  sum_pur$ci_hi <- ci[,"hi"]

  overall_rate  <- with(sum_pur, sum(defaults)/sum(n))

  # 4) Chuẩn hoá label hiển thị và trục Y
  to_title <- function(s) gsub("_"," ", tools::toTitleCase(s))

  sum_pur <- sum_pur %>%
    arrange(rate) %>%
    mutate(
      # Nhãn trục Y: Tên + (n = ...)
      y_lab = paste0(
        to_title(as.character(purpose_top)),
        "  (n=", formatC(n, format="d", big.mark=","), ")"
      ),
      # Text hiển thị bên phải mỗi điểm: "14.70% · Exp 1.1%"
      ann_text = paste0(
        percent(rate, accuracy = 0.01),
        " · ",
        "Exp ", percent(exp, accuracy = 0.1)
      )
    )

  # 5) Thông số cho layout
  xmax  <- max(sum_pur$ci_hi, na.rm = TRUE)
  # khoảng đẩy nhãn sang phải so với điểm
  x_pad <- 0.05 * xmax
  # nới phải để khỏi cắt nhãn
  x_limit_right <- xmax * 1.30

  gp <- ggplot(sum_pur, aes(x = rate, y = fct_rev(y_lab))) +

    # thanh lollipop
    geom_segment(
      aes(x = 0, xend = rate, yend = fct_rev(y_lab)),
      linewidth = 0.9, color = "grey70"
    ) +

    # CI ngang
    geom_errorbarh(
      aes(xmin = ci_lo, xmax = ci_hi),
      height = 0.16,
      linewidth = 0.6,
      color = "grey45"
    ) +

    # baseline toàn tập
    geom_vline(
      xintercept = overall_rate,
      linetype   = "dashed",
      linewidth  = 0.6,
      color      = "grey35"
    ) +

    # điểm trung tâm
    geom_point(
      aes(fill = rate),
      shape  = 21,
      size   = 3.8,
      color  = "grey20",
      stroke = 0.5
    ) +

    # nhãn duy nhất (rate + exposure), hơi nhích lên một chút để tránh đè vào điểm
    geom_text(
      aes(x = pmin(rate + x_pad, x_limit_right * 0.98),
          label = ann_text),
      hjust = 0,
      vjust = -0.3,
      size  = 10,
      color = "grey15"
    ) +

    # trục & màu
    scale_x_continuous(
      breaks  = pretty(c(0, xmax), n = 4),
      labels  = percent_format(accuracy = 1),
      limits  = c(0, x_limit_right),
      expand  = expansion(mult = c(0, 0.02))
    ) +
    scale_fill_viridis(
      name      = "Default rate",
      option    = "C",
      direction = -1
    ) +
    guides(
      fill = guide_colorbar(
        title.position   = "top",
        direction        = "horizontal",
        barwidth         = unit(120, "pt"),
        barheight        = unit(8, "pt"),
        ticks            = FALSE,
        label.position   = "bottom"
      )
    ) +

    labs(
      title = "Default theo mục đích vay-Lollipop 95% CI (Wilson)",
      subtitle = paste0(
        "Đường gạch: mức toàn tập = ",
        percent(overall_rate, accuracy = 0.01)
      ),
      x = "Default rate",
      y = NULL,
      caption = "Màu = Default rate (viridis); nhãn bên phải hiển thị Default% và Exposure%."
    ) +

    theme_classic(base_size = 30) +
    theme(
      plot.title.position = "plot",
      plot.title          = element_text(face = "bold",size=30),
      axis.title.x        = element_text(margin = margin(t = 6)),
      axis.text.y         = element_text(size = 26),
      panel.grid.major.x  = element_line(color = "grey90", linewidth = 0.35),
      panel.grid.minor.x  = element_blank(),
      legend.position     = "top",
      legend.justification= "left",
      legend.text         = element_text(size = 26),
      legend.title        = element_text(size = 26),
      plot.margin         = margin(t = 6, r = 24, b = 6, l = 6)
    ) +
    coord_cartesian(clip = "off")

  gp
}

Đoạn code trên dùng các gói dplyr, ggplot2, scales, và viridis để xử lý và trực quan hóa dữ liệu vỡ nợ theo mục đích vay. Trước tiên, biến purpose được chuẩn hóa, gộp thành top 8 nhóm phổ biến và nhóm “other”, sau đó tính số quan sát (n), số vỡ nợ (defaults), tỷ lệ vỡ nợ (rate) và tỷ trọng (exp). Hàm wilson_ci() được viết thủ công để tính khoảng tin cậy Wilson 95% cho từng nhóm.
Dữ liệu được sắp xếp và gắn nhãn hiển thị tỷ lệ và tỷ trọng. Biểu đồ Lollipop chart được tạo bằng ggplot2 với các lớp: geom_segment (thanh ngang), geom_errorbarh (CI), geom_point (điểm chính), geom_text (nhãn tỷ lệ) và geom_vline (đường trung bình toàn tập). Màu sắc tỷ lệ vỡ nợ được mã hóa bằng thang viridis, giúp trực quan mức độ rủi ro. Kết quả cho thấy “Small business” có tỷ lệ vỡ nợ cao nhất (≈29,54%), trong khi “Car” thấp nhất (≈14,65%), phản ánh chênh lệch rủi ro đáng kể giữa các mục đích vay. ### Biểu đồ Cột chồng: Khả năng Phân tách Rủi ro (Separability) theo Gánh nặng LTI

library(dplyr)
library(ggplot2)
library(scales)
library(forcats)
library(viridis)

# 1. Chuẩn bị dữ liệu: Dùng df_binned
df_sep_plot <- df_binned %>%
  filter(!is.na(Default), !is.na(lti_band)) %>%
  mutate(
    # Chuyển Default thành Factor để trực quan hóa
    default_status = factor(Default, levels = c(0, 1), labels = c("0: Current (Không vỡ nợ)", "1: Default (Vỡ nợ)")),
    # Sắp xếp LTI từ thấp đến cao
    lti_band_ordered = forcats::fct_relevel(lti_band, levels(df_binned$lti_band))
  )

# 2. Vẽ Biểu đồ Phân tách - Biểu đồ cột chồng
p_separability <- ggplot(df_sep_plot, aes(x = lti_band_ordered)) +
  
  # Layer 1: Vẽ Histogram (cột) chồng nhau cho Tỷ lệ %
  geom_bar(aes(fill = default_status), position = "fill", color = "white", linewidth = 0.2) +
  
  # Layer 2: Thêm nhãn Tỷ lệ % vỡ nợ (Default Rate)
  geom_text(
    stat = "count", 
    aes(y = after_stat(count) / tapply(after_stat(count), after_stat(x), sum)[after_stat(x)],
        label = scales::percent(after_stat(count) / tapply(after_stat(count), after_stat(x), sum)[after_stat(x)], 0.1)),
    position = position_stack(vjust = 0.5), # Đặt nhãn ở giữa
    color = "white", size = 3.5
  ) +
  
  # Layer 3: Thang màu tùy chỉnh (Xanh lá cho Current, Đỏ/Cam cho Default)
  scale_fill_manual(
    values = c("0: Current (Không vỡ nợ)" = "#1E8449", "1: Default (Vỡ nợ)" = "#C0392B"), 
    name = "Tình trạng khoản vay"
  ) +
  
  # Layer 4: Trục Y (Hiển thị dưới dạng %)
  scale_y_continuous(labels = scales::percent) +
  
  # Layer 5: Thiết lập Tiêu đề và Nhãn
  labs(
    title = "Khả năng Phân tách Rủi ro (Separability) theo Gánh nặng LTI",
    subtitle = "Phân tích tỷ lệ % vỡ nợ (Màu Đỏ) trong từng dải LTI. Kiểm định Giả thuyết H2.",
    x = "Nhóm LTI (Loan-to-Income)",
    y = "Tỷ lệ % khách hàng"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(face="bold"), legend.position = "bottom")

print(p_separability)

Khối code này tạo Biểu đồ Cột chồng 100% (Stacked Bar Chart) để trực quan hóa Khả năng Phân tách Rủi ro (Separability) của biến LTI band đối với biến mục tiêu (Default). Về mặt kỹ thuật và thống kê, mã chuyển biến Default thành kiểu Factor, sử dụng geom_bar(position = “fill”) để tính tỷ lệ phần trăm (thay vì số đếm tuyệt đối) của vỡ nợ (Màu Đỏ) trong từng dải LTI. Ý nghĩa thống kê của biểu đồ là kiểm tra trực quan mối quan hệ đơn điệu giữa biến dự báo và rủi ro. Kết quả trực quan xác nhận khả năng phân tách rủi ro tuyệt vời của LTI: Tỷ lệ vỡ nợ (Màu Đỏ) tăng đều và mạnh mẽ theo từng dải LTI, từ \(\approx 11.7\%\) ở dải thấp nhất (\(\le 5\%\)) lên \(\approx 28.6\%\) ở dải cao nhất (\(>30\%\)). Phân tích này chứng minh LTI band là một biến phân loại rủi ro rất mạnh và có thể dễ dàng tách biệt các nhóm rủi ro khác nhau.

Biểu đồ Mật độ điều kiện: Phân phối Khoản vay theo Thái cực LTI

library(dplyr); library(ggplot2); library(scales); library(forcats); library(viridis)

# 1. Chuẩn bị dữ liệu: Lọc chỉ lấy 2 nhóm LTI cực đoan
df_density <- df_binned %>%
  filter(lti_band %in% c("≤5%", ">30%")) %>%
  mutate(
    # Tạo nhãn và màu sắc: Low Risk (Green) vs High Risk (Red)
    risk_extreme = factor(lti_band, levels = c("≤5%", ">30%"), 
                          labels = c("Low Risk (LTI ≤ 5%)", "High Risk (LTI > 30%)"))
  )

# 2. Tính toán các điểm Median cho từng nhóm
median_loans <- df_density %>%
  group_by(risk_extreme) %>%
  summarise(
    median_loan = median(loan_amnt, na.rm = TRUE),
    .groups = 'drop'
  )

# 3. Tính toán Median toàn danh mục (cho đường tham chiếu)
overall_median <- median(df_binned$loan_amnt, na.rm = TRUE)

# 4. Vẽ Biểu đồ Mật độ (Density Plot)
p_density_cond <- ggplot(df_density, aes(x = loan_amnt, fill = risk_extreme)) +
  
  # Layer 1 & 2: Vẽ các đường mật độ (Distribution)
  geom_density(alpha = 0.5, linewidth = 1) + 
  
  # Layer 3: Đường tham chiếu Median Toàn danh mục
  geom_vline(xintercept = overall_median, linetype = "dotted", color = "grey50", linewidth = 1) +
  
  # Layer 4 & 5: Đường tham chiếu Median cho từng nhóm Cực đoan
  geom_vline(data = median_loans, aes(xintercept = median_loan, color = risk_extreme), 
             linetype = "dashed", linewidth = 1) + 
  
  # Layer 6: Thiết lập Trục và Thang màu
 scale_x_continuous(labels = comma_format(scale = 0.001, suffix = "K"), limits = c(0, 40000)) +
  scale_fill_manual(values = c("Low Risk (LTI ≤ 5%)" = "#1E8449", "High Risk (LTI > 30%)" = "#C0392B"), name = "Nhóm Rủi ro LTI") +
  scale_color_manual(values = c("Low Risk (LTI ≤ 5%)" = "#1E8449", "High Risk (LTI > 30%)" = "#C0392B"), name = "Nhóm Rủi ro LTI") +
  
  # Layer 7: Thiết lập Tiêu đề
  labs(
    title = "Phân phối Quy mô Khoản vay theo Thái cực Rủi ro LTI",
    subtitle = "So sánh Loan Amount giữa khách hàng Rủi ro Thấp nhất (LTI ≤ 5%) và Cao nhất (LTI > 30%).",
    x = "Số tiền vay (USD)",
    y = "Mật độ"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(face="bold"), legend.position = "top")

print(p_density_cond)

Khối code này tạo Biểu đồ Mật độ điều kiện (Conditional Density Plot) nhằm so sánh sự phân phối của Số tiền vay (loan_amnt) giữa hai nhóm rủi ro LTI cực đoan: Low Risk (\(\le 5\%\)) và High Risk (\(>30\%\)).
Về mặt kỹ thuật và thống kê, mã lọc dữ liệu chỉ lấy hai nhóm này và tính Trung vị (Median) cho từng nhóm, sau đó sử dụng geom_density để trực quan hóa mật độ.Kết quả trực quan xác nhận sự khác biệt rõ rệt về quy mô khoản vay giữa hai thái cực rủi ro:
Low Risk (Màu Xanh lá): Đường mật độ tập trung cao nhất ở mức khoản vay thấp (dưới $10K), và Trung vị của nhóm này cũng nằm ở mức thấp (đường đứt nét màu xanh).
High Risk (Màu Đỏ): Đường mật độ trải rộng hơn nhiều và có các đỉnh ở mức khoản vay cao hơn, cho thấy nhóm rủi ro cao nhất có xu hướng vay số tiền lớn hơn và có độ biến thiên lớn về quy mô khoản vay. Trung vị của nhóm này nằm ở mức cao hơn đáng kể (đường đứt nét màu đỏ), vượt xa trung vị của nhóm Low Risk.
Ý nghĩa thống kê: Biểu đồ này chứng minh rằng rủi ro cao không chỉ liên quan đến thu nhập thấp, mà còn liên quan đến quy mô khoản vay lớn hơn, đặc biệt là trong bối cảnh gánh nặng nợ cao (LTI cao).

Biểu đồ phân phối mật đồ gánh nặng

library(dplyr); library(ggplot2); library(scales); library(forcats); library(viridis)

# 1. Chuẩn bị dữ liệu: Gom nhóm Purpose và lọc outlier LTI
df_lti_density <- df_binned %>%
  filter(!is.na(lti), !is.na(purpose), lti < quantile(lti, 0.99)) %>%
  mutate(
    # Lấy 5 mục đích chính (tạo thành các đường cong riêng)
    purpose_group = fct_lump_n(purpose, n = 5, other_level = "Other Purposes")
  )

# 2. Tính toán Median LTI toàn danh mục (đường tham chiếu)
overall_median_lti <- median(df_lti_density$lti, na.rm = TRUE)

# 3. Vẽ Biểu đồ Mật độ Có điều kiện
p_lti_density <- ggplot(df_lti_density, aes(x = lti, fill = purpose_group, color = purpose_group)) +
  
  # Layer 1: Vẽ các đường mật độ (Density)
  geom_density(alpha = 0.4, linewidth = 1) + 
  
  # Layer 2: Đường tham chiếu Median LTI Toàn danh mục
  geom_vline(xintercept = overall_median_lti, linetype = "dashed", color = "black", linewidth = 0.8) +
  
  # Layer 3: Thêm nhãn cho đường Median Toàn danh mục
  annotate("text", x = overall_median_lti + 0.05, y = Inf, vjust = 1.5, color = "black", size = 3.5, 
           label = paste("Median Toàn danh mục:", scales::percent(overall_median_lti, 0.1))) +

  # Layer 4: Thiết lập Trục và Thang màu
  scale_x_continuous(labels = percent) +
  scale_y_continuous(labels = comma) +
  scale_fill_viridis_d(option = "Plasma", name = "Mục đích vay") +
  scale_color_viridis_d(option = "Plasma", name = "Mục đích vay") +
  
  # Layer 5: Thiết lập Tiêu đề
  labs(
    title = "Phân phối Mật độ Gánh nặng LTI theo Mục đích vay",
    subtitle = "So sánh hình dạng và vị trí của chỉ số LTI giữa các mục đích vay chính.",
    x = "Tỷ lệ LTI (Loan-to-Income)",
    y = "Mật độ"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(face="bold"), legend.position = "top")

print(p_lti_density)

Khối code này tạo Biểu đồ Mật độ có điều kiện (geom_density) để so sánh sự phân phối của Tỷ lệ LTI giữa các mục đích vay chính. Về mặt kỹ thuật và thống kê, mã lọc các outlier LTI cực đoan (trên phân vị 99%) để làm sạch biểu đồ, gom nhóm các mục đích vay nhỏ thành “Other Purposes” (fct_lump_n), và vẽ biểu đồ với các đường mật độ chồng lên nhau (alpha=0.4). Mã tính và thêm đường tham chiếu Trung vị Toàn danh mục (19.5%). Kết quả trực quan xác nhận sự khác biệt đáng kể về gánh nặng nợ giữa các nhóm:
* Gánh nặng cao nhất: Mục đích debt_consolidation (Màu Tím đậm) có đường mật độ dịch chuyển mạnh nhất sang phải (về phía LTI cao hơn) và tập trung cao ở vùng LTI cao (trên 20%), cho thấy nhóm này có gánh nặng nợ lớn nhất, vượt xa Trung vị Toàn danh mục.
* Gánh nặng tthấp nhất: Các mục đích other và home_improvement (Màu Xanh lá/Vàng) có mật độ tập trung cao nhất ở vùng LTI thấp (dưới 15%).
Ý nghĩa Thống kê: Sự khác biệt về hình dạng và vị trí của các đường mật độ chứng minh rằng Mục đích vay là một biến có giá trị trong việc phân khúc mức độ gánh nặng nợ và, do đó, là một chỉ báo rủi ro quan trọng.

Biểu đồ Boxplot: Phân phối rủi ro theo gánh nặng LTI

library(dplyr); library(ggplot2); library(scales); library(forcats)
library(viridis)

# 1. Chuẩn bị dữ liệu: Lấy sample và thêm nhiễu (jitter) vào biến Default
# Tạo một cột y_jittered để biểu diễn 0/1 dưới dạng boxplot
df_plot_default <- df_binned %>%
  filter(!is.na(Default), !is.na(lti_band)) %>%
  mutate(
    # Thêm nhiễu nhẹ để các điểm 0/1 không bị chồng lên nhau
    default_jittered = Default + runif(n(), -0.1, 0.1) 
  )

# 2. Tính toán điểm Median Rủi ro (đường chấm)
median_risk <- df_plot_default %>%
  group_by(lti_band) %>%
  summarise(
    median_rate = mean(Default), # Median của biến 0/1 chính là Tỷ lệ Vỡ nợ (Mean)
    .groups = 'drop'
  )

# 3. Vẽ Boxplot Rủi ro
p_risk_boxplot <- ggplot(df_plot_default, aes(x = lti_band, y = default_jittered, fill = lti_band)) +
  
  # Layer 1: Vẽ Boxplot cho Default (thể hiện sự phân phối rủi ro)
  geom_boxplot(alpha = 0.6, outlier.shape = NA) +
  
  # Layer 2: Vẽ đường tham chiếu Tỷ lệ Vỡ nợ Toàn danh mục
  geom_hline(yintercept = mean(df_plot_default$Default), linetype = "dashed", color = "red", linewidth = 0.8) +
  
  # Layer 3: Vẽ điểm Tỷ lệ Vỡ nợ Median (Median Rate)
  geom_point(data = median_risk, aes(y = median_rate), 
             shape = 23, fill = "black", color = "white", size = 3) +
  
  # Layer 4: Thêm nhãn Tỷ lệ Vỡ nợ (%) cho mỗi Boxplot
  geom_text(data = median_risk, aes(y = median_rate, 
                                   label = scales::percent(median_rate, 0.1)),
            vjust = -1.5, size = 3.5, color = "black") +
  
  # Layer 5: Thiết lập Trục (chỉ hiện 0/1)
  scale_y_continuous(breaks = c(0, 1), labels = c("Không Vỡ nợ (0)", "Vỡ nợ (1)")) +
  scale_fill_viridis_d(option = "D", name = "Nhóm LTI") +
  
  # Layer 6: Thiết lập Tiêu đề và Theme
  labs(
    title = "Phân phối Rủi ro (Default) theo Gánh nặng LTI",
    subtitle = "Mỗi Boxplot cho thấy sự phân tách (separability) giữa 0 và 1 theo dải LTI. Đường đỏ = Mức nền.",
    x = "Nhóm LTI (Loan-to-Income)",
    y = "Tình trạng Vỡ nợ",
    caption = "Hộp càng thấp/hẹp, nhóm càng an toàn."
  ) +
  theme_minimal() +
  theme(legend.position = "none", plot.title = element_text(face="bold"))

print(p_risk_boxplot)

Khối code này tạo Biểu đồ Boxplot Rủi ro để trực quan hóa khả năng phân tách rủi ro của biến LTI band trên biến mục tiêu nhị phân (Default). Về mặt kỹ thuật, mã xử lý biến nhị phân bằng cách thêm nhiễu nhẹ (jittering) (default_jittered) để Boxplot có thể hiển thị phân phối, và sử dụng Trung bình (mean(Default)) để tính Tỷ lệ Vỡ nợ cho mỗi dải LTI. Kết quả trực quan xác nhận khả năng phân tách rủi ro cực kỳ cao của LTI: Tỷ lệ Vỡ nợ tăng đều từ mức thấp nhất 11.7% (nhóm \(\le 5\%\)) lên mức cao nhất 28.6% (nhóm \(>30\%\)), thể hiện mối quan hệ tỷ lệ thuận mạnh mẽ. Sự dịch chuyển vị trí rõ rệt của các hộp Boxplot ra xa mức 0 (Không Vỡ nợ) chứng minh LTI band là một biến phân loại rủi ro rất mạnh, có khả năng phân biệt hiệu quả các nhóm rủi ro thấp, trung bình và cao, làm cơ sở vững chắc cho mô hình. # CHƯƠNG 2: PHÂN TÍCH BCTC CTCP VÀNG BẠC ĐÁ QUÝ PHÚ NHUẬN PNJPNJ

Bài phân tích này thực hiện đánh giá định lượng về hiệu suất tài chính và hoạt động kinh doanh của CTCP Vàng bạc Đá quý Phú Nhuận (PNJ), dựa trên dữ liệu Báo cáo Tài chính Quý 1/2017 – Quý 2/2025.
Mục tiêu chính là phân tích cơ cấu tài chính và các nhân tố thúc đẩy lợi nhuận của PNJ, tập trung trả lời bốn câu hỏi cốt lõi:
Động lực nào chi phối xu hướng tăng trưởng của công ty?
Việc quản trị vốn lưu động và sử dụng đòn bẩy tài chính đã hiệu quả như thế nào?
Các yếu tố cấu thành khả năng sinh lời trên vốn chủ sở hữu (ROE) là gì (Phân tích Du Pont)?
Chất lượng lợi nhuận có tương xứng với dòng tiền từ hoạt động kinh doanh (CFO) không?

Tải dữ liệu thô và kiểm tra tiền phân tích

data_path2 <- if (exists("params") && !is.null(params$data_path2)) params$data_path2 else
  "BCTC_PNJ.xlsx"
if (!file.exists(data_path2)) stop("Không tìm thấy file: ", data_path2)
message("[INFO] Đọc dữ liệu từ: ", data_path2)  # hoặc cat("[INFO] ...\n")
#------

data_path2 <- if (exists("params") && !is.null(params$data_path2)) params$data_path2 else 
  "BCTC_PNJ.xlsx"
#-----
if (!file.exists(data_path2)) stop("Không tìm thấy file: ", data_path2)
df_raw <- readxl::read_excel(path = data_path2, sheet = 1)
dplyr::glimpse(df_raw)

## Rows: 17
## Columns: 35
## $ Indicator <chr> "ta", "equity", "tl", "ca", "cl", "inv_net", "ar", "ap", "re…
## $ `Q1 2017` <dbl> 3.595808e+12, 1.749066e+12, 1.846742e+12, 3.106734e+12, 1.78…
## $ `Q2 2017` <dbl> 3.561007e+12, 1.735442e+12, 1.825565e+12, 3.055472e+12, 1.77…
## $ `Q3 2017` <dbl> 4.063461e+12, 2.728456e+12, 1.335006e+12, 3.480268e+12, 1.31…
## $ `Q4 2017` <dbl> 4.492513e+12, 2.950301e+12, 1.542211e+12, 3.896409e+12, 1.52…
## $ `Q1 2018` <dbl> 4.688396e+12, 3.198904e+12, 1.489493e+12, 3.872677e+12, 1.40…
## $ `Q2 2018` <dbl> 5.173682e+12, 3.366712e+12, 1.806970e+12, 4.233071e+12, 1.74…
## $ `Q3 2018` <dbl> 5.481233e+12, 3.382261e+12, 2.098972e+12, 4.491050e+12, 2.07…
## $ `Q4 2018` <dbl> 6.303185e+12, 3.745313e+12, 2.557872e+12, 5.280216e+12, 2.54…
## $ `Q1 2019` <dbl> 6.152152e+12, 4.040218e+12, 2.111934e+12, 5.065795e+12, 2.09…
## $ `Q2 2019` <dbl> 6.526498e+12, 4.092534e+12, 2.433964e+12, 5.333230e+12, 2.42…
## $ `Q3 2019` <dbl> 7.625552e+12, 4.120300e+12, 3.505252e+12, 6.396513e+12, 3.49…
## $ `Q4 2019` <dbl> 8.600312e+12, 4.574035e+12, 4.026277e+12, 7.330560e+12, 4.01…
## $ `Q1 2020` <dbl> 8.389422e+12, 4.985528e+12, 3.403894e+12, 7.097183e+12, 3.39…
## $ `Q2 2020` <dbl> 8.162334e+12, 4.701173e+12, 3.461161e+12, 6.863106e+12, 3.45…
## $ `Q3 2020` <dbl> 8.091042e+12, 4.903259e+12, 3.187783e+12, 6.786222e+12, 3.17…
## $ `Q4 2020` <dbl> 8.483371e+12, 5.241862e+12, 3.241509e+12, 7.144154e+12, 3.23…
## $ `Q1 2021` <dbl> 8.136212e+12, 5.752949e+12, 2.383263e+12, 6.821148e+12, 2.37…
## $ `Q2 2021` <dbl> 9.175654e+12, 5.718861e+12, 3.456793e+12, 7.870492e+12, 3.44…
## $ `Q3 2021` <dbl> 9168881692481, 5559332736154, 3609548956327, 7887307408212, …
## $ `Q4 2021` <dbl> 1.054664e+13, 6.016510e+12, 4.530132e+12, 9.220118e+12, 4.52…
## $ `Q1 2022` <dbl> 1.141021e+13, 8.008233e+12, 3.401972e+12, 1.009283e+13, 3.32…
## $ `Q2 2022` <dbl> 1.101943e+13, 8.029628e+12, 2.989799e+12, 9.707247e+12, 2.98…
## $ `Q3 2022` <dbl> 1.246321e+13, 8.318218e+12, 4.144991e+12, 1.115783e+13, 4.13…
## $ `Q4 2022` <dbl> 1.332110e+13, 8.587800e+12, 4.733304e+12, 1.195792e+13, 4.72…
## $ `Q1 2023` <dbl> 1.283083e+13, 9.192688e+12, 3.638140e+12, 1.143193e+13, 3.62…
## $ `Q2 2023` <dbl> 1.349285e+13, 9.194055e+12, 4.298800e+12, 1.208256e+13, 4.28…
## $ `Q3 2023` <dbl> 1.305516e+13, 9.436999e+12, 3.618164e+12, 1.161576e+13, 3.60…
## $ `Q4 2023` <dbl> 1.442994e+13, 9.806566e+12, 4.623377e+12, 1.296011e+13, 4.61…
## $ `Q1 2024` <dbl> 1.296893e+13, 1.047438e+13, 2.494545e+12, 1.153522e+13, 2.48…
## $ `Q2 2024` <dbl> 1.296771e+13, 1.071481e+13, 2.252909e+12, 1.153787e+13, 2.24…
## $ `Q3 2024` <dbl> 1.496249e+13, 1.052443e+13, 4.438053e+12, 1.353531e+13, 4.42…
## $ `Q4 2024` <dbl> 1.720773e+13, 1.125531e+13, 5.952424e+12, 1.569260e+13, 5.94…
## $ `Q1 2025` <dbl> 1.741932e+13, 1.173029e+13, 5.689026e+12, 1.594966e+13, 5.67…
## $ `Q2 2025` <dbl> 1.715378e+13, 1.196938e+13, 5.184392e+12, 1.571393e+13, 5.17…

Code này đọc và kiểm tra cấu trúc một file Excel chứa báo cáo tài chính.

• Tạo biến data_path2 để lưu đường dẫn file. Điều này giúp Ta chỉ cần sửa một nơi khi đường dẫn thay đổi.
• Dùng file.exists() để xác nhận file tồn tại. Nếu không, stop() dừng chương trình và báo lỗi để ngăn các vấn đề phát sinh sau này.
• Đọc dữ liệu từ sheet đầu tiên bằng hàm readxl::read_excel. Kết quả được lưu vào biến df_raw, cho biết đây là dữ liệu thô.
• Dùng dplyr::glimpse() để xem tóm tắt cấu trúc dữ liệu. Hàm này hiển thị thông tin rõ ràng khi có nhiều cột.

Kết quả cho thấy dữ liệu có cấu trúc dạng rộng (wide format).

• Bảng dữ liệu có 17 dòng và 35 cột.
• Mỗi dòng là một chỉ số tài chính, xác định bởi cột Indicator.
• Các cột còn lại, như Q1 2017, chứa giá trị của các chỉ số theo từng quý.

Tái cấu trúc dữ liệu

Đoạn code này thực hiện việc tái định hình dữ liệu Báo cáo Tài chính sang định dạng chuỗi thời gian tiêu chuẩn: * Tách ngữ cảnh và chuyển vị: Dùng column_to_rownames để đặt tên các chỉ tiêu BCTC làm nhãn hàng, sau đó dùng lệnh t() để hoán đổi vị trí (Transpose). Kết quả: Các kỳ báo cáo (Quý/Năm) trở thành hàng (quan sát) và các chỉ số tài chính trở thành cột (biến số). * Ép kiểu toàn diện: Lệnh lapply(…, as.numeric) chuyển đổi tất cả các cột thành kiểu số (numeric). Đây là bước kiểm soát chất lượng quan trọng, đảm bảo dữ liệu sẵn sàng cho mọi phép tính và mô hình thống kê. * Phục hồi Nhãn hàng: Gán lại tên hàng (là các quý) để duy trì ngữ cảnh thời gian cho mỗi quan sát. Cấu trúc dữ liệu mới đã sẵn sàng cho phân tích: Mỗi hàng là một thời điểm; mỗi cột là một biến số tài chính định lượng.

Chuẩn hóa đơn vị

# Đơn vị (Ngàn tỷ)
df_numeric <- df_numeric / 1e12

# In 5 dòng đầu tiên của dữ liệu đã chuẩn hóa để kiểm tra
print(t(head(df_numeric, 1)))

##              Q1 2017
## ta       3.595808173
## equity   1.749065808
## tl       1.846742365
## ca       3.106733885
## cl       1.784913164
## inv_net  2.846786699
## ar       0.040388070
## ap       0.323535911
## revenue  3.130971058
## cogs    -2.580113338
## gp       0.550857720
## ebit     0.311240665
## ebt      0.311075206
## ni       0.248739097
## cfo      0.204604720
## cfi     -0.007393507
## cff     -0.225897519

Code này chuẩn hóa đơn vị dữ liệu.

• Chia toàn bộ dữ liệu cho 1 nghìn tỷ (1e12).
• Thao tác này đổi đơn vị các giá trị từ Đồng sang Nghìn tỷ Đồng.
• R thực hiện phép tính này trên mọi ô dữ liệu cùng lúc. Điều này an toàn vì Ta đã đảm bảo tất cả các cột đều là số.

Feature Engineering

df_analysis <- df_numeric

# cogs nên dương (magnitude)
if ("cogs" %in% names(df_analysis)) {
  df_analysis$cogs <- abs(df_analysis$cogs)
  message("[INFO] Đã chuẩn hóa 'cogs' sang giá trị dương.")
} else {
  warning("[WARN] Không tìm thấy cột 'cogs'.")
}

## [INFO] Đã chuẩn hóa 'cogs' sang giá trị dương.

# Các biến nên dương tuyệt đối
positive_cols <- c('ta', 'equity', 'tl', 'ca', 'cl', 'inv_net', 'ar', 'ap', 'revenue')
for (col in positive_cols) {
  if (col %in% names(df_analysis)) {
    df_analysis[[col]] <- abs(df_analysis[[col]])
  } else {
    warning(paste("[WARN] Không tìm thấy cột '", col, "' để chuẩn hóa dương.", sep=""))
  }
}
message("[INFO] Hoàn tất chuẩn hóa nền tảng.")

## [INFO] Hoàn tất chuẩn hóa nền tảng.

Code này làm sạch dữ liệu bằng cách áp dụng các quy tắc tài chính.

• Tạo một bản sao dữ liệu tên df_analysis. Hành động này bảo vệ dữ liệu gốc của Ta.
• Chuyển đổi giá trị của cột cogs thành số dương bằng hàm abs(). Dữ liệu kế toán thường ghi cogs là số âm, nhưng phân tích cần độ lớn thực của chi phí.
• Định nghĩa một danh sách các cột phải có giá trị dương, ví dụ như ta (tổng tài sản).
• Code dùng một vòng lặp để áp dụng hàm abs() cho tất cả các cột trong danh sách đó.

Các hàm tiện ích

# Hàm tiện ích dùng xuyên suốt
avg_lag1 <- function(x) {
  (x + dplyr::lag(x, n = 1)) / 2
}

Code này tạo ra các hàm để tiện tính toán trong quá trình phân tích

• Tạo hàm avg_lag1 để tính trung bình của một giá trị tại kỳ hiện tại và kỳ trước đó.
• Hàm dùng dplyr::lag() để lấy giá trị của kỳ trước, sau đó tính trung bình cộng của giá trị hiện tại và giá trị đó.

Mục đích:

• Ta không thể so sánh một số liệu đo lường trong một kỳ (như lợi nhuận) với một số liệu đo tại một thời điểm (như tài sản).
• Giải pháp là Ta phải dùng giá trị tài sản trung bình trong kỳ đó.
• Hàm avg_lag1 tính toán chính xác giá trị trung bình này, giúp các tỷ số tài chính chính xác hơn.

Biên lợi nhuận

# 1) Khả năng sinh lời — Biên lợi nhuận
if (all(c("gp","ebit","ni","revenue") %in% names(df_analysis))) {
  df_analysis <- df_analysis %>%
    dplyr::mutate(
      gross_margin = gp / revenue,
      operating_margin = ebit / revenue,
      net_margin = ni / revenue
    )
  message("[INFO] Đã tính biên lợi nhuận (gross/operating/net).")
} else {
  warning("[WARN] Thiếu một trong gp/ebit/ni/revenue — bỏ qua margins.")
}

## [INFO] Đã tính biên lợi nhuận (gross/operating/net).

Code này tính toán các biên lợi nhuận để đánh giá khả năng sinh lời của công ty.

• Ta dùng all() để đảm bảo các cột gp, ebit, ni, và revenue đều tồn tại.
• Ta dùng dplyr::mutate để tạo ba cột mới: gross_margin, operating_margin, và net_margin.

Mỗi biên lợi nhuận đo lường hiệu quả ở một giai đoạn khác nhau

• Biên lợi nhuận gộp: Đo lường hiệu quả sản xuất.
• Biên lợi nhuận hoạt động: Đo lường hiệu quả kinh doanh cốt lõi.
• Biên lợi nhuận ròng: Đo lường lợi nhuận cuối cùng.

Biến trung bình

# Biến trung bình cho ROA/ROE/Efficiency
need <- c("ta","equity","inv_net","ar")
if (any(need %in% names(df_analysis))) {
  df_analysis <- df_analysis %>%
    dplyr::mutate(
      avg_ta = avg_lag1(ta),
      avg_equity = avg_lag1(equity),
      avg_inv_net = avg_lag1(inv_net),
      avg_ar = avg_lag1(ar)
    )
  message("[INFO] Đã tạo các biến trung bình: avg_ta, avg_equity, avg_inv_net, avg_ar.")
}

## [INFO] Đã tạo các biến trung bình: avg_ta, avg_equity, avg_inv_net, avg_ar.

Code này tạo ra các cột dữ liệu trung bình để chuẩn bị cho các phép tính tỷ số.

• Áp dụng hàm avg_lag1 đã tạo trước đó.
• Hàm này tạo ra các cột mới như avg_ta và avg_equity, chứa giá trị trung bình của chỉ số gốc và chỉ số kỳ trước.
• Hành động này giải quyết vấn đề so sánh giữa hai loại dữ liệu tài chính.
• Ta không thể chia một chỉ số của cả kỳ (lợi nhuận) cho một chỉ số tại một thời điểm (tài sản).
• Các cột trung bình mới này làm cho các phép tính tỷ số của Ta chính xác về mặt logic.

Tính ROA & ROE

# ROA & ROE
if (all(c("ni","avg_ta","avg_equity") %in% names(df_analysis))) {
  df_analysis <- df_analysis %>%
    dplyr::mutate(
      roa = ni / avg_ta,
      roe = ni / avg_equity
    )
  message("[INFO] Đã tính ROA, ROE.")
} else {
  warning("[WARN] Thiếu ni/avg_ta/avg_equity — bỏ qua ROA/ROE.")
}

## [INFO] Đã tính ROA, ROE.

Code này tính Tỷ suất sinh lời trên tài sản (ROA) và Tỷ suất sinh lời trên vốn chủ sở hữu (ROE).

• Ta dùng mutate để tạo hai cột mới: roa và roe.
• Công thức sử dụng lợi nhuận ròng (ni) và các giá trị tài sản (avg_ta) hoặc vốn chủ sở hữu (avg_equity) trung bình.
• ROA đo lường hiệu quả Ta sử dụng toàn bộ tài sản để sinh lời.
• ROE đo lường lợi nhuận Ta tạo ra cho mỗi đồng vốn của cổ đông.

ROE có thể được phân tích thành ba yếu tố chính.

• Biên lợi nhuận: Khả năng sinh lời trên doanh thu.
• Vòng quay tài sản: Hiệu quả sử dụng tài sản để tạo doanh thu.
• Đòn bẩy tài chính: Mức độ sử dụng nợ để tài trợ cho tài sản.

Tính các giá trị YOY (theo quý)

# 2) Tăng trưởng YOY (theo quý)
if (all(c("revenue","gp","ebit","ni") %in% names(df_analysis))) {
  df_analysis <- df_analysis %>%
    dplyr::mutate(
      revenue_yoy = (revenue / dplyr::lag(revenue, 4)) - 1,
      gp_yoy      = (gp       / dplyr::lag(gp, 4)) - 1,
      ebit_yoy    = (ebit     / dplyr::lag(ebit, 4)) - 1,
      ni_yoy      = (ni       / dplyr::lag(ni, 4)) - 1
    )
  message("[INFO] Đã tính tăng trưởng YOY (revenue/gp/ebit/ni).")
}

## [INFO] Đã tính tăng trưởng YOY (revenue/gp/ebit/ni).

Code này tính toán tăng trưởng so với cùng kỳ năm trước (Year-over-Year).

• Ta dùng dplyr::lag(x, 4) để lấy dữ liệu từ 4 quý trước.
• Con số 4 được chọn vì dữ liệu là theo quý.
• Công thức (hiện tại / quá khứ) - 1 chuyển đổi tỷ lệ thành phần trăm tăng trưởng.

Phép tính này loại bỏ yếu tố mùa vụ.

• So sánh quý 4 với quý 3 thường không chính xác. Doanh thu quý 4 thường cao hơn do các yếu tố mùa vụ.
• So sánh quý 4 năm nay với quý 4 năm ngoái cho thấy mức tăng trưởng thực sự của công ty.
• Công thức tính toán là: (Giá trị hiện tại / Giá trị cùng kỳ năm trước) - 1.

Tỷ số Vòng quay

# 3) Hiệu quả — Vòng quay
if (all(c("revenue","avg_ta","cogs","avg_inv_net","avg_ar") %in% names(df_analysis))) {
  df_analysis <- df_analysis %>%
    dplyr::mutate(
      asset_turnover       = revenue / avg_ta,
      inventory_turnover   = cogs / avg_inv_net,
      receivables_turnover = revenue / avg_ar
    )
  message("[INFO] Đã tính asset/inventory/receivables turnover.")
}

## [INFO] Đã tính asset/inventory/receivables turnover.

Code này tính các tỷ số vòng quay để đo lường hiệu quả hoạt động.

• Vòng quay tổng tài sản: Ta tính bằng revenue chia cho avg_ta. Chỉ số này đo lường Ta tạo ra bao nhiêu doanh thu từ mỗi đồng tài sản.
• Vòng quay hàng tồn kho: Ta tính bằng cogs chia cho avg_inv_net. Ta phải dùng cogs (giá vốn) vì hàng tồn kho được ghi nhận theo giá vốn, không phải giá bán.
• Vòng quay các khoản phải thu: Ta tính bằng revenue chia cho avg_ar. Chỉ số này đo tốc độ Ta thu tiền từ khách hàng.

Chỉ số Hiệu quả

# 3) Hiệu quả — DSO/DIO/DPO/CCC
if (all(c("ar","revenue","inv_net","cogs","ap") %in% names(df_analysis))) {
  df_analysis <- df_analysis %>%
    dplyr::mutate(
      dso = 365 * ar / revenue,
      dio = 365 * inv_net / cogs,
      dpo = 365 * ap / cogs,
      ccc = dso + dio - dpo
    )
  message("[INFO] Đã tính DSO/DIO/DPO/CCC.")
}

## [INFO] Đã tính DSO/DIO/DPO/CCC.

Code này tính toán chu kỳ chuyển đổi tiền mặt.

• Tính Số ngày phải thu bình quân (DSO). Công thức này cho biết Ta mất trung bình bao nhiêu ngày để thu tiền từ khách hàng.
• Tính Số ngày tồn kho bình quân (DIO). Nó đo lường thời gian trung bình để Ta bán hết hàng trong kho.
• Tính Số ngày phải trả bình quân (DPO). Nó cho thấy Ta có bao nhiêu ngày để trả tiền cho nhà cung cấp.
• Tính Chu kỳ chuyển đổi tiền mặt (CCC). Đây là thời gian ròng mà tiền của Ta bị kẹt trong hoạt động kinh doanh, tính bằng công thức DSO + DIO - DPO

Đòn bẩy tài chính

# 4) Đòn bẩy tài chính
if (all(c("tl","equity","ta") %in% names(df_analysis))) {
  df_analysis <- df_analysis %>%
    dplyr::mutate(
      debt_to_equity    = tl / equity,
      financial_leverage = ta / equity
    )
  message("[INFO] Đã tính debt_to_equity, financial_leverage.")
}

## [INFO] Đã tính debt_to_equity, financial_leverage.

Code này tính toán các tỷ số đòn bẩy tài chính để đánh giá rủi ro.

• Tính tỷ số Nợ trên Vốn chủ sở hữu (debt_to_equity). Công thức này cho biết Ta có bao nhiêu đồng nợ cho mỗi đồng vốn chủ sở hữu.
• Tính Đòn bẩy tài chính (financial_leverage). Nó đo lường Ta kiểm soát bao nhiêu đồng tài sản cho mỗi đồng vốn chủ sở hữu.

Các công thức này dựa trên phương trình kế toán:

• Tài sản = Nợ + Vốn chủ sở hữu.
• Đòn bẩy tài chính = Tỷ số Nợ trên Vốn chủ sở hữu + 1.

Thanh khoản

# 5) Thanh khoản
if (all(c("ca","cl","inv_net") %in% names(df_analysis))) {
  df_analysis <- df_analysis %>%
    dplyr::mutate(
      current_ratio = ca / cl,
      quick_ratio   = (ca - inv_net) / cl
    )
  message("[INFO] Đã tính current_ratio, quick_ratio.")
}

## [INFO] Đã tính current_ratio, quick_ratio.

Code này tính toán các tỷ số thanh khoản để đánh giá khả năng trả nợ ngắn hạn.

• Ta tính Tỷ số thanh khoản hiện thời (current_ratio). Công thức này cho biết Ta có bao nhiêu đồng tài sản ngắn hạn để trả cho mỗi đồng nợ ngắn hạn.
• Ta tính Tỷ số thanh khoản nhanh (quick_ratio). Tỷ số này loại bỏ hàng tồn kho. Hàng tồn kho là tài sản khó chuyển thành tiền nhất. Do đó, đây là một bài kiểm tra nghiêm ngặt hơn về khả năng trả nợ của Ta.

So sánh hai tỷ số này để đánh giá rủi ro. Nếu current_ratio cao nhưng quick_ratio thấp, thanh khoản của Ta phụ thuộc vào việc bán hàng tồn kho.

Dòng tiền

# 6) Dòng tiền 
if (all(c("cfo","ni","cfi","revenue") %in% names(df_analysis))) {
  df_analysis <- df_analysis %>%
    dplyr::mutate(
      cfo_ni     = cfo / ni,
      capex_proxy = -pmin(cfi, 0, na.rm = TRUE),
      fcf_proxy   = cfo - capex_proxy,
      fcf_margin  = fcf_proxy / revenue,
      fcf_margin_proxy_crude = (cfo + cfi) / revenue
    )
  message("[INFO] Đã tính CFO/NI, CAPEX proxy, FCF & FCF margin.")
}

## [INFO] Đã tính CFO/NI, CAPEX proxy, FCF & FCF margin.

Code này tính các chỉ số dòng tiền, bao gồm Dòng tiền tự do (FCF).

• Tính tỷ số cfo / ni. Tỷ số này đo lường chất lượng lợi nhuận bằng cách so sánh tiền mặt thực tế với lợi nhuận trên giấy tờ.
• Ước tính Chi tiêu Vốn (CAPEX) bằng công thức -pmin(cfi, 0). Công thức này lọc ra các khoản chi đầu tư âm từ dòng tiền đầu tư (cfi) và đổi nó thành số dương.
• Tính Dòng tiền tự do (FCF) bằng cfo - capex_proxy. Đây là lượng tiền còn lại sau khi trừ chi phí hoạt động và đầu tư.
• Tính Biên Dòng tiền tự do (fcf_margin). Chỉ số này cho biết Ta tạo ra bao nhiêu tiền mặt tự do từ mỗi đồng doanh thu.

EBITDA Margin

#  Conditional — EBITDA margin
if ("ebitda" %in% names(df_analysis) && "revenue" %in% names(df_analysis)) {
  df_analysis <- df_analysis %>% dplyr::mutate(ebitda_margin = ebitda / revenue)
  message("[INFO] Đã tính ebitda_margin.")
} else {
  warning("[WARN] Thiếu ebitda hoặc revenue — bỏ qua ebitda_margin.")
}

Code này tính Biên lợi nhuận EBITDA.

• Dùng mutate để tạo cột ebitda_margin bằng ebitda chia cho revenue.
• Hành động này cho Ta một thước đo lợi nhuận.
• Thước đo này loại bỏ ảnh hưởng từ các quyết định về tài chính, thuế và khấu hao.

So sánh các biên lợi nhuận khác nhau.

• Khoảng cách giữa biên EBITDA và biên hoạt động thể hiện tác động của chi phí khấu hao.
• Khoảng cách giữa biên hoạt động và biên ròng thể hiện tác động của lãi vay và thuế.

ROIC

#  Conditional — ROIC (cần ic & tax_rate)
if (all(c("ic","tax_rate","ebit") %in% names(df_analysis))) {
  df_analysis <- df_analysis %>% dplyr::mutate(
    avg_ic = avg_lag1(ic),
    roic   = ebit * (1 - tax_rate) / avg_ic
  )
  message("[INFO] Đã tính ROIC.")
} else {
  warning("[WARN] Thiếu ic/tax_rate/ebit — bỏ qua ROIC.")
}

Code này tính Tỷ suất sinh lời trên Vốn đầu tư (ROIC) để Ta đo lường khả năng tạo ra giá trị.

• Ta tính Lợi nhuận Hoạt động Ròng Sau Thuế (NOPAT). Công thức là EBIT * (1 - tax_rate).
• Ta tính Vốn đầu tư trung bình (avg_ic). Việc này đảm bảo Ta so sánh lợi nhuận trong kỳ với vốn hoạt động trong kỳ.
• Ta chia NOPAT cho vốn đầu tư trung bình để có ROIC.

Ta phải so sánh ROIC với Chi phí vốn bình quân gia quyền (WACC).

• ROIC > WACC: Công ty tạo ra giá trị.
• ROIC < WACC: Công ty phá hủy giá trị.

Net Debt / EBITDA

#  Conditional — Net Debt / EBITDA
req <- c("st_debt","lt_debt","cash_eq","ebitda")
if (all(req %in% names(df_analysis))) {
  df_analysis <- df_analysis %>% dplyr::mutate(
    net_debt        = st_debt + lt_debt - cash_eq,
    net_debt_ebitda = net_debt / ebitda
  )
  message("[INFO] Đã tính net_debt và net_debt_ebitda.")
} else {
  warning("[WARN] Thiếu một trong st_debt/lt_debt/cash_eq/ebitda — bỏ qua NetDebt/EBITDA.")
}

Code này tính toán tỷ số Nợ ròng trên EBITDA. Ta dùng nó để đánh giá khả năng trả nợ.

• Tính Nợ ròng (net_debt) bằng tổng nợ trừ đi tiền mặt. Công thức này phản ánh gánh nặng nợ thực tế vì Ta có thể dùng tiền mặt để trả nợ ngay lập tức.
• Chia Nợ ròng cho EBITDA. Kết quả cho thấy Ta cần bao nhiêu năm lợi nhuận hoạt động để trả hết nợ.

Tỷ số bảo đảm trả lãi vay

#  Conditional — Interest Coverage
if (all(c("ebit","int_exp") %in% names(df_analysis))) {
  df_analysis <- df_analysis %>% dplyr::mutate(interest_coverage = ebit / abs(int_exp))
  message("[INFO] Đã tính interest_coverage.")
} else {
  warning("[WARN] Thiếu ebit hoặc int_exp — bỏ qua interest_coverage.")
}

Code này tính Tỷ số Bảm đảm Trả lãi vay. Dùng nó để đánh giá mức độ an toàn của công ty.

• Tính tỷ số này bằng ebit chia cho abs(int_exp).
• Phải dùng ebit (lợi nhuận trước lãi vay). Đây là khoản lợi nhuận có sẵn để Ta trả lãi.
• Dùng abs() vì chi phí lãi vay thường là số âm. Hàm này lấy giá trị dương của nó.

Tỷ số tiền mặt

#  Conditional — Cash Ratio
if (all(c("cash_eq","cl") %in% names(df_analysis))) {
  df_analysis <- df_analysis %>% dplyr::mutate(cash_ratio = cash_eq / cl)
  message("[INFO] Đã tính cash_ratio.")
} else {
  warning("[WARN] Thiếu cash_eq hoặc cl — bỏ qua cash_ratio.")
}

Code này tính Tỷ số Tiền mặt. Ta dùng nó để đánh giá khả năng trả nợ trong kịch bản xấu nhất.

• Tính tỷ số này bằng cách chia tiền mặt (cash_eq) cho nợ ngắn hạn (cl).
• Nó chỉ dùng tài sản có thể sử dụng ngay lập tức, bỏ qua các khoản phải thu và hàng tồn kho.

Tỷ số này là bài kiểm tra thanh khoản nghiêm ngặt nhất. Nó nằm trong một hệ thống gồm ba cấp độ:

• Tỷ số thanh khoản hiện thời: Đánh giá chung.
• Tỷ số thanh khoản nhanh: Loại bỏ hàng tồn kho.
• Tỷ số Tiền mặt: Chỉ dùng tiền mặt.

Hàm an toàn

#  Dọn dẹp — thay Inf bằng NA
df_analysis[] <- lapply(df_analysis, function(x) {
  x[is.infinite(x)] <- NA
  x
})
message("[INFO] Hoàn tất tạo biến và dọn dẹp giá trị Inf.")

## [INFO] Hoàn tất tạo biến và dọn dẹp giá trị Inf.

Code này dọn dẹp các giá trị Inf (vô cùng) trong dữ liệu.

• Giá trị Inf xuất hiện khi Ta thực hiện phép chia cho 0.
• Giá trị này gây lỗi cho các phân tích thống kê và thuật toán học máy.
• Code này quét qua toàn bộ dữ liệu của Ta. Nó tìm tất cả các giá trị Inf và thay thế chúng bằng NA.

Từ điển biến

# 1) Tự đảm bảo có hàm ptbl()
if (!exists("ptbl")) {
  suppressPackageStartupMessages(library(kableExtra))
  ptbl <- function(x, caption = NULL, digits = 3) {
    stopifnot(is.data.frame(x) || inherits(x, "matrix"))
    fmt <- if (knitr::is_latex_output()) "latex" else "html"
    kb <- knitr::kable(x, format = fmt, caption = caption, digits = digits,
                       booktabs = TRUE, longtable = TRUE, linesep = "")
    if (knitr::is_latex_output()) {
      kb |>
        kableExtra::kable_styling(
          latex_options = c("hold_position", "repeat_header", "scale_down"),
          position = "center", font_size = 10
        )
    } else {
      kb |>
        kableExtra::kable_styling(full_width = FALSE,
                                  bootstrap_options = c("striped","condensed")) |>
        kableExtra::scroll_box(width = "100%")
    }
  }
}

# 2) Tạo từ điển biến 
suppressPackageStartupMessages(library(tibble))
dict <- tibble::tribble(
  ~Biến,                ~Mô_tả,                         ~Nhóm,
  "revenue",            "Doanh thu thuần",              "Input",
  "gp",                 "Lợi nhuận gộp",                "Input",
  "ebit",               "EBIT",                         "Input",
  "ni",                 "LN ròng",                      "Input",
  "ta",                 "Tổng tài sản",                 "Input",
  "equity",             "Vốn chủ sở hữu",               "Input",
  "cogs",               "Giá vốn hàng bán",             "Input",
  "inv_net",            "HTK ròng",                     "Input",
  "ar",                 "Phải thu KH",                  "Input",
  "ap",                 "Phải trả NB",                  "Input",
  "ca",                 "TS ngắn hạn",                  "Input",
  "cl",                 "Nợ ngắn hạn",                  "Input",
  "cfo",                "LCTT HĐKD",                    "Input",
  "cfi",                "LCTT HĐĐT",                    "Input",
  "ebitda",             "EBITDA",                       "Optional",
  "st_debt",            "Nợ vay ngắn hạn",              "Optional",
  "lt_debt",            "Nợ vay dài hạn",               "Optional",
  "cash_eq",            "Tiền & TĐT",                   "Optional",
  "tax_rate",           "Thuế suất hiệu dụng",          "Optional",
  "gross_margin",       "Biên gộp",                     "Derived",
  "operating_margin",   "Biên EBIT",                    "Derived",
  "net_margin",         "Biên ròng",                    "Derived",
  "avg_ta",             "TS bình quân",                 "Derived",
  "avg_equity",         "VCSH bình quân",               "Derived",
  "roa",                "ROA",                          "Derived",
  "roe",                "ROE",                          "Derived",
  "asset_turnover",     "Vòng quay TS",                 "Derived",
  "inventory_turnover", "Vòng quay HTK",                "Derived",
  "receivables_turnover","Vòng quay PT",                "Derived",
  "dso",                "DSO (ngày)",                   "Derived",
  "dio",                "DIO (ngày)",                   "Derived",
  "dpo",                "DPO (ngày)",                   "Derived",
  "ccc",                "Chu kỳ tiền (CCC)",            "Derived",
  "debt_to_equity",     "Nợ/VCSH",                      "Derived",
  "current_ratio",      "Current ratio",                "Derived",
  "quick_ratio",        "Quick ratio",                  "Derived",
  "cfo_ni",             "CFO/NI",                       "Derived",
  "fcf_margin",         "Biên FCF",                     "Derived"
)

# 3) In bảng ngay trong chunk
ptbl(dict, caption = "Từ điển biến", digits = 2)

Từ điển biến

Biến	Mô_tả	Nhóm
revenue	Doanh thu thuần	Input
gp	Lợi nhuận gộp	Input
ebit	EBIT	Input
ni	LN ròng	Input
ta	Tổng tài sản	Input
equity	Vốn chủ sở hữu	Input
cogs	Giá vốn hàng bán	Input
inv_net	HTK ròng	Input
ar	Phải thu KH	Input
ap	Phải trả NB	Input
ca	TS ngắn hạn	Input
cl	Nợ ngắn hạn	Input
cfo	LCTT HĐKD	Input
cfi	LCTT HĐĐT	Input
ebitda	EBITDA	Optional
st_debt	Nợ vay ngắn hạn	Optional
lt_debt	Nợ vay dài hạn	Optional
cash_eq	Tiền & TĐT	Optional
tax_rate	Thuế suất hiệu dụng	Optional
gross_margin	Biên gộp	Derived
operating_margin	Biên EBIT	Derived
net_margin	Biên ròng	Derived
avg_ta	TS bình quân	Derived
avg_equity	VCSH bình quân	Derived
roa	ROA	Derived
roe	ROE	Derived
asset_turnover	Vòng quay TS	Derived
inventory_turnover	Vòng quay HTK	Derived
receivables_turnover	Vòng quay PT	Derived
dso	DSO (ngày)	Derived
dio	DIO (ngày)	Derived
dpo	DPO (ngày)	Derived
ccc	Chu kỳ tiền (CCC)	Derived
debt_to_equity	Nợ/VCSH	Derived
current_ratio	Current ratio	Derived
quick_ratio	Quick ratio	Derived
cfo_ni	CFO/NI	Derived
fcf_margin	Biên FCF	Derived

Code này tạo một từ điển dữ liệu để Ta ghi lại ý nghĩa của tất cả các biến.

• Dùng tibble::tribble để tạo từ điển vì cú pháp của nó rõ ràng. Một hàm khác in từ điển ra thành một bảng chuyên nghiệp.
• Phân loại các biến thành ba nhóm: Input (dữ liệu gốc), Optional (dữ liệu có thể thiếu), và Derived (các biến Ta đã tạo ra).
• Cấu trúc này ghi lại toàn bộ quy trình phân tích của Ta, từ nguyên liệu thô đến sản phẩm cuối cùng.

Ổn định từ điển

# 1) Tự đảm bảo có hàm ptbl()
if (!exists("ptbl")) {
  suppressPackageStartupMessages(library(kableExtra))
  ptbl <- function(x, caption = NULL, digits = 3) {
    stopifnot(is.data.frame(x) || inherits(x, "matrix"))
    fmt <- if (knitr::is_latex_output()) "latex" else "html"
    kb <- knitr::kable(x, format = fmt, caption = caption, digits = digits,
                       booktabs = TRUE, longtable = TRUE, linesep = "")
    if (knitr::is_latex_output()) {
      kb |>
        kableExtra::kable_styling(
          latex_options = c("hold_position", "repeat_header", "scale_down"),
          position = "center", font_size = 10
        )
    } else {
      kb |>
        kableExtra::kable_styling(full_width = FALSE,
                                  bootstrap_options = c("striped","condensed")) |>
        kableExtra::scroll_box(width = "100%")
    }
  }
}

# 2) Tự đảm bảo có từ điển 
if (!exists("dict")) {
  suppressPackageStartupMessages(library(tibble))
  dict <- tibble::tribble(
    ~Biến, ~Mô_tả, ~Nhóm,
    "revenue","Doanh thu thuần","Input",
    "gp","Lợi nhuận gộp","Input",
    "ebit","EBIT","Input",
    "ni","LN ròng","Input",
    "ta","Tổng tài sản","Input",
    "equity","Vốn chủ sở hữu","Input",
    "cogs","Giá vốn hàng bán","Input",
    "inv_net","HTK ròng","Input",
    "ar","Phải thu KH","Input",
    "ap","Phải trả NB","Input",
    "ca","TS ngắn hạn","Input",
    "cl","Nợ ngắn hạn","Input",
    "cfo","LCTT HĐKD","Input",
    "cfi","LCTT HĐĐT","Input",
    "gross_margin","Biên gộp","Derived",
    "operating_margin","Biên EBIT","Derived",
    "net_margin","Biên ròng","Derived",
    "roa","ROA","Derived",
    "roe","ROE","Derived"
  )
}

# 3) Xác định dữ liệu tham chiếu
suppressPackageStartupMessages(library(dplyr))
data_ref <- if (exists("df_analysis")) df_analysis else if (exists(
"df_numeric")) df_numeric else df_raw
vars_avail <- intersect(dict$Biến, names(data_ref))

# 4) Lọc & in
if (length(vars_avail) == 0) {
  message("[INFO] Không thấy biến nào trong dữ liệu khớp từ điển.")
} else {
  dict_in_use <- dict |>
    dplyr::filter(Biến %in% vars_avail) |>
    dplyr::arrange(Nhóm, Biến) |>
    dplyr::mutate(`Đã có` = "\u2713")
  ptbl(dict_in_use, caption = "Danh mục biến dùng trong phân tích")
}

Danh mục biến dùng trong phân tích

Biến	Mô_tả	Nhóm	Đã có
asset_turnover	Vòng quay TS	Derived	✓
avg_equity	VCSH bình quân	Derived	✓
avg_ta	TS bình quân	Derived	✓
ccc	Chu kỳ tiền (CCC)	Derived	✓
cfo_ni	CFO/NI	Derived	✓
current_ratio	Current ratio	Derived	✓
debt_to_equity	Nợ/VCSH	Derived	✓
dio	DIO (ngày)	Derived	✓
dpo	DPO (ngày)	Derived	✓
dso	DSO (ngày)	Derived	✓
fcf_margin	Biên FCF	Derived	✓
gross_margin	Biên gộp	Derived	✓
inventory_turnover	Vòng quay HTK	Derived	✓
net_margin	Biên ròng	Derived	✓
operating_margin	Biên EBIT	Derived	✓
quick_ratio	Quick ratio	Derived	✓
receivables_turnover	Vòng quay PT	Derived	✓
roa	ROA	Derived	✓
roe	ROE	Derived	✓
ap	Phải trả NB	Input	✓
ar	Phải thu KH	Input	✓
ca	TS ngắn hạn	Input	✓
cfi	LCTT HĐĐT	Input	✓
cfo	LCTT HĐKD	Input	✓
cl	Nợ ngắn hạn	Input	✓
cogs	Giá vốn hàng bán	Input	✓
ebit	EBIT	Input	✓
equity	Vốn chủ sở hữu	Input	✓
gp	Lợi nhuận gộp	Input	✓
inv_net	HTK ròng	Input	✓
ni	LN ròng	Input	✓
revenue	Doanh thu thuần	Input	✓
ta	Tổng tài sản	Input	✓

Code này tạo một từ điển dữ liệu tự động. Nó kiểm tra dữ liệu của Ta và tạo một báo cáo về các biến hiện có.

• Code tự động tìm phiên bản dữ liệu đầy đủ nhất mà Ta có (df_analysis, df_numeric, hoặc df_raw).
• Nó so sánh từ điển gốc với dữ liệu này để tìm ra các biến thực sự tồn tại.
• Bảng kết quả hiển thị các biến đã tồn tại trong bộ dữ liệu, được sắp xếp và có dấu xác nhận.

Thống kê mô tả

# 1) Tự đảm bảo có hàm ptbl()
if (!exists("ptbl")) {
  suppressPackageStartupMessages(library(kableExtra))
  ptbl <- function(x, caption = NULL, digits = 3) {
    stopifnot(is.data.frame(x) || inherits(x, "matrix"))
    fmt <- if (knitr::is_latex_output()) "latex" else "html"
    kb <- knitr::kable(x, format = fmt, caption = caption, digits = digits,
                       booktabs = TRUE, longtable = TRUE, linesep = "")
    if (knitr::is_latex_output()) {
      kb |>
        kableExtra::kable_styling(
          latex_options = c("hold_position", "repeat_header", "scale_down"),
          position = "center", font_size = 10
        )
    } else {
      kb |>
        kableExtra::kable_styling(full_width = FALSE,
                                  bootstrap_options = c("striped","condensed")) |>
        kableExtra::scroll_box(width = "100%")
    }
  }
}

# 2) Tự đảm bảo có dict & data_ref
if (!exists("dict")) {
  suppressPackageStartupMessages(library(tibble))
  dict <- tibble::tribble(
    ~Biến, ~Mô_tả, ~Nhóm,
    "revenue","Doanh thu thuần","Input",
    "gp","Lợi nhuận gộp","Input",
    "ebit","EBIT","Input",
    "ni","LN ròng","Input",
    "ta","Tổng tài sản","Input",
    "equity","Vốn chủ sở hữu","Input",
    "cogs","Giá vốn hàng bán","Input",
    "inv_net","HTK ròng","Input",
    "ar","Phải thu KH","Input",
    "ap","Phải trả NB","Input",
    "ca","TS ngắn hạn","Input",
    "cl","Nợ ngắn hạn","Input"
  )
}

data_ref <- if (exists("df_analysis")) df_analysis else if (
  exists("df_numeric")) df_numeric else df_raw
vars_avail <- intersect(dict$Biến, names(data_ref))
num_vars <- vars_avail[sapply(data_ref[vars_avail], is.numeric)]

# 3) Tính & in
if (length(num_vars) > 0) {
  M <- sapply(data_ref[num_vars], function(x) c(
    n = sum(!is.na(x)),
    mean = mean(x, na.rm = TRUE),
    sd = stats::sd(x, na.rm = TRUE),
    min = suppressWarnings(min(x, na.rm = TRUE)),
    max = suppressWarnings(max(x, na.rm = TRUE))
  ))
  summ <- data.frame(Biến = colnames(M), t(M), row.names = NULL, check.names = FALSE)
  ptbl(summ, caption = "Thống kê nhanh các biến số (đơn giản)")
} else {
  message("[INFO] Không có biến số nào để tóm tắt.")
}

Thống kê nhanh các biến số (đơn giản)

Biến	n	mean	sd	min	max
revenue	34	5.933	2.842	0.877	12.594
gp	34	1.103	0.520	0.156	2.149
ebit	34	0.445	0.269	-0.193	0.941
ni	34	0.352	0.216	-0.160	0.749
ta	34	9.739	4.092	3.561	17.419
equity	34	6.465	3.110	1.735	11.969
cogs	34	4.830	2.337	0.721	10.445
inv_net	34	7.449	3.105	2.847	13.709
ar	34	0.069	0.034	0.027	0.200
ap	34	0.385	0.173	0.128	0.823
ca	34	8.535	3.851	3.055	15.950
cl	34	3.254	1.230	1.313	5.942
cfo	34	0.044	0.707	-1.573	1.887
cfi	34	-0.058	0.282	-0.965	0.767
gross_margin	34	0.186	0.016	0.156	0.219
operating_margin	34	0.067	0.055	-0.221	0.114
net_margin	34	0.053	0.045	-0.182	0.090
avg_ta	33	9.720	3.953	3.578	17.314
avg_equity	33	6.453	3.015	1.742	11.850
roa	33	0.038	0.019	-0.017	0.073
roe	33	0.059	0.030	-0.028	0.110
asset_turnover	33	0.627	0.170	0.096	0.924
inventory_turnover	33	0.669	0.192	0.096	1.039
receivables_turnover	33	91.904	41.495	24.624	164.907
dso	34	5.170	3.535	2.225	18.570
dio	34	668.030	573.493	324.351	3801.405
dpo	34	35.762	25.327	8.158	157.608
ccc	34	637.439	553.836	305.407	3655.352
debt_to_equity	34	0.567	0.194	0.210	1.056
current_ratio	34	2.618	0.730	1.723	5.143
quick_ratio	34	0.319	0.218	0.078	0.814
cfo_ni	34	-0.095	1.663	-3.509	2.565
fcf_margin	34	-0.031	0.120	-0.400	0.168

Code này tính toán và trình bày một bảng thống kê mô tả cho dữ liệu.

• Nó tự động tìm các biến dạng số hiện có trong dữ liệu của Ta.
• Nó dùng hàm sapply để tính toán số lượng (n), trung bình (mean), độ lệch chuẩn (sd), giá trị nhỏ nhất (min), và lớn nhất (max) cho mỗi biến.
• Bảng kết quả giúp Ta nhanh chóng hiểu được xu hướng trung tâm, mức độ biến động và phạm vi giá trị của từng biến. Bảng thống kê này cho Ta những thông tin quan trọng.
• n: Cho biết mức độ đầy đủ của dữ liệu.
• mean: Cho thấy giá trị trung tâm.
• sd: Đo lường sự biến động. sd lớn cho thấy dữ liệu phân tán rộng.
• min/max: Xác định khoảng giá trị. Chúng có thể giúp Ta phát hiện các điểm bất thường, ví dụ như một quý bị lỗ khi min là số âm.

Cấu hình vẽ biểu đồ

Code này thiết lập một môi trường để Ta vẽ các biểu đồ chuyên nghiệp và nhất quán.

• Ta tạo các công tắc điều khiển. Các biến use_plotly và save_plots cho phép Ta bật chế độ tương tác hoặc tự động lưu file. Ta thay đổi chúng ở một nơi duy nhất.
• Ta tạo các hàm định dạng. Chúng tự động định dạng các con số và tỷ lệ phần trăm trên biểu đồ, đảm bảo tính nhất quán.
• Ta xây dựng một theme tùy chỉnh. Hàm theme_pro định nghĩa phông chữ, màu sắc và bố cục. Ta dùng theme_set để áp dụng nó làm mặc định cho mọi biểu đồ.
• Ta chuẩn bị một data frame mới. Ta chuyển thông tin thời gian từ rownames thành một cột Date thực sự. ggplot2 yêu cầu cột Date này để vẽ biểu đồ chuỗi thời gian chính xác.

# Chuẩn bị dữ liệu cho ggplot: rownames -> Quarter, parse Date từ "Q1 2017"
df_plot <- df_analysis %>%
  tibble::rownames_to_column(var = "Quarter") %>%
  mutate(
    # Tạo cột Date an toàn từ các định dạng phổ biến: "Q1 2017", "2017 Q1", "2017Q1"
    Date = {
      d1 <- suppressWarnings(as.Date(zoo::as.yearqtr(
      Quarter, format = "Q%q %Y")))
      
      d2 <- suppressWarnings(as.Date(zoo::as.yearqtr(
      stringr::str_replace(Quarter, "^(\\d{4})Q(\\d)$", "Q\\2 \\1"),
                                                    format = "Q%q %Y")))
      
      d3 <- suppressWarnings(as.Date(zoo::as.yearqtr(
      stringr::str_replace(Quarter, "^(\\d{4})\\s*Q(\\d)$", "Q\\2 \\1"),
                                                    format = "Q%q %Y")))
      
      dplyr::coalesce(d1, d2, d3)
    },
    Year = factor(stringr::str_extract(Quarter, "\\d{4}$")), # Cũng sửa ở đây
    qtr  = as.integer(((as.integer(format(Date, "%m")) - 1) %/% 3) + 1),
    yr   = as.integer(format(Date, "%Y"))
  ) %>%
  arrange(Date)

message("Đã tạo 'df_plot' sẵn sàng cho trực quan hóa.")

theme_pro <- function() {
  theme_minimal(base_family = "sans", base_size = 12) +
  theme(
    plot.title = element_text(face = "bold", size = 16, color = "#222222"),
    plot.subtitle = element_text(size = 12, color = "#555555", margin = margin(b = 10)),
    plot.caption = element_text(size = 9, color = "#888888", face = "italic", hjust = 0),
    axis.title = element_text(face = "bold", size = 11, color = "#333333"),
    axis.text = element_text(size = 10, color = "#444444"),
    legend.position = "bottom",
    legend.title = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.major.y = element_line(color = "#E5E7EB", linetype = "dashed"),
    plot.background = element_rect(fill = "#FAFAFA", color = NA),
    panel.background = element_rect(fill = "#FAFAFA", color = NA),
    strip.text = element_text(face = "bold", size = 12, hjust = 0, color = "#333333"),
    plot.margin = margin(15, 15, 15, 15)
  )
}

theme_set(theme_pro())

# Breaks theo năm + xoay nhãn
date_break_years <- function(n = 1) list(
  scale_x_date(date_breaks = paste0(n, " year"), date_labels = "%Y"),
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
)
add_xdate <- function(p, n = 1) { if (inherits(p, 'gg')) p + date_break_years(n) else p }

# Hệ số trục phụ an toàn
safe_coef <- function(lhs, rhs) {
  a <- suppressWarnings(max(lhs, na.rm = TRUE));
  b <- suppressWarnings(max(rhs, na.rm = TRUE));
  if (!is.finite(a) || !is.finite(b) || b == 0) return(1)
  a / b
}

cap_default <- "Đơn vị: Ngàn Tỷ Đồng đối với số tuyệt đối; % với tỷ lệ. Dữ liệu đã được chuẩn hóa."
num_fmt  <- scales::label_number(big.mark = ".", decimal.mark = ",")
pct_fmt  <- scales::label_percent(accuracy = 0.1, big.mark = ".", decimal.mark = ",")
qty_fmt  <- scales::label_number(accuracy = 0.1, big.mark = ",", decimal.mark = ",")

Code này thiết lập một môi trường để vẽ các biểu đồ chuyên nghiệp và nhất quán.

• Ta chuẩn bị một data frame mới. Ta chuyển thông tin thời gian từ rownames thành một cột Date thực sự. ggplot2 yêu cầu cột Date này để vẽ biểu đồ chuỗi thời gian chính xác.
• Ta xây dựng một theme tùy chỉnh. Hàm theme_pro định nghĩa phông chữ, màu sắc và bố cục. Ta dùng theme_set để áp dụng nó làm mặc định cho mọi biểu đồ.

Phân tích Thống kê và Trực quan hóa Dữ liệu từ BCTC PNJ

Phân tích Tăng trưởng & Quy mô

Doanh thu & Tăng trưởng YoY

if (all(c("revenue","revenue_yoy","Date") %in% names(df_plot))) {

  # Tính hệ số quy đổi cho trục phụ
  s <- safe_coef(df_plot$revenue, df_plot$revenue_yoy)

  # Bảng màu
  col_rev        <- "#2563EB"  # xanh royal cho cột
  col_yoy        <- "#F9F316"  # đường YoY
  col_yoy_trend  <- "#FB923C"  # xu hướng
  zero_line      <- 0 * s

  # (MỚI) Độ rộng cột ~ 80% khoảng cách giữa các mốc thời gian
  bw <- max(10, as.numeric(median(diff(sort(unique(df_plot$Date))))) * 0.8)

  # Điểm cuối để gắn nhãn % YoY
  last_df <- df_plot %>%
    dplyr::filter(!is.na(revenue_yoy), !is.na(Date)) %>%
    dplyr::slice_tail(n = 1) %>%
    dplyr::mutate(
      y_scaled  = revenue_yoy * s,
      yoy_label = scales::percent(revenue_yoy, accuracy = 0.1)
    )

  p1 <- ggplot(df_plot, aes(x = Date)) +
    # --- CỘT DOANH THU (to hơn) ---
    geom_col(aes(y = revenue, fill = "Revenue"),
             width = bw, alpha = 0.95, color = NA) +
    # --- ĐƯỜNG 0% của YoY ---
    geom_hline(yintercept = zero_line, linetype = "dotted", color = "#9CA3AF") +

    # --- Hiệu ứng glow cho YoY ---
    geom_line(aes(y = revenue_yoy * s),
              linewidth = 2.2, color = col_yoy, alpha = 0.12, 
              show.legend = FALSE) +

    # --- Đường & điểm YoY ---
    geom_line(aes(y = revenue_yoy * s, color = "YoY"),
              linewidth = 1.25, lineend = "round") +
    geom_point(aes(y = revenue_yoy * s, color = "YoY"), size = 2.1, stroke = 0.2) +

    # --- Xu hướng LOESS của YoY ---
    geom_smooth(aes(y = revenue_yoy * s, color = "YoY (LOESS)"),
                method = "loess", se = FALSE, linewidth = 0.9, span = 0.5) +

    # --- Nhãn % tại điểm cuối ---
    geom_text(data = last_df,
              aes(x = Date, y = y_scaled, label = yoy_label),
              vjust = -0.8, size = 3.3, color = col_yoy) +

    # --- Thang & trục ---
    scale_y_continuous(
      name = "Doanh thu (Ngàn Tỷ)", labels = num_fmt,
      expand = expansion(mult = c(0.02, 0.1)),
      sec.axis = sec_axis(~ . / s, name = "Tăng trưởng YoY", labels = pct_fmt)
    ) +
    scale_fill_manual(values = c("Revenue" = col_rev)) +
    scale_color_manual(values = c("YoY" = col_yoy, "YoY (LOESS)" = col_yoy_trend)) +

    labs(
      title = "Doanh thu & Tăng trưởng YoY",
      x = NULL, caption = cap_default
    ) +
    guides(fill = guide_legend(order = 1), color = guide_legend(order = 2)) +
    theme_pro() +
    theme(
      legend.box = "horizontal",
      plot.title.position = "plot"
    )

  print(maybe_plotly(add_xdate(p1)))
  maybe_save(add_xdate(p1), "01_revenue_yoy_dual_axis.png") }

Phân tích Code:

• Dùng geom_col để vẽ các cột doanh thu trên trục chính bên trái.
• Dùng geom_line và geom_smooth để vẽ các đường tăng trưởng trên trục phụ bên phải.
• Kỹ thuật trục phụ hoạt động bằng cách nhân dữ liệu trục phụ với một hệ số s, sau đó dùng sec_axis(~ . / s) để hiển thị lại giá trị gốc.
• Dùng geom_smooth với phương pháp loess để làm nổi bật xu hướng dài hạn, loại bỏ các biến động nhiễu trong ngắn hạn.

Kết quả:

• 2017-2020: Doanh thu tăng trưởng ổn định.
• 2021: Tăng trưởng chững lại.
• 2022: Doanh thu và tốc độ tăng trưởng bùng nổ mạnh mẽ.
• 2023-2025: Tốc độ tăng trưởng quay về mức bình thường khi quy mô công ty đã lớn hơn.

Phân rã Tăng trưởng

cols <- c("revenue_yoy","gp_yoy","ebit_yoy","ni_yoy")
avail <- intersect(cols, names(df_plot))

if (length(avail) >= 2) {
  df_long <- df_plot |>
    dplyr::select(Date, dplyr::all_of(avail)) |>
    tidyr::pivot_longer(-Date, names_to = "metric", values_to = "value") |>
    dplyr::mutate(
      metric_lab = dplyr::recode(
        metric,
        "revenue_yoy" = "Doanh thu (YoY)",
        "gp_yoy"      = "Lợi nhuận gộp (YoY)",
        "ebit_yoy"    = "EBIT (YoY)",
        "ni_yoy"      = "Lợi nhuận ròng (YoY)"
      )
    )

  # Màu sắc đậm nhưng dịu, phân biệt tốt
  pal_yoy <- c(
    "Doanh thu (YoY)"      = "#2563EB",
    "Lợi nhuận gộp (YoY)"  = "#10B981",
    "EBIT (YoY)"           = "#F59E0B",
    "Lợi nhuận ròng (YoY)" = "#EF4444"
  )

  # Nhãn ở điểm gần nhất hiện có
  last_pts <- df_long |>
    dplyr::filter(!is.na(value)) |>
    dplyr::group_by(metric_lab) |>
    dplyr::slice_max(order_by = Date, n = 1, with_ties = FALSE) |>
    dplyr::ungroup() |>
    dplyr::mutate(lbl = scales::percent(value, accuracy = 0.1))

  has_repel <- requireNamespace("ggrepel", quietly = TRUE)

  p2 <- ggplot(df_long, aes(Date, value, color = metric_lab, group = metric_lab)) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "#9CA3AF") +
    geom_line(linewidth = 1.05, lineend = "round") +
    geom_point(size = 1.9, stroke = 0.4, shape = 21, fill = "white", show.legend = FALSE) +
    {
      if (has_repel) {
        ggrepel::geom_text_repel(
          data = last_pts, aes(label = lbl),
          size = 3.2, segment.color = "#D1D5DB",
          box.padding = 0.25, point.padding = 0.2,
          max.overlaps = Inf, seed = 1
        )
      } else {
        geom_text(data = last_pts, aes(label = lbl), vjust = -0.6, size = 3.2)
      }
    } +
    scale_color_manual(values = pal_yoy) +
    scale_y_continuous(labels = pct_fmt, breaks = scales::breaks_pretty(n = 4)) +
    labs(
      title = "Phân rã tăng trưởng YoY (tách theo chỉ tiêu)",
      subtitle = "Mỗi ô một chỉ tiêu, trục Y riêng giúp nhìn xu hướng rõ ràng hơn.",
      x = NULL, y = "Tăng trưởng YoY", caption = cap_default
    ) +
    theme_pro() +
    theme(legend.position = "none") +
    facet_wrap(~ metric_lab, ncol = 2, scales = "free_y")

  print(maybe_plotly(add_xdate(p2)))
  maybe_save(add_xdate(p2), "02_growth_decomposition_facets.png")
}

Phân tích Code:

• Dùng pivot_longer để chuyển dữ liệu từ dạng rộng sang dạng dài. ggplot2 yêu cầu định dạng này để tạo biểu đồ.
• Dùng facet_wrap để tạo một biểu đồ nhỏ riêng cho mỗi chỉ tiêu.
• Dùng scales = "free_y". Đây là tham số quan trọng nhất. Nó cho phép mỗi biểu đồ nhỏ có trục tung riêng, giúp các biến động nhỏ không bị ẩn đi.

Kết quả:

• Cả bốn chỉ tiêu đều có chung một mẫu hình tăng trưởng, với một đỉnh đột biến vào năm 2022.
• Tuy nhiên, biên độ tăng trưởng của EBIT và Lợi nhuận ròng thấp hơn và biến động mạnh hơn so với Doanh thu.
• Điều này cho thấy chi phí hoạt động đã tăng. Tăng trưởng doanh thu không được chuyển đổi hiệu quả thành lợi nhuận ròng.

Cấu trúc Tài sản

if (all(c("ta","ca","Date") %in% names(df_plot))) {

  # Chuẩn bị dữ liệu + đặt tên hiển thị đẹp
  df_area <- df_plot |>
    dplyr::transmute(
      Date,
      `Tài sản ngắn hạn (CA)`      = ca,
      `Tài sản dài hạn (TA-CA)`    = ta - ca
    ) |>
    tidyr::pivot_longer(-Date, names_to = "Khoản mục", values_to = "Giá trị") |>
    dplyr::mutate(`Khoản mục` = factor(
      `Khoản mục`,
      levels = c("Tài sản ngắn hạn (CA)", "Tài sản dài hạn (TA-CA)")
    ))

  # Nhãn % tại quý gần nhất (đặt giữa mỗi lớp)
  labels_last <- df_area |>
    dplyr::filter(Date == max(Date, na.rm = TRUE)) |>
    dplyr::arrange(`Khoản mục`) |>
    dplyr::mutate(
      Tổng   = sum(`Giá trị`, na.rm = TRUE),
      Tỷ_trọng = dplyr::if_else(is.finite(Tổng) & Tổng != 0, `Giá trị`/Tổng, NA_real_),
      y_mid  = cumsum(`Giá trị`) - `Giá trị`/2
    ) |>
    dplyr::filter(!is.na(Tỷ_trọng)) |>
    dplyr::mutate(lbl = scales::percent(Tỷ_trọng, accuracy = 0.1))

  # Bảng màu: teal + violet (tươi, dễ phân biệt)
  pal_asset <- c(
    "Tài sản ngắn hạn (CA)"   = "#06B6D4",  # teal
    "Tài sản dài hạn (TA-CA)" = "#8B5CF6"   # violet
  )

  p3 <- ggplot(df_area, aes(Date, `Giá trị`, fill = `Khoản mục`, group = `Khoản mục`)) +
    # nền nhẹ giúp nổi khối
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = Inf,
             fill = "#F8FAFC", alpha = 0.6) +
    # stacked area
    geom_area(alpha = 0.92, color = NA) +
    # nhãn % trong ô ở quý gần nhất (màu chữ trắng, không viền)
    geom_label(
      data = labels_last,
      aes(x = max(df_area$Date, na.rm = TRUE), y = y_mid, label = lbl, fill = `Khoản mục`),
      color = "white", size = 3.2, label.size = 0, label.padding = grid::unit(0.12, "lines"),
      inherit.aes = FALSE
    ) +
    scale_fill_manual(values = pal_asset) +
    scale_y_continuous(labels = num_fmt, expand = expansion(mult = c(0.02, 0.12))) +
    labs(
      title = "Cấu trúc Tài sản",
      subtitle = "Stacked area: CA vs. (TA−CA). Nhãn hiển thị tỷ trọng tại quý gần nhất.",
      x = NULL, y = "Ngàn Tỷ", caption = cap_default
    ) +
    theme_pro() +
    theme(
      legend.position = "top",
      legend.box = "horizontal",
      plot.title.position = "plot"
    )

  print(maybe_plotly(add_xdate(p3)))
  maybe_save(add_xdate(p3), "03_asset_structure.png")
}

Phân tích Code:

• Tính tài sản dài hạn bằng ta - ca, sau đó dùng pivot_longer để chuyển dữ liệu sang dạng dài. Định dạng này là yêu cầu để ggplot2 vẽ biểu đồ xếp chồng.
• Tính toán vị trí của các nhãn phần trăm. Công thức này đặt chúng vào giữa mỗi lớp diện tích tại điểm dữ liệu cuối cùng.
• Dùng geom_area để vẽ biểu đồ và geom_label để thêm các nhãn đã tính toán. Biểu đồ này cho thấy:

Kết quả:

• Tổng tài sản của công ty tăng trưởng đều đặn.
• Tài sản ngắn hạn chiếm phần lớn, khoảng 91.6% tại điểm cuối. Cấu trúc này ổn định qua thời gian.
• Có thể kết luận đây là một doanh nghiệp thâm dụng vốn lưu động, ví dụ như ngành bán lẻ, phụ thuộc vào việc quản lý hàng tồn kho và các khoản phải thu.

Nhóm chỉ số liên quan tới ROE vs ROA

ROE vs ROA

if (all(c("roe","roa","Date") %in% names(df_plot))) {

  # 1) Dữ liệu gốc
  df_lr <- df_plot |>
    dplyr::arrange(Date) |>
    dplyr::transmute(Date, ROE = roe, ROA = roa)

  # 2) Chèn điểm giao nhau (gap = 0) giữa hai quý nếu có đổi dấu
  add_crossings <- function(d){
    d <- d[order(d$Date), ]
    out <- d[1, , drop = FALSE]
    if (nrow(d) >= 2) {
      for (i in seq_len(nrow(d) - 1)) {
        r1 <- d[i, ]; r2 <- d[i + 1, ]
        g1 <- r1$ROE - r1$ROA; g2 <- r2$ROE - r2$ROA
        if (is.finite(g1) && is.finite(g2) && g1 * g2 < 0) {
          t <- abs(g1) / (abs(g1) + abs(g2))                    # t ∈ (0,1)
          cross <- r1
          cross$Date <- r1$Date + as.difftime(as.numeric(r2$Date - r1$Date) * t, units = "days")
          cross$ROE  <- r1$ROE + (r2$ROE - r1$ROE) * t
          cross$ROA  <- r1$ROA + (r2$ROA - r1$ROA) * t
          out <- rbind(out, cross, r2)
        } else {
          out <- rbind(out, r2)
        }
      }
    }
    out
  }

  rib <- add_crossings(df_lr) |>
    dplyr::filter(!is.na(ROE), !is.na(ROA)) |>
    dplyr::mutate(
      gap  = ROE - ROA,
      ymin = pmin(ROE, ROA),
      ymax = pmax(ROE, ROA),
      sign = ifelse(gap >= 0, "Đòn bẩy dương (ROE > ROA)", "Đòn bẩy âm (ROE < ROA)")
    ) |>
    dplyr::arrange(Date) |>
    dplyr::mutate(                           # chia polygon theo đoạn cùng dấu
      sg  = ifelse(gap >= 0, 1L, -1L),
      grp = cumsum(dplyr::coalesce(sg != dplyr::lag(sg), FALSE))
    )

  # Nhãn Δ tại quý gần nhất có đủ cặp
  last_pair <- rib |>
    dplyr::slice_tail(n = 1) |>
    dplyr::mutate(
      y_mid     = (ROE + ROA) / 2,
      label_gap = paste0(ifelse(gap >= 0, "+", ""),
                         scales::number(gap * 100, accuracy = 0.1, decimal.mark = ","),
                         " điểm %"),
      col_gap   = ifelse(gap >= 0, "#10B981", "#EF4444")
    )

  col_line <- c("ROE" = "#F59E0B", "ROA" = "#2563EB")
  col_fill <- c("Đòn bẩy dương (ROE > ROA)" = "#10B981",
                "Đòn bẩy âm (ROE < ROA)"   = "#EF4444")

  # Zoom trục Y ~30% (thu hẹp phạm vi để phóng lớn)
  rng <- range(c(df_lr$ROE, df_lr$ROA), na.rm = TRUE)
  mid <- mean(rng); half <- diff(rng) / 2
  ylim_zoom <- c(mid - half * 0.70, mid + half * 0.70)

  p6 <- ggplot() +
    annotate("rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=Inf, fill="#F8FAFC", alpha=.6) +
    geom_ribbon(data = rib,
                aes(Date, ymin = ymin, ymax = ymax, fill = sign, group = grp),
                alpha = .25, show.legend = TRUE) +
    geom_line(data = df_lr, aes(Date, ROE), linewidth = 2, color = col_line["ROE"], alpha = .12) +
    geom_line(data = df_lr, aes(Date, ROE, color = "ROE"), linewidth = 1.25, lineend = "round") +
    geom_point(data = df_lr, aes(Date, ROE, color = "ROE"), size = 2, stroke = .3) +
    geom_line(data = df_lr, aes(Date, ROA), linewidth = 2, color = col_line["ROA"], alpha = .12) +
    geom_line(data = df_lr, aes(Date, ROA, color = "ROA"), linewidth = 1.25, lineend = "round") +
    geom_point(data = df_lr, aes(Date, ROA, color = "ROA"), size = 2, stroke = .3) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "#9CA3AF") +
    geom_label(
      data = last_pair,
      aes(Date, y_mid, label = label_gap),
      inherit.aes = FALSE, size = 2, label.size = 0, fill = "white",
      color = last_pair$col_gap[1], fontface = "bold"
    ) +
    scale_y_continuous(labels = pct_fmt, expand = expansion(mult = c(.06, .18))) +
    scale_color_manual(values = col_line, name = NULL) +
    scale_fill_manual(values = col_fill,  name = NULL) +
    scale_x_date(expand = expansion(add = c(30, 150)), date_breaks = "1 year", date_labels = "%Y") +
    coord_cartesian(ylim = ylim_zoom, clip = "off") +
    labs(
      title = "ROE vs ROA",
      x = NULL, y = "Tỷ lệ", caption = cap_default
    ) +
    theme_pro() +
    theme(legend.position = "top", legend.box = "horizontal",
          plot.title.position = "plot", axis.text.x = element_text(angle = 45, hjust = 1))

  print(maybe_plotly(p6))
  maybe_save(p6, "06_roe_roa.png")
}

Phân tích Code:

• Dùng geom_ribbon để tô màu vùng giữa hai đường ROE và ROA.
• Dùng hàm add_crossings. Hàm này chèn thêm các điểm dữ liệu tại đúng vị trí hai đường cắt nhau. Việc này đảm bảo vùng tô màu chính xác.
• Tạo một biến nhóm grp. Biến này giúp geom_ribbon vẽ các vùng màu riêng biệt cho từng giai đoạn đòn bẩy dương và âm.
• Các lớp khác như geom_line và geom_point được vẽ lên trên để làm rõ xu hướng của từng chỉ số.

Kết quả:

• Vùng màu xanh lá: ROE > ROA. Việc dùng nợ đã khuếch đại lợi nhuận cho cổ đông (đòn bẩy dương).
• Vùng màu đỏ: ROE < ROA. Chi phí lãi vay cao hơn lợi nhuận tài sản tạo ra, làm giảm lợi nhuận của cổ đông (đòn bẩy âm).
• Phân tích:
  • Công ty chủ yếu sử dụng đòn bẩy hiệu quả.
  • Vào cuối năm 2021, cả hai chỉ số đều âm. Quan trọng là ROE giảm sâu hơn ROA, cho thấy gánh nặng lãi vay đã làm tình hình tồi tệ hơn.
  • Khoảng cách giữa ROE và ROA đang có xu hướng thu hẹp.

DuPont 3 thành tố

if (all(c("roe","net_margin","asset_turnover","financial_leverage","Date") %in% names(df_plot))) {

  # Hệ số quy đổi an toàn cho trục phụ
  s <- suppressWarnings(
    max(df_plot$roe, na.rm = TRUE) /
      max(c(df_plot$net_margin, df_plot$asset_turnover, df_plot$financial_leverage), na.rm = TRUE)
  )
  if (!is.finite(s) || s <= 0) s <- 1

  uq  <- sort(unique(df_plot$Date))
  step_days <- if (length(uq) > 1) median(diff(uq)) else 90
  bar_w <- as.numeric(step_days) * 0.75  # ~75% khoảng cách giữa 2 quý

  # Long-format cho các thành tố
  comps <- df_plot %>%
    dplyr::select(Date, net_margin, asset_turnover, financial_leverage) %>%
    tidyr::pivot_longer(-Date, names_to = "component", values_to = "value") %>%
    dplyr::mutate(
      component_lab = dplyr::recode(
        component,
        net_margin        = "Net margin",
        asset_turnover    = "Asset turnover",
        financial_leverage= "Financial leverage"
      ),
      value_scaled = value * s
    )

  # Nhãn ở quý gần nhất
  last_labels <- comps %>%
    dplyr::filter(!is.na(value)) %>%
    dplyr::group_by(component_lab) %>%
    dplyr::slice_max(order_by = Date, n = 1, with_ties = FALSE) %>%
    dplyr::ungroup() %>%
    dplyr::mutate(
      label = dplyr::case_when(
        component_lab == "Net margin"        ~ scales::percent(value, accuracy = 0.1),
        component_lab == "Asset turnover"    ~ scales::number(value, accuracy = 0.01, suffix = "x",
                                                             big.mark = ".", decimal.mark = ","),
        component_lab == "Financial leverage"~ scales::number(value, accuracy = 0.01, suffix = "x",
                                                             big.mark = ".", decimal.mark = ","),
        TRUE ~ as.character(value)
      )
    )

  # Bảng màu
  col_bar  <- "#22C55E"  # ROE (emerald)
  col_line <- c(
    "Net margin"         = "#0EA5E9",
    "Asset turnover"     = "#F59E0B",
    "Financial leverage" = "#8B5CF6"
  )

  has_repel <- requireNamespace("ggrepel", quietly = TRUE)

  p7 <- ggplot() +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = Inf,
             fill = "#F8FAFC", alpha = 0.6) +

    # CỘT ROE
    geom_col(data = dplyr::filter(df_plot, !is.na(roe)),
             aes(Date, roe), width = bar_w,
             fill = col_bar, alpha = 0.92, color = NA, na.rm = TRUE) +

    # GLOW + ĐƯỜNG CÁC THÀNH TỐ (đã scale sang trục trái)
    geom_line(data = comps, aes(Date, value_scaled, color = component_lab),
              linewidth = 2.3, alpha = 0.10, show.legend = FALSE) +
    geom_line(data = comps, aes(Date, value_scaled, color = component_lab),
              linewidth = 1.25, lineend = "round") +
    geom_point(data = comps, aes(Date, value_scaled, color = component_lab),
               size = 1.9, stroke = 0.3, show.legend = FALSE) +

    # NHÃN Ở QUÝ GẦN NHẤT
    {
      if (has_repel) {
        ggrepel::geom_text_repel(
          data = last_labels,
          aes(Date, value_scaled, label = label, color = component_lab),
          size = 3.1, segment.color = "#D1D5DB",
          box.padding = 0.25, point.padding = 0.2,
          max.overlaps = Inf, seed = 1, show.legend = FALSE
        )
      } else {
        geom_text(
          data = last_labels,
          aes(Date, value_scaled, label = label, color = component_lab),
          vjust = -0.6, size = 3.1, show.legend = FALSE
        )
      }
    } +

    geom_hline(yintercept = 0, linetype = "dotted", color = "#9CA3AF") +

    scale_y_continuous(
      name = "ROE",
      labels = pct_fmt,
      expand = expansion(mult = c(0.05, 0.18)),
      sec.axis = sec_axis(~ . / s, name = "Thành tố DuPont")
    ) +
    scale_color_manual(values = col_line, name = NULL) +
    labs(
      title = "Phân rã DuPont (3 thành tố)",
      x = NULL, caption = cap_default
    ) +
    theme_pro() +
    theme(
      legend.position = "top",
      legend.box = "horizontal",
      plot.title.position = "plot"
    )

  print(maybe_plotly(add_xdate(p7))) }

Phân tích Code:

• Dùng pivot_longer để tái cấu trúc dữ liệu, chuẩn bị các thành tố DuPont cho việc vẽ biểu đồ.
• Vẽ ROE dưới dạng các cột trên trục tung chính bên trái.
• Vẽ ba thành tố DuPont (Biên lợi nhuận ròng, Vòng quay tài sản, Đòn bẩy tài chính) dưới dạng các đường trên trục tung phụ bên phải. Kỹ thuật này cho phép Ta so sánh các biến có thang đo khác nhau trên cùng một biểu đồ.

Kết quả:

• Các cột ROE có xu hướng biến động cùng chiều với đường Vòng quay tài sản.
• Biên lợi nhuận ròng và Đòn bẩy tài chính tương đối ổn định.
• Điều này cho thấy sự thay đổi trong ROE của công ty chủ yếu được quyết định bởi hiệu quả sử dụng tài sản, không phải bởi khả năng sinh lời hay việc sử dụng nợ.

ROA Decomposition: Scatter (AT vs Net margin)

if (all(c("asset_turnover","net_margin","Date") %in% names(df_plot))) {

  # Tính median để vẽ trục tham chiếu chiến lược
  med_x <- median(df_plot$asset_turnover, na.rm = TRUE)
  med_y <- median(df_plot$net_margin,     na.rm = TRUE)

  # Breaks/labels cho thanh màu thời gian
  t_min <- suppressWarnings(min(df_plot$Date, na.rm = TRUE))
  t_max <- suppressWarnings(max(df_plot$Date, na.rm = TRUE))
  yr_seq <- try(seq(as.Date(cut(t_min, "year")), as.Date(cut(t_max, "year")), by = "2 years"), silent = TRUE)
  if (inherits(yr_seq, "try-error") || length(yr_seq) < 2) {
    yr_seq <- unique(as.Date(c(t_min, t_max)))
  }

  # Bảng màu gradient (đẹp, dễ đọc)
  time_cols <- c("#2563EB", "#8B5CF6", "#F59E0B", "#EF4444")  # blue → violet → amber → coral

  p8 <- ggplot(df_plot, aes(x = asset_turnover, y = net_margin)) +
    # Nền cho vùng Net margin âm
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0,
             fill = "#FEE2E2", alpha = 0.10) +

    # Halo trắng (lớp dưới) để điểm nổi bật trên nền
    geom_point(color = "white", size = 3.6, alpha = 0.8, na.rm = TRUE) +
    # Điểm màu theo thời gian (lớp trên)
    geom_point(aes(color = as.numeric(Date)), size = 2.2, alpha = 0.95, na.rm = TRUE) +

    # Đường LOESS mượt (xu hướng phi tuyến)
    geom_smooth(aes(color = NULL), method = "loess", se = FALSE, linewidth = 1.1, span = 0.6, color = "#111827") +
    # Đường hồi quy tuyến tính tổng quát (trade-off margin vs turnover)
    geom_smooth(aes(color = NULL), method = "lm", se = FALSE, linewidth = 0.9, linetype = "dashed", color = "#6B7280") +

    # Đường tham chiếu theo median (phân mảnh chiến lược)
    geom_vline(xintercept = med_x, linetype = "dotted", color = "#9CA3AF") +
    geom_hline(yintercept = med_y, linetype = "dotted", color = "#9CA3AF") +

    # Tỷ lệ & thang màu
    scale_x_continuous(labels = scales::number_format(accuracy = 0.01, suffix = "x"),
                       breaks = scales::breaks_pretty(n = 5), expand = expansion(mult = c(0.05, 0.08))) +
    scale_y_continuous(labels = pct_fmt,
                       breaks = scales::breaks_pretty(n = 5), expand = expansion(mult = c(0.05, 0.10))) +
    scale_color_gradientn(
      colours = time_cols,
      breaks  = as.numeric(yr_seq),
      labels  = format(yr_seq, "%Y"),
      name    = "Thời gian"
    ) +

    labs(
      title = "Phân rã ROA: Asset Turnover vs Net Margin",
      x = "Asset turnover (lần)",
      y = "Net margin",
      caption = cap_default
    ) +
    theme_pro() +
    theme(
      legend.position = "right",
      legend.box = "vertical",
      plot.title.position = "plot"
    ) +
    coord_cartesian(clip = "off")

  # Không dùng add_xdate vì trục X không phải Date
  print(maybe_plotly(p8))}

Phân tích Code:

• Tạo một biểu đồ phân tán (geom_point) với Asset turnover trên trục hoành và Net margin trên trục tung.
• Mã hóa thời gian vào màu sắc của các điểm. Màu xanh/tím là quá khứ, màu cam/đỏ là hiện tại. Điều này biến một biểu đồ tĩnh thành một chuỗi thời gian.
• Vẽ hai đường xu hướng: một đường hồi quy tuyến tính (lm) thể hiện mối quan hệ tổng thể, và một đường loess thể hiện xu hướng phi tuyến cục bộ.
• Thêm các đường tham chiếu tại giá trị trung vị của mỗi trục, chia biểu đồ thành bốn góc phần tư chiến lược.

Kết quả:

• Mỗi điểm là một quý, màu sắc cho thấy dòng thời gian.
• Các điểm dữ liệu di chuyển từ trái sang phải theo thời gian, cho thấy Vòng quay tài sản ngày càng cải thiện.
• Mối quan hệ tổng thể (đường gạch nối) là dương: khi công ty bán hàng hiệu quả hơn (turnover tăng), biên lợi nhuận cũng có xu hướng tăng.
• Tuy nhiên, đường xu hướng cục bộ (đường đậm) cho thấy một sự đánh đổi: ở mức Vòng quay tài sản rất cao (trên 0.6x), Biên lợi nhuận bắt đầu đi ngang hoặc giảm nhẹ. Điều này có thể do công ty phải giảm giá hoặc tăng chi phí để đạt được doanh số cao hơn.
• Hầu hết các điểm gần đây (màu cam/đỏ) đều nằm ở góc phần tư trên bên phải, thể hiện một giai đoạn hoạt động hiệu quả cả về biên lợi nhuận và vòng quay tài sản.

ROIC vs EBIT margin

df_plot2 <- df_plot
if (!("roic" %in% names(df_plot2)) || all(is.na(df_plot2$roic))) {
  n <- nrow(df_plot2)
  base_ic <- if ("ta" %in% names(df_plot2)) df_plot2$ta else rep(NA_real_, n)
  cash    <- if ("cash_eq" %in% names(df_plot2)) dplyr::coalesce(df_plot2$cash_eq, 0) else 0
  nibcl   <- if (all(c("cl","ap") %in% names(df_plot2))) dplyr::coalesce(df_plot2$cl - df_plot2$ap, 0) else 0
  ic_proxy <- if ("ic" %in% names(df_plot2)) df_plot2$ic else base_ic - cash - nibcl
  tax_eff <- if ("tax_rate" %in% names(df_plot2)) pmin(pmax(df_plot2$tax_rate, 0), 1) else rep(0.20, n)
  avg_ic  <- avg_lag1(ic_proxy)
  roic_proxy <- ifelse(is.finite(avg_ic) & avg_ic != 0 & "ebit" %in% names(df_plot2),
                       df_plot2$ebit * (1 - tax_eff) / avg_ic, NA_real_)
  df_plot2$roic <- if ("roic" %in% names(df_plot2)) dplyr::coalesce(df_plot2$roic, roic_proxy) else roic_proxy
}

# --- Vẽ ---
cols <- intersect(c("roic","operating_margin"), names(df_plot2))
if (length(cols) >= 1 && "Date" %in% names(df_plot2)) {

  rename_map <- c(roic = "ROIC", operating_margin = "EBIT margin")
  df_sel <- df_plot2 %>%
    dplyr::select(Date, dplyr::all_of(cols)) %>%
    dplyr::rename_with(~ rename_map[.x] %||% .x, .cols = -Date) %>%
    tidyr::pivot_longer(-Date, names_to = "metric", values_to = "value")

  # Nhãn quý gần nhất (dạng %)
  last_labels <- df_sel %>%
    dplyr::filter(!is.na(value)) %>%
    dplyr::group_by(metric) %>%
    dplyr::slice_max(order_by = Date, n = 1, with_ties = FALSE) %>%
    dplyr::ungroup() %>%
    dplyr::mutate(lbl = scales::percent(value, accuracy = 0.1))

  pal_line <- c("ROIC" = "#10B981", "EBIT margin" = "#2563EB", "WACC" = "#EF4444")
  has_repel <- requireNamespace("ggrepel", quietly = TRUE)

  has_wacc <- "wacc" %in% names(df_plot2)
  rib <- if (has_wacc && "roic" %in% names(df_plot2)) {
    df_plot2 %>% dplyr::select(Date, roic, wacc) %>%
      dplyr::filter(!is.na(roic), !is.na(wacc)) %>%
      dplyr::mutate(ymin = pmin(roic, wacc), ymax = pmax(roic, wacc))
  } else NULL

  p9 <- ggplot() +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = Inf,
             fill = "#F8FAFC", alpha = 0.6) +
    { if (!is.null(rib)) list(
        geom_ribbon(data = rib, aes(Date, ymin = ymin, ymax = ymax, fill = "Spread"),
                    alpha = 0.15, show.legend = FALSE),
        scale_fill_manual(values = c("Spread" = "#10B981"))
      ) else NULL } +
    geom_line(data = df_sel, aes(Date, value, color = metric),
              linewidth = 2.2, alpha = 0.10, show.legend = FALSE) +
    geom_line(data = df_sel, aes(Date, value, color = metric),
              linewidth = 1.25, lineend = "round") +
    geom_point(data = df_sel, aes(Date, value, color = metric),
               size = 1.9, stroke = 0.3, show.legend = FALSE) +
    { if (has_wacc) list(
        geom_line(data = df_plot2, aes(Date, wacc, color = "WACC"),
                  linewidth = 1.1, linetype = "longdash"),
        geom_point(data = df_plot2, aes(Date, wacc, color = "WACC"),
                   size = 1.6, alpha = 0.9, show.legend = FALSE)
      ) else NULL } +
    geom_hline(yintercept = 0, linetype = "dotted", color = "#9CA3AF") +

    # ===== NHÃN % =====
    {
      if (has_repel) {
        ggrepel::geom_label_repel(
          data = last_labels,
          aes(Date, value, label = lbl, color = metric),
          fill = "white", label.size = 0, fontface = "bold",
          size = 3.1, segment.color = "#D1D5DB",
          box.padding = 0.25, point.padding = 0.2,
          max.overlaps = Inf, seed = 1, show.legend = FALSE
        )
      } else {
        geom_label(
          data = last_labels,
          aes(Date, value, label = lbl, color = metric),
          fill = "white", label.size = 0, fontface = "bold",
          size = 3.1, show.legend = FALSE
        )
      }
    } +

    scale_color_manual(values = pal_line, name = NULL) +
    scale_y_continuous(labels = pct_fmt,
                       breaks = scales::breaks_pretty(n = 5),
                       expand = expansion(mult = c(0.05, 0.16))) +
    labs(
      title = "ROIC & EBIT margin",
      x = NULL, y = "Tỷ lệ", caption = cap_default
    ) +
    theme_pro() +
    theme(legend.position = "top", legend.box = "horizontal", plot.title.position = "plot")

  print(maybe_plotly(add_xdate(p9)))}

Phân tích Code:

• Tự động ước tính ROIC nếu dữ liệu gốc bị thiếu. Điều này làm cho quy trình của Ta linh hoạt hơn.
• Dùng geom_line và geom_point để vẽ xu hướng của hai chỉ số này theo thời gian.

Kết quả:

• Đường EBIT Margin nằm trên đường ROIC. Khoảng cách này cho thấy Vòng quay vốn đầu tư của PNJ nhỏ hơn 1.
• Ta cần hơn 1 đồng vốn đầu tư để tạo ra 1 đồng doanh thu. Đây là đặc điểm của ngành bán lẻ hoặc phân phối.
• Cả hai chỉ số đều âm vào cuối năm 2021. Điều này xác nhận công ty đã lỗ từ hoạt động kinh doanh cốt lõi trong giai đoạn đó.

Hiệu quả & Vốn lưu động

Phân rã CCC

#  cột
need <- c("Date","ccc","dso","dio","dpo")
if (all(need %in% names(df_plot))) {

  # ===== Chuẩn bị dữ liệu + MA 4 quý =====
  df_ccc <- df_plot |> dplyr::arrange(Date) |>
    dplyr::mutate(ccc_ma4 = zoo::rollmean(ccc, k = 4, fill = NA, align = "right"))
  df_dso <- df_plot |> dplyr::arrange(Date) |>
    dplyr::transmute(Date, DSO = dso, MA4 = zoo::rollmean(dso, k = 4, fill = NA, align = "right"))
  df_dio <- df_plot |> dplyr::arrange(Date) |>
    dplyr::transmute(Date, DIO = dio, MA4 = zoo::rollmean(dio, k = 4, fill = NA, align = "right"))
  df_dpo <- df_plot |> dplyr::arrange(Date) |>
    dplyr::transmute(Date, DPO = dpo, MA4 = zoo::rollmean(dpo, k = 4, fill = NA, align = "right"))

  # ===== Màu sắc nhất quán =====
  col_ccc <- "#2563EB"   # royal blue
  col_ma  <- "#F59E0B"   # amber
  col_dso <- "#0EA5E9"   # sky
  col_dio <- "#8B5CF6"   # violet
  col_dpo <- "#EF4444"   # coral

  # ===== Helper panel gọn =====
  base_panel <- list(
    theme_pro(),
    theme(
      legend.position = "none",
      plot.title      = element_text(size = 11, face = "bold"),
      axis.text.x     = element_text(size = 8, angle = 45, hjust = 1),
      axis.text.y     = element_text(size = 8)
    ),
    scale_y_continuous(
      labels = scales::label_number(accuracy = 1, big.mark = ".", decimal.mark = ","),
      breaks = scales::breaks_pretty(n = 4),
      expand = expansion(mult = c(0.05, 0.10))
    )
  )

  # ===== Panel A: CCC tổng =====
  p_ccc <- ggplot(df_ccc, aes(Date, ccc)) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf,
             fill = "#F8FAFC", alpha = 0.65) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "#9CA3AF") +
    geom_line(linewidth = 2.2, color = col_ccc, alpha = 0.12) +           # halo
    geom_line(color = col_ccc, linewidth = 1.15, lineend = "round") +
    geom_line(aes(y = ccc_ma4), color = col_ma, linewidth = 0.95, linetype = "longdash", na.rm = TRUE) +
    labs(title = "CCC", x = NULL, y = "Ngày", caption = cap_default) +
    base_panel

  # ===== Panel B: DSO =====
  p_dso <- ggplot(df_dso, aes(Date, DSO)) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf,
             fill = "#F8FAFC", alpha = 0.65) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "#9CA3AF") +
    geom_line(linewidth = 2.2, color = col_dso, alpha = 0.12) +
    geom_line(color = col_dso, linewidth = 1.10) +
    geom_line(aes(y = MA4), color = col_ma, linewidth = 0.9, linetype = "longdash", na.rm = TRUE) +
    labs(title = "DSO (Kỳ thu tiền)", x = NULL, y = "Ngày") +
    base_panel

  # ===== Panel C: DIO =====
  p_dio <- ggplot(df_dio, aes(Date, DIO)) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf,
             fill = "#F8FAFC", alpha = 0.65) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "#9CA3AF") +
    geom_line(linewidth = 2.2, color = col_dio, alpha = 0.12) +
    geom_line(color = col_dio, linewidth = 1.10) +
    geom_line(aes(y = MA4), color = col_ma, linewidth = 0.9, linetype = "longdash", na.rm = TRUE) +
    labs(title = "DIO (Kỳ tồn kho)", x = NULL, y = "Ngày") +
    base_panel

  # ===== Panel D: DPO =====
  p_dpo <- ggplot(df_dpo, aes(Date, DPO)) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf,
             fill = "#F8FAFC", alpha = 0.65) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "#9CA3AF") +
    geom_line(linewidth = 2.2, color = col_dpo, alpha = 0.12) +
    geom_line(color = col_dpo, linewidth = 1.10) +
    geom_line(aes(y = MA4), color = col_ma, linewidth = 0.9, linetype = "longdash", na.rm = TRUE) +
    labs(title = "DPO (Kỳ thanh toán)", x = NULL, y = "Ngày") +
    base_panel

  # ===== Gộp 4 panel thành 1 hình =====
  if (requireNamespace("patchwork", quietly = TRUE)) {
    combo <- (p_ccc + p_dso) / (p_dio + p_dpo) +
      patchwork::plot_annotation(
        title = "Phân rã CCC — Đường nét đứt là MA4",
        theme = theme_pro() + theme(plot.title = element_text(size = 14, face = "bold"))
      )
    print(maybe_plotly(combo))
    maybe_save(combo, "11_ccc_4panel.png", w = 10, h = 7)
  } else if (requireNamespace("gridExtra", quietly = TRUE)) {
    gridExtra::grid.arrange(p_ccc, p_dso, p_dio, p_dpo, ncol = 2,
                            top = grid::textGrob("Phân rã CCC — Đường nét đứt là MA4",
                                                 gp = grid::gpar(fontsize = 14, fontface = "bold")))
  } else {
    # Fallback đơn giản nếu thiếu cả patchwork & gridExtra: in tuần tự
    print(p_ccc); print(p_dso); print(p_dio); print(p_dpo)
  }
}

Phân tích Code:

• Tính toán trung bình trượt 4 kỳ (zoo::rollmean) cho mỗi chỉ số. Đường trung bình trượt giúp làm mịn dữ liệu và làm nổi bật xu hướng dài hạn.
• Tạo ra bốn biểu đồ riêng biệt cho CCC, DSO, DIO và DPO.
• Dùng gói patchwork để ghép bốn biểu đồ này thành một lưới 2x2. patchwork cung cấp một cú pháp đơn giản (+ và /) để sắp xếp các biểu đồ.

Kết quả:

• Biểu đồ cho thấy CCC (ô trên bên trái) tăng đột biến vào cuối năm 2021.
• Nhìn vào các biểu đồ thành phần, Ta có thể xác định nguyên nhân. DIO (ô dưới bên trái) cũng có một đỉnh cực lớn tại cùng thời điểm.
• Điều này cho thấy vấn đề của công ty là do quản lý hàng tồn kho kém hiệu quả, không phải do thu tiền từ khách hàng (DSO).

Tỷ lệ Vònng quay

cols <- c("asset_turnover","inventory_turnover","receivables_turnover")
avail <- intersect(cols, names(df_plot))

if (length(avail) >= 2 && "Date" %in% names(df_plot)) {

  rename_map <- c(
    asset_turnover       = "T.sản",
    inventory_turnover   = "T.kho",
    receivables_turnover = "P.thu"
  )

  df_long <- df_plot |>
    dplyr::select(Date, dplyr::all_of(avail)) |>
    dplyr::rename_with(~ rename_map[.x] %||% .x, .cols = -Date) |>
    tidyr::pivot_longer(-Date, names_to = "Chỉ số", values_to = "Giá trị") |>
    dplyr::arrange(Date)

  df_ma <- df_long |>
    dplyr::group_by(`Chỉ số`) |>
    dplyr::mutate(`MA 4 quý` = zoo::rollmean(`Giá trị`, k = 4, fill = NA, align = "right")) |>
    dplyr::ungroup()

  last_lbls <- df_ma |>
    dplyr::filter(!is.na(`Giá trị`)) |>
    dplyr::group_by(`Chỉ số`) |>
    dplyr::slice_max(order_by = Date, n = 1, with_ties = FALSE) |>
    dplyr::ungroup() |>
    dplyr::mutate(
      x_lbl = as.Date(Date) + 90, 
      y_lbl = `Giá trị`,
      lbl   = paste0(scales::number(`Giá trị`, accuracy = 0.1, big.mark = ".", decimal.mark = ","), "x")
    )

  
  col_map <- c(
    "T.sản" = "#10B981",  # emerald
    "T.kho" = "#8B5CF6",  # violet
    "P.thu KH"= "#0EA5E9"   # sky
  )
  col_ma <- "#F59E0B"

  p12 <- ggplot(df_ma, aes(x = Date)) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf,
             fill = "#F8FAFC", alpha = 0.65) +
    geom_line(aes(y = `Giá trị`, color = `Chỉ số`), linewidth = 2.3, 
              alpha = 0.12, show.legend = FALSE) +
    geom_line(aes(y = `Giá trị`, color = `Chỉ số`), linewidth = 1.15, lineend = "round") +
    geom_line(aes(y = `MA 4 quý`), color = col_ma, linewidth = 0.95, 
    linetype = "longdash", na.rm = TRUE) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "#9CA3AF") +
    geom_segment(data = last_lbls,
                 aes(x = Date, xend = x_lbl, y = y_lbl, yend = y_lbl, color = `Chỉ số`),
                 linewidth = 0.6, alpha = 0.6, show.legend = FALSE) +
    geom_label(data = last_lbls,
               aes(x = x_lbl, y = y_lbl, label = lbl, color = `Chỉ số`),
               fill = "white", label.size = 0, size = 3.1, show.legend = FALSE) +
    scale_color_manual(values = col_map, name = NULL) +
    scale_y_continuous(
      labels = scales::label_number(accuracy = 0.1, big.mark = ".", decimal.mark = ",", suffix = "x"),
      breaks = scales::breaks_pretty(n = 3),
      expand = expansion(mult = c(0.05, 0.12)),
      guide  = guide_axis(check.overlap = TRUE)
    ) +
    scale_x_date(date_labels = "%Y", date_breaks = "1 year",
                 expand = expansion(mult = c(0.02, 0.32))) +
    coord_cartesian(clip = "off") +
    labs(
      title = "Vòng quay",
      subtitle = "Dường nét đứt là MA4",
      x = NULL, y = NULL, caption = cap_default
    ) +
    theme_pro() +
    theme(
      legend.position    = "top",
      strip.placement    = "outside",
      strip.background   = element_rect(fill = "#EEF2FF", color = NA),
      strip.text         = element_text(face = "bold", margin = margin(b = 4, t = 4)),
      axis.text.x        = element_text(angle = 45, hjust = 1),
      axis.text.y        = element_text(margin = margin(r = 8)),
      panel.spacing.y    = grid::unit(16, "pt"),
      plot.margin        = margin(10, 46, 10, 24)
    ) +
    facet_wrap(~ `Chỉ số`, ncol = 1, scales = "free_y", strip.position = "left")

  print(maybe_plotly(p12))
  maybe_save(p12, "12_turnover_ratios.png")
}

Phân tích Code:

• Dùng pivot_longer để tái cấu trúc dữ liệu. ggplot2 yêu cầu định dạng này để tạo biểu đồ.
• Tính trung bình trượt 4 kỳ (zoo::rollmean). Đường trung bình trượt giúp làm mịn dữ liệu và làm nổi bật xu hướng.
• Dùng facet_wrap để tạo một biểu đồ nhỏ cho mỗi chỉ số.
• Dùng scales = "free_y". Tham số này là quan trọng nhất. Nó cho phép mỗi biểu đồ có trục tung riêng vì các chỉ số có thang đo rất khác nhau.

Kết quả:

• Vòng quay các khoản phải thu rất cao, khoảng 100 lần. Điều này cho thấy công ty thu tiền từ khách hàng rất nhanh.
• Vòng quay hàng tồn kho và tài sản thấp hơn nhiều, chỉ khoảng 0.5x đến 1.0x.
• Sự khác biệt lớn về thang đo này là lý do chính để sử dụng các biểu đồ riêng biệt. Nếu vẽ chung trên một trục, các chỉ số vòng quay thấp sẽ trông như một đường thẳng.

Dòng tiền

Waterfall CFO / CFI / CFF

if (all(c("cfo","cfi","cff","Date") %in% names(df_plot))) {

  # Data 
  df_cf <- df_plot |>
    dplyr::transmute(
      Date,
      `CFO (HĐKD)`      = cfo,
      `CFI (Đầu tư)`    = cfi,
      `CFF (Tài chính)` = cff
    ) |>
    tidyr::pivot_longer(-Date, names_to = "Dòng tiền", values_to = "Giá trị") |>
    dplyr::mutate(
      `Dòng tiền` = factor(`Dòng tiền`,
                           levels = c("CFI (Đầu tư)", "CFF (Tài chính)", "CFO (HĐKD)"))
    )

  # ΔCash theo kỳ
  net <- df_cf |>
    dplyr::group_by(Date) |>
    dplyr::summarise(net = sum(`Giá trị`, na.rm = TRUE), .groups = "drop")

  # ---- ĐỘ RỘNG CỘT THEO QUÝ  ----
  uq_dates <- sort(unique(df_cf$Date))
  step_days <- if (length(uq_dates) > 1) median(diff(uq_dates)) else 90
  bar_w <- as.numeric(step_days) * 0.75   # ~ 75% độ rộng mỗi quý

  # Nhãn ΔCash: chỉ hiển thị điểm nổi bật (top 20% | kỳ cuối)  tránh rối
  net <- net |>
    dplyr::mutate(
      keep   = abs(net) >= stats::quantile(abs(net), 0.8, na.rm = TRUE) | Date == max(Date),
      offset = 0.06 * max(abs(net), na.rm = TRUE),
      y_lbl  = net + ifelse(net >= 0, offset, -offset),
      lbl    = scales::number(net, accuracy = 0.1, big.mark = ".", decimal.mark = ","),
      col    = ifelse(net >= 0, "pos", "neg")
    )

  # Bảng màu
  fill_map <- c(
    "CFO (HĐKD)"      = "#16A34A",  # green-600
    "CFI (Đầu tư)"    = "#EF4444",  # red-500
    "CFF (Tài chính)" = "#6366F1"   # indigo-500
  )
  lbl_col  <- c(pos = "#059669", neg = "#DC2626")

  p18 <- ggplot(df_cf, aes(x = Date, y = `Giá trị`)) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf,
             fill = "#F8FAFC", alpha = 0.65) +
    geom_hline(yintercept = 0, color = "#9CA3AF", linetype = "dashed") +
    geom_col(aes(fill = `Dòng tiền`), width = bar_w, color = NA, alpha = 0.92) +
    geom_line(data = net, aes(x = Date, y = net),
              color = "#111827", linewidth = 1.1, inherit.aes = FALSE) +
    geom_point(data = net, aes(x = Date, y = net),
               color = "#111827", size = 1.9, inherit.aes = FALSE) +
    geom_label(
      data = subset(net, keep),
      aes(x = Date, y = y_lbl, label = lbl, color = col),
      fill = "white", label.size = 0, size = 3.1, fontface = "bold",
      show.legend = FALSE, inherit.aes = FALSE
    ) +
    scale_fill_manual(values = fill_map, name = NULL) +
    scale_color_manual(values = lbl_col, guide = "none") +
    scale_y_continuous(labels = num_fmt, expand = expansion(mult = c(0.08, 0.12))) +
    scale_x_date(date_labels = "%Y", date_breaks = "1 year",
                 expand = expansion(mult = c(0.02, 0.04))) +
    coord_cartesian(clip = "off") +
    labs(
      title = "Thác nước Dòng tiền (CFO / CFI / CFF)",
      subtitle = "Cột theo quý (rộng 75% kỳ)|Đường đen: ΔCash|Nhãn chỉ hiển thị kỳ nổi bật.",
      x = NULL, y = "Ngàn Tỷ", caption = cap_default
    ) +
    theme_pro() +
    theme(
      legend.position = "top",
      axis.text.x = element_text(angle = 45, hjust = 1),
      panel.grid.major.x = element_blank(),
      plot.margin = margin(10, 28, 10, 10)
    )

  print(maybe_plotly(p18))  # nếu plotly làm nhãn trùng, đổi tạm thành: print(p18)
  maybe_save(p18, "18_cashflow_waterfall.png")
}

Phân tích Code:

• Dùng factor để sắp xếp thứ tự các dòng tiền. Điều này đảm bảo các cột xếp chồng được vẽ một cách logic.
• Tính toán dòng tiền ròng (net) bằng cách cộng ba dòng tiền thành phần.
• Kết hợp geom_col để vẽ các thành phần và geom_line để vẽ dòng tiền ròng. Kỹ thuật này cho thấy cả nguyên nhân và kết quả.

Kết quả:

*   Đầu năm 2022, dòng tiền ròng đạt đỉnh do CFO rất lớn.
*   Cuối năm 2022, dòng tiền ròng giảm mạnh do CFO và CFI cùng âm.

Kết luận: Mô hình kinh doanh của công ty lành mạnh, với dòng tiền chính đến từ hoạt động kinh doanh.

Chất lượng Lợi nhuận (CFO vs NI)

if (all(c("Date","cfo","ni") %in% names(df_plot))) {

  # hệ số trục phụ an toàn
  s <- safe_coef(df_plot$cfo, df_plot$ni)

  # độ rộng cột theo khoảng cách quý
  uq  <- sort(unique(df_plot$Date))
  step_days <- if (length(uq) > 1) median(diff(uq)) else 90
  bar_w <- as.numeric(step_days) * 0.75

  df_cfo <- df_plot |>
    dplyr::select(Date, cfo, ni) |>
    dplyr::mutate(
      `CFO ký hiệu` = dplyr::case_when(cfo > 0 ~ "CFO (+)", cfo < 0 ~ "CFO (-)", TRUE ~ "CFO (0)"),
      diverge = dplyr::if_else(sign(cfo) * sign(ni) == -1, TRUE, FALSE, missing = FALSE)
    )

  # điểm cuối cho nhãn NI (trục phụ)
  last_ni <- df_cfo |>
    dplyr::filter(!is.na(ni)) |>
    dplyr::slice_tail(n = 1) |>
    dplyr::mutate(
      x_lbl = Date + step_days * 0.7,
      y_lbl = ni * s,
      lbl   = paste0("NI: ", scales::number(ni, accuracy = 0.1, big.mark=".", decimal.mark=",")))
  
  fill_map <- c("CFO (+)"="#16A34A","CFO (-)"="#EF4444","CFO (0)"="#A3A3A3")
  line_col <- "#1D4ED8"  # xanh dương đậm cho NI

  p19 <- ggplot(df_cfo, aes(x = Date)) +
    # nền hồng các quý phân kỳ CFO vs NI
    geom_rect(
      data = dplyr::filter(df_cfo, diverge),
      aes(xmin = Date - step_days*0.5, xmax = Date + step_days*0.5,
          ymin = -Inf, ymax = Inf),
      inherit.aes = FALSE, fill = "#FDE2E2", alpha = 0.35
    ) +
    # đường zero
    geom_hline(yintercept = 0, linetype = "dashed", color = "#9CA3AF") +
    # cột CFO
    geom_col(aes(y = cfo, fill = `CFO ký hiệu`), width = bar_w, alpha = 0.92, color = NA) +
    # đường & điểm NI (được scale lên trục trái rồi trả về bằng sec.axis)
    geom_line(aes(y = ni * s, color = "NI (trục phụ)", group = 1), linewidth = 1.25) +
    geom_point(aes(y = ni * s, color = "NI (trục phụ)"), size = 1.9) +
    # nhãn NI ở kỳ cuối
    geom_label(
      data = last_ni,
      aes(x = x_lbl, y = y_lbl, label = lbl),
      inherit.aes = FALSE, size = 3.1, label.size = 0,
      fill = "white", color = line_col, fontface = "bold"
    ) +
    scale_fill_manual(values = fill_map, name = NULL) +
    scale_color_manual(values = c("NI" = line_col), name = NULL) +
    scale_y_continuous(
      name = "CFO (Ngàn Tỷ)", labels = num_fmt,
      sec.axis = sec_axis(~ . / s, name = "Net Income (Ngàn Tỷ)", labels = num_fmt),
      expand = expansion(mult = c(0.08, 0.12))
    ) +
    scale_x_date(date_breaks = "1 year", date_labels = "%Y",
                 expand = expansion(mult = c(0.02, 0.04))) +
    coord_cartesian(clip = "off") +
    labs(
      title = "Chất lượng Lợi nhuận: CFO vs Net Income",
      x = NULL, caption = cap_default
    ) +
    theme_pro() +
    theme(
      legend.position = "top",
      axis.text.x = element_text(angle = 45, hjust = 1),
      panel.grid.major.x = element_blank(),
      plot.margin = margin(10, 28, 10, 10)
    )

  print(maybe_plotly(p19))  # Nếu tooltip gây rối, dùng print(p19)
  maybe_save(p19, "19_quality_of_earnings.png")
}

Phân tích Code:

• Tạo một biến diverge để xác định các kỳ có sự phân kỳ giữa CFO và Lợi nhuận ròng. Logic này dùng tích của hàm sign() để phát hiện các trường hợp hai chỉ số trái dấu.
• Dùng geom_rect để tô màu nền cho các kỳ có sự phân kỳ. Đây là một cảnh báo trực quan mạnh mẽ.
• Dùng geom_col cho CFO và tô màu xanh/đỏ tùy theo giá trị dương/âm.
• Dùng geom_line và kỹ thuật trục phụ để vẽ Lợi nhuận ròng.

Kết quả:

• Nhìn chung, xu hướng của CFO và Lợi nhuận ròng đồng pha với nhau.
• Tuy nhiên, các vùng màu hồng cho thấy có nhiều kỳ công ty báo lãi nhưng dòng tiền hoạt động lại âm. Đây là một cảnh báo về chất lượng lợi nhuận.
• Ngay cả khi cùng dương, CFO thường thấp hơn Lợi nhuận ròng, cho thấy khả năng chuyển đổi lợi nhuận thành tiền mặt cần được cải thiện.

Tỷ lệ CFO/NI

if (all(c("Date","cfo_ni") %in% names(df_plot))) {

  # kích thước cột mốc theo khoảng cách quý (dùng cho nhãn lệch)
  uq  <- sort(unique(df_plot$Date))
  step_days <- if (length(uq) > 1) median(diff(uq)) else 90

  # dải nền
  min_y <- suppressWarnings(min(df_plot$cfo_ni, na.rm = TRUE))
  max_y <- suppressWarnings(max(df_plot$cfo_ni, na.rm = TRUE))
  max_y <- if (!is.finite(max_y)) 2 else max(max_y, 2)

  # điểm cuối để gắn nhãn
  last_pt <- df_plot |>
    dplyr::filter(!is.na(cfo_ni)) |>
    dplyr::slice_tail(n = 1) |>
    dplyr::mutate(
      x_lbl = Date + step_days*0.7,
      y_lbl = cfo_ni,
      lbl   = paste0(scales::number(cfo_ni, accuracy = 0.01, big.mark=".", decimal.mark=","), "×")
    )

  line_col  <- "#0EA5E9"  # cyan-600
  loess_col <- "#F59E0B"  # amber-500

  p20 <- ggplot(df_plot, aes(x = Date, y = cfo_ni, group = 1)) +
    # nền: yếu (<1) đỏ nhạt, khỏe (>=1) xanh nhạt
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = 1, ymax = Inf,
             fill = "#ECFDF5", alpha = 0.45) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = 1,
             fill = "#FDE2E2", alpha = 0.45) +
    # nếu có phần <0, tô đỏ đậm hơn
    { if (is.finite(min_y) && min_y < 0)
        annotate("rect", xmin = -Inf, xmax = Inf, ymin = min_y, ymax = 0,
                 fill = "#FECACA", alpha = 0.45) } +
    # đường tham chiếu
    geom_hline(yintercept = 0, color = "#9CA3AF", linetype = "dotted") +
    geom_hline(yintercept = 1, color = "#6B7280", linetype = "dashed") +
    geom_hline(yintercept = 2, color = "#CBD5E1", linetype = "dotted") +
    # chuỗi chính + LOESS xu hướng
    geom_line(color = line_col, linewidth = 1.2) +
    geom_point(color = line_col, size = 1.9) +
    geom_smooth(method = "loess", se = FALSE, color = loess_col, linewidth = 1) +
 
    geom_label(
      data = last_pt,
      aes(x = x_lbl, y = y_lbl, label = lbl),
      inherit.aes = FALSE, size = 3.1, label.size = 0,
      fill = "white", color = line_col, fontface = "bold"
    ) +
    scale_y_continuous(labels = qty_fmt, expand = expansion(mult = c(0.08, 0.12))) +
    scale_x_date(date_breaks = "1 year", date_labels = "%Y",
                 expand = expansion(mult = c(0.02, 0.04))) +
    coord_cartesian(clip = "off") +
    labs(
      title = "Tỷ lệ CFO / NI",
      subtitle = "Nền đỏ: Tỷ lệ < 1X | Đường vàng: LOESS xu hướng.",
      x = NULL, y = "Lần", caption = cap_default
    ) +
    theme_pro() +
    theme(
      panel.grid.major.x = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1),
      plot.margin = margin(10, 28, 10, 10)
    )

  print(maybe_plotly(p20))   
  maybe_save(p20, "20_cfo_ni_ratio.png")
}

Phân tích Code:

• Dùng annotate("rect", ...) để tạo các vùng nền màu xanh/hồng/đỏ. Các vùng này cung cấp một bối cảnh trực quan, giúp người xem đánh giá ngay lập tức mức độ “khỏe mạnh” của chỉ số.
• Dùng geom_hline để vẽ các đường tham chiếu tại các ngưỡng quan trọng (0 và 1).
• Kết hợp geom_line để vẽ dữ liệu gốc và geom_smooth để vẽ một đường xu hướng dài hạn. Đường xu hướng giúp Ta nhìn xuyên qua các biến động ngắn hạn.

Kết quả:

• Khả năng chuyển đổi lợi nhuận thành tiền mặt của công ty rất không ổn định, thể hiện qua sự biến động mạnh của đường màu xanh.
• Đường xu hướng dài hạn (màu vàng) cho thấy chất lượng lợi nhuận đã cải thiện trong giai đoạn 2019-2022 nhưng đang có dấu hiệu suy giảm trở lại.
• Phần lớn thời gian, tỷ lệ này nằm dưới mức 1 (vùng màu hồng). Điều này xác nhận rằng chất lượng lợi nhuận là một điểm yếu của công ty.

Nội suy tuyến tính

Ma trận Tương quan (Heatmap)

key_cols <- intersect(
  c("revenue_yoy","net_margin","roe","ccc","debt_to_equity","fcf_margin","asset_turnover"),
  names(df_plot)
)

if (length(key_cols) >= 3) {
  sub <- df_plot |>
    dplyr::select(dplyr::all_of(key_cols)) |>
    dplyr::mutate(dplyr::across(dplyr::everything(), as.numeric))

  cm <- cor(sub, use = "pairwise.complete.obs")
  cm[is.na(cm)] <- 0

  # Sắp xếp theo phân cụm dựa trên |cor|
  ord <- hclust(as.dist(1 - abs(cm)))$order
  cm_ord <- cm[ord, ord]

  # Chỉ hiển thị tam giác dưới
  cm_ord[upper.tri(cm_ord, diag = FALSE)] <- NA
  dfc <- as.data.frame(as.table(cm_ord))
  names(dfc) <- c("Var1","Var2","Corr")
  dfc <- dplyr::filter(dfc, !is.na(Corr))
  dfc$Var1 <- factor(dfc$Var1, levels = rownames(cm_ord))
  dfc$Var2 <- factor(dfc$Var2, levels = colnames(cm_ord))

  p22 <- ggplot(dfc, aes(x = Var2, y = Var1, fill = Corr)) +
    # Ô vuông + lưới trắng mảnh
    geom_tile(width = .95, height = .95, color = "white", linewidth = .6) +
    # Nhãn số; màu chữ tự động: trắng khi |Corr| lớn
    geom_text(
      aes(
        label = scales::number(Corr, accuracy = 0.01, decimal.mark = ","),
        color = ifelse(abs(Corr) >= 0.5, "hi", "lo")
      ),
      size = 3.2, fontface = "bold"
    ) +
    # Palette phân kỳ (đỏ–trắng–xanh), cố định [-1,1]
    scale_fill_gradient2(
      low = "#EF4444", mid = "#F8FAFC", high = "#22C55E",
      midpoint = 0, limits = c(-1, 1), name = "Tương quan"
    ) +
    scale_color_manual(values = c(hi = "white", lo = "#374151"), guide = "none") +
    coord_fixed() +
    labs(
      title = "Ma trận tương quan ",
      x = NULL, y = NULL, caption = cap_default
    ) +
    theme_pro() +
    theme(
      legend.position = "right",
      axis.text.x = element_text(angle = 45, hjust = 1),
      panel.grid = element_blank(),
      plot.margin = margin(10, 20, 10, 10)
    )

  print(p22)                  # hoặc maybe_plotly(p22) nếu muốn tương tác
  maybe_save(p22, "22_correlation_heatmap.png")
}

Phân tích Code:

• Dùng hàm cor() để tính ma trận tương quan giữa các biến.
• Sử dụng kỹ thuật Phân cụm theo thứ bậc (hclust) để sắp xếp lại ma trận. Việc này nhóm các biến có tương quan cao lại gần nhau, làm nổi bật các cụm quan hệ.
• Dùng geom_tile để vẽ bản đồ nhiệt.
• Dùng scale_fill_gradient2, một thang màu phân kỳ. Màu xanh cho tương quan thuận, màu đỏ cho tương quan nghịch, và màu trắng cho giá trị gần 0.

Kết quả:

• Cụm Lợi nhuận & Hiệu quả: ROE có tương quan thuận rất mạnh với Biên lợi nhuận ròng và Vòng quay tài sản. Điều này xác nhận kết quả từ phân tích DuPont.
• Mối quan hệ Nghịch chiều mạnh nhất: Chu kỳ chuyển đổi tiền mặt (CCC) có tương quan nghịch rất mạnh với Biên lợi nhuận ròng và Vòng quay tài sản. Khi công ty quản lý vốn lưu động hiệu quả hơn (CCC thấp hơn), lợi nhuận và hiệu quả của họ cao hơn.
• Tăng trưởng và Đòn bẩy: Tăng trưởng doanh thu không có tương quan mạnh với lợi nhuận. Đòn bẩy tài chính có tương quan nghịch với các chỉ số hiệu quả.

Phân rã STL (Revenue hoặc NI)

series_name <- if ("revenue" %in% names(df_plot)) "revenue" else if ("ni" %in% names(df_plot)) "ni" else NA_character_

if (!is.na(series_name) && sum(!is.na(df_plot[[series_name]])) >= 8) {
  y0 <- suppressWarnings(min(df_plot$yr, na.rm = TRUE))
  q0 <- df_plot$qtr[which.min(df_plot$Date)]
  if (!is.finite(y0) || is.na(y0)) y0 <- 2000
  if (is.na(q0)) q0 <- 1

  ts_obj <- stats::ts(df_plot[[series_name]], frequency = 4, start = c(y0, q0))
  fit <- try(stats::stl(ts_obj, s.window = "periodic", robust = TRUE), silent = TRUE)

  if (!inherits(fit, "try-error")) {
    comp_wide <- dplyr::tibble(
      Date      = df_plot$Date,
      observed  = as.numeric(fit$time.series[, "trend"] + fit$time.series[, "seasonal"] + fit$time.series[, "remainder"]),
      trend     = as.numeric(fit$time.series[, "trend"]),
      seasonal  = as.numeric(fit$time.series[, "seasonal"]),
      remainder = as.numeric(fit$time.series[, "remainder"])
    )

    # dữ liệu dài + panel
    comp <- comp_wide |>
      tidyr::pivot_longer(-Date, names_to = "name", values_to = "value") |>
      dplyr::mutate(
        panel = dplyr::recode(name,
          observed  = "Observed & Trend",
          trend     = "Observed & Trend",
          seasonal  = "Seasonal",
          remainder = "Remainder"
        ),
        sign_season = dplyr::case_when(name == "seasonal" & value >= 0 ~ "Dương",
                                       name == "seasonal" & value <  0 ~ "Âm",
                                       TRUE ~ NA_character_)
      )

    # nhãn Trend kỳ gần nhất
    last_trend <- comp |>
      dplyr::filter(name == "trend", !is.na(value)) |>
      dplyr::slice_tail(n = 1) |>
      dplyr::mutate(lbl = scales::number(value, accuracy = 0.1, big.mark=".", decimal.mark=","))

    # màu pro
    col_obs <- "#0EA5E9"  # cyan
    col_trd <- "#F59E0B"  # amber
    col_rem <- "#8B5CF6"  # violet
    fill_pos <- "#DCFCE7" # green-100
    fill_neg <- "#FEE2E2" # red-100

    p23 <- ggplot(comp, aes(x = Date)) +
      # Seasonal ribbon (theo dấu)
      geom_ribbon(
        data = dplyr::filter(comp, name == "seasonal"),
        aes(ymin = pmin(value, 0), ymax = pmax(value, 0), fill = sign_season),
        alpha = 0.6, color = NA
      ) +
      # đường zero cho Seasonal & Remainder
      geom_hline(
        data = subset(comp, panel %in% c("Seasonal","Remainder")) |> dplyr::distinct(panel),
        aes(yintercept = 0), linetype = "dashed", color = "#9CA3AF"
      ) +
      # Observed
      geom_line(
        data = dplyr::filter(comp, name == "observed"),
        aes(y = value, color = "Observed"), linewidth = 0.9
      ) +
      # Trend
      geom_line(
        data = dplyr::filter(comp, name == "trend"),
        aes(y = value, color = "Trend"), linewidth = 1.25
      ) +
      # Remainder (điểm + đường mảnh)
      geom_line(
        data = dplyr::filter(comp, name == "remainder"),
        aes(y = value, color = "Remainder"), linewidth = 0.7, alpha = 0.8
      ) +
      geom_point(
        data = dplyr::filter(comp, name == "remainder"),
        aes(y = value, color = "Remainder"), size = 1.2, alpha = 0.9
      ) +
      # nhãn Trend
      geom_label(
        data = last_trend,
        aes(x = Date + 20, y = value, label = lbl),
        inherit.aes = FALSE, size = 3.1, label.size = 0,
        fill = "white", color = col_trd, fontface = "bold"
      ) +
      facet_wrap(~ panel, ncol = 1, scales = "free_y") +
      scale_fill_manual(values = c("Dương" = fill_pos, "Âm" = fill_neg), guide = "none") +
      scale_color_manual(values = c("Observed" = col_obs, "Trend" = col_trd, "Remainder" = col_rem), name = NULL) +
      scale_y_continuous(labels = num_fmt, breaks = scales::breaks_pretty(n = 3),
                         expand = expansion(mult = c(0.08, 0.12))) +
      scale_x_date(date_breaks = "1 year", date_labels = "%Y",
                   expand = expansion(mult = c(0.02, 0.04))) +
      coord_cartesian(clip = "off") +
      labs(
        title = paste0("Phân rã STL: ", toupper(series_name)),
        x = NULL, y = NULL, caption = cap_default
      ) +
      theme_pro() +
      theme(
        legend.position = "top",
        panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1),
        plot.margin = margin(10, 24, 10, 10)
      )

    print(p23)                       # hoặc maybe_plotly(p23) nếu muốn tương tác
    maybe_save(p23, "23_stl_decomposition.png", w = 9, h = 8)
  }
}

Phân tích Code:

• Chuyển đổi dữ liệu thành một đối tượng chuỗi thời gian (ts) với tần suất 4 kỳ mỗi năm. Đây là yêu cầu của hàm stl.
• Áp dụng thuật toán stl để phân rã chuỗi thời gian thành ba thành phần: xu hướng, mùa vụ và phần dư.
• Dùng facet_wrap để vẽ mỗi thành phần trên một biểu đồ nhỏ riêng biệt.
• Dùng scales = "free_y". Điều này là cần thiết vì mỗi thành phần có một thang đo khác nhau.

Kết quả:

• Xu hướng (Trend): Đường màu vàng trong panel đầu tiên cho thấy sự tăng trưởng dài hạn đã được làm mịn, loại bỏ các biến động theo quý.
• Mùa vụ (Seasonal): Panel cuối cùng cho thấy một mẫu hình lặp lại hàng năm, với các quý mạnh và quý yếu rõ ràng.
• Phần dư (Remainder): Panel ở giữa làm nổi bật các sự kiện bất thường. Cú sụt giảm mạnh vào cuối năm 2021 là một sự kiện ngoại lệ, không phải do xu hướng hay mùa vụ.

Scatter Matrix (P&L chính)

cols <- intersect(c("revenue","gp","ebit","ni"), names(df_plot))
if (length(cols) >= 3) {

  # Dữ liệu & nhãn cột
  df_pairs <- df_plot[, cols, drop = FALSE]
  vi_labels <- c(revenue = "Doanh thu", gp = "Ln gộp", ebit = "EBIT", ni = "Ln ròng")
  colnames(df_pairs) <- vi_labels[cols]

  # Panel trên: hệ số tương quan có nền theo độ mạnh (đỏ âm, xanh dương)
  upper_fun <- function(data, mapping, ...){
    x <- as.numeric(rlang::eval_tidy(mapping$x, data))
    y <- as.numeric(rlang::eval_tidy(mapping$y, data))
    c <- suppressWarnings(stats::cor(x, y, use = "pairwise.complete.obs"))
    c <- ifelse(is.finite(c), c, NA_real_)

    fill_fun <- scales::col_numeric(
      palette  = c("#EF4444", "#F8FAFC", "#22C55E"),
      domain   = c(-1, 1),
      na.color = "#E5E7EB"
    )
    bg  <- fill_fun(c)
    txt <- if (!is.na(c) && abs(c) >= 0.5) "white" else "#111827"

    ggplot() +
      annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, fill = bg) +
      annotate("text", x = 0.5, y = 0.5,
               label = scales::number(c, accuracy = 0.01, decimal.mark = ","),
               colour = txt, size = 5, fontface = "bold") +
      theme_void()
  }

  # Panel dưới: scatter + đường hồi quy
  lower_fun <- function(data, mapping, ...){
    ggplot(data = data, mapping = mapping) +
      geom_point(alpha = 0.6, size = 1.6, color = "#0EA5E9") +
      geom_smooth(method = "lm", se = FALSE, linewidth = 0.9, color = "#F59E0B") +
      theme_minimal(base_size = 10) +
      theme(panel.grid.minor = element_blank(),
            panel.grid.major = element_line(color = "#ECEFF4"))
  }

  # Đường chéo: mật độ (density)
  diag_fun <- GGally::wrap("densityDiag", alpha = 0.6, fill = "#8B5CF6", adjust = 1)

  p24 <- GGally::ggpairs(
    df_pairs,
    upper = list(continuous = upper_fun),
    lower = list(continuous = lower_fun),
    diag  = list(continuous = diag_fun)
  ) +
    theme_pro() +
    theme(
      strip.text       = element_text(face = "bold", size = 10),
      panel.spacing    = unit(0.6, "lines"),
      axis.text.x      = element_text(size = 8, angle = 45, hjust = 1),
      axis.text.y      = element_text(size = 8),
      plot.title       = element_text(size = 14),
      plot.caption     = element_text(size = 9),
      panel.grid.major = element_blank()
    ) +
    labs(
      title   = "Ma trận phân tán ",
      caption = cap_default )
  print(p24)
  if (save_plots) ggsave(file.path('figs', '24_pairs_plot.png'), p24, width = 9, height = 9, dpi = 300)
}

Phân tích Code:

• Dùng hàm GGally::ggpairs để tạo một lưới biểu đồ N x N.
• Tùy chỉnh từng khu vực của ma trận:
  • Tam giác trên: Hiển thị hệ số tương quan dưới dạng bản đồ nhiệt.
  • Tam giác dưới: Hiển thị biểu đồ phân tán với một đường hồi quy tuyến tính.
  • Đường chéo: Hiển thị biểu đồ mật độ để xem phân phối của từng biến.

Kết quả:

• Phân phối (Đường chéo): Các biến đều có phân phối lệch phải. Điều này có nghĩa là có một vài quý với doanh thu và lợi nhuận rất cao.
• Tương quan (Tam giác trên): Tất cả các chỉ số lợi nhuận đều có tương quan thuận rất mạnh với nhau. Khi doanh thu tăng, các loại lợi nhuận khác cũng tăng theo.
• Mối quan hệ Tuyến tính (Tam giác dưới): Các điểm dữ liệu nằm rất sát đường hồi quy. Điều này xác nhận mối quan hệ giữa chúng là tuyến tính mạnh mẽ, không có các điểm ngoại lệ rõ rệt.

Doanh thu Rolling 4Q

if (all(c("Date","revenue") %in% names(df_plot))) {
  df_roll <- df_plot |>
    dplyr::mutate(
      revenue_4Q = zoo::rollsum(revenue, k = 4, fill = NA, align = "right"),
      sign_4Q    = dplyr::case_when(
        is.na(revenue_4Q) ~ NA_character_,
        revenue_4Q >= 0   ~ "≥ 0",
        TRUE              ~ "< 0"
      )
    )

  if (sum(!is.na(df_roll$revenue_4Q)) >= 2) {
    # bước thời gian để canh nhãn
    uq <- sort(unique(df_roll$Date))
    step_days <- if (length(uq) > 1) median(diff(uq)) else 90

    last_pt <- df_roll |>
      dplyr::filter(!is.na(revenue_4Q)) |>
      dplyr::slice_tail(n = 1) |>
      dplyr::mutate(
        x_lbl = Date + step_days*0.7,
        lbl   = paste0("4Q: ", scales::number(revenue_4Q, accuracy = 0.1, big.mark=".", decimal.mark=","))
      )

    col_line  <- "#0EA5E9"  # cyan-600
    col_loess <- "#F59E0B"  # amber-500

    p26 <- ggplot(df_roll, aes(x = Date)) +
      # ribbon phân cực theo dấu của rolling 4Q
      geom_ribbon(
        aes(ymin = pmin(revenue_4Q, 0), ymax = pmax(revenue_4Q, 0), fill = sign_4Q),
        alpha = 0.55, color = NA
      ) +
      # đường zero mảnh
      geom_hline(yintercept = 0, linetype = "dotted", color = "#9CA3AF") +
      # đường rolling 4Q + điểm
      geom_line(aes(y = revenue_4Q), color = col_line, linewidth = 1.4) +
      geom_point(aes(y = revenue_4Q), color = col_line, size = 1.7) +
      # xu hướng loess
      geom_smooth(aes(y = revenue_4Q), method = "loess", se = FALSE, color = col_loess, linewidth = 1) +
      # nhãn kỳ mới nhất
      geom_label(
        data = last_pt,
        aes(x = x_lbl, y = revenue_4Q, label = lbl),
        inherit.aes = FALSE, size = 3.1, label.size = 0,
        fill = "white", color = col_line, fontface = "bold"
      ) +
      scale_fill_manual(values = c("≥ 0" = "#DCFCE7", "< 0" = "#EF4444"), guide = "none") +
      scale_y_continuous(labels = num_fmt, expand = expansion(mult = c(0.08, 0.12))) +
      scale_x_date(date_breaks = "1 year", date_labels = "%Y",
                   expand = expansion(mult = c(0.02, 0.04))) +
      coord_cartesian(clip = "off") +
      labs(
        title = "Doanh thu Rolling 4Q (TTM)",
        subtitle = "Ribbon xanh/đỏ theo dấu|Đường vàng là LOESS xu hướng",
        x = NULL, y = "Ngàn Tỷ", caption = cap_default
      ) +
      theme_pro() +
      theme(
        panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1),
        plot.margin = margin(10, 24, 10, 10)
      )

    print(maybe_plotly(p26))   # hoặc print(p26) nếu cần tĩnh tuyệt đối
    maybe_save(p26, "26_revenue_rolling4Q.png")
  }
}

Phân tích Code:

• Dùng zoo::rollsum(..., k = 4, ...) để tính tổng doanh thu của 4 quý gần nhất. Thao tác này tạo ra dữ liệu Trailing Twelve Months (TTM).
• Dùng geom_ribbon để tô màu vùng dưới đường doanh thu TTM. Kỹ thuật này tạo cảm giác về quy mô và làm cho sự tăng trưởng trở nên ấn tượng hơn.
• Kết hợp geom_line (dữ liệu TTM) và geom_smooth (xu hướng của TTM). geom_smooth hoạt động như một lớp làm mịn thứ hai, giúp xác định sự thay đổi trong tốc độ tăng trưởng.

Kết quả:

• Biểu đồ này mượt mà hơn nhiều so với dữ liệu quý gốc. Các biến động mùa vụ đã được loại bỏ, giúp Ta tập trung vào xu hướng tăng trưởng cốt lõi.
• Xu hướng tăng trưởng rõ ràng từ năm 2018.
• Có một giai đoạn tăng tốc mạnh mẽ từ giữa năm 2022 đến 2024.
• Dữ liệu TTM là một tiêu chuẩn trong tài chính. Nó cập nhật hơn dữ liệu năm nhưng ổn định hơn dữ liệu quý.

Biên lợi nhuận

cols <- intersect(c("gross_margin","operating_margin","net_margin"), names(df_plot))
if (length(cols) >= 2) {

  dfm <- df_plot |>
    dplyr::select(Date, dplyr::all_of(cols)) |>
    tidyr::pivot_longer(-Date, names_to = "Metric", values_to = "Value") |>
    dplyr::mutate(
      Metric = dplyr::recode(Metric,
        gross_margin     = "Biên gộp",
        operating_margin = "Biên EBIT",
        net_margin       = "Biên ròng"
      )
    )

  uq <- sort(unique(dfm$Date))
  step_days <- if (length(uq) > 1) stats::median(diff(uq)) else 90

  last_pts <- dfm |>
    dplyr::filter(!is.na(Value)) |>
    dplyr::group_by(Metric) |>
    dplyr::slice_tail(n = 1) |>
    dplyr::ungroup() |>
    dplyr::mutate(
      x_lbl = Date + step_days*0.9,
      lbl   = paste0(Metric, ": ",
                     scales::percent(Value, accuracy = 0.1, big.mark=".", decimal.mark=","))
    )

  pal_margin <- c(
    "Biên ròng" = "#8B5CF6",
    "Biên EBIT" = "#10B981",
    "Biên gộp"  = "#0EA5E9"
  )

  base <- ggplot(dfm, aes(x = Date, y = Value, color = Metric, group = Metric)) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0, ymax = Inf,
             fill = "#EEFDF3", alpha = 0.30) +
    annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0,
             fill = "#F5EEFF", alpha = 0.30) +
    geom_line(linewidth = 2.2, alpha = 0.10, lineend = "round") +
    geom_line(linewidth = 1.25, lineend = "round") +
    # ==== LOESS ====
    geom_smooth(se = FALSE, method = "loess", linewidth = 1, linetype = "dashed") +
    # ===============
    scale_color_manual(values = pal_margin) +
    scale_y_continuous(labels = pct_fmt, expand = expansion(mult = c(0.08, 0.12))) +
    scale_x_date(date_breaks = "1 year", date_labels = "%Y",
                 expand = expansion(mult = c(0.02, 0.18))) +
    labs(
      title = "Biên lợi nhuận (LOESS)",
      subtitle = "Nét đứt là đường Xu Hướng",
      x = NULL, y = "Tỷ lệ", caption = cap_default
    ) +
    theme_pro() +
    theme(
      panel.grid.major.x = element_blank(),
      axis.text.x = element_text(angle = 45, hjust = 1),
      plot.margin = margin(10, 28, 10, 10)
    )

  if (requireNamespace("ggrepel", quietly = TRUE)) {
    p27 <- base +
      ggrepel::geom_label_repel(
        data = last_pts,
        aes(x = x_lbl, y = Value, label = lbl, color = Metric),
        inherit.aes = FALSE,
        size = 3.1, label.size = 0, fill = "white", fontface = "bold",
        nudge_x = 0, nudge_y = 0,
        min.segment.length = 0, segment.alpha = 0.4,
        direction = "y", box.padding = 0.25, point.padding = 0.3,
        max.overlaps = Inf
      ) +
      coord_cartesian(clip = "off")
  } else {
    message("[INFO] ggrepel chưa cài đặt: dùng geom_label thông thường.")
    p27 <- base +
      geom_label(
        data = last_pts,
        aes(x = x_lbl, y = Value, label = lbl, color = Metric),
        inherit.aes = FALSE, size = 3.1, label.size = 0,
        fill = "white", fontface = "bold"
      ) +
      coord_cartesian(clip = "off")
  }

  print(maybe_plotly(p27))
  maybe_save(p27, "27_margins_loess.png")
}

Phân tích Code:

• Dùng pivot_longer để chuyển dữ liệu sang dạng dài, giúp ggplot có thể tự động vẽ và tô màu ba đường riêng biệt.
• Kết hợp geom_line (dữ liệu gốc) và geom_smooth (xu hướng). Việc dùng đường nét đứt cho xu hướng giúp phân biệt rõ ràng hai loại thông tin.
• Dùng ggrepel để dán nhãn trực tiếp lên cuối mỗi đường. Kỹ thuật này hiệu quả hơn một chú thích riêng biệt vì người xem không cần phải tra cứu.

Kết quả:

• Có một sự phân cấp rõ ràng: Biên gộp > Biên EBIT > Biên ròng. Khoảng cách giữa các đường đại diện cho các lớp chi phí.
• Biên gộp tương đối ổn định. Điều này cho thấy công ty kiểm soát giá vốn tốt.
• Biên EBIT và Biên ròng biến động mạnh hơn, với một cú sụt giảm nghiêm trọng vào cuối năm 2021.

Kết luận: Vì Biên gộp ổn định, sự biến động của các biên lợi nhuận thấp hơn phải đến từ chi phí hoạt động. Các đường xu hướng cho thấy về dài hạn, các biên lợi nhuận tương đối ổn định.

Quỹ đạo Margin và CCC

req <- c("Date","net_margin","ccc")
if (all(req %in% names(df_plot))) {

  # ===== Data & TTM for size =====
  df_b3 <- df_plot |>
    dplyr::select(Date, net_margin, ccc, revenue) |>
    dplyr::arrange(Date) |>
    dplyr::mutate(
      revenue_ttm = if ("revenue" %in% names(df_plot))
        zoo::rollsum(revenue, k = 4, fill = NA, align = "right") else NA_real_,
      Qtr   = zoo::as.yearqtr(Date),
      label_q = format(Qtr, "Q%q %Y")
    )

  # ===== Medians for regimes =====
  med_nm  <- stats::median(df_b3$net_margin, na.rm = TRUE)
  med_ccc <- stats::median(df_b3$ccc,        na.rm = TRUE)

  # ===== Start/End points for labels =====
  start_pt <- df_b3 |>
    dplyr::filter(!is.na(net_margin), !is.na(ccc)) |>
    dplyr::slice_head(n = 1) |>
    dplyr::mutate(lbl = paste0("Bắt đầu\n", label_q))
  end_pt   <- df_b3 |>
    dplyr::filter(!is.na(net_margin), !is.na(ccc)) |>
    dplyr::slice_tail(n = 1) |>
    dplyr::mutate(lbl = paste0("Hiện tại\n", label_q))

  # ===== Spans for smart nudges (x reversed) =====
  rng_x <- range(df_b3$ccc, na.rm = TRUE);  x_span <- diff(rng_x);  if (!is.finite(x_span) || x_span == 0) x_span <- 1
  rng_y <- range(df_b3$net_margin, na.rm = TRUE); y_span <- diff(rng_y); if (!is.finite(y_span) || y_span == 0) y_span <- 0.1
  nudge_right <- -0.10 * x_span   # x reversed => right = negative nudge
  nudge_left  <-  0.10 * x_span

  time_cols <- c("#2563EB", "#8B5CF6", "#F59E0B", "#EF4444")

  pB3 <- ggplot(df_b3, aes(x = ccc, y = net_margin)) +
    # Regime background (4 quadrants)
    annotate("rect", xmin = -Inf, xmax = med_ccc, ymin = med_nm,  ymax = Inf, fill = "#E8F5E9", alpha = .65) +
    annotate("rect", xmin =  med_ccc, xmax =  Inf,  ymin = med_nm,  ymax = Inf, fill = "#FFF8E1", alpha = .65) +
    annotate("rect", xmin = -Inf, xmax = med_ccc, ymin = -Inf,    ymax = med_nm, fill = "#E3F2FD", alpha = .65) +
    annotate("rect", xmin =  med_ccc, xmax =  Inf,  ymin = -Inf,    ymax = med_nm, fill = "#FDE2E2", alpha = .65) +

    geom_vline(xintercept = med_ccc, linetype = "dotted", color = "#94A3B8") +
    geom_hline(yintercept = med_nm,  linetype = "dotted", color = "#94A3B8") +

    # Trajectory with arrow
    geom_path(linewidth = .9, color = "#374151",
              arrow = grid::arrow(length = grid::unit(3, "mm"), type = "closed"), na.rm = TRUE) +

    # Halo + points
    geom_point(color = "white", size = 4.0, alpha = 0.8, na.rm = TRUE) +
    geom_point(aes(color = as.numeric(Date), size = revenue_ttm), alpha = 0.95, na.rm = TRUE) +

    # ===== Non-overlapping labels (ggrepel) =====
    ggrepel::geom_label_repel(
      data = start_pt, aes(label = lbl), seed = 2,
      nudge_x = nudge_left,  nudge_y = -0.06 * y_span,
      force = 2, max.time = 1.5, box.padding = 0.35, point.padding = 0.35,
      min.segment.length = 0, direction = "both",
      segment.color = "#CBD5E1", segment.size = 0.3,
      fill = "white", color = "#111827", size = 3.2, label.size = 0, show.legend = FALSE
    ) +
    ggrepel::geom_label_repel(
      data = end_pt, aes(label = lbl), seed = 3,
      nudge_x = nudge_right, nudge_y =  0.06 * y_span,
      force = 2, max.time = 1.5, box.padding = 0.35, point.padding = 0.35,
      min.segment.length = 0, direction = "both",
      segment.color = "#CBD5E1", segment.size = 0.3,
      fill = "white", color = "#111827", size = 3.2, label.size = 0, show.legend = FALSE
    ) +

    # Scales
    scale_size_continuous(range = c(1.4, 5.5), name = "Doanh thu (TTM)",
                          labels = scales::label_number(big.mark = ".", decimal.mark = ",")) +
    scale_color_gradientn(colours = time_cols, name = "Thời gian",
                          breaks = as.numeric(unique(as.Date(cut(df_b3$Date, "year")))),
                          labels = format(unique(as.Date(cut(df_b3$Date, "year"))), "%Y")) +

    # Extra right padding for labels (x reversed: pad on second mult)
    scale_x_reverse(labels = scales::label_number(accuracy = 1, big.mark = ".", decimal.mark = ","),
                    expand = expansion(mult = c(0.06, 0.32)),
                    name = "CCC (ngày)") +
    scale_y_continuous(labels = pct_fmt, expand = expansion(mult = c(0.06, 0.16)),
                       name = "Net margin") +

    labs(title = "Regime Map: Net margin vs CCC (Quỹ đạo theo thời gian)",
         caption = cap_default) +
    theme_pro() +
    theme(legend.position = "right", plot.title.position = "plot") +
    coord_cartesian(clip = "off")

  print(pB3)
  maybe_save(pB3, "B3_regime_map_margin_ccc_nolabeloverlap.png", w = 9, h = 6)
}

Phân tích Code:

• Dùng giá trị trung vị để chia biểu đồ thành bốn vùng chiến lược, hay còn gọi là các “chế độ” kinh doanh.
• Đảo ngược trục CCC. Việc này đặt vùng lý tưởng (CCC thấp, Net Margin cao) vào góc trên bên phải, giúp Ta dễ đọc biểu đồ.
• Mã hóa bốn chiều dữ liệu vào biểu đồ: Vị trí X là hiệu quả (CCC), vị trí Y là lợi nhuận (Net Margin), màu sắc là thời gian, và kích thước là quy mô doanh thu.

Kết quả:

• Đường quỹ đạo cho thấy hành trình của công ty Ta.
• Biểu đồ xác nhận cuộc khủng hoảng vào cuối năm 2021. Điểm dữ liệu tại thời điểm đó rơi vào vùng tệ nhất (CCC cao, Net Margin âm).
• Sau đó, công ty đã phục hồi và quay trở lại vùng lý tưởng. Các điểm dữ liệu gần đây cho thấy sự cải thiện.
• Kích thước các điểm lớn dần. Điều này cho thấy quy mô doanh thu tăng mà không phải hy sinh hiệu quả hay lợi nhuận.

Phân rã ROE dạng vòng

need <- c("Date","roe","net_margin","asset_turnover","financial_leverage")
if (all(need %in% names(df_plot))) {

  # === 1) Tính đóng góp DuPont (xấp xỉ tuyến tính quanh t-1) ===
  d <- df_plot |>
    dplyr::select(Date, roe, net_margin, asset_turnover, financial_leverage) |>
    dplyr::arrange(Date) |>
    dplyr::mutate(
      nm_lag = dplyr::lag(net_margin),
      at_lag = dplyr::lag(asset_turnover),
      fl_lag = dplyr::lag(financial_leverage),
      roe_lag= dplyr::lag(roe),
      d_nm = net_margin - nm_lag,
      d_at = asset_turnover - at_lag,
      d_fl = financial_leverage - fl_lag,
      contrib_nm = at_lag * fl_lag * d_nm,
      contrib_at = nm_lag * fl_lag * d_at,
      contrib_fl = nm_lag * at_lag * d_fl,
      delta_roe  = roe - roe_lag,
      approx_delta = contrib_nm + contrib_at + contrib_fl,
      resid = delta_roe - approx_delta
    )

  # === 2) Dài + thang điểm phần trăm (điểm %) ===
  d_long <- d |>
    dplyr::select(Date, contrib_nm, contrib_at, contrib_fl) |>
    tidyr::pivot_longer(-Date, names_to = "component", values_to = "value") |>
    dplyr::mutate(
      component = dplyr::recode(component,
        contrib_nm = "Net margin",
        contrib_at = "Asset turnover",
        contrib_fl = "Financial leverage"
      ),
      value_pp = value * 100,
      Date_f   = factor(Date, levels = unique(Date))
    ) |>
    dplyr::filter(!is.na(value_pp))

  # Tổng (≈ ΔROE, bỏ phần dư) để chấm một vòng mảnh phía ngoài
  total_pp <- d_long |>
    dplyr::group_by(Date, Date_f) |>
    dplyr::summarise(total_pp = sum(value_pp, na.rm = TRUE), .groups = "drop")

  # Bảng màu
  pal_bridge <- c(
    "Net margin" = "#0EA5E9",      # sky
    "Asset turnover" = "#F59E0B",  # amber
    "Financial leverage" = "#8B5CF6" # violet
  )

  # === 3) Polar stream (hiếm gặp): geom_col xếp lớp + coord_polar ===
  pB4_radial <- ggplot(d_long, aes(x = Date_f, y = value_pp, fill = component)) +
    # Cột xếp lớp (cả dương lẫn âm) — độ rộng vừa đủ để tạo dòng chảy mượt
    geom_col(width = 0.98, alpha = 0.95, color = NA, position = "stack") +

    # Điểm mảnh biểu thị tổng (≈ ΔROE) mỗi kỳ
    geom_point(data = total_pp, aes(x = Date_f, y = total_pp, color = "Tổng (≈ΔROE)"),
               size = 1.6, inherit.aes = FALSE) +

    scale_fill_manual(values = pal_bridge, name = "Đóng góp ΔROE") +
    scale_color_manual(values = c("Tổng (≈ΔROE)" = "#111827"), name = NULL) +

    # Vòng tròn: theo trục thời gian (theta = x). Bắt đầu ở đỉnh, quay thuận chiều kim đồng hồ.
    coord_polar(theta = "x", start = pi/2, direction = -1, clip = "off") +

    labs(
      title = "Radial ROE Bridge",
      subtitle = "Mỗi nan quạt là một quý",
      x = NULL, y = NULL,
      caption = "Đơn vị: điểm %. ΔROE xấp xỉ tổng đóng góp (bỏ phần dư nhỏ)."
    ) +
    theme_pro() +
    theme(
      legend.position = "top",
      axis.text = element_blank(),
      panel.grid = element_blank(),
      plot.margin = margin(20, 20, 20, 20)
    )

  print(pB4_radial)
  maybe_save(pB4_radial, "B4_radial_roe_bridge.png", w = 9, h = 9)
}

Phân tích Code:

• Tính toán mức độ đóng góp của từng thành phần DuPont vào sự thay đổi của ROE theo từng quý.
• Dùng pivot_longer để chuẩn bị dữ liệu cho việc vẽ biểu đồ.
• Tạo một biểu đồ cột xếp chồng tiêu chuẩn. Sau đó, Ta dùng coord_polar để biến đổi biểu đồ này thành một dạng tròn. Trục thời gian trở thành góc, và giá trị đóng góp trở thành bán kính.

Kết quả:

• Mỗi nan quạt là một quý. Thời gian chạy theo chiều kim đồng hồ.
• Nan quạt hướng ra ngoài cho thấy đóng góp dương vào sự thay đổi của ROE. Nan quạt hướng vào trong cho thấy đóng góp âm.
• Màu cam (Vòng quay tài sản) và màu xanh (Biên lợi nhuận ròng) là các nan quạt lớn nhất.
• Điều này có nghĩa là sự thay đổi trong ROE chủ yếu đến từ sự thay đổi trong hiệu quả sử dụng tài sản và khả năng sinh lời.
• Màu tím (Đòn bẩy tài chính) có đóng góp không đáng kể vào sự thay đổi này.