> su.mod1<- lm(lnY ~ X1+X2+X3+X4)
> 
> summary(su.mod1)

Call:
lm(formula = lnY ~ X1 + X2 + X3 + X4)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.43500 -0.17591 -0.02091  0.18400  0.56192 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3.851948   0.266258  14.467  < 2e-16 ***
X1          0.083684   0.028833   2.902  0.00554 ** 
X2          0.012665   0.002315   5.471 1.51e-06 ***
X3          0.015632   0.002100   7.443 1.37e-09 ***
X4          0.032161   0.051465   0.625  0.53493    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.2509 on 49 degrees of freedom
Multiple R-squared: 0.7592,     Adjusted R-squared: 0.7396 
F-statistic: 38.62 on 4 and 49 DF,  p-value: 1.388e-14 

> 
> su.mod2 <- step(su.mod1)
Start:  AIC=-144.59
lnY ~ X1 + X2 + X3 + X4

       Df Sum of Sq    RSS     AIC
- X4    1    0.0246 3.1085 -146.16
<none>              3.0840 -144.59
- X1    1    0.5302 3.6141 -138.02
- X2    1    1.8839 4.9678 -120.84
- X3    1    3.4862 6.5702 -105.75

Step:  AIC=-146.16
lnY ~ X1 + X2 + X3

       Df Sum of Sq    RSS      AIC
<none>              3.1085 -146.161
- X1    1    1.2040 4.3125 -130.483
- X2    1    2.6724 5.7810 -114.658
- X3    1    6.3341 9.4427  -88.162
> 
> summary(su.mod2)

Call:
lm(formula = lnY ~ X1 + X2 + X3)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46992 -0.17954 -0.03071  0.17883  0.59253 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3.766176   0.226758  16.609  < 2e-16 ***
X1          0.095458   0.021692   4.401 5.67e-05 ***
X2          0.013340   0.002035   6.556 2.97e-08 ***
X3          0.016452   0.001630  10.094 1.17e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.2493 on 50 degrees of freedom
Multiple R-squared: 0.7573,     Adjusted R-squared: 0.7427 
F-statistic:    52 on 3 and 50 DF,  p-value: 2.124e-15 

> 
> install.packages("leaps")     # Must have "set mirror" in R
Warning in install.packages("leaps") :
  'lib = "\\ad.ufl.edu/clas/home/stat/Documents/R/R-2.15.0/library"' is not writable
Error in install.packages("leaps") : unable to install packages
> library(leaps)
> 
> allpossreg <- regsubsets(lnY ~ X1+X2+X3+X4,nbest=6,data=surgunit1)
> aprout <- summary(allpossreg)
> 
> with(aprout,round(cbind(which,rsq,adjr2,cp,bic),3))
  (Intercept) X1 X2 X3 X4   rsq adjr2      cp     bic
1           1  0  0  1  0 0.428 0.417  66.489 -22.146
1           1  0  0  0  1 0.422 0.410  67.715 -21.581
1           1  0  1  0  0 0.221 0.206 108.556  -5.498
1           1  1  0  0  0 0.061 0.043 141.164   4.602
2           1  0  1  1  0 0.663 0.650  20.520 -46.814
2           1  0  0  1  1 0.599 0.584  33.504 -37.443
2           1  1  0  1  0 0.549 0.531  43.852 -30.989
2           1  0  1  0  1 0.483 0.463  57.215 -23.654
2           1  1  0  0  1 0.430 0.408  67.972 -18.397
2           1  1  1  0  0 0.263 0.234 102.031  -4.493
3           1  1  1  1  0 0.757 0.743   3.391 -60.502
3           1  0  1  1  1 0.718 0.701  11.424 -52.365
3           1  1  0  1  1 0.612 0.589  32.932 -35.186
3           1  1  1  0  1 0.487 0.456  58.392 -20.089
4           1  1  1  1  1 0.759 0.740   5.000 -56.942
> 
> 
> # These Stepwise Methods are based on Model Criteria, not individual regression coefficients
> # direction="both" begins like backward and works down
> # Criteria:  k=2 uses AIC (default)   k=log(length(y))  uses BIC 
> 
> reg.full <- lm(lnY ~ X1+X2+X3+X4)
> reg.null <- lm(lnY ~ 1)
> 
> summary(reg.full)

Call:
lm(formula = lnY ~ X1 + X2 + X3 + X4)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.43500 -0.17591 -0.02091  0.18400  0.56192 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3.851948   0.266258  14.467  < 2e-16 ***
X1          0.083684   0.028833   2.902  0.00554 ** 
X2          0.012665   0.002315   5.471 1.51e-06 ***
X3          0.015632   0.002100   7.443 1.37e-09 ***
X4          0.032161   0.051465   0.625  0.53493    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.2509 on 49 degrees of freedom
Multiple R-squared: 0.7592,     Adjusted R-squared: 0.7396 
F-statistic: 38.62 on 4 and 49 DF,  p-value: 1.388e-14 

