collinearity.pca <- function(x, nlin=dim(x)[2], lin.combs=F,
                                 digits=3, standardize=F) {
      if(!is.data.frame(x) & !is.matrix(x)) {
        cat("\nError in collinearity.pca(): first (x) argument must be a data frame or matrix.\n")
        return(NULL)
      }
      xx <- as.matrix(x)
      vars.unique <- unlist(lapply(apply(xx, 2, unique), length))
      if(any(vars.unique==1)) {
        cat("The following variables take on only one value and will not be used:\n")
        print(names(vars.unique)[vars.unique==1])
        xx <- xx[,names(vars.unique)[vars.unique!=1]]
      }
      if(standardize) xx <- apply(xx, 2, scale)
      p <- dim(xx)[2];  nlin <- min(p,nlin)   # number of predictors
      eig <- eigen(var(xx))   # eigendecomposition of the covariance matrix
      labels <- paste("PC",1:p,sep="") # eig is badly labeled, correct this:
      names(eig$values) <- labels           # relabel the variances
      dimnames(eig$vectors)[[2]] <- labels  # relabel the lin. combs.
      dimnames(eig$vectors)[[1]] <- colnames(x)  # relabel the lin. combs.
      ord <- (p:1)[1:nlin]    # pick the last nlin in reverse order
      eig$values  <- round(eig$values[ord], digits)   # round for readability
      eig$vectors <- round(eig$vectors[,ord], digits) # dito
      names(eig) <- c("variances","coefficients")     # e'vals=variances
      if(lin.combs) eig$lin.combs <- xx %*% eig$coefficients
      return(eig)
    }