library(leaps)

regress.table <- function(covariates,y){
	n <- length(y)
	p <- dim(covariates)[2]
	num.models <- ((2^p)-1)
results <- regsubsets(covariates,y, nbest=num.models)
coeffs <- coef(results,id=seq(1,num.models))
coeff.matrix <- matrix(NA,num.models,(p+1))
dimnames(coeff.matrix) <- list(seq(1,num.models),c("(Intercept)",dimnames(covariates)[[2]]))
for(i in 1:length(coeffs)){
	for(j in 1:(length(coeffs[[i]]))){
		tmp <- coeffs[[i]]
		variable <- (names(tmp))[j]
	coeff.matrix[i,variable] <- tmp[j]
	}
}
return(coeff.matrix)
}


n <- 2000

# Simple example
x1 <- 2 + rnorm(n,0,1)
x2 <- 2.5 + 2*x1 + rnorm(n,0,1)
x3 <- 3.6 -1.2*x1 + rnorm(n,0,1)
x4 <- 1.5 -2*x2 + 2.5*x3 + rnorm(n,0,1)
covariates <- cbind(x1,x2,x3)

table1 <- regress.table(covariates,x4)


# Not so simple example
x1 <- rnorm(n,0,2)
x2 <- rnorm(n,0,2)
x3 <- 1 + 2*x1 + 1.5*x2 + rnorm(n,0,1)
x4 <- -1 + -2*x1 + rnorm(n,0,1)
x5 <- 2 + 3*x2 -1*x4 + rnorm(n,0,1)

covariates <- cbind(x1,x2,x3,x4)
table2 <- regress.table(covariates,x5)
round(table2,2)



summary(results)