Kiểm định Fisher Exact – Monte Carlo

* Kiểm định Fisher Exact trong SPSS
* Mở dữ liệu
GET FILE='/data/vietlod.sav'.

* Tạo biến upuniv từ biến degree
RECODE degree (6=1) (1 thru 5=0) INTO upuniv.
VALUE LABELS upuniv 0 'Không có bằng sau đại học' 1 'Có bằng sau đại học'.

* Tạo biến agegroup (nếu chưa có)
RECODE age (LOWEST thru 25=1) (26 thru 35=2) (36 thru 45=3) 
           (46 thru 55=4) (56 thru HIGHEST=5) INTO agegroup.
VALUE LABELS agegroup 1 '≤25 tuổi' 2 '26-35 tuổi' 3 '36-45 tuổi' 
                     4 '46-55 tuổi' 5 '≥56 tuổi'.

* Kiểm tra phân phối dữ liệu
FREQUENCIES VARIABLES=upuniv agegroup.

* Thực hiện kiểm định chi bình phương thông thường
CROSSTABS
  /TABLES=agegroup BY upuniv
  /FORMAT=AVALUE TABLES
  /STATISTICS=CHISQ
  /CELLS=COUNT EXPECTED ROW COLUMN TOTAL
  /COUNT ROUND CELL.

* Thực hiện kiểm định Fisher Exact
CROSSTABS
  /TABLES=agegroup BY upuniv
  /FORMAT=AVALUE TABLES
  /STATISTICS=CHISQ
  /CELLS=COUNT EXPECTED ROW COLUMN TOTAL
  /COUNT ROUND CELL
  /EXACT FISHER.

* Kiểm định Fisher Exact với Monte Carlo (nếu cần)
CROSSTABS
  /TABLES=agegroup BY upuniv
  /FORMAT=AVALUE TABLES
  /STATISTICS=CHISQ
  /CELLS=COUNT EXPECTED ROW COLUMN TOTAL
  /COUNT ROUND CELL
  /EXACT FISHER MONTE CARLO CIN(95) SAMPLES(10000).

* Phân tích chi tiết cho từng nhóm tuổi
SPLIT FILE SEPARATE BY agegroup.
FREQUENCIES VARIABLES=upuniv.
SPLIT FILE OFF.

* Kiểm định Fisher Exact trong Stata
* Mở dữ liệu
use "vietlod.dta", clear

* Tạo biến upuniv từ biến degree
gen upuniv = (degree == 6)
label define upuniv_lab 0 "Không có bằng sau đại học" 1 "Có bằng sau đại học"
label values upuniv upuniv_lab

* Tạo biến agegroup (nếu chưa có)
gen agegroup = 1 if age <= 25 replace agegroup = 2 if age >= 26 & age <= 35 replace agegroup = 3 if age >= 36 & age <= 45 replace agegroup = 4 if age >= 46 & age <= 55 replace agegroup = 5 if age >= 56

label define age_lab 1 "≤25 tuổi" 2 "26-35 tuổi" 3 "36-45 tuổi" ///
                     4 "46-55 tuổi" 5 "≥56 tuổi"
label values agegroup age_lab

* Kiểm tra phân phối dữ liệu
tab upuniv
tab agegroup

* Tạo bảng chéo cơ bản
tab agegroup upuniv

* Kiểm định chi bình phương thông thường
tab agegroup upuniv, chi2

* Kiểm định Fisher Exact
tab agegroup upuniv, exact

* Kiểm định Fisher Exact với chi tiết
tab agegroup upuniv, exact all

