SE.comparisons <- 
function(x, y, Nboot=10000, digits=3) {   # x must be a matrix or dataframe, even if it's just a column
  X <- cbind(icept=1,as.matrix(x));  Y <- y;  N <- nrow(X);  p <- ncol(X)
  if(nrow(X) != length(y)) { cat("Number of X rows and length of y are different!\n");  return() }
  ## X <- scale(x);  Y <- scale(y)
  ## mod.lm <- summary(lm(Y ~ . - 1, data=as.data.frame(X)))
  mod.lm <- summary(lm(Y ~ ., data=as.data.frame(X)))
  slopes.bs <- rbind(sapply(1:Nboot, function(i) { sel <- sample(nrow(X), replace=T);  if(i%%1000==0) cat(i,"")
                                                   ## Xb <- scale(X[sel,],scale=F);  Yb <- Y[sel]
                                                   Xb <- X[sel,];  Yb <- Y[sel]
                                                   solve(t(Xb)%*%Xb)%*%t(Xb)%*%Yb } )); cat("\n") # rbind() for matrix
  slopes.bs.SE <- apply(slopes.bs, 1, sd)
  # sandwich estimator, corrected for design:
  X3 <- X%*%solve(t(X)%*%X)
  Hdiag <- diag(X%*%solve(t(X)%*%X)%*%t(X))
  # SE.sandw <- sqrt(apply(resid(mod.lm)^2 * X3^2, 2, sum) * N/(N-p))  # global correction
  SE.sandw <- sqrt(apply(resid(mod.lm)^2 * X3^2 / (1-Hdiag), 2, sum))  # local correction
  # root of RAV sample estimate:
  X.adj <- sapply(1:ncol(X), function(j) resid(lm(X[,j] ~ ., as.data.frame(X[,-j]))))
  X.adj[,1] <- resid(lm(X[,1] ~ .-1, data=as.data.frame(X[,-1])))
  colnames(X.adj) <- colnames(X)
  y.res <- resid(lm(y ~ X))
  root.RAV <- sapply(colnames(X.adj), function(j) sqrt(mean(y.res^2 * X.adj[,j]^2)/mean(y.res^2)/mean(X.adj[,j]^2)) )
  # List:
  round(cbind("Slope"=mod.lm$coeff[,"Estimate"],
              "SE.conv"=mod.lm$coeff[,"Std. Error"],
              "SE.boot"=slopes.bs.SE,
              "SE.sandw"=SE.sandw,
              "boot/conv"=slopes.bs.SE/mod.lm$coeff[,"Std. Error"],
              "sandw/conv"=SE.sandw/mod.lm$coeff[,"Std. Error"],
              "sandw/boot"=SE.sandw/slopes.bs.SE,
              "root.RAV"=root.RAV,
              "t.conv"=mod.lm$coeff[,"t value"],
              "t.boot"=mod.lm$coeff[,"t value"] * mod.lm$coeff[,"Std. Error"] / slopes.bs.SE,
              "t.sand"=mod.lm$coeff[,"t value"] * mod.lm$coeff[,"Std. Error"] / SE.sandw,
              NULL), digits)
}


## SE.comparisons <- function(x, y, Nboot=10000, digits=3) {   # x must be a matrix or dataframe or a column
##   X <- cbind(icept=1,as.matrix(x));  Y <- y;  N <- nrow(X);  p <- ncol(X)
##   if(nrow(X) != length(y)) { cat("Number of X rows and length of y are different!\n");  return() }
##   mod.lm <- summary(lm(Y ~ ., data=as.data.frame(X)))
##   slopes.bs <- rbind(sapply(1:Nboot, function(i) { sel <- sample(nrow(X), replace=T);  if(i%%1000==0) cat(i,"")
##                                                    ## Xb <- scale(X[sel,],scale=F);  Yb <- Y[sel]
##                                                    Xb <- X[sel,];  Yb <- Y[sel]
##                                                    solve(t(Xb)%*%Xb)%*%t(Xb)%*%Yb } )); cat("\n") # rbind() for matrix
##   slopes.bs.SE <- apply(slopes.bs, 1, sd)
##   # sandwich estimator, corrected for design:
##   X3 <- X%*%solve(t(X)%*%X)
##   Hdiag <- diag(X%*%solve(t(X)%*%X)%*%t(X))
##   # SE.sandw <- sqrt(apply(resid(mod.lm)^2 * X3^2, 2, sum) * N/(N-p))  # global correction
##   SE.sandw <- sqrt(apply(resid(mod.lm)^2 * X3^2 / (1-Hdiag), 2, sum))  # local correction
##   # root of RAV sample estimate:
##   X.adj <- sapply(1:ncol(X), function(j) resid(lm(X[,j] ~ X[,-j])))
##   colnames(X.adj) <- colnames(X)
##   y.res <- resid(lm(y ~ X))
##   root.RAV <- sapply(colnames(X.adj), function(j) sqrt(mean(y.res^2 * X.adj[,j]^2)/mean(y.res^2)/mean(X.adj[,j]^2)) )
##   # List:
##   round(cbind("Slope"=mod.lm$coeff[,"Estimate"],
##               "SE.conv"=mod.lm$coeff[,"Std. Error"],
##               "SE.boot"=slopes.bs.SE,
##               "SE.sandw"=SE.sandw,
##               "boot/conv"=slopes.bs.SE/mod.lm$coeff[,"Std. Error"],
##               "sandw/conv"=SE.sandw/mod.lm$coeff[,"Std. Error"],
##               "sandw/boot"=SE.sandw/slopes.bs.SE,
##               "t.conv"=mod.lm$coeff[,"t value"],
##               "t.boot"=mod.lm$coeff[,"t value"] * mod.lm$coeff[,"Std. Error"] / slopes.bs.SE,
##               "t.sand"=mod.lm$coeff[,"t value"] * mod.lm$coeff[,"Std. Error"] / SE.sandw,
##               "root.RAV"=root.RAV,
##               NULL), digits)
## }