重抽样与模拟

先查看各命令展开的示例,再按需展开其他选项及 R 实现。改写代码前,请阅读适用条件及差异说明。

bayesboot

使用贝叶斯自助Dirichlet权重估计斜率

Stata

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

R 包:stats

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

固定流程。独立指数权重归一化为总和一后服从Dirichlet(1,…,1);统一缩放不影响加权最小二乘系数。

bootstrap

在固定分层内自助重抽样回归系数

Stata

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

R 包:boot

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

另有 10 个示例及替代实现

对两个回归斜率之差进行自助估计

Stata

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

R 包: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)

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

通过重抽样完整聚类自助估计均值

Stata

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

R 包: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) }) }

固定流程。对聚类标识重抽样,并保留被抽中聚类中的全部观测;三个聚类仅用于演示。

对样本均值进行自助重抽样

Stata

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

R 包:boot

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

计算偏差校正加速自助置信区间

Stata

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

R 包:boot

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

对样本中位数进行自助重抽样

Stata

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

R 包:boot

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

对下四分位数进行自助重抽样

Stata

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

R 包:boot

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

对上四分位数进行自助重抽样

Stata

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

R 包:boot

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

对样本标准差进行自助重抽样

Stata

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

R 包:boot

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

通过重抽样观测对回归斜率进行自助估计

Stata

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

R 包:boot

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

对样本方差进行自助重抽样

Stata

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

R 包:boot

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

bsample

抽取非参数观测自助样本

Stata

bsample

R 包:base

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

抽取分层观测自助样本

Stata

bsample, strata(g)

R 包:base

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

bstat

汇总已保存的自助估计值、偏差与标准误

Stata

bstat

R 包:stats

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

软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

estat

计算正态自助置信区间

Stata

estat bootstrap, normal

R 包:boot

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

软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

另有 1 个示例及替代实现

计算百分位自助置信区间

Stata

estat bootstrap, percentile

R 包:boot

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

软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

jackknife

对回归斜率进行刀切估计

Stata

jackknife _b[x]: regress y x

R 包: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))) }

固定流程。删除单个观测的刀切法为确定性流程;方差等于(n-1)/n乘以重复估计值离均差的平方和。

另有 4 个示例及替代实现

枚举删除一个聚类的刀切重抽样

Stata

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

R 包:base

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

固定流程。返回重复次数并保留删除估计值;聚类大小不等时的方差校正需明确选择刀切设计。

生成样本均值的刀切伪值

Stata

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

R 包:base

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

固定流程。删除单个观测的刀切法为确定性流程;方差等于(n-1)/n乘以重复估计值离均差的平方和。

对样本均值进行刀切估计

Stata

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

R 包: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))) }

固定流程。删除单个观测的刀切法为确定性流程;方差等于(n-1)/n乘以重复估计值离均差的平方和。

对样本标准差进行刀切估计

Stata

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

R 包: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))) }

固定流程。删除单个观测的刀切法为确定性流程;方差等于(n-1)/n乘以重复估计值离均差的平方和。

permute

在固定分层内部置换结果变量

Stata

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

R 包: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"]) })

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

生成回归斜率的置换分布

Stata

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

R 包:stats

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

生成相关系数的置换分布

Stata

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

R 包:stats

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

固定流程。软件间随机序列和重抽样默认值不同;固定R种子使此流程可复现。

simulate

模拟正态样本均值的抽样分布

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 包:stats

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

固定流程。mymean是在Stata中定义的程序,抽取40个normal(2,1)观测并返回均值;此例为固定流程。