# To read your data for example
 library(leaps)
 
################################################################
# This is an example of applying the R square, Adjusted R      #   #    square, PRESS and AIC Criteria to selection model         #
################################################################

 data <- read.table("fromage.txt",header=T)
 
# Selection variables by Cp and R^2 criteria
 X=data[,c(3,4,5,7,8,9,10,11,12,13,14,15)] # Matrixe of the covariate(free) variables
Y<-data$lfines  # The response Variable 

 # How we can get the CP
 cp<-leaps(x=X, y=Y,method="Cp",nbest=2,names=colnames(X))
 # For adjusted R^2, put the option method=adjr2
 cp$which

 
# The valuse of interest are (Cp-p), so we should caculate
 cpmp<-cp$Cp-cp$size  # Cp-p
 cpmp
abscpmp<-abs(cp$Cp-cp$size) # The absolute values of Cp-p
abscpmp
 # Creat a table of models and values of Cp-p
 res.cp<-cbind(cp$which,cpmp,abscpmp)
 res.cp
 # Then ordeh the table based on small to grand of Cp-p
 ordre<-order(abscpmp)
 ordre
 res.cp.c<-res.cp[ordre,] # Orderd by contain the models based on Cp values.
 res.cp.c
 # Calculate the other statistics for 4 selected models
 nbvar<-12   # Important nomber of variables
 lignes<-c(2,3,4,5) # Les chiffres correspondent aux 4 modèles #que l'on souhaite étudier plus en détail
 press<-NULL
 aic<-NULL
 r2a<-NULL
 #Add the PRESS residuals and AIC table
 var1<- as.matrix(X[,as.logical(res.cp.c[lignes[1],1:nbvar])] )
 mod1<-lm(Y~var1)
 summary(mod1)

 # PRESS residuals for 4 selected models
 press[1]<-sum((mod1$residuals/(1-hatvalues(mod1)))^2)
 # AIC of the model
 aic[1]<-AIC(mod1)
 AIC(mod1)
 s.mod1<-summary(mod1)
 r2a[1]<-s.mod1$adj.r.squared
 r2a[1]
 var2<- as.matrix(X[,as.logical(res.cp.c[lignes[2],1:nbvar])] )
 mod2<-lm(Y~var2)
 press[2]<-sum((mod2$residuals/(1-hatvalues(mod2)))^2)
 press[2]
 aic[2]<-AIC(mod2)
 AIC(mod2)
 s.mod2<-summary(mod2)
 summary(mod2)

 r2a[2]<-s.mod2$adj.r.squared
 r2a[2]
 var3<- as.matrix(X[,as.logical(res.cp.c[lignes[3],1:nbvar])] )
 mod3<-lm(Y~var3)
 press[3]<-sum((mod3$residuals/(1-hatvalues(mod3)))^2)
 press[3]
 aic[3]<-AIC(mod3)
 aic[3]
 s.mod3<-summary(mod1)
 r2a[3]<-s.mod3$adj.r.squared
 r2a[3]
 
 var4<- as.matrix(X[,as.logical(res.cp.c[lignes[4],1:nbvar])] )
 mod4<-lm(Y~var4)
 press[4]<-sum((mod4$residuals/(1-hatvalues(mod4)))^2)
 aic[4]<-AIC(mod4)
 s.mod4<-summary(mod4)
 r2a[4]<-s.mod4$adj.r.squared
 cbind(r2a,press)
 #Table of results selected modeles
 resume<-cbind(res.cp.c[1:4,],press,aic,r2a)
 resume