Phân tích phương sai một chiều

* Phân tích phương sai một chiều 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'.

* Tạo biến logarit thu nhập nếu chưa có
COMPUTE lnearn = LN(earn).
VARIABLE LABELS lnearn 'Logarit tự nhiên của thu nhập'.
EXECUTE.

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

* Kiểm tra phân phối chuẩn cho từng nhóm
EXAMINE VARIABLES=lnearn BY structure
  /PLOT=HISTOGRAM NORMPLOT
  /STATISTICS=DESCRIPTIVES
  /CINTERVAL 95
  /MISSING LISTWISE
  /NOTOTAL.

* Phân tích phương sai một chiều
ONEWAY lnearn BY structure
  /STATISTICS DESCRIPTIVES HOMOGENEITY
  /PLOT MEANS
  /POSTHOC=TUKEY GAMES ALPHA(0.05)
  /MISSING ANALYSIS.

* Kiểm định robust (Welch) khi phương sai không đồng nhất
ONEWAY lnearn BY structure
  /STATISTICS DESCRIPTIVES HOMOGENEITY WELCH
  /POSTHOC=GAMES ALPHA(0.05)
  /MISSING ANALYSIS.

* Kiểm tra điểm dị biệt
EXAMINE VARIABLES=lnearn BY structure
  /PLOT=BOXPLOT
  /STATISTICS=EXTREME
  /CINTERVAL 95
  /MISSING LISTWISE
  /NOTOTAL.

* Tính effect size (eta-squared)
COMPUTE eta_squared = SS_between / SS_total.
* Chú thích: Cần tính thủ công từ bảng ANOVA

* Xuất kết quả
DISPLAY "Phân tích phương sai một chiều hoàn tất".
DISPLAY "Kiểm tra giả thiết và sử dụng post-hoc phù hợp".

* Phân tích phương sai một chiều 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 "vietlod.dta", clear

* Tạo biến logarit thu nhập
gen lnearn = ln(earn)
label variable lnearn "Logarit tự nhiên của thu nhập"

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

* Thống kê mô tả theo nhóm
by structure, sort: summarize lnearn, detail
tabstat lnearn, by(structure) statistics(n mean sd min max)

* Kiểm tra phân phối chuẩn
by structure: swilk lnearn
histogram lnearn, by(structure) normal

* Kiểm tra điểm dị biệt
by structure: summarize lnearn, detail
graph box lnearn, over(structure)

* Phân tích phương sai một chiều
oneway lnearn structure, tabulate

* Kiểm định đồng nhất phương sai
robvar lnearn, by(structure)

* Kiểm định Welch (robust)
oneway lnearn structure, welch

* Post-hoc tests
pwmean lnearn, over(structure) mcompare(tukey) effects
pwmean lnearn, over(structure) mcompare(bonferroni) effects

* Biểu đồ means
graph hbar (mean) lnearn, over(structure) ///
    title("Thu nhập trung bình theo lĩnh vực") ///
    ytitle("Logarit thu nhập") ///
    blabel(bar, format(%9.2f))

* Tính effect size
anova lnearn structure
estat esize

* Kết luận
display "Phân tích phương sai một chiều hoàn tất"
display "Sử dụng Welch test nếu phương sai không đồng nhất"

# Phân tích phương sai một chiều 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(car)) install.packages("car")
if (!require(multcomp)) install.packages("multcomp")
if (!require(onewaytests)) install.packages("onewaytests")

library(foreign)
library(dplyr)
library(ggplot2)
library(car)
library(multcomp)
library(onewaytests)

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

# Tạo biến logarit thu nhập
data$lnearn <- log(data$earn)

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

# Chuyển structure thành factor
data$structure <- as.factor(data$structure)
levels(data$structure) <- c("Dịch vụ", "Công nghiệp", "Nông nghiệp")

# Thống kê mô tả theo nhóm
summary_stats <- data %>%
  group_by(structure) %>%
  summarise(
    n = n(),
    mean = mean(lnearn, na.rm = TRUE),
    sd = sd(lnearn, na.rm = TRUE),
    min = min(lnearn, na.rm = TRUE),
    max = max(lnearn, na.rm = TRUE),
    se = sd/sqrt(n)
  )

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

