انواع روش های رگرسیونی — راهنمای جامع

در دیگر نوشتارهای مجله فرادرس با مفهوم رگرسیون خطی و همچنین نحوه محاسبات آن آشنا شده‌‌اید. در این نوشتار با انواع روش های رگرسیونی خطی (Linear) و غیرخطی (NonLinear) آشنا خواهیم شد و کدهایی مربوط به پیاده‌سازی آن‌ها را در زبان برنامه‌نویسی R، فرا می‌گیریم. در این بین ابتدا با مفاهیم اولیه رگرسیون و سپس دستورات و کدهای R که تکنیک‌های مختلف رگرسیونی را اجرا می‌کنند، آشنا می‌شویم. از آنجایی که تکنیک‌های رگرسیونی برای داده‌های کیفی و کمی، به شکلی جداگانه‌ و متفاوت اجرا می‌شوند، اطلاع از نحوه اجرا این تکنیک‌های رگرسیونی در این حالت‌ها بسیار ضروری است، بخصوص برای کسانی که می‌خواهند در حوزه «علم داده» (Data Science) فعالیت کنند.

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

انواع روش های رگرسیونی

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

1. رگرسیون خطی (Linear Regression)
2. رگرسیون لجستیک (Logistic Regression)
3. رگرسیون چندکی (Quantile Regression)
4. رگرسیون ستیغی (Ridge Regression)
5. رگرسیون لاسو (Lasso Regression)
6. رگرسیون شبکه الاستیک (Elastic Net Regression)
7. رگرسیون مولفه‌های اصلی (Principle Component Regression)
8. رگرسیون کمترین مربعات جزئی (Partial Least Square (PLS) Regression)
9. رگرسیون بردار پشتیبان (Support Vector Regression)
10. رگرسیون ترتیبی (Ordinal Regression)
11. رگرسیون پواسون (Poisson Regression)
12. رگرسیون دوجمله‌ای منفی (Negative Binomial Regression)
13. رگرسیون شبه پواسن (Quasi Poisson Regression)
14. رگرسیون کاکس (Cox Regression)

آگاهی از این روش‌ها به یک «دانشمند داده» (Data Scientist) کمک می‌کند که بهترین روش و الگو را برای تحلیل داده‌های خود به کار ببرد و در نتیجه مدل‌های ساخته شده از بیشترین کارایی و دقت برخوردار شوند.

هر یک از روش‌های رگرسیونی، پیش‌فرض‌های مخصوص خود را دارد که برحسب ویژگی و مشخصات «متغیرهای توصیفی» (Explanatory Variables) و «متغیر پاسخ» (Response Variable) تعیین می‌شوند. توجه داشته باشید که گاهی به متغیرهای توصیفی، «متغیرهای مستقل» (Independent Variable) و به متغیر پاسخ، «متغیر وابسته» (Dependent Variable) می‌گویند. البته اصطلاح متغیرهای «رگرسور» (Regressor) یا «پیشگو» (Predictor) نیز گاهی برای متغیرهای توصیفی به کار می‌رود.

فرضیه‌های صحت مدل رگرسیونی

قبل از اجرای مدل رگرسیونی باید بوسیله رسم نمودار یا محاسبه ضریب همبستگی و نظایر آن، وجود وابستگی بین متغیرها توصیفی با متغیر وابسته را مورد بررسی قرار دهیم و در صورتی که نتایج حاصل از این ابزارها، نشانگر وجود رابطه (رابطه خطی یا رابطه غیر خطی) بین این دو دسته متغیر باشد، به دنبال ایجاد مدل رگرسیونی باشیم.

ولی در اینجا هم کار تمام نشده است. پس از محاسبات و برآورد پارامترهای مدل باید فرضیه‌هایی که مدل رگرسیونی براساس آن استوار شده است نیز آزمون شود. بیشتر این گونه آزمون‌ها براساس باقی‌مانده‌های مدل صورت می‌گیرد در نتیجه باید آن‌ها را پس از پیدا کردن مدل و محاسبه باقی‌مانده‌ها (Residuals) مورد سنجش قرار داد.

این آزمون‌ها در فهرست زیر معرفی و توضیح داده شده‌اند. در ادامه با کدهایی از زبان R آشنا می‌شویم که این گونه آزمون‌ها را برایمان در مدل رگرسیونی اجرا می‌کنند.

• بررسی نقاط پرت (Outlier): از آنجایی که وجود نقاط دور افتاده یا پرت (Outlier)، باعث می‌شود که برآورد پارامترهای مدل رگرسیونی به درستی صورت نگیرد، قبل از اجرای محاسبات مربوط به مدل رگرسیونی باید از عدم چنین داده‌های اطمینان حاصل پیدا کنیم و اگر به این گونه مشاهدات برخوردیم، آن‌ها را از مدل خارج کرده، سپس دوباره مدل‌سازی را انجام دهیم. نقاط پرت در اینجا به مشاهداتی اشاره دارد که باقی‌مانده‌های خیلی بزرگی نسبت به بقیه نقاط دارند.
• هم‌خطی و هم‌خطی چندگانه: اگر متغیرهای توصیفی نسبت به هم وابستگی داشته باشند، پارامترهای مدل، واریانس بزرگی پیدا کرده، در نتیجه قابل اعتماد نیستند. در این حالت باید مدل رگرسیونی را از جهت وجود مشکل هم‌خطی (Colinearity)  یا  هم‌خطی چندگانه (Multicolinearity) مورد بررسی قرار داد.
• ناهم‌واریانسی (Heteroscedasticity): باقی‌مانده‌های حاصل از مدل رگرسیونی باید دارای واریانس ثابت باشند. معمولا این موضوع به متغیر وابسته نیز بر می‌گردد. اگر با تغییر مقادیر متغیرهای توصیفی، واریانس متغیر وابسته تغییر کند، با مشکل ناهم‌واریانس مواجه هستیم. وزن‌دهی به مشاهدات یا تبدیل روی متغیر وابسته، روش‌هایی برای ثابت کردن واریانس متغیر پاسخ محسوب می‌شوند.
• بیش‌برازش و کم‌برازش (Under and Over Fitting): به کارگیری متغیرهای توصیفی زیاد در کاهش میزان خطای مدل رگرسیونی موثر است ولی در این صورت ممکن است مدل دچار مشکل «بیش‌برازش» (Overfitting) شود. به این معنی که مدل قادر به پیش‌بینی دقیق مقادیر متغیر پاسخ برای مشاهداتی است براساس آن مدل ساخته شده ولی برای داده‌های جدید، توانایی پیش‌بینی مناسب را ندارد. عکس این حالت که کم‌برازش (Underfitting) نامیده می‌شود، ناکافی بودن متغیرهای توصیفی در مدل است، بطوری که واریانس مقادیر پیش‌بینی بسیار زیاد خواهد بود. به این ترتیب باید از مدل رگرسیونی استفاده کرد که تعداد متغیرهای آن مناسب بوده و چنین مسائلی را بوجود نیاورد.

۱. رگرسیون خطی

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

زمانی که فقط یک متغیر مستقل وجود داشته باشد، مدل رگرسیونی خطی را ساده (Simple Regression) می‌نامند و اگر بیش از یک متغیر مستقل (توصیفی) وجود داشته باشد، رگرسیون را چندگانه (Multiple Regression) می‌گویند.

در اینجا، پارامترهای مدل، بوسیله کمینه‌سازی «مجموع مربعات خطا» (Sum of Square Error) صورت می‌گیرد که به این تکنیک، «رگرسیون عادی» (Ordinary Least Square) یا OLS نیز گفته می‌شود.

برای برازش و اجرای رگرسیون خطی، از مجموعه داده swiss که در بسته datasets از زبان R قرار دارد، استفاده می‌کنیم. ابتدا به این مجموعه داده نگاهی می‌اندازیم. برای دسترسی به این مجموعه کافی است کتابخانه datasets را بارگذاری کنید. البته از اینجا نیز این مجموعه داده با قالب فشرده را می‌توانید دریافت کنید. توجه داشته باشید که پس از خارج کردن این فایل از حالت فشرده، به یک فایل csv برخورد خواهید کرد که با دستور read.csv در R قابل بارگذاری است.

مجموعه داده swiss شامل 47 مشاهده و شش متغیر است که مربوط به ویژگی‌های اقتصادی و رشد جمعیت در 47 ناحیه مختلف فرانسوی زبان کشور سوئیس هستند. متغیرهای مربوط به باروری (Fertility)، درصد جمعیت کشاورزان از کل (Agriculture)، درصد نفرات برتر در آزمون نظامی (Examination)،  درصد افراد با تحصیلات عالی (Education)، درصد کاتولیک‌ها (Catholic) و درصد فوت نوزادان کمتر از یک سال (Infant.Mortality) در سال 1888 در این فایل، ثبت شده‌اند.

1library(datasets)
2head(swiss)
3             Fertility Agriculture Examination Education Catholic Infant.Mortality
4Courtelary        80.2        17.0          15        12     9.96             22.2
5Delemont          83.1        45.1           6         9    84.84             22.2
6Franches-Mnt      92.5        39.7           5         5    93.40             20.2
7Moutier           85.8        36.5          12         7    33.77             20.3
8Neuveville        76.9        43.5          17        15     5.16             20.6
9Porrentruy        76.1        35.3           9         7    90.57             26.6

به منظور انجام محاسبات و برآورد پارامترهای رگرسیون OLS از تابع lm در زبان‌ برنامه‌نویسی R استفاده می‌شود. در ادامه نمونه‌ای از کدهای مربوطه را مشاهده می‌کنید.

1model = lm(Fertility ~ .,data = swiss)
2lm_coeff = model$coefficients
3lm_coeff
4summary(model)

