Kiểm định Wilcoxon-Mann-Whitney

* Kiểm định Wilcoxon-Mann-Whitney trong SPSS
* Tác giả: Thực hành kinh tế lượng
* Mục đích: So sánh thu nhập theo giới tính

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

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

* Thực hiện kiểm định Wilcoxon-Mann-Whitney
NPAR TESTS
  /M-W= earn BY gender(0 1)
  /STATISTICS DESCRIPTIVES QUARTILES
  /MISSING ANALYSIS.

* Thống kê mô tả theo nhóm
EXAMINE VARIABLES=earn BY gender
  /PLOT BOXPLOT HISTOGRAM
  /STATISTICS DESCRIPTIVES
  /CINTERVAL 95
  /MISSING LISTWISE
  /NOTOTAL.

* Kiểm định Wilcoxon-Mann-Whitney trong Stata
* Tác giả: Thực hành kinh tế lượng
* Mục đích: So sánh thu nhập theo giới tính

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

* Kiểm tra phân phối bằng histogram
histogram earn, by(gender) normal ///
    title("Phân phối thu nhập theo giới tính") ///
    subtitle("Dữ liệu lao động Việt Nam") ///
    note("Nguồn: Khảo sát lao động quốc gia")

* Thống kê mô tả theo nhóm
tabstat earn, by(gender) statistics(n mean median sd min max p25 p75) ///
    columns(statistics) format(%9.2f)

* Thực hiện kiểm định Wilcoxon-Mann-Whitney (ranksum)
ranksum earn, by(gender)

* Kiểm định với các tùy chọn bổ sung
ranksum earn, by(gender) porder

* Tạo biểu đồ box plot để so sánh
graph box earn, over(gender) ///
    title("So sánh thu nhập theo giới tính") ///
    subtitle("Kiểm định Wilcoxon-Mann-Whitney") ///
    ytitle("Thu nhập (triệu đồng/năm)") ///
    note("Nguồn: Dữ liệu lao động Việt Nam")

# Kiểm định Wilcoxon-Mann-Whitney trong R
# Tác giả: Thực hành kinh tế lượng
# Mục đích: So sánh thu nhập theo giới tính

# Tải thư viện cần thiết
library(haven)      # Đọc file SPSS
library(ggplot2)    # Vẽ biểu đồ
library(dplyr)      # Xử lý dữ liệu
library(coin)       # Kiểm định phi tham số

# Đọc dữ liệu từ file SPSS
data <- read_sav("vietlod.sav") # Kiểm tra cấu trúc dữ liệu str(data) summary(data[c("earn", "gender")]) # Thống kê mô tả theo nhóm data %>%
  group_by(gender) %>%
  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),
    q25 = quantile(earn, 0.25, na.rm = TRUE),
    q75 = quantile(earn, 0.75, na.rm = TRUE)
  )

# Tạo histogram so sánh
ggplot(data, aes(x = earn, fill = factor(gender))) +
  geom_histogram(alpha = 0.7, position = "identity", bins = 30) +
  facet_wrap(~gender, labeller = labeller(gender = c("0" = "Nữ", "1" = "Nam"))) +
  labs(
    title = "Phân phối thu nhập theo giới tính",
    x = "Thu nhập (triệu đồng/năm)",
    y = "Tần suất",
    fill = "Giới tính"
  ) +
  theme_minimal()

# Thực hiện kiểm định Wilcoxon-Mann-Whitney
# Phương pháp 1: Sử dụng hàm wilcox.test()
result1 <- wilcox.test(earn ~ gender, data = data, 
                       alternative = "two.sided")
print(result1)

# Phương pháp 2: Sử dụng thư viện coin (chính xác hơn)
result2 <- wilcox_test(earn ~ factor(gender), data = data,
                       distribution = "exact")
print(result2)

# Tạo box plot so sánh
ggplot(data, aes(x = factor(gender), y = earn)) +
  geom_boxplot() +
  labs(
    title = "So sánh thu nhập theo giới tính",
    x = "Giới tính (0 = Nữ, 1 = Nam)",
    y = "Thu nhập (triệu đồng/năm)"
  ) +
  theme_minimal()

# Tính effect size (r)
z_score <- qnorm(result1$p.value/2)
n <- nrow(data)
effect_size <- abs(z_score) / sqrt(n)
cat("Effect size (r):", effect_size, "\n")

# Kiểm định Wilcoxon-Mann-Whitney trong Python
# Tác giả: Thực hành kinh tế lượng
# Mục đích: So sánh thu nhập theo giới tính

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

# Thiết lập hiển thị tiếng Việt
plt.rcParams['font.family'] = 'DejaVu Sans'
sns.set_style("whitegrid")

# Đọc dữ liệu từ file SPSS
data, meta = pyreadstat.read_sav("vietlod.sav")

# Kiểm tra thông tin cơ bản
print("Thông tin dữ liệu:")
print(data[['earn', 'gender']].info())
print("\nThống kê mô tả:")
print(data[['earn', 'gender']].describe())

# Thống kê mô tả theo nhóm
print("\nThống kê theo giới tính:")
stats_by_gender = data.groupby('gender')['earn'].agg([
    'count', 'mean', 'median', 'std', 'min', 'max',
    lambda x: x.quantile(0.25),
    lambda x: x.quantile(0.75)
]).round(2)
stats_by_gender.columns = ['Số lượng', 'Trung bình', 'Trung vị', 'Độ lệch chuẩn', 
                          'Min', 'Max', 'Q25', 'Q75']
print(stats_by_gender)

# Tạo histogram so sánh
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
data[data['gender'] == 0]['earn'].hist(alpha=0.7, bins=30, label='Nữ', color='pink')
data[data['gender'] == 1]['earn'].hist(alpha=0.7, bins=30, label='Nam', color='lightblue')
plt.xlabel('Thu nhập (triệu đồng/năm)')
plt.ylabel('Tần suất')
plt.title('Phân phối thu nhập theo giới tính')
plt.legend()

plt.subplot(1, 2, 2)
sns.boxplot(data=data, x='gender', y='earn')
plt.xlabel('Giới tính (0 = Nữ, 1 = Nam)')
plt.ylabel('Thu nhập (triệu đồng/năm)')
plt.title('Box plot thu nhập theo giới tính')

plt.tight_layout()
plt.show()

# Tách dữ liệu theo nhóm
earn_female = data[data['gender'] == 0]['earn'].dropna()
earn_male = data[data['gender'] == 1]['earn'].dropna()

# Thực hiện kiểm định Wilcoxon-Mann-Whitney
statistic, p_value = stats.mannwhitneyu(earn_male, earn_female, 
                                       alternative='two-sided')

print(f"\nKết quả kiểm định Wilcoxon-Mann-Whitney:")
print(f"Thống kê U: {statistic}")
print(f"P-value: {p_value:.6f}")

# Tính z-score
n1, n2 = len(earn_male), len(earn_female)
mean_u = n1 * n2 / 2
std_u = np.sqrt(n1 * n2 * (n1 + n2 + 1) / 12)
z_score = (statistic - mean_u) / std_u

print(f"Z-score: {z_score:.3f}")

# Diễn giải kết quả
alpha = 0.05
if p_value < alpha:
    print(f"\nKết luận: Với α = {alpha}, p-value = {p_value:.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 theo giới tính.") else: print(f"\nKết luận: Với α = {alpha}, p-value = {p_value:.6f} >= {alpha}")
    print("Không đủ bằng chứng để bác bỏ giả thuyết H0.")

# Tính effect size
effect_size = abs(z_score) / np.sqrt(n1 + n2)
print(f"Effect size (r): {effect_size:.4f}")

# Diễn giải effect size
if effect_size < 0.1:
    effect_interpretation = "nhỏ"
elif effect_size < 0.3:
    effect_interpretation = "trung bình"
else:
    effect_interpretation = "lớn"
print(f"Độ lớn hiệu ứng: {effect_interpretation}")