/*************************************************************************/ /* DASP : Distributive Analysis Stata Package (version 1.5) */ /*************************************************************************/ /* Conceived and programmed by Dr. Araar Abdelkrim (2008) */ /* Universite Laval, Quebec, Canada */ /* email : aabd@ecn.ulaval.ca */ /* Phone : 1 418 656 7507 */ /*************************************************************************/ /* module : difgt */ /*************************************************************************/ ** EDITED BY SEAN HIGGINS: ** 07-03-2016 Replace inefficient while loop through observations ** Formatting and remove table of results ** Calculate and return p-value of test of difference = 0 #delim ; /*****************************************************/ /* Density function : ici fw=Hweight*Hsize */ /*****************************************************/ cap program drop dens; program define dens, rclass; args fw x xval; qui su `x' [aw=`fw'], detail; local tmp = (`r(p75)'-`r(p25)')/1.34; local tmp = (`tmp'<`r(sd)')*`tmp'+(`tmp'>=`r(sd)')*`r(sd)' ; local h = 0.9*`tmp'*_N^(-1.0/5.0); tempvar s1 s2; gen `s1' = sum( `fw' *exp(-0.5* ( ((`xval'-`x')/`h')^2 ) )); gen `s2' = sum( `fw' ); return scalar den = `s1'[_N]/( `h'* sqrt(2*c(pi)) * `s2'[_N] ); end; /***************************************/ /* Non parametric regression */ /***************************************/ cap program drop npe; program define npe, rclass; args fw x y xval; qui su `x' [aw=`fw'], detail; local tmp = (`r(p75)'-`r(p25)')/1.34 ; local tmp = (`tmp'<`r(sd)')*`tmp'+(`tmp'>=`r(sd)')*`r(sd)'; local h = 0.9*`tmp'*_N^(-1.0/5.0); tempvar s1 s2; gen `s1' = sum( `fw'*exp(-0.5* ( ((`xval'-`x')/`h')^2 ) ) *`y'); gen `s2' = sum( `fw'*exp(-0.5* ( ((`xval'-`x')/`h')^2 ) )); return scalar nes = `s1'[_N]/`s2'[_N]; end; /***************************************/ /* Quantile */ /***************************************/ cap program drop quant; program define quant, rclass; args fw yyy xval; preserve; sort `yyy'; qui cap drop if `yyy'>=. | `fw'>=.; tempvar ww qp pc; qui gen `ww'=sum(`fw'); qui gen `pc'=`ww'/`ww'[_N]; qui gen `qp' = `yyy' ; qui sum `yyy' [aw=`fw']; /* added by Sean to replace inefficient while loop through observations */ tempvar ordering ; gen `ordering' = _n; count if (`pc'<`xval'); if r(N) { ; // i.e. positive number of obs with `pc'<`xval' ; summ `ordering' if (`pc'<`xval'); local ar = r(max); } ; else { ; local ar = 0; }; local i = `ar' + 1; if (`i'> 1) local qnt=`qp'[`ar']+((`qp'[`i']-`qp'[`ar'])/(`pc'[`i']-`pc'[`ar']))*(`pc'[`i']-`pc'[`ar']); if (`i'==1) local qnt=(max(0,`qp'[`i'])/(`pc'[`i']))*(`pc'[`i']); return scalar qnt = `qnt'; restore; end; /***************************************/ /* FGT */ /***************************************/ capture program drop fgt; program define fgt, rclass; syntax varlist(min=1 max=1) [, FWeight(string) HGroup(string) GNumber(integer 1) PLine(real 0) ALpha(real 0) type(string) ]; tokenize `varlist'; tempvar fw ga hy; gen `hy' = `1'; gen `fw'=1; if ("`fweight'"~="") qui replace `fw'=`fw'*`fweight'; if ("`hgroup'" ~="") qui replace `fw'=`fw'*(`hgroup'==`gnumber'); gen `ga' = 0; if (`alpha'==0) qui replace `ga' = (`pline'>`hy'); if (`alpha'~=0) qui replace `ga' = ((`pline'-`hy'))^`alpha' if (`pline'>`hy'); qui sum `ga' [aweight= `fw']; local fgt = r(mean); if ("`type'" == "nor") local fgt = `fgt' /(`pline'^`alpha'); return scalar fgt = `fgt'; end; capture program drop ceqdifgt2d; /* note deleted ceqdifgt2 since that was only for those coming from different files */ program define ceqdifgt2d, rclass; version 9.2; syntax namelist [, FWEIGHT1(string) FWEIGHT2(string) HS1(string) HS2(string) AL(real 0) PL1(string) PL2(string) OPL1(string) PROP1(real 0.5) PERC1(real 0.4) OPL2(string) PROP2(real 0.5) PERC2(real 0.4) TYpe(string) INDex(string) CONF(string) LEVEL(real 95) ITER(real 1000) /* added this to be able to adjust iterations but it makes no difference */ ]; tokenize `namelist'; tempvar num1 num2 snum1 snum2 hsy1 hsy2; qui gen `num1'=0; qui gen `num2'=0; if ("`opl1'"=="mean") {; qui sum `1' [aw=`fweight1'], meanonly; local pl1=`prop1'*`r(mean)'; local mu = `r(mean)'; if (`al'==0) {; qui dens `fweight1' `1' `pl1'; local pprim1 = `r(den)'; }; if (`al'>=1) {; local al1=`al'-1; qui fgt `1', fweight(`fweight1') pline(`pl1') alpha(`al1'); local pprim1 = `al'*`r(fgt)'; }; }; if ("`opl1'"=="quantile") {; qui quant `fweight1' `1' `perc1'; local pl1=`prop1'*`r(qnt)'; local qnt1= `r(qnt)'; if (`al'==0){; qui dens `fweight1' `1' `pl1'; local pprim1 = `r(den)'; qui dens `fweight1' `1' `qnt1'; local fqp1 = `r(den)'; }; if (`al'>=1){; local al1=`al'-1; qui fgt `1', fweight(`fweight1') pline(`pl1') alpha(`al1'); local pprim1 = `al'*`r(fgt)'; qui dens `fweight1' `1' `qnt1'; local fqp1 = `r(den)'; }; }; fgt `1', fweight(`fweight1') pline(`pl1') alpha(`al') type(`type'); local est1 = `r(fgt)'; cap drop `num1' `snum1'; tempvar num1 snum1; qui gen `num1'=0; qui gen `snum1'=0; if (`al' == 0) qui replace `num1' = `hs1'*(`pl1'> `1'); if (`al' ~= 0) qui replace `num1' = `hs1'*(`pl1'-`1')^`al' if (`pl1'>`1'); if ((`"`opl1'"'=="") & (`"`type'"'=="nor") & `al'!=0) /* " */ qui replace `num1' = `num1'/(`pl1'^`al') if (`pl1'>`1'); if ("`opl1'"=="mean") qui replace `snum1' = `pprim1'*`hs1'*(`prop1'*`1'- `pl1'); if ("`opl1'"=="quantile") qui replace `snum1' = -`pprim1'*`hs1'* `prop1' * ((`qnt1'>`1')-`perc1')/`fqp1' ; if (`"`opl1'"'=="" ) {; // " ; qui svy: ratio (eq1: `num1'/`hs1'); cap drop matrix _vv; matrix _vv=e(V); local std1 = el(_vv,1,1)^0.5; }; if ("`opl1'"~="") {; cap drop `hsy1'; if ("`opl1'"=="mean") qui gen `hsy1' = `prop1'*`hs1'*`1'; if ("`opl1'"=="quantile") qui gen `hsy1'= -`hs1'* `prop1' *((`qnt1'>`1')-`perc1')/`fqp1' + `hs1'*`pl1'; qui svy: mean `num1' `snum1' `hsy1' `hs1'; if ("`type'"!="nor" | `al'==0 ) qui nlcom (eq1: _b[`num1']/_b[`hs1']+_b[`snum1']/_b[`hs1']), iterate(`iter'); if ("`type'"=="nor" & `al'!=0) qui nlcom (eq1: (_b[`num1']/_b[`hs1']+_b[`snum1']/_b[`hs1'])/(_b[`hsy1']/_b[`hs1'])^`al' ), iterate(`iter'); cap drop matrix _vv; matrix _vv=r(V); local std1 = el(_vv,1,1)^0.5; }; if ("`opl2'"=="mean") {; qui sum `2' [aw=`fweight2'], meanonly; local pl2=`prop2'*`r(mean)'; local mu = `r(mean)'; if (`al'==0) {; qui dens `fweight2' `2' `pl2'; local pprim2 = `r(den)'; }; if (`al'>=1) {; local al2=`al'-1; qui fgt `2', fweight(`fweight2') pline(`pl2') alpha(`al2'); local pprim2 = `al'*`r(fgt)'; }; }; if ("`opl2'"=="quantile") {; qui quant `fweight2' `2' `perc2'; local pl2=`prop2'*`r(qnt)'; local qnt2= `r(qnt)'; if (`al'==0) {; qui dens `fweight2' `2' `pl2'; local pprim2 = `r(den)'; qui dens `fweight2' `2' `qnt2'; local fqp2 = `r(den)'; }; if (`al'>=1){; local al1=`al'-1; qui fgt `2', fweight(`fweight2') pline(`pl2') alpha(`al1'); local pprim2 = `al'*`r(fgt)'; qui dens `fweight2' `2' `qnt2'; local fqp2 = `r(den)'; }; }; fgt `2', fweight(`fweight2') pline(`pl2') alpha(`al') type(`type'); local est2 = `r(fgt)'; cap drop `num2' `snum2'; tempvar num2 snum2; qui gen `num2'=0; qui gen `snum2'=0; if (`al' == 0) qui replace `num2' = `hs2'*(`pl2'> `2'); if (`al' ~= 0) qui replace `num2' = `hs2'*(`pl2'-`2')^`al' if (`pl2'>`2'); if ((`"`opl2'"'=="") & (`"`type'"'=="nor") & `al'!=0) /* " */ qui replace `num2' = `num2'/(`pl2'^`al') if (`pl2'>`2'); if ("`opl2'"=="mean") qui replace `snum2' = `pprim2'*`hs2'*(`prop2'*`2'- `pl2'); if ("`opl2'"=="quantile") qui replace `snum2' = -`pprim2'*`hs2'* `prop2' * ((`qnt2'>`2')-`perc2')/`fqp2'; if (`"`opl2'"'=="" ) {; // " ; qui svy: ratio (eq2: `num2'/`hs2'); cap drop matrix _vv; matrix _vv=e(V); local std2 = el(_vv,1,1)^0.5; }; if ("`opl2'"~="") {; cap drop `hsy2'; if ("`opl2'"=="mean") qui gen `hsy2' = `prop2'*`hs2'*`2'; if ("`opl2'"=="quantile") qui gen `hsy2'= -`hs2'* `prop2' *((`qnt2'>`2')-`perc2')/`fqp2' + `hs1'*`pl2'; qui svy: mean `num2' `snum2' `hsy2' `hs2'; if ("`type'"!="nor" | `al'==0 ) qui nlcom (eq2: _b[`num2']/_b[`hs2']+_b[`snum2']/_b[`hs2']), iterate(`iter'); if ("`type'"=="nor" & `al'!=0) qui nlcom (eq2: (_b[`num2']/_b[`hs2']+_b[`snum2']/_b[`hs2'])/(_b[`hsy2']/_b[`hs2'])^`al'), iterate(`iter'); cap drop matrix _vv; matrix _vv=r(V); local std2 = el(_vv,1,1)^0.5; }; local est3=`est2'-`est1'; local equa1=""; if (`"`opl1'"'=="") {; // " ; local equa1 = "(_b[`num1']/_b[`hs1'])"; local smean1 = "`num1' `hs1' "; }; if (`"`opl1'"'!="") {; // " ; if ("`type'"!="nor" | `al'==0) local equa1 = "(_b[`num1']/_b[`hs1']+_b[`snum1']/_b[`hs1'] )"; if ("`type'"=="nor" & `al'!=0) local equa1 = "(_b[`num1']/_b[`hs1'])/((_b[`hsy1']/_b[`hs1'])^`al')"; local smean1 = "`num1' `snum1' `hsy1' `hs1' "; }; local equa2=""; if (`"`opl2'"'=="") {; // " ; local equa2 = "(_b[`num2']/_b[`hs2'])"; local smean2 = "`num2' `hs2' "; }; if (`"`opl2'"'!="") {; // " ; if ("`type'"!="nor" | `al'==0) local equa2 = "(_b[`num2']/_b[`hs2']+_b[`snum2']/_b[`hs2'] )"; if ("`type'"=="nor" & `al'!=0) local equa2 = "(_b[`num2']/_b[`hs2'])/((_b[`hsy2']/_b[`hs2'])^`al')"; local smean2 = "`num2' `snum2' `hsy2' `hs2' "; }; qui svy: mean `smean1' `smean2'; qui nlcom (`equa1'-`equa2'), iterate(`iter'); cap drop matrix _vv; matrix _vv=r(V); local std3 = el(_vv,1,1)^0.5; qui svydes; local fr=`r(N_units)'-`r(N_strata)'; scalar _p = ttail(`fr',abs(`est3'/`std3')); /* added by Sean: calculate p-value */ local lvl=(100-`level')/100; if ("`conf'"=="ts") local lvl=`lvl'/2; local tt=invttail(`fr',`lvl'); return scalar std1 = `std1'; return scalar est1 = `est1'; return scalar lb1 = `est1' - `tt'*`std1'; return scalar ub1 = `est1' + `tt'*`std1'; return scalar pl1 = `pl1'; return scalar df = `fr'; qui count; return scalar nobs = `r(N)'; return scalar std2 = `std2'; return scalar est2 = `est2'; return scalar lb2 = `est2' - `tt'*`std2'; return scalar ub2 = `est2' + `tt'*`std2'; return scalar pl2 = `pl2'; return scalar std3 = `std3'; return scalar est3 = `est3'; return scalar lb3 = `est3' - `tt'*`std3'; return scalar ub3 = `est3' + `tt'*`std3'; /* added by Sean */ return scalar p = _p; cap scalar drop _p; end; capture program drop ceqdifgt; program define ceqdifgt, eclass; version 9.2; syntax namelist(min=2 max=2) [, FILE1(string) FILE2(string) HSize1(string) HSize2(string) ALpha(real 0) COND1(string) COND2(string) PLINE1(real 10000) OPL1(string) PROP1(real 50) PERC1(real 0.4) PLINE2(real 10000) OPL2(string) PROP2(real 50) PERC2(real 0.4) type(string) LEVEL(real 95) CONF(string) DEC(int 6) BOOT(string) NREP(string) TEST(string) ]; global indica=3; tokenize `namelist'; if ("`prop1'"!="") local prop1=`prop1'/100; if ("`prop2'"!="") local prop2=`prop2'/100; if ("`type'"=="") local type="nor"; if ("`conf'"=="") local conf="ts"; if ("`opl1'"=="median") {; local opl1 = "quantile"; local perc1 = 0.5; }; if ("`opl2'"=="median") {; local opl2 = "quantile"; local perc2 = 0.5; }; preserve; local indep = ( (`"`file1'"'=="" & `"`file2'"'=="") | (`"`file1'"'==`"`file2'"') ); // " ; // if they come from the same data file if ( (`"`file1'"'=="" & `"`file2'"'=="") | (`"`file1'"'==`"`file2'"') ) {; // " ; if (`"`file1'"'~="") use `"`file1'"', replace; // " ; tempvar cd1 ths1; qui gen `ths1'=1; if ( "`hsize1'"!="") qui replace `ths1'=`hsize1'; if ( "`cond1'"!="") {; gen `cd1'=`cond1'; qui replace `ths1'=`ths1'*`cd1'; qui sum `cd1'; if (`r(sum)'==0) {; dis as error " With condition(s) of distribution_1, the number of observations is 0."; exit; }; }; tempvar cd2 ths2; qui gen `ths2'=1; if ( "`hsize2'"!="") qui replace `ths2'=`hsize2'; if ( "`cond2'"!="") {; gen `cd2'=`cond2'; qui replace `ths2'=`ths2'*`cd2'; qui sum `cd2i'; if (`r(sum)'==0) {; dis as error " With condition(s) of distribution_2, the number of observations is 0."; exit; }; }; qui svyset ; if ( "`r(settings)'"==", clear") qui svyset _n, vce(linearized); local hweight=""; cap qui svy: total `1'; local hweight=`"`e(wvar)'"'; // " ; cap ereturn clear; tempvar fw1 fw2; gen `fw1'=`ths1'; gen `fw2'=`ths2'; if ("`hweight'"!="") qui replace `fw1'=`fw1'*`hweight'; if ("`hweight'"!="") qui replace `fw2'=`fw2'*`hweight'; ceqdifgt2d `1' `2' , fweight1(`fw1') fweight2(`fw2') hs1(`ths1') hs2(`ths2') pl1(`pline1') pl2(`pline2') opl1(`opl1') prop1(`prop1') perc1(`perc1') opl2(`opl2') prop2(`prop2') perc2(`perc2') al(`alpha') type(`type') level(`level') conf(`conf') ; matrix _res_d1 =(`r(est1)',`r(std1)',`r(lb1)',`r(ub1)' ,`r(pl1)'); matrix _res_d2 =(`r(est2)',`r(std2)',`r(lb2)',`r(ub2)' ,`r(pl2)'); matrix _res_di =(`r(est3)',`r(std3)',`r(lb3)',`r(ub3)'); scalar _p = r(p); /* _dasp_dif_table_ifgt `1' `1', name1("`1'") name2("`2'") pl1(`r(pl1)') pl2(`r(pl2)') m1(`r(est1)') m2(`r(est2)') s1(`r(std1)') s2(`r(std2)') df1(`r(df)') df2(`r(df)') dif(`r(est3)') sdif(`r(std3)') level(`level') conf(`conf') indep(`indep') test(`test'); */ }; // if they come from different data files; not applicable (deleted) restore; ereturn matrix d1 = _res_d1; ereturn matrix d2 = _res_d2; ereturn matrix di = _res_di; ereturn scalar p = _p; cap matrix drop _res_d1; cap matrix drop _res_d2; cap matrix drop _res_di; cap scalar drop _p; end;