library(leaps) library(MASS) library(lmtest) y=c(2.30, 2.00, 2.31, 2.00, 2.71, 1.90, 1.90, 2.10, 2.31, 2.31, 2.52, 1.81, 2.91, 2.52, 2.23, 2.34, 1.94, 1.53, 2.33, 2.24, 2.04, 2.13, 1.85, 2.34, 2.44, 2.16, 2.26, 2.76, 1.86, 2.25, 1.85, 1.76, 2.06, 2.47, 2.06, 2.26, 2.07, 2.08, 2.18, 2.17, 2.00, 1.88, 2.68, 2.18, 2.28, 2.09, 2.49, 2.60, 2.20, 2.49, 2.09, 2.10, 2.30, 2.50) x1=c(6.7, 5.1, 7.4, 6.5, 7.8, 5.8, 5.7, 3.7, 6.0, 3.7, 6.3, 6.7, 5.8, 5.8, 7.7, 7.4, 6.0, 3.7, 7.3, 5.6, 5.2, 3.4, 6.7, 5.8, 6.3, 5.8, 5.2, 11.2, 5.2, 5.8, 3.2, 8.7, 5.0, 5.8, 5.4, 5.3, 2.6, 4.3, 4.8, 5.4, 5.2, 3.6, 8.8, 6.5, 3.4, 6.5, 4.5, 4.8, 5.1, 3.9, 6.6, 6.4, 6.4, 8.8) x2=c(62, 59, 57, 73, 65, 38, 46, 68, 67, 76, 84, 51, 96, 83, 62, 74, 85, 51, 68, 57, 52, 83, 26, 67, 59, 61, 52, 76, 54, 76, 64, 45, 59, 72, 58, 51, 74, 8, 61, 52, 49, 28, 86, 56, 77, 40, 73, 86, 67, 82, 77, 85, 59, 78) x3=c(81, 66, 83, 41, 115, 72, 63, 81, 93, 94, 83, 43, 114, 88, 67, 68, 28, 41, 74, 87, 76, 53, 68, 86, 100, 73, 86, 90, 56, 59, 65, 23, 73, 93, 70, 99, 86, 119, 76, 88, 72, 99, 88, 77, 93, 84, 106, 101, 77, 103, 46, 40, 85, 72) x4=c(2.6, 1.7, 2.2, 2.0, 4.3, 1.4, 1.9, 2.6, 2.5, 2.4, 4.1, 1.9, 4.0, 4.0, 3.4, 2.4, 3.0, 1.6, 3.6, 3.0, 2.9, 1.1, 2.1, 3.4, 3.0, 3.5, 2.5, 5.6, 2.7, 2.6, 0.7, 2.5, 3.5, 3.3, 2.6, 2.6, 2.1, 2.9, 2.5, 1.8, 1.8, 1.3, 6.4, 2.9, 1.5, 3.0, 3.1, 4.1, 2.9, 4.6, 2.0, 1.1, 2.3, 3.2) x=cbind(x1, x2, x3, x4) mydata=data.frame(y, x1, x2, x3, x4) pairs(x) rr=cor(x) rr vif=diag(solve(rr)) vif eigv=eigen(rr)$values ci=sqrt(eigv[1]/eigv) ci ls=regsubsets(y ~ x1 + x2 + x3 + x4, data=mydata) rs=summary(ls) rs cp=rs$cp cp adjr2=rs$adjr2 adjr2 fit=lm(y ~ x1 + x2 + x3 + x4, data=mydata) step=stepAIC(fit, direction="both")