Linear regression
Thanh Bui
6/18/2021
#Loading packages
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.3 v purrr 0.3.4
## v tibble 3.1.1 v dplyr 1.0.5
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(leaps) #Mallow's Cp
library(BMA) #Bayesian Model## Loading required package: survival## Loading required package: robustbase##
## Attaching package: 'robustbase'## The following object is masked from 'package:survival':
##
## heart## Loading required package: inline## Loading required package: rrcov## Scalable Robust Estimators with High Breakdown Point (version 1.5-5)#Backward regression ##x~.= “x is explained by everything else”
FitALL=lm(mpg~., data=mtcars)
step(FitALL,direction = "backward",test="F", k=log(nrow(mtcars)))## Start: AIC=87.02
## mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## - cyl 1 0.0799 147.57 83.572 0.0114 0.91609
## - vs 1 0.1601 147.66 83.590 0.0228 0.88142
## - carb 1 0.4067 147.90 83.643 0.0579 0.81218
## - gear 1 1.3531 148.85 83.847 0.1926 0.66521
## - drat 1 1.6270 149.12 83.906 0.2317 0.63528
## - disp 1 3.9167 151.41 84.394 0.5576 0.46349
## - hp 1 6.8399 154.33 85.006 0.9739 0.33496
## - qsec 1 8.8641 156.36 85.423 1.2621 0.27394
## - am 1 10.5467 158.04 85.765 1.5016 0.23399
## <none> 147.49 87.021
## - wt 1 27.0144 174.51 88.937 3.8463 0.06325 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=83.57
## mpg ~ disp + hp + drat + wt + qsec + vs + am + gear + carb
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## - vs 1 0.2685 147.84 80.165 0.0400 0.84326
## - carb 1 0.5201 148.09 80.219 0.0775 0.78326
## - gear 1 1.8211 149.40 80.499 0.2715 0.60754
## - drat 1 1.9826 149.56 80.534 0.2956 0.59214
## - disp 1 3.9009 151.47 80.942 0.5815 0.45381
## - hp 1 7.3632 154.94 81.665 1.0977 0.30615
## - qsec 1 10.0933 157.67 82.224 1.5047 0.23292
## - am 1 11.8359 159.41 82.575 1.7645 0.19768
## <none> 147.57 83.572
## - wt 1 27.0280 174.60 85.488 4.0293 0.05716 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=80.16
## mpg ~ disp + hp + drat + wt + qsec + am + gear + carb
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## - carb 1 0.6855 148.53 76.847 0.1066 0.74696
## - gear 1 2.1437 149.99 77.160 0.3335 0.56922
## - drat 1 2.2139 150.06 77.175 0.3444 0.56301
## - disp 1 3.6467 151.49 77.479 0.5673 0.45897
## - hp 1 7.1060 154.95 78.201 1.1055 0.30399
## - am 1 11.5694 159.41 79.110 1.7999 0.19283
## - qsec 1 15.6830 163.53 79.925 2.4398 0.13195
## <none> 147.84 80.165
## - wt 1 27.3799 175.22 82.136 4.2595 0.05049 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=76.85
## mpg ~ disp + hp + drat + wt + qsec + am + gear
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## - gear 1 1.565 150.09 73.717 0.2529 0.619641
## - drat 1 1.932 150.46 73.795 0.3122 0.581508
## - disp 1 10.110 158.64 75.489 1.6337 0.213420
## - am 1 12.323 160.85 75.932 1.9913 0.171042
## - hp 1 14.826 163.35 76.426 2.3956 0.134763
## <none> 148.53 76.847
## - qsec 1 26.408 174.94 78.618 4.2672 0.049815 *
## - wt 1 69.127 217.66 85.610 11.1699 0.002717 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=73.72
## mpg ~ disp + hp + drat + wt + qsec + am
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## - drat 1 3.345 153.44 70.956 0.5571 0.462401
## - disp 1 8.545 158.64 72.023 1.4233 0.244054
## - hp 1 13.285 163.38 72.965 2.2127 0.149381
## <none> 150.09 73.717
## - am 1 20.036 170.13 74.261 3.3372 0.079692 .
## - qsec 1 25.574 175.67 75.286 4.2598 0.049551 *
## - wt 1 67.572 217.66 82.146 11.2550 0.002536 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=70.96
## mpg ~ disp + hp + wt + qsec + am
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## - disp 1 6.629 160.07 68.844 1.1232 0.298972
## - hp 1 12.572 166.01 70.011 2.1303 0.156387
## <none> 153.44 70.956
## - qsec 1 26.470 179.91 72.583 4.4853 0.043908 *
## - am 1 32.198 185.63 73.586 5.4559 0.027488 *
## - wt 1 69.043 222.48 79.380 11.6993 0.002075 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=68.84
## mpg ~ hp + wt + qsec + am
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## - hp 1 9.219 169.29 67.170 1.5551 0.223088
## <none> 160.07 68.844
## - qsec 1 20.225 180.29 69.186 3.4115 0.075731 .
## - am 1 25.993 186.06 70.193 4.3845 0.045791 *
## - wt 1 78.494 238.56 78.147 13.2403 0.001141 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=67.17
## mpg ~ wt + qsec + am
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 169.29 67.170
## - am 1 26.178 195.46 68.306 4.3298 0.0467155 *
## - qsec 1 109.034 278.32 79.614 18.0343 0.0002162 ***
## - wt 1 183.347 352.63 87.187 30.3258 6.953e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Call:
## lm(formula = mpg ~ wt + qsec + am, data = mtcars)
##
## Coefficients:
## (Intercept) wt qsec am
## 9.618 -3.917 1.226 2.936#Forward regression
FitStart=lm(mpg~1, data = mtcars)
step(FitStart,mpg~ cyl+disp+hp+drat+wt+qsec+vs+am+gear+carb ,direction="forward",test="F")## Start: AIC=115.94
## mpg ~ 1
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + wt 1 847.73 278.32 73.217 91.3753 1.294e-10 ***
## + cyl 1 817.71 308.33 76.494 79.5610 6.113e-10 ***
## + disp 1 808.89 317.16 77.397 76.5127 9.380e-10 ***
## + hp 1 678.37 447.67 88.427 45.4598 1.788e-07 ***
## + drat 1 522.48 603.57 97.988 25.9696 1.776e-05 ***
## + vs 1 496.53 629.52 99.335 23.6622 3.416e-05 ***
## + am 1 405.15 720.90 103.672 16.8603 0.000285 ***
## + carb 1 341.78 784.27 106.369 13.0736 0.001084 **
## + gear 1 259.75 866.30 109.552 8.9951 0.005401 **
## + qsec 1 197.39 928.66 111.776 6.3767 0.017082 *
## <none> 1126.05 115.943
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=73.22
## mpg ~ wt
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + cyl 1 87.150 191.17 63.198 13.2203 0.001064 **
## + hp 1 83.274 195.05 63.840 12.3813 0.001451 **
## + qsec 1 82.858 195.46 63.908 12.2933 0.001500 **
## + vs 1 54.228 224.09 68.283 7.0177 0.012926 *
## + carb 1 44.602 233.72 69.628 5.5343 0.025646 *
## + disp 1 31.639 246.68 71.356 3.7195 0.063620 .
## <none> 278.32 73.217
## + drat 1 9.081 269.24 74.156 0.9781 0.330854
## + gear 1 1.137 277.19 75.086 0.1189 0.732668
## + am 1 0.002 278.32 75.217 0.0002 0.987915
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Step: AIC=63.2
## mpg ~ wt + cyl
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## + hp 1 14.5514 176.62 62.665 2.3069 0.1400
## + carb 1 13.7724 177.40 62.805 2.1738 0.1515
## <none> 191.17 63.198
## + qsec 1 10.5674 180.60 63.378 1.6383 0.2111
## + gear 1 3.0281 188.14 64.687 0.4507 0.5075
## + disp 1 2.6796 188.49 64.746 0.3980 0.5332
## + vs 1 0.7059 190.47 65.080 0.1038 0.7497
## + am 1 0.1249 191.05 65.177 0.0183 0.8933
## + drat 1 0.0010 191.17 65.198 0.0001 0.9903
##
## Step: AIC=62.66
## mpg ~ wt + cyl + hp
##
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 176.62 62.665
## + am 1 6.6228 170.00 63.442 1.0519 0.3142
## + disp 1 6.1762 170.44 63.526 0.9784 0.3314
## + carb 1 2.5187 174.10 64.205 0.3906 0.5372
## + drat 1 2.2453 174.38 64.255 0.3477 0.5603
## + qsec 1 1.4010 175.22 64.410 0.2159 0.6459
## + gear 1 0.8558 175.76 64.509 0.1315 0.7197
## + vs 1 0.0599 176.56 64.654 0.0092 0.9245##
## Call:
## lm(formula = mpg ~ wt + cyl + hp, data = mtcars)
##
## Coefficients:
## (Intercept) wt cyl hp
## 38.75179 -3.16697 -0.94162 -0.01804#Stepwise regression
step(FitStart,scope=formula(FitALL),direction="both",test="F",trace=0)##
## Call:
## lm(formula = mpg ~ wt + cyl + hp, data = mtcars)
##
## Coefficients:
## (Intercept) wt cyl hp
## 38.75179 -3.16697 -0.94162 -0.01804#Mallow’s Cp
a= regsubsets(mpg~.,nbest=3, data=mtcars)
plot(a, scale="Cp")
#Bayesian Model Average
X=mtcars%>%mutate(mpg=NULL)
Y=mtcars$mpg
search=bicreg(X,Y, strict = FALSE, OR=20)
summary(search)##
## Call:
## bicreg(x = X, y = Y, strict = FALSE, OR = 20)
##
##
## 86 models were selected
## Best 5 models (cumulative posterior probability = 0.2457 ):
##
## p!=0 EV SD model 1 model 2 model 3
## Intercept 100.0 26.3181386 13.324900 9.61778 39.68626 37.22727
## cyl 37.8 -0.4285797 0.685863 . -1.50779 .
## disp 12.5 0.0004793 0.005408 . . .
## hp 41.1 -0.0104394 0.015618 . . -0.03177
## drat 15.4 0.1988568 0.739629 . . .
## wt 98.8 -3.5386668 1.121556 -3.91650 -3.19097 -3.87783
## qsec 48.8 0.4506877 0.568549 1.22589 . .
## vs 10.6 0.0878027 0.645424 . . .
## am 40.1 1.0835917 1.688212 2.93584 . .
## gear 13.1 0.0964727 0.514582 . . .
## carb 22.9 -0.1294021 0.342176 . . .
##
## nVar 3 2 2
## r2 0.850 0.830 0.827
## BIC -50.23818 -49.81447 -49.17255
## post prob 0.071 0.057 0.042
## model 4 model 5
## Intercept 19.74622 38.75179
## cyl . -0.94162
## disp . .
## hp . -0.01804
## drat . .
## wt -5.04798 -3.16697
## qsec 0.92920 .
## vs . .
## am . .
## gear . .
## carb . .
##
## nVar 2 3
## r2 0.826 0.843
## BIC -49.10426 -48.88168
## post prob 0.040 0.036imageplot.bma(search)