*! Version 1.0, Dirk Enzmann (20-Aug-2026) *! *! Design-based (survey) H/SE/CI for Loevinger's H (e.g. to use with pweight) program define loevh2_svy, rclass byable(recall) version 16.0 local _origframe = c(frame) local _origchanged = c(changed) syntax varlist(min=2 max=2 numeric) [if] [in] /// [, Table Level(cilevel) Compare meta(string asis) _svytest] // Tokenize `varlist' into clean `1'/`2' variable-name macros right // away, so that the missing-byvar skip-forwarding branch below (which // returns r(var1)/r(var2) using `1'/`2' BEFORE the later, pre-existing // "tokenize `varlist'" call further down the program) never picks up // a stray trailing comma or other leftover token from the raw command // line's positional-macro parse. The later tokenize call (kept as-is, // for `v1'/`v2') is now redundant but harmless. tokenize `varlist' // See: https://www.statalist.org/forums/forum/general-stata-discussion/general/152134 capture which fre if ( _rc ) { capture window stopbox rusure "This program requires -fre-. Do you want to install -fre-?" if ( _rc ) exit 111 ssc install fre } if "`compare'"=="compare" & !_by() { di as error "option compare requires {help by}" exit 498 } // meta() no longer requires -compare- (nor -by-): see loevh2.ado's // own changelog note for the full rationale (identical here). // Require an active svyset design with at least one of pweight/ // iweight, cluster(), or strata() actually set -- i.e., a genuine // complex-survey design, not merely an "svyset _n" no-op. If none of // these three is set, loevh2_svy has nothing to offer beyond loevh2/ // loevh2_boot (which are faster and better documented for the // unweighted/simple-pweight case), so it errors out and points the // user there instead. // // The svyset design's own descriptive fields (weight type/expression, // cluster/strata variable names) are captured here too, right after // the design is confirmed active, for later use (regardless of // -compare-) as purely informational string columns whenever meta() // is specified -- see loevh2_meta_save.ado's header for the full // rationale on why the raw data/design itself is never replayed. local _svy_wexp "" local _svy_fpc "" if "`_svytest'" == "" { capture svyset if _rc { di as error "no svyset design is active -- loevh2_svy requires " /// "{help svyset} with pweight, cluster(), and/or " /// "strata() to be set first. For unweighted data, " /// "or pweight without clustering/stratification, " /// "use {help loevh2} or {help loevh2_boot} instead." exit 119 } local _svy_wtype "`r(wtype)'" local _svy_wexp "`r(wexp)'" local _svy_clust "`r(su1)'" local _svy_strata "`r(strata1)'" local _svy_fpc "`r(fpc1)'" local _svy_has_pw = inlist("`_svy_wtype'", "pweight", "iweight") local _svy_has_cl = ("`_svy_clust'" != "") local _svy_has_st = ("`_svy_strata'" != "") if !`_svy_has_pw' & !`_svy_has_cl' & !`_svy_has_st' { di as error "the active svyset design has no pweight/iweight, " /// "no cluster(), and no strata() -- loevh2_svy " /// "requires at least one of these to be genuinely " /// "set. For unweighted, non-clustered, non-" /// "stratified data, use {help loevh2} or " /// "{help loevh2_boot} instead." exit 119 } } marksample touse // Build a second sample marker, `tousegrp', that reflects only the // -in- restriction (if any), but NOT the user's -if- condition, // exactly as loevh2/loevh2_boot do, so a by-group's missing by-value // can be detected independent of any -if- that might otherwise zero // out `touse' for an entire missing-value by-group. tempvar tousegrp mark `tousegrp' `in' markout `tousegrp' if _by() { local _byvars_miss = 0 foreach _bv of local _byvars { qui count if `tousegrp' & missing(`_bv') if r(N) > 0 local _byvars_miss = 1 } // Also proactively check whether `touse', though not entirely // zero, still leaves ZERO observations with non-missing analysis // variables (v1/v2) for this by-group -- exactly mirroring // loevh2.ado's own analogous fix (see its changelog/comments for // the full rationale). This prevents ever reaching the binary- // variable-check/tabulate/matrix logic below with zero usable // observations for this by-group. local _byvars_empty = 0 if !`_byvars_miss' { qui count if `touse' if r(N) == 0 local _byvars_empty = 1 else { tokenize `varlist' qui count if `touse' & !missing(`1') & !missing(`2') if r(N) == 0 local _byvars_empty = 1 } } if `_byvars_miss' | `_byvars_empty' { di as txt _n "Note: by-group with missing value of " /// `"`_byvars'"' " (or no observations selected " /// "by if/in) skipped (H is not estimated for " /// "missing values of the by-variable)." return local error "missing_byvar" // Retrieve and forward the last VALID by-group's results // (r(lastgroup), r(loevh), r(se), r(N), r(overlap), etc.) from // the "svyressubsamp" frame, UNCONDITIONALLY (regardless of // -compare-) whenever this is the true last by-group of the // whole by-sequence -- exactly mirroring loevh2.ado's own // equivalent behavior (see its "ressubsamp" frame handling), // so that r() is never left populated with only // r(error)="missing_byvar" and nothing else just because the // by-sequence happens to end on a missing/excluded group. The // "svyressubsamp" frame itself is now created/posted for every // valid by-group regardless of -compare- (see the main-flow // posting block further below), so it is always available // here to retrieve from, whether or not -compare- was // requested. Only the loevh2_svy_compare() call and its // compare-specific returns (Hbar, chi2, df, p_chi2, H_SE_N) // remain conditional on -compare-. // // NOTE: meta() is deliberately NOT triggered in this skip- // forwarding branch, even though the forwarded results // originate from the last VALID group -- that valid group // already triggered its own meta() save when IT was // processed, so saving again here would create a duplicate // row in the meta() output file. if _bylastcall() { capture confirm frame svyressubsamp if _rc==0 { if "`compare'"=="compare" { di _n _n as txt "Categories of by: variable(s):" fre `_byvars' `if' `in' } frame change svyressubsamp lab val sample sample local _last_n = _N if `_last_n' > 0 { local _lg_h = h[`_last_n'] local _lg_se = se[`_last_n'] local _lg_N = n[`_last_n'] local _lg_sval = sample[`_last_n'] local _lg_lab : label (sample) `_lg_sval' local _lg_ov = ov[`_last_n'] local _lg_ov_se = ov_se[`_last_n'] local _lg_ov_lb = ov_lb[`_last_n'] local _lg_ov_ub = ov_ub[`_last_n'] local _lg_sample = sample[`_last_n'] } if "`compare'"=="compare" { loevh2_svy_compare local pchi2 = `loevh2_pchi2' local df = `loevh2_df' local chi2 = `loevh2_chi2' local Hbar = `loevh2_Hbar' local Hbar_se = `loevh2_Hbar_se' local Hbar_N = `loevh2_Hbar_N' local HSEN = "`loevh2_HSEN'" if `"`meta'"' != `""' & "`HSEN'"!="" { local HSENLAB = `"`loevh2_HSEN_labels'"' loevh2_meta_save, metaspec(`meta') /// source(loevh2_svy) suffix("_svy") /// setype("svy, nlcom (delta method)") /// hsen(`HSEN') labels(`"`HSENLAB'"') /// svywtype(`"`_svy_wtype'"') svywexp(`"`_svy_wexp'"') /// svycluster(`"`_svy_clust'"') svystrata(`"`_svy_strata'"') /// svyfpc(`"`_svy_fpc'"') /// origframe(`"`_origframe'"') origchanged(`_origchanged') } } frame change `_origframe' frame drop svyressubsamp if `_last_n' > 0 { tempname _crit scalar `_crit' = abs(invnormal(`=(1-`level'/100)/2')) return scalar loevh = `_lg_h' return scalar se = `_lg_se' return scalar lb = `_lg_h' - scalar(`_crit')*`_lg_se' return scalar ub = `_lg_h' + scalar(`_crit')*`_lg_se' return scalar level = `level' return scalar N = `_lg_N' return scalar overlap = `_lg_ov' return scalar se_overlap = `_lg_ov_se' return scalar lb_overlap = `_lg_ov_lb' return scalar ub_overlap = `_lg_ov_ub' } if "`compare'"=="compare" { return scalar Hbar = `Hbar' return scalar Hbar_se = `Hbar_se' return scalar chi2 = `chi2' return scalar df = `df' return scalar p_chi2 = `pchi2' return scalar Hbar_N = `Hbar_N' if "`HSEN'" != "" return matrix H_SE_N = `HSEN' } if `_last_n' > 0 { return local var1 "`1'" return local var2 "`2'" return local se_type "svy, nlcom (delta method)" return local group "`_byvars'" return local lastgroup "`_lg_lab'" } } } exit } } tempname crit scalar `crit' = abs(invnormal(`=(1-`level'/100)/2')) // Check that variables are binary (0/1) foreach var of varlist `varlist' { qui levelsof `var' if `touse', local(levels) if `: word count `levels'' != 2 { di as error "Error: `var' must be a binary (0/1) variable" exit 450 } foreach val of local levels { if !inlist(`val', 0, 1) { di as error "Error: `var' must contain only values 0 and 1" exit 450 } } } tokenize `varlist' local v1 "`1'" local v2 "`2'" // Structural pre-check for a degenerate (or non-2x2) table, using a // plain, non-svy tabulate purely as a diagnostic -- never used for // the actual H/SE computation, which always comes from svy: mean + // nlcom below. This mirrors loevh2's own degenerate-table guard // exactly (any cell or margin equal to zero makes H undefined). // // CRITICAL: _T is a fixed (non-tempname) matrix name, which // persists across successive by-group calls (byable(recall) re-runs // this whole program once per group) unless explicitly cleared. If // this by-group's touse-restricted subsample yields ZERO // observations to tabulate, "tabulate ... matcell(_T)" silently // returns rc=0 WITHOUT creating or modifying _T at all, leaving it // holding the STALE matrix from the previous by-group. Dropping _T // first (exactly mirroring loevh2.ado's own analogous fix for its // T matrix) ensures a genuine "no observations" case leaves _T // undefined, so rowsof(_T)/colsof(_T) below correctly return 0, // correctly tripping the `wrongshape' guard instead of silently // reusing a stale previous-group table. capture matrix drop _T qui tabulate `v1' `v2' if `touse', matcell(_T) // _T may not exist at all here (rather than merely having the wrong // shape) whenever the touse-restricted tabulate above genuinely // matched zero observations -- capture matrix drop _T just above // guarantees this is a clean "not found" rather than a stale // leftover matrix, but rowsof()/colsof() would themselves raise // r(111) on a nonexistent matrix, so its existence must be checked // explicitly first and treated as an (rowsof=0, colsof=0) wrong- // shape case (mirrors loevh2.ado's own analogous fix). capture confirm matrix _T local _Texists = (_rc==0) local wrongshape = !`_Texists' | (rowsof(_T)!=2 | colsof(_T)!=2) if `wrongshape' { tempname _Ndeg scalar `_Ndeg' = 0 if `_Texists' { forvalues _r=1/`=rowsof(_T)' { forvalues _c=1/`=colsof(_T)' { scalar `_Ndeg' = `_Ndeg' + _T[`_r',`_c'] } } } di as error "Warning: degenerate 2x2 table (a variable has " /// "only one distinct value in this (sub)sample) " /// "-- H and its SE/CI are undefined; returning " /// "missing values" return scalar loevh = . return scalar se = . return scalar lb = . return scalar ub = . return scalar N = `_Ndeg' return local error "degenerate" exit } tempname _a2 _b2 _c2 _d2 _e2 _f2 scalar `_a2' = _T[2,2] scalar `_b2' = max(_T[1,2], _T[2,1]) scalar `_c2' = min(_T[1,2], _T[2,1]) scalar `_d2' = _T[1,1] scalar `_e2' = `_a2' + `_b2' scalar `_f2' = `_a2' + `_c2' tempname _Ndeg scalar `_Ndeg' = `_a2' + `_b2' + `_c2' + `_d2' local degenerate = (`_a2'==0 | `_b2'==0 | `_c2'==0 | `_d2'==0 | /// `_e2'==0 | `_f2'==0 | (`_Ndeg'-`_e2')==0 | /// (`_Ndeg'-`_f2')==0) if `degenerate' { di as error "Warning: degenerate 2x2 table (a zero cell or " /// "margin) -- H and its SE/CI are undefined; " /// "returning missing values" return scalar loevh = . return scalar se = . return scalar lb = . return scalar ub = . return scalar N = `_Ndeg' return local error "degenerate" exit } // Build the four joint-cell-membership indicators for the FULL // svyset-declared dataset in memory, so that svy, subpop(): mean can // use the full design information (PSUs/strata) while restricting // point estimation to the current touse subpopulation. Outside // touse==1 (excluded by if/in, or a different by-group, or missing // var1/var2 -- marksample already excludes missing values), the // indicators are set to 0, following the standard svy subpop() // convention (never left missing). // // NOTE: the branch-selection ("swap") decision itself is NOT made // here (see the block AFTER svy: mean below for why -- in short: it // must be based on the WEIGHTED margins, not on the unweighted `_T' // matrix, which is only ever used above for the structural // degenerate-table diagnostic). tempvar raw11 raw10 raw01 raw00 qui gen byte `raw11' = 0 qui gen byte `raw10' = 0 qui gen byte `raw01' = 0 qui gen byte `raw00' = 0 qui replace `raw11' = 1 if `touse' & `v1'==1 & `v2'==1 qui replace `raw10' = 1 if `touse' & `v1'==1 & `v2'==0 qui replace `raw01' = 1 if `touse' & `v1'==0 & `v2'==1 qui replace `raw00' = 1 if `touse' & `v1'==0 & `v2'==0 // svy, subpop(): mean on the full design, restricted to the current // touse subpopulation -- the statistically correct way to handle // if/in/by-group restriction under a complex survey design (naively // filtering with "svy: mean ... if ..." can bias variance estimation // by silently dropping PSU/strata information; subpop() retains the // full design while estimating only for the subpopulation). capture qui svy, subpop(`touse'): mean `raw11' `raw10' `raw01' `raw00' if _rc { di as error "svy: mean failed (rc=" _rc ")" return scalar loevh = . return scalar se = . return scalar lb = . return scalar ub = . return scalar N = `_Ndeg' return local error "svy_failed" exit } // Save the design/header e() scalars now, immediately after svy: // mean, before nlcom's own "post" (below) replaces e() with the // nlcom estimation results -- needed for the cross-tab display's // header block. tempname _e_Nstrata _e_Npsu _e_N _e_Nsub _e_Npop _e_Nsubpop _e_dfr scalar `_e_Nstrata' = e(N_strata) scalar `_e_Npsu' = e(N_psu) scalar `_e_N' = e(N) scalar `_e_Nsub' = e(N_sub) scalar `_e_Npop' = e(N_pop) scalar `_e_Nsubpop' = e(N_subpop) scalar `_e_dfr' = e(df_r) // Cross-tab cell proportions (x100, as percentages) and their // design-based SEs, straight from the svy: mean estimates just // computed -- always in the fixed v1 x v2 orientation (never // swapped), for the display table and r(overlap)/r(lb_overlap)/ // r(ub_overlap) below. tempname _p11 _p10 _p01 _p00 tempname _se11 _se10 _se01 _se00 scalar `_p11' = _b[`raw11'] scalar `_p10' = _b[`raw10'] scalar `_p01' = _b[`raw01'] scalar `_p00' = _b[`raw00'] scalar `_se11' = _se[`raw11'] scalar `_se10' = _se[`raw10'] scalar `_se01' = _se[`raw01'] scalar `_se00' = _se[`raw00'] // Branch-selection ("swap") decision -- WEIGHTED margins. // // The branch-selection logic must match loevh2's own Mata branch // EXACTLY -- and that branch condition is NOT based on comparing the // raw off-diagonal cell counts, but on comparing the two MARGINS: // loevh2.ado's Mata block branches on colsum(T)[1,2] (= P(var2=1) // margin) versus rowsum(T)[2,1] (= P(var1=1) margin): // // if colsum(T)[1,2] < rowsum(T)[2,1] // i.e. p2 < p1 // H uses T[1,2] (var1=0,var2=1) in the numerator // else // i.e. p1 <= p2 // H uses T[2,1] (var1=1,var2=0) in the numerator // // i.e. the numerator ("cell10" role below) is always the off- // diagonal cell belonging to the variable with the SMALLER marginal // probability of being 1. // // CRITICAL FIX (version 1.2): this decision MUST be made using the // same (WEIGHTED, design-based) margins that H itself is computed // from -- i.e. `_b[raw10]'+`_b[raw11]' (weighted p1) versus // `_b[raw01]'+`_b[raw11]' (weighted p2), taken from the svy: mean // estimates just obtained above -- NOT the unweighted `_T' matrix // used earlier purely for the structural degenerate-table // diagnostic. An earlier version of this program computed the swap // decision from the unweighted `_T' margins before svy: mean was // even run. Whenever weighting shifts the ordering of the two // marginal probabilities relative to their UNWEIGHTED ordering -- // which can and does happen whenever the two margins are close to a // tie, exactly the boundary condition where the branch choice // matters most -- the unweighted-margin-based decision picks the // WRONG branch (equivalently, the wrong off-diagonal cell/margin // pair) for the weighted nlcom H formula below, silently producing // an incorrect H, SE, and CI with no error or warning. This was // discovered via the ISRD dataset (test_isrd_bh.do/.log): unweighted // margins gave p1margin=304 > p2margin=299 (swap=1), while the // WEIGHTED margins gave weighted-p1=.20725 < weighted-p2=.20935 // (swap=0) -- the two disagree, and using the unweighted decision // produced H=.34736 instead of the correct H=.35183 (matching an // independent pweight-based bootstrap recomputation and manual // verification to 5 decimal places). See // verify_h_formula.py/verify_h_formula2.py/verify_h_formula3.py for // the full numerical diagnosis. local _p1margin_w = `_p10' + `_p11' local _p2margin_w = `_p01' + `_p11' local _swap = (`_p1margin_w' > `_p2margin_w') local branch = cond(`_swap', "pi01", "pi10") // nlcom's numerator cell reference (cell10 role): _b[raw10] if // p1<=p2 (loevh2's "else" branch), _b[raw01] if p1>p2 (loevh2's // "if" branch) -- matching loevh2's own margin-based swap exactly, // now correctly computed from the WEIGHTED margins above. local _cell10ref = cond(`_swap', "`raw01'", "`raw10'") local _cell01ref = cond(`_swap', "`raw10'", "`raw01'") // r(overlap): the "both = 1" (v1==1 & v2==1) proportion, i.e. raw11. // Point estimate and delta-method SE come directly from the svy: // mean estimates already computed above (a direct linear parameter, // no nlcom needed for these two). The CI, however, is obtained via // nlcom's delta method applied to logit(raw11) and back-transformed // with invlogit() -- an asymmetric, boundary-respecting CI, which is // the natural and well-behaved choice for a bounded proportion in // (0,1) (unlike H itself, which can be negative and for which an // earlier Monte Carlo simulation found the analogous logit-CI // construction to badly under-cover for small H -- see the removed- // `logit'-option note at the top of this file; that finding does // NOT apply here, since raw11 is a genuine bounded proportion, the // textbook use case for a logit-transformed Wald CI). This nlcom // call is done FIRST (before H's own nlcom below), with its results // captured into plain scalars immediately, so that H's subsequent // nlcom ..., post call (which replaces the currently posted // estimation results) does not need to be undone or re-run. tempname overlap overlap_se overlap_lb overlap_ub scalar `overlap' = `_p11' scalar `overlap_se' = `_se11' // NOTE: no `post' here -- nlcom's r(b)/r(V) are populated regardless // of `post', and using `post' would replace the currently-posted // svy: mean estimation results (the 4 raw-cell coefficients) with a // single-coefficient set, breaking the subsequent `_b[raw10]' etc. // references in H's own nlcom call below. tempname _ov_logit_h _ov_logit_se capture qui nlcom (overlap_logit: logit(_b[`raw11'])) if _rc { di as error "nlcom (logit CI for overlap) failed (rc=" _rc ") " /// "-- overlap CI is undefined (this can happen at, " /// "or very near, a boundary cell proportion)" scalar `overlap_lb' = . scalar `overlap_ub' = . } else { scalar `_ov_logit_h' = r(b)[1,1] scalar `_ov_logit_se' = sqrt(r(V)[1,1]) scalar `overlap_lb' = invlogit(`_ov_logit_h' - `crit'*`_ov_logit_se') scalar `overlap_ub' = invlogit(`_ov_logit_h' + `crit'*`_ov_logit_se') } tempname h se lb ub N scalar `N' = `_Ndeg' capture qui nlcom (H: 1 - _b[`_cell10ref']/((_b[`raw11']+_b[`_cell10ref'])* /// (1-(_b[`raw11']+_b[`_cell01ref'])))), post if _rc { di as error "nlcom failed (rc=" _rc ") -- H/SE/CI are undefined " /// "(this can happen at, or very near, a boundary " /// "cell proportion, which makes the Jacobian " /// "singular)" return scalar loevh = . return scalar se = . return scalar lb = . return scalar ub = . return scalar N = `N' return scalar overlap = `overlap' return scalar se_overlap = `overlap_se' return scalar lb_overlap = `overlap_lb' return scalar ub_overlap = `overlap_ub' return local error "nlcom_failed" exit } if _rc==0 { scalar `h' = _b[H] scalar `se' = _se[H] scalar `lb' = `h' - `crit'*`se' scalar `ub' = `h' + `crit'*`se' } // Display: survey design summary header, then -- ONLY if `table' is // specified -- the cross-tab of design-based cell percentages, // obtained by simply wrapping Stata's own svy: tabulate (subpop()'d // to the current touse restriction). This reproduces svy: tabulate's // exact row/column variable name and value-label display, including // its automatic word-wrapping of arbitrarily long variable labels, // rather than reimplementing that logic by hand. The 4 cell // proportions/SEs already extracted above (from the svy: mean call // used for H itself) are NOT used for this display -- svy: tabulate // computes its own (numerically identical) design-based cell // percentages directly -- so the table body and the H/SE/CI results // below it remain guaranteed consistent, just via two separate (but // equivalent) svy estimations instead of one. // // Mirrors loevh2's own `table' option (default: no cross-tab shown // at all; with `table', the full cross-tab is displayed). When // `table' is NOT specified, the design summary header is instead // built by hand from the e()-scalars already saved above (off the // svy: mean call), since svy: tabulate itself is not run in this // case; this keeps the design summary (Number of strata/PSUs/obs/ // population size/subpop size/design df) always visible regardless // of `table'. Full variable names are always used (no abbreviation) // in the H results table below, matching svy: tabulate's own // unabbreviated variable-name display in the `table' case. local v1a "`v1'" local v2a "`v2'" local cil = string(`level') di _n as txt "Survey: Cross tabulation for Loevinger's H" if "`table'"=="table" { svy, subpop(`touse'): tabulate `v1' `v2', cell percent pearson // Overlap line: placed here, between the cross-tab (above) and // the H results table (below), in a narrow two-line format // (variable names on their own first line, values -- shown as // PERCENTAGES, i.e. x100, matching the cross-tab's own cell- // percentage scale -- on a second, narrower line). di _n as txt "Overlap (" as res "`v1a'" as txt "=1 & " as res "`v2a'" as txt "=1):" di as txt " " as res %7.2f `=scalar(`overlap')*100' /// as txt " Std. err. " as res %5.2f `=scalar(`overlap_se')*100' /// as txt " [" as res "`cil'" as txt "% CI (logit): " /// as res %5.2f `=scalar(`overlap_lb')*100' " " %5.2f `=scalar(`overlap_ub')*100' as txt "]" } else { di _n as txt "Number of strata = " as res %8.0fc scalar(`_e_Nstrata') /// as txt _col(36) "Number of obs" as txt " = " as res %12.0gc scalar(`_e_N') di as txt "Number of PSUs = " as res %8.0fc scalar(`_e_Npsu') /// as txt _col(36) "Population size" as txt " = " as res %12.0gc scalar(`_e_Npop') di as txt _col(36) "Subpop. no. obs" as txt " = " as res %12.0gc scalar(`_e_Nsub') di as txt _col(36) "Subpop. size" as txt " = " as res %12.0gc scalar(`_e_Nsubpop') di as txt _col(36) "Design df" as txt " = " as res %12.0gc scalar(`_e_dfr') } local vars "`v1a' `v2a'" local colsv = max(length("`vars'"), 17) + 2 local inc = `colsv' - 19 local dig = ceil(log10(max(1,scalar(`N')))) + 2 local nd = 8 - `dig' local spaces = cond(length("`cil'") == 2, " ", " ") tempname z p scalar `z' = `h'/`se' scalar `p' = 2*(1-normal(abs(`z'))) di _n in smcl in gr _col(`=58+`inc'+`nd'') /* */ "Number of obs = " as res %`dig'.0fc `N' /* */ _n _n in smcl in gr _col(`=23 + `inc'') "Loevinger (svy)" /* */ _n in smcl in gr /* */ " Variables" _col(`colsv') " {c |}" _col(`=18+`inc'') /* */ _col(`=25+`inc'') "H Coeff" /* */ _col(`=34+`inc'') "Std. err." /* */ _col(`=51+`inc'') "z" /* */ _col(`=58+`inc'') "p" /* */ _col(`=59+`inc'') `"`spaces'[`=strsubdp("`level'")'% conf. interval]"'/* */ _n "{hline `colsv'}{c +}{hline 61}" di in smcl in gr /* */ " `vars'" _col(`colsv') " {c |}" _col(`=18+`inc'') as res /* */ _col(`=25+`inc'') %7.0g `h' /* */ _col(`=36+`inc'') %7.0g `se' /* */ _col(`=45+`inc'') %7.2f `z' /* */ _col(`=51+`inc'') %7.3f `p' /* */ _col(`=62+`inc'') %7.0g `lb' /* */ _col(`=75+`inc'') %7.0g `ub' // by-group label, for r(lastgroup) and the compare frame // // Each by-variable can be either numeric (in which case its value // label, if any, is looked up via "summarize" + ": label") or a // string (in which case -summarize- does not apply -- it silently // produces no r(min), leaving the label empty, which caused blank // sub-sample labels whenever the by-variable was a string). For a // string by-variable, its own text value (for this by-group) is // used directly as the label. local sgroupslab "" if _by() { foreach _bv of local _byvars { capture confirm string variable `_bv' if !_rc { qui levelsof `_bv' if `touse', local(_bv_lvls) clean local _bv_lab : word 1 of `_bv_lvls' } else { qui summarize `_bv' if `touse', meanonly local _bv_val = r(min) local _bv_lab : label (`_bv') `_bv_val' } local sgroupslab = trim("`sgroupslab' `_bv_lab'") } } // The "svyressubsamp" frame (posting each valid by-group's h, se, // n, and overlap statistics, together with its label) is now built // for EVERY valid by-group, UNCONDITIONALLY (regardless of // -compare-) -- mirroring loevh2.ado's own "ressubsamp" frame // handling -- so that the last valid by-group's results can always // be retrieved and forwarded to r() (r(lastgroup), r(loevh), // r(se), r(N), r(overlap), etc.) even when -compare- is not // requested and the by-sequence happens to end on a // missing/excluded group (see the missing-byvar branch above). Only // the loevh2_svy_compare() call itself, and its compare-specific // returns (Hbar, chi2, df, p_chi2, H_SE_N), remain conditional on // -compare-. if _by() { if _byindex()==1 { capture frame drop svyressubsamp } capture confirm frame svyressubsamp if _rc { frame create svyressubsamp sample h se n ov ov_se ov_lb ov_ub } local sample = _byindex() frame post svyressubsamp (`sample') (`h') (`se') (`N') /// (`overlap') (`overlap_se') (`overlap_lb') (`overlap_ub') frame svyressubsamp: lab def sample `sample' "`sgroupslab'", modify if _bylastcall() { if "`compare'"=="compare" { di _n _n as txt "Categories of by: variable(s):" fre `_byvars' `if' `in' } frame change svyressubsamp lab val sample sample if "`compare'"=="compare" { loevh2_svy_compare local pchi2 = `loevh2_pchi2' local df = `loevh2_df' local chi2 = `loevh2_chi2' local Hbar = `loevh2_Hbar' local Hbar_se = `loevh2_Hbar_se' local Hbar_N = `loevh2_Hbar_N' local HSEN = "`loevh2_HSEN'" if `"`meta'"' != `""' & "`HSEN'"!="" { local HSENLAB = `"`loevh2_HSEN_labels'"' loevh2_meta_save, metaspec(`meta') /// source(loevh2_svy) suffix("_svy") /// setype("svy, nlcom (delta method)") /// hsen(`HSEN') labels(`"`HSENLAB'"') /// svywtype(`"`_svy_wtype'"') svywexp(`"`_svy_wexp'"') /// svycluster(`"`_svy_clust'"') svystrata(`"`_svy_strata'"') /// svyfpc(`"`_svy_fpc'"') /// origframe(`"`_origframe'"') origchanged(`_origchanged') } } frame change `_origframe' frame drop svyressubsamp } } // meta() single-row save: fires whenever meta() is specified and we // are NOT in the by:+compare pooling branch just above (which, when // it fires on the last by-group, already saves ALL rows -- including // this group's -- via the HSEN matrix). See loevh2.ado's own // equivalent block for the full rationale (identical logic here): // (a) a plain, non-by: call -> label = "var1_var2" // (b) by: without compare -> label = sgroupslab, // one row saved per by-group as each is processed // (c) by: + compare, for every by-group EXCEPT the last one. if `"`meta'"' != `""' & "`compare'"!="compare" { local _metalab = cond(_by(), "`sgroupslab'", "`v1'_`v2'") loevh2_meta_save, metaspec(`meta') source(loevh2_svy) /// suffix("_svy") setype("svy, nlcom (delta method)") /// label(`"`_metalab'"') hval(`=`h'') seval(`=`se'') nval(`=`N'') /// svywtype(`"`_svy_wtype'"') svywexp(`"`_svy_wexp'"') /// svycluster(`"`_svy_clust'"') svystrata(`"`_svy_strata'"') /// svyfpc(`"`_svy_fpc'"') /// origframe(`"`_origframe'"') origchanged(`_origchanged') } // Returns return scalar loevh = `h' return scalar se = `se' return scalar lb = `lb' return scalar ub = `ub' return scalar level = `level' return scalar N = `N' return scalar overlap = `overlap' return scalar se_overlap = `overlap_se' return scalar lb_overlap = `overlap_lb' return scalar ub_overlap = `overlap_ub' if _bylastcall() & "`compare'"=="compare" { return scalar Hbar = `Hbar' return scalar Hbar_se = `Hbar_se' return scalar chi2 = `chi2' return scalar df = `df' return scalar p_chi2 = `pchi2' return scalar Hbar_N = `Hbar_N' if "`HSEN'" != "" return matrix H_SE_N = `HSEN' } return local var1 "`v1'" return local var2 "`v2'" return local se_type "svy, nlcom (delta method)" if _by() return local group "`_byvars'" if _by() return local lastgroup "`sgroupslab'" end * ------------------------------------------------------------------------------ program define loevh2_svy_compare version 16.0 tempname wsum wrsum wr2sum Hbar Hbar_se Hbar_N chi2 df pchi2 scalar `Hbar' = . scalar `Hbar_se' = . scalar `Hbar_N' = . scalar `chi2' = . scalar `df' = . scalar `pchi2' = . tempvar ok qui gen byte `ok' = !missing(h) & !missing(se) qui count if `ok' local k = r(N) if `k' > 1 { mata: loevh2_svy_hs = st_data(., ("h", "se"), "`ok'") mata: st_numscalar("`wsum'", sum(1 :/ loevh2_svy_hs[.,2]:^2)) mata: st_numscalar("`wrsum'", sum(loevh2_svy_hs[.,1] :/ loevh2_svy_hs[.,2]:^2)) mata: st_numscalar("`wr2sum'", sum(loevh2_svy_hs[.,1]:^2 :/ loevh2_svy_hs[.,2]:^2)) mata: mata drop loevh2_svy_hs scalar `Hbar' = `wrsum' / `wsum' // Pooled (inverse-variance-weighted meta-analysis) SE of the // weighted average H -- SE(Hbar) = sqrt(1 / sum(1/se_i^2)) -- // plus the total N summed across all valid sub-samples, both // shown alongside Hbar on the "Weighted average H" summary line // below. scalar `Hbar_se' = sqrt(1/`wsum') qui summarize n if `ok', meanonly scalar `Hbar_N' = r(sum) scalar `chi2' = `wr2sum' - (`wrsum')^2 / `wsum' scalar `df' = `k' - 1 scalar `pchi2' = chi2tail(`df', `chi2') local hlwidth = 30 local maxlablen = length("Weighted average H") local N = _N forvalues obs = 1/`N' { local sval = sample[`obs'] local slab : label (sample) `sval' if length("`slab'") > `maxlablen' local maxlablen = length("`slab'") } local labcol = `maxlablen' + 3 di _n as txt "Test of equality of H's across " as res `N' /// as txt " sub-samples:" di _n in gr /// _col(2) "Sub-sample" /// _col(`labcol') "{c |}" /// _col(`=`labcol'+3') "H Coeff" /// _col(`=`labcol'+13') "Std. err." /// _col(`=`labcol'+32') "N" di in gr "{hline `=`labcol'-1'}{c +}{hline `=`hlwidth'+2'}" capture matrix drop __loevh2_svy_HSEN matrix __loevh2_svy_HSEN = J(`N', 3, .) local HSEN "__loevh2_svy_HSEN" matrix colnames `HSEN' = H SE N local HSEN_rownames "" local HSEN_labels_full "" forvalues obs = 1/`N' { local hval = h[`obs'] local seval = se[`obs'] local nval = n[`obs'] local sval = sample[`obs'] local slab : label (sample) `sval' if missing(`hval') | missing(`seval') continue di in gr _col(2) "`slab'" /// _col(`labcol') "{c |}" as res /// _col(`=`labcol'+3') %7.4f `hval' /// _col(`=`labcol'+15') %7.4f `seval' /// _col(`=`labcol'+20') %11.0fc `nval' matrix `HSEN'[`obs',1] = `hval' matrix `HSEN'[`obs',2] = `seval' matrix `HSEN'[`obs',3] = `nval' // Row names are simple, guaranteed-short, collision-free // "group N" identifiers (N = this row's position, 1-based) -- // NOT sanitized copies of the descriptive sub-sample label // (see loevh2.ado's loevh2_compare for the full rationale). // The full, untruncated label is carried forward separately // via `HSEN_labels_full' for loevh2_meta_save's numeric, // value-labelled `study' variable. local HSEN_rownames `"`HSEN_rownames' "group `obs'""' local HSEN_labels_full `"`HSEN_labels_full' `"`slab'"'"' } matrix rownames `HSEN' = `HSEN_rownames' di in gr "{hline `=`labcol'-1'}{c +}{hline `=`hlwidth'+2'}" di in gr _col(2) "Weighted average H" /// _col(`labcol') "{c |}" as res /// _col(`=`labcol'+3') %7.4f scalar(`Hbar') /// _col(`=`labcol'+15') %7.4f scalar(`Hbar_se') /// _col(`=`labcol'+20') %11.0fc scalar(`Hbar_N') di in gr "{hline `=`labcol'-1'}{c BT}{hline `=`hlwidth'+2'}" di _n as txt " Pearson chi2(" as res %1.0f scalar(`df') as text ") = " /// as res %6.4f scalar(`chi2') /// as text " Pr = " as res %5.3f scalar(`pchi2') c_local loevh2_HSEN "`HSEN'" c_local loevh2_HSEN_labels `"`HSEN_labels_full'"' } else { di as error "Note: fewer than 2 valid sub-samples; " /// "cannot test equality of H's" c_local loevh2_HSEN "" c_local loevh2_HSEN_labels "" } c_local loevh2_chi2 = `chi2' c_local loevh2_df = `df' c_local loevh2_pchi2 = `pchi2' c_local loevh2_Hbar = `Hbar' c_local loevh2_Hbar_se = `Hbar_se' c_local loevh2_Hbar_N = `Hbar_N' end