Phân tích ANOVA lặp một chiều

* Phân tích ANOVA đo lường lặp lại một chiều trong SPSS
* Kiểm tra thống kê mô tả
DESCRIPTIVES VARIABLES=DT2011 DT2012 DT2013
  /STATISTICS=MEAN STDDEV MIN MAX.

* Kiểm tra tính chuẩn của dữ liệu
EXAMINE VARIABLES=DT2011 DT2012 DT2013
  /PLOT BOXPLOT STEMLEAF NPPLOT
  /COMPARE GROUPS
  /STATISTICS DESCRIPTIVES
  /CINTERVAL 95
  /MISSING LISTWISE
  /NOTOTAL.

* Thực hiện phân tích ANOVA đo lường lặp lại
GLM DT2011 DT2012 DT2013
  /WSFACTOR=Year 3 Polynomial
  /MEASURE=DT
  /METHOD=SSTYPE(3)
  /PLOT=PROFILE(Year)
  /EMMEANS=TABLES(Year) COMPARE ADJ(BONFERRONI)
  /PRINT=DESCRIPTIVE ETASQ OPOWER
  /CRITERIA=ALPHA(.05)
  /WSDESIGN=Year.

* Kiểm định sphericity riêng biệt (nếu cần)
RELIABILITY
  /VARIABLES=DT2011 DT2012 DT2013
  /SCALE('ALL VARIABLES') ALL
  /MODEL=ALPHA
  /STATISTICS=DESCRIPTIVE SCALE CORR COV
  /SUMMARY=TOTAL.

* Phân tích ANOVA đo lường lặp lại một chiều trong Stata
* Mở dữ liệu
use "anova_lap_1_chieu.dta", clear

* Tạo cấu trúc dữ liệu dài (long format)
reshape long DT, i(id) j(year)

* Đặt nhãn cho biến
label define year_lab 2011 "Năm 2011" 2012 "Năm 2012" 2013 "Năm 2013"
label values year year_lab

* Thống kê mô tả
bysort year: summarize DT

* Kiểm tra phân phối chuẩn
bysort year: swilk DT

* Thực hiện ANOVA đo lường lặp lại
anova DT year id, repeated(year)

* Kiểm định sphericity (sử dụng epsilon)
estat sphericity

* So sánh từng cặp với điều chỉnh Bonferroni
pwcompare year, mcompare(bonferroni) effects

* Vẽ đồ thị profile plot
egen mean_DT = mean(DT), by(year)
twoway line mean_DT year, ///
    title("Doanh thu trung bình theo năm") ///
    xtitle("Năm") ///
    ytitle("Doanh thu trung bình (triệu đồng)") ///
    xlabel(2011 2012 2013) ///
    lwidth(medium) ///
    msymbol(circle) ///
    msize(medium)

* Kiểm tra hiệu ứng kích thước
estat esize

* Tính toán confidence intervals
margins year
marginsplot, ///
    title("Doanh thu trung bình với khoảng tin cậy 95%") ///
    xtitle("Năm") ///
    ytitle("Doanh thu dự đoán (triệu đồng)")

# Phân tích ANOVA đo lường lặp lại một chiều trong R
# Cài đặt và load các package cần thiết
if (!require(car)) install.packages("car")
if (!require(ez)) install.packages("ez")
if (!require(reshape2)) install.packages("reshape2")
if (!require(ggplot2)) install.packages("ggplot2")
if (!require(emmeans)) install.packages("emmeans")

library(car)
library(ez)
library(reshape2)
library(ggplot2)
library(emmeans)

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

# Tạo cấu trúc dữ liệu dài
data_long <- melt(data, 
                  id.vars = "id",
                  measure.vars = c("DT2011", "DT2012", "DT2013"),
                  variable.name = "year",
                  value.name = "DT")

# Đặt lại nhãn cho biến year
data_long$year <- factor(data_long$year, 
                        levels = c("DT2011", "DT2012", "DT2013"),
                        labels = c("2011", "2012", "2013"))

# Thống kê mô tả
aggregate(DT ~ year, data = data_long, 
          FUN = function(x) c(mean = mean(x), 
                             sd = sd(x), 
                             n = length(x)))

# Kiểm tra phân phối chuẩn
by(data_long$DT, data_long$year, shapiro.test)

# Thực hiện ANOVA đo lường lặp lại
anova_result <- ezANOVA(data = data_long,
                        dv = DT,
                        wid = id,
                        within = year,
                        detailed = TRUE,
                        type = 3)

print(anova_result)

# Kiểm định sphericity (Mauchly's test)
print(anova_result$`Mauchly's Test for Sphericity`)

# Điều chỉnh sphericity nếu cần
if (anova_result$`Mauchly's Test for Sphericity`$p < 0.05) {
  print("Vi phạm giả thiết sphericity - sử dụng điều chỉnh")
  print(anova_result$`Sphericity Corrections`)
}

# So sánh từng cặp
model <- aov(DT ~ year + Error(id/year), data = data_long)
emm <- emmeans(model, ~ year)
pairs(emm, adjust = "bonferroni")

# Vẽ đồ thị profile plot
ggplot(data_long, aes(x = year, y = DT, group = 1)) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  stat_summary(fun = mean, geom = "line", aes(group = 1), size = 1) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1) +
  labs(title = "Doanh thu trung bình theo năm",
       x = "Năm",
       y = "Doanh thu trung bình (triệu đồng)") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

