Linear Regression in R

Let’s build artificial data:

> x = 1:100
> y = 3*x + rnorm(100, 0, 0.1)

Now, y has 3*x + error (which follows independent normal distribution.)

> y
  [1]   3.084463   5.921660   9.043298  11.863854  14.975701  17.980051  20.812897  23.950546  27.013619  29.879748  33.029797
 [12]  35.933785  39.032856  42.021709  45.150473  47.988470  51.005696  53.921992  56.956936  60.154561  62.870397  66.056555
 [23]  68.818562  71.878614  75.037485  77.971538  81.212420  84.075885  86.924858  89.879004  92.987709  96.250259  98.866730
 [34] 101.816094 104.976485 108.138008 111.180993 114.016497 117.019846 119.892135 123.005147 125.958263 129.165459 131.892357
 [45] 134.985222 137.992608 141.158806 144.089719 146.924594 149.944683 153.130148 156.113614 158.992470 162.068371 165.093236
 [56] 168.134131 171.038950 174.138382 176.801168 180.132722 183.137894 185.884679 188.922193 192.098572 195.171722 198.037255
 [67] 201.071104 204.023956 206.970078 209.960827 212.839931 215.842587 218.939796 221.790648 225.065676 228.100543 231.090998
 [78] 233.996619 236.942362 239.952409 242.805987 246.144036 248.835554 251.792460 255.217650 257.822964 260.989347 264.009189
 [89] 267.065940 270.070240 272.876032 275.924888 279.106800 282.153525 284.942045 288.031658 290.881470 293.930968 296.895243
[100] 299.968869

Build a model:

> m = lm(y ~ x)
> m

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x
 -1.168e-05    3.000e+00

This shows that yhat = 3.0 * x – 1.168e-05. Let’s check the goodness of this model:

> summary(m)

Call:
lm(formula = y ~ x)

Residuals:
      Min        1Q    Median        3Q       Max
-0.204330 -0.074342 -0.004509  0.079480  0.252437

Coefficients:
              Estimate Std. Error  t value Pr(>|t|)
(Intercept) -1.168e-05  2.239e-02   -0.001        1
x            3.000e+00  3.849e-04 7793.923   —
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1111 on 98 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1
F-statistic: 6.075e+07 on 1 and 98 DF,  p-value: < 2.2e-16

Adjusted R-squared says that 100% of variation is explained by this linear regression, and we can see that t value of x has p < 0.05 meaning that x is statistically meaningful in this model. It’s easy to draw some diagrams:

> par(mfrow=c(2,2))
> plot(m)


Residual VS Fitted looks fine (though it would have been much nicer if the line were totally flat). Normal Q-Q plot shows us that residual follows a normal distribution.

Let’s make an intentional mistake. In the following we will try to model y = ax + bz + c while y is actually 3x + error:

> z = runif(100, 0, 50)
> m2 = lm (y ~ x + z)
> m2

Call:
lm(formula = y ~ x + z)

Coefficients:
(Intercept)            x            z
  0.0052540    2.9999384   -0.0002374

> summary(m2)

Call:
lm(formula = y ~ x + z)

Residuals:
      Min        1Q    Median        3Q       Max
-0.199865 -0.075759 -0.007949  0.077534  0.256803

Coefficients:
              Estimate Std. Error  t value Pr(>|t|)
(Intercept)  0.0052540  0.0284124    0.185    0.854
x            2.9999384  0.0003872 7747.345   z           -0.0002374  0.0007824   -0.303    0.762
—
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1116 on 97 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1
F-statistic: 3.009e+07 on 2 and 97 DF,  p-value: < 2.2e-16  

As one can see, z’s p value is 0.762 which is larger than 0.05 meaning that the value is not useful in modeling. Let’s compare m and m2 using ANOVA:

> anova(m, m2)
Analysis of Variance Table

Model 1: y ~ x
Model 2: y ~ x + z
  Res.Df    RSS Df Sum of Sq     F Pr(>F)
1     98 1.2098
2     97 1.2087  1 0.0011466 0.092 0.7623
This shows that the model 2 is not too different from model 1 (as p value 0.7623 is larger than 0.05, they’re not statistically different model.) Therefore, it’s fine to think that z is meaningless. (Recall that p value of z shows that z is not meaningful.)
We may find y ~ x starting from y ~ x + z automatically:
> step(m2, scope=list(lower=y ~ 1, direction=”backward”))
Start:  AIC=-435.57
y ~ x + z

       Df Sum of Sq    RSS     AIC
- z     1         0      1 -437.47
                   1 -435.57
- x     1    747884 747885  895.98

Step:  AIC=-437.47
y ~ x

       Df Sum of Sq    RSS     AIC
                   1 -437.47
- x     1    749891 749892  894.25

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x
 -1.168e-05    3.000e+00

Final model from step was y ~ x. When we removed z, AIC became -437.47 while AIC of y ~ x + z was -435.57 (smaller AIC is better).

Let’s perform boxcox transformation. For that, let’s change x to sqrt(x):

> x = 1:100
> y = 3 * x + rnorm(100, 0, 0.1)
> x = sqrt(x)
> m3 = lm(y ~ x)
> summary(m3)

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max
-17.535 -14.181  -4.509  11.042  60.449

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -93.8145     5.1147  -18.34   x            36.5330     0.7197   50.76   —
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.75 on 98 degrees of freedom
Multiple R-squared: 0.9634,	Adjusted R-squared: 0.963
F-statistic:  2576 on 1 and 98 DF,  p-value: < 2.2e-16  

This model looks good, but boxcox transformation may make it better:

> library(MASS)
> boxcox(m3)
> bc = boxcox(m3)

> bc$x[which.max(bc$y)]
[1] 0.5050505

bc$x and bc$y are points in the above plot. Thus, we can see that lambda = 0.5050505 has max. likelihood. To apply that:

> m4 = lm(I(y^0.5050505) ~ x)
> summary(m4)

Call:
lm(formula = I(y^0.5050505) ~ x)

Residuals:
      Min        1Q    Median        3Q       Max
-0.020409 -0.008784 -0.001066  0.006676  0.091507

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0916328  0.0044541  -20.57   x            1.7903556  0.0006268 2856.43   —
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01458 on 98 degrees of freedom
Multiple R-squared:     1,	Adjusted R-squared:     1
F-statistic: 8.159e+06 on 1 and 98 DF,  p-value: < 2.2e-16  

Now we have much better Adjusted R-square. Or we may do forward-backward selection:

> m5 = lm(y ~ 1)
> step(m5, list(lower=y~1, upper=y~x+z), direction=”both”)
Start:  AIC=894.27
y ~ 1

       Df Sum of Sq    RSS    AIC
+ x     1    722550  27484 565.62
              750034 894.27
+ z     1      2015 748019 896.00

Step:  AIC=565.62
y ~ x

       Df Sum of Sq    RSS    AIC
               27484 565.62
+ z     1        45  27439 567.46
- x     1    722550 750034 894.27

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x
     -93.81        36.53

But this is not a silver bullet; it’s better to find better feature in the first place.
References)
1. R Cookbook by Paul Teetor. Copyright 2011 Paul Teetor, 978-0-596-80915-7.
2. 구자용, 박현진, 최대우, 김성수, 데이터마이닝, KNOU Press.

Similar Posts:

Post a Comment

Your email is never published nor shared. Required fields are marked *