"mat.regress" <-
function(m,x,y,n.obs=NULL)  {
 #a function to extract subsets of variables (a and b) from a correlation matrix m or data set m
  #and find the multiple correlation beta weights + R2 of the a set predicting the b set
  #seriously rewritten, March 24, 2009 to make much simpler
  #minor additons, October, 20, 2009 to allow for print and summary function
   cl <- match.call()
  if(!is.matrix(m)) m <- as.matrix(m)
  if(dim(m)[1]!=dim(m)[2]) {n.obs=dim(m)[1]
                    C <- cov(m,use="pairwise")
                    m <- cov2cor(C)}  else {
                    C <- m
                    m <- cov2cor(m)}
                   
 
        nm <- dim(m)[1]
        xy <- c(x,y)
        numx <- length(x)
     	numy <- length(y)
        nxy <- numx+numy
     	a.matrix <- m[x,x]
     	ac.matrix <- C[x,x]
     	b.matrix <- m[x,y]
     	bc.matrix <- C[x,y]
     	
        beta <- solve(a.matrix,b.matrix)       #solve the equation bY~aX
        beta <- as.matrix(beta)
     	if (length(y) >1 ) { if(is.null(rownames(beta))) {rownames(beta) <- colnames(m)[x]}
     	                      if(is.null(colnames(beta))) {colnames(beta) <- colnames(m)[y]}
     	 
     	R2 <- colSums(beta * b.matrix) } else { colnames(beta) <- colnames(m)[1]
     	 R2 <- sum(beta * b.matrix)
     	 names(beta) <- colnames(m)[x]
     	 names(R2) <- colnames(m)[y]}
     	if(!is.null(n.obs)) {k<- length(x)
     	                     
     	                     uniq <- (1-smc(a.matrix))
     	                     se.beta <- list() 
     	                     for (i in 1:length(y)) {
     	                     df <- n.obs-k-1
     	                     se.beta[[i]] <- (sqrt((1-R2[i])/(df))*sqrt(1/uniq))}
     	                     se <- matrix(unlist(se.beta),ncol=length(y))
     	                     colnames(se) <- colnames(beta)
     	                     rownames(se) <- rownames(beta)
     	                     tvalue <- beta/se
     	                     se <- t(t(se) * sqrt(diag(C)[y]))/sqrt(diag(ac.matrix))
     	                     
     	                prob <- 2*(1- pt(abs(tvalue),df))
     	                     SE2 <- 4*R2*(1-R2)^2*(df^2)/((n.obs^2-1)*(n.obs+3))
     	                     SE =sqrt(SE2)
     	                     F <- R2*df/(k*(1-R2))
     	                     pF <- 1 - pf(F,k,df)
     	                     shrunkenR2 <- 1-(1-R2)*(n.obs-1)/df    	                     }
     	
     	beta <-  t(t(beta) * sqrt(diag(C)[y]))/sqrt(diag(ac.matrix)) #this puts the betas into the raw units
        
     	if(is.null(n.obs)) {mat.regress <- list(beta=beta,R=sqrt(R2),R2=R2,Call = cl)} else {
     	              mat.regress <- list(beta=beta,se=se,t=tvalue,Probability = prob,R=sqrt(R2),R2=R2,shrunkenR2 = shrunkenR2,seR2 = SE,F=F,probF=pF,df=c(k,df),Call = cl)}
     	class(mat.regress) <- c("psych","mat.regress")
     	return(mat.regress)
     	}
#modified July 12,2007 to allow for NA in the overall matrix
#modified July 9, 2008 to give statistical tests
#modified yet again August 15 , 2008 to convert covariances to correlations