########################################################################
# NAME: Tom Loughin #
# DATE: 08-26-2013 #
# PURPOSE: Simulate random effects for logistic regression and compare #
# curve with average random effect to average curve #
# #
# NOTES: #
########################################################################
set.seed(3789022)
beta1 <- 0.2
beta0 <- -100*beta1
b <- rnorm(n = 50, mean = 0, sd = 1.5)
x <- seq(from = 70, to = 130, by = .5)
probs <- matrix(data = 0, ncol = 50, nrow = length(x))
x11(height = 7, width = 12)
# pdf(file = "c:\\figures\\Figure6.8color.pdf", width = 12, height = 6, colormodel = "cmyk", pointsize = 14) # Create plot for book
curve(expr = plogis(beta0 + beta1*x), from = 70, to = 130, col = "red", lwd = 2, ylab = "P(Y = 1)")
for(i in c(1:50)){
curve(expr = plogis((beta0 + b[i]) + beta1*x), from = 70, to = 130,
col = "black", lwd = 1, lty = "dotted", add = TRUE)
probs[,i] <- plogis((beta0 + b[i]) + beta1*x)
}
curve(expr = plogis(beta0 + beta1*x), from = 70, to = 130, col = "red", lwd = 3, add = TRUE)
marg <- apply(X = probs, MARGIN = 1, FUN = mean)
lines(x = x, y = marg, lty = "dashed", lwd = 3, col = "blue")
legend(x = 115, y = 0.3, legend = c("Subject-specific", "Marginal", "Individual"), bty = "n",
lty = c("solid", "dashed", "dotted"), lwd = c(3,3,1), col = c("red", "blue", "black"))
# dev.off() # Create plot for book
# Black-and-white version of plot
# pdf(file = "c:\\figures\\Figure6.8BW.pdf", width = 12, height = 6, colormodel = "cmyk", pointsize = 14) # Create plot for book
curve(expr = plogis(beta0 + beta1*x), from = 70, to = 130, col = "red", lwd = 2, ylab = "P(Y = 1)")
for(i in c(1:50)){
curve(expr = plogis((beta0 + b[i]) + beta1*x), from = 70, to = 130,
col = "black", lwd = 1, lty = "dotted", add = TRUE)
probs[,i] <- plogis((beta0 + b[i]) + beta1*x)
}
curve(expr = plogis(beta0 + beta1*x), from = 70, to = 130, col = "black", lwd = 3, add = TRUE)
marg <- apply(X = probs, MARGIN = 1, FUN = mean)
lines(x = x, y = marg, lty = "dashed", lwd = 3, col = "black")
legend(x = 115, y = 0.3, legend = c("Subject-specific", "Marginal", "Individual"), bty = "n",
lty = c("solid", "dashed", "dotted"), lwd = c(3,3,1), col = c("black", "black", "black"))
# dev.off() # Create plot for book