# Kiểm tra phân phối chuẩn
shapiro_results <- data %>%
  group_by(structure) %>%
  summarise(
    shapiro_statistic = shapiro.test(lnearn)$statistic,
    shapiro_p_value = shapiro.test(lnearn)$p.value
  )

print("Kiểm định Shapiro-Wilk:")
print(shapiro_results)

# Vẽ biểu đồ kiểm tra phân phối
ggplot(data, aes(x = lnearn)) +
  geom_histogram(bins = 30, alpha = 0.7, fill = "skyblue") +
  facet_wrap(~structure) +
  geom_density(alpha = 0.5, color = "red") +
  labs(title = "Phân phối logarit thu nhập theo lĩnh vực",
       x = "Logarit thu nhập",
       y = "Tần số") +
  theme_minimal()

# Boxplot để kiểm tra outliers
ggplot(data, aes(x = structure, y = lnearn, fill = structure)) +
  geom_boxplot(alpha = 0.7) +
  labs(title = "Boxplot logarit thu nhập theo lĩnh vực",
       x = "Lĩnh vực",
       y = "Logarit thu nhập") +
  theme_minimal() +
  theme(legend.position = "none")

# Phân tích phương sai một chiều
aov_model <- aov(lnearn ~ structure, data = data)
aov_summary <- summary(aov_model)
print("Kết quả ANOVA:")
print(aov_summary)

# Kiểm định đồng nhất phương sai (Levene)
levene_test <- leveneTest(lnearn ~ structure, data = data)
print("Kiểm định Levene:")
print(levene_test)

# Kiểm định Welch (robust)
welch_test <- oneway.test(lnearn ~ structure, data = data)
print("Kiểm định Welch:")
print(welch_test)

# Post-hoc tests
# Tukey HSD (khi phương sai đồng nhất)
tukey_results <- TukeyHSD(aov_model)
print("Kiểm định Tukey HSD:")
print(tukey_results)

# Games-Howell (khi phương sai không đồng nhất)
games_howell <- paircomp(lnearn ~ structure, data = data)
print("Kiểm định Games-Howell:")
print(games_howell)

# Tính effect size (eta-squared)
aov_summary_table <- summary(aov_model)[[1]]
eta_squared <- aov_summary_table[1, "Sum Sq"] / sum(aov_summary_table[, "Sum Sq"])
cat("Effect size (eta-squared):", eta_squared, "\n")

# Biểu đồ means với confidence intervals
ggplot(summary_stats, aes(x = structure, y = mean, fill = structure)) +
  geom_col(alpha = 0.7) +
  geom_errorbar(aes(ymin = mean - 1.96*se, ymax = mean + 1.96*se), 
                width = 0.2) +
  labs(title = "Thu nhập trung bình theo lĩnh vực (95% CI)",
       x = "Lĩnh vực",
       y = "Logarit thu nhập trung bình") +
  theme_minimal() +
  theme(legend.position = "none")

# Kiểm tra residuals
residuals <- residuals(aov_model)
fitted_values <- fitted(aov_model)

# Q-Q plot cho residuals
qqnorm(residuals)
qqline(residuals)

# Residuals vs Fitted
plot(fitted_values, residuals,
     xlab = "Fitted values", ylab = "Residuals",
     main = "Residuals vs Fitted Values")
abline(h = 0, col = "red")

# Kết luận
alpha <- 0.05
f_statistic <- aov_summary[[1]][1, "F value"]
p_value <- aov_summary[[1]][1, "Pr(>F)"]

