Kiểm định dấu hạng Wilcoxon

* Kiểm định McNemar 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/mcnemar.sav'.

* Thực hiện kiểm định McNemar
NPAR TESTS
  /MCNEMAR = Before WITH After (PAIRED)
  /MISSING ANALYSIS
  /STATISTICS DESCRIPTIVES.

* Tạo bảng chéo để xem phân bố dữ liệu
CROSSTABS
  /TABLES=Before BY After
  /FORMAT=AVALUE TABLES
  /STATISTICS=MCNEMAR
  /CELLS=COUNT EXPECTED ROW COLUMN TOTAL.

* Kiểm định McNemar 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 sang Stata format)
use "mcnemar.dta", clear

* Mô tả dữ liệu
describe Before After
summarize Before After

* Tạo bảng chéo
tabulate Before After

* Thực hiện kiểm định McNemar
mcci a b c d
* Trong đó: a, b, c, d là các giá trị từ bảng chéo 2x2

* Hoặc sử dụng lệnh symmetry cho kiểm định McNemar
symmetry Before After

* Kiểm định chính xác cho mẫu nhỏ
mcc Before After, exact

* Xuất kết quả
display "Kiểm định McNemar hoàn tất"
display "Kiểm tra p-value để đưa ra kết luận"

# Kiểm định McNemar 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(exact2x2)) install.packages("exact2x2")
library(foreign)
library(exact2x2)

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

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

# Tạo bảng chéo
table_cross <- table(data$Before, data$After)
print("Bảng chéo:")
print(table_cross)

# Thực hiện kiểm định McNemar
mcnemar_result <- mcnemar.test(table_cross)
print("Kết quả kiểm định McNemar:")
print(mcnemar_result)

# Kiểm định chính xác với correction
mcnemar_exact <- mcnemar.test(table_cross, correct = FALSE)
print("Kiểm định McNemar không hiệu chính:")
print(mcnemar_exact)

# Tính toán thêm các thống kê
n_discord <- table_cross[1,2] + table_cross[2,1]
cat("Số quan sát không concordant:", n_discord, "\n")

# Vẽ biểu đồ minh họa (tùy chọn)
par(mfrow=c(1,2))
barplot(c(sum(data$Before), sum(data$After)), 
        names.arg = c("Trước", "Sau"),
        main = "So sánh tỷ lệ hút thuốc",
        ylab = "Số lượng")

# Biểu đồ mosaic
mosaicplot(table_cross, main = "Bảng chéo McNemar")

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

import pandas as pd
import numpy as np
from scipy.stats import mcnemar
from statsmodels.stats.contingency_tables import mcnemar as mcnemar_sm
import matplotlib.pyplot as plt
import seaborn as sns

# Đọc dữ liệu từ file SPSS
try:
    data = pd.read_spss('mcnemar.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)
    data = pd.DataFrame({
        'Before': np.random.choice([0, 1], 50, p=[0.5, 0.5]),
        'After': np.random.choice([0, 1], 50, p=[0.6, 0.4])
    })

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

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

# Thực hiện kiểm định McNemar với scipy
statistic, p_value = mcnemar(contingency_table.values, exact=True)
print(f"\nKết quả kiểm định McNemar (scipy):")
print(f"Test statistic: {statistic}")
print(f"P-value: {p_value}")

# Thực hiện kiểm định McNemar với statsmodels (nhiều tùy chọn hơn)
result_sm = mcnemar_sm(contingency_table.values, exact=True)
print(f"\nKết quả kiểm định McNemar (statsmodels):")
print(f"Test statistic: {result_sm.statistic}")
print(f"P-value: {result_sm.pvalue}")

# Tính toán các chỉ số bổ sung
b = contingency_table.iloc[0, 1]  # False -> True
c = contingency_table.iloc[1, 0]  # True -> False
total_changes = b + c

print(f"\nThống kê bổ sung:")
print(f"Số người chuyển từ 0 sang 1: {b}")
print(f"Số người chuyển từ 1 sang 0: {c}")
print(f"Tổng số thay đổi: {total_changes}")

# Vẽ biểu đồ minh họa
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Biểu đồ so sánh trước và sau
before_sum = data['Before'].sum()
after_sum = data['After'].sum()
axes[0].bar(['Trước can thiệp', 'Sau can thiệp'], 
           [before_sum, after_sum], 
           color=['lightcoral', 'lightblue'])
axes[0].set_title('So sánh tỷ lệ hút thuốc')
axes[0].set_ylabel('Số lượng người hút thuốc')

# Heatmap của bảng chéo
sns.heatmap(contingency_table, annot=True, fmt='d', 
           cmap='Blues', ax=axes[1])
axes[1].set_title('Bảng chéo Before vs After')
axes[1].set_xlabel('After')
axes[1].set_ylabel('Before')

plt.tight_layout()
plt.show()

# Kết luận
alpha = 0.05
if p_value < alpha:
    print(f"\nKết luận: Với p-value = {p_value:.4f} < α = {alpha}") print("Bác bỏ giả thuyết H0. Có sự thay đổi có ý nghĩa thống kê.") else: print(f"\nKết luận: Với p-value = {p_value:.4f} >= α = {alpha}")
    print("Không bác bỏ giả thuyết H0. Không có sự thay đổi có ý nghĩa thống kê.")