Replicating Stata Command bsample, cluster(nr) into R

I am replicating a paper by Clément de Chaisemartin: Two-Way Fixed Effects Estimators with Heterogeneous Treatment Effects.

I face a problem with Stata command bsample - Sampling with replacement.

I cannot understand what bsample, cluster(nr) Stata command is doing in for loop and what its results are.

I need to convert bsample, cluster(nr) Stata code into R.

sample data view

Stata code below:

set seed 1
quietly {
forvalue i=1/200 {
    preserve
    bsample, cluster(nr)
    xtreg lwage d81-d87 union, fe robust
    scalar betafe=_b[union]
    reg diff_lwage diff_union d82-d87, cluster(nr)
    scalar betafd=_b[diff_union]
    did_multiplegt lwage nr year union
    scalar did_m=e(effect_0)
    matrix B=B\(betafe-betafd,betafe-did_m,betafd-did_m)
    restore
}
}

The State command bsample, cluster(nr) is doing a clustered bootstrap - it resamples cluster ids and then takes all the observations in each resampled cluster. While there would be lots of ways to do this in R, one is:

ids <- na.omit(unique(data$nr))
for(i in 1:nBootstraps){
  samp <- sample(ids, length(ids), replace=TRUE)
  samps <- lapply(samp, function(x)data[which(nr == x), ])
  bs_data <- do.call(rbind, samps)
  ... other stuff ...
}

Where ... other stuff ... are the other functions executed on the bootstrapped data.