رگرسیون

رگرسیون خطی ساده

فرض کنید 2 بردار y و x به حجم n وجود دارد و به نظر می‌رسد که بین این دو نمونه رابطه‌ای خطی وجود دارد. با استفاده از تابع ()lm رگرسیون موجود بین این 2 بردار را بدست می‌آوریم.

بطور مثال فرض کنید دو بردار x و y زیر موجودند:

(set.seed(50

(x<-rnorm(1000,0,1

(y<-rnorm(1000,3,1

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

y_i=a+bx_i+e_i

 (lm(y~x

را اجرا می‌کنیم که براورد عرض از مبدا (intercept) و شیب خط رگرسیون (که در مدل رگرسیونی فوق به ترتیب مقادیر a و b) هستند به صورت زیر نشان داده می‌شود:

:Call

(lm(formula = y ~ x

 :Coefficients

Intercept)          x) 

      0.008005-        3.018635

رگرسیون خطی چندگانه

فرض کنید متغیر پاسخ (y) و چند متغیر پیشگو (برای مثال u، v و w) موجود است و تصور می‌کنیم که بین این متغیرهای پیشگو و متغیر پاسخ ارتباط خطی وجود دارد.

بطور مثال فرض کنید دو بردار u، v، w و y زیر موجودند:

(set.seed(500

(u<-rnorm(1000,1,1

(v<-rnorm(1000,3,1

(w<-rnorm(1000,1,1

(y<-rnorm(1000,10,1

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

y_i=a₀+a₁u_i+ a₂v_i + a₃w_i + e_i

 (lm(y~u+v+w

بدست آوردن آماره‌های رگرسیون

فرض کنید درباره رگرسیونی که بدست آورده‌ایم، می‌خواهیم اطلاعاتی چون R₂، آماره F، فواصل اطمینان برای ضرائب رگرسیون، مانده‌ها، جدول ANOVA و... را بدست آوریم.

برای این منظور پس از اجرای رگرسیون و نامگذاری دلخواه آن (مثلاً برای راحتی کار نام متغیر را m در نظر می‌گیریم)، برای

-جدول ANOVA، تابع (anova(m؛

-ضرائب مدل، (coefficient(m یا (coef(m؛

-فواصل اطمینان ضرائب مدل، (confint(m؛

-مجموع مربعات خطا، (deviance(m؛

-بردار اثرات متعامد، (effects(m؛

-بردار مقادیر برازیده شده متغیر پاسخ (y، fitted(m؛

- مانده‌ها، (residuals(m یا (resid(m؛

-خلاصه آماره‌ها چون، R₂، F و خطای استاندارد مانده‌ها، (summary(m؛

-ماتریس واریانس – کوواریانس پارامترهای اصلی، (vcov(m.

تشخیص یک مدل رگرسیون خطی

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

برای این منظور ابتدا تابع plot را برای مدل مورد نظر اجرا می‌کنیم که برخی نمودارهای را جهت تشخیص این موضوع ارائه می‌کند.

بعد از نصب پکیج car می‌توان از تابع ()outlierTest برای بررسی مانده‌های مدل استفاده کرد.

با در نظر گرفتن مثال زیر به ارائه مختصری از کار فوق می‌پردازیم.

فرض کنید:

(set.seed(1500

(x<-rnorm(1000,0,1

(y<-rnorm(1000,00,1

(m<- lm(y~x

پس از اجرای تایع:

(plot(m

به طور مثال نمودار اول بیانگر پراکنش مقادیر برازیده شده در برابر مانده‌هاست. همچنین دومین نمودار بیانگر Q-Q plot است که نشان می‌دهد آیا مانده‌ها از تابع توزیع نرمال پیروی می‌کنند یا خیر.

 پیشگویی مقادیر جدید

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

برای مثال اگر:

(set.seed(2500

(x<-rnorm(1000,0,1

(y<-rnorm(1000, 0,1

(m<- lm(y~x

داده‌های مدل رگرسیونی باشد و

(u<-rnorm(1000,0,1

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

(pred<-data.frame(u

ساخته و سپس با استفاده از دستور:

(predict(m,pred

مقادیر را پیشگویی می‌کنیم.

ساخت بافت‌نگار

برای این منظور از تابع ()hist استفاده می‌کنیم. آرگومان مربوط به این تابع نیز مقادیر عددی را در بر می‌گیرد.

برای مثال، فرض کنید نمودار جعبه‌ای مربوط به نمونه‌ای از اعداد نرمال با میانگین 0 و واریانس 1 به حجم1000 را می‌خواهیم رسم کنیم:

(set.seed(100

(x<-rnorm(1000,0,1

(hist(x

می‌توان تعداد ستون‌هایی که در بافت‌نگار مورد نظر است را انتخاب کرد. همچنین می‌توان عنوان نمودار و نام متغیر سطر افقی را می‌توان به ترتیب با استفاده از پارامترهای main و xlab تغییر داد.

به‌طور مثال می‌توان عنوان نمودار و نام متغیر سطر  افقی را به random normal و values تغییر داد. همچنین تعداد ستون‌ها را 100 در نظر گرفت. بنابراین خواهیم داشت:

("hist(x,100,main="random normal", xlab="valuse

برای رسم بافت‌نگار همچنین می‌توانید به تابع histogram از پکیج lattice مراجعه کنید.

ساخت نمودار جعبه‌ای

برای این منظور کافی است از تابع boxplot استفاده کرد. آرگومان مربوط به متغیری که می‌خواهیم نمودار جعبه‌ای آن ترسیم شود از نوع عددی (numeric) است.

برای مثال، فرض کنید نمودار جعبه‌ای مربوط به نمونه‌ای از اعداد نرمال با میانگین 0 و واریانس 1 به حجم1000 را می‌خواهیم رسم کنیم:

(set.seed(100

(x<-rnorm(1000,0,1

(boxplot(x

سطر اول دستور به منظور یکسان بودن نمونه اعداد تصادفی انتخاب شده از جامعه نرمال با میانگین 0 و واریانس 1 است. شکل نمودار جعبه‌ای به صورت زیر است:

در توضیح شکل باید گفت که:

- خط پررنگی که در وسط جعبه وجود دارد، میانه جامعه است.

-در جعبه موجود سطر پایین چارک اول مشاهدات و سطر بالا، چارک سوم مشاهدات است.

-ارتفاع کادری که حول جعبه است نشان‌دهنده دامنه مشاهدات است.

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

ساخت نمودار میله‌ای

برای ساخت نمودار میله‌ای از دستور ()barplot استفاده می‌کنیم. برای این منظور مقادیری را که می‌خواهیم نمودار میله‌ای آن رسم شوند، با یک بردار مشخص کرده و آنرا داخل آرگومان تابع قرار می‌دهیم. مثلاً اگر بخواهیم نمودار میله‌ای را برای بردار

(c(3,4,4.5,3.5

رسم کنیم، کافی است دستور زیر را اجرا کنیم:

barplot(c(3,4,4.5,3.5))

در صورتی که بخواهیم نام مقادیر را برای سطر افقی مشخص کنیم،‌ می‌توان پارامتر names.arg را برای هر کدام از این مقادیر مشخص کرد. در مثال قبل اگر بخواهیم عنوان هر یک از مقادیر بردار را مشخص کنیم (با مقادیر A، B، C و D) کافی است دستور زیر را اجرا کنیم:

barplot(c(3,4,4.5,3.5),names.arg=c("A","B","C","D"))

در صورتی که بخواهیم برای محور عمودی یا افقی برچسب مشخص کنیم، به ترتیب پارامتر ylab و xlab را مقداردهی می‌کنیم. در مثال قبل اگر بخواهیم برای محور افقی و محور عمودی به ترتیب برچسب‌های Names و Values را مشخص کنیم خواهیم داشت:

barplot(c(3,4,4.5,3.5),names.arg=c("A","B","C","D"),xlab="Names",ylab="Values")

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

فرض کنید در مجموعه داده airquality می‌خواهیم برای مقادیر متغیر Temp با توجه به مقادیر گسسته Month نمودار میله‌ای رسم کنیم. با توجه به دستور ()tapply و با توجه به شاخص آماری میانگین خواهیم داشت:

tapply(airquality$Temp,airquality$month,mean)

حال نمودار میله‌ای را پس از بدست آمدن میانگین مقادیر متغیر Temp با توجه به مقادیر گسسته Month رسم می‌کنیم که دستور آن به صورت زیر است:

barplot(tapply(airquality$Temp,airquality$Month,mean))

ساخت نمودار پراکنش

فرض کنید برای مشاهدات زوجی (x₁,y₁)، (x₂,y₂)، ... و (x_n,y_n) می‌خواهیم نمودار پراکنش رسم کنیم. برای این منظور از دستور plot استفاده می‌کنیم. اگر دو بردار هم‌اندازه باشند (مثلاً بردار x و y) آنگاه باید دستور را به صورت (plot(x,y اجرا کرد. در ضورتی که مجموعه داده‌ای (datafram) حاوی دو ستون داشته باشیم آنگاه کافی است نام مجموعه داده را در تابع فوق نوشت و دستور را اجرا کرد یعنی (plot(dataframe. در صورتی که مجموعه داده بیش از دو ستون داشته باشد،‌ در محیط گرافیکی، نمودار پراکنش دو به دوی این ستون‌ها و به صورت ماتریسی نمایش داده خواهد شد.

به طور مثال اگر:

(x<-rnorm(1000,0,1

(y<-rnorm(1000,0,1

(z<-rnorm(1000,0,1

(dfr1<-data.frame(x,y

(dfr2<-data.frame(x,y,z

در این صورت اگر بخواهیم نمودار پراکنش x و y یا نمودار پراکنش مجموعه داده dfr1 یا مجموعه نمودار پراکنش تمامی ستون‌های مجموعه داده dfr2 را مشاهده کنیم به ترتیب باید (plot(x,y)، plot(dfr1 و (plot(drf2 را اجرا کنیم.