# Bootstrapping in Binary Response Data with Few Clusters and Within-Cluster Correlation

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

# Introduction

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^*$:

$y^*_{ig}=\beta(z_g+z_{ig})+u_{ig}$

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:

$b=\sqrt{\frac{27z_g^2}{\pi^2}}$

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

# Simulation

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
clear

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)
end


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


# Results

Rejection rates are:

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

# Questions

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?

# EDIT

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

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.

• 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? – Roberto Liebscher Jun 30 '16 at 11:01
• 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). – BeingQuisitive Jul 1 '16 at 13:57
• 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? – Roberto Liebscher Jul 2 '16 at 10:06
• 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. – BeingQuisitive Jul 3 '16 at 20:20
• 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. – Roberto Liebscher Jul 4 '16 at 7:36