if (p_value < alpha) {
  cat("\nKết luận: Với F =", f_statistic, "và p-value =", p_value, 
      "< α =", alpha, "\n") cat("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.\n") } else { cat("\nKết luận: Với F =", f_statistic, "và p-value =", 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ê về thu nhập giữa các lĩnh vực.\n")
}

cat("Phân tích phương sai một chiều hoàn tất!\n")

# Phân tích phương sai một chiều 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
from statsmodels.stats.anova import anova_lm
from statsmodels.formula.api import ols
import warnings
warnings.filterwarnings('ignore')

# Đọ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.35, 0.35, 0.3])
    
    # 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.lognormal(3.5, 0.8)
        elif structure[i] == 2:  # Công nghiệp
            earn[i] = np.random.lognormal(3.2, 0.7)
        else:  # Nông nghiệp
            earn[i] = np.random.lognormal(2.8, 0.6)
    
    data = pd.DataFrame({
        'earn': earn,
        'structure': structure
    })

# Tạo biến logarit thu nhập
data['lnearn'] = np.log(data['earn'])

# Chuyển đổi structure thành categorical
data['structure'] = data['structure'].astype('category')
data['structure'] = data['structure'].cat.rename_categories({1: 'Dịch vụ', 2: 'Công nghiệp', 3: 'Nông nghiệp'})

# 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')['lnearn'].agg([
    'count', 'mean', 'std', 'min', 'max'
]).round(4)
print(summary_stats)

# Kiểm tra phân phối chuẩn cho từng nhóm
print("\nKiểm định Shapiro-Wilk cho từng nhóm:")
for group in data['structure'].cat.categories:
    group_data = data[data['structure'] == group]['lnearn']
    if len(group_data) > 3:
        shapiro_stat, shapiro_p = stats.shapiro(group_data.sample(min(5000, len(group_data))))
        print(f"{group}: Shapiro-Wilk statistic = {shapiro_stat:.4f}, p-value = {shapiro_p:.4f}")

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

# Histogram
for i, group in enumerate(data['structure'].cat.categories):
    group_data = data[data['structure'] == group]['lnearn']
    axes[0, 0].hist(group_data, alpha=0.7, label=group, bins=30)
axes[0, 0].set_title('Histogram logarit thu nhập theo lĩnh vực')
axes[0, 0].set_xlabel('Logarit thu nhập')
axes[0, 0].set_ylabel('Tần số')
axes[0, 0].legend()

# Boxplot
data.boxplot('lnearn', by='structure', ax=axes[0, 1])
axes[0, 1].set_title('Boxplot logarit thu nhập theo lĩnh vực')
axes[0, 1].set_xlabel('Lĩnh vực')
axes[0, 1].set_ylabel('Logarit thu nhập')

# Q-Q plot
stats.probplot(data['lnearn'], dist="norm", plot=axes[1, 0])
axes[1, 0].set_title('Q-Q plot cho logarit thu nhập')

# Violin plot
data_pivot = [data[data['structure'] == group]['lnearn'] for group in data['structure'].cat.categories]
axes[1, 1].violinplot(data_pivot, positions=[1, 2, 3])
axes[1, 1].set_title('Violin plot logarit thu nhập theo lĩnh vực')
axes[1, 1].set_ylabel('Logarit thu nhập')
axes[1, 1].set_xticks([1, 2, 3])
axes[1, 1].set_xticklabels(data['structure'].cat.categories)

plt.tight_layout()
plt.show()

# Phân tích phương sai một chiều
print("\nPhân tích phương sai một chiều:")

# Cách 1: Sử dụng scipy.stats
groups = [data[data['structure'] == group]['lnearn'] for group in data['structure'].cat.categories]
f_stat, p_value = stats.f_oneway(*groups)
print(f"F-statistic: {f_stat:.6f}")
print(f"P-value: {p_value:.6f}")

# Cách 2: Sử dụng statsmodels (chi tiết hơn)
model = ols('lnearn ~ C(structure)', data=data).fit()
anova_table = anova_lm(model, typ=2)
print("\nBảng ANOVA chi tiết:")
print(anova_table)

# Kiểm định đồng nhất phương sai (Levene)
levene_stat, levene_p = stats.levene(*groups)
print(f"\nKiểm định Levene:")
print(f"Levene statistic: {levene_stat:.6f}")
print(f"P-value: {levene_p:.6f}")

