# Answer set.seed(123321) nrep <- 1000 b0 <- 0 bs <- c(0, 0.5, 1, 1.5) ns <- c(50, 100, 200, 400) ovresult <- NULL for(k in 1:length(bs)){ b1 <- bs[k] for(j in 1:length(ns)) { n <- ns[j] result <- matrix(NA, nrep, 3) for(i in 1:nrep) { x <- rnorm(n, 0, 1) x <- rnorm(n, 0, 1) fy <- b0 + b1*x p <- exp(fy)/(1 + exp(fy)) y <- rep(NA, n) for(w in 1:n) y[w] <- rbinom(1, 1, p[w]) dat <- data.frame(x = x, y = y) out <- glm(y ~ x, data = dat, family = binomial(link = "logit")) sumout <- summary(out) result[i,] <- sumout[["coefficients"]][2, c(1, 2, 4)] } result <- cbind(b1, n, result) ovresult <- rbind(ovresult, result) } } sig <- ovresult[,5] < 0.05 ovresult <- cbind(ovresult, sig) colnames(ovresult) <- c("b", "n", "est", "se", "p", "sig") aggresult <- aggregate(cbind(est, se, sig) ~ n + b, data = ovresult, FUN = mean) # average aggsdresult <- aggregate(est ~ n + b, data = ovresult, FUN = sd) # sd biasest <- aggresult[,"est"] - aggresult[,"b"] relbiasest <- biasest / aggresult[,"b"] biasse <- aggresult[,"se"] - aggsdresult[,"est"] relbiasse <- biasse / aggsdresult[,"est"] power <- aggresult[,"sig"] tableresult <- data.frame(aggresult[,c(1, 2)], biasest, relbiasest, biasse, relbiasse, power) tableresult equal <- function(a, b) { abs(a - b) < 0.000001 } putLine <- function(mat, lty, pch, col) { points(mat, pch = pch, col = col) lines(mat, lty = lty, col = col) } plotResult <- function(mat, targetcol, aty, ylab, poslegend, hline) { lty <- 1:4 pch <- c(20, 19, 18, 17) col <- c("black", "red", "darkgreen", "orange") legendtxt <- c("b = 0", "b = 0.5", "b = 1", "b = 1.5") plot(c(0,0), type = "n", xlim=c(40, 400), ylim=range(aty), xlab="Sample Size", ylab=ylab, main = "My Simulation", axes = FALSE) axis(1, at=c(50, 100, 200, 400)) axis(2, at=aty) box() putLine(mat[equal(mat[,"b"], 0), c("n", targetcol)], lty = lty[1], pch = pch[1], col = col[1]) putLine(mat[equal(mat[,"b"], 0.5), c("n", targetcol)], lty = lty[2], pch = pch[2], col = col[2]) putLine(mat[equal(mat[,"b"], 1), c("n", targetcol)], lty = lty[3], pch = pch[3], col = col[3]) putLine(mat[equal(mat[,"b"], 1.5), c("n", targetcol)], lty = lty[4], pch = pch[4], col = col[4]) abline(h = hline, lty = 2) legend(x = poslegend[1], y = poslegend[2], legend = legendtxt, pch = pch, lty = lty, col = col) } par(mfrow=c(3,1)) plotResult(tableresult, "relbiasest", seq(-0.02, 0.1, 0.02), ylab = "Relative Bias", poslegend = c(300, 0.1), hline = 0) plotResult(tableresult, "relbiasse", seq(-0.1, 0.02, 0.02), ylab = "Relative SE Bias", poslegend = c(300, -0.05), hline = 0) plotResult(tableresult, "power", seq(0, 1, 0.2), ylab = "Power", poslegend = c(300, 0.7), hline = 0.8)