set.seed(1234)
p  <- 450000
n  <- 200
np <- n*p

## Generate the data
tt <- system.time( B <- matrix(rbeta(np,0.5,0.5),nrow=p,ncol=n) )

## Generate the covariate
age <- round(rnorm(n,mean=45,sd=10))
sex <- rbinom(n,size=1,prob=0.5)
bmi <- sex*rnorm(n,mean=23,sd=3) + (1-sex)*rnorm(n,mean=28,sd=5)
t2d <- rep(0:1,each=n/2)

## Linear regression
## B~t2d+age+sex+bmi


## Solution 0: for loop ----
mod <- NULL
tt0 <- system.time( 
for(j in 1:10000){
    mod <- rbind(mod,lm(B[j,]~t2d+age+sex+bmi)$coefficients)
})
  

## Solution 1: sapply
tt1 <- system.time( mod <- sapply(1:10000,function(j) lm(B[j,]~t2d+age+sex+bmi)$coefficients) )

## Solution 2: mclapply
library(parallel)
tt2a <- system.time( mod <- mclapply(1:10000,function(j) lm(B[j,]~t2d+age+sex+bmi)$coefficients,mc.cores=4) )
tt2b <- system.time( mod <- mclapply(1:10000,function(j) lm(B[j,]~t2d+age+sex+bmi)$coefficients,mc.cores=10) )

## Solution 3: Matrix lm
tt3 <- system.time( mod <- lm(t(B[1:10000,])~t2d+age+sex+bmi)$coefficients )

## Solution 3: matrix manipulation
tt4 <- system.time({
  x      <- cbind(Intercept=1,T2D=t2d,Age=age,Sex=sex,BMI=bmi)
  xTx    <- t(x)%*%x
  G      <- solve(xTx)%*%t(x)
  Betas  <- B[1:10000,]%*%t(G)
})

## Full Dataset
tt3f <- system.time( mod <- lm(t(B)~t2d+age+sex+bmi)$coefficients )

## Solution 3: matrix manipulation
tt4f <- system.time({
  x      <- cbind(Intercept=1,T2D=t2d,Age=age,Sex=sex,BMI=bmi)
  xTx    <- t(x)%*%x
  G      <- solve(xTx)%*%t(x)
  Betas  <- B%*%t(G)
})