set.seed(798)
n<-200
p<-3
Sig<-rWishart(1,p+2,diag(1,p)/(p+2))[,,1]
library(MASS)
x<-mvrnorm(n,mu=rep(0,p),Sigma=Sig)

be<-c(2,2,2,0,-2)

xm<-cbind(rep(1,n),x,x[,1]*x[,2])
Y<-xm %*% be + rnorm(n)

x1<-x[,1]
x2<-x[,2]
x3<-x[,3]

fit0<-lm(Y~1)
fit1<-lm(Y~x1)
fit2<-lm(Y~x2)
fit3<-lm(Y~x3)
fit12<-lm(Y~x1+x2)
fit13<-lm(Y~x1+x3)
fit23<-lm(Y~x2+x3)
fit123<-lm(Y~x1+x2+x3)
fit12i<-lm(Y~x1*x2)
fit13i<-lm(Y~x1*x3)
fit23i<-lm(Y~x2*x3)
fit123i<-lm(Y~x1*x2*x3)

criteria.eval<-function(fit.obj,nv,bigsig.hat){
	cvec<-rep(0,5)
	SSRes<-sum(residuals(fit.obj)^2)
	p<-length(coef(fit.obj))
	cvec[1]<-summary(fit.obj)$r.squared
	cvec[2]<-summary(fit.obj)$adj.r.squared
	cvec[3]<-SSRes/bigsig.hat^2-n+2*p

	#AIC in R computes
	#n*log(sum(residuals(fit.obj)^2)/n)+
      #2*(length(coef(fit.obj))+1)+n*log(2*pi)+n
	cvec[4]<-AIC(fit.obj)

	#BIC in R computes
	#n*log(sum(residuals(fit.obj)^2)/n)+
      #log(n)*(length(coef(fit.obj))+1)+n*log(2*pi)+n
	cvec[5]<-BIC(fit.obj)

	return(cvec)
}

bigs.hat<-summary(fit123i)$sigma



cvals<-matrix(0,nrow=12,ncol=5)
cvals[1,]<-criteria.eval(fit0,n,bigs.hat)
cvals[2,]<-criteria.eval(fit1,n,bigs.hat)
cvals[3,]<-criteria.eval(fit2,n,bigs.hat)
cvals[4,]<-criteria.eval(fit3,n,bigs.hat)
cvals[5,]<-criteria.eval(fit12,n,bigs.hat)
cvals[6,]<-criteria.eval(fit13,n,bigs.hat)
cvals[7,]<-criteria.eval(fit23,n,bigs.hat)
cvals[8,]<-criteria.eval(fit123,n,bigs.hat)
cvals[9,]<-criteria.eval(fit12i,n,bigs.hat)
cvals[10,]<-criteria.eval(fit13i,n,bigs.hat)
cvals[11,]<-criteria.eval(fit23i,n,bigs.hat)
cvals[12,]<-criteria.eval(fit123i,n,bigs.hat)

Criteria<-data.frame(cvals)
names(Criteria)<-c('Rsq','Adj.Rsq','Cp','AIC','BIC')

rownames(Criteria)<-c('1','x1','x2','x3','x1+x2','x1+x3','x2+x3','x1+x2+x3',
                    'x1*x2','x1*x3','x2*x3','x1*x2*x3')