Kiểm định Chi bình phương

* Kiểm định chi bình phương trong SPSS
* Tác giả: Hướng dẫn thực hành kinh tế lượng
* Mục đích: Phân tích mối quan hệ giữa giới tính và loại hình doanh nghiệp

* Mở file dữ liệu
GET FILE='/data/vietlod.sav'.

* Thực hiện kiểm định chi bình phương
CROSSTABS
  /TABLES=gender BY section
  /FORMAT=AVALUE TABLES
  /STATISTICS=CHISQ PHI CC
  /CELLS=COUNT EXPECTED ROW COLUMN TOTAL
  /COUNT ROUND CELL
  /BARCHART.

* Thống kê mô tả cho các biến
FREQUENCIES VARIABLES=gender section
  /ORDER=ANALYSIS.

* Tạo bảng chéo chi tiết hơn
CROSSTABS
  /TABLES=gender BY section
  /FORMAT=AVALUE TABLES
  /STATISTICS=CHISQ PHI CC LAMBDA CORR
  /CELLS=COUNT EXPECTED ROW COLUMN TOTAL RESID SRESID
  /COUNT ROUND CELL.

* Xuất kết quả ra file
OUTPUT EXPORT
  /CONTENTS EXPORT=VISIBLE LAYERS=PRINTSETTING MODELVIEWS=PRINTSETTING
  /HTML DOCUMENTFILE='chi_square_results.html'.

* Kiểm định chi bình phương trong Stata
* Tác giả: Hướng dẫn thực hành kinh tế lượng
* Mục đích: Phân tích mối quan hệ giữa giới tính và loại hình doanh nghiệp

* Mở file dữ liệu
use "/data/vietlod.dta", clear

* Thống kê mô tả
tab gender, missing
tab section, missing
tab gender section, missing

* Tạo bảng chéo với thống kê chi bình phương
tab gender section, chi2 V exact

* Tạo bảng chéo chi tiết với tần suất kỳ vọng
tab gender section, chi2 expected

* Tạo bảng chéo với phần trăm
tab gender section, row col cell

* Kiểm định chi bình phương với correction
tab gender section, chi2 exact

* Tính Cramer's V
tab gender section, V

* Vẽ biểu đồ thanh
graph bar, over(section) over(gender) ///
    title("Phân bố giới tính theo loại hình doanh nghiệp") ///
    ylabel(, angle(0)) ///
    legend(title("Giới tính"))

* Lưu biểu đồ
graph export "gender_section_bar.png", replace

* Kiểm tra điều kiện áp dụng
tab gender section, expected chi2
* Kiểm tra xem có ô nào có tần suất kỳ vọng < 5 không

* Tính các thống kê bổ sung
tab gender section, gamma taub

# Kiểm định chi bình phương trong R
# Tác giả: Hướng dẫn thực hành kinh tế lượng
# Mục đích: Phân tích mối quan hệ giữa giới tính và loại hình doanh nghiệp

# Cài đặt và load các package cần thiết
library(haven)
library(dplyr)
library(ggplot2)
library(vcd)
library(corrplot)
library(gmodels)

# Đọc dữ liệu
data <- read_sav("/data/vietlod.sav")

# Kiểm tra cấu trúc dữ liệu
str(data)
head(data[c("gender", "section")])

# Thống kê mô tả
table(data$gender, useNA = "always")
table(data$section, useNA = "always")

# Tạo bảng chéo
cross_table <- table(data$gender, data$section)
print(cross_table)

# Thực hiện kiểm định chi bình phương
chi_test <- chisq.test(cross_table)
print(chi_test)

# Tần suất kỳ vọng
print("Tần suất kỳ vọng:")
print(chi_test$expected)

# Residuals
print("Residuals:")
print(chi_test$residuals)

# Standardized residuals
print("Standardized residuals:")
print(chi_test$stdres)

# Tính Cramer's V
cramers_v <- sqrt(chi_test$statistic / (sum(cross_table) * (min(dim(cross_table)) - 1)))
print(paste("Cramer's V:", round(cramers_v, 3)))

# Tạo bảng chéo chi tiết với CrossTable
CrossTable(data$gender, data$section, 
           expected = TRUE, 
           prop.r = TRUE, 
           prop.c = TRUE, 
           prop.t = TRUE, 
           prop.chisq = FALSE,
           chisq = TRUE,
           fisher = TRUE,
           format = "SPSS")

