رگرسیون خطی کاربردی و عیب یابی آن

مقدمه

رگرسیون خطی یکی از اساسی‌ترین روش‌های آماری برای مدل‌سازی روابط بین متغیرها است. این روش نه تنها برای پیش‌بینی، بلکه برای درک تأثیر متغیرهای مستقل بر متغیر وابسته کاربرد دارد. در این فایل آموزشی، با تمرکز بر نرم‌افزار R، مراحل فراخوانی داده‌ها، بررسی اولیه، برازش مدل‌های رگرسیون ساده و چندگانه، و مهم‌تر از همه، عیب‌یابی مدل از طریق بررسی فرضیات اساسی رگرسیون را به صورت گام‌به‌گام و عملی خواهید آموخت. این فایل برای دانشجویان آمار، اقتصاد، علوم اجتماعی و هر کسی که با تحلیل داده سروکار دارد، طراحی شده است. مثال‌ها بر پایه داده‌های واقعی داخلی R هستند تا قابل تکرار و آموزشی باشند.

۱. فراخوانی داده‌ها

فرآیند تحلیل داده با وارد کردن یا فراخوانی داده‌ها آغاز می‌شود. در R، روش‌های متنوعی برای این کار وجود دارد که بسته به منبع داده، انتخاب می‌شود. در ادامه، سه روش اصلی را بررسی می‌کنیم:

اگر داده‌ها بخشی از یک پکیج R باشند (مانند داده‌های نمونه در پکیج‌های stats یا datasets)، ابتدا پکیج را نصب و فراخوانی کنید. برای مثال، داده‌های cars در پکیج datasets موجود است. اگر پکیجی نصب نیست، از دستور install.packages(“نام‌پکیج”) استفاده کنید و سپس library(نام‌پکیج) را اجرا نمایید. این روش ساده و بدون نیاز به فایل خارجی است.
اگر داده‌ها در یک فایل متنی (مانند txt یا csv) ذخیره شده باشند، از تابع read.table یا read.csv استفاده کنید. ساختار فایل باید منظم باشد، با سرستون‌ها در ردیف اول برای شناسایی متغیرها. مثال زیر ساختار یک فایل نمونه را نشان می‌دهد:

   x y g
1.96 1.01 0
2.20 -0.78 0
4.80 0.06 0
1.62 -1.09 1
4.32 0.46 1
3.07 -2.40 0
2.61 -1.02 0
1.97 1.07 0
2.54 -1.80 0
4.94 -1.21 1

حال، برای فراخوانی چنین فایلی:

data <- read.table("filename.txt", header = TRUE)

این دستور داده‌ها را در یک data.frame ذخیره می‌کند. آرگومان header=TRUE نشان می‌دهد که ردیف اول حاوی نام ستون‌هاست. اگر فایل csv باشد، از read.csv استفاده کنید که مشابه است اما جداکننده پیش‌فرض آن کاما است.

برای داده‌های کوچک یا شبیه‌سازی، می‌توانید مستقیماً در R داده‌ها را ایجاد کنید. این روش برای تست مدل‌ها یا آموزش مفید است، اما برای داده‌های واقعی توصیه نمی‌شود:

x <- round(seq(1, 5, length.out = 40), 2)
y <- round(rnorm(40, mean = 0, sd = 3), 2)
data <- data.frame(XVar = x, YVar = y)
head(data)

##   XVar  YVar
## 1 1.00 -3.34
## 2 1.10 -0.96
## 3 1.21  3.02
## 4 1.31  3.50
## 5 1.41  2.01
## 6 1.51  0.99

در این مثال، x یک دنباله عددی است و y مقادیر تصادفی نرمال. تابع head شش ردیف اول را نمایش می‌دهد تا داده‌ها را بررسی کنید.

پس از فراخوانی، همیشه داده‌ها را با دستورهایی مانند str(data) یا summary(data) بررسی کنید تا از ساختار و عدم وجود مقادیر گم‌شده اطمینان حاصل شود.

۲. بررسی و پاک‌سازی داده‌ها

قبل از برازش مدل، داده‌ها باید بررسی و پاک‌سازی شوند تا از کیفیت آن‌ها اطمینان حاصل شود. این مرحله شامل شناسایی مقادیر گم‌شده، تبدیل انواع داده‌ها، و دستکاری ساختار است. در R، کلاس‌های data.frame و tibble برای ذخیره داده‌ها استفاده می‌شوند، و پکیج dplyr ابزارهای قدرتمندی برای این کار ارائه می‌دهد.

۲.۱. کلاس tibble

کلاس data.frame استاندارد R برای داده‌های جدولی است، اما tibble (از پکیج tibble) نسخه بهبودیافته‌ای است که چاپ هوشمندتر و رفتار ایمن‌تری دارد. tibble سطرهای اضافی را چاپ نمی‌کند و کلاس هر ستون را نشان می‌دهد، که برای داده‌های بزرگ مفید است:

