Re: st: Bootstrapping with 'bca' option will not give BCa confidence interval

[Stas Kolenikov <skolenik@gmail.com>]

I would suspect a fundamental difficulty here. BCa CIs rely on
jackknife to compute the acceleration constant -a- (see Methods and
Formulas in [R] bootstrap), but jackknife may fail to deliver enough
(or in fact any) variability in the sample median(s), producing a zero
in the denominator of -a-, and hence a missing value overall. I am not
sure what regularity conditions are needed for consistency of the BCa
CIs; it may well be that the difference in two medians is not smooth
enough to satisfy these conditions. Never assume that the bootstrap
works in each and every case; it does for the means of i.i.d. data
with two moments, but that's about the complete list of the situations
where the bootstrap is guaranteed to work.

Note that your -gettoken- way of parsing the input variables
completely obscures what the program does. You could've used -args-,
or better yet -syntax varlist(numeric min=2 max=2), or two required
options like -outcome(varname) group(varname)-. While what you have is
a technically valid solution, it is difficult to digest and handle.

On Mon, Apr 16, 2012 at 12:04 PM, Philip Jones <pjones8@uwo.ca> wrote:
> Hello,
>
> I am trying to use bootstrapping to calculate the difference between two
> medians, as in the following program. Unfortunately, even though the code
> works correctly, I cannot seem to get BCa confidence intervals, although I can
> get every other type of confidence interval.
>
> Note that this problem seems to be the same as experienced in the following
> Statalist discussion from 2005
> (http://www.stata.com/statalist/archive/2005-12/msg00094.html)
> but it doesn't seem to have been resolved as far as I can see.
>
> My code is:
>
> ------------
> capture program drop mydiff3
> program mydiff3, rclass
>
> version 12
> gettoken number 0 : 0, parse(" ,")
> gettoken grp 0 : 0, parse(" ,")
>
>
>     quietly summarize `number' if `grp'==0, detail
>     local median_group_hi = r(p50)
>     quietly summarize `number' if `grp'==1, detail
>     local median_group_lo = r(p50)
>     return scalar difference = `median_group_lo'-`median_group_hi'
>
> end
>
>
> bootstrap r(difference), bca reps(1000) seed(12345) nodots saving(bs, replace) ///
> nowarn: mydiff3 dose group
> estat bootstrap, all
> ------------
>
> This is used on a dataset with real numbers as the "number" variable, and an
> indicator (0/1) variable ("grp").


-- 
Stas Kolenikov, also found at http://stas.kolenikov.name
Small print: I use this email account for mailing lists only.

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/