Phân tích ANOVA 2 chiều

* Phân tích ANOVA hai chiều trong SPSS
* Kiểm tra giả thiết phân phối chuẩn
EXAMINE VARIABLES=lnearn BY gender section
  /PLOT BOXPLOT STEMLEAF NPPLOT
  /COMPARE GROUPS
  /STATISTICS DESCRIPTIVES
  /CINTERVAL 95
  /MISSING LISTWISE
  /NOTOTAL.

* Thực hiện phân tích ANOVA hai chiều
UNIANOVA lnearn BY gender section
  /METHOD=SSTYPE(3)
  /INTERCEPT=INCLUDE
  /PLOT=PROFILE(section*gender)
  /EMMEANS=TABLES(gender)
  /EMMEANS=TABLES(section)
  /EMMEANS=TABLES(gender*section)
  /PRINT=ETASQ DESCRIPTIVE HOMOGENEITY
  /CRITERIA=ALPHA(.05)
  /DESIGN=gender section gender*section
  /POSTHOC=section(TUKEY).

* Phân tích ANOVA hai chiều trong Stata
* Mở dữ liệu
use "vietlod.dta", clear

* Tạo biến logarit thu nhập
gen lnearn = ln(earn)

* Kiểm tra thống kê mô tả theo nhóm
bysort gender section: summarize lnearn

* Thực hiện phân tích ANOVA hai chiều
anova lnearn i.gender##i.section

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

* Tính toán estimated marginal means
margins gender
margins section
margins gender#section

* Vẽ đồ thị tương tác
marginsplot, by(gender) ///
    title("Tương tác giữa giới tính và loại hình doanh nghiệp") ///
    xtitle("Loại hình doanh nghiệp") ///
    ytitle("Log thu nhập dự đoán")

* Kiểm định Post Hoc cho biến section
pwcompare section, mcompare(tukey) effects

# Phân tích ANOVA hai chiều trong R
# Cài đặt các package cần thiết
if (!require(car)) install.packages("car")
if (!require(emmeans)) install.packages("emmeans")
if (!require(ggplot2)) install.packages("ggplot2")

library(car)
library(emmeans)
library(ggplot2)

# Đọc dữ liệu
data <- read.csv("vietlod.csv")

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

# Chuyển đổi biến thành factor
data$gender <- as.factor(data$gender)
data$section <- as.factor(data$section)

# Thống kê mô tả theo nhóm
aggregate(lnearn ~ gender + section, data = data, 
          FUN = function(x) c(mean = mean(x), sd = sd(x), n = length(x)))

# Thực hiện ANOVA hai chiều
model <- aov(lnearn ~ gender * section, data = data)
summary(model)

# Kiểm định đồng nhất phương sai (Levene test)
leveneTest(lnearn ~ gender * section, data = data)

# Estimated marginal means
emm <- emmeans(model, ~ gender * section)
summary(emm)

# Đồ thị tương tác
emmip(model, gender ~ section, CIs = TRUE) +
  labs(title = "Tương tác giữa giới tính và loại hình doanh nghiệp",
       x = "Loại hình doanh nghiệp",
       y = "Log thu nhập dự đoán")

# Kiểm định Post Hoc
pairs(emmeans(model, ~ section), adjust = "tukey")

# Kiểm tra phân phối chuẩn của residuals
qqnorm(residuals(model))
qqline(residuals(model))

# Kiểm định Shapiro-Wilk cho residuals
shapiro.test(residuals(model))

# Phân tích ANOVA hai chiều trong Python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from statsmodels.stats.anova import anova_lm
from statsmodels.formula.api import ols
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import warnings
warnings.filterwarnings('ignore')

# Đọc dữ liệu
data = pd.read_csv('vietlod.csv')

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

# Loại bỏ giá trị missing
data = data.dropna(subset=['lnearn', 'gender', 'section'])

# Thống kê mô tả theo nhóm
desc_stats = data.groupby(['gender', 'section'])['lnearn'].agg(['mean', 'std', 'count'])
print("Thống kê mô tả theo nhóm:")
print(desc_stats)

# Tạo model ANOVA hai chiều
model = ols('lnearn ~ C(gender) + C(section) + C(gender):C(section)', data=data).fit()

# Bảng ANOVA
anova_table = anova_lm(model, typ=2)
print("\nBảng ANOVA hai chiều:")
print(anova_table)

# Kiểm định đồng nhất phương sai (Levene test)
from scipy.stats import levene

groups = []
for gender in data['gender'].unique():
    for section in data['section'].unique():
        group_data = data[(data['gender'] == gender) & (data['section'] == section)]['lnearn']
        if len(group_data) > 0:
            groups.append(group_data)

levene_stat, levene_p = levene(*groups)
print(f"\nKiểm định Levene: F = {levene_stat:.4f}, p = {levene_p:.4f}")

# Đồ thị tương tác
plt.figure(figsize=(10, 6))
for gender in data['gender'].unique():
    subset = data[data['gender'] == gender]
    means = subset.groupby('section')['lnearn'].mean()
    plt.plot(means.index, means.values, marker='o', label=f'Giới tính: {gender}')

plt.xlabel('Loại hình doanh nghiệp')
plt.ylabel('Log thu nhập trung bình')
plt.title('Tương tác giữa giới tính và loại hình doanh nghiệp')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# Kiểm định Post Hoc (Tukey HSD) cho biến section
tukey_result = pairwise_tukeyhsd(data['lnearn'], data['section'], alpha=0.05)
print("\nKiểm định Post Hoc (Tukey HSD) cho loại hình doanh nghiệp:")
print(tukey_result)

# Kiểm tra phân phối chuẩn của residuals
residuals = model.resid
shapiro_stat, shapiro_p = stats.shapiro(residuals[:5000])  # Giới hạn 5000 quan sát
print(f"\nKiểm định Shapiro-Wilk cho residuals: W = {shapiro_stat:.4f}, p = {shapiro_p:.4f}")

# Q-Q plot cho residuals
plt.figure(figsize=(8, 6))
stats.probplot(residuals, dist="norm", plot=plt)
plt.title("Q-Q Plot của Residuals")
plt.show()