get.data <- function(n, k){ X <- matrix(rnorm(n*k, mean =0, sd= 1), nrow=n, ncol = k) epsilon <- rnorm(n, mean =0, sd = 1) indexp <- rowSums(X[,(k/2):k])/sqrt(k/2) p <- exp(indexp)/ (1+exp(indexp)) u <- runif(n, min=0, max =1) D <- ifelse(p>u, 1, 0) y1 <- 10 + rowSums(X)/sqrt(k) + epsilon Y <- D*y1 cbind(Y,D,X) }