.Rhistory

gcdnet?
)
install.packages('data.table')
library(data.table)
library(MASS)
library(Matrix)
set.seed(0808)
n<-1000
q<-2000
q10<-50
q01<-100
q11<-50 #the number of active mediators
p1<-10 #the number of exposures in one toxicant class
rho1<-0.2
covA<-matrix(rho1,p1,p1); diag(covA)<-1
A<-mvrnorm(n,rep(0,p1),covA)
gamma<-c(rep(0.2,p1/2),rep(0,p1/2))
ers<-A%*%gamma + rnorm(n,0,1)
length(ers)
rho2<-.5
covM<-matrix(0,q,q)
block_size = 20; times <- 1:block_size #simulate 50 biological pathways for mediators, each contains 20 mediators
nblock = 50
start = 11 + 40*(c(1:nblock)-1); end = start + block_size-1
for (i in c(1:nblock)) {
covM[start[i]:end[i], start[i]:end[i]] <- rho2
}
diag(covM) <- 1
Sigma11<-matrix(c(0.5,0.2,0.2,0.5),2,2)
Sigma01<-matrix(c(0,0,0,0.5),2,2)
Sigma10<-matrix(c(0.5,0,0,0),2,2)
param <- matrix(0,q,2)
mu<-c(0,0)
idx11<-c(c(start[7]:end[7])[c(1:10)], c(start[17]:end[17])[c(1:10)], c(start[27]:end[27])[c(1:10)], c(start[37]:end[37])[c(1:10)], c(start[47]:end[47])[c(1:10)]) #index for active mediators
range<-1:q
range10<-range[!range %in% idx11]
idx10<-sample(range10, q10, replace = FALSE)
range01<-range[!range %in% c(idx11,idx10)]
idx01<-sample(range01, q01, replace = FALSE)
param[idx11,] <- mvrnorm(q11,mu,Sigma11)
param[idx10,] <- mvrnorm(q10,mu,Sigma10)
param[idx01,] <- mvrnorm(q01,mu,Sigma01)
beta_m <- param[,1]; alpha_a <- param[,2]
E<-mvrnorm(n,rep(0, q),covM)
M<-ers %*% alpha_a + E
beta_a<-0.5
e<-matrix(rnorm(n,0,1),n,1)
Y<-beta_a*ers + M %*% beta_m + e
dim(Y)
hist(Y)
dim(E)
dim(M)
dim(e)
dim(E)
hist(ers)
dim(A)
data<-cbind(A,M,Y)
dim(dat)
dim(data)
write.csv(data,"sample.data.csv")
names(data)
colnames(data)
head(data)