# R for simulation study ### Generate a data set and analyze it using multiple regression library(MASS) set.seed(123321) n <- 100 miv <- rep(0, 3) siv <- matrix(0.5, 3, 3) diag(siv) <- 1 dat <- mvrnorm(n, miv, siv) y <- 0.3*dat[,1] + 0.4*dat[,2] + 0.2*dat[,3] + rnorm(n, 0, sqrt(0.45)) dat <- cbind(dat, y) colnames(dat) <- c("x1", "x2", "x3", "y") dat <- as.data.frame(dat) out <- lm(y ~ x1 + x2 + x3, data = dat) sumout <- summary(out) names(sumout) sumout[["coefficients"]] out2 <- lm(y ~ x1 + x2, data = dat) anova(out2, out) ### Generate multiple data sets nrep <- 3 n <- 100 miv <- rep(0, 3) siv <- matrix(0.5, 3, 3) diag(siv) <- 1 dat.l <- list() set.seed(123321) for(i in 1:nrep) { dat <- mvrnorm(n, miv, siv) y <- 0.3*dat[,1] + 0.4*dat[,2] + 0.2*dat[,3] + rnorm(n, 0, sqrt(0.45)) dat <- cbind(dat, y) colnames(dat) <- c("x1", "x2", "x3", "y") dat.l[[i]] <- as.data.frame(dat) } ### Analyze them using multiple regression by for loop resultcoef <- matrix(NA, nrep, 4) for(i in 1:nrep) { out <- lm(y ~ x1 + x2 + x3, data = dat.l[[i]]) sumout <- summary(out) resultcoef[i,] <- sumout\$coefficients[,1] } apply(resultcoef, 2, mean) ### sapply function sapply(dat.l, nrow) sapply(dat.l, colMeans) dvstat <- function(dat) c(mean(dat[,4]), sd(dat[,4])) sapply(dat.l, dvstat) sapply(dat.l, function(dat) c(mean(dat[,4]), sd(dat[,4]))) analysisfun <- function(dat) { out <- lm(y ~ x1 + x2 + x3, data = dat) sumout <- summary(out) sumout\$coefficients[,1] } sapply(dat.l, analysisfun) ### Time it ########## Data Generation ########## nrep2 <- 10000 miv <- rep(0, 3) siv <- matrix(0.5, 3, 3) diag(siv) <- 1 dat.l2 <- list() set.seed(123321) for(i in 1:nrep2) { dat <- mvrnorm(n, miv, siv) y <- 0.3*dat[,1] + 0.4*dat[,2] + 0.2*dat[,3] + rnorm(n, 0, sqrt(0.45)) dat <- cbind(dat, y) colnames(dat) <- c("x1", "x2", "x3", "y") dat.l2[[i]] <- as.data.frame(dat) } ########## For loop ########## system.time({ resultcoef <- matrix(NA, nrep2, 4) for(i in 1:nrep2) { out <- lm(y ~ x1 + x2 + x3, data = dat.l2[[i]]) sumout <- summary(out) resultcoef[i,] <- sumout\$coefficients[,1] } }) ########## sapply ########## system.time(resultcoef2 <- sapply(dat.l2, analysisfun)) ### lapply lapply(dat.l, colMeans) lapply(dat.l, cor) lapply(dat.l, cor, method = "spearman") output.l <- lapply(dat.l, lm, formula = y ~ x1 + x2 + x3) output.l sumoutput.l <- lapply(output.l, summary) sumoutput.l extractcoef <- function(out) out[["coefficients"]][,2] coef.l <- lapply(sumoutput.l, extractcoef) coef.l do.call(rbind, coef.l) coef.l <- lapply(sumoutput.l, "[[", "coefficients") sapply(coef.l, "[", 1:4, 1) residual.l <- lapply(sumoutput.l, "[[", "residuals") ### mapply output1.l <- lapply(dat.l, lm, formula = y ~ x1 + x2 + x3) output2.l <- lapply(dat.l, lm, formula = y ~ x1 + x2) mapply(anova, output1.l, output2.l, SIMPLIFY = FALSE) diffrsquared <- function(out1, out2, adjust = FALSE) { name <- "r.squared" if(adjust) name <- "adj.r.squared" rsquared1 <- summary(out1)[[name]] rsquared2 <- summary(out2)[[name]] abs(rsquared1 - rsquared2) } mapply(diffrsquared, output1.l, output2.l) mapply(diffrsquared, output1.l, output2.l, MoreArgs = list(adjust = TRUE)) corpredict <- function(out1, out2) { cor(predict(out1), predict(out2)) } mapply(corpredict, output1.l, output2.l) ### Parallel Processing output.l <- lapply(dat.l, lm, formula = y ~ x1 + x2 + x3) library(parallel) output.l <- mclapply(dat.l, lm, formula = y ~ x1 + x2 + x3, mc.cores = 2) cl <- makeCluster(rep("localhost", 2), type = "PSOCK") output.l <- parLapply(cl, dat.l, lm, formula = y ~ x1 + x2 + x3) stopCluster(cl) detectCores() ### Exercises g1 <- cbind(1, rnorm(100, 0, 1)) g2 <- cbind(2, rnorm(100, 0.5, 1)) dat <- data.frame(rbind(g1, g2)) colnames(dat) <- c("group", "y") outeq <- t.test(y ~ group, data = dat, var.equal = TRUE) outneq <- t.test(y ~ group, data = dat) cieq <- outeq[["conf.int"]] cineq <- outneq[["conf.int"]] widtheq <- cieq[2] - cieq[1] widthneq <- cineq[2] - cineq[1] diffwidth <- widtheq - widthneq