# Vẽ biểu đồ thanh
# Chuyển đổi dữ liệu
data_plot <- data %>%
  count(gender, section) %>%
  mutate(
    gender = factor(gender, levels = c(1, 2), labels = c("Nam", "Nữ")),
    section = factor(section, levels = c(1, 2, 3, 4), 
                    labels = c("DNNN", "DNTN", "FDI", "HTX-hộ cá thể"))
  )

# Tạo biểu đồ thanh
ggplot(data_plot, aes(x = section, y = n, fill = gender)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Phân bố giới tính theo loại hình doanh nghiệp",
       x = "Loại hình doanh nghiệp",
       y = "Số lượng",
       fill = "Giới tính") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Lưu biểu đồ
ggsave("gender_section_bar.png", width = 10, height = 6, dpi = 300)

# Vẽ heatmap
library(reshape2)
cross_table_pct <- prop.table(cross_table, margin = 1) * 100
cross_table_melt <- melt(cross_table_pct)

ggplot(cross_table_melt, aes(Var2, Var1, fill = value)) +
  geom_tile() +
  scale_fill_gradient(low = "white", high = "steelblue") +
  labs(title = "Phần trăm theo hàng - Giới tính vs Loại hình doanh nghiệp",
       x = "Loại hình doanh nghiệp",
       y = "Giới tính",
       fill = "Phần trăm (%)") +
  theme_minimal()

# Kiểm tra điều kiện áp dụng
min_expected <- min(chi_test$expected) print(paste("Tần suất kỳ vọng nhỏ nhất:", round(min_expected, 2))) if(min_expected >= 5) {
  print("Điều kiện tần suất kỳ vọng ≥ 5 được thỏa mãn")
} else {
  print("Cảnh báo: Có ô có tần suất kỳ vọng < 5")
}

# Kết luận
alpha <- 0.05
if(chi_test$p.value < alpha) {
  print(paste("Kết luận: Với α =", alpha, ", p-value =", round(chi_test$p.value, 6)))
  print("Bác bỏ H0. Có mối quan hệ giữa giới tính và loại hình doanh nghiệp")
} else {
  print(paste("Kết luận: Với α =", alpha, ", p-value =", round(chi_test$p.value, 6)))
  print("Chưa có đủ bằng chứng để bác bỏ H0")
}

# Kiểm định chi bình phương trong Python
# Tác giả: Hướng dẫn thực hành kinh tế lượng
# Mục đích: Phân tích mối quan hệ giữa giới tính và loại hình doanh nghiệp

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import chi2_contingency
from scipy.stats import chi2
import pyreadstat

# Đọc dữ liệu SPSS
data, meta = pyreadstat.read_sav('/data/vietlod.sav')

# Kiểm tra cấu trúc dữ liệu
print(data.info())
print(data[['gender', 'section']].head())

# Thống kê mô tả
print("Thống kê mô tả:")
print(data['gender'].value_counts())
print(data['section'].value_counts())

# Tạo bảng chéo
cross_table = pd.crosstab(data['gender'], data['section'], margins=True)
print("\nBảng chéo:")
print(cross_table)

# Thực hiện kiểm định chi bình phương
chi2_stat, p_value, dof, expected = chi2_contingency(cross_table.iloc[:-1, :-1])

print(f"\nKết quả kiểm định chi bình phương:")
print(f"Chi-square statistic: {chi2_stat:.3f}")
print(f"P-value: {p_value:.6f}")
print(f"Degrees of freedom: {dof}")

# Tần suất kỳ vọng
print("\nTần suất kỳ vọng:")
expected_df = pd.DataFrame(expected, 
                          index=cross_table.index[:-1], 
                          columns=cross_table.columns[:-1])
print(expected_df)

# Tính Cramer's V
n = cross_table.iloc[:-1, :-1].sum().sum()
cramers_v = np.sqrt(chi2_stat / (n * (min(cross_table.shape) - 2)))
print(f"\nCramer's V: {cramers_v:.3f}")