همانطور که مشخص است رابطه بین متغیر وابسته (Fertility) با بقیه متغیرها که نقش متغیرهای مستقل را ایفا می‌کنند به صورت زیر نوشته شده است. واضح است که مجموعه داده نیز swiss نام دارد. توجه داشته باشید که در اینجا بقیه متغیرها را با نماد «.» مشخص کرده‌ایم.

$$\large lm (Fertility \sim . , data = swiss)$$

حاصل اجرای این کد به صورت زیر است:

1   (Intercept)      Agriculture      Examination        Education         Catholic 
2      66.9151817       -0.1721140       -0.2580082       -0.8709401        0.1041153 
3Infant.Mortality 
4       1.0770481 
5
6lm(formula = Fertility ~ ., data = swiss)
7
8Residuals:
9     Min       1Q   Median       3Q      Max 
10-15.2743  -5.2617   0.5032   4.1198  15.3213 
11
12Coefficients:
13                 Estimate Std. Error t value Pr(>|t|)    
14(Intercept)      66.91518   10.70604   6.250 1.91e-07 ***
15Agriculture      -0.17211    0.07030  -2.448  0.01873 *  
16Examination      -0.25801    0.25388  -1.016  0.31546    
17Education        -0.87094    0.18303  -4.758 2.43e-05 ***
18Catholic          0.10412    0.03526   2.953  0.00519 ** 
19Infant.Mortality  1.07705    0.38172   2.822  0.00734 ** 
20---
21Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
22
23Residual standard error: 7.165 on 41 degrees of freedom
24Multiple R-squared:  0.7067, Adjusted R-squared:  0.671 
25F-statistic: 19.76 on 5 and 41 DF,  p-value: 5.594e-10