> 
> forward.reg <- step(reg.null,direction="forward",scope=list(upper=reg.full,lower=reg.null))
Start:  AIC=-75.7
lnY ~ 1

       Df Sum of Sq     RSS      AIC
+ X3    1    5.4762  7.3316 -103.827
+ X4    1    5.3990  7.4087 -103.262
+ X2    1    2.8285  9.9792  -87.178
+ X1    1    0.7763 12.0315  -77.079
<none>              12.8077  -75.703

Step:  AIC=-103.83
lnY ~ X3

       Df Sum of Sq    RSS     AIC
+ X2    1    3.0191 4.3125 -130.48
+ X4    1    2.2019 5.1297 -121.11
+ X1    1    1.5506 5.7810 -114.66
<none>              7.3316 -103.83

Step:  AIC=-130.48
lnY ~ X3 + X2

       Df Sum of Sq    RSS     AIC
+ X1    1   1.20395 3.1085 -146.16
+ X4    1   0.69836 3.6141 -138.02
<none>              4.3125 -130.48

Step:  AIC=-146.16
lnY ~ X3 + X2 + X1

       Df Sum of Sq    RSS     AIC
<none>              3.1085 -146.16
+ X4    1  0.024578 3.0840 -144.59
> 
> summary(forward.reg)

Call:
lm(formula = lnY ~ X3 + X2 + X1)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46992 -0.17954 -0.03071  0.17883  0.59253 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3.766176   0.226758  16.609  < 2e-16 ***
X3          0.016452   0.001630  10.094 1.17e-13 ***
X2          0.013340   0.002035   6.556 2.97e-08 ***
X1          0.095458   0.021692   4.401 5.67e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.2493 on 50 degrees of freedom
Multiple R-squared: 0.7573,     Adjusted R-squared: 0.7427 
F-statistic:    52 on 3 and 50 DF,  p-value: 2.124e-15 

> 
> backward.reg <- step(reg.full,direction="backward",k=log(length(lnY)))
Start:  AIC=-134.64
lnY ~ X1 + X2 + X3 + X4

       Df Sum of Sq    RSS      AIC
- X4    1    0.0246 3.1085 -138.205
<none>              3.0840 -134.645
- X1    1    0.5302 3.6141 -130.067
- X2    1    1.8839 4.9678 -112.888
- X3    1    3.4862 6.5702  -97.792

Step:  AIC=-138.2
lnY ~ X1 + X2 + X3

       Df Sum of Sq    RSS      AIC
<none>              3.1085 -138.205
- X1    1    1.2040 4.3125 -124.516
- X2    1    2.6724 5.7810 -108.691
- X3    1    6.3341 9.4427  -82.195
> 
> summary(backward.reg)

Call:
lm(formula = lnY ~ X1 + X2 + X3)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46992 -0.17954 -0.03071  0.17883  0.59253 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3.766176   0.226758  16.609  < 2e-16 ***
X1          0.095458   0.021692   4.401 5.67e-05 ***
X2          0.013340   0.002035   6.556 2.97e-08 ***
X3          0.016452   0.001630  10.094 1.17e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.2493 on 50 degrees of freedom
Multiple R-squared: 0.7573,     Adjusted R-squared: 0.7427 
F-statistic:    52 on 3 and 50 DF,  p-value: 2.124e-15 

> 
> stepwise.reg <- step(reg.full,direction="both")
Start:  AIC=-144.59
lnY ~ X1 + X2 + X3 + X4

       Df Sum of Sq    RSS     AIC
- X4    1    0.0246 3.1085 -146.16
<none>              3.0840 -144.59
- X1    1    0.5302 3.6141 -138.02
- X2    1    1.8839 4.9678 -120.84
- X3    1    3.4862 6.5702 -105.75

Step:  AIC=-146.16
lnY ~ X1 + X2 + X3

       Df Sum of Sq    RSS      AIC
<none>              3.1085 -146.161
+ X4    1    0.0246 3.0840 -144.590
- X1    1    1.2040 4.3125 -130.483
- X2    1    2.6724 5.7810 -114.658
- X3    1    6.3341 9.4427  -88.162
> 
> summary(stepwise.reg)

Call:
lm(formula = lnY ~ X1 + X2 + X3)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.46992 -0.17954 -0.03071  0.17883  0.59253 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 3.766176   0.226758  16.609  < 2e-16 ***
X1          0.095458   0.021692   4.401 5.67e-05 ***
X2          0.013340   0.002035   6.556 2.97e-08 ***
X3          0.016452   0.001630  10.094 1.17e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.2493 on 50 degrees of freedom
Multiple R-squared: 0.7573,     Adjusted R-squared: 0.7427 
F-statistic:    52 on 3 and 50 DF,  p-value: 2.124e-15