library(pscl)
# simulating data -------------------------------------------
set.seed(1839)
n <- 1500
x <- runif(n)
y <- cut(x + rnorm(n), 2, c("disagree", "agree"))

# fit model -------------------------------------------------
model <- glm(y ~ x, family = binomial("logit"))
pR2(model)[["McFadden"]]

# bootstrap -------------------------------------------------
number_of_resamples <- 5000
bootstrap_pr2 <- sapply(1:number_of_resamples, function(j) {
  bootstrap_cases <- sample(c(1:n), n, TRUE)
  pR2(glm(y[bootstrap_cases] ~ x[bootstrap_cases],
          family = binomial("logit")))[["McFadden"]]
})

# report 95% CI ---------------------------------------------
mean(bootstrap_pr2) - 1.96 * sd(bootstrap_pr2)
mean(bootstrap_pr2) + 1.96 * sd(bootstrap_pr2)