Statalist The Stata Listserver


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

st: Lowess confidence bands


From   "Austin Nichols" <austinnichols@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   st: Lowess confidence bands
Date   Wed, 7 Mar 2007 20:41:21 -0500

Paswel Phiri Marenya --
You may want to set reps() higher than 10.  I only set it low so the
demonstration would run quickly.  You should try 1000 or 2000 or more
and make sure the results are not sensitive to the number of reps. You
should read -help simulate- and -h bsample- and all of the other
wonderful help files in the Stata directories on your machine for
details of usage.  As for statistical validity, Stas has raised
concerns (in several threads, I believe) about the bootstrap and
kernel-based estimators' convergence and how big N has to be "close"
to \infty.  He is much more knowledgable about such things than I, so
there is reason to worry...  Bradley Anderson pointed you to this
reference (I have not read it):
John Fox (2000).  Nonparametric Simple Regression: Smooothing
Scatterplots.  Quantitative Applications in the Social Sciences,
07-130.  Sage Publications.

As always with free advice, caveat emptor.

On 3/7/07, Paswel Phiri Marenya <ppm3@cornell.edu> wrote:
Dear Austin:
This is very helpful.

I have tried to use these lines and get the message:
''insufficient obervations to compute bootstrap standard errors''

 One question. My data has 445 rows, with GPS code identifiers. How do I
modify the pr eachone, bsample commands to suit this situation. Dioes
this have something to do with the number of rows? or the way x is
generated here?

Thanks
Paswel


> Paswel Phiri Marenya asked:
>> > How do I get confidence bands around my -lowess- graphs?
>
> This example may be a useful starting point
> (note I use -locpoly- rather than -lowess-, but the idea is the same):
>
> net inst st0053_3, replace
> webuse nlswork, clear
> gen work=wks_w>0 if !mi(wks_w)
> keep work age id
> save /out, replace
>
> cap pr drop _all
> pr eachone, rclass
> drop _all
> use /out
> bsample, cluster(idcode)
> gen x=_n in 20/40
> locpoly work age, gen(sm) at(x) nog
> forv a=20/40 {
>  return scalar a`a'=sm[`a']
> }
> end
>
> gen x=_n in 20/40
> locpoly work age, gen(sm) at(x) nog
> qui reg work age
> local n=e(N)
> local rownames
> local rnames
> forv a=20/40 {
>  mat a=nullmat(a) , `=sm[`a']'
>  local rownames `rnames' a`a'
>  local rnames "`rnames' a`a'=r(a`a')"
> }
> simulate `rnames', reps(10) seed(12345): eachone
> bstat, stat(a) n(`n')
> mat b=e(b)
> mat v=e(V)
> drop _all
> set obs 21
> gen coef=.
> gen ll=.
> gen ul=.
> forv i=1/21 {
>  qui replace coef=b[1,`i'] in `i'
>  qui replace ll=b[1,`i']-1.96*sqrt(v[`i',`i']) in `i'
>  qui replace ul=b[1,`i']+1.96*sqrt(v[`i',`i']) in `i'
> }
> g a=_n+19
> la var a "Age in years"
> tw rarea ll ul a || line coef a, leg(lab(1 "CI") lab(2 "Working"))
> *
> *   For searches and help try:
> *   http://www.stata.com/support/faqs/res/findit.html
> *   http://www.stata.com/support/statalist/faq
> *   http://www.ats.ucla.edu/stat/stata/
>


*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/

*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index