Resampling and simulation

Start with the expanded example for each command. Open the variants for options and alternative R approaches. Read the qualifications before adapting the code.

bayesboot

Use Bayesian-bootstrap Dirichlet weights for a slope

Stata

bayesboot _b[x], reps(199): regress y x

R Package:stats

replicate(199,{ weights <- stats::rexp(nrow(data)); unname(stats::coef(stats::lm(y~x,data=data,weights=weights))["x"]) })

Fixed workflow. Independent exponential weights normalized to sum one are Dirichlet(1,…,1); common scaling does not change weighted least-squares coefficients.

bootstrap

Bootstrap regression coefficients within fixed strata

Stata

bootstrap _b, reps(199) strata(g): regress y x

R Package:boot

boot::boot(data,statistic=function(d,i) stats::coef(stats::lm(y~x,data=d[i,])),R=199,strata=data$g)

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

10 more examples and alternatives

Bootstrap a difference between regression slopes

Stata

bootstrap difference=(_b[x]-_b[z]), reps(199): regress y x z

R Package:boot

boot::boot(data,statistic=function(d,i) unname(diff(stats::coef(stats::lm(y~x+z,data=d[i,]))[c("z","x")])),R=199)

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

Bootstrap a mean by resampling complete clusters

Stata

bootstrap mean=r(mean), reps(199) cluster(g): summarize y

R Package:base

{ ids <- unique(data$g); replicate(199,{ picked <- sample(ids,length(ids),replace=TRUE); sample_data <- do.call(rbind,lapply(picked,function(id) data[data$g==id,,drop=FALSE])); mean(sample_data$y) }) }

Fixed workflow. Resamples cluster identifiers, retaining every observation in each selected cluster. Three clusters are only a demonstration.

Bootstrap the sample mean

Stata

bootstrap mean=r(mean), reps(199): summarize y

R Package:boot

boot::boot(data,statistic=function(d,i) mean(d$y[i]),R=199)

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

Compute a bias-corrected accelerated bootstrap interval

Stata

bootstrap mean=r(mean), reps(499) bca: summarize y

R Package:boot

{ boot_fit <- boot::boot(data,function(d,i) mean(d$y[i]),R=499); boot::boot.ci(boot_fit,type="bca") }

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

Bootstrap the sample median

Stata

bootstrap median=r(p50), reps(199): summarize y, detail

R Package:boot

boot::boot(data,statistic=function(d,i) stats::median(d$y[i]),R=199)

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

Bootstrap the lower quartile

Stata

bootstrap q1=r(p25), reps(199): summarize y, detail

R Package:boot

boot::boot(data,statistic=function(d,i) unname(stats::quantile(d$y[i],0.25,type=2)),R=199)

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

Bootstrap the upper quartile

Stata

bootstrap q3=r(p75), reps(199): summarize y, detail

R Package:boot

boot::boot(data,statistic=function(d,i) unname(stats::quantile(d$y[i],0.75,type=2)),R=199)

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

Bootstrap the sample standard deviation

Stata

bootstrap sd=r(sd), reps(199): summarize y

R Package:boot

boot::boot(data,statistic=function(d,i) stats::sd(d$y[i]),R=199)

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

Bootstrap a regression slope by resampling rows

Stata

bootstrap slope=_b[x], reps(199): regress y x

R Package:boot

boot::boot(data,statistic=function(d,i) unname(stats::coef(stats::lm(y~x,data=d[i,]))["x"]),R=199)

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

Bootstrap the sample variance

Stata

bootstrap variance=r(Var), reps(199): summarize y

R Package:boot

boot::boot(data,statistic=function(d,i) stats::var(d$y[i]),R=199)

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

bsample

Draw a nonparametric row-bootstrap sample

Stata

bsample

R Package:base

data[sample.int(nrow(data),nrow(data),replace=TRUE),,drop=FALSE]

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

Draw a stratified row-bootstrap sample

Stata

bsample, strata(g)

R Package:base

do.call(rbind,lapply(split(data,data$g),function(d) d[sample.int(nrow(d),nrow(d),replace=TRUE),,drop=FALSE]))

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

bstat

Summarize stored bootstrap estimates, bias and standard error

Stata

bstat

R Package:stats

c(estimate=boot_fit$t0,bias=mean(boot_fit$t)-boot_fit$t0,se=stats::sd(boot_fit$t))

Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

estat

Compute a normal bootstrap interval

Stata

estat bootstrap, normal

R Package:boot

boot::boot.ci(boot_fit,type="norm")

Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

1 more examples and alternatives

Compute a percentile bootstrap interval

Stata

estat bootstrap, percentile

R Package:boot

boot::boot.ci(boot_fit,type="perc")

Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

jackknife

Jackknife a regression slope

Stata

jackknife _b[x]: regress y x

R Package:stats

{ n <- nrow(data); deleted <- vapply(seq_len(n),function(i) unname(stats::coef(stats::lm(y~x,data=data[-i,]))["x"]),numeric(1)); c(estimate=unname(stats::coef(stats::lm(y~x,data=data))["x"]),se=sqrt((n-1)/n*sum((deleted-mean(deleted))^2))) }

Fixed workflow. The delete-one jackknife is deterministic; its variance uses (n-1)/n times the squared deviations of the replicate estimates.

4 more examples and alternatives

Enumerate delete-one-cluster jackknife replicates

Stata

jackknife mean=r(mean), rclass cluster(g): summarize y

R Package:base

{ groups <- unique(data$g); deleted <- vapply(groups,function(id) mean(data$y[data$g!=id]),numeric(1)); length(groups) }

Fixed workflow. Returns replicate count and retains deleted estimates. Unequal cluster-size variance corrections need an explicitly chosen jackknife design.

Generate jackknife pseudovalues for the sample mean

Stata

jackknife mean=r(mean), rclass keep: summarize y

R Package:base

{ n <- nrow(data); deleted <- vapply(seq_len(n),function(i) mean(data$y[-i]),numeric(1)); n*mean(data$y)-(n-1)*deleted }

Fixed workflow. The delete-one jackknife is deterministic; its variance uses (n-1)/n times the squared deviations of the replicate estimates.

Jackknife the sample mean

Stata

jackknife mean=r(mean), rclass: summarize y

R Package:stats

{ n <- nrow(data); deleted <- vapply(seq_len(n),function(i) mean(data$y[-i]),numeric(1)); c(estimate=mean(data$y),se=sqrt((n-1)/n*sum((deleted-mean(deleted))^2))) }

Fixed workflow. The delete-one jackknife is deterministic; its variance uses (n-1)/n times the squared deviations of the replicate estimates.

Jackknife the sample standard deviation

Stata

jackknife sd=r(sd), rclass: summarize y

R Package:stats

{ n <- nrow(data); deleted <- vapply(seq_len(n),function(i) stats::sd(data$y[-i]),numeric(1)); c(estimate=stats::sd(data$y),se=sqrt((n-1)/n*sum((deleted-mean(deleted))^2))) }

Fixed workflow. The delete-one jackknife is deterministic; its variance uses (n-1)/n times the squared deviations of the replicate estimates.

permute

Permute outcomes within fixed strata

Stata

permute y _b[x], reps(199) strata(g): regress y x

R Package:stats

replicate(199,{ d <- data; d$y <- unsplit(lapply(split(data$y,data$g),sample),data$g); unname(stats::coef(stats::lm(y~x,data=d))["x"]) })

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

Generate a permutation distribution for a regression slope

Stata

permute y _b[x], reps(199): regress y x

R Package:stats

replicate(199,{ d <- data; d$y <- sample(d$y); unname(stats::coef(stats::lm(y~x,data=d))["x"]) })

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

Generate a permutation distribution for a correlation

Stata

permute y r(rho), reps(199): correlate y x

R Package:stats

replicate(199,stats::cor(sample(data$y),data$x))

Fixed workflow. Random streams and resampling defaults differ across software; a fixed R seed makes this workflow reproducible.

simulate

Simulate the sampling distribution of a normal sample mean

Stata

program define mymean, rclass
    preserve
    clear
    set obs 40
    generate y = rnormal(2,1)
    summarize y
    local estimate = r(mean)
    restore
    return scalar mean = `estimate'
end
simulate mean=r(mean), reps(199): mymean

R Package:stats

replicate(199,mean(stats::rnorm(40,mean=2,sd=1)))

Fixed workflow. mymean is a user-defined Stata program drawing 40 normal(2,1) observations and returning their mean; this is a fixed workflow.