# Kiểm định Welch (robust)
welch_stat, welch_p = stats.alexandergovern(*groups)
print(f"\nKiểm định Alexander

# Kiểm định Welch (robust) - thay thế khi phương sai không đồng nhất
try:
    # Sử dụng Alexander-Govern test (tương tự Welch)
    welch_stat, welch_p = stats.alexandergovern(*groups)
    print(f"\nKiểm định Alexander-Govern (tương tự Welch):")
    print(f"Statistic: {welch_stat:.6f}")
    print(f"P-value: {welch_p:.6f}")
except:
    # Nếu không có, sử dụng welch_anova từ pingouin
    print(f"\nKiểm định Welch không khả dụng với scipy hiện tại")

# Post-hoc tests
print("\nKiểm định post-hoc:")

# Tukey HSD
tukey_results = sp.posthoc_tukey(data, val_col='lnearn', group_col='structure')
print("Tukey HSD:")
print(tukey_results)

# Games-Howell (robust khi phương sai không đồng nhất)
games_howell_results = sp.posthoc_gameshowell(data, val_col='lnearn', group_col='structure')
print("\nGames-Howell:")
print(games_howell_results)

# Tính effect size (eta-squared)
ss_between = anova_table['sum_sq'][0]
ss_total = anova_table['sum_sq'].sum()
eta_squared = ss_between / ss_total
print(f"\nEffect size (eta-squared): {eta_squared:.6f}")

# Phân loại effect size
if eta_squared >= 0.14:
    effect_size_interpretation = "Large effect"
elif eta_squared >= 0.06:
    effect_size_interpretation = "Medium effect"
else:
    effect_size_interpretation = "Small effect"

print(f"Diễn giải effect size: {effect_size_interpretation}")

# Biểu đồ means với confidence intervals
plt.figure(figsize=(10, 6))
means = summary_stats['mean']
errors = summary_stats['std'] / np.sqrt(summary_stats['count']) * 1.96

plt.bar(range(len(means)), means, yerr=errors, 
        capsize=5, alpha=0.7, color=['skyblue', 'lightgreen', 'lightcoral'])
plt.xlabel('Lĩnh vực')
plt.ylabel('Logarit thu nhập trung bình')
plt.title('Thu nhập trung bình theo lĩnh vực (95% CI)')
plt.xticks(range(len(means)), data['structure'].cat.categories)
plt.grid(axis='y', alpha=0.3)

# Thêm giá trị lên đầu cột
for i, (mean, error) in enumerate(zip(means, errors)):
    plt.text(i, mean + error + 0.05, f'{mean:.3f}', 
             ha='center', va='bottom', fontweight='bold')

plt.tight_layout()
plt.show()

# Kiểm tra residuals
residuals = model.resid
fitted_values = model.fittedvalues

# Residuals vs Fitted
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.scatter(fitted_values, residuals, alpha=0.6)
plt.axhline(y=0, color='red', linestyle='--')
plt.xlabel('Fitted Values')
plt.ylabel('Residuals')
plt.title('Residuals vs Fitted Values')

# Q-Q plot của residuals
plt.subplot(1, 2, 2)
stats.probplot(residuals, dist="norm", plot=plt)
plt.title('Q-Q Plot của Residuals')

plt.tight_layout()
plt.show()

# Tạo bảng tóm tắt kết quả
results_summary = pd.DataFrame({
    'Lĩnh vực': data['structure'].cat.categories,
    'N': summary_stats['count'],
    'Mean': summary_stats['mean'],
    'Std': summary_stats['std'],
    'Min': summary_stats['min'],
    'Max': summary_stats['max']
})

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

# Kết luận
alpha = 0.05
if p_value < alpha:
    print(f"\nKết luận: Với F = {f_stat:.6f} và 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 giữa các lĩnh vực.")
    
    if levene_p < alpha:
        print("Lưu ý: Phương sai không đồng nhất (Levene p < 0.05), nên sử dụng Games-Howell post-hoc.") else: print("Phương sai đồng nhất (Levene p >= 0.05), có thể sử dụng Tukey HSD post-hoc.")
        
else:
    print(f"\nKết luận: Với F = {f_stat:.6f} và p-value = {p_value:.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.")

print(f"\nEffect size: {eta_squared:.6f} ({effect_size_interpretation})")
print("Phân tích phương sai một chiều hoàn tất!")