* Tính toán exact test cho từng cặp so sánh
forvalues i = 1/5 {
    forvalues j = 1/1 {
        tab agegroup upuniv if agegroup == `i', exact
    }
}

* Vẽ đồ thị phân tích
graph bar (count), over(upuniv) over(agegroup) ///
    title("Phân bố trình độ học vấn theo nhóm tuổi") ///
    ytitle("Số lượng") ///
    legend(label(1 "Không có bằng sau ĐH") label(2 "Có bằng sau ĐH"))

* Tính tỷ lệ phần trăm
tab agegroup upuniv, row col

* Kiểm tra điều kiện sử dụng Fisher Exact
tab agegroup upuniv, expected

* Xuất kết quả ra file
outreg2 using "fisher_exact_results.xls", replace ///
    title("Kết quả kiểm định Fisher Exact") ///
    addstat("Chi-square", r(chi2), "Fisher Exact p-value", r(p_exact))

# Kiểm định Fisher Exact trong R
# Cài đặt và load các package cần thiết
if (!require(foreign)) install.packages("foreign")
if (!require(gmodels)) install.packages("gmodels")
if (!require(exact2x2)) install.packages("exact2x2")
if (!require(DescTools)) install.packages("DescTools")

library(foreign)
library(gmodels)
library(exact2x2)
library(DescTools)

# Đọc dữ liệu SPSS
data <- read.spss("vietlod.sav", to.data.frame = TRUE)

# Tạo biến upuniv từ biến degree
data$upuniv <- ifelse(data$degree == 6, 1, 0)
data$upuniv <- factor(data$upuniv, 
                     levels = c(0, 1),
                     labels = c("Không có bằng sau đại học", "Có bằng sau đại học"))

# Tạo biến agegroup (nếu chưa có)
data$agegroup <- cut(data$age, 
                    breaks = c(0, 25, 35, 45, 55, 100),
                    labels = c("≤25 tuổi", "26-35 tuổi", "36-45 tuổi", 
                              "46-55 tuổi", "≥56 tuổi"),
                    include.lowest = TRUE)

# Kiểm tra phân phối dữ liệu
table(data$upuniv)
table(data$agegroup)

# Tạo bảng chéo
contingency_table <- table(data$agegroup, data$upuniv)
print("Bảng chéo:")
print(contingency_table)

# Kiểm định chi bình phương thông thường
chi_square_test <- chisq.test(contingency_table)
print("Kết quả kiểm định chi bình phương:")
print(chi_square_test)

# Kiểm tra tần suất kỳ vọng
print("Tần suất kỳ vọng:")
print(chi_square_test$expected)

# Kiểm định Fisher Exact
fisher_test <- fisher.test(contingency_table)
print("Kết quả kiểm định Fisher Exact:")
print(fisher_test)

# Kiểm định Fisher Exact với simulate.p.value (tương tự Monte Carlo)
fisher_mc <- fisher.test(contingency_table, simulate.p.value = TRUE, B = 10000)
print("Kết quả kiểm định Fisher Exact (Monte Carlo):")
print(fisher_mc)

# Sử dụng package exact2x2 cho kiểm định chính xác hơn
if (ncol(contingency_table) == 2) {
    exact_test <- exact2x2::exact2x2(contingency_table)
    print("Kết quả kiểm định exact2x2:")
    print(exact_test)
}

# Phân tích chi tiết
CrossTable(data$agegroup, data$upuniv, 
           expected = TRUE, prop.r = TRUE, prop.c = TRUE, prop.t = TRUE,
           fisher = TRUE, chisq = TRUE, format = "SAS")

# Vẽ đồ thị
library(ggplot2)
library(dplyr)

# Tạo data frame cho visualization
plot_data <- data %>%
    group_by(agegroup, upuniv) %>%
    summarise(count = n(), .groups = "drop") %>%
    group_by(agegroup) %>%
    mutate(prop = count / sum(count))

# Đồ thị thanh
ggplot(plot_data, aes(x = agegroup, y = count, fill = upuniv)) +
    geom_bar(stat = "identity", position = "stack") +
    labs(title = "Phân bố trình độ học vấn theo nhóm tuổi",
         x = "Nhóm tuổi",
         y = "Số lượng",
         fill = "Trình độ học vấn") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Đồ thị tỷ lệ
ggplot(plot_data, aes(x = agegroup, y = prop, fill = upuniv)) +
    geom_bar(stat = "identity", position = "fill") +
    labs(title = "Tỷ lệ trình độ học vấn theo nhóm tuổi",
         x = "Nhóm tuổi",
         y = "Tỷ lệ",
         fill = "Trình độ học vấn") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    scale_y_continuous(labels = scales::percent)

# So sánh kết quả
comparison_results <- data.frame(
    Method = c("Chi-square", "Fisher Exact", "Fisher MC"),
    P_value = c(chi_square_test$p.value, fisher_test$p.value, fisher_mc$p.value),
    Statistic = c(chi_square_test$statistic, NA, NA)
)

print("So sánh kết quả các phương pháp:")
print(comparison_results)

# Kiểm tra điều kiện sử dụng Fisher Exact
min_expected <- min(chi_square_test$expected)
prop_cells_below_5 <- sum(chi_square_test$expected < 5) / length(chi_square_test$expected)

print(paste("Tần suất kỳ vọng nhỏ nhất:", round(min_expected, 2)))
print(paste("Tỷ lệ ô có tần suất kỳ vọng < 5:", round(prop_cells_below_5 * 100, 1), "%"))

if (min_expected < 5 | prop_cells_below_5 > 0.2) {
    print("Nên sử dụng Fisher Exact test")
} else {
    print("Có thể sử dụng Chi-square test")
}

# Kiểm định Fisher Exact trong Python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from scipy.stats import chi2_contingency, fisher_exact
import warnings
warnings.filterwarnings('ignore')

# Đọc dữ liệu
# Nếu từ file SPSS
try:
    data = pd.read_spss('vietlod.sav')
except:
    # Nếu từ file CSV
    data = pd.read_csv('vietlod.csv')

# Tạo biến upuniv từ biến degree
data['upuniv'] = (data['degree'] == 6).astype(int)
data['upuniv'] = data['upuniv'].map({0: 'Không có bằng sau đại học', 
                                    1: 'Có bằng sau đại học'})

# Tạo biến agegroup (nếu chưa có)
data['agegroup'] = pd.cut(data['age'], 
                         bins=[0, 25, 35, 45, 55, 100],
                         labels=['≤25 tuổi', '26-35 tuổi', '36-45 tuổi', 
                                '46-55 tuổi', '≥56 tuổi'],
                         include_lowest=True)

# Kiểm tra phân phối dữ liệu
print("Phân phối biến upuniv:")
print(data['upuniv'].value_counts())
print("\nPhân phối biến agegroup:")
print(data['agegroup'].value_counts())

# Tạo bảng chéo
contingency_table = pd.crosstab(data['agegroup'], data['upuniv'])
print("\nBảng chéo:")
print(contingency_table)

# Kiểm định chi bình phương thông thường
chi2, p_chi2, dof, expected = chi2_contingency(contingency_table)
print(f"\nKết quả kiểm định chi bình phương:")
print(f"Chi-square statistic: {chi2:.4f}")
print(f"P-value: {p_chi2:.6f}")
print(f"Degrees of freedom: {dof}")

print(f"\nTần suất kỳ vọng:")
expected_df = pd.DataFrame(expected, 
                          index=contingency_table.index,
                          columns=contingency_table.columns)
print(expected_df)

# Kiểm tra điều kiện sử dụng Fisher Exact
min_expected = expected.min()
prop_cells_below_5 = (expected < 5).sum() / expected.size

print(f"\nTần suất kỳ vọng nhỏ nhất: {min_expected:.2f}")
print(f"Tỷ lệ ô có tần suất kỳ vọng < 5: {prop_cells_below_5*100:.1f}%") # Kiểm định Fisher Exact # Lưu ý: scipy.stats.fisher_exact chỉ làm việc với bảng 2x2 # Với bảng lớn hơn, cần sử dụng thư viện khác def fisher_exact_2xk(table): """ Thực hiện Fisher Exact test cho bảng 2xk """ try: from scipy.stats import chi2_contingency # Sử dụng Monte Carlo simulation _, p_value, _, _ = chi2_contingency(table, lambda_=None) return p_value except: return None # Nếu bảng chỉ có 2 hàng, có thể sử dụng Fisher exact if contingency_table.shape[0] == 2: # Chuyển đổi thành bảng 2x2 nếu cần if contingency_table.shape[1] == 2: odds_ratio, p_fisher = fisher_exact(contingency_table) print(f"\nKết quả kiểm định Fisher Exact (2x2):") print(f"Odds ratio: {odds_ratio:.4f}") print(f"P-value: {p_fisher:.6f}") else: print(f"\nBảng {contingency_table.shape[0]}x{contingency_table.shape[1]} - cần phương pháp khác") # Sử dụng thư viện rpy2 để gọi R's fisher.test (nếu có) try: import rpy2.robjects as ro from rpy2.robjects.packages import importr from rpy2.robjects import pandas2ri pandas2ri.activate() stats_r = importr('stats') # Chuyển đổi bảng chéo sang R r_table = pandas2ri.py2rpy(contingency_table) # Thực hiện Fisher Exact test fisher_result = stats_r.fisher_test(r_table) p_fisher_r = fisher_result.rx2('p.value')[0] print(f"\nKết quả Fisher Exact (sử dụng R):") print(f"P-value: {p_fisher_r:.6f}") except ImportError: print("\nrpy2 không có sẵn - không thể sử dụng R's fisher.test") # Simulation-based Fisher Exact (Monte Carlo) def monte_carlo_fisher_exact(table, n_simulations=10000): """ Thực hiện Fisher Exact test bằng Monte Carlo simulation """ from scipy.stats import hypergeom # Tính toán p-value bằng simulation row_sums = table.sum(axis=1) col_sums = table.sum(axis=0) total = table.sum().sum() observed_stat = chi2_contingency(table)[0] count = 0 for _ in range(n_simulations): # Tạo bảng ngẫu nhiên với cùng marginal totals simulated_table = np.random.multivariate_hypergeometric( row_sums, col_sums.sum()) # Tính chi-square statistic sim_stat = chi2_contingency(simulated_table.reshape(table.shape))[0] if sim_stat >= observed_stat:
            count += 1
    
    return count / n_simulations

# Thực hiện Monte Carlo simulation
try:
    p_mc = monte_carlo_fisher_exact(contingency_table.values)
    print(f"\nKết quả Monte Carlo Fisher Exact:")
    print(f"P-value: {p_mc:.6f}")
except:
    print("\nKhông thể thực hiện Monte Carlo simulation")

# Vẽ đồ thị
plt.style.use('default')

# Đồ thị 1: Heatmap
plt.figure(figsize=(10, 6))
sns.heatmap(contingency_table, annot=True, fmt='d', cmap='Blues')
plt.title('Bảng chéo: Nhóm tuổi vs Trình độ học vấn')
plt.tight_layout()
plt.show()

# Đồ thị 2: Stacked bar chart
plt.figure(figsize=(12, 6))
contingency_table.plot(kind='bar', stacked=True, color=['lightblue', 'lightcoral'])
plt.title('Phân bố trình độ học vấn theo nhóm tuổi')
plt.xlabel('Nhóm tuổi')
plt.ylabel('Số lượng')
plt.legend(title='Trình độ học vấn')
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

# Đồ thị 3: Proportion chart
plt.figure(figsize=(12, 6))
prop_table = contingency_table.div(contingency_table.sum(axis=1), axis=0)
prop_table.plot(kind='bar', stacked=True, color=['lightblue', 'lightcoral'])
plt.title('Tỷ lệ trình độ học vấn theo nhóm tuổi')
plt.xlabel('Nhóm tuổi')
plt.ylabel('Tỷ lệ')
plt.legend(title='Trình độ học vấn')
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()

# So sánh kết quả
results_comparison = pd.DataFrame({
    'Method': ['Chi-square', 'Fisher Exact (if 2x2)', 'Monte Carlo'],
    'P_value': [p_chi2, 
               p_fisher if 'p_fisher' in locals() else np.nan,
               p_mc if 'p_mc' in locals() else np.nan],
    'Statistic': [chi2, 
                 odds_ratio if 'odds_ratio' in locals() else np.nan,
                 np.nan]
})

print("\nSo sánh kết quả các phương pháp:")
print(results_comparison)

# Kết luận
print("\n" + "="*50)
print("KẾT LUẬN VÀ KHUYẾN NGHỊ")
print("="*50)

if min_expected < 5 or prop_cells_below_5 > 0.2:
    print("✓ Nên sử dụng Fisher Exact test")
    print("✗ Không nên sử dụng Chi-square test")
    
    if 'p_fisher' in locals():
        conclusion = "có" if p_fisher < 0.05 else "không có"
        print(f"\nKết luận: {conclusion} mối quan hệ có ý nghĩa thống kê")
        print(f"P-value (Fisher): {p_fisher:.6f}")
    
else:
    print("✓ Có thể sử dụng Chi-square test")
    print("✓ Fisher Exact test cũng có thể sử dụng (chính xác hơn)")
    
    conclusion = "có" if p_chi2 < 0.05 else "không có"
    print(f"\nKết luận: {conclusion} mối quan hệ có ý nghĩa thống kê")
    print(f"P-value (Chi-square): {p_chi2:.6f}")