RE: Re: st: ci for a standard deviation

["Seed, Paul" <paul.seed@kcl.ac.uk>]

Richard Goldstein is interested in a sensitivity analysis for a power calculation,
And in particular in variation in the SD (estimated from a sample).

It is important to use a method that allows for a non-symmetrical
Confidence Interval, as the sampling distribution for as SD is positively skew, 
with the upper limit (which matters most) much further from the estimate 
than the lower limit.  The method based on the chi-sq distribution catches this, 
the ordinary bootstrap (N), which assumes a Normally distributed sampling distribution, 
 does not.

However, the BC (bias-corrected) and in particular the BCA (bias corrected & accelerated)
 bootstrap does.  The percentile (P) bootstrap is overoptimistic, because it ignores resampling 
bias but it is at least asymmetric.

The following code demonstrates.  I have increased the replications for greater stability.
I have also introduced a second variable (mpg_10) with much greater skew to exaggerate the problem.
Of course, if Richard's variables are as badly behaved as mpg_10 (skew 7.3, kurtosis 58), 
methods based on the Normal distribution (t-test, linear regression, ANOVA etc. will not work 
well for the power calculation or for the final analysis.


*************CODE BEGINS*************
	sysuse auto, clear
	set seed 49891
	
	su mpg, detail
	local lb = r(sd)*sqrt((r(N)-1)/invchi2((r(N)-1), 1-(.05/2)))
	local ub = r(sd)*sqrt((r(N)-1)/invchi2((r(N)-1), .05/2))

	qui bootstrap r(sd) , reps(1000)  bca : sum mpg 
	estat bootstrap, all
	di "Normal formula: " %5.4f `lb' " to " %5.4f `ub'
	set seed 49891
	cap gen mpg_10 = (mpg/20)^10
	su mpg_10, detail
	local lb = r(sd)*sqrt((r(N)-1)/invchi2((r(N)-1), 1-(.05/2)))
	local ub = r(sd)*sqrt((r(N)-1)/invchi2((r(N)-1), .05/2))

	qui bootstrap r(sd),  reps(1000) bca : sum mpg_10 
	estat bootstrap, all
	di "Normal formula: " %5.4f `lb' " to " %5.4f `ub'
**************CODE ENDS**************

> Date: Wed, 26 Feb 2014 08:35:51 -0500
> From: Richard Goldstein <richgold@ix.netcom.com>
> Subject: Re: st: ci for a standard deviation
> 
> Steve,
> 
> thanks for this
> 
> I note, however, that (1) the bootstrap CI is narrower than the normal
> CI; (2) since I am using this as a "sensitivity" analysis for a sample
> size estimation, I prefer the broader CI :-)
> 
> Best,
> Rich
> 
> On 2/25/14, 7:11 PM, Steve Samuels wrote:
> > Alan's advice is valid only for normally distributed data. For a general
> approach, use
> > the bootstrap or jackknife:
> >
> > *************CODE BEGINS*************
> > sysuse auto, clear
> > set seed 49891
> > bootstrap r(sd) : sum mpg
> > estat bootstrap, all
> > **************CODE ENDS**************
> >
> > Steve
> > sjsamuels@gmail.com
> >
> >
> > On Feb 25, 2014, at 10:08 AM, Richard Goldstein <richgold@ix.netcom.com>
> wrote:
> >
> > Alan,
> >
> > thank you
> >
> > Best,
> > Rich
> >
> > On 2/25/14, 9:57 AM, Alan Neustadtl wrote:
> >> Rich,
> >>
> >>
> >> How about:
> >> Assuming sd=15 and n=1,000
> >>
> >> di 15*sqrt((1000-1)/invchi2((1000-1), .05/2))
> >> di 15*sqrt((1000-1)/invchi2((1000-1), 1-(.05/2)))
> >>
> >> Best,
> >> Alan
> >>
> >> On Tue, Feb 25, 2014 at 9:49 AM, Richard Goldstein
> >> <richgold@ix.netcom.com> wrote:
> >>> Hi,
> >>>
> >>> I want to calculate a CI for a standard deviation; anyone have Stata
> >>> code for this?
> >>>
> >>> And, no, it is not the CI given by sdtest (that is a CI for the mean for
> >>> some reason)
> >>>
> >>> Rich


Paul T Seed, Senior Lecturer in Medical Statistics, 
Division of Women's Health, King's College London
Women's Health Academic Centre, King's Health Partners 
(+44) (0) 20 7188 3642.


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