Beware: This is (almost) a cross-post to a thread I started on the Statalist but that has not received much attention so far.


I am learning about the problems when conducting hypothesis tests on a cluster-sample with very few clusters (<30) but considerable within-cluster correlation. So far, I read the work of Cameron/Gelbach/Miller "Bootstrap-Based Improvements for Inference with Clustered Errors (Review of Economics and Statistics 90, 414–427) [Working Paper here] as well as Cameron and Miller's "Practitioner’s Guide to Cluster-Robust Inference" (Journal of Human Resources 50, 317–370) [ Preprint here]. From what I learned if the number of clusters is small one would reject the null hypothesis too often and bootstraping procedures (possibly with asymptotic refinement) can help to overcome this problem. In their simulations a wild cluster bootstrap t procedure works best with rejection rates very close to the nominal 5%. My questions are about transferring their ideas to binary response models. To provide a basis for discussion I did some simulations in Stata to find out about the rejection rates of the null hypothesis.

Data Generating Process

I assumed the following data generating process for the "latent" variable $y^*$:


where $z_g$ is a standard random normal variable constant for any group $g$ and $z_{ig}$ is an independent random draw from the standard normal. To induce group dependent errors as well as heteroskedasticity I assume $u_{ig}$ to be a random draw from a logistic distribution with parameter $b$ defined as follows:


If $y^* > 0$ the binary response variable y takes on value 1, otherwise it is 0.


I tried to learn about the problem with simulated data. I simulated 499 datasets according to the data generating process described above, estimated a logit model and counted how often $H_0: \beta = 1$ is rejected. I assumed 15 groups, each comprising 30 observations. I estimated logit models:

  1. with no adjustment,
  2. a cluster-robust estimate for the variance-covariance matrix,
  3. a cluster bootstrap estimate for the variance-covariance matrix,
  4. a clustered bootstrap on the wald statistics.

Here is the Stata programme I wrote:

