library(fixest)
library(haven)
library(knitr)
library(kableExtra)
library(magrittr)
auto<-read_dta("http://www.stata-press.com/data/r16/auto.dta")
lm1<-feols(mpg~weight+gear_ratio+foreign,data=auto)New summary function:
summary.bboot1 = function(object, level = 0.95,lmobject, ...) {
alpha = 1 - level
sd.v=apply(object,2,sd)
z.v=coef(lmobject)/sd.v
bias.v<-apply(object,2, mean)-coef(lmobject)
pval.v=2*(1-pt(abs(z.v),Inf , lower.tail = TRUE))
lb.q=apply(object,2,function(x) quantile(x,alpha/2,type=2)) # Looks like Stata uses quantile(...,type=2)?
ub.q=apply(object,2,function(x) quantile(x,1-alpha/2,type=2))
lb.n=coef(lmobject)-qt(p = 1-alpha/2,df=Inf,lower.tail = TRUE)*sd.v
ub.n=coef(lmobject)+qt(p = 1-alpha/2,df=Inf,lower.tail = TRUE)*sd.v
out<-data.frame("Observed Coef."=coef(lmobject),
"bias" = bias.v,
"Bootstrap_S.E."=sd.v,
"z"=z.v,
"p.val"=pval.v,
"LB"=lb.n,
"UB"=ub.n,
"CI type"="Normal Based")
out2<-data.frame("Observed Coef."=NA,
"bias" = NA,
"Bootstrap_S.E."=NA,
"z"=NA,
"p.val"=NA,
"LB"=lb.q,
"UB"=ub.q,
"CI type"="Qunatile")
rownames(out2)<-rownames(out)
out<-rbind(out,out2)
out<-out[order(row.names(out)),]
row.names(out)
out}Stata code:
use http://www.stata-press.com/data/r16/auto, clear
bootstrap, reps(100) seed(1) ///
saving(C:\R\boot.dta, every(1) double replace): ///
regress mpg weight gear foreign
estat bootstrap, all##
##
## . use http://www.stata-press.com/data/r16/auto, cl(1978 Automobile Data)
##
## . bootstrap, reps(100) seed(1) ///
## > saving(C:\R\boot.dta, every(1) double replace): ///
## > regress mpg weight gear foreign
## (running regress on estimation sample)
##
## Bootstrap replications (100)
## ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
## .................................................. 50
## .................................................. 100
##
## Linear regression Number of obs = 74
## Replications = 100
## Wald chi2(3) = 167.13
## Prob > chi2 = 0.0000
## R-squared = 0.6670
## Adj R-squared = 0.6527
## Root MSE = 3.4096
##
## ------------------------------------------------------------------------------
## | Observed Bootstrap Normal-based
## mpg | Coef. Std. Err. z P>|z| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## weight | -.006139 .0006063 -10.13 0.000 -.0073273 -.0049507
## gear_ratio | 1.457113 1.367917 1.07 0.287 -1.223954 4.138181
## foreign | -2.221682 1.169727 -1.90 0.058 -4.514305 .070942
## _cons | 36.10135 5.20581 6.93 0.000 25.89815 46.30455
## ------------------------------------------------------------------------------
##
## .
## . estat bootstrap, all
##
## Linear regression Number of obs = 74
## Replications = 100
##
## ------------------------------------------------------------------------------
## | Observed Bootstrap
## mpg | Coef. Bias Std. Err. [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## weight | -.00613903 .0001897 .00060628 -.0073273 -.0049507 (N)
## | -.0070678 -.004665 (P)
## | -.0075629 -.0050886 (BC)
## gear_ratio | 1.4571134 .2023129 1.3679166 -1.223954 4.138181 (N)
## | -.7289374 4.986959 (P)
## | -.7289374 4.986959 (BC)
## foreign | -2.2216815 .2298969 1.1697273 -4.514305 .070942 (N)
## | -4.298888 1.484208 (P)
## | -4.646578 .1750952 (BC)
## _cons | 36.101353 -1.284148 5.2058097 25.89815 46.30455 (N)
## | 24.38359 43.59073 (P)
## | 25.85769 49.81502 (BC)
## ------------------------------------------------------------------------------
## (N) normal confidence interval
## (P) percentile confidence interval
## (BC) bias-corrected confidence intervalReplicate the same output in R:
stata_bb<-read_dta("c:/R/boot.dta")
stata_bb<-stata_bb[,c(4,1,2,3)] # Move cons to first col
summary_df<-summary.bboot1(object=stata_bb,level = 0.95, lmobject = lm1)
summary_df[c(7,8,5,6,3,4,1,2),] %>% # Reorder rows to match stata's order
kbl() %>%
kable_classic_2(full_width = F)| Observed.Coef. | bias | Bootstrap_S.E. | z | p.val | LB | UB | CI.type | |
|---|---|---|---|---|---|---|---|---|
| weight | -0.006139 | 0.000190 | 0.000606 | -10.12571 | 0.000000 | -0.007327 | -0.004951 | Normal Based |
| weight1 | -0.007068 | -0.004665 | Qunatile | |||||
| gear_ratio | 1.457113 | 0.202313 | 1.367917 | 1.06521 | 0.286783 | -1.223954 | 4.138181 | Normal Based |
| gear_ratio1 | -0.728937 | 4.986959 | Qunatile | |||||
| foreign | -2.221682 | 0.229897 | 1.169727 | -1.89932 | 0.057523 | -4.514305 | 0.070942 | Normal Based |
| foreign1 | -4.298888 | 1.484208 | Qunatile | |||||
| (Intercept) | 36.101353 | -1.284148 | 5.205810 | 6.93482 | 0.000000 | 25.898153 | 46.304552 | Normal Based |
| (Intercept)1 | 24.383594 | 43.590725 | Qunatile |