به نام خدا

رگرسیون خطی

فهرست مطالب

کتاب‌خانه‌ها
رگرسیون خطی ساده
رگرسیون خطی چندتایی
روابط متقابل(interaction terms)
تبدیل غیرخطی پیش‌بین‌ها

کتاب‌خانه‌ها

در این بخش از دیتافریم‌های موجود در کتابخانه MASS استفاده می‌کنیم.
دیتافریم Boston : داده‌های جمع‌آوری شده از میانگین قیمت خانه‌های ۵۰۶ محله در اطراف بوستون.
برای دیدن اطلاعات این دیتافریم و توضیح هر یک از متغیر‌ها کافیست مستندات آن را مطالعه کنید.

install.packages("MASS", repos='http://cran.us.r-project.org')

library(MASS)

fix(Boston)
names(Boston)
?Boston

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

برای شروع با استفاده از تابع $\texttt{lm()}$ یک مدل ساده‌ی رگرسیون خطی را با medv به عنوان پاسخ(response) و lstat به عنوان پیش‌بین(predictor) بر داده‌ها منطبق می کنیم. سینتکس این کد به صورت $\texttt{lm(response~predictor, data)}$ است.

lm.fit=lm(medv~lstat ,data=Boston)
attach(Boston)
lm.fit=lm(medv~lstat)

با تایپ کردن lm.fit برخی اطلاعات اساسی در مورد مدل نشان داده می‌شوند.

lm.fit

Call:
lm(formula = medv ~ lstat)

Coefficients:
(Intercept)        lstat  
      34.55        -0.95

برای اطلاعات جزئی‌تر از $\texttt{summary(lm.fit)}$ استفاده ‌می‌کنیم.

summary(lm.fit)

Call:
lm(formula = medv ~ lstat)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.168  -3.990  -1.318   2.034  24.500 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 34.55384    0.56263   61.41   <2e-16 ***
lstat       -0.95005    0.03873  -24.53   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.216 on 504 degrees of freedom
Multiple R-squared:  0.5441,	Adjusted R-squared:  0.5432 
F-statistic: 601.6 on 1 and 504 DF,  p-value: < 2.2e-16

با استفاده از تابع $\texttt{names()}$ می‌توانیم بقیه اطلاعاتی که در lm.fit ذخیره‌ شده‌اند ببینیم.

names(lm.fit)

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

lm.fit$coefficients

اما استفاده از توابع استخراج کننده مانند $\texttt{coef()}$ برای دسترسی به آن‌ها ساده‌تر است:

coef(lm.fit)

برای بدست آوردن یک بازه‌ی اطمینان برای تقریب ضرایب می‌توانیم از دستور زیر استفاده کنیم:

confint(lm.fit)

تابع $\texttt{predict()}$ برای ساختن بازه‌های اطمینان و بازه‌های پیش‌بینی برای پیش‌بینی medv برای یک مقدار lstat ورودی است.

predict(lm.fit,data.frame(lstat=(c(5,10,15))), interval ="confidence")

predict(lm.fit,data.frame(lstat=(c(5,10,15))), interval ="prediction")

حال با استفاده از $\texttt{plot()}$ و $\texttt{abline()}$، نمودار medv و lstat را به همراه رگرسیون کوچکترین مربعات رسم می‌کنیم:

plot(lstat ,medv)
abline(lm.fit)

به نظر می‌رسد رابطه‌ی بین lstat و medv در اینجا غیرخطی است. جلوتر به این مساله خواهیم پرداخت.
با اعمال کردن $\texttt{plot()}$ مستقیما بر روی خروجی $\texttt{lm()}$ چهار نمودارِ تشخیصی رسم می‌شوند. در واقع این دستور یک نمودار را نمایش می‌دهد و با زدن دکمه‌ی enter می‌توان نمودار بعدی را مشاهده کرد. با استفاده از تابع $\texttt{par()}$ می‌توان هر چهار نمودار را با هم دید.

## divides the plotting region into a 2 × 2 grid of panels
par(mfrow=c(2,2)) 
plot(lm.fit)

می‌توانیم انحراف‌‌های برازش یک رگرسیون خطی را با استفاده از تابع $\texttt{residuals()}$ محاسبه کنیم.

plot(predict(lm.fit), residuals(lm.fit))
plot(predict(lm.fit), rstudent(lm.fit))

با مشاهده نمودارهای انحراف می‌توان به غیرخطی بودن رابطه پی برد.
می‌توانیم با تابع $\texttt{hatvalues()}$ آماره‌های leverage را برای هر تعداد پیش‌بین محاسبه کنیم. تابع $\texttt{which.max()}$ اندیس بزرگترین عنصر یک بردار را مشخص می‌کند. در اینجا نشان می‌دهد که کدام مشاهده بزرگترین آماره‌ی leverage را داشته است:

