Hồi quy tuyến tính đa biến

* Hồi quy tuyến tính đa biến 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=read age mem iq
  /STATISTICS=MEAN STDDEV MIN MAX.

* Ma trận tương quan
CORRELATIONS
  /VARIABLES=read age mem iq
  /PRINT=TWOTAIL NOSIG
  /MISSING=PAIRWISE.

* Hồi quy tuyến tính đa biến
REGRESSION
  /MISSING LISTWISE
  /STATISTICS COEFF OUTS R ANOVA COLLIN
  /CRITERIA=PIN(.05) POUT(.10)
  /NOORIGIN 
  /DEPENDENT read
  /METHOD=ENTER age mem iq
  /SCATTERPLOT=(*ZRESID ,*ZPRED)
  /RESIDUALS DURBIN HISTOGRAM(ZRESID) NORMPROB(ZRESID).

* Kiểm tra điểm dị biệt
REGRESSION
  /MISSING LISTWISE
  /STATISTICS COEFF OUTS R ANOVA
  /CRITERIA=PIN(.05) POUT(.10)
  /NOORIGIN 
  /DEPENDENT read
  /METHOD=ENTER age mem iq
  /SAVE COOK LEVER DFBETA DFFIT
  /RESIDUALS OUTLIERS(COOK) OUTLIERS(LEVER).

* Kiểm tra tính tuyến tính
GRAPH
  /SCATTERPLOT(BIVAR)=age WITH read
  /MISSING=LISTWISE.

GRAPH
  /SCATTERPLOT(BIVAR)=mem WITH read
  /MISSING=LISTWISE.

GRAPH
  /SCATTERPLOT(BIVAR)=iq WITH read
  /MISSING=LISTWISE.

* Biến đổi dữ liệu nếu cần
COMPUTE ln_read = LN(read).
COMPUTE sqrt_read = SQRT(read).
EXECUTE.

* Hồi quy tuyến tính đa biến trong Stata
* Sử dụng dữ liệu về khả năng đọc hiểu

* Tải và mô tả dữ liệu
use "ols_data.dta", clear
describe read age mem iq
summarize read age mem iq, detail

* Ma trận tương quan
corr read age mem iq
pwcorr read age mem iq, sig

* Hồi quy tuyến tính đa biến
regress read age mem iq

* Kiểm tra các giả thuyết
* 1. Kiểm tra đa cộng tuyến (VIF)
estat vif

* 2. Kiểm tra phương sai đồng nhất (Breusch-Pagan test)
estat hettest

* 3. Kiểm tra tự tương quan (Durbin-Watson)
estat dwatson

* 4. Kiểm tra phân phối chuẩn của phần dư
predict residuals, residuals
histogram residuals, normal
qnorm residuals

* 5. Kiểm tra tính tuyến tính
scatter read age || lfit read age
scatter read mem || lfit read mem
scatter read iq || lfit read iq

* 6. Kiểm tra điểm dị biệt
predict leverage, leverage
predict cooksd, cooksd
predict dfbeta_age, dfbeta(age)
predict dfbeta_mem, dfbeta(mem)
predict dfbeta_iq, dfbeta(iq)

* Liệt kê các quan sát có Cook's distance > 4/n
scalar threshold = 4/_N
list if cooksd > threshold

* Robust standard errors nếu có phương sai không đồng nhất
regress read age mem iq, robust

* Weighted least squares nếu cần
regress read age mem iq [aweight=1/variance_weight]

* Stepwise regression
stepwise, pr(0.05) pe(0.10): regress read age mem iq

* Thống kê mô tả theo nhóm nếu có biến categorical
* egen age_group = cut(age), group(3)
* by age_group: summarize read

# Hồi quy tuyến tính đa biến với R
# Sử dụng các gói car, lmtest và ggplot2

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

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

# Thống kê mô tả
summary(data[c("read", "age", "mem", "iq")])

# Ma trận tương quan
cor_matrix <- cor(data[c("read", "age", "mem", "iq")])
print(cor_matrix)
corrplot(cor_matrix, method = "circle")

# Hồi quy tuyến tính đa biến
model <- lm(read ~ age + mem + iq, data = data)
summary(model)

# Kiểm tra các giả thuyết
# 1. Kiểm tra đa cộng tuyến
vif_values <- vif(model)
print(vif_values)

# 2. Kiểm tra phương sai đồng nhất
# Breusch-Pagan test
bptest(model)

# 3. Kiểm tra tự tương quan
# Durbin-Watson test
dwtest(model)

# 4. Kiểm tra phân phối chuẩn của phần dư
# Shapiro-Wilk test
shapiro.test(residuals(model))

# 5. Đồ thị chẩn đoán
par(mfrow = c(2, 2))
plot(model)

# 6. Kiểm tra tính tuyến tính
# Scatter plots
ggplot(data, aes(x = age, y = read)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(title = "Mối quan hệ giữa Age và Read")

ggplot(data, aes(x = mem, y = read)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(title = "Mối quan hệ giữa Memory và Read")

ggplot(data, aes(x = iq, y = read)) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(title = "Mối quan hệ giữa IQ và Read")

# 7. Kiểm tra điểm dị biệt
# Cook's distance
cooksd <- cooks.distance(model)
plot(cooksd, pch = "*", cex = 2, main = "Cook's Distance")
abline(h = 4/length(cooksd), col = "red")

# Leverage
leverage <- hatvalues(model)
plot(leverage, pch = "*", cex = 2, main = "Leverage")

# Standardized residuals
std_residuals <- rstandard(model)
plot(std_residuals, pch = "*", cex = 2, main = "Standardized Residuals")

# 8. Robust standard errors (nếu cần)
library(sandwich)
coeftest(model, vcov = vcovHC(model, type = "HC0"))

# 9. Weighted Least Squares (nếu cần)
# weights <- 1/fitted(lm(abs(residuals(model)) ~ fitted(model)))^2
# model_wls <- lm(read ~ age + mem + iq, data = data, weights = weights)

# 10. Stepwise regression
model_stepwise <- step(model, direction = "both")
summary(model_stepwise)

# 11. Tạo bảng kết quả đẹp
stargazer(model, type = "text", 
          title = "Kết quả hồi quy tuyến tính đa biến",
          dep.var.labels = "Khả năng đọc hiểu",
          covariate.labels = c("Tuổi", "Khả năng ghi nhớ", "Chỉ số IQ"))

# 12. Dự báo
new_data <- data.frame(age = 10, mem = 50, iq = 100)
prediction <- predict(model, new_data, interval = "confidence")
print(prediction)

# Hồi quy tuyến tính đa biến với Python
# Sử dụng scikit-learn, 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 sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score, mean_squared_error
from sklearn.preprocessing import StandardScaler
from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan, linear_rainbow
from statsmodels.stats.stattools import durbin_watson

# Đọc dữ liệu
data = pd.read_csv('ols_data.csv')
print(data.info())
print(data.describe())

# Kiểm tra missing values
print(data.isnull().sum())

# Thống kê mô tả
variables = ['read', 'age', 'mem', 'iq']
desc_stats = data[variables].describe()
print(desc_stats)

# Ma trận tương quan
correlation_matrix = data[variables].corr()
print(correlation_matrix)

# Vẽ heatmap 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')
plt.show()

# Scatter plots để kiểm tra tính tuyến tính
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

sns.scatterplot(data=data, x='age', y='read', ax=axes[0])
sns.regplot(data=data, x='age', y='read', ax=axes[0], scatter=False)
axes[0].set_title('Age vs Read')

sns.scatterplot(data=data, x='mem', y='read', ax=axes[1])
sns.regplot(data=data, x='mem', y='read', ax=axes[1], scatter=False)
axes[1].set_title('Memory vs Read')

sns.scatterplot(data=data, x='iq', y='read', ax=axes[2])
sns.regplot(data=data, x='iq', y='read', ax=axes[2], scatter=False)
axes[2].set_title('IQ vs Read')

plt.tight_layout()
plt.show()

# Chuẩn bị dữ liệu
X = data[['age', 'mem', 'iq']]
y = data['read']

# Hồi quy với scikit-learn
model_sklearn = LinearRegression()
model_sklearn.fit(X, y)

# Kết quả cơ bản
print("Sklearn Results:")
print(f"R² Score: {model_sklearn.score(X, y):.4f}")
print(f"Intercept: {model_sklearn.intercept_:.4f}")
print("Coefficients:")
for i, coef in enumerate(model_sklearn.coef_):
    print(f"  {X.columns[i]}: {coef:.4f}")

# Hồi quy với statsmodels (chi tiết hơn)
X_with_const = sm.add_constant(X)
model_sm = sm.OLS(y, X_with_const).fit()
print("\nStatsmodels Results:")
print(model_sm.summary())

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

# 1. Kiểm tra đa cộng tuyến (VIF)
print("\n1. Variance Inflation Factor (VIF):")
vif_data = pd.DataFrame()
vif_data["Variable"] = X.columns
vif_data["VIF"] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(vif_data)

# 2. Kiểm tra phương sai đồng nhất (Breusch-Pagan test)
print("\n2. Breusch-Pagan Test (Homoscedasticity):")
bp_test = het_breuschpagan(model_sm.resid, X_with_const)
bp_labels = ['LM Statistic', 'LM-Test p-value', 'F-Statistic', 'F-Test p-value']
print(dict(zip(bp_labels, bp_test)))

# 3. Kiểm tra tự tương quan (Durbin-Watson)
print("\n3. Durbin-Watson Test:")
dw_statistic = durbin_watson(model_sm.resid)
print(f"Durbin-Watson statistic: {dw_statistic:.4f}")

# 4. Kiểm tra phân phối chuẩn của phần dư
print("\n4. Shapiro-Wilk Test (Normality):")
shapiro_stat, shapiro_p = stats.shapiro(model_sm.resid)
print(f"Shapiro-Wilk statistic: {shapiro_stat:.4f}, p-value: {shapiro_p:.4f}")

# 5. Đồ thị chẩn đoán
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Residuals vs Fitted
fitted_values = model_sm.fittedvalues
residuals = model_sm.resid

axes[0, 0].scatter(fitted_values, 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_values, 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()

# 6. Kiểm tra điểm dị biệt
print("\n6. Outlier Detection:")

# Cook's distance
influence = model_sm.get_influence()
cook_distance = influence.cooks_distance[0]
leverage = influence.hat_matrix_diag
dfbetas = influence.dfbetas

print(f"Max Cook's Distance: {cook_distance.max():.4f}")
print(f"Threshold (4/n): {4/len(data):.4f}")

# Vẽ đồ thị Cook's distance
plt.figure(figsize=(10, 6))
plt.subplot(1, 2, 1)
plt.scatter(range(len(cook_distance)), cook_distance)
plt.axhline(y=4/len(data), color='red', linestyle='--', label='Threshold')
plt.xlabel('Observation')
plt.ylabel("Cook's Distance")
plt.title("Cook's Distance")
plt.legend()

# Vẽ đồ thị Leverage
plt.subplot(1, 2, 2)
plt.scatter(range(len(leverage)), leverage)
plt.axhline(y=2*X.shape[1]/len(data), color='red', linestyle='--', label='Threshold')
plt.xlabel('Observation')
plt.ylabel('Leverage')
plt.title('Leverage')
plt.legend()

plt.tight_layout()
plt.show()

# 7. Dự báo
print("\n7. Prediction Example:")
new_data = pd.DataFrame({
    'age': [10],
    'mem': [50], 
    'iq': [100]
})

prediction_sklearn = model_sklearn.predict(new_data)
print(f"Sklearn prediction: {prediction_sklearn[0]:.2f}")

new_data_with_const = sm.add_constant(new_data)
prediction_sm = model_sm.predict(new_data_with_const)
confidence_interval = model_sm.get_prediction(new_data_with_const).conf_int()
print(f"Statsmodels prediction: {prediction_sm[0]:.2f}")
print(f"95% Confidence Interval: [{confidence_interval[0,0]:.2f}, {confidence_interval[0,1]:.2f}]")

# 8. Model comparison và feature selection
from sklearn.feature_selection import f_regression

print("\n8. Feature Importance:")
f_stats, f_pvalues = f_regression(X, y)
feature_importance = pd.DataFrame({
    'Feature': X.columns,
    'F-statistic': f_stats,
    'p-value': f_pvalues
}).sort_values('F-statistic', ascending=False)

print(feature_importance)

# 9. Robust regression (nếu cần)
print("\n9. Robust Regression Results:")
robust_model = sm.RLM(y, X_with_const).fit()
print(robust_model.summary())