## first import the dataset using the "File" and "Import Data"

attach(Case1202)

## variables are "bsal" "sal77" "fsex" "senior" "age" "educ" "exper" 

## Sleuth uses log transformed salary variable, squared terms 
## interaction terms, so create these

## log transformed salary
log.bsal_log(bsal)

## squared variables
senior2_senior^2
age2_age^2
educ2_educ^2
exper2_exper^2

## interaction terms
senior.i.age_senior*age
senior.i.educ_senior*educ
senior.i.exper_senior*exper
age.i.educ_age*educ
age.i.exper_age*exper
educ.i.exper_educ*exper

out <- select.cp.aic.bic.keep.fun(y=log.bsal,x=cbind(senior,age,educ,exper,senior2,age2,educ2,exper2,senior.i.age,senior.i.educ,senior.i.exper,age.i.educ,age.i.exper,educ.i.exper),keep=1:4)

> names(out)
[1] "model" "p"     "Cp"    "AIC"   "BIC"   "R2"    "R2adj"

The top 25 models:

  > out$model
 [1] senior,age,educ,exper,age.i.educ,age.i.exper                               
 [2] senior,age,educ,exper,exper2,age.i.educ                                    
 [3] senior,age,educ,exper,age.i.exper,educ.i.exper                             
 [4] senior,age,educ,exper,senior.i.educ,age.i.educ,age.i.exper                 
 [5] senior,age,educ,exper,age2,age.i.educ,age.i.exper                          
 [6] senior,age,educ,exper,senior2,age.i.educ,age.i.exper                       
 [7] senior,age,educ,exper,age2,age.i.exper,educ.i.exper                        
 [8] senior,age,educ,exper,exper2,senior.i.educ,age.i.educ                      
 [9] senior,age,educ,exper,age.i.exper                                          
[10] senior,age,educ,exper,educ2,age.i.educ,age.i.exper                         
[11] senior,age,educ,exper,age.i.educ,age.i.exper,educ.i.exper                  
[12] senior,age,educ,exper,exper2,age.i.educ,age.i.exper                        
[13] senior,age,educ,exper,exper2,educ.i.exper                                  
[14] senior,age,educ,exper,senior.i.age,age.i.educ,age.i.exper                  
[15] senior,age,educ,exper,senior.i.exper,age.i.educ,age.i.exper                
[16] senior,age,educ,exper,age2,age.i.educ                                      
[17] senior,age,educ,exper,exper2                                               
[18] senior,age,educ,exper,age2,age.i.exper                                     
[19] senior,age,educ,exper,age.i.educ                                           
[20] senior,age,educ,exper,senior2,senior.i.educ,age.i.educ,age.i.exper         
[21] senior,age,educ,exper,educ2,age.i.exper                                    
[22] senior,age,educ,exper,senior.i.educ,age.i.educ,age.i.exper,educ.i.exper    
[23] senior,age,educ,exper,age2,senior.i.educ,age.i.educ,age.i.exper            
[24] senior,age,educ,exper,exper2,senior.i.educ,age.i.educ,age.i.exper          
[25] senior,age,educ,exper,educ2,senior.i.educ,age.i.educ,age.i.exper  

> range(out$BIC)
[1] -401.402 -361.240
> range(out$Cp)
[1]  1 15
> range(out$AIC)
[1] -419.131 -378.610

postscript(file="comparecrit.ps")
par(mfrow=c(1,3))
plot(p,out$BIC,ylim=c(-425,-360),ylab="BIC or AIC criterion",type="n",sub="A=AIC, B=BIC")
mtext(side=3,line=1,"BIC,AIC vs. p",cex=1.2)
points(p,jitter(out$BIC,factor=2),pch="B", cex=.5)
points(p,jitter(out$AIC,factor=2),pch="A", cex=.5)
text(-360,2,"Note p at min(BIC) and min(AIC)")
plot(p,out$R2,cex=.5,main="R2")
plot(p,out$R2adj,cex=.5,main="R2 adjusted")
dev.off()