Hồi quy tuyến tính giản đơn

* Thực hiện hồi quy tuyến tính giản đơn
REGRESSION
  /MISSING LISTWISE
  /STATISTICS COEFF OUTS CI(95) R ANOVA COLLIN TOL CHANGE ZPP
  /CRITERIA=PIN(.05) POUT(.10)
  /NOORIGIN 
  /DEPENDENT y
  /METHOD=ENTER x
  /SCATTERPLOT=(*ZRESID ,*ZPRED)
  /RESIDUALS NORMPROB(ZRESID)
  /RESIDUALS HISTOGRAM(ZRESID)
  /SAVE PRED COOK RESID ZRESID ZPRED.

* Kiểm tra thống kê mô tả
DESCRIPTIVES VARIABLES=x y
  /STATISTICS=MEAN STDDEV MIN MAX.

* Tạo đồ thị phân tán
GRAPH
  /SCATTERPLOT(BIVAR)=x WITH y
  /MISSING=LISTWISE.

* Đọc dữ liệu
use "simple-ols.dta", clear

* Thống kê mô tả
summarize x y

* Tạo đồ thị phân tán
scatter y x, title("Mối quan hệ giữa kinh nghiệm và thu nhập") ///
    xtitle("Số năm kinh nghiệm") ytitle("Thu nhập (triệu đồng/năm)")

* Thực hiện hồi quy tuyến tính giản đơn
regress y x

* Kiểm tra giả thiết
* Durbin-Watson test
estat dwatson

* Kiểm tra phân phối chuẩn của phần dư
predict residuals, residuals
histogram residuals, normal title("Phân phối phần dư")
qnorm residuals

* Kiểm tra phương sai đồng nhất
predict fitted, xb
scatter residuals fitted, title("Đồ thị phần dư vs giá trị dự báo")

* Đồ thị regression line
twoway (scatter y x) (lfit y x), ///
    title("Đường hồi quy tuyến tính") ///
    legend(label(1 "Dữ liệu quan sát") label(2 "Đường hồi quy"))

* Kiểm định tính ổn định của mô hình
ovtest

* Tính confidence intervals
regress y x, level(95)

# Load thư viện cần thiết
library(ggplot2)
library(car)
library(lmtest)

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

# Thống kê mô tả
summary(data$x)
summary(data$y)

# Tạo đồ thị phân tán
ggplot(data, aes(x = x, y = y)) +
  geom_point() +
  geom_smooth(method = "lm", se = TRUE) +
  labs(title = "Mối quan hệ giữa kinh nghiệm và thu nhập",
       x = "Số năm kinh nghiệm",
       y = "Thu nhập (triệu đồng/năm)")

# Thực hiện hồi quy tuyến tính giản đơn
model <- lm(y ~ x, data = data)
summary(model)

# Tính confidence intervals
confint(model)

# Kiểm tra giả thiết
# Durbin-Watson test
dwtest(model)

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

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

# Histogram phần dư
hist(residuals(model), breaks = 10, main = "Histogram phần dư", 
     xlab = "Phần dư", ylab = "Tần suất")

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

# Kiểm tra phương sai đồng nhất
bptest(model)

# Đồ thị phần dư vs fitted values
plot(fitted(model), residuals(model),
     main = "Phần dư vs Giá trị dự báo",
     xlab = "Giá trị dự báo", ylab = "Phần dư")
abline(h = 0, col = "red")

# Tính R-squared
summary(model)$r.squared
summary(model)$adj.r.squared

# Dự báo với dữ liệu mới
new_data <- data.frame(x = c(5, 10, 15))
predictions <- predict(model, new_data, interval = "prediction")
print(predictions)

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
import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.stattools import durbin_watson

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

# Thống kê mô tả
print("Thống kê mô tả:")
print(data[['x', 'y']].describe())

# Tạo đồ thị phân tán
plt.figure(figsize=(10, 6))
plt.scatter(data['x'], data['y'], alpha=0.6)
plt.plot(np.unique(data['x']), np.poly1d(np.polyfit(data['x'], data['y'], 1))(np.unique(data['x'])), 'r--')
plt.xlabel('Số năm kinh nghiệm')
plt.ylabel('Thu nhập (triệu đồng/năm)')
plt.title('Mối quan hệ giữa kinh nghiệm và thu nhập')
plt.grid(True, alpha=0.3)
plt.show()

# Chuẩn bị dữ liệu
X = data[['x']]
y = data['y']

# Thực hiện hồi quy tuyến tính với statsmodels
X_with_const = sm.add_constant(X)
model = sm.OLS(y, X_with_const).fit()
print("\nKết quả hồi quy tuyến tính:")
print(model.summary())

# Tính R-squared
r_squared = model.rsquared
adj_r_squared = model.rsquared_adj
print(f"\nR-squared: {r_squared:.4f}")
print(f"Adjusted R-squared: {adj_r_squared:.4f}")

# Kiểm tra giả thiết
residuals = model.resid
fitted_values = model.fittedvalues

# Durbin-Watson test
dw_stat = durbin_watson(residuals)
print(f"\nDurbin-Watson statistic: {dw_stat:.4f}")

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

# Kiểm tra phương sai đồng nhất
bp_stat, bp_p, _, _ = het_breuschpagan(residuals, X_with_const)
print(f"Breusch-Pagan test: statistic={bp_stat:.4f}, p-value={bp_p:.4f}")

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

# Histogram phần dư
axes[0, 0].hist(residuals, bins=10, edgecolor='black', alpha=0.7)
axes[0, 0].set_title('Histogram phần dư')
axes[0, 0].set_xlabel('Phần dư')
axes[0, 0].set_ylabel('Tần suất')

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

# Phần dư vs fitted values
axes[1, 0].scatter(fitted_values, residuals, alpha=0.6)
axes[1, 0].axhline(y=0, color='red', linestyle='--')
axes[1, 0].set_title('Phần dư vs Giá trị dự báo')
axes[1, 0].set_xlabel('Giá trị dự báo')
axes[1, 0].set_ylabel('Phần dư')

# Leverage plot
leverage = model.get_influence().hat_matrix_diag
axes[1, 1].scatter(leverage, residuals, alpha=0.6)
axes[1, 1].set_title('Leverage Plot')
axes[1, 1].set_xlabel('Leverage')
axes[1, 1].set_ylabel('Phần dư')

plt.tight_layout()
plt.show()

# Dự báo với dữ liệu mới
new_x = np.array([5, 10, 15]).reshape(-1, 1)
new_x_with_const = sm.add_constant(new_x)
predictions = model.predict(new_x_with_const)
pred_intervals = model.get_prediction(new_x_with_const).conf_int()

print("\nDự báo cho giá trị mới:")
for i, x_val in enumerate([5, 10, 15]):
    print(f"Kinh nghiệm {x_val} năm: {predictions[i]:.2f} triệu đồng")
    print(f"  Confidence interval: [{pred_intervals[i, 0]:.2f}, {pred_intervals[i, 1]:.2f}]")

# Hiển thị phương trình hồi quy
intercept = model.params[0]
slope = model.params[1]
print(f"\nPhương trình hồi quy: y = {intercept:.3f} + {slope:.3f}*x")