Skip to content

Instantly share code, notes, and snippets.

@briatte
Forked from dmarcelinobr/myboot.R
Last active December 18, 2015 21:39
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save briatte/5849358 to your computer and use it in GitHub Desktop.
Daniel Marcelino's bootstrap function for regression coefficients, made parallelizable + bootstrapped CIs
library(boot)
r2 = function(f, d, i) {
d = d[i, ]
f = lm(f, data = d)
return(summary(f)$r.square)
}
br <- boot(COI, r2, 10000, formula = nb.publis ~ nb.liens)
br
plot(br)
boot.ci(br, .95, "bca")
require(ggplot2)
require(reshape)
qog = read.csv("http://www.qogdata.pol.gu.se/data/qog_std_cs_20dec13.csv", sep = ";")
mod1 <- function(qog) lm(wdi_fr ~ bl_asy25f + gid_wip, data = qog)[[1]]
mod1.sds <- boot.lm(qog, "mod1", 100, hist = TRUE)
yhat = with(mod1.sds, data.frame(estimation))
names(yhat) = mod1.sds$names
yhat$x = row.names(yhat)
yhat = melt(yhat)
# ugly scaling
qplot(data = yhat, x = value, geom = "density") +
facet_grid(variable ~ ., scales = "free")
boot.lm <- function(data, stat, nreps, hist = TRUE, quiet = FALSE) {
require(foreach)
estimates <- get(stat)(data)
len <- length(estimates)
container <- matrix(NA, ncol = len , nrow = nreps)
nobs <- nrow(data)
if(!quiet) p <- txtProgressBar()
foreach(i = 1:nreps) %dopar% {
resample <- data[sample(rownames(data), nobs, replace = TRUE), ]
container[i, ] <- get(stat)(resample)
if(!quiet) setTxtProgressBar(p, i / nreps)
}
close(p)
sds <- apply(container, 2, sd)
if(hist) {
mfrow = c(1,1)
frame()
if(len <= 3) par(mfrow=c(len,1))
if((len > 3) & (len <= 6)) par(mfrow = c(3,2))
for(j in 1:len) hist(container[,j],
main=paste("Estimates for ", names(estimates)[j]), xlab="")
}
print(rbind(estimates, sds))
return(list(names = names(estimates), estimation = container, sds=sds))
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment