Kiểm định Kruskal-Wallis

* Kiểm định Kruskal-Wallis trong SPSS Syntax
* Tác giả: Thực hành kinh tế lượng
* Ngày: 2025

* Mở file dữ liệu
GET FILE='/data/2-way-anova.sav'.

* Mô tả dữ liệu theo nhóm
DESCRIPTIVES VARIABLES=earn
  /STATISTICS=MEAN STDDEV MIN MAX
  /SORT=MEAN (A).

* Kiểm tra phân phối bằng histogram
GRAPH
  /HISTOGRAM(NORMAL)=earn
  /PANEL ROWVAR=structure ROWOP=CROSS.

* Thực hiện kiểm định Kruskal-Wallis
NPAR TESTS
  /K-W=earn BY structure(1 3)
  /STATISTICS DESCRIPTIVES QUARTILES
  /MISSING ANALYSIS.

* Thống kê mô tả chi tiết theo nhóm
MEANS TABLES=earn BY structure
  /CELLS=MEAN COUNT STDDEV MEDIAN
  /STATISTICS ANOVA.

* Kiểm định post-hoc (nếu cần)
* Sử dụng Mann-Whitney U cho từng cặp
NPAR TESTS
  /M-W=earn BY structure(1 2)
  /MISSING ANALYSIS.

NPAR TESTS
  /M-W=earn BY structure(1 3)
  /MISSING ANALYSIS.

NPAR TESTS
  /M-W=earn BY structure(2 3)
  /MISSING ANALYSIS.

* Xuất kết quả
DISPLAY "Kiểm định Kruskal-Wallis hoàn tất".

* Kiểm định Kruskal-Wallis trong Stata
* Tác giả: Thực hành kinh tế lượng
* Ngày: 2025

* Mở file dữ liệu (chuyển đổi từ SPSS)
use "kruskal_wallis.dta", clear

* Mô tả dữ liệu
describe earn structure
summarize earn, detail

* Thống kê mô tả theo nhóm
by structure, sort: summarize earn, detail

* Kiểm tra phân phối bằng histogram
histogram earn, by(structure) normal ///
    title("Phân phối thu nhập theo lĩnh vực") ///
    xtitle("Thu nhập (triệu đồng)") ///
    ytitle("Mật độ")

* Boxplot để so sánh
graph box earn, over(structure) ///
    title("So sánh thu nhập theo lĩnh vực") ///
    ytitle("Thu nhập (triệu đồng)")

* Thực hiện kiểm định Kruskal-Wallis
kwallis earn, by(structure)

* Kiểm định post-hoc bằng Dunn's test
* Cần cài đặt package dunntest
* ssc install dunntest
dunntest earn, by(structure) ma(bonferroni)

* Thống kê chi tiết
tabstat earn, by(structure) statistics(n mean sd min max median q)

* Kiểm tra tính đồng nhất phương sai (để so sánh)
robvar earn, by(structure)

* Xuất kết quả
display "Kiểm định Kruskal-Wallis hoàn tất"
display "Mean ranks được sử dụng để so sánh các nhóm"

# Kiểm định Kruskal-Wallis trong R
# Tác giả: Thực hành kinh tế lượng
# Ngày: 2025

# Cài đặt và load các package cần thiết
if (!require(foreign)) install.packages("foreign")
if (!require(dplyr)) install.packages("dplyr")
if (!require(ggplot2)) install.packages("ggplot2")
if (!require(FSA)) install.packages("FSA")
if (!require(dunn.test)) install.packages("dunn.test")

library(foreign)
library(dplyr)
library(ggplot2)
library(FSA)
library(dunn.test)

# Đọc dữ liệu từ file SPSS
data <- read.spss("2-way-anova.sav", to.data.frame = TRUE)

# Kiểm tra cấu trúc dữ liệu
str(data)
head(data)

# Thống kê mô tả theo nhóm
summary_stats <- data %>%
  group_by(structure) %>%
  summarise(
    n = n(),
    mean = mean(earn, na.rm = TRUE),
    median = median(earn, na.rm = TRUE),
    sd = sd(earn, na.rm = TRUE),
    min = min(earn, na.rm = TRUE),
    max = max(earn, na.rm = TRUE),
    q1 = quantile(earn, 0.25, na.rm = TRUE),
    q3 = quantile(earn, 0.75, na.rm = TRUE)
  )

print("Thống kê mô tả theo lĩnh vực:")
print(summary_stats)

# Vẽ boxplot
ggplot(data, aes(x = factor(structure), y = earn)) +
  geom_boxplot(fill = "lightblue", alpha = 0.7) +
  labs(title = "Phân phối thu nhập theo lĩnh vực kinh tế",
       x = "Lĩnh vực (1=Dịch vụ, 2=Công nghiệp, 3=Nông nghiệp)",
       y = "Thu nhập (triệu đồng)") +
  theme_minimal()

# Vẽ histogram
ggplot(data, aes(x = earn)) +
  geom_histogram(bins = 30, fill = "skyblue", alpha = 0.7) +
  facet_wrap(~structure, ncol = 3) +
  labs(title = "Histogram thu nhập theo lĩnh vực",
       x = "Thu nhập (triệu đồng)",
       y = "Tần số") +
  theme_minimal()

# Thực hiện kiểm định Kruskal-Wallis
kruskal_result <- kruskal.test(earn ~ structure, data = data)
print("Kết quả kiểm định Kruskal-Wallis:")
print(kruskal_result)

# Kiểm định post-hoc bằng Dunn's test
dunn_result <- dunn.test(data$earn, data$structure, 
                        method = "bonferroni", alpha = 0.05)
print("Kết quả kiểm định post-hoc (Dunn's test):")

# Tính effect size (eta-squared)
eta_squared <- kruskal_result$statistic / (nrow(data) - 1)
cat("Effect size (eta-squared):", eta_squared, "\n")

# Kiểm định Levene để kiểm tra homogeneity of variance
levene_result <- car::leveneTest(earn ~ factor(structure), data = data)
print("Kiểm định Levene (homogeneity of variance):")
print(levene_result)

# Tạo bảng kết quả đẹp
results_table <- data.frame(
  Statistic = c("Chi-square", "df", "p-value", "Effect size"),
  Value = c(kruskal_result$statistic, 
           kruskal_result$parameter, 
           kruskal_result$p.value,
           eta_squared)
)

print("Bảng tóm tắt kết quả:")
print(results_table)

# Kết luận
alpha <- 0.05
if (kruskal_result$p.value < alpha) {
  cat("\nKết luận: Với p-value =", kruskal_result$p.value, 
      "< α =", alpha, "\n") cat("Bác bỏ giả thuyết H0. Có sự khác biệt có ý nghĩa thống kê.\n") } else { cat("\nKết luận: Với p-value =", kruskal_result$p.value, ">= α =", alpha, "\n")
  cat("Không bác bỏ giả thuyết H0. Không có sự khác biệt có ý nghĩa thống kê.\n")
}

# Kiểm định Kruskal-Wallis trong Python
# Tác giả: Thực hành kinh tế lượng
# Ngày: 2025

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import scikit_posthocs as sp

# Đọc dữ liệu từ file SPSS
try:
    data = pd.read_spss('2-way-anova.sav')
except:
    # Nếu không có pyreadstat, tạo dữ liệu mẫu
    print("Tạo dữ liệu mẫu...")
    np.random.seed(42)
    n = 7287
    structure = np.random.choice([1, 2, 3], n, p=[0.4, 0.35, 0.25])
    
    # Tạo thu nhập với sự khác biệt giữa các nhóm
    earn = np.zeros(n)
    for i in range(n):
        if structure[i] == 1:  # Dịch vụ
            earn[i] = np.random.gamma(2, 20) + 10
        elif structure[i] == 2:  # Công nghiệp
            earn[i] = np.random.gamma(2, 15) + 5
        else:  # Nông nghiệp
            earn[i] = np.random.gamma(2, 10) + 2
    
    data = pd.DataFrame({
        'earn': earn,
        'structure': structure
    })

# Kiểm tra dữ liệu
print("Thông tin dữ liệu:")
print(data.info())
print("\nMẫu dữ liệu:")
print(data.head(10))

# Thống kê mô tả theo nhóm
print("\nThống kê mô tả theo lĩnh vực:")
summary_stats = data.groupby('structure')['earn'].agg([
    'count', 'mean', 'median', 'std', 'min', 'max',
    lambda x: x.quantile(0.25),  # Q1
    lambda x: x.quantile(0.75)   # Q3
]).round(2)

summary_stats.columns = ['n', 'mean', 'median', 'std', 'min', 'max', 'Q1', 'Q3']
print(summary_stats)

# Vẽ biểu đồ
fig, axes = plt.subplots(2, 2, figsize=(15, 10))

# Boxplot
axes[0, 0].boxplot([data[data['structure'] == i]['earn'] for i in [1, 2, 3]], 
                   labels=['Dịch vụ', 'Công nghiệp', 'Nông nghiệp'])
axes[0, 0].set_title('Boxplot thu nhập theo lĩnh vực')
axes[0, 0].set_ylabel('Thu nhập (triệu đồng)')

# Histogram
for i, sector in enumerate(['Dịch vụ', 'Công nghiệp', 'Nông nghiệp']):
    sector_data = data[data['structure'] == i+1]['earn']
    axes[0, 1].hist(sector_data, alpha=0.7, label=sector, bins=30)
axes[0, 1].set_title('Histogram thu nhập theo lĩnh vực')
axes[0, 1].set_xlabel('Thu nhập (triệu đồng)')
axes[0, 1].set_ylabel('Tần số')
axes[0, 1].legend()

# Violin plot
data_for_violin = [data[data['structure'] == i]['earn'] for i in [1, 2, 3]]
axes[1, 0].violinplot(data_for_violin, positions=[1, 2, 3])
axes[1, 0].set_title('Violin plot thu nhập theo lĩnh vực')
axes[1, 0].set_ylabel('Thu nhập (triệu đồng)')
axes[1, 0].set_xticks([1, 2, 3])
axes[1, 0].set_xticklabels(['Dịch vụ', 'Công nghiệp', 'Nông nghiệp'])

# QQ plot để kiểm tra phân phối chuẩn
stats.probplot(data['earn'], dist="norm", plot=axes[1, 1])
axes[1, 1].set_title('Q-Q Plot để kiểm tra phân phối chuẩn')

plt.tight_layout()
plt.show()

# Kiểm tra giả thiết về phân phối chuẩn
print("\nKiểm định Shapiro-Wilk cho từng nhóm:")
for i in [1, 2, 3]:
    sector_data = data[data['structure'] == i]['earn']
    if len(sector_data) > 3:  # Shapiro-Wilk cần ít nhất 3 quan sát
        shapiro_stat, shapiro_p = stats.shapiro(sector_data[:5000])  # Giới hạn mẫu
        print(f"Lĩnh vực {i}: Shapiro-Wilk statistic = {shapiro_stat:.4f}, p-value = {shapiro_p:.4f}")

# Thực hiện kiểm định Kruskal-Wallis
group1 = data[data['structure'] == 1]['earn']
group2 = data[data['structure'] == 2]['earn']
group3 = data[data['structure'] == 3]['earn']

kruskal_stat, kruskal_p = stats.kruskal(group1, group2, group3)

print(f"\nKết quả kiểm định Kruskal-Wallis:")
print(f"Test statistic (H): {kruskal_stat:.6f}")
print(f"P-value: {kruskal_p:.6f}")
print(f"Degrees of freedom: {len([group1, group2, group3]) - 1}")

# Tính effect size (eta-squared)
n_total = len(data)
eta_squared = (kruskal_stat - len([group1, group2, group3]) + 1) / (n_total - len([group1, group2, group3]))
print(f"Effect size (eta-squared): {eta_squared:.6f}")

# Kiểm định post-hoc bằng Dunn's test
print("\nKiểm định post-hoc (Dunn's test):")
dunn_results = sp.posthoc_dunn(data, val_col='earn', group_col='structure', 
                               p_adjust='bonferroni')
print(dunn_results)

# Tính mean ranks
from scipy.stats import rankdata
ranks = rankdata(data['earn'])
data['ranks'] = ranks

mean_ranks = data.groupby('structure')['ranks'].mean()
print("\nMean ranks:")
for i in [1, 2, 3]:
    print(f"Lĩnh vực {i}: {mean_ranks[i]:.2f}")

# Tạo bảng kết quả tóm tắt
results_summary = pd.DataFrame({
    'Lĩnh vực': ['Dịch vụ', 'Công nghiệp', 'Nông nghiệp'],
    'N': [len(group1), len(group2), len(group3)],
    'Mean': [group1.mean(), group2.mean(), group3.mean()],
    'Median': [group1.median(), group2.median(), group3.median()],
    'Mean Rank': [mean_ranks[1], mean_ranks[2], mean_ranks[3]]
})

print("\nBảng tóm tắt kết quả:")
print(results_summary.round(2))

# Kết luận
alpha = 0.05
if kruskal_p < alpha:
    print(f"\nKết luận: Với p-value = {kruskal_p:.6f} < α = {alpha}") print("Bác bỏ giả thuyết H0. Có sự khác biệt có ý nghĩa thống kê về thu nhập giữa các lĩnh vực.") else: print(f"\nKết luận: Với p-value = {kruskal_p:.6f} >= α = {alpha}")
    print("Không bác bỏ giả thuyết H0. Không có sự khác biệt có ý nghĩa thống kê về thu nhập giữa các lĩnh vực.")

# Lưu kết quả
results_dict = {
    'kruskal_statistic': kruskal_stat,
    'p_value': kruskal_p,
    'effect_size': eta_squared,
    'mean_ranks': mean_ranks.to_dict(),
    'summary_stats': summary_stats.to_dict()
}

print("\nKiểm định Kruskal-Wallis hoàn tất!")