# Boxplot để kiểm tra phân phối
ggplot(data_long, aes(x = year, y = DT)) +
  geom_boxplot(fill = "lightblue", alpha = 0.7) +
  geom_jitter(width = 0.2, alpha = 0.5) +
  labs(title = "Phân phối doanh thu theo năm",
       x = "Năm",
       y = "Doanh thu (triệu đồng)") +
  theme_minimal()

# Tính effect size
eta_squared <- anova_result$ANOVA$SSn[1] / 
               (anova_result$ANOVA$SSn[1] + anova_result$ANOVA$SSd[1])
print(paste("Eta squared:", round(eta_squared, 3)))

# Phân tích ANOVA đo lường lặp lại một 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 scipy.stats import f_oneway
from statsmodels.stats.anova import AnovaRM
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import warnings
warnings.filterwarnings('ignore')

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

# Tạo cấu trúc dữ liệu dài
data_long = pd.melt(data, 
                    id_vars=['id'],
                    value_vars=['DT2011', 'DT2012', 'DT2013'],
                    var_name='year',
                    value_name='DT')

# Đặt lại nhãn cho biến year
data_long['year'] = data_long['year'].map({
    'DT2011': '2011',
    'DT2012': '2012', 
    'DT2013': '2013'
})

# Thống kê mô tả
desc_stats = data_long.groupby('year')['DT'].agg(['mean', 'std', 'count'])
print("Thống kê mô tả:")
print(desc_stats)

# Kiểm tra phân phối chuẩn
print("\nKiểm định Shapiro-Wilk cho từng năm:")
for year in data_long['year'].unique():
    subset = data_long[data_long['year'] == year]['DT']
    stat, p_value = stats.shapiro(subset)
    print(f"Năm {year}: W = {stat:.4f}, p = {p_value:.4f}")

# Thực hiện ANOVA đo lường lặp lại
anova_rm = AnovaRM(data_long, 'DT', 'id', within=['year'])
anova_result = anova_rm.fit()

print("\nKết quả ANOVA đo lường lặp lại:")
print(anova_result.summary())

# Kiểm định sphericity (Mauchly's test) - approximation
def mauchly_test(data_wide):
    """
    Xấp xỉ kiểm định Mauchly cho sphericity
    """
    # Tính ma trận hiệp phương sai
    cov_matrix = np.cov(data_wide.T)
    
    # Tính determinant và trace
    det_cov = np.linalg.det(cov_matrix)
    trace_cov = np.trace(cov_matrix)
    
    # Tính epsilon (Greenhouse-Geisser correction)
    k = len(data_wide.columns)
    epsilon = (trace_cov / k) ** 2 / np.mean(np.diag(cov_matrix))
    
    return epsilon

# Áp dụng kiểm định sphericity
data_wide = data.set_index('id')[['DT2011', 'DT2012', 'DT2013']]
epsilon = mauchly_test(data_wide)
print(f"\nGreenhouse-Geisser epsilon: {epsilon:.4f}")

# So sánh từng cặp (Post-hoc test)
print("\nSo sánh từng cặp (Tukey HSD):")
tukey_result = pairwise_tukeyhsd(data_long['DT'], data_long['year'], alpha=0.05)
print(tukey_result)

# Vẽ đồ thị profile plot
plt.figure(figsize=(10, 6))

# Tính trung bình và standard error
means = data_long.groupby('year')['DT'].mean()
sems = data_long.groupby('year')['DT'].sem()

plt.errorbar(means.index, means.values, yerr=sems.values, 
             marker='o', linewidth=2, markersize=8, capsize=5)
plt.xlabel('Năm')
plt.ylabel('Doanh thu trung bình (triệu đồng)')
plt.title('Doanh thu trung bình theo năm với khoảng tin cậy 95%')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Boxplot
plt.figure(figsize=(10, 6))
sns.boxplot(data=data_long, x='year', y='DT')
plt.xlabel('Năm')
plt.ylabel('Doanh thu (triệu đồng)')
plt.title('Phân phối doanh thu theo năm')
plt.tight_layout()
plt.show()

# Tính effect size (eta squared)
def eta_squared(anova_result):
    """Tính eta squared từ kết quả ANOVA"""
    ss_effect = anova_result.anova_table['sum_sq'].iloc[0]
    ss_total = anova_result.anova_table['sum_sq'].sum()
    return ss_effect / ss_total

eta_sq = eta_squared(anova_result)
print(f"\nEta squared: {eta_sq:.3f}")

# Phân loại effect size
if eta_sq < 0.01:
    effect_size = "nhỏ"
elif eta_sq < 0.06:
    effect_size = "trung bình"
elif eta_sq < 0.14:
    effect_size = "lớn"
else:
    effect_size = "rất lớn"

print(f"Kích thước hiệu ứng: {effect_size}")

# Ma trận tương quan giữa các thời điểm
print("\nMa trận tương quan giữa các năm:")
correlation_matrix = data_wide.corr()
print(correlation_matrix.round(3))

# Heatmap cho ma trận tương quan
plt.figure(figsize=(8, 6))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', center=0)
plt.title('Ma trận tương quan giữa các năm')
plt.tight_layout()
plt.show()