Phân tích hiệp phương sai – ANCOVA

* Phân tích hiệp phương sai (ANCOVA) với SPSS
* Tác giả: Hướng dẫn thực hành kinh tế lượng

* Thống kê mô tả cơ bản
DESCRIPTIVES VARIABLES=After Before BY Group
  /STATISTICS=MEAN STDDEV MIN MAX.

* Phân tích hiệp phương sai chính
UNIANOVA After BY Group WITH Before
  /METHOD=SSTYPE(3)
  /INTERCEPT=INCLUDE
  /PRINT=DESCRIPTIVE PARAMETER HOMOGENEITY
  /CRITERIA=ALPHA(.05)
  /DESIGN=Group Before.

* So sánh từng cặp với điều chỉnh Bonferroni
UNIANOVA After BY Group WITH Before
  /EMMEANS=TABLES(Group) WITH(Before=MEAN) COMPARE ADJ(BONFERRONI)
  /PRINT=DESCRIPTIVE PARAMETER.

* Kiểm tra giả thuyết đồng nhất hệ số (tương tác)
UNIANOVA After BY Group WITH Before
  /METHOD=SSTYPE(3)
  /INTERCEPT=INCLUDE
  /PRINT=PARAMETER
  /CRITERIA=ALPHA(.05)
  /DESIGN=Group Before Group*Before.

* Kiểm tra tính tuyến tính bằng đồ thị scatter
GRAPH
  /SCATTERPLOT(BIVAR)=Before WITH After BY Group
  /MISSING=LISTWISE.

* Kiểm tra phân phối chuẩn của residuals
UNIANOVA After BY Group WITH Before
  /SAVE=RESID
  /PRINT=DESCRIPTIVE.

EXAMINE VARIABLES=RES_1
  /PLOT NPPLOT
  /STATISTICS DESCRIPTIVES
  /CINTERVAL 95
  /MISSING LISTWISE
  /NOTOTAL.

* Phân tích hiệp phương sai (ANCOVA) trong Stata
* Sử dụng General Linear Model

* Tải và mô tả dữ liệu
use "ancova_data.dta", clear
describe After Before Group
summarize After Before, detail

* Thống kê mô tả theo nhóm
bysort Group: summarize After Before

* Kiểm tra tương quan giữa các biến
corr After Before
bysort Group: corr After Before

* Thực hiện ANCOVA
anova After Group Before

* Mô hình tương đương với regress
regress After i.Group Before

* So sánh các nhóm (post-hoc)
margins i.Group, at(Before=`r(mean)')
marginsplot

* Kiểm tra giả thuyết đồng nhất hệ số (interaction)
regress After i.Group##c.Before
test 2.Group#c.Before 3.Group#c.Before

* Kiểm tra đồng nhất phương sai (Levene test)
robvar After, by(Group)

* Vẽ đồ thị scatter để kiểm tra tính tuyến tính
scatter After Before, by(Group) || lfit After Before, by(Group)

* Kiểm tra phân phối chuẩn của residuals
regress After i.Group Before
predict residuals, residuals
histogram residuals, normal
qnorm residuals

* Kiểm tra điểm dị biệt
regress After i.Group Before
predict leverage, leverage
predict cook, cooksd
list if cook > 4/e(N)

* Estimated marginal means
margins Group, at(Before=`r(mean)')

# Phân tích hiệp phương sai (ANCOVA) với R
# Sử dụng các gói car, emmeans và ggplot2

# Tải các thư viện cần thiết
library(car)
library(emmeans)
library(ggplot2)
library(dplyr)
library(broom)

# Đọc dữ liệu
data <- read.csv("ancova_data.csv")
data$Group <- as.factor(data$Group) str(data) # Thống kê mô tả summary(data) data %>% 
  group_by(Group) %>%
  summarise(
    n = n(),
    mean_after = mean(After),
    sd_after = sd(After),
    mean_before = mean(Before),
    sd_before = sd(Before)
  )

# Vẽ đồ thị khám phá dữ liệu
ggplot(data, aes(x = Before, y = After, color = Group)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(title = "Mối quan hệ giữa Before và After theo Group")

# Boxplot để kiểm tra outliers
ggplot(data, aes(x = Group, y = After)) +
  geom_boxplot() +
  labs(title = "Phân phối After theo Group")

# Thực hiện ANCOVA
ancova_model <- aov(After ~ Group + Before, data = data)
summary(ancova_model)

# Mô hình tương đương với lm
lm_model <- lm(After ~ Group + Before, data = data)
summary(lm_model)

# Kiểm tra các giả thuyết
# 1. Đồng nhất phương sai (Levene test)
leveneTest(After ~ Group, data = data)

# 2. Phân phối chuẩn của residuals
shapiro.test(residuals(ancova_model))

# 3. Tính tuyến tính - vẽ đồ thị residuals vs fitted
plot(ancova_model, which = 1)

# 4. Q-Q plot cho phân phối chuẩn
plot(ancova_model, which = 2)

# 5. Kiểm tra điểm dị biệt
plot(ancova_model, which = 4)

# Kiểm tra giả thuyết đồng nhất hệ số (interaction)
interaction_model <- lm(After ~ Group * Before, data = data)
anova(lm_model, interaction_model)

# Estimated Marginal Means
emm <- emmeans(ancova_model, ~ Group)
print(emm)

# So sánh từng cặp với điều chỉnh Bonferroni
pairwise <- pairs(emm, adjust = "bonferroni")
print(pairwise)

# Vẽ đồ thị estimated marginal means
plot(emm) + 
  labs(title = "Estimated Marginal Means với 95% CI")

# Tính effect size (eta squared)
effectsize::eta_squared(ancova_model)

# Báo cáo kết quả chi tiết
anova_table <- anova(ancova_model)
print(anova_table)

# Confidence intervals cho các hệ số
confint(lm_model)

# Phân tích hiệp phương sai (ANCOVA) với Python
# Sử dụng scipy, statsmodels và seaborn

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
from statsmodels.stats.multicomp import pairwise_tukeyhsd
from statsmodels.stats.diagnostic import het_breuschpagan
import statsmodels.api as sm

# Đọc dữ liệu
data = pd.read_csv('ancova_data.csv')
data['Group'] = data['Group'].astype('category')
print(data.info())
print(data.describe())

# Thống kê mô tả theo nhóm
desc_stats = data.groupby('Group').agg({
    'After': ['count', 'mean', 'std'],
    'Before': ['mean', 'std']
}).round(3)
print("Thống kê mô tả theo nhóm:")
print(desc_stats)

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Scatter plot cho mối quan hệ Before vs After
sns.scatterplot(data=data, x='Before', y='After', hue='Group', ax=axes[0,0])
sns.regplot(data=data, x='Before', y='After', scatter=False, ax=axes[0,0])
axes[0,0].set_title('Mối quan hệ Before vs After')

# Box plot cho After theo Group
sns.boxplot(data=data, x='Group', y='After', ax=axes[0,1])
axes[0,1].set_title('Phân phối After theo Group')

# Histogram cho Before
sns.histplot(data=data, x='Before', hue='Group', ax=axes[1,0])
axes[1,0].set_title('Phân phối Before theo Group')

# Correlation matrix
corr_matrix = data[['After', 'Before']].corr()
sns.heatmap(corr_matrix, annot=True, ax=axes[1,1])
axes[1,1].set_title('Ma trận tương quan')

plt.tight_layout()
plt.show()

# Thực hiện ANCOVA
model = ols('After ~ C(Group) + Before', data=data).fit()
anova_table = anova_lm(model, typ=2)
print("\nKết quả ANCOVA:")
print(anova_table)

# Tóm tắt mô hình
print("\nTóm tắt mô hình:")
print(model.summary())

# Kiểm tra các giả thuyết
print("\n=== KIỂM TRA GIẢ THUYẾT ===")

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

groups = [group['After'].values for name, group in data.groupby('Group')]
levene_stat, levene_p = levene(*groups)
print(f"Levene test: statistic = {levene_stat:.4f}, p-value = {levene_p:.4f}")

# 2. Kiểm tra phân phối chuẩn của residuals
residuals = model.resid
shapiro_stat, shapiro_p = stats.shapiro(residuals)
print(f"Shapiro-Wilk test: statistic = {shapiro_stat:.4f}, p-value = {shapiro_p:.4f}")

# 3. Vẽ đồ thị diagnostic
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Residuals vs Fitted
fitted = model.fittedvalues
axes[0,0].scatter(fitted, residuals)
axes[0,0].axhline(y=0, color='red', linestyle='--')
axes[0,0].set_xlabel('Fitted values')
axes[0,0].set_ylabel('Residuals')
axes[0,0].set_title('Residuals vs Fitted')

# Q-Q plot
stats.probplot(residuals, dist="norm", plot=axes[0,1])
axes[0,1].set_title('Q-Q Plot')

# Histogram của residuals
axes[1,0].hist(residuals, bins=20, edgecolor='black')
axes[1,0].set_xlabel('Residuals')
axes[1,0].set_ylabel('Frequency')
axes[1,0].set_title('Histogram of Residuals')

# Scale-Location plot
axes[1,1].scatter(fitted, np.sqrt(np.abs(residuals)))
axes[1,1].set_xlabel('Fitted values')
axes[1,1].set_ylabel('√|Residuals|')
axes[1,1].set_title('Scale-Location Plot')

plt.tight_layout()
plt.show()

# 4. Kiểm tra tương tác (đồng nhất hệ số)
interaction_model = ols('After ~ C(Group) * Before', data=data).fit()
interaction_anova = anova_lm(interaction_model, typ=2)
print("\nKiểm tra tương tác:")
print(interaction_anova)

# So sánh mô hình có và không có tương tác
print("\nSo sánh mô hình:")
print(anova_lm(model, interaction_model))

# Estimated Marginal Means (thủ công)
mean_before = data['Before'].mean()
group_means = []

for group in data['Group'].cat.categories:
    # Tính predicted value cho mỗi group tại mean của Before
    temp_data = pd.DataFrame({
        'Group': [group],
        'Before': [mean_before]
    })
    temp_data['Group'] = temp_data['Group'].astype('category')
    predicted = model.predict(temp_data)
    group_means.append(predicted[0])

emm_df = pd.DataFrame({
    'Group': data['Group'].cat.categories,
    'EMM': group_means
})
print("\nEstimated Marginal Means:")
print(emm_df)

# Post-hoc comparisons (Bonferroni)
# Tính confidence intervals cho EMM
conf_int = model.conf_int()
print("\nConfidence Intervals:")
print(conf_int)

# Effect size (eta squared)
ss_effect = anova_table.loc['C(Group)', 'sum_sq']
ss_total = anova_table['sum_sq'].sum()
eta_squared = ss_effect / ss_total
print(f"\nEta squared cho Group: {eta_squared:.4f}")

# Tóm tắt kết quả
print("\n=== TÓM TẮT KẾT QUẢ ===")
print(f"F-statistic cho Group: {anova_table.loc['C(Group)', 'F']:.4f}")
print(f"p-value cho Group: {anova_table.loc['C(Group)', 'PR(>F)']:.4f}")
print(f"R-squared: {model.rsquared:.4f}")
print(f"Adjusted R-squared: {model.rsquared_adj:.4f}")