simhet <- function(n,p,nreps,sdpow) { bh1s <- vector(length=nreps) ses <- vector(length=nreps) for (i in 1:nreps) { x <- matrix(rnorm(n*p),ncol=p) meany <- x %*% rep(1,p) sds <- abs(meany)^sdpow y <- meany + rnorm(n,sd=sds) lmout <- lm(y ~ x) bh1s[i] <- coef(lmout)[2] ses[i] <- sqrt(vcov(lmout)[2,2]) } mean(abs(bh1s - 1.0) < 1.96*ses) }