{
#Suppose you have your data in X and Y
#Y being the vector dependent variables and X being the matrix containing the
#independent variables.

#Y=scan("UBN_precip.txt")
#X=matrix(scan("UBN_predictors.txt"),ncol=7,byrow=T)

# Linear regression of a single independent variable X..

N=length(Y)	# of data points..

# to select the best subset or 'best' model, fit models with all possible
# combinations, compute a quantity PRESS or Cp and, select the model
# with the minimum PRESS

library(MPV)
library(leaps)

combs = leaps(X,Y, nbest=N)     #so that it computes all combinations..
combos = combs$which

ncombos = length(combos[,1])

xpress=1:ncombos
xmse = 1:ncombos

for(i in 1:ncombos){

zz=glm(Y ~ X[,combos[i,]])

xpress[i]=PRESS(zz)

xmse[i] = sum((zz$res)^2) / (N - length(zz$coef))

}

# best model based on MSE
zc=order(xmse)[1]   #best model for PRESS
bestmodelm = lsfit(X[,combos[zc,]], Y)

#best model based on PRESS
zc=order(xpress)[1]   #best model for PRESS
bestmodelp = lsfit(X[,combos[zc,]], Y)

# using AIC
X1 = as.data.frame(X)
zm = glm(Y ~ ., data=X1)
zms = step(zm)