Bookmark and Share

Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down on April 23, and its replacement, statalist.org is already up and running.


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: proportion confidence intervals <0 and >1 in MI SVY


From   Austin Nichols <austinnichols@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: proportion confidence intervals <0 and >1 in MI SVY
Date   Tue, 11 May 2010 14:42:57 -0400

Steve--

I think -nlcom- is harder to integrate with the poster's mi/svy setup, but try:

clear all
program ilogit
version 11
tempname v b d
mat `v'=e(V_mi)
mat `b'=e(b_mi)
mat `d'=e(df_mi)
forv i=1/`=colsof(`b')' {
loc l: word `i' of `c'
loc df=`d'[1,`i']
loc pr=`b'[1,`i']
loc xb=ln(`pr'/(1-`pr'))
loc se=sqrt(`v'[`i',`i'])/`pr'/(1-`pr')
loc bound=invttail(`df',.025)*`se'
loc ll = exp(`xb'-`bound')
loc ul = exp(`xb'+`bound')
loc ll = `ll'/(1+`ll')
loc ul = `ul'/(1+`ul')
di as res `l',`pr' " (" `ll' "," `ul' ")"
}
end
webuse mheart1s20
g acat=age-mod(age,20)
g w=100+_mi_id
mi svyset [pw=w]
mi estimate:svy:proportion acat
ilogit acat
mi extract 0, clear
svy:tab acat, ci

* or set up as estimation cmd to use e.g. -estout- (on SSC)
clear all
program misvypr, eclass
version 11
if replay() {
 syntax [anything]
 eret di
 }
else {
syntax varname [, f(string) Percent ]
if "`f'"=="" loc f `f' "%7.4f"
if "`percent'"!="" loc m 100*
qui mi estimate:svy:proportion `varlist'
tempname v b d o o1 o2 o3
mat `v'=e(V_mi)
mat `b'=e(b_mi)
mat `d'=e(df_mi)
loc c: colnames `b'
forv i=1/`=colsof(`b')' {
loc l: word `i' of `c'
loc df=`d'[1,`i']
loc pr=`b'[1,`i']
loc xb=ln(`pr'/(1-`pr'))
loc se=sqrt(`v'[`i',`i'])/`pr'/(1-`pr')
loc bound=invttail(`df',.025)*`se'
loc ll = exp(`xb'-`bound')
loc ul = exp(`xb'+`bound')
loc ll = `ll'/(1+`ll')
loc ul = `ul'/(1+`ul')
di as res `l',`f' `m'`pr' " (" `f' `m'`ll' "," `f' `m'`ul' ")"
mat `o1'=nullmat(`o1'),`m'`pr'
mat `o2'=nullmat(`o2'),`m'`ll'
mat `o3'=nullmat(`o3'),`m'`ul'
loc c1 `c1' pr:`l'
loc c2 `c2' lb:`l'
loc c3 `c3' ub:`l'
}
mat colnames `o1'=`c1'
mat colnames `o2'=`c2'
mat colnames `o3'=`c3'
mat `o'=`o1',`o2',`o3'
eret post `o'
eret local depvar "`varlist'"
eret local cmd "misvypr"
eret local properties "b"
 }
end
webuse mheart1s20
g acat=age-mod(age,20)
g w=100+_mi_id
mi svyset [pw=w]
misvypr acat
est sto pr
esttab pr, unstack
mi extract 0, clear
svy:tab acat, ci


On Tue, May 11, 2010 at 11:39 AM, Steve Samuels <sjsamuels@gmail.com> wrote:
> Here's the corrected version. I also added the display of all the
> proportions from -svy: prop-.   Sorry for the error.
>
> Steve
>
> ************CODE BEGINS***********
> capture program drop _all
> program backlogit
> nlcom log($parm/(1-$parm))
> scalar lparm = el(r(b),1,1)
> scalar se  =  sqrt(el(r(V),1,1))
> scalar bound = invttail(e(df_r),.025)*se
> scalar ll = exp(lparm -bound)
> scalar ul = exp(lparm +bound)
> scalar ll = ll/(1+ll)
> scalar ul = ul/(1+ul)
> di %9.4f $parm   %9.4f ll  %9.4f ul
> end
>
> sysuse auto,clear
> svyset _n
> svy: prop rep78
> levelsof rep78, local(levels)
> foreach x of local  levels{
> global parm _b[`x']
> backlogit
> }
> ********CODE ENDS*********************
>
>
> On Tue, May 11, 2010 at 11:18 AM, Steve Samuels <sjsamuels@gmail.com> wrote:
>> There's an error someplace in this program, so please disregard for
>> the time being.
>>
>> Steve
>>
>>> On Tue, May 11, 2010 at 8:44 AM, Paul Shattuck <pshattuck@wustl.edu> wrote:
>>>> Hi,
>>>>
>>>> I am analyzing imputed survey data using the MI SVY estimation commands. I am estimating proportions and requesting confidence intervals (syntax below). However, many of the confidence intervals are falling outside of 0,1 bounds. The tab confidence intervals available in plain SVY (without MI) are generated using a logistic transform that prevents this problem. However, the tab procedure is not available in MI estimation.
>>>>
>>>> Can anyone suggest a work-around that will let me obtain proportion confidence intervals in MI SVY estimation that remain within the 0,1 bounds?
>>>>
>>>> mi estimate: svy, subpop(if var1==1): proportion var2
>>>>
>>>> Thank you,
>>>>
>>>> Paul

*
*   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/


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