پیاده سازی رگرسیون چندگانه | گام به گام و به زبان ساده

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

برای آشنایی بیشتر با مدل رگرسیون خطی، بهتر است مطالب دیگر مجله فرادرس با عنوان‌های رگرسیون خطی — مفهوم و محاسبات به زبان ساده و رگرسیون خطی چندگانه (Multiple Linear Regression) — به زبان ساده را مطالعه کنید. همچنین خواندن مطالب هم خطی در مدل رگرسیونی — به زبان ساده و انواع روش های رگرسیونی — راهنمای جامع نیز خالی از لطف نیست.

پیاده سازی رگرسیون چندگانه

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

فیلم آموزش درس رگرسیون 1 – رگرسیون خطی در فرادرس

کلیک کنید

مقدمات و تعریف مدل رگرسیون چندگانه

فرض کنید $$X$$ یک ماتریس با $$n$$ سطر و $$k$$ ستون است. در این صورت آن را به شکل زیر در نظر می‌گیریم. این مقادیر همان مشاهدات برای هر یک از ابعاد یا متغیرها محسوب می‌شود. بنابراین ماتریس $$X$$ را ماتریس مشاهدات در نظر می‌گیریم.

$$\large {\displaystyle X = \left(\begin{array}{c}X'_1\\X'_2\\ \vdots \\ X'_n \end{array}\right) }$$

واضح است که $$X_i$$ یک بردار ستونی به صورت $$k$$ سطر و یک ستون است. به همین منظور $$X_i'$$، ترانهاده آن محسوب می‌شود که رابطه بالا را به یک ماتریس تبدیل می‌کند. توجه داشته باشید که اندیس $$i$$‌ بیانگر مشاهده $$i$$ام است.

نکته: همانطور که گفته شد، $$X$$ یک ماتریس $$n \times k$$ است که $$n$$ تعداد مشاهدات و $$k$$ تعداد متغیرهای مستقل است.

متغیر وابسته یا $$Y$$ را هم به صورت یک بردار با $$n$$ سطر در نظر بگیرید. در اینصورت $$y = (y_1, \ldots, y_n)'$$ یک بردار ستونی خواهد بود.

طبق مدل OLS یا همان «کمترین مربعات عادی» (Ordinary Least Square)، برآورد حاصل کمینه‌سازی خطا به صورت زیر خواهد بود.

$$\large {\displaystyle \widehat{\beta} = (X'X)^{-1}X'y }$$

از طرفی ماتریس کوواریانس-واریانس برآوردگرهای OLS نیز به شکل زیر خواهد بود.

$$\large {\displaystyle \widehat{\sigma} (X'X)^{-1} }$$

که در آن $$\widehat{\sigma}$$ به وسیله رابطه زیر برآورد می‌شود.

$$\large {\displaystyle \widehat{ \sigma} = \dfrac{ \sum_{i = 1}^n \widehat{u}_i}{n - k}}$$

برآورد خطا ($$\widehat{u}_i$$) نیز بر طبق رابطه زیر حاصل می‌شود.

$$\large {\displaystyle \widehat{u}_i = y_i - X_i' \widehat{\beta} }$$

نکته: مدل مطرح شده در قسمت بالا، با این شرط نوشته شده است که مدل رگرسیونی دارای عرض از مبدا (Intercept) نیست. اگر مدل به همراه عرض از مبدا (مقدار ثابت به عنوان میانگین) در مدل لحاظ شود، باید ماتریس $$X$$ را به صورت زیر در نظر بگیریم. به این ترتیب خواهیم داشت:

$$\large {\displaystyle X = \left(\begin{array}{c}1 & X'_1\\ 1 &X'_2\\ \vdots \\1 & X'_n \end{array}\right) }$$

به یاد داشته باشید که هنگام پیاده سازی رگرسیون چندگانه با کدهای R، به این موضوع توجه خواهیم داشت.

برای مثال مدل رگرسیون زیر را با دو متغیر مستقل و عبارت عرض از مبدا در نظر بگیرید.

$$\large {\displaystyle y_i = \beta_0 + \beta_1 x_{1,i} + \beta_2 x_{2,i} + u_{i}\;\;, \;\;\;\;\; i = 1 , 2 , \ldots ,n }$$

در این صورت مدل رگرسیونی را به صورت ماتریس و برداری به شکل زیر خواهیم نوشت:

$$\large {\displaystyle y  = X \beta + u }$$

در رابطه بالا، بردارها و ماتریس‌ها به شکل زیر معرفی می‌شوند. ابتدا متغیرهای مستقل و وابسته را معرفی می‌کنیم.

$$\large {\displaystyle X = \left(\begin{array}{c} 1 & x_{1,1} & x_{1,2} \\ 1 & x_{2,1} & x_{2,2} \\ \vdots& \vdots & \vdots \\1 & x_{n,1} & x_{n,2} \end{array} \right) }$$

و

$$\large {\displaystyle y = \left(\begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array}\right) }$$

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

$$\large {\displaystyle \beta = \left(\begin{array}{c}\beta_0 \\ \beta_1\\ \beta_2 \end{array}\right) }$$

و

$$\large {\displaystyle u = \left( \begin{array}{c} u_1 \\ \vdots \\ u_n \end{array} \right) }$$

به این ترتیب با اجرای محاسباتی که برای برآورد پارامترها مدل رگرسیونی بیان شد، می‌توان مدل خطی را ساخت و خطای مدل را هم محاسبه کرد.

کدهای پیاده‌سازی محاسبات رگرسیون چندگانه

پیاده‌سازی رگرسیون چندگانه (تک متغیره و چندگانه) در زبان برنامه‌نویسی R به استفاده از تابع lm و معرفی پارامترهای مربوطه قابل اجرا است. کافی است که متغیر وابسته و مستقل را از یک چارچوب داده مشخص کرده و به کمک عملگر  معرفی مدل (~)، نوع رابطه بین متغیر مستقل و وابسته را مشخص کنید. برای آشنایی بیشتر با این شیوه اجرای رگرسیون خطی، بهتر است به نوشتار دیگری از مجله فرادرس با عنوان روش‌های رگرسیون در R — کاربرد در یادگیری ماشین (قسمت اول) مراجعه فرمایید.

فیلم آموزش درس رگرسیون 2 – جامع و کاربردی در فرادرس

کلیک کنید

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

1library(ggplot2)
2####################################
3#     Regression Function          #
4####################################
5
6linear_regression <- function(data, dep, indep, intercept = TRUE) {
7y <- as.matrix(data[, dep])
8x <- as.matrix(data[, indep])
9if (intercept == TRUE) {
10x <- cbind(1, x)
11}
12beta <- c(solve(crossprod(x)) %*% crossprod(x, y))
13fits <- x %*% beta
14resids <- y - fits
15sigma2 <- sum(resids^2)/(length(resids)-ncol(x))
16se <- sqrt(diag(sigma2 * solve(crossprod(x))))
17names(beta) <- colnames(x)
18return_obj <- list(beta = beta, se = se,
19                   residuals = c(resids), fitted = c(fits),
20                   sigma = sqrt(sigma2), dep = dep, indep = indep,
21                   intercept = intercept, y = c(y))
22
23class(return_obj) <- "linear_regression"
24}
25

فیلم آموزش برنامه نویسی R و نرم افزار RStudio – مقدماتی در فرادرس

کلیک کنید

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

پارامتر اول (data)، یک چارچوب داده (Dataframe) است که ستون‌های مربوط به متغیر وابسته (dep) و مستقل (indep) نیز با ذکر اسامی آن‌ها در چارچوب data مشخص می‌شود. برای آشنایی با نحوه ورود پارامترها به مثالی که در ادامه مشاهده خواهید کرد، توجه کنید. در انتها نیز پارامتر منطقی intercept، با مقدار پیش‌فرض TRUE، نوع مدل رگرسیونی را با توجه به حضور یا عدم عرض از مبداء، تعیین می‌کند.

در خطوط بعدی، با استفاده از تابع as.matrix متغیرهای مستقل به یک ماتریس، مطابق با فرمول‌های بالا تبدیل می‌شوند. همچنین به منظور امکان اجرای محاسبات ماتریسی و برداری، متغیر وابسته نیز به شکل یک ماتریس در می‌آید. توجه داشته باشید که با حضور عرض از مبدا، به ماتریس متغیرهای مستقل، یک ستون با مقادیر ۱ اضافه می‌شود.

در قسمت بعدی نیز با توجه به فرمول گفته شده، برآورد پارامترها در متغیر beta قرار داده می‌شود. این متغیر، از نوع بردار است که هر یک از مقادیر $$\beta_0$$ تا $$\beta_2$$ را در خود جای خواهد داد.

نکته: تابع crossprod همان محاسبه ضرب ترانهاده یک ماتریس در ماتریس دیگر را نشان می‌دهد. بنابراین (crossprod(x,y همان محاسبه $$X'y$$ است. در صورتی که پارامتر دوم در این تابع مشخص نشده باشد، حاصل ضرب ترانهاده یک ماتریس در خودش حاصل می‌شود، در نتیجه (crossprod(x نیز مربوط به بدست آوردن $$X'X$$‌ است. برای بدست آوردن معکوس این ماتریس نیز از تابع solve استفاده شده است.

متغیر fit، شامل مقادیر برازش شده بوده که نتایج را به صورت یک بردار ذخیره می‌کند. عملگر %*%، همان ضرب ماتریسی معمولی است. باقی‌مانده‌های مدل رگرسیونی نیز در resids و براساس تفاضل مقدار واقعی و برازش شده توسط مدل، حاصل شده است.

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

در انتهای کد نیز خروجی به کمک دستور return طبق یک لیست (list) به صورت مقادیر پارامترها (beta)،  خطای استاندارد (se)، بردار سطری باقی‌مانده‌ها (residuals) بردار سطری مقادیر برازش شده (fitted)، انحراف معیار متغیر وابسته (sigma) که توسط جذر گیری از واریانس بدست آمده و همچنین نام متغیرهای مستقل (indp) و وابسته (dep) به همراه وضعیت حضور عرض از مبدا (intercept) و مقادیر بردار سطری متغیر وابسته (y)، مشخص شده است.

تا اینجا انجام محاسبات مربوط به پیاده سازی رگرسیون خطی چندگانه به پایان رسیده است. ولی می‌خواهیم با ایجاد دو کلاس برای چاپ خروجی کار را به پایان ببریم.

نوشتن دو روش برای کلاس رگرسیون خطی به گونه‌ای است که اگر obj کلاس باشد رگرسیون خطی خروجی‌ها را به صورتی که در ادامه می‌خوانید ارائه دهد.

1. چاپ (obj) چاپ جدول بر روی کنسول یک جدول با ردیف 1) نام از متغیرها، 2) ضرایب برآورد شده و 3) خطای استاندارد برآوردگرها.
2. رسم نمودار (obj) به صورت یک نمودار پراکندگی از مشاهدات مرتب شده که در محور افقی شماره یا اندیس مشاهده و در محور عمودی نیز مقدار واقعی و برآورد شده برای متغیر وابسته نمایش داده شود. این کار را به کمک کتابخانه ggplot2 اجرا خواهیم کرد.

ابتدا با یک مثال برای تابع linear_regression پرداخته، سپس کلاس‌های print و plot را تعریف می‌کنیم. به این ترتیب پیاده سازی رگرسیون چندگانه را به شکل دلخواه انجام داده‌ایم.

فیلم آموزش همبستگی و رگرسیون خطی در اس پی اس اس SPSS در فرادرس

کلیک کنید

در کتابخانه ggplot2، یک چارچوب داده به نام economics وجود دارد که از آن برای اجرای مدل رگرسیون خطی استفاده می‌کنیم. بهتر است از دستور زیر به منظور نمایش محتویات چند سطر اول این فایل اطلاعاتی کمک بگیرم.

1> head(economics)
2# A tibble: 6 x 6
3  date         pce    pop psavert uempmed unemploy
4  <date>     <dbl>  <dbl>   <dbl>   <dbl>    <dbl>
51 1967-07-01  507. 198712    12.6     4.5     2944
62 1967-08-01  510. 198911    12.6     4.7     2945
73 1967-09-01  516. 199113    11.9     4.6     2958
84 1967-10-01  512. 199311    12.9     4.9     3143
95 1967-11-01  517. 199498    12.8     4.7     3066
106 1967-12-01  525. 199657    11.8     4.8     3018

مشخص است که این مجموعه، شامل ۶ متغیر بوده که میزان رشد اقتصادی را در سال‌های مختلف بر حسب پارامترهای pce، pop، psavert  uempmed و unemploy مشخص می‌کند. می‌خواهیم مقدار pce را به کمک بقیه متغیرها، طبق یک مدل خطی، برآورد و پیش‌بینی کنیم.

دستورات زیر به این منظور نوشته شده‌اند. به وضعیت معرفی متغیرهای مستقل (indep) و وابسته (dep) در تابع linear_regression توجه کنید.

1####################################
2#     Regression Example           #
3####################################
4
5obj <- linear_regression(ggplot2::economics,
6                         dep = "pce",
7                         indep = c("unemploy", "pop", "psavert"))
8
9
10class(obj)
11
12###################################
13#     Customization Output        #
14#   Method Print and Plot         #
15###################################
16
17
18print <- function(object) {
19  cat("  \t Intercept ",object$indep,"\nCoeffcient  ", object$beta, "\nSe \t   ", object$se, "\n")
20
21}
22
23plot <- function(object) {
24  ggplot(as.data.frame(cbind(object$fitted,object$y)),
25         aes(x=(1:length(object$fitted)))) +
26    geom_point(aes(y = object$fitted),color="blue")+
27    geom_point(aes(y = object$y),color="red")+
28    labs(x= "Indexs", y="Values")
29}
30

کلاسی که برای  obj معرفی شده، همان linear_regression است که توسط فرمان (class(obj مشخص شده است.

برای نمایش خروجی به شکلی که فقط مقادیر مورد نظر از نتایج پیاده سازی رگرسیون چندگانه ظاهر شود، کلاس  print را تغییر داده و معرفی کرده‌ایم. همانطور که مشخص است، با بهره‌گیری از تابع cat، نتایج به به دلخواه توسط فرمان print ظاهر کرده‌ایم.

در قسمت انتهایی به کمک کلاسی که از تابع plot ایجاد کرده‌ایم، نمودار را به کمک کتابخانه ggplot2‌ و با دستور ggplot و ملحقات آن (مانند معرفی محورها، نحوه ترسیم نقاط و برچسب‌ها) انجام داده‌ایم. به این ترتیب خروجی دستور plot برای obj به گونه‌ای متفاوت از دستور plot عادی عمل خواهد کرد و عملیات پیاده سازی رگرسیون چندگانه در زبان R به پایان می‌رسد.

کافی است ابتدا تابع را اجرا کرده، سپس، کلاس‌های تعریف شده را با کلید‌های میانبر ctrl+R در حافظه قرار دهید. به این ترتیب با فراخوانی قسمت کد مربوط به مثال، نتایج زیر حاصل خواهد شد.

نکته: با اجرای هر یک از دستورات print و plot، مطابق با کد مورد نظر، خروجی‌ها به صورت زیر ساخته می‌شوند. واضح است که این خروجی‌ها، مطابق با درخواست شما در کد برنامه، تعیین شده‌‌اند.

1####################################
2#     Getting Output New format    #
3####################################
4
5
6print(obj)
7plot(obj)

همانطور که در ادامه مشاهده می‌کنید، شکل خروجی توسط فرمان print تغییر یافته و فقط «پارامترها برآورد شده» (Coefficient) و «خطای برآورد» (SE) نمایش داده شده است.

به این ترتیب خروجی حاصل از اجرای دستور (plot(obj دیگر به شیوه عادی در R اجرا نشده و مطابق با کلاسی که معرفی کرده‌اید، نمودار ترسیم و به کمک کتابخانه ggplot2 رسم می‌شود.

خلاصه و نتیجه‌گیری

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