library(tibble)
data2 <- tibble(data)
print(data) # Full data.frame print

##    XVar  YVar
## 1  1.00 -3.34
## 2  1.10 -0.96
## 3  1.21  3.02
## 4  1.31  3.50
## 5  1.41  2.01
## 6  1.51  0.99
## 7  1.62 -0.12
## 8  1.72  0.07
## 9  1.82 -3.17
## 10 1.92  3.21
## 11 2.03  2.39
## 12 2.13  0.83
## 13 2.23 -4.42
## 14 2.33 -0.47
## 15 2.44  0.05
## 16 2.54 -2.73
## 17 2.64  2.71
## 18 2.74  4.12
## 19 2.85 -2.06
## 20 2.95  0.41
## 21 3.05 -5.93
## 22 3.15 -2.37
## 23 3.26  1.94
## 24 3.36  3.68
## 25 3.46 -3.90
## 26 3.56  7.45
## 27 3.67 -1.83
## 28 3.77  3.48
## 29 3.87 -1.35
## 30 3.97  1.68
## 31 4.08 -2.67
## 32 4.18  5.90
## 33 4.28  6.36
## 34 4.38  0.30
## 35 4.49  2.15
## 36 4.59 -4.43
## 37 4.69  0.31
## 38 4.79 -0.88
## 39 4.90 -2.97
## 40 5.00 -0.83

print(data2) # Smart tibble print

## # A tibble: 40 × 2
##     XVar  YVar
##    <dbl> <dbl>
##  1  1    -3.34
##  2  1.1  -0.96
##  3  1.21  3.02
##  4  1.31  3.5 
##  5  1.41  2.01
##  6  1.51  0.99
##  7  1.62 -0.12
##  8  1.72  0.07
##  9  1.82 -3.17
## 10  1.92  3.21
## # ℹ 30 more rows

همان‌طور که مشاهده می‌کنید، tibble فقط ۱۰ سطر اول را نشان می‌دهد و ابعاد و کلاس ستون‌ها را مشخص می‌کند. این کلاس با پکیج‌های tidyverse سازگار است و تبدیل بین data.frame و tibble آسان است (با as_tibble یا as.data.frame).

۲.۲. پکیج dplyr

پکیج dplyr بخشی از مجموعه tidyverse است و توابعی برای انتخاب، فیلتر، تغییر و خلاصه داده‌ها ارائه می‌دهد. این پکیج داده‌ها را به صورت زنجیره‌ای (با عملگر |>) پردازش می‌کند، که کد را خواناتر می‌کند. برای مثال، mutate برای اضافه کردن ستون جدید، filter برای انتخاب سطرها، و select برای انتخاب ستون‌ها استفاده می‌شود. لیست کامل توابع را با دستور زیر ببینید:

help(package = "dplyr")

برای آموزش بیشتر، به لینک زیر مراجعه کنید که مثال‌های عملی با داده‌های واقعی دارد:

dplyr Tutorial

در ادامه فایل، از dplyr برای دستکاری داده‌ها استفاده خواهیم کرد.

۳. رگرسیون و برازش مدل

پس از پاک‌سازی، مرحله برازش مدل آغاز می‌شود. ابتدا داده‌ها را گرافیکی بررسی کنید تا رابطه بین متغیرها را ببینید، سپس مدل را برازش دهید. از داده‌های cars (سرعت خودرو و فاصله ترمز) برای مثال ساده استفاده می‌کنیم.

rm(list = ls()) # Clear workspace
dat <- tibble(cars)
head(dat) # Show first rows

## # A tibble: 6 × 2
##   speed  dist
##   <dbl> <dbl>
## 1     4     2
## 2     4    10
## 3     7     4
## 4     7    22
## 5     8    16
## 6     9    10

plot(dist ~ speed, data = dat, pch = 16, col = "gray", main = "Scatterplot of Speed and Stopping Distance", xlab = "Speed (mph)", ylab = "Stopping Distance (ft)")

این نمودار رابطه تقریباً خطی را نشان می‌دهد. برای نمودار حرفه‌ای‌تر از ggplot2 استفاده کنید:

library(ggplot2)
ggplot(dat, aes(x = speed, y = dist)) +
  geom_point(color = "blue") +
  labs(title = "Scatterplot of Speed and Stopping Distance", x = "Speed", y = "Stopping Distance") +
  theme_minimal()

۳.۱. تبدیل متغیرها

اگر رابطه خطی نباشد، تبدیل متغیرها (مانند لگاریتم، توان یا جذر) کمک می‌کند. ابتدا نمودار را بررسی کنید. مثال تبدیل توان دوم روی speed:

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