براساس خروجی‌ها و مقدار R-squared که برابر با 0٫7067 است،‌ می‌توان گفت که حدود ۷۰ درصد تغییرات متغیر وابسته، برحسب مدل بیان شده. همچنین در ستون (|Pr(>|t که همان «مقدار احتمال» (p-Value) است می‌توان همه پارامترها به جز پارامتر Examination را معنی‌دار (Significant) دانست. معنی‌دار بودن از لحاظ آماری در اینجا به معنی مخالف صفر بودن پارامتر است. به این ترتیب می‌توان حضور آن را در مدل، معنی‌دار در نظر گرفت.

نکته: متغیر یا پارامترهایی که از لحاظ آماری، معنی‌دار (مخالف صفر) هستند با علامت * در خروجی تابع lm مشخص شده‌اند.

۲. رگرسیون لجستیک

در رگرسیون لجستیک (Logistic Regression)، متغیر وابسته، به صورت دو دویی (Binary) است. به این معنی که مقادیر آن به دو طبقه صفر و یک دسته‌بندی شده‌اند. البته زمانی که از رگرسیون چند جمله‌ای لجستیک (Multinomial Logistic Regression) استفاده می‌کنید، ممکن است تعداد سطوح متغیر طبقه‌ای بیشتر از دو باشد. در این حال مدل رگرسیون لجستیک به شکل زیر نوشته می‌شود.

$$\large P(Y = 1) = \dfrac{1}{1 + e^{ -(\beta_0 + \beta_1X_1 + \beta_2X_2 + \ldots + \beta_kX_k)}}$$

واضح است که در این مدل رگرسیونی، خطاها، دارای توزیع نرمال نیستند و متغیر وابسته دارای توزیع دو یا چند جمله‌ای است در نتیجه نمی‌توان از مدل رگرسیون ساده یا خطی استفاده کرد.

نکته: معمولا از این شیوه یا مدل رگرسیونی، برای طبقه‌بندی کردن مشاهدات جدید برحسب مقادیر قبلی استفاده می‌کنند و به نوع «یادگیری نظارت شده» (Supervised Learning) محسوب می‌شود. به این ترتیب اگر مقدار $$P(Y=1)$$ از یک مقدار آستانه (مثلا 0٫5) بیشتر باشد، آن مشاهده را در گروه ۱ طبقه‌بندی می‌کنیم.

برای اجرای رگرسیون لجستیک، از مجموعه داده‌های سرطان سینه (case2002) که در بسته یا کتابخانه Sleuth2 قرار دارد، استفاده می‌کنیم. این مجموعه داده شامل 147 مشاهده و 7 متغیر است. نمونه‌ای از اطلاعات مربوط به این داده‌ها را در ادامه مشاهده می‌کنید. اگر دسترسی به این کتابخانه برایتان مقدور نیست از اینجا نیز می‌توانید این فایل را دریافت و از حالت فشرده خارج کنید. یک فایل به نام case2002 از نوع csv در اختیارتان قرار خواهد گرفت که با دستور read.csv در R قابل بارگذاری است.

1> head(case2002)
2          LC   FM   SS     BK AG YR CD
31 LungCancer Male  Low   Bird 37 19 12
42 LungCancer Male  Low   Bird 41 22 15
53 LungCancer Male High NoBird 43 19 15
64 LungCancer Male  Low   Bird 46 24 15
75 LungCancer Male  Low   Bird 49 31 20
86 LungCancer Male High NoBird 51 24 15

قرار است براساس متغیر مصرف سیگار روزانه (CD)، تشخیص دهیم که فرد با چه احتمالی دچار سرطان سینه (LC) خواهد شد.

1model <- glm(LC~CD,data = case2002, family = "binomial")
2summary(model)

نتیجه اجرای این کد، محاسبه احتمال داشتن سرطان براساس تعداد سیگارهای کشیده شده در روز است.

1glm(formula = LC ~ CD, family = "binomial", data = case2002)
2
3Deviance Residuals: 
4    Min       1Q   Median       3Q      Max  
5-1.5148  -0.9688  -0.7166   1.3449   1.8603  
6
7Coefficients:
8            Estimate Std. Error z value Pr(>|z|)    
9(Intercept) -1.53541    0.37707  -4.072 4.66e-05 ***
10CD           0.05113    0.01939   2.637  0.00836 ** 
11---
12Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
13
14(Dispersion parameter for binomial family taken to be 1)
15
16    Null deviance: 187.14  on 146  degrees of freedom
17Residual deviance: 179.62  on 145  degrees of freedom
18AIC: 183.62
19
20Number of Fisher Scoring iterations: 4

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

1        1         2         3         4         5         6         7         8         9        10 
20.2845787 0.3168095 0.3168095 0.3168095 0.3745313 0.3168095 0.3745313 0.3745313 0.2642256 0.4360573 
3       11        12        13        14        15        16        17        18        19        20 
40.6247508 0.4360573 0.3745313 0.4360573 0.3745313 0.4360573 0.4996171 0.3745313 0.3168095 0.3168095 
5       21        22        23        24        25        26        27        28        29        30 
60.3168095 0.3745313 0.3745313 0.3168095 0.6247508 0.3745313 0.2642256 0.4996171 0.3745313 0.3168095 
7       31        32        33        34        35        36        37        38        39        40 
80.3279782 0.2006837 0.3745313 0.2642256 0.4996171 0.4360573 0.2642256 0.3168095 0.3168095 0.4360573 
9       41        42        43        44        45        46        47        48        49        50 
100.3745313 0.3745313 0.3745313 0.2264199 0.3745313 0.3168095 0.1772028 0.3745313 0.4360573 0.1926074 
11       51        52        53        54        55        56        57        58        59        60 
120.3745313 0.4360573 0.2448299 0.3168095 0.4360573 0.3745313 0.6247508 0.2642256 0.2642256 0.2090110

از آنجایی که تعداد مشاهدات زیاد است، فقط برای ۶۰ نفر اول، نتایج را نشان داده‌ایم. حال داشتن سرطان را برای افراد سیگاری برحسب مقدار احتمال ۵۰٪ انجام می‌دهیم.

1data$prediction <- model$fitted.values>0.5
2data$prediction

به این ترتیب خروجی به صورت زیر در آمده و نشان می‌دهد که چه کسانی به احتمال زیاد مبتلا به سرطان خواهند شد. این امر با مقدار TRUE برای هر فرد مشخص شده است.

1    1     2     3     4     5     6     7     8     9    10    11    12    13    14    15    16 
2FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE 
3   17    18    19    20    21    22    23    24    25    26    27    28    29    30    31    32 
4FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
5   33    34    35    36    37    38    39    40    41    42    43    44    45    46    47    48 
6FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 
7   49    50    51    52    53    54    55    56    57    58    59    60    61    62    63    64 
8FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE

۳. رگرسیون چندکی

رگرسیون چندکی (Quantile Regression) را می‌توان جایگزینی برای مدل رگرسیون خطی در نظر گرفت که نسبت به مشکلات حاصل از «نقاط دورافتاده یا پرت» (Outlier)، «چولگی زیاد» (High Skewness) و همچنین «ناهم‌واریانسی» (Heteroscedasticity) مقاوم است.

در رگرسیون خطی، میانگین متغیر وابسته به شرط مشاهدات برآورد می‌شود. در حقیقیت مدل رگرسیون خطی با در نظر گرفتن ماتریس $$X$$ و بردارهای $$\beta$$ و $$\alpha$$ به صورت زیر نوشته می‌شود.

$$\large E(Y|x) = \alpha + \beta X$$

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

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

استفاده از چندک‌ها (Quantile) که نسبت به شرایط گفته شده، مقاوم‌تر هستند، می‌تواند مدل رگرسیون کامل‌تر و دقیق‌تری را ارائه کند. در رگرسیون چندکی به جای برآورد میانگین متغیر وابسته، از برآورد چندک‌های آن (مانند صدک، دهک یا چارک) به شرط متغیر مستقل، کمک گرفته می‌شود.

نکته: در رگرسیون چندکی، باید متغیر وابسته از نوع عددی و از نوع مقیاس (Scale) با مقادیر پیوسته باشد تا امکان محاسبه چندک‌ها وجود داشته باشد.

اگر $$Q_{Y|X}(p)$$ را چندک $$p$$ام متغیر وابسته $$Y$$ به شرط $$X$$ در نظر بگیریم، آنگاه مدل رگرسیون به صورت زیر خواهد بود.

$$\large Q_{Y|X}(p) = X \beta_{p}$$

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

نکته: توجه داشته باشید که کمینه‌سازی خطا طبق رابطه زیر صورت خواهد گرفت.

$$\large p ( \sum |e_i|) + (1-p) \sum |e_i|$$

به این ترتیب اگر مقدار $$p=0.5$$ باشد، چندک به میانه تبدیل شده و مدل رگرسیونی را «رگرسیون میانه» (Median Regression) می‌نامند.

فرض کنید معادله رگرسیونی برای صدک ۲۵ام به صورت زیر بدست آمده باشد:

$$\large Q_{y|x}(0.25) = 5 + 7x$$

در این صورت هر واحد اضافه شدن مقدار متغیر مستقل (x) باعث افزایش ۷ واحدی مقدار چندک ۲۵ام خواهد شد.

برای استفاده از رگرسیون چندکی در زبان برنامه‌نویسی R از بسته یا کتابخانه quantreg استفاده خواهیم کرد. در کد زیر نحوه نصب و بارگذاری این کتابخانه مشخص شده است.

1install.packages("quantreg")
2library(quantreg)

به کمک تابع rq امکان اجرای رگرسیون چندکی فراهم می‌شود. کدی که در زیر مشاهده می‌کنید به اجرای رگرسیون چندکی روی مجموعه داده swiss پرداخته است. در اینجا Fertility متغیر وابسته و بقیه متغیرها، مستقل محسوب شده‌اند.

1model1 = rq(Fertility~.,data = swiss,tau = 0.25)
2summary(model1)

مقدار $$p$$‌ در تابع rq به کمک پارامتر tau مشخص می‌شود. واضح است که در مدل اجرا شده، از صدک ۲۵ام (چارک اول) استفاد شده است. خروجی به صورت زیر خواهد بود.

1tau: [1] 0.25
2Coefficients:
3                 coefficients lower bd upper bd
4(Intercept)      76.63132      2.12518 93.99111
5Agriculture      -0.18242     -0.44407  0.10603
6Examination      -0.53411     -0.91580  0.63449
7Education        -0.82689     -1.25865 -0.50734
8Catholic          0.06116      0.00420  0.22848
9Infant.Mortality  0.69341     -0.10562  2.36095

اگر مقدار tau را برابر با ۰٫۵ در نظر بگیریم، رگرسیون میانه را تولید خواهیم کرد.

1model2 = rq(Fertility~.,data = swiss,tau = 0.5)
2summary(model2)

خروجی به صورت زیر محاسبه خواهد شد.

1tau: [1] 0.5
2Coefficients:
3                 coefficients lower bd upper bd
4(Intercept)      63.49087     38.04597 87.66320
5Agriculture      -0.20222     -0.32091 -0.05780
6Examination      -0.45678     -1.04305  0.34613
7Education        -0.79138     -1.25182 -0.06436
8Catholic          0.10385      0.01947  0.15534
9Infant.Mortality  1.45550      0.87146  2.21101

در این قسمت، رگرسیون چندکی را براساس صدک‌های پنجم تا ۹۵ام با فاصله ۵ درصدی اجرا خواهیم کرد.

1model3 = rq(Fertility~.,data = swiss, tau = seq(0.05,0.95,by = 0.05))
2quantplot = summary(model3)
3quantplot
4plot(quantplot)

خروجی به صورت زیر خواهد بود. البته به علت اینکه ۱۹ صدک‌ مختلف در این کد استفاده شده است، از نمایش همه خروجی‌ها صرف نظر کرده‌ایم و فقط صدک پنجم، پنجاهم و صدک ۹۵ را نمایش خواهیم دارد.

1Call: rq(formula = Fertility ~ ., tau = seq(0.05, 0.95, by = 0.05), 
2    data = swiss)
3tau: [1] 0.05
4Coefficients:
5                 coefficients lower bd  upper bd 
6(Intercept)       51.09373    -86.11344 137.77856
7Agriculture       -0.07513     -0.61030   0.58050
8Examination        0.25953     -1.09496   1.41414
9Education         -1.22152     -2.59070  -0.63000
10Catholic           0.14562     -0.54350   0.21595
11Infant.Mortality   0.77762     -0.93163   4.32563
12...
13Call: rq(formula = Fertility ~ ., tau = seq(0.05, 0.95, by = 0.05), 
14    data = swiss)
15tau: [1] 0.5
16Coefficients:
17                 coefficients lower bd upper bd
18(Intercept)      63.49087     38.04597 87.66320
19Agriculture      -0.20222     -0.32091 -0.05780
20Examination      -0.45678     -1.04305  0.34613
21Education        -0.79138     -1.25182 -0.06436
22Catholic          0.10385      0.01947  0.15534
23Infant.Mortality  1.45550      0.87146  2.21101
24...
25Call: rq(formula = Fertility ~ ., tau = seq(0.05, 0.95, by = 0.05), 
26    data = swiss)
27tau: [1] 0.95
28Coefficients:
29                 coefficients lower bd  upper bd 
30(Intercept)       75.46363     72.06874 101.45314
31Agriculture       -0.08391     -0.86499   0.04293
32Examination       -0.81276     -1.91404   1.38522
33Education         -0.21915     -2.89058   2.24446
34Catholic           0.09827     -0.10995   0.23662
35Infant.Mortality   1.05260     -0.75085   1.18204

نمودارهای حاصل به شکل زیر ترسیم می‌شوند.

در این نمودارها، رگرسیون چندکی به ازاء هر یک از متغیرها مستقل ترسیم شده است. محور افقی، مقادیر چندک‌ها (از صدک پنجم تا صدک ۹۵ام) را نشان می‌دهد. محور عمودی نیز مقدار پارامتر را مشخص کرده است. خطوط سیاه رنگ چندک‌‌ها و فاصله اطمینان آن‌ها نیز با رنگ خاکستری نشان داده شده. رگرسیون خطی ساده نیز با استفاده از خط قرمز و فاصله اطمینان ۹۵٪ برای پارامتر آن نیز با خط‌چین قرمز در نمودار مشخص است.

مثلا برای صدک پنجم، مقدار عرض از مبدا (Intercept)، تقریبا برابر با 51 است و ضریب متغیر تحصیلات (Education) نیز در این حالت ۱٫۲۲- است. ولی این ضرایب برای چندک ۹۵ام تغییر کرده و به 7۵٫46 برای عرض از مبدا و ۰٫۲۱۹ برای تحصیلات رسیده است.

۴. رگرسیون ستیغی

قبل از بررسی «رگرسیون ستیغی» (Ridge Regression)، بهتر است در مورد «قاعده‌سازی» (Regularization) توضیح دهیم. همانطور که قبلا اشاره شد، بیش‌برازش و کم‌برازش، از مسائلی است که ممکن است که در رگرسیون چندگانه رخ دهد. یکی از راه‌های جلوگیری از این مشکلات، قاعده‌سازی است. به این معنی که مدل رگرسیونی را با توجه به تعداد پارامترهای آن جریمه کرد تا تعداد آن‌ها به یک مقدار بهینه برسد. به این ترتیب پیچیدگی مدل کاهش یافته، بدون آنکه از کارایی آن کاسته شود.

قاعده‌سازی بخصوص در موارد زیر کارساز است:

• تعداد زیاد متغیرهای توصیفی
• زیاد بودن تعداد متغیرها نسبت به تعداد مشاهدات
• هم‌خطی یا هم‌خطی چندگانه در بین متغیرهای توصیفی

در رگرسیون ستیغی از تابع زیان درجه ۲ استفاده می‌شود. به این ترتیب مقدار جریمه (Penalty) برای مدل رگرسیونی، به صورت مجموع مربعات ضرایب مشخص می‌شود.

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

$$\large y_{i} = \beta _{0}+\beta _{1}x_{i1} + \cdots +\beta _{p}x_{ip} + \varepsilon _{i} , \qquad i = 1 , \ldots , n$$

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

$$\large \operatorname{argmin} || y - \widehat{y} ||_2^2 = \operatorname{argmin} \sum [y_i - (\beta_0 + \beta_1x_1 + \beta_2x_2 + \cdots + \beta_px_p)]^2$$

توجه داشته باشید که منظور از  $$\operatorname{argmin}$$ مقدارهایی از $$\beta$$ است که تابع مورد نظر را کمینه می‌کنند.

برای برآورد کردن پارامترهای رگرسیونی در روش ستیغی، قیدی روی پارامترها وجود دارد که به صورت زیر نوشته می‌شود.

$$\large \beta_0^2 + \beta_1^2 + \cdots + \beta_p^2 \leq C^2$$

این محدودیت، مشخص می‌کند که باید مجموع مربعات پارامترها از مقدار ثابت یا آستانه‌ای کمتر باشند. به این ترتیب شیوه برآورد پارامترها به صورت زیر در خواهد آمد. مشخص است که بین وجود پارامترهای $$\beta$$ و صفر شدن آن‌ها در بخش قید، توازنی برقرار شده و تعداد پارامترها و متغیرهای مربوطه، بهینه می‌شود.

$$\large \operatorname{argmin} || y - X \beta ||_2^2 + \lambda || \beta ||_2^2$$

پارامتر $$\lambda$$ در اینجا میزان جریمه (Regularization Parameter) نامیده می‌شود.

نکته: توجه داشته باشید که قاعده‌سازی فقط برای پارامترهای $$\beta_1$$ تا $$\beta_n$$ صورت می‌گیرد و عرض از مبدا $$\beta_0$$‌ از این موضوع مستثنی است.

برآورد پارامترهای مدل رگرسیون ستیغی با توجه به قیدی که ذکر شده، به صورت زیر خواهد بود.

$$\large \widehat{\beta}^{ridge} =(X^TX+\lambda I)^{-1}X^Ty$$

نکته: مشخص است که اگر هیچ جریمه‌ای برای مدل در نظر نگیریم، به این معنی که مقدار $$\lambda$$ را صفر انتخاب کنیم، روش برآورد در رگرسیون ستیغی به شیوه OLS تبدیل خواهد شد. در نتیجه می‌توان روش OLS را حالت خاصی از روش ستیغی در نظر گرفت.

در تعریف پارامترها در مدل ستیغی، می‌توان دید که مقدار بزرگ برای $$\lambda$$ باعث مشکل کم‌برازش (Underfitting) می‌شود. بنابراین تعیین مقدار صحیح برای لاندا ($$\lambda$$) ضروری است. برای مشخص کردن مقدار مناسب بهتر است نموداری برحسب مقادیر پارامترها و لاندا ترسیم کرد و زمانی که برآوردگرهای ثابتی بوجود آمده‌اند، کوچکترین مقدار لاندا را در نظر گرفت.

برای پیاده‌سازی رگرسیون ستیغی در زبان برنامه‌نویسی R از مجموعه داده سوئیس استفاده کرده و از کتابخانه glment کمک می‌گیریم.

همچنین تابع cv.glmnet نیز با استفاده از اعتبارسنجی متقابل مدل، مقدار مناسب برای لاندا را مشخص می‌کند. حوزه تغییرات لاندا در این مثال $$10^4$$ تا $$10^{-1}$$‌ است که با گام‌های ۰٫۱، کاهش می‌یابد.

1library(glmnet)
2X = swiss[,-1]
3y = swiss[,1]
4
5set.seed(123) #Setting the seed to get similar results.
6model = cv.glmnet(as.matrix(X),y,alpha = 0,lambda = 10^seq(4,-1,-0.1))
7model
8summary(model)

خروجی به صورت زیر خواهد بود.

1Call:  cv.glmnet(x = as.matrix(X), y = y, lambda = 10^seq(4, -1, -0.1),      alpha = 0) 
2Measure: Mean-Squared Error 
3    Lambda Measure     SE Nonzero
4min  1.585   63.41  8.941       5
51se  7.943   70.64 12.626       5
6
7           Length Class  Mode     
8lambda     51     -none- numeric  
9cvm        51     -none- numeric  
10cvsd       51     -none- numeric  
11cvup       51     -none- numeric  
12cvlo       51     -none- numeric  
13nzero      51     -none- numeric  
14call        5     -none- call     
15name        1     -none- character
16glmnet.fit 12     elnet  list     
17lambda.min  1     -none- numeric  
18lambda.1se  1     -none- numeric

به کمک کد زیر، کمترین مقدار لاندا را با lambda.min بدست آورده و مدل مناسب را مشخص می‌کنیم.

1best_lambda = model$lambda.min
2ridge_coeff = predict(model,s = best_lambda,type = "coefficients")
3ridge_coeff

ضرایب حاصل از اجرای مدل رگرسیون ستیغی در این حالت در ادامه دیده می‌شود.

16 x 1 sparse Matrix of class "dgCMatrix"
2                           1
3(Intercept)      64.92994664
4Agriculture      -0.13619967
5Examination      -0.31024840
6Education        -0.75679979
7Catholic          0.08978917
8Infant.Mortality  1.09527837

۵. رگرسیون لاسو

اصطلاح لاسو (Lasso) مخفف عبارت (Least Absolute Shrinkage and Selection Operator) یا «عمل‌گر گزینش و انقباض کمترین قدرمطلق» است. در این مدل، نحوه قاعده‌سازی براساس تابع زیان قدر مطلق انجام می‌شود. در نتیجه تابع هدف در «رگرسیون لاسو» (Lasso Regression) به صورت زیر نوشته می‌شود.

$$\sum_{i=1}^n(y_i-\beta_0-\sum_{j=1}^px_{ij}\beta_j)^2+\lambda \sum_j|\beta_j|$$

نکته: در اینجا نیز به مانند رگرسیون ستیغی، فرض نرمال بودن باقی‌مانده‌ها وجود ندارد. همچنین مقدار ثابت یا عرض از مبدا هم در قاعد‌ه‌سازی دخیل نمی‌شود.

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

1X = swiss[,-1]
2y = swiss[,1]
3
4set.seed(123)
5model = cv.glmnet(as.matrix(X),y,alpha = 1,lambda = 10^seq(4,-1,-0.1))
6#By default standardize = TRUE
7set.seed(123)
8model = cv.glmnet(as.matrix(X),y,alpha = 1,lambda = 10^seq(4,-1,-0.1))
9model
10summary(model)

باز هم از cv.glmnet برای استفاده از اعتبارسنجی متقابل برای تعیین مقدار مناسب برای پارامتر لاندا استفاده شده است. توجه داشته باشید که برای استفاده از رگرسیون لاسو، مقادیر متغیرها باید استاندارد (Standardize) شده باشند. به طور پیش‌فرض این کار در تابع cv.glmnet صورت می‌گیرد.

خروجی به صورت زیر ظاهر خواهد شد.

1> model
2Call:  cv.glmnet(x = as.matrix(X), y = y, lambda = 10^seq(4, -1, -0.1),      alpha = 1) 
3Measure: Mean-Squared Error 
4    Lambda Measure    SE Nonzero
5min 0.1259   66.05  9.00       5
61se 1.5849   74.24 12.53       4
7> summary(model)
8           Length Class  Mode     
9lambda     51     -none- numeric  
10cvm        51     -none- numeric  
11cvsd       51     -none- numeric  
12cvup       51     -none- numeric  
13cvlo       51     -none- numeric  
14nzero      51     -none- numeric  
15call        5     -none- call     
16name        1     -none- character
17glmnet.fit 12     elnet  list     
18lambda.min  1     -none- numeric  
19lambda.1se  1     -none- numeric

به منظور پیدا کردن بهترین مقدار لاندا از کد زیر استفاده کرده‌ایم.

1best_lambda = model$lambda.min
2lasso_coeff = predict(model,s = best_lambda,type = "coefficients")
3lasso_coeff

در نتیجه مقادیر ضرایب بر این اساس به صورت زیر درخواهد آمد.

16 x 1 sparse Matrix of class "dgCMatrix"
2                           1
3(Intercept)      65.46374579
4Agriculture      -0.14994107
5Examination      -0.24310141
6Education        -0.83632674
7Catholic          0.09913931
8Infant.Mortality  1.07238898

از آنجایی که نتایج حاصل از رگرسیون لاسو با رگرسیون ستیغی نزدیک هستند، سوالی که مطرح می‌شود این است که کدام یک نسبت به دیگری ارجح است؟

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

۶. رگرسیون شبکه الاستیک

«رگرسیون شبکه الاستیک» (Elastic Net Regression)، با ترکیب رگرسیون لاسو و رگرسیون ستیغی، بر معایب آن‌ها غلبه کرده و جایگزین مطمئن برای آن‌ها است. به این ترتیب اگر با مدلی مواجه هستید که متغیرهای توصیفی آن با یکدیگر همبستگی دارند، بهتر است از رگرسیون شبکه الاستیک استفاده کنید.

به این ترتیب یک قاعده‌سازی مرتبه ۱ و ۲ روی مدل همزمان اعمال می‌شود. در نتیجه تابع هدف در رگرسیون شبکه الاستیک به صورت زیر نوشته خواهد شد.

$$\large \min(\sum \epsilon^2 + \lambda_1\sum\beta_i+\lambda_2\sum |\beta_i|)$$

که با در نظر گرفتن مدل رگرسیون خطی چندگانه، می‌توان آن را به صورت زیر نیز در نظر گرفت.

$$\large \min( \sum y_i - ( \beta_0 + \beta_1 X_1 + \beta_2X_2 + \ldots + \beta_kX_k)^2 + \lambda_1 \sum \beta_i^2 + \lambda_2\sum|\beta_i|)$$

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

فایل داده swiss (سوئیس) را در نظر بگیرید. متغیر پاسخ Y براساس ستون اول این مجموعه داده و بقیه به عنوان متغیر X یا توصیفی در نظر گرفته شده‌اند.

مدل رگرسیونی شبکه الاستیک توسط کدهای زیر در نرم‌افزار R پیاده‌سازی شده است.

1library(glmnet)
2X = swiss[,-1]
3y = swiss[,1]
4model = cv.glmnet(as.matrix(X),y,alpha = 0.5,lambda = 10^seq(4,-1,-0.1))
5#Taking the best lambda
6model
7best_lambda = model$lambda.min
8en_coeff = predict(model,s = best_lambda,type = "coefficients")
9en_coeff

مقادیر alpha و lambda از قبل برای تشکیل مدل در نظر گفته شده‌اند که همان پارامترهای اصلی برای مدل رگرسیون سه تیغی و لاسو هستند. حاصل خروجی به صورت زیر خواهد بود.

1> model
2Call:  cv.glmnet(x = as.matrix(X), y = y, lambda = 10^seq(4, -1, -0.1),      alpha = 0.5) 
3
4Measure: Mean-Squared Error 
5
6    Lambda Measure    SE Nonzero
7min  0.398   64.34 15.02       5
81se  3.981   75.70 18.39       4
9> en_coeff
106 x 1 sparse Matrix of class "dgCMatrix"
11                           1
12(Intercept)      64.05919123
13Agriculture      -0.12604667
14Examination      -0.25484517
15Education        -0.77825368
16Catholic          0.09140641
17Infant.Mortality  1.07581457

نتایج حاصل از برآورد پارامترهای مدل را با دو شیوه رگرسیون لاسو و ستیغی مقایسه کنید. به نظر می‌رسد که نتیجه برآورد در رگرسیون شبکه الاستیک به آن‌ها بسیار نزدیک است.

۷. رگرسیون مولفه‌های اصلی

زمانی که هم‌خطی یا هم‌خطی چندگانه در مدل رگرسیونی وجود داشته باشد، بهتر است از مدل رگرسیون مولفه‌های اصلی (Principle Component Regression) که به اختصار با PCR نشان داده می‌شود، استفاده کنیم.

رگرسیون مولفه‌های اصلی در دو گام اجرا می‌شود.

1. استخراج مولفه‌های اصلی براساس متغیرهای توصیفی
2. اجرای رگرسیون براساس مولفه‌های ایجاد شده به عنوان متغیرهای مستقل با متغیر پاسخ

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

بخش اول در محاسبات مربوط به رگرسیون مولفه‌های اصلی، تعیین «بارهای عاملی» (Loads) است که به کمک آن مولفه‌ها ایجاد می‌شوند. هر مولفه (مثل $$U_i$$) به صورت زیر تشکیل می‌شود.

$$\large U_i =  \beta_1X_1 + \beta_2X_2 + \ldots + \beta_pX_p$$

با شرط اینکه بارها ($$\beta_i$$) در شرط زیر صدق کنند.

$$\large \sum_{i=1}^p \beta_i^2 =  1$$

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

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

همچنین از آنجایی که مقدار $$p$$ را می‌توان کمتر یا مساوی با $$k$$ انتخاب کرد، کاهش بعد مسئله نیز از مزایای استفاده از PCR محسوب می‌شود. در نتیجه می‌توان به جای استفاده از مدل با ۱۰ متغیر توصیفی، فقط با ۲ یا ۳ مولفه، مدل رگرسیونی را ایجاد کرد بطوریکه کمترین میزان اطلاعات در مورد متغیر وابسته، در مدل از بین رفته یا نادیده گرفته شده باشد.

نکته: باید این موضوع را در نظر بگیرم که استفاده از PCR، روشی برای تعیین ویژگی‌های موثر در مدل رگرسیونی نیست بلکه با بهره‌گیری از آن، مولفه‌های جدیدی ایجاد می‌شود که بیشترین توصیف یا سهم تغییرات برای متغیر وابسته را در خود دارند. در نتیجه نمی‌توان گفت که کدام متغیر توصیفی، بیشترین نقش را در مدل رگرسیونی PCR دارد.

به منظور پیاده‌سازی رگرسیون PCR از کتابخانه pls در زبان برنامه‌نویسی و محاسباتی R استفاده می‌کنیم. پس ابتدا کتابخانه‌ را نصب و سپس بارگذاری می‌کنیم.

1install.packages("pls")
2library(pls)

در اینجا از مجموعه داده longley استفاده می‌کنیم که در R به طور پیش‌فرض وجود دارد. به جز ستون سال (Year) بقیه متغیرها را استفاده خواهیم کرد.

1data1 = longley[,colnames(longley) != "Year"]
2View(data1)

در ادامه این داده‌ها را مشاهده می‌کنید. قرار است تعداد افراد شاغل (Employed) را براساس متغیرهای دیگر نظیر «جمعیت» (Population)، «نیروی نظامی» (Armed.Forces)، «تعداد افراد بیکار» (Unemployed)، «شاخص درآمد ناخالص ملی» (GNP) و «شاخص درآمد ضمنی ناخالص ملی» (GNP.deflator)، مدل‌سازی کنیم.

1> longley
2     GNP.deflator     GNP Unemployed Armed.Forces Population Year Employed
31947         83.0 234.289      235.6        159.0    107.608 1947   60.323
41948         88.5 259.426      232.5        145.6    108.632 1948   61.122
51949         88.2 258.054      368.2        161.6    109.773 1949   60.171
61950         89.5 284.599      335.1        165.0    110.929 1950   61.187
71951         96.2 328.975      209.9        309.9    112.075 1951   63.221
81952         98.1 346.999      193.2        359.4    113.270 1952   63.639
91953         99.0 365.385      187.0        354.7    115.094 1953   64.989
101954        100.0 363.112      357.8        335.0    116.219 1954   63.761
111955        101.2 397.469      290.4        304.8    117.388 1955   66.019
121956        104.6 419.180      282.2        285.7    118.734 1956   67.857
131957        108.4 442.769      293.6        279.8    120.445 1957   68.169
141958        110.8 444.546      468.1        263.7    121.950 1958   66.513
151959        112.6 482.704      381.3        255.2    123.366 1959   68.655
161960        114.2 502.601      393.1        251.4    125.368 1960   69.564
171961        115.7 518.173      480.6        257.2    127.852 1961   69.331
181962        116.9 554.894      400.7        282.7    130.081 1962   70.551

اگر دسترسی به این مجموعه داده برایتان مقدور نیست، می‌توانید فایل longley را از اینجا با قالب فشرده از اینجا دریافت کنید و پس از خارج کردن از حالت فشرده با دستور read.csv در R‌ بارگذاری کنید.

از آنجایی که بعضی از این متغیرها به یکدیگر وابستگی دارند (مثلا درآمد ناخالص ملی با تعداد افراد بیکار و همچنین درآمد ضمنی ناخالص ملی در رابطه است)، بهتر است از مدل PCR استفاده کنیم. کد زیر به این منظور نوشته شده است. همانطور که مشاهده می‌کنید، مدل ~.Employed نشانگر ارتباط بین متغیر Employed با بقیه متغیرها است.

1pcr_model <- pcr(Employed~., data = data1, scale = TRUE, validation = "CV")
2summary(pcr_model)

خروجی به صورت زیر خواهد بود.

1Data:  X dimension: 16 5 
2 Y dimension: 16 1
3Fit method: svdpc
4Number of components considered: 5
5
6VALIDATION: RMSEP
7Cross-validated using 10 random segments.
8       (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps
9CV           3.627    1.194    1.118   0.5555   0.6514   0.5954
10adjCV        3.627    1.186    1.111   0.5489   0.6381   0.5819
11
12TRAINING: % variance explained
13          1 comps  2 comps  3 comps  4 comps  5 comps
14X           72.19    95.70    99.68    99.98   100.00
15Employed    90.42    91.89    98.32    98.33    98.74

همانطور که مشاهده می‌کنید، حدود 95٫70 درصد تغییرات متغیر وابسته توسط دو مولفه اول و دوم در مدل رگرسیونی PCR پوشش داده می‌شوند. به عنوان شاخص برای اعتبار مدل از RMSEP یا «ریشه میانگین مربعات خطای پیش‌بینی» (Root Mean Square Error Of Prediction) استفاده شده است. همچنین محاسبه اعتبارسنجی متقابل (Crossvalidation) برای سنجش کارایی مدل روی داده‌ها براساس هر یک از مدل‌ها (تعداد مولفه‌ها) صورت گرفته است.

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

1validationplot(pcr_model,val.type = "MSEP")
2
3validationplot(pcr_model,val.type = "R2")
4
5pred = predict(pcr_model,data1,ncomp = 3)

خروجی به صورت زیر خواهد بود.

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

از طرفی نمودار $$R^2$$ نیز نشان می‌دهد که با افزایش تعداد مولفه‌ها، ضریب تعیین افزایش داشته ولی از مولفه سوم به بعد باز هم تغییرات آن چشمگیر نیست. بنابراین اگر سه مولفه اول در مدل به کار گرفته شوند، بهترین دقت و کمترین خطا را به همراه کمترین تعداد مولفه‌ها خواهیم داشت.

حاصل پیش‌بینی‌ها توسط این سه مولفه در ادامه قابل مشاهده است.

1> pred
2, , 3 comps
3
4     Employed
51947 60.14949
61948 61.52509
71949 60.14052
81950 61.28339
91951 63.57905
101952 64.23811
111953 65.08944
121954 63.66381
131955 65.47474
141956 66.68995
151957 67.77152
161958 66.59554
171959 68.67805
181960 69.47458
191961 69.41667
201962 71.30205

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

1firstcomp1=pcr_model$fitted.values[,,1]
2firstcomp2=pcr_model$fitted.values[,,2]
3firstcomp3=pcr_model$fitted.values[,,3]
4
5plot(firstcomp1,data1$Employed,col=1,xlab="Comp Values",ylab="Real Value")
6lines(firstcomp2,data1$Employed,type='p',col=5)
7lines(firstcomp3,data1$Employed,type='p',col=2)
8legend(60,70,legend=c("comp1","comp2","comp3"),col=c(1,5,2),lty=1)

۸. رگرسیون کمترین مربعات جزئی

زمانی که بین متغیرهای توصیفی، وابستگی شدید وجود داشته باشد، به جای رگرسیون مولفه‌های اصلی (PCR) بهتر است از رگرسیون کمترین مربعات جزئی (PLS) یا Partial Least Square Regression استفاده شود. همچنین زمانی که تعداد متغیرهای توصیفی زیاد هستند و می‌خواهیم موثرترین متغیرها در مدل حضور داشته باشند، از رگرسیون کمترین مربعات جزئی (PLS) استفاده می‌کنیم.

هم در روش رگرسیون مولفه‌های اصلی و هم کمترین مربعات جزئی، متغیر جدیدی به عنوان متغیر پیش‌گو ساخته می‌شود که به آن مولفه (Component) گفته می‌شود. این متغیر جدید، ترکیب خطی از متغیرهای توصیفی است. ولی تفاوت در این است که در تحلیل رگرسیون PCR، مولفه‌ها براساس توصیف واریانس کل متغیرهای توصیفی (پیشگو-Predictor) تولید می‌شوند بدون آنکه به مقایر متغیر پاسخ توجه شود. در حالیکه در PLS با در نظر گرفتن متغیر پاسخ و متغیرهای پیشگو، مولفه‌ها تولید می‌شوند و در نهایت مدلی ایجاد می‌شود که با کمترین عوامل، بهترین برازش را دارد.

در کد زیر براساس مجموعه داده vehicles که مشخصات خودروها را ثبت کرده، می‌خواهیم قیمت خودرو (ستون سیزدهم) با نام متغیر price را برحسب متغیرهای دیگر مدل‌بندی کنیم. به منظور دسترسی به این مجموعه داده می‌توانید فایل vehicles را از اینجا به قالب فشرده دریافت کنید و پس از خارج کردن از حالت فشرده، به کمک دستور read.csv آن را در در برنامه R فراخوانی و بارگذاری کنید.

ابتدا بعضی از مقادیر این مجموعه داده را نمایش می‌دهیم. این مجموعه داده، شامل 30 مشاهده و 16 متغیر است. برای ایجاد مدل رگرسیونی PLS از متغیر سیزدهم یعنی قیمت خودرو (price) به عنوان متغیر وابسته استفاده می‌کنیم. همچنین بقیه متغیرها به عنوان متغیرهای پیشگو به کار خواهند رفت.

1> head(vehicles)
2          diesel turbo two.doors hatchback wheel.base length width height curb.weight eng.size
3alfaromeo      0     0         1         1       94.5  171.2  65.5   52.4        2823      152
4audi           0     0         0         0      105.8  192.7  71.4   55.7        2844      136
5bmw            0     0         1         0      101.2  176.8  64.8   54.3        2395      108
6chevrolet      0     0         0         0       94.5  158.8  63.6   52.0        1909       90
7dodge1         0     1         1         1       93.7  157.3  63.8   50.8        2128       98
8dodge2         0     0         0         1       93.7  157.3  63.8   50.6        1967       90
9          horsepower peak.rpm price symbol city.mpg highway.mpg
10alfaromeo        154     5000 16500      1       19          26
11audi             110     5500 17710      1       19          25
12bmw              101     5800 16430      2       23          29
13chevrolet         70     5400  6575      0       38          43
14dodge1           102     5500  7957      1       24          30
15dodge2            68     5500  6229      1       31          38

توجه داشته باشید که برای اجرای رگرسیون PLS باید از کتابخانه plsdepot استفاده کنید. مجموعه داده vehicles نیز در این کتابخانه قرار دارد.

1library(plsdepot)
2data(vehicles)
3pls.model = plsreg1(vehicles[, c(1:12,14:16)], vehicles[, 13], comps = 3)
4# summary
5summary(pls.model)
6# Predicted values
7pls.model$y.pred
8# loads of factors
9pls.model$x.loads
10# R-Square
11pls.model$R2

همانطور که مشاهده می‌کنید، رگرسیون PLS‌ از سه مولفه استفاده کرده است زیرا در کد عبارت comps=3 به کار رفته است. به این ترتیب خروجی به صورت زیر خواهد بود.

1> summary(pls.model)
2          Length Class  Mode   
3x.scores  90     -none- numeric
4x.loads   45     -none- numeric
5y.scores  90     -none- numeric
6y.loads    3     -none- numeric
7cor.xyt   48     -none- numeric
8raw.wgs   45     -none- numeric
9mod.wgs   45     -none- numeric
10std.coefs 15     -none- numeric
11reg.coefs 16     -none- numeric
12R2         3     -none- numeric
13R2Xy      48     -none- numeric
14y.pred    30     -none- numeric
15resid     30     -none- numeric
16T2        93     -none- numeric
17Q2        15     -none- numeric
18y         30     -none- numeric
19> # Predicted values
20> pls.model$y.pred
21         1          2          3          4          5          6          7          8          9 
2218186.8186 17880.1674 14547.1022  5990.3904  8805.9785  4357.6998  5545.7037 10140.0942  7306.1875 
23        10         11         12         13         14         15         16         17         18 
2436512.4736  5221.2953 25798.8095 17472.6262  6653.6606  6809.4727  6977.2829  4357.6998 17427.1386 
25        19         20         21         22         23         24         25         26         27 
2633931.7662 14254.0017  8285.7662   784.7039 10960.9806 15410.8928 10749.3651 11546.7997 14782.5522 
27        28         29         30 
2813053.9396 18398.4627 19615.1680 
29> # loads of factors
30> pls.model$x.loads
31                     p1         p2           p3
32diesel       0.09271495 -0.3890200  0.678689303
33turbo        0.18536097 -0.2436168 -0.086493899
34two.doors   -0.10124057  0.3717346 -0.039073060
35hatchback   -0.14373397  0.2481944 -0.480716370
36wheel.base   0.32361239 -0.2729353 -0.018798612
37length       0.35328772 -0.1794076  0.022031790
38width        0.33579957 -0.1300531  0.008125624
39height       0.20723718 -0.3762250 -0.099462195
40curb.weight  0.37681297 -0.0261901  0.117117365
41eng.size     0.34243554  0.1890635  0.252856936
42horsepower   0.29507625  0.3616392 -0.081524076
43peak.rpm    -0.11417298  0.3743210 -0.273045978
44symbol      -0.14375163  0.4326949  0.131893643
45city.mpg    -0.31971248 -0.2441165  0.352607126
46highway.mpg -0.34605016 -0.1329582  0.317595249
47> # R-Square
48> pls.model$R2
49        t1         t2         t3 
500.70474894 0.11708109 0.09872169

در بخش اول، متغیرها و مقادیرشان معرفی شده‌اند. در بخش دوم، مقدارهای پیش‌بینی شده توسط مدل برای متغیر وابسته با دستور pls.model$y.pred ظاهر شده است. همچنین بارهای عاملی برای تشکیل مولفه‌ها، براساس سه مولفه اول مشخص شده. به این ترتیب ضریب هر یک از متغیرهای توصیفی که در اینجا همان عامل‌ها نامیده می‌شوند، در تشکیل مدل رگرسیونی، دیده می‌شوند.

در انتها نیز ضریب $$R^2$$ برای هر یک از مولفه‌ها مشخص شده است. همانطور که انتظار داریم، اولین مولفه، بیشترین همبستگی را با متغیر وابسته دارد. به همین ترتیب برای مولفه‌های دوم و سوم، ضریب همبستگی کاهش یافته است.

۹. رگرسیون بردار پشتیبان

به کمک رگرسیون بردار پشتیبان (Support Vector Regression)  که گاهی آن را با SVR نیز نشان می‌دهند، می‌توان مدل‌های خطی و غیرخطی را ایجاد و پارامترهای آن را محاسبه کرد. این کار توسط به کارگیری یک «تابع هسته» (kernel) غیرخطی (مانند چندجمله‌ای) حاصل می‌شود. محاسبه پارامترهای این تابع به شکل است که خطا کمینه شود بطوری که فاصله بین صفحاتی که عمل جداسازی بین دسته‌ها را ایجاد می‌کنند، بیشینه شود.

کتابخانه e1071 از زبان برنامه‌نویسی R، امکان اجرای رگرسیون بردار پشتیبان را می‌دهد. در کدی که در ادامه مشاهده می‌کنید، با استفاده از داده‌های شبیه‌سازی شده، یک مدل رگرسیونی بردار پشتیبان ایجاد می‌شود.

متغیر مستقل (x) براساس دنباله‌ای از مقادیر ۰ تا ۵ با افزایش ۰٫۰۵ تولید شده است. متغیر وابسته (y) نیز براساس رابطه زیر ساخته شده است.

$$\large y = \log(x) + z, \; \; \; z \sim N(mean = x , standard deviation = 0.2)$$

پارامترها و کد دستوری برای اجرای این مدل رگرسیونی در ادامه قابل مشاهده است.

1svm.default(x = x, y = y)
2Parameters:
3   SVM-Type:  eps-regression 
4 SVM-Kernel:  radial 
5       cost:  1 
6      gamma:  1 
7    epsilon:  0.1 
8
9Number of Support Vectors:  71

کدی که در ادامه مشاهده می‌کنید، برای اجرای رگرسیون SVM برای داده‌های شبیه‌سازی شده به کار می‌رود.  پس از اجرای مدل، مقادیر پیش‌بینی شده نیز محاسبه شده و به همراه مقادیر x و y در یک نمودار ترسیم می‌شود.

1library(e1071)
2# create data
3x <- seq(0.1, 5, by = 0.05)
4y <- log(x) + rnorm(x, sd = 0.2)
5
6# estimate model and predict input values
7m   <- svm(x, y)
8new <- predict(m, x)
9
10# visualize
11plot(x, y)
12points(x, log(x), col = 2)
13points(x, new, col = 4)

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

نقطه‌های مشکی رنگ، بیانگر نقاط x و y است. نقطه‌های قرمز رنگ هم نتیجه ارتباط بین نقاط x با (log(x را نشان می‌دهد. قرار است ارتباط بین y با (log(x را به صورت خطی نمایش دهیم. مدل برازش شده توسط رگرسیون SVM نیز با نقطه‌های آبی رنگ روی نمودار دیده می‌شود که انطباق مناسبی با نقطه‌های قرمز رنگ دارد. البته در مقادیر نزدیک به صفر این مطابقت از بین رفته است زیرا تغییرات لگاریتم در این قسمت شدید خواهد بود.

۱۰. رگرسیون ترتیبی

اگر متغیر پاسخ یا وابسته به صورت ترتیبی (Ordinal) یا رتبه‌ای (Rank) باشد، بهتر است از مدل رگرسیون ترتیبی (Ordinal Regression) استفاده شود. به این ترتیب اگر با استفاده از طیف لیکرت (Likert Scale)، با پاسخ‌های مربوط به پرسشنامه مواجه هستید، بهتر است برای مدل‌سازی بین این متغیر پاسخ با متغیرهای مستقل یا توصیفی، از رگرسیون ترتیبی استفاده کنید. برای مثال ممکن است براساس یک پرسش‌نامه که پاسخ‌های آن با مقیاس یا طیف لیکرت ۵ سطحی تشکیل شده بخواهید میزان بهبود درد بیماران نسبت به دوزهای مختلف یک دارو را بسنجید. چنین مدلی به کمک رگرسیون ترتیبی ساخته خواهد شد.

به منظور پیاده‌سازی این مدل از رگرسیون با مجموعه داده wine که میزان تندی سرکه را مشخص کرده است کار خواهیم کرد. این مجموعه داده، شامل ۷۲ مشاهده است که براساس قضاوت ۹ داور میزان تندی با متغیر رتبه‌ای response از ۰ تا ۱۰۰ مشخص و درجه‌بندی شده است. همچنین متغیر rating نیز مقادیر response را به ۵ طبقه به صورت ترتیبی، گروه‌بندی کرده است. دما (temp) با دو سطح گرم و سرد و همچنین نحوه چیدن انگور با متغیر contact با دوسطح no و yes ثبت شده. نوع بطری (bottle) نیز نحوه نگهداری سرکه را با توجه به هشت وضعیت مختلف تعیین می‌کند.

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

1> head(wine)
2  response rating temp contact bottle judge
31       36      2 cold      no      1     1
42       48      3 cold      no      2     1
53       47      3 cold     yes      3     1
64       67      4 cold     yes      4     1
75       77      4 warm      no      5     1
86       60      4 warm      no      6     1

در اینجا، مدل ارائه شده طبق کد زیر به بررسی رابطه بین متغیر وابسته rating به عنوان متغیر عامل یا ترتیبی با داورها (judge) و دما (temp) پرداخته‌ایم.

1library(ordinal)
2o.model <- clm( rating~ temp*judge, data = wine)
3summary(o.model)

خروجی به صورت زیر خواهد بود.

1formula: rating ~ temp * judge
2data:    wine
3
4 link  threshold nobs logLik AIC    niter max.grad cond.H 
5 logit flexible  72   -73.79 189.58 6(0)  3.99e-07 4.7e+02
6
7Coefficients:
8                  Estimate Std. Error z value Pr(>|z|)   
9tempwarm         3.683e+00  1.465e+00   2.514  0.01194 * 
10judge2          -4.132e+00  1.643e+00  -2.515  0.01189 * 
11judge3          -1.357e+00  1.360e+00  -0.998  0.31838   
12judge4          -1.357e+00  1.360e+00  -0.998  0.31838   
13judge5          -3.433e-16  1.390e+00   0.000  1.00000   
14judge6          -1.357e+00  1.360e+00  -0.998  0.31838   
15judge7          -4.720e+00  1.579e+00  -2.990  0.00279 **
16judge8          -2.152e+00  1.403e+00  -1.534  0.12495   
17judge9          -2.995e+00  1.549e+00  -1.933  0.05320 . 
18tempwarm:judge2  1.703e+00  2.218e+00   0.768  0.44259   
19tempwarm:judge3  1.357e+00  1.893e+00   0.717  0.47359   
20tempwarm:judge4 -2.086e+00  2.022e+00  -1.031  0.30240   
21tempwarm:judge5 -3.663e+00  1.970e+00  -1.859  0.06301 . 
22tempwarm:judge6 -1.402e-01  1.950e+00  -0.072  0.94269   
23tempwarm:judge7 -3.202e-01  2.094e+00  -0.153  0.87850   
24tempwarm:judge8 -8.934e-01  1.956e+00  -0.457  0.64786   
25tempwarm:judge9  1.057e-01  2.102e+00   0.050  0.95990   
26---
27Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
28
29Threshold coefficients:
30    Estimate Std. Error z value
311|2   -4.656      1.271  -3.663
322|3   -1.302      1.013  -1.285
333|4    1.341      1.013   1.325
344|5    3.507      1.171   2.994

۱۱. رگرسیون پواسون

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

متغیر وابسته در این مدل باید شرایط زیر را داشته باشد:

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

به منظور اجرای مدل رگرسیون پواسون از یک مجموعه داده که مربوط به پارگی دو نوع نخ است استفاده می‌کنیم. ابتدا به داده‌های warpbreaks نگاهی می‌اندازیم. این مجموعه داده شامل ۵۴ مشاهده و ۳ متغیر است.

• تعداد پارگی‌ها (breaks) که مقادیر آن به صورت شمارشی است.
• نوع نخ (wool) که با دو مقدار A و B‌ مشخص شده است.
• میزان کشش (tension) که در سه سطح کم (L)، متوسط (M) و زیاد (H) طبقه‌بندی شده است.

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

1> head(warpbreaks)
2  breaks wool tension
31     26    A       L
42     30    A       L
53     54    A       L
64     25    A       L
75     70    A       L
86     52    A       L

کدی که در ادامه مشاهده می‌کنید، مربوط به مجموعه داده استحکام نخ‌ها است. تعداد پارگی‌های نخ‌ها (breaks) برحسب نوع نخ (wool)، کشش نخ (tension) و اثر متقابل این دو براساس یک مدل رگرسیونی پواسون، مورد محاسبه قرار گرفته و پارامترهای آن بدست آمده است.

1pos.model<-glm(breaks~wool*tension, data = warpbreaks, family=poisson)
2summary(pos.model)

از آنجایی که در مدل عبارت wool*tension دیده می‌شود، اثرات اصلی و متقابل بین این دو متغیر در مدل نیز حضور خواهند داشت. خروجی به صورت زیر خواهد بود.

1glm(formula = breaks ~ wool * tension, family = poisson, data = warpbreaks)
2
3Deviance Residuals: 
4    Min       1Q   Median       3Q      Max  
5-3.3383  -1.4844  -0.1291   1.1725   3.5153  
6
7Coefficients:
8               Estimate Std. Error z value Pr(>|z|)    
9(Intercept)     3.79674    0.04994  76.030  < 2e-16 ***
10woolB          -0.45663    0.08019  -5.694 1.24e-08 ***
11tensionM       -0.61868    0.08440  -7.330 2.30e-13 ***
12tensionH       -0.59580    0.08378  -7.112 1.15e-12 ***
13woolB:tensionM  0.63818    0.12215   5.224 1.75e-07 ***
14woolB:tensionH  0.18836    0.12990   1.450    0.147    
15---
16Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
17
18(Dispersion parameter for poisson family taken to be 1)
19
20    Null deviance: 297.37  on 53  degrees of freedom
21Residual deviance: 182.31  on 48  degrees of freedom
22AIC: 468.97
23
24Number of Fisher Scoring iterations: 4

همانطور که مشخص است همه متغیرها (یا ضرایب) در مدل از لحاظ آماری معنی‌دار بوده و فقط نوع نخ B با میزان کشش سطح H در سطح آزمون ۵٪ معنی‌دار نیستند.

از طرفی دیده می‌شود که میانگین مقدار متغیر پاسخ (Intercept) برابر با 3٫79674 است در حالیکه مقدار پراکندگی (Dispersion) برای توزیع پواسن عدد ۱ در نظر گرفته شده است. در نتیجه با مشکل «کم‌پراکنش» (Underdispersion) مواجه هستیم و بهتر است با مدل رگرسیون دوجمله‌ای منفی یا شبه پواسن مدل را ایجاد کنیم. این مدل در ادامه معرفی و به کار گرفته شده است.

۱۲. رگرسیون دوجمله‌ای منفی

برای فراگیری رگرسیون دو جمله‌ای منفی (Negative Binomial Regression) و دستورات مرتبط با آن در R به بررسی یک مجموعه داده به نام quine‌ می‌پردازیم که مربوط به تعداد روزهای غیبت دانش‌آموزان با توجه به قومیت، سن، جنس و سرعت یادگیری آن‌ها است.

این مجموعه داده شامل 146 مشاهده و ۵ متغیر است که از بین مدارس چهار منطقه متفاوت جمع‌آوری شده است.

1. قومیت (Eth) که شامل دو مقدار بومی (A) و غیربومی (N) است.
2. جنسیت (Sex) که شامل دو مقدار دختر (F) و پسر (M) است.
3. سن (Age) که به طبقه‌های F0 تا F3 تقسیم‌بندی شده است.
4. سرعت یادگیری (Lrn) که شامل دو سطح متوسط (AL) و کند (SL) می‌شود.
5. تعداد روزهای غیبت در مدرسه در طول سال (Days) که یک متغیر شمارشی است و به عنوان متغیر وابسته به کار رفته است.

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

1library(MASS)
2> head(quine)
3  Eth Sex Age Lrn Days
41   A   M  F0  SL    2
52   A   M  F0  SL   11
63   A   M  F0  SL   14
74   A   M  F0  AL    5
85   A   M  F0  AL    5
96   A   M  F0  AL   13

قرار است به کمک متغیرهای توصیفی که در ستون های اول تا چهارم قرار دارند، مدلی برای پیش‌بینی متغیر پنجم یعنی تعداد روزهای غیبت در مدرسه (Days)، ایجاد کنیم. کد زیر به این منظور نوشته شده است.

1nb.model <- glm.nb(Days ~ Sex/(Age + Eth*Lrn), data = quine)
2summary(nb.model)

همانطور که مشاهده می‌کنید در اینجا از تابع glm.nb استفاده کرده‌ایم تا مدل رگرسیون دوجمله‌ای منفی را به کار ببریم. در اینجا متغیر وابسته Days و متغیرهای مستقل Sex, Age, Eth و Lrn هستند که البته از متغیر Sex برای تفکیک مدل‌ها استفاده شده است.

خروجی به صورت زیر خواهد بود.

1glm.nb(formula = Days ~ Sex/(Age + Eth * Lrn), data = quine, 
2    init.theta = 1.597990733, link = log)
3
4Deviance Residuals: 
5    Min       1Q   Median       3Q      Max  
6-2.8950  -0.8827  -0.2299   0.5669   2.1071  
7
8Coefficients:
9                Estimate Std. Error z value Pr(>|z|)    
10(Intercept)      3.01919    0.29706  10.163  < 2e-16 ***
11SexM            -0.47541    0.39550  -1.202 0.229355    
12SexF:AgeF1      -0.70887    0.32321  -2.193 0.028290 *  
13SexM:AgeF1      -0.72373    0.33048  -2.190 0.028527 *  
14SexF:AgeF2      -0.61486    0.37141  -1.655 0.097826 .  
15SexM:AgeF2       0.62842    0.27366   2.296 0.021655 *  
16SexF:AgeF3      -0.34235    0.32717  -1.046 0.295388    
17SexM:AgeF3       1.15084    0.31385   3.667 0.000246 ***
18SexF:EthN       -0.07312    0.26539  -0.276 0.782908    
19SexM:EthN       -0.67899    0.25631  -2.649 0.008072 ** 
20SexF:LrnSL       0.94358    0.32246   2.926 0.003432 ** 
21SexM:LrnSL       0.23891    0.33553   0.712 0.476443    
22SexF:EthN:LrnSL -1.35849    0.37719  -3.602 0.000316 ***
23SexM:EthN:LrnSL  0.76142    0.44134   1.725 0.084483 .  
24---
25Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
26
27(Dispersion parameter for Negative Binomial(1.598) family taken to be 1)
28
29    Null deviance: 234.56  on 145  degrees of freedom
30Residual deviance: 167.56  on 132  degrees of freedom
31AIC: 1093
32
33Number of Fisher Scoring iterations: 1
34
35
36              Theta:  1.598 
37          Std. Err.:  0.213 
38
39 2 x log-likelihood:  -1063.025

همانطور که مشاهده می‌کنید برای هر متغیر مستقل، در دو وضعیت پسر و دختر (M و F)، ضرایب به طور مجزا محاسبه شده‌اند. توجه داشته باشید که عرض از مبدا (Intercept) برای گروه دخترها محاسبه شده. برای پسر‌ها باید مقدار محاسبه شده در ضریب SexM به مقدار Intercept اضافه شده تا عرض از مبدا برای مدل رگرسیون دوجمله‌ای منفی در گروه پسرها، حاصل شود.

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

1pos.model<-glm(breaks~wool*tension, data = warpbreaks, family=quasipoisson)
2summary(pos.model

خروجی به صورت زیر خواهد بود.

1glm(formula = breaks ~ wool * tension, family = quasipoisson, 
2    data = warpbreaks)
3
4Deviance Residuals: 
5    Min       1Q   Median       3Q      Max  
6-3.3383  -1.4844  -0.1291   1.1725   3.5153  
7
8Coefficients:
9               Estimate Std. Error t value Pr(>|t|)    
10(Intercept)     3.79674    0.09688  39.189  < 2e-16 ***
11woolB          -0.45663    0.15558  -2.935 0.005105 ** 
12tensionM       -0.61868    0.16374  -3.778 0.000436 ***
13tensionH       -0.59580    0.16253  -3.666 0.000616 ***
14woolB:tensionM  0.63818    0.23699   2.693 0.009727 ** 
15woolB:tensionH  0.18836    0.25201   0.747 0.458436    
16---
17Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
18
19(Dispersion parameter for quasipoisson family taken to be 3.76389)
20
21    Null deviance: 297.37  on 53  degrees of freedom
22Residual deviance: 182.31  on 48  degrees of freedom
23AIC: NA
24
25Number of Fisher Scoring iterations: 4

همانطور که مشاهده می‌کنید در قسمت Dispersion parameter، مقدار واریانس 3٫76389 محاسبه شده که با میانگین مدل (Intercept) با مقدار 3٫79674 تقریبا برابر است و مشکل بیش‌پراکندگی یا کم‌پراکندگی بوجود نیامده است.

۱۳. رگرسیون شبه پواسون

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

به منظور پیاده‌سازی این مدل رگرسیونی در زبان برنامه‌نویسی R باز هم از تابع glm‌ استفاده می‌کنیم ولی مدل را در قسمت Family برابر با quasipoisson قرار می‌دهیم. مجموعه داده به کار رفته همان quine است که در رگرسیون دو جمله‌ای منفی به کار بردیم.

1qs.pos.model <- glm(Days ~ Sex/(Age + Eth*Lrn), data = quine, family = "quasipoisson")

در کد بالا، مدل رگرسیون شبه پواسن را براساس متغیرهای وابسته و مستقل به شکل زیر اجرا کرده‌ایم.

اثر متغیرهای طبقه‌ای سن، به همراه اثرات اصلی (Main Effect) و اثرات متقابل (Interaction) بین متغیرهای قومیت (Eth) و سرعت یادگیری (Lrn) به تفکیک جنسیت (Sex) برای متغیر روز‌های غیبت (Days) مدل‌سازی شده است.

1glm(formula = Days ~ Sex/(Age + Eth * Lrn), family = "quasipoisson", 
2    data = quine)
3
4Deviance Residuals: 
5    Min       1Q   Median       3Q      Max  
6-7.3924  -2.4923  -0.5745   1.5507   8.5631  
7
8Coefficients:
9                Estimate Std. Error t value Pr(>|t|)    
10(Intercept)      2.89964    0.28113  10.314  < 2e-16 ***
11SexM            -0.39728    0.38999  -1.019 0.310204    
12SexF:AgeF1      -0.48941    0.31188  -1.569 0.118988    
13SexM:AgeF1      -0.82267    0.40664  -2.023 0.045081 *  
14SexF:AgeF2      -0.19579    0.34387  -0.569 0.570071    
15SexM:AgeF2       0.62626    0.26615   2.353 0.020096 *  
16SexF:AgeF3      -0.21627    0.31725  -0.682 0.496619    
17SexM:AgeF3       1.01691    0.30175   3.370 0.000986 ***
18SexF:EthN       -0.08595    0.27661  -0.311 0.756503    
19SexM:EthN       -0.49065    0.24160  -2.031 0.044278 *  
20SexF:LrnSL       0.73625    0.30990   2.376 0.018950 *  
21SexM:LrnSL       0.33779    0.29888   1.130 0.260451    
22SexF:EthN:LrnSL -1.27702    0.41039  -3.112 0.002280 ** 
23SexM:EthN:LrnSL  0.59021    0.40165   1.469 0.144083    
24---
25Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
26
27(Dispersion parameter for quasipoisson family taken to be 10.53364)
28
29    Null deviance: 2073.5  on 145  degrees of freedom
30Residual deviance: 1393.5  on 132  degrees of freedom
31AIC: NA
32
33Number of Fisher Scoring iterations: 5

ضرایبی که با علامت * مشخص شده‌اند، از لحاظ آماری معنی‌دار بوده و حضورشان در مدل ضروری است.

۱۴. رگرسیون کاکس

از «رگرسیون کاکس» (Cox Regression) برای داده‌های وابسته به زمان استفاده می‌شود. برای درک بیشتر این موضوع به مثال‌های زیر توجه کنید:

1. زمان طی شده از وقتی که کاربر وارد حساب خود شده تا زمانی که خارج می‌شود.
2. طول دوره درمان بیماری سرطان تا مرگ بیمار.
3. زمان بین اولین و دومین حمله قلبی.

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

چنین تحلیل‌هایی را به عنوان «تحلیل بقا» (Survival Analysis) می‌شناسیم. همانطور که مشخص است متغیر وابسته از دو بخش تشکیل شده. در بخش اول که مقداری پیوسته است، زمان یا طول دوره رخداد یک پیشامد ثبت شده و در بخش دوم،‌ رخداد یا عدم رخداد آن پیشامد مشخص شده است.

کد زیر به منظور بررسی رابطه بین طول عمر بیماران سرطانی (SurvObj) و عوامل دیگر مانند سن (age)، جنس (sex) و امتیاز کارایی کارنوفکسی ( Karnofsky performance score) که با ph.karno قابل مشاهده است، به همراه میزان کاهش وزن (wt.loss) تهیه شده است. برای دسترسی به این مجموعه داده باید کتابخانه survival را راه‌اندازی کنید.

نکته: برای دسترسی مستقیم به فایل مربوط به این داده‌ها می‌توانید اطلاعات و مقادیر متغیرها را از اینجا در قالب یک فایل فشرده دریافت کنید و بعد از خارج کردن فایل از حالت زیپ، آن را با دستور read.csv بازیابی کنید.

1library(survival)
2# Lung Cancer Data
3# status: 2=death
4lung$SurvObj <- with(lung, Surv(time, status == 2))
5cox.reg <- coxph(SurvObj ~ age + sex + ph.karno + wt.loss, data = lung)
6cox.reg