# Tạo bảng chéo với phần trăm
print("\nBảng chéo với phần trăm:")
cross_table_pct = pd.crosstab(data['gender'], data['section'], 
                             normalize='index') * 100
print(cross_table_pct.round(1))

# Vẽ biểu đồ thanh
plt.figure(figsize=(12, 8))

# Chuẩn bị dữ liệu cho biểu đồ
data_plot = data.copy()
data_plot['gender'] = data_plot['gender'].map({1: 'Nam', 2: 'Nữ'})
data_plot['section'] = data_plot['section'].map({
    1: 'DNNN', 2: 'DNTN', 3: 'FDI', 4: 'HTX-hộ cá thể'
})

# Tạo biểu đồ thanh
plt.subplot(2, 2, 1)
sns.countplot(data=data_plot, x='section', hue='gender')
plt.title('Phân bố giới tính theo loại hình doanh nghiệp')
plt.xlabel('Loại hình doanh nghiệp')
plt.ylabel('Số lượng')
plt.xticks(rotation=45)
plt.legend(title='Giới tính')

# Tạo heatmap
plt.subplot(2, 2, 2)
sns.heatmap(cross_table_pct, annot=True, fmt='.1f', cmap='Blues')
plt.title('Phần trăm theo hàng (%)')
plt.xlabel('Loại hình doanh nghiệp')
plt.ylabel('Giới tính')

# Biểu đồ residuals
residuals = (cross_table.iloc[:-1, :-1] - expected) / np.sqrt(expected)
plt.subplot(2, 2, 3)
sns.heatmap(residuals, annot=True, fmt='.2f', cmap='RdBu_r', center=0)
plt.title('Standardized Residuals')
plt.xlabel('Loại hình doanh nghiệp')
plt.ylabel('Giới tính')

# Biểu đồ tần suất kỳ vọng vs quan sát
plt.subplot(2, 2, 4)
observed = cross_table.iloc[:-1, :-1].values.flatten()
expected_flat = expected.flatten()
plt.scatter(expected_flat, observed, alpha=0.7)
plt.plot([0, max(expected_flat)], [0, max(expected_flat)], 'r--', label='y=x')
plt.xlabel('Tần suất kỳ vọng')
plt.ylabel('Tần suất quan sát')
plt.title('Tần suất kỳ vọng vs Quan sát')
plt.legend()

plt.tight_layout()
plt.show()

# Kiểm tra điều kiện áp dụng
min_expected = expected.min()
print(f"\nTần suất kỳ vọng nhỏ nhất: {min_expected:.2f}")

if min_expected >= 5:
    print("Điều kiện tần suất kỳ vọng ≥ 5 được thỏa mãn")
else:
    print("Cảnh báo: Có ô có tần suất kỳ vọng < 5")

# Kết luận
alpha = 0.05
if p_value < alpha:
    print(f"\nKết luận: Với α = {alpha}, p-value = {p_value:.6f}")
    print("Bác bỏ H₀. Có mối quan hệ giữa giới tính và loại hình doanh nghiệp")
    print(f"Độ mạnh mối quan hệ (Cramer's V): {cramers_v:.3f}")
else:
    print(f"\nKết luận: Với α = {alpha}, p-value = {p_value:.6f}")
    print("Chưa có đủ bằng chứng để bác bỏ H₀")

# Tạo báo cáo tóm tắt
print("\n" + "="*50)
print("TÓM TẮT KẾT QUẢ KIỂM ĐỊNH CHI BÌNH PHƯƠNG")
print("="*50)
print(f"Cỡ mẫu: {n}")
print(f"Chi-square statistic: {chi2_stat:.3f}")
print(f"Degrees of freedom: {dof}")
print(f"P-value: {p_value:.6f}")
print(f"Cramer's V: {cramers_v:.3f}")
print(f"Tần suất kỳ vọng nhỏ nhất: {min_expected:.2f}")
print("\nKết luận:")
if p_value < alpha:
    print("- Có mối quan hệ có ý nghĩa thống kê")
    if cramers_v < 0.3:
        strength = "yếu"
    elif cramers_v < 0.5:
        strength = "trung bình"
    else:
        strength = "mạnh"
    print(f"- Độ mạnh mối quan hệ: {strength}")
else:
    print("- Không có mối quan hệ có ý nghĩa thống kê")