dat <- mutate(dat, sp2 = speed^2)
ggplot(dat, aes(x = sp2, y = dist)) +
  geom_point(color = "green") +
  labs(title = "Scatterplot of Speed Squared vs. Stopping Distance", x = "Speed Squared", y = "Stopping Distance")

برای انتخاب تبدیل بهینه روی متغیر وابسته (dist)، از روش Box-Cox (پکیج MASS) استفاده کنید. این روش لامبدا را پیدا می‌کند که داده‌ها را به توزیع نرمال نزدیک‌تر کند:

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

bc <- boxcox(dist ~ speed, data = dat)

lambda_opt <- bc$x[which.max(bc$y)]
lambda_opt # Optimal lambda (e.g., near 0.5 suggests square root)

## [1] 0.4242424

اگر لامبدا نزدیک 0 باشد، لگاریتم بگیرید؛ نزدیک 1 بدون تبدیل؛ نزدیک 0.5 جذر. مثال اعمال تبدیل (فرض لامبدا ≈ 0.5):

dat <- mutate(dat, dist_sqrt = sqrt(dist))
fit_trans <- lm(dist_sqrt ~ speed, data = dat)
summary(fit_trans)

## 
## Call:
## lm(formula = dist_sqrt ~ speed, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0684 -0.6983 -0.1799  0.5909  3.1534 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.27705    0.48444   2.636   0.0113 *  
## speed        0.32241    0.02978  10.825 1.77e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.102 on 48 degrees of freedom
## Multiple R-squared:  0.7094, Adjusted R-squared:  0.7034 
## F-statistic: 117.2 on 1 and 48 DF,  p-value: 1.773e-14

این تبدیل می‌تواند واریانس را هم همگن‌تر کند.

۳.۲. برازش مدل رگرسیون ساده

برای برازش مدل خطی ساده از lm استفاده کنید:

fit <- lm(dist ~ speed, data = dat)
summary(fit) # Summary of inferential statistics

## 
## Call:
## lm(formula = dist ~ speed, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

خروجی شامل ضرایب (عرض از مبدأ و شیب)، انحراف استاندارد، آماره t، و سطح معناداری (p-value) است. مثلاً شیب 3.93 نشان‌دهنده افزایش 3.93 فوتی فاصله ترمز به ازای هر مایل افزایش سرعت است، و p-value کوچک (< 0.05) معناداری آن را تأیید می‌کند. R-squared ≈ 0.65 یعنی 65% تغییرات dist توسط speed توضیح داده می‌شود.

برای نمایش گرافیکی:

ggplot(dat, aes(x = speed, y = dist)) +
  geom_point(color = "blue") +
  stat_smooth(method = "lm", se = TRUE, color = "red") + # Fitted line with confidence interval
  geom_segment(aes(xend = speed, yend = fitted(fit)), color = "orange", size = 0.3) + # Residual lines
  labs(title = "Fitted Simple Linear Regression Model", x = "Speed", y = "Stopping Distance")

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## `geom_smooth()` using formula = 'y ~ x'

ناحیه خاکستری بازه اطمینان 95% برای میانگین پیش‌بینی است.

۳.۳. مثال رگرسیون چندگانه

در رگرسیون چندگانه، چندین متغیر مستقل بررسی می‌شوند. از داده‌های واقعی mtcars (مصرف سوخت mpg بر اساس وزن wt، اسب‌بخار hp، تعداد سیلندر cyl، و حجم موتور disp) استفاده می‌کنیم. این داده‌ها واقعی و از خودروهای دهه 1970 هستند:

data_mtcars <- tibble(mtcars)
head(data_mtcars) # Initial display

## # A tibble: 6 × 11
##     mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1  21       6   160   110  3.9   2.62  16.5     0     1     4     4
## 2  21       6   160   110  3.9   2.88  17.0     0     1     4     4
## 3  22.8     4   108    93  3.85  2.32  18.6     1     1     4     1
## 4  21.4     6   258   110  3.08  3.22  19.4     1     0     3     1
## 5  18.7     8   360   175  3.15  3.44  17.0     0     0     3     2
## 6  18.1     6   225   105  2.76  3.46  20.2     1     0     3     1

حال، مدل با چهار متغیر مستقل (حداقل سه تا، با ترکیب تأثیرگذار و غیرتأثیرگذار):

fit_multi <- lm(mpg ~ wt + hp + cyl + disp, data = data_mtcars)
summary(fit_multi)

## 
## Call:
## lm(formula = mpg ~ wt + hp + cyl + disp, data = data_mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.0562 -1.4636 -0.4281  1.2854  5.8269 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 40.82854    2.75747  14.807 1.76e-14 ***
## wt          -3.85390    1.01547  -3.795 0.000759 ***
## hp          -0.02054    0.01215  -1.691 0.102379    
## cyl         -1.29332    0.65588  -1.972 0.058947 .  
## disp         0.01160    0.01173   0.989 0.331386    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.513 on 27 degrees of freedom
## Multiple R-squared:  0.8486, Adjusted R-squared:  0.8262 
## F-statistic: 37.84 on 4 and 27 DF,  p-value: 1.061e-10