plot(hatvalues (lm.fit))
which.max(hatvalues (lm.fit))

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

برای برازش یک مدل با رگرسیون خطی چندتایی با استفاده از کمترین مربعات، مجددا از تابع $\texttt{lm()}$ استفاده می‌کنیم.

## The syntax lm(y∼x1+x2+x3) is used to fit a model with three predictors, x1, x2, and x3.
lm.fit=lm(medv~lstat+age,data=Boston) 
summary(lm.fit)

Call:
lm(formula = medv ~ lstat + age, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.981  -3.978  -1.283   1.968  23.158 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 33.22276    0.73085  45.458  < 2e-16 ***
lstat       -1.03207    0.04819 -21.416  < 2e-16 ***
age          0.03454    0.01223   2.826  0.00491 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.173 on 503 degrees of freedom
Multiple R-squared:  0.5513,	Adjusted R-squared:  0.5495 
F-statistic:   309 on 2 and 503 DF,  p-value: < 2.2e-16

مجموعه داده‌ی Boston شامل ۱۳ متغیر می‌شود. برای سادگی نوشتار، می‌توان از سینکس زیر برای پیاده‌سازی رگرسیون با استفاده از تمامی این متغیرها استفاده کرد:

lm.fit=lm(medv~.,data=Boston) 
summary(lm.fit)

Call:
lm(formula = medv ~ ., data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.595  -2.730  -0.518   1.777  26.199 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
crim        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
zn           4.642e-02  1.373e-02   3.382 0.000778 ***
indus        2.056e-02  6.150e-02   0.334 0.738288    
chas         2.687e+00  8.616e-01   3.118 0.001925 ** 
nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
age          6.922e-04  1.321e-02   0.052 0.958229    
dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
tax         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
black        9.312e-03  2.686e-03   3.467 0.000573 ***
lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406,	Adjusted R-squared:  0.7338 
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

ما می‌توانیم به تک تک اجزای شیء summary با نام‌شان دسترسی پیدا کنیم ( دستور ?summary.lm را تایپ کنید تا موارد ممکن را ببینید).

## gives us R^2
summary(lm.fit)$r.sq

## gives us RSE
summary(lm.fit)$sigma

مثالی از اجرای رگرسیون با همه‌ی متغیرها به غیر از یکی:

lm.fit1=lm(medv~.-age,data=Boston)
summary(lm.fit1)

Call:
lm(formula = medv ~ . - age, data = Boston)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.6054  -2.7313  -0.5188   1.7601  26.2243 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  36.436927   5.080119   7.172 2.72e-12 ***
crim         -0.108006   0.032832  -3.290 0.001075 ** 
zn            0.046334   0.013613   3.404 0.000719 ***
indus         0.020562   0.061433   0.335 0.737989    
chas          2.689026   0.859598   3.128 0.001863 ** 
nox         -17.713540   3.679308  -4.814 1.97e-06 ***
rm            3.814394   0.408480   9.338  < 2e-16 ***
dis          -1.478612   0.190611  -7.757 5.03e-14 ***
rad           0.305786   0.066089   4.627 4.75e-06 ***
tax          -0.012329   0.003755  -3.283 0.001099 ** 
ptratio      -0.952211   0.130294  -7.308 1.10e-12 ***
black         0.009321   0.002678   3.481 0.000544 ***
lstat        -0.523852   0.047625 -10.999  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.74 on 493 degrees of freedom
Multiple R-squared:  0.7406,	Adjusted R-squared:  0.7343 
F-statistic: 117.3 on 12 and 493 DF,  p-value: < 2.2e-16

هم‌چنین از تابع $\texttt{update()}$ می‌توان استفاده کرد:

lm.fit1=update(lm.fit, ~.-age)

روابط متقابل(interaction terms)

دستورِ lstat:black در R یک تقابل بین lstat و black ایجاد می‌کند. دستورِ lstat*age به طور همزمان lstat و age رابطه lstat×age را به عنوان پیش‌بین قرار می‌دهد (در واقع خلاصه شده‌ی عبارت lstat+age+lstat:age است)

summary(lm(medv~lstat*age,data=Boston))

Call:
lm(formula = medv ~ lstat * age, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.806  -4.045  -1.333   2.085  27.552 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 36.0885359  1.4698355  24.553  < 2e-16 ***
lstat       -1.3921168  0.1674555  -8.313 8.78e-16 ***
age         -0.0007209  0.0198792  -0.036   0.9711    
lstat:age    0.0041560  0.0018518   2.244   0.0252 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.149 on 502 degrees of freedom
Multiple R-squared:  0.5557,	Adjusted R-squared:  0.5531 
F-statistic: 209.3 on 3 and 502 DF,  p-value: < 2.2e-16

تبدیل غیرخطی پیش‌بین‌ها

تابع $\texttt{lm()}$ می‌تواند با تبدیل‌های غیرخطی پیش‌بین‌ها استفاده شود. به طور مثال ما می‌توانیم با داشتن یک پیش‌بینِ X، یک پیش‌بینِ X² را با استفاده از $\text{I(X^2)}$ ساخت. به طور مثال:

lm.fit2=lm(medv~lstat+I(lstat^2))
summary(lm.fit2)

Call:
lm(formula = medv ~ lstat + I(lstat^2))

Residuals:
     Min       1Q   Median       3Q      Max 
-15.2834  -3.8313  -0.5295   2.3095  25.4148 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 42.862007   0.872084   49.15   <2e-16 ***
lstat       -2.332821   0.123803  -18.84   <2e-16 ***
I(lstat^2)   0.043547   0.003745   11.63   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.524 on 503 degrees of freedom
Multiple R-squared:  0.6407,	Adjusted R-squared:  0.6393 
F-statistic: 448.5 on 2 and 503 DF,  p-value: < 2.2e-16

یک p-value غیرصفر مرتبط با یک عبارت درجه دوم، نشان می‌دهد که مدل بهتر شده است. برای نشان دادن برتری بردازش درجه دوم به بردازش خطی به بیانی دیگر، از تابع $\texttt{anova()}$ استفاده می‌کنیم.

lm.fit=lm(medv~lstat)
anova(lm.fit ,lm.fit2)

تابع $\texttt{anova()}$ در این جا یک آزمون فرض اجرا می‌کند که دو مدل را مقایسه کند. فرض صفر این است که دو مدل به یک اندازه بر داده منطبق می‌شوند و فرض دیگر این است که مدل دوم بهتر است.
در اینجا F-statistic، عدد ۱۳۵ است و p-value مربوطه تقریبا صفر است. این به روشنی نشان می‌دهد که مدلی که شامل پیش‌بین‌های ‌‌lstat و lstat² می‌شود به مراتب برتر از مدلی‌‌ست که فقط پیش‌بینِ lstat را دارد. همان‌طور که پیش‌تر هم شواهدی از رابطه‌ای غیرخطی بین medv و lstat دیده بودیم.

با استفاده از تابع $\texttt{poly()}$ می‌توان چندجمله‌ای‌هایی با درجات بالاتر را نیز به سادگی در $\texttt{lm()}$ ساخت. به طور مثال:

## produces a fifth-order polynomial fit
lm.fit5=lm(medv~poly(lstat ,5))
summary(lm.fit5)

Call:
lm(formula = medv ~ poly(lstat, 5))

Residuals:
     Min       1Q   Median       3Q      Max 
-13.5433  -3.1039  -0.7052   2.0844  27.1153 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       22.5328     0.2318  97.197  < 2e-16 ***
poly(lstat, 5)1 -152.4595     5.2148 -29.236  < 2e-16 ***
poly(lstat, 5)2   64.2272     5.2148  12.316  < 2e-16 ***
poly(lstat, 5)3  -27.0511     5.2148  -5.187 3.10e-07 ***
poly(lstat, 5)4   25.4517     5.2148   4.881 1.42e-06 ***
poly(lstat, 5)5  -19.2524     5.2148  -3.692 0.000247 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.215 on 500 degrees of freedom
Multiple R-squared:  0.6817,	Adjusted R-squared:  0.6785 
F-statistic: 214.2 on 5 and 500 DF,  p-value: < 2.2e-16

مثالی از تبدیل غیرچندجمله‌ای:

summary(lm(medv~log(rm),data=Boston))

Call:
lm(formula = medv ~ log(rm), data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-19.487  -2.875  -0.104   2.837  39.816 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -76.488      5.028  -15.21   <2e-16 ***
log(rm)       54.055      2.739   19.73   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.915 on 504 degrees of freedom
Multiple R-squared:  0.4358,	Adjusted R-squared:  0.4347 
F-statistic: 389.3 on 1 and 504 DF,  p-value: < 2.2e-16

	2.5 %	97.5 %
(Intercept)	33.448457	35.6592247
lstat	-1.026148	-0.8739505

fit	lwr	upr
29.80359	29.00741	30.59978
25.05335	24.47413	25.63256
20.30310	19.73159	20.87461

fit	lwr	upr
29.80359	17.565675	42.04151
25.05335	12.827626	37.27907
20.30310	8.077742	32.52846

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
504	19472.38	NA	NA	NA	NA
503	15347.24	1	4125.138	135.1998	7.630116e-28