Simulating Confidence Intervals
[This article was first published on The PolStat Feed, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
As I am finishing up my thesis I have recently been plotting effects from many models. An important aspect of this is to show the uncertainty surrounding different estimates and effects. Following a paper by Gary King, Michael Tomz and Jason Wittenberg it has become very popular among quants in political science to simulate the uncertainty rather than relying on the delta method. They take advantage of the fact that the central limit theorem show that with a large enough sample and bounded variance it is possible to simulate parameter values by drawing from a multivariate normal distribution. The means are set equal to the estimated parameters and variance equal to the variance-covariance matrix from the model.
Stata users have access to the excellent Clarify module from Gary King and colleagues which implements their approrach, however since I am doing most of my work in R this is not an option. King and friends have implemented much of the same functionality in the Zelig package for R, but I find extracting exactly what I want from objects created by Zelig to be a chore. Luckily it is very easy to implement this on your own in R. The approach advocated by King and colleagues follows a 5 step process:
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
- Simulate the parameters
- define the values at which the variables are held constant
- Calculate the systematic component of model for each round of simulated parameters
- Use the systematic component to calculate your quantity of interest
- repeat step 1-4 a 1000 times, or until you have the desired degree of accuracy
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(MASS) | |
library(ggplot2) | |
#generate the fake data | |
covvar <- matrix(c(1,.6,.6,.6,1,.6,.6,.6,1),nr = 3) | |
means <- c(20,25,30) | |
data <- as.data.frame(mvrnorm(100,means,covvar)) | |
#dichotomize the dependent variable | |
data$V1 <- ifelse(data$V1 > 20,1,0) | |
#estimate the logit model | |
model1 <- glm(V1 ~ V2 + V3, data = data, family = "binomial") | |
summary(model1) | |
#get the components necessary for the simulation | |
beta <- coef(model1) | |
covvar <- vcov(model1) | |
#define the values at which the variables are held constant, V3 is allowed to vary with V2 at its mean | |
x <- cbind(1,mean(data$V2), seq(27,33,0.1)) | |
#simulate parameter values from the multivariate normal distribution, 1000 draws | |
beta.sim <- mvrnorm(1000, beta,covvar) | |
#use some matrix algebra to calculate to systematic component for each draw | |
xb <- x %*% t(beta.sim) | |
#for each draw get prediced probability | |
p <- 1/(1 + exp(-xb)) | |
#get the .5 and .95 confidence intervals for each value in the sequence for V3 | |
predicted.prob <- apply(p,1,quantile, probs = c(.05,.5,.95)) | |
#collect the results in a data frame for plotting in ggplot2 | |
dat <- data.frame(x = seq(27,33,0.1), lower = predicted.prob[1,], middle = predicted.prob[2,], upper = predicted.prob[3,]) | |
#plot the results with ggplot2 | |
plot <- ggplot(dat, aes(x = x)) | |
plot <- plot + geom_line(aes(y = lower), linetype = 2) | |
plot <- plot + geom_line(aes(y = middle), linetype = 1) | |
plot <- plot + geom_line(aes(y = upper), linetype = 2) | |
plot <- plot + theme_bw() | |
plot |

To leave a comment for the author, please follow the link and comment on their blog: The PolStat Feed.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.