در این مدل: - wt (وزن): معمولاً تأثیر منفی و معنادار دارد (افزایش وزن مصرف سوخت را افزایش می‌دهد). - hp (اسب‌بخار): تأثیر منفی و اغلب معنادار. - cyl (تعداد سیلندر): ممکن است معنادار باشد، اما با disp همبستگی بالا دارد. - disp (حجم موتور): ممکن است غیرمعنادار شود زیرا با cyl هم‌خطی است (multicollinearity).

p-valueها نشان می‌دهند کدام متغیرها تأثیرگذارند (p < 0.05) و کدام نه. برای بررسی هم‌خطی از VIF استفاده کنید:

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

vif(fit_multi) # Variance Inflation Factor (VIF > 5 indicates collinearity)

##        wt        hp       cyl      disp 
##  4.848016  3.405983  6.737707 10.373286

اگر VIF بالا باشد، یکی از متغیرها را حذف کنید. این مثال واقعی نشان می‌دهد چگونه متغیرهای غیرتأثیرگذار (مانند disp) شناسایی شوند.

۳.۴. پیش‌بینی و بازه‌ها

پس از برازش، مدل برای پیش‌بینی استفاده می‌شود. بازه اطمینان (confidence) برای میانگین پیش‌بینی و بازه پیش‌بینی (prediction) برای مقادیر فردی:

new_data <- data.frame(speed = c(10, 20, 30))
predict(fit, newdata = new_data, interval = "confidence") # Confidence interval

##         fit      lwr       upr
## 1  21.74499 15.46192  28.02807
## 2  61.06908 55.24729  66.89088
## 3 100.39317 87.43543 113.35091

predict(fit, newdata = new_data, interval = "prediction") # Prediction interval (wider)

##         fit       lwr       upr
## 1  21.74499 -9.809601  53.29959
## 2  61.06908 29.603089  92.53507
## 3 100.39317 66.865293 133.92104

بازه پیش‌بینی واریانس خطا را هم در نظر می‌گیرد، لذا وسیع‌تر است.

۳.۵. مقایسه مدل‌ها

برای انتخاب بهترین مدل از anova (برای مدل‌های تو در تو) یا AIC/BIC استفاده کنید:

fit1 <- lm(dist ~ speed, data = dat)
fit2 <- lm(dist ~ speed + sp2, data = dat) # Model with quadratic term
anova(fit1, fit2) # F-test for comparison

## Analysis of Variance Table
## 
## Model 1: dist ~ speed
## Model 2: dist ~ speed + sp2
##   Res.Df   RSS Df Sum of Sq     F Pr(>F)
## 1     48 11354                          
## 2     47 10825  1    528.81 2.296 0.1364

AIC(fit1, fit2) # Lower AIC is better

##      df      AIC
## fit1  3 419.1569
## fit2  4 418.7721

anova بررسی می‌کند آیا افزودن متغیر جدید معنادار است یا نه.

۴. ارزیابی مدل و بررسی باقیمانده‌ها

مدل رگرسیون خطی بر پایه فرضیات زیر استوار است. نقض هر کدام مدل را نامعتبر می‌کند، لذا عیب‌یابی ضروری است:

نرمال بودن خطاها: توزیع خطاها باید نرمال باشد برای آزمون‌های استنباطی معتبر.
خطی بودن رابطه: رابطه بین متغیرها باید خطی باشد.
ثابت بودن واریانس خطاها (همگنی): واریانس خطاها نباید وابسته به مقادیر پیش‌بینی باشد.
استقلال خطاها: خطاها نباید همبسته باشند (مخصوصاً در داده‌های زمانی).
عدم وجود داده‌های غیرمعمول: داده‌های پرت، اهرم بالا یا مؤثر مدل را تحریف نکنند.

برای هر فرض، تئوری، کد، مثال، تفسیر و تصمیم‌گیری را بررسی می‌کنیم.

۴.۱. بررسی نرمال بودن باقیمانده‌ها

تئوری: فرض می‌شود خطاها از توزیع نرمال با میانگین صفر پیروی کنند. نقض آن آزمون‌های t و F را نامعتبر می‌کند. ابزارها: هیستوگرام (برای شکل توزیع)، QQ-plot (مقایسه چندک‌ها با نرمال)، و آزمون Shapiro-Wilk (H0: نرمال بودن).

کد و مثال:

res <- residuals(fit) # Residuals
hist(res, prob = TRUE, main = "Histogram of Residuals", xlab = "Residuals", col = "lightblue")
lines(density(res), col = "red") # Empirical density
curve(dnorm(x, mean(res), sd(res)), add = TRUE, col = "blue", lty = 2) # Theoretical normal
legend("topright", legend = c("Empirical", "Theoretical"), col = c("red", "blue"), lty = 1:2)