program define mysimu, rclass
    version 13.1

    syntax [, groups(integer 15) obspgroup(integer 30) alpha(real 0.05) bootstraprep(integer 499)]

    tempname cmat vmat cmat2 vmat2 cmat3 vmat3 cmat4 vmat4 b se w rej se_clust w_clust rej_clust se_boot w_boot rej_boot 
    tempvar group newgroup helpvar helpvar2 z_g helpvar3 e_g parb u x ystar y
    //Set the number of observations
    set obs `groups'
    gen `group' = _n    
    expand `obspgroup'
    bysort `group': gen `helpvar' = _n

    //Cluster-Level Regressor Correlation
    gen `helpvar2' = rnormal() if `helpvar' == 1
    bysort `group': egen `z_g' = min(`helpvar2')

    //Cluster-Level Error Correlation
    gen `helpvar3' = runiform() if `helpvar' == 1   
    bysort `group': egen `e_g' = min(`helpvar3')

    //Data Generating Process
    generate `parb' = sqrt(27*`z_g'^2/_pi^2)
    generate `u' = -`parb'*log(1/uniform() - 1)
    generate `x' = rnormal()+`z_g'
    generate `ystar' = `x'+`u'
    generate `y' = `ystar' > 0

    //Simple Logit
    logit `y' `x'
    matrix `cmat' = e(b)
    matrix `vmat' = e(V)
    return scalar b = `cmat'[1,1]
    return scalar se = sqrt(`vmat'[1,1])
    return scalar w = (`cmat'[1,1]-1)^2/`vmat'[1,1]
    return scalar rej = (`cmat'[1,1]-1)^2/`vmat'[1,1] >  invchi2(1, 1-`alpha')

    //Clustered SE
    logit `y' `x', cluster(`group') 
    matrix `cmat2' = e(b)
    matrix `vmat2' = e(V)
    return scalar b_clust = `cmat2'[1,1]
    return scalar se_clust = sqrt(`vmat2'[1,1])
    return scalar w_clust = (`cmat2'[1,1]-1)^2/`vmat2'[1,1]
    return scalar rej_clust = (`cmat2'[1,1]-1)^2/`vmat2'[1,1] >  invchi2(1, 1-`alpha')

    //Pairs clustered bootstrap se
    logit `y' `x', vce(boot, reps(`bootstraprep') cluster(`group'))
    matrix `cmat3' = e(b)
    matrix `vmat3' = e(V)
    return scalar b_boot = `cmat3'[1,1] 
    return scalar se_boot = sqrt(`vmat3'[1,1])
    return scalar w_boot = (`cmat3'[1,1]-1)^2/`vmat3'[1,1]
    return scalar rej_boot = (`cmat3'[1,1]-1)^2/`vmat3'[1,1] >  invchi2(1, 1-`alpha')

    //Pairs clustered bootstrap wald test
    tempfile dataset
    logit `y' `x', cluster(`group')
    matrix `cmat4' = e(b)
    local theta = `cmat4'[1,1]
    matrix `vmat4' = e(V)
    local w = (`cmat4'[1,1]-1)^2/`vmat4'[1,1]
    bootstrap wstar=((_b[`x']-`theta')^2/(_se[`x'])^2), reps(`bootstraprep') cluster(`group') idcluster(`newgroup') saving(`dataset', replace): logit `y' `x', cluster(`newgroup')
    use `dataset', clear
    sum wstar, d
    return scalar rej_boot2 = `w' > r(p95)

The programme simulates one dataset at a time and estimates the four different versions mentioned above and saves the information on whether $H_0$ is rejected in r(rej*) at 5%. To get 499 datasets and 499 times the information on r(rej*) I run:

simulate rej=r(rej) rej_clust=r(rej_clust) rej_boot=r(rej_boot) rej_boot2=r(rej_boot2) , reps(499) seed(1342): mysimu


Rejection rates are:

  1. 53%
  2. 24%
  3. 21%
  4. 19%


  1. The rejection rate for approach 4) is smaller than for approaches 2) and 3) in accordance with the results in Cameron [3]. However, the difference is small relative to their findings. Does someone have an idea why this is?
  2. As I read in Cameron/Trivedis "Microeconometrics Using Stata", the wild bootstrap procedure (not used here but "best" approach in their paper [2]) is for linear models only ("For linear regression, a wild bootstrap accomodates the more realistic assumptions that ..."). Is there a non-linear counterpart that may help in getting closer to the 5% rejection rate in this case?
  3. In general, is there a "state-of-the-art"-approach to handle the problem of few clusters when modelling a binary response? If so, how would the (bootstrap?) procedure look like?


In their practitioners guide Cameron and Miller point out that:

"If cluster-specific effects are present then the pairs cluster bootstrap must be adapted to account for the following complication. Suppose cluster 3 appears twice in a bootstrap resample. Then if clusters in the bootstrap resample are identified from the original cluster-identifier, the two occurrences of cluster 3 will be incorrectly treated as one large cluster rather than two distinct clusters. In Stata, the bootstrap option idcluster ensures that distinct identifiers are used in each bootstrap resample. Examples are regress y x i.id_clu, vce(boot, cluster(id_clu) idcluster(newid) reps(400) seed(10101)) and, more simply, xtreg y x, vce(boot, reps(400) seed(10101)) , as in this latter case Stata automatically accounts for this complication."

Because of this I changed the code so that each cluster-draw in the bootstrap is identified as individual cluster using the idcluster option. Also, I decided to drop the initial seed value from the program as it is also defined within the simulate command. I revised the rejection rate above accordingly.


1 Answer 1


There is an extension of the wild bootstrap called the "score bootstrap" developed by Kline and Santos (2012) (working paper here). Whereas the wild works for OLS, the score method works additionally for ML models such as logit/probit and 2SLS and GMM models. The user-written Stata command boottest can calculate p-values using the score method after an initial estimation.

  • $\begingroup$ Thanks for this valuable answer. I downloaded boottest from ssc and tried to redo the analysis using the "score bootstrap". Unfortunately, Stata (version 13.1) came up with an error message: "<istmt>: 3499 boottestStataVersion() not found" after typing "boottest, h0(1) cluster(group) bootcluster(group) ". Do you know what to do about it? $\endgroup$ Jun 30, 2016 at 11:01
  • 1
    $\begingroup$ That's a function in the provided mlib so it seems like Stata can't find the mlib. You can see the package contents here. Make sure the mlib is installed correctly and that it is found in the mlib search path (outputted from mata: mata mlib query). You can try rebuilding the mlib search list with mata: mata mlib index. Also ensure your ado-search path picks up the install folder/root (see $S_ADO). $\endgroup$ Jul 1, 2016 at 13:57
  • $\begingroup$ If I understand correctly it is about telling Stata to call a function from a library (like with library() in R). I typed mata: mata mlib index and found lboottest within the list. To see which functions are in there I typed mata: mata describe using lboottest and the function that comes closest is boottest_stata() but boottestStataVersion() was not in the list. I tried mata: mata mlib add lboottest boottestStataVersion() giving an error that the function could not be found. Have I missed something here? Does the package work for you? $\endgroup$ Jul 2, 2016 at 10:06
  • 1
    $\begingroup$ boottest.ado has some code at the top that rebuilds the mlib (via boottest.mata) if one's Stata version is different than the Stata version the mlib was created with (11.2). Maybe that rebuilding isn't working well (it does assume you install into the PLUS folder, which is usually <Stata program folder>/plus/b/). You could try debugging the rebuilding, or alternatively disable the rebuild (with a clean install it shouldn't be necessary). To display, uninstall boottest, install boottest from ssc, close Stata, edit boottest.ado to comment out lines 21-30, open Stata and try example. $\endgroup$ Jul 3, 2016 at 20:20
  • $\begingroup$ I tried to uninstall the package and encountered an error: criterion matches more than one package. Luckily enough I found your conversation with Robert Picard here which helped solving this issue. I manually deleted the package, cleaned stata.trk, reinstalled boottest from ssc and it worked. The simulation above now yields a rejection rate close to 10% when using the score bootstrap. This is indeed a large improvement. $\endgroup$ Jul 4, 2016 at 7:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.