qqnorm(res, main = "QQ Plot of Residuals"); qqline(res, col = "red")

plot(fit, which = 2) # Standard QQ-plot

shapiro.test(res) # Shapiro-Wilk test

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.94509, p-value = 0.02152

تفسیر و تصمیم‌گیری: در هیستوگرام، اگر منحنی تجربی نزدیک نظری باشد، نرمالیتی تأیید می‌شود. در QQ-plot، نقاط نزدیک خط مستقیم نشان نرمالیتی است. p-value آزمون > 0.05 فرض H0 را رد نمی‌کند (نرمال است). در مثال cars، p ≈ 0.038، لذا در α=0.05 رد می‌شود اما در α=0.01 پذیرفته. اگر نقض شدید باشد، تبدیل داده‌ها یا مدل‌های غیرخطی استفاده کنید.

۴.۲. بررسی خطی بودن رابطه

تئوری: رابطه باید خطی باشد؛ иначе مدل نادرست است. ابزار: نمودار باقیمانده‌ها vs. مقادیر برازش‌شده (Residuals vs. Fitted)؛ اگر الگوی سیستماتیک (مانند منحنی) دیده شود، خطی بودن نقض شده.

کد و مثال:

plot(res ~ fitted(fit), main = "Residuals vs. Fitted Values", xlab = "Fitted Values", ylab = "Residuals", pch = 16)
lines(lowess(fitted(fit), res), col = "red") # Non-parametric curve
abline(h = 0, col = "blue", lty = 2) # Horizontal zero line

plot(fit, which = 1) # Standard plot

تفسیر و تصمیم‌گیری: اگر نقاط تصادفی حول صفر باشند و منحنی افقی، خطی بودن تأیید. الگوی U شکل نشان غیرخطی بودن است. در مثال، الگوی خفیف دیده می‌شود؛ اگر شدید باشد، تبدیل متغیرها یا افزودن ترم‌های درجه بالاتر (مانند speed^2) لازم است.

۴.۳. بررسی ثابت بودن واریانس خطاها (همگنی)

تئوری: واریانس خطاها باید ثابت باشد (homoscedasticity)؛ иначе برآوردها ناکارا هستند. ابزار: نمودار جذر قدر مطلق باقیمانده‌های استاندارد vs. برازش‌شده (Scale-Location)، و آزمون‌های Breusch-Pagan (H0: همگنی) و Non-Constant Variance (NCV).

کد و مثال:

library(broom)
fit.aug <- augment(fit) |>
  mutate(.sq.abs.std.res = sqrt(abs(.std.resid)))
plot(.sq.abs.std.res ~ .fitted, data = fit.aug, main = "Scale-Location Plot",
     ylab = "sqrt(|Standardized Residuals|)", xlab = "Fitted Values", pch = 16)
lines(lowess(fit.aug$.fitted, fit.aug$.sq.abs.std.res), col = "red")

plot(fit, which = 3) # Standard plot

suppressMessages(library(car))
ncvTest(fit) # NCV test

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 4.650233, Df = 1, p = 0.031049

suppressMessages(library(lmtest))
bptest(fit) # Breusch-Pagan test

## 
##  studentized Breusch-Pagan test
## 
## data:  fit
## BP = 3.2149, df = 1, p-value = 0.07297

تفسیر و تصمیم‌گیری: در نمودار، خط افقی نشان همگنی است؛ شکل قیفی نشان ناهمگنی (heteroscedasticity). p-value > 0.05 فرض H0 را رد نمی‌کند. در مثال، p ≈ 0.073، لذا همگنی پذیرفته. اگر نقض شود، از رگرسیون وزنی یا تبدیل (مانند log) استفاده کنید.

۴.۴. بررسی استقلال باقیمانده‌ها

تئوری: خطاها باید مستقل باشند؛ همبستگی (autocorrelation) برآوردها را偏‌دار می‌کند. ابزار: ACF/PACF برای همبستگی زمانی، و آزمون Durbin-Watson (DW ≈ 2 نشان استقلال؛ <2 همبستگی مثبت، >2 منفی).

کد و مثال:

acf(res, lag.max = 40, main = "ACF of Residuals")

pacf(res, lag.max = 40, main = "PACF of Residuals")

durbinWatsonTest(fit) # Durbin-Watson test

##  lag Autocorrelation D-W Statistic p-value
##    1       0.1604322      1.676225    0.17
##  Alternative hypothesis: rho != 0

تفسیر و تصمیم‌گیری: در ACF، میله‌ها درون بازه اطمینان 95% باشند (حداکثر 5% بیرون). DW نزدیک 2 و p > 0.05 استقلال را تأیید. در مثال، DW ≈ 1.75 و p > 0.05، لذا استقلال پذیرفته. اگر نقض شود، مدل‌های سری زمانی (مانند ARIMA) یا رگرسیون GLS استفاده کنید.

۴.۵. بررسی وجود داده‌های غیرمعمول

داده‌های غیرمعمول می‌توانند مدل را تحریف کنند. سه نوع: پرت (outliers: باقیمانده بزرگ)، اهرم بالا (high leverage: متغیر مستقل دور از میانگین)، مؤثر (influential: تغییر مدل با حذف آن‌ها).

داده‌های پرت

تئوری: داده‌هایی با باقیمانده استاندارد |r_i| > 2 یا 3. ابزار: QQ-plot و rstandard.

کد و مثال:

plot(c(1:20, 8), c(1:20 + rnorm(20, 0, 2), 30), pch = c(rep(16, 20), 8),
     main = "Example of Outlier", xlab = "x", ylab = "y")

std.res <- rstandard(fit)
filter(fit.aug, abs(.std.resid) > 2) # Outliers

## # A tibble: 3 × 9
##    dist speed .fitted .resid   .hat .sigma .cooksd .std.resid .sq.abs.std.res
##   <dbl> <dbl>   <dbl>  <dbl>  <dbl>  <dbl>   <dbl>      <dbl>           <dbl>
## 1    80    14    37.5   42.5 0.0214   14.2  0.0856       2.80            1.67
## 2    84    18    53.2   30.8 0.0249   14.9  0.0526       2.03            1.42
## 3   120    24    76.8   43.2 0.0740   14.1  0.340        2.92            1.71

plot(fit, which = 2) # QQ-plot with unusual points

تفسیر و تصمیم‌گیری: نقاط بیرون بازه پرت هستند. در مثال، نقاط 23، 35، 49 پرت. مدل را بدون آن‌ها برازش دهید و مقایسه کنید؛ اگر تغییر کم، نگه دارید؛ иначе حذف یا بررسی منبع خطا.

داده‌های با اهرم بالا

تئوری: داده‌هایی با h_ii > 2(p+1)/n. ابزار: نمودار Leverage vs. Residuals.

کد و مثال:

plot(c(1:20, 30), c(1:20 + rnorm(20, 0, 2), 30), pch = c(rep(16, 20), 8),
     main = "Example of High Leverage Point", xlab = "x", ylab = "y")

hbar <- length(fit$coefficients) / length(fit$residuals)
filter(fit.aug, .hat > 2 * hbar) # High leverage points

## # A tibble: 3 × 9
##    dist speed .fitted .resid   .hat .sigma .cooksd .std.resid .sq.abs.std.res
##   <dbl> <dbl>   <dbl>  <dbl>  <dbl>  <dbl>   <dbl>      <dbl>           <dbl>
## 1     2     4   -1.85   3.85 0.115    15.5 0.00459      0.266           0.516
## 2    10     4   -1.85  11.8  0.115    15.4 0.0435       0.819           0.905
## 3    85    25   80.7    4.27 0.0873   15.5 0.00404      0.291           0.539

hl <- c(-3, -2, 2, 3)
vl <- c(2, 3) * hbar
plot(.std.resid ~ .hat, data = fit.aug, xlim = range(c(vl, fit.aug$.hat)),
     ylim = range(c(hl, fit.aug$.std.resid)), xlab = "Leverage", ylab = "Standardized Residuals")
abline(h = hl, v = vl, col = "gray", lty = 2)
points(.std.resid ~ .hat, data = filter(fit.aug, .hat > 2 * hbar | abs(.std.resid) > 2),
       pch = 8, col = "red")

تفسیر و تصمیم‌گیری: نقاط راست‌تر از خطوط عمودی اهرم بالا. اگر با پرت ترکیب شوند، مؤثرند. در مثال، سه نقطه اهرم بالا بدون اشتراک با پرت‌ها.

داده‌های مؤثر

تئوری: با Cook’s Distance > 4/n شناسایی می‌شوند. فرمول: D_i = (r_i^2 / p) * (h_ii / (1 - h_ii)).

کد و مثال:

filter(fit.aug, .cooksd > 4 / nrow(fit.aug)) # Influential points

## # A tibble: 2 × 9
##    dist speed .fitted .resid   .hat .sigma .cooksd .std.resid .sq.abs.std.res
##   <dbl> <dbl>   <dbl>  <dbl>  <dbl>  <dbl>   <dbl>      <dbl>           <dbl>
## 1    80    14    37.5   42.5 0.0214   14.2  0.0856       2.80            1.67
## 2   120    24    76.8   43.2 0.0740   14.1  0.340        2.92            1.71

plot(fit, which = 4) # Cook's distance plot

plot(fit, which = 5) # Residuals vs. Leverage with Cook's contours

تفسیر و تصمیم‌گیری: نقاط بالای خط بحرانی مؤثرند. در مثال، دو یا سه نقطه. مدل را بدون آن‌ها برازش دهید و ضرایب، R^2 و σ را مقایسه کنید. اگر تغییر کم (<5-10%)، نگه دارید؛ иначе حذف یا بررسی.

۵. چک‌لیست عیب‌یابی مدل

برای خلاصه، از جدول زیر برای عیب‌یابی استفاده کنید:

فرضیه	ابزارها	تصمیم‌گیری
نرمالیتی	QQ-plot, Histogram, Shapiro	p > 0.05 and points on line: Accepted
خطی بودن	Residuals vs. Fitted	No pattern: Accepted
همگنی	Scale-Location, BP/NCV	p > 0.05 and horizontal line: Accepted
استقلال	ACF/PACF, DW	DW ≈ 2 and p > 0.05: Accepted
غیرمعمول	Outliers/Leverage/Cook’s	No critical points or small change after removal: Accepted

چک‌لیست گام‌به‌گام تحلیل رگرسیون خطی چندگانه (برای دانشجویان)

در این بخش، چک‌لیست کامل و دقیق تحلیل رگرسیون خطی چندگانه ارائه شده است. لطفاً تمام مراحل را به ترتیب انجام دهید و در گزارش خود، برای هر مرحله کد R، خروجی مربوطه، نمودارها و تفسیر کامل را بیاورید.

۱. بررسی اولیه داده‌ها

داده را فراخوانی کنید و نام متغیرها را بررسی کنید.
با استفاده از str()، summary() و head() ساختار داده، نوع متغیرها، خلاصه آماری و چند سطر اول را مشاهده کنید.
بررسی کنید آیا مقادیر گم‌شده (NA) وجود دارد؟ اگر وجود داشت، گزارش دهید و در صورت نیاز (مثلاً با na.omit()) آن‌ها را حذف کنید.

۲. تحلیل اکتشافی اولیه (Exploratory Data Analysis)

هیستوگرام متغیر پاسخ (متغیر وابسته) را رسم کنید و توزیع آن را تفسیر کنید (آیا نرمال به نظر می‌رسد؟ آیا نیاز به تبدیل دارد؟).
برای هر متغیر مستقل پیوسته، نمودار پراکنش (scatter plot) در مقابل متغیر پاسخ رسم کنید (plot(y ~ x_i) یا با ggplot2).
ضریب همبستگی پیرسون بین متغیر پاسخ و هر متغیر مستقل را محاسبه کنید (cor.test(y, x_i)) و معناداری آن را گزارش دهید (p-value < 0.05 نشان‌دهنده رابطه معنادار است).
ماتریس همبستگی تمام متغیرهای عددی را رسم کنید (corrplot یا pairs()) و همبستگی‌های قوی (بیش از ۰.۷ یا کمتر از -۰.۷) را شناسایی کنید (نشانه احتمالی هم‌خطی).

۳. برازش مدل اولیه رگرسیون چندگانه

مدل اولیه را با تمام متغیرهای اختصاص داده شده برازش دهید:
fit <- lm(response ~ pred1 + pred2 + pred3 + ..., data = dat)
خروجی summary(fit) را گزارش دهید.
تفسیر کنید:
- ضرایب هر متغیر (شیب و عرض از مبدأ).
- معناداری هر ضریب (p-value < 0.05 → معنادار).
- R-squared و Adjusted R-squared (چه درصدی از تغییرات پاسخ توسط مدل توضیح داده می‌شود؟).
- آماره F کلی مدل (آیا مدل کلی معنادار است؟).

۴. بررسی هم‌خطی (Multicollinearity)

فاکتور تورم واریانس (VIF) را برای متغیرهای مستقل محاسبه کنید (vif(fit) از پکیج car).
اگر VIF برای هر متغیر بیشتر از ۵ (یا ۱۰ در برخی منابع) باشد، هم‌خطی وجود دارد. پیشنهاد دهید کدام متغیر حذف یا ترکیب شود.

۵. بررسی فرضیات مدل رگرسیون

۵.۱. خطی بودن رابطه

نمودار باقیمانده‌ها در مقابل مقادیر برازش‌شده را رسم کنید (plot(fit, which = 1) یا دستی).
اگر الگویی سیستماتیک (مثل منحنی یا U شکل) دیده شود، رابطه خطی نیست.

۵.۲. نرمال بودن باقیمانده‌ها

هیستوگرام باقیمانده‌ها را رسم کنید.
نمودار QQ باقیمانده‌ها را رسم کنید (plot(fit, which = 2) یا qqnorm(residuals(fit))).
آزمون Shapiro-Wilk را اجرا کنید (shapiro.test(residuals(fit))).
تفسیر: اگر نقاط در QQ-plot روی خط باشند و p-value > 0.05، نرمالیتی قابل قبول است.

۵.۳. ثابت بودن واریانس باقیمانده‌ها (همگنی واریانس)

نمودار Scale-Location را رسم کنید (plot(fit, which = 3)).
آزمون Breusch-Pagan و/یا NCV را اجرا کنید (bptest(fit) و ncvTest(fit)).
اگر خط در نمودار افقی باشد و p-value > 0.05، واریانس ثابت است.

۵.۴. استقلال باقیمانده‌ها

نمودار ACF باقیمانده‌ها را رسم کنید (acf(residuals(fit))).
آزمون Durbin-Watson را اجرا کنید (durbinWatsonTest(fit)).
اگر DW نزدیک به ۲ باشد و p-value > 0.05، استقلال تأیید می‌شود.

۵.۵. بررسی داده‌های غیرمعمول

داده‌های پرت: باقیمانده‌های استاندارد شده (rstandard(fit)) — مقادیر مطلق بیشتر از ۲ یا ۳ را گزارش دهید.
داده‌های با اهرم بالا: مقادیر hat (hatvalues(fit)) — بیشتر از ۲(k+1)/n را گزارش دهید (k = تعداد متغیرهای مستقل).
داده‌های مؤثر: فاصله کوک (cooks.distance(fit)) — مقادیر بیشتر از ۴/n را گزارش دهید.
نمودارهای تشخیصی (plot(fit, which = 4) و which = 5) را رسم و تفسیر کنید.

۶. اصلاح مدل (در صورت نیاز)

اگر مشکلی (مثل غیرنرمالیتی یا ناهمگنی) مشاهده شد، تبدیل مناسب اعمال کنید (مثل log، sqrt، یا Box-Cox).
مدل اصلاح‌شده را برازش دهید و دوباره فرضیات را بررسی کنید.
اگر داده مؤثر یا پرت قوی وجود داشت، مدل را بدون آن‌ها برازش دهید و نتایج (ضرایب، R-squared) را مقایسه کنید.

۷. مدل نهایی و تفسیر

مدل نهایی خود را (اصلاح‌شده یا اولیه) معرفی کنید.
ضرایب معنادار را تفسیر کنید (از نظر عملی و آماری).
قدرت توضیحی مدل (Adjusted R-squared) را گزارش دهید.
نتیجه‌گیری کنید که آیا مدل قابل اعتماد است یا خیر.

۸. پیش‌بینی

برای ۲-۳ مجموعه مقدار جدید از متغیرهای مستقل، پیش‌بینی کنید (با predict(fit, newdata, interval = "confidence") و "prediction").
بازه‌ها را تفسیر کنید.

این چک‌لیست را دقیقاً دنبال کنید و هر مرحله را با کد، خروجی و تفسیر در گزارش خود بگنجانید. گزارش باید کامل، منظم و حرفه‌ای باشد. موفق باشید!

رگرسیون خطی کاربردی و عیب یابی آن

حیدرعلی مردانی فرد

زمستان ۱۴۰۱ (به روزرسانی: آذر ۱۴۰۴)

مقدمه

۱. فراخوانی داده‌ها

۲. بررسی و پاک‌سازی داده‌ها

۲.۱. کلاس tibble

۲.۲. پکیج dplyr

۳. رگرسیون و برازش مدل

۳.۱. تبدیل متغیرها

۳.۲. برازش مدل رگرسیون ساده

۳.۳. مثال رگرسیون چندگانه

۳.۴. پیش‌بینی و بازه‌ها

۳.۵. مقایسه مدل‌ها

۴. ارزیابی مدل و بررسی باقیمانده‌ها

۴.۱. بررسی نرمال بودن باقیمانده‌ها

۴.۲. بررسی خطی بودن رابطه

۴.۳. بررسی ثابت بودن واریانس خطاها (همگنی)

۴.۴. بررسی استقلال باقیمانده‌ها

۴.۵. بررسی وجود داده‌های غیرمعمول

داده‌های پرت

داده‌های با اهرم بالا

داده‌های مؤثر

۵. چک‌لیست عیب‌یابی مدل

چک‌لیست گام‌به‌گام تحلیل رگرسیون خطی چندگانه (برای دانشجویان)

۱. بررسی اولیه داده‌ها

۲. تحلیل اکتشافی اولیه (Exploratory Data Analysis)

۳. برازش مدل اولیه رگرسیون چندگانه

۴. بررسی هم‌خطی (Multicollinearity)

۵. بررسی فرضیات مدل رگرسیون

۵.۱. خطی بودن رابطه

۵.۲. نرمال بودن باقیمانده‌ها

۵.۳. ثابت بودن واریانس باقیمانده‌ها (همگنی واریانس)

۵.۴. استقلال باقیمانده‌ها

۵.۵. بررسی داده‌های غیرمعمول

۶. اصلاح مدل (در صورت نیاز)

۷. مدل نهایی و تفسیر

۸. پیش‌بینی