*! svylet.ado v1.6 - 27aug2026
*! Wald omnibus F-test + Bonferroni pairwise comparisons + Compact Letter
*! Display (CLD) for means, totals, proportions, and ratios under survey
*! design.
*! Adds r(letra_1..k) and r(nombre_categoria_1..k) so the CLD letters can
*! be read programmatically after the command runs, not only displayed.
*! Author: Andres Talavera Cuya. Affiliation stated for identification
*! purposes only -- this software is not an official product of INEI
*! and INEI bears no responsibility for it. Developed independently by
*! the author; views, methods, and results are the author's own and do
*! not necessarily reflect the position of INEI.
*! Distributed under the GNU General Public License v3
*! (https://www.gnu.org/licenses/gpl-3.0.txt).
*!
*! v1.5 -- FIX real (encontrado corriendo revision_pre_ssc.do en Stata
*! real, StataNow 19.5, antes de mandar el paquete a SSC): over() STRING
*! (ej. una variable armada con decode()) hacia que "over() debe tener
*! al menos 2 categorias en la muestra" saltara SIEMPRE, sin importar
*! cuantas categorias reales tuviera la variable. Causa raiz: la linea
*! "markout `touse' `over' `denominator'" (sin la opcion strok) hace que
*! Stata excluya el 100% de las observaciones en cuanto CUALQUIER
*! variable de su lista es de tipo string -- no solo las de string vacio,
*! TODAS -- confirmado aislando el markout solo (fuera de svylet.ado) con
*! datos donde la variable string no tenia NINGUN valor vacio. touse
*! quedaba en 0 para las 74 observaciones de auto.dta, "levelsof origin
*! if touse" devolvia vacio, y el chequeo de k<2 saltaba por una causa
*! que no tenia nada que ver con el numero real de categorias. Afecta a
*! CUALQUIER llamada con over() string -- el camino numerico (el uso mas
*! comun) nunca estuvo afectado, porque strok no cambia el manejo de
*! variables numericas. Corregido agregando ", strok" al markout. Este
*! bug estaba presente desde v1.0; el ejemplo de over() string agregado
*! al help en esta misma sesion (svylet.sthlp) nunca se habia corrido en
*! Stata real hasta ahora -- ver revision_pre_ssc.do en el repo, seccion 5.
*!
*! v1.6 -- FIX (encontrado revisando tabulados de produccion reales,
*! tipo_actividad/Pecuaria): filas con las 4 categorias de over() con
*! estimacion presente (ninguna realmente faltante) igual mostraban el F
*! omnibus en blanco, porque 1-3 de esas 4 categorias tenian proporcion
*! exactamente 0% (varianza degenerada). Causa raiz: el contraste R del
*! F omnibus (svylet_core(), Mata) se armaba con las k categorias
*! COMPLETAS -- una sola con V[i,i] missing o <=0 propaga "." por
*! R*V*R' y lusolve() devuelve missing para TODO el estadistico, aunque
*! el resto de las k-1 categorias tuvieran variables perfectamente
*! calculables. Mismo criterio que ya se aplicaba a los pares (Pmat: una
*! categoria degenerada solo invalida SUS pares, no todos) -- ahora el F
*! omnibus se calcula sobre el SUBCONJUNTO de categorias con varianza
*! definida (minimo 2, no se puede testear una igualdad con menos). Con
*! 0 categorias degeneradas el resultado es identico al de v1.5 (mismo
*! contraste, misma formula). Replicado en tsvy.ado (tsvy_core(), que
*! mantiene su propia copia de este motor -- ver nota ahi sobre por que
*! no puede llamar a svylet_core() entre archivos .ado separados de SSC).
*!
*! v1.5 (parte 2, mismo dia, corriendo el FIX de arriba en Stata real):
*! con touse corregido, over() string sigue sin funcionar -- ahora falla
*! en el "svy: mean/total/proportion/ratio ..., over(`over')" interno con
*! el error NATIVO de Stata "string variables not allowed in varlist;
*! invalid over() option" (r(109)). Esto es una limitacion DURA de
*! svy: mismo (over() debe ser numerico, sin excepcion) -- no algo que
*! svylet pueda evitar sin re-mapear el string a codigos numericos por
*! dentro (no intentado; agrega complejidad no validada en Stata real
*! justo antes de una entrega). La documentacion previa ("varname may be
*! numeric or string") era incorrecta desde antes de esta sesion -- ver
*! svylet.sthlp/svylet_es.sthlp. Agregado un chequeo temprano (justo
*! despues de syntax) que detecta over() string y corta con un mensaje
*! claro sugiriendo encode(), en vez de dejar que el usuario llegue al
*! error crudo de Stata mas abajo.
*!
*! v1.1 -- auditoria: dos cambios de comportamiento respecto a v1.0:
*!   1. boot(): la reasignacion de pseudo-grupo bajo H0 (prepivoting,
*!      Beran 1988 / Hall & Wilson 1991) ahora se hace a nivel de
*!      CONGLOMERADO (`newpsu' tras bsample), no a nivel de observacion.
*!      En v1.0, un sorteo por FILA + reparto en bloques de filas
*!      contiguas partia un mismo conglomerado entre varios pseudo-
*!      grupos, destruyendo la correlacion intra-conglomerado que el
*!      propio diseno (svyset) asume al calcular varianza -- eso sesga
*!      el p-valor bootstrap hacia el rechazo (ver simulacion en
*!      sim/simulacion_bootstrap.py: con v1.0 el error tipo I bajo H0
*!      queda inflado por encima del alpha nominal; con v1.1 queda
*!      centrado en el nominal). Respaldo: Rao & Wu (1988, JASA),
*!      Rust & Rao (1996, Stat Methods Med Res), Wolter (2007,
*!      "Introduction to Variance Estimation" 2nd ed, cap. 6),
*!      Field & Welsh (2007, JRSS-B), Hall & Wilson (1991, Biometrics),
*!      Canty & Davison (1999, The Statistician), Davison & Hinkley
*!      (1997, "Bootstrap Methods and Their Application", cap. 3 y 9).
*!   2. refgroup() ELIMINADO de la sintaxis. Estaba declarado desde v1.0
*!      pero nunca se usaba en el programa ni en el motor Mata -- el
*!      contraste F omnibus siempre corria contra el grupo 1. Al revisar
*!      la matematica: el estadistico de Wald para "todos los grupos
*!      iguales" es INVARIANTE al grupo de referencia elegido, porque
*!      cualquier conjunto de k-1 contrastes independientes que genere
*!      el mismo subespacio da el mismo Wald (R2 = A*R para A invertible
*!      => (R2 b)'(R2 V R2')^-1(R2 b) = (R b)'(R V R')^-1(R b)). Por eso
*!      no se "implementa" refgroup(): no hay nada que implementar que
*!      cambie el resultado. Dejarlo en la sintaxis como no-op silencioso
*!      es peor que quitarlo -- un usuario que lo pase ahora obtiene un
*!      error claro de Stata ("option refgroup() not allowed") en vez de
*!      creer que cambio algo cuando no lo hizo.
*!   3. Nota agregada: level() solo tiene efecto para stat(proportion).
*!      Si se pasa level() distinto del default junto con mean/total, se
*!      muestra un aviso (no es un error -- level() simplemente no aplica
*!      ahi) en vez de ignorarlo en silencio.
*!
*! v1.2 -- agrega stat(ratio): svy: ratio numerador/denominador, over().
*!   El denominador se pasa por la opcion denominator(varname), no por
*!   texto "num/den" en varlist -- svylet arma la sintaxis "num/den" que
*!   svy: ratio espera internamente. e(b)/e(V) de svy: ratio tienen la
*!   MISMA forma que mean/total (una sola ecuacion, k columnas, una por
*!   categoria de over()) -- por eso _svylet_seleccionar no necesito
*!   cambios: su rama "else" (ya usada por mean/total) sirve igual para
*!   ratio, sin logica nueva de seleccion de ecuacion.
*!
*! v1.3 -- FIX: r(n_ponderado)/r(n_sin_ponderar) devolvian los valores
*!   INVERTIDOS desde v1.0 (e(_N) mapeado a "ponderado" y e(_N_subp) a
*!   "sin ponderar" -- al reves). Encontrado en produccion via tsvy:
*!   en la columna N_SIN_PON del frame acumulador aparecia un numero del
*!   orden de millones (identico a ESTIMA en stat(total) con varname
*!   constante =1, que por construccion ES el tamano de poblacion
*!   ponderado) en vez de un conteo de casos de muestra. Ver nota junto a
*!   "tempname Nmat Nsubpmat Rtable" en el cuerpo del programa. Afecta a
*!   TODO llamador de svylet que use r(n_ponderado)/r(n_sin_ponderar) --
*!   en particular tsvy (columnas N_PONDERA/N_SIN_PON del frame
*!   acumulador). No afecta r(b)/r(V) ni el test F/CLD -- esas cantidades
*!   nunca pasaron por Nmat/Nsubpmat.
*!
*! v1.4 -- agrega ref(): comparaciones "cada categoria de over() contra
*!   una categoria base fija", Bonferroni sobre k-1 comparaciones -- una
*!   FAMILIA DE HIPOTESIS DISTINTA de la que responde el F omnibus + CLD
*!   existente (que compara TODOS los pares, Bonferroni sobre
*!   k(k-1)/2). Ver svylet.sthlp, seccion Remarks, para la justificacion
*!   y las referencias (Dunn 1961; Dunnett 1955, 1964; Hsu 1996) sobre
*!   por que estas dos preguntas dan, legitimamente, resultados
*!   distintos para el MISMO dato -- encontrado al comparar las letras
*!   CLD de tsvy contra un script de referencia que solo comparaba cada
*!   anio contra el anio mas reciente (Bonferroni sobre 3 comparaciones,
*!   no sobre 6): ninguna de las dos salidas estaba mal, respondian
*!   preguntas distintas.
*!
*! See svylet.sthlp for full documentation, options, examples, and
*! references.

program define svylet, rclass
    version 14
    syntax varname(numeric) [if] [in], OVER(varname) ///
        STAT(string) [ALPHA(real 0.05) LEVEL(integer 1) ///
        DENOMINATOR(varname numeric) BOOT(integer 0) BSEED(integer -1) ///
        REF(string)]

    * over() DEBE ser numerico -- svy: mean/total/proportion/ratio no
    * aceptan un over() string en absoluto ("string variables not
    * allowed in varlist; invalid over() option", r(109)), confirmado en
    * Stata real (StataNow 19.5). Rechazar ACA con un mensaje claro, antes
    * de marksample/markout, en vez de dejar que el usuario llegue al
    * error crudo de Stata mas abajo -- sugiere encode() como salida.
    capture confirm string variable `over'
    if !_rc {
        di as err "over(`over') es una variable string -- svy: `stat' no acepta over() string."
        di as err "Solucion: encode `over', gen(nombre_numerico) y use over(nombre_numerico)."
        exit 109
    }

    local stat = lower("`stat'")
    if !inlist("`stat'", "mean", "proportion", "total", "ratio") {
        di as err "stat() debe ser mean, proportion, total o ratio"
        exit 198
    }
    if "`stat'" != "proportion" & `level' != 1 {
        di as text "Nota: level() solo tiene efecto con stat(proportion) -- se ignora para `stat'."
    }
    if "`stat'" == "ratio" & "`denominator'" == "" {
        di as err "stat(ratio) requiere denominator(varname) -- el numerador es `varlist'"
        exit 198
    }
    if "`stat'" != "ratio" & "`denominator'" != "" {
        di as text "Nota: denominator() solo tiene efecto con stat(ratio) -- se ignora para `stat'."
    }

    marksample touse
    * strok: SIN esta opcion, markout excluye TODAS las observaciones
    * (no solo las de string vacio) en cuanto `over' es string -- no un
    * comportamiento documentado de forma obvia, confirmado en Stata real
    * (ver changelog v1.4.1 arriba): "markout touse over" sobre una
    * variable string (ej. over() con decode()) dejaba touse=0 para el
    * 100% de las observaciones, y el error resultante ("over() debe
    * tener al menos 2 categorias") no tenia nada que ver con la causa
    * real. strok le dice a markout que una variable string no-missing
    * (incluida "") se trate como valida -- las observaciones con string
    * REALMENTE vacio en `over' siguen quedando excluidas mas abajo,
    * porque levelsof nunca las cuenta como una categoria (una categoria
    * "" no es una categoria valida de over() para el test).
    markout `touse' `over' `denominator', strok

    quietly levelsof `over' if `touse', local(niveles_over)
    local k : word count `niveles_over'
    if `k' < 2 {
        di as err "over() debe tener al menos 2 categorias en la muestra"
        exit 198
    }

    * -- v1.4 -- ref(): comparaciones "cada categoria vs una categoria
    * base", Bonferroni sobre k-1 comparaciones -- PREGUNTA DISTINTA de
    * la que responde el F omnibus + CLD de mas abajo (que compara TODOS
    * los pares, Bonferroni sobre k(k-1)/2). Ver help para la cita de
    * literatura sobre por que estas dos familias de hipotesis dan,
    * legitimamente, resultados distintos (Dunn 1961; Dunnett 1955).
    local ref_idx = 0
    if "`ref'" != "" {
        local i = 0
        foreach v of local niveles_over {
            local i = `i' + 1
            if `v' == `ref' local ref_idx = `i'
        }
        if `ref_idx' == 0 {
            di as err "svylet: ref(`ref') no es un valor observado de over(`over') en la muestra."
            di as err "  valores observados: `niveles_over'"
            exit 198
        }
    }

    * -- Correr el comando svy correspondiente ------------------------------
    if "`stat'" == "mean" {
        svy: mean `varlist' if `touse', over(`over')
    }
    else if "`stat'" == "total" {
        svy: total `varlist' if `touse', over(`over')
    }
    else if "`stat'" == "ratio" {
        svy: ratio `varlist'/`denominator' if `touse', over(`over')
    }
    else {
        svy: proportion `varlist' if `touse', over(`over')
    }

    tempname B V
    matrix `B' = e(b)
    matrix `V' = e(V)
    local df_r = e(df_r)
    if "`df_r'" == "" {
        di as err "e(df_r) no disponible tras el svy -- revisar version de Stata/diseno"
        exit 498
    }
    * Capturar N, N_subp, y r(table) INMEDIATAMENTE, antes de que
    * cualquier otro comando pueda tocar e()/r() -- mismo criterio que
    * ya establecimos: guardar apenas esta disponible, no confiar en que
    * sobreviva hasta mas adelante en el programa.
    *
    * v1.2 -- FIX: `Nmat' (e(_N)) y `Nsubpmat' (e(_N_subp)) estaban
    * INVERTIDOS al exponerse como r(n_ponderado)/r(n_sin_ponderar) desde
    * v1.0. Confirmado con datos reales (ENA, stat(total) con varname
    * constante =1): para NIVEL=NACIONAL/ANIO=2026, e(_N) devolvia
    * 20783 (el N de casos SIN ponderar, coincide con el conteo de filas
    * de la muestra ese anio) y e(_N_subp) devolvia 2057264.91... --
    * EXACTAMENTE igual a r(b) (la estimacion PONDERADA, ya que
    * varname=1 hace que total(varname) sea el tamano de poblacion
    * ponderado). O sea e(_N) es el tamano SIN ponderar y e(_N_subp) es
    * el PONDERADO -- al reves de lo que el nombre "_N_subp" sugiere y de
    * lo que este mismo archivo asumia. Los nombres de macro (`Nmat'/
    * `Nsubpmat') se mantienen como estaban (mapeados 1 a 1 a e(_N)/
    * e(_N_subp) respectivamente, para no confundir de donde viene cada
    * uno) -- el swap se hace mas abajo, en el return matrix, que es
    * donde se decide el significado (ponderado/sin ponderar) que ve
    * quien llama.
    tempname Nmat Nsubpmat Rtable
    capture matrix `Nmat' = e(_N)
    capture matrix `Nsubpmat' = e(_N_subp)
    capture matrix `Rtable' = r(table)

    * -- Seleccionar las filas relevantes de b/V (subrutina reutilizable, ----
    *    la necesita tambien el loop de boot() mas abajo) --------------------
    tempname b_sel V_sel
    _svylet_seleccionar `B' `V' `k' `level' "`stat'" "`varlist'" `b_sel' `V_sel'
    local idx_usar "`r(idx_usar)'"

    * -- v1.4 -- ref(): k-1 contrastes de Wald "categoria vs ref_idx",
    * Bonferroni sobre k-1 (no sobre k(k-1)/2 como el CLD). Usa el MISMO
    * b_sel/V_sel/df_r ya calculados -- ningun svy: adicional.
    tempname p_vsref p_vsref_raw
    matrix `p_vsref'     = J(`k', 1, .)
    matrix `p_vsref_raw' = J(`k', 1, .)
    if `ref_idx' > 0 {
        mata: svylet_vsref("`b_sel'", "`V_sel'", `df_r', `ref_idx', `k', "`p_vsref'", "`p_vsref_raw'")
    }

    * N ponderado / sin ponderar por categoria, en el mismo orden que
    * b_sel/V_sel (usa los mismos indices de columna que ya selecciono
    * _svylet_seleccionar). Si e(_N)/e(_N_subp) no existieran por algun
    * motivo (comando svy distinto sin esa info), quedan en missing en
    * vez de cortar la ejecucion -- esto es informativo, no critico para
    * el test en si.
    tempname Nout Nsubpout
    matrix `Nout' = J(`k', 1, .)
    matrix `Nsubpout' = J(`k', 1, .)

    * Limite inferior/superior del IC por categoria: se leen DIRECTO de
    * r(table) (filas "ll"/"ul"), NO se reconstruyen con una formula
    * simetrica propia. r(table) ya trae el IC correcto para cada
    * stat() -- simetrico para mean/total, en escala logit para
    * proportion (respeta [0,1]) -- construirlo de nuevo a mano
    * arriesgaba reimplementar mal esa transformacion. Usar lo que Stata
    * ya calculo es mas robusto que reimplementarlo.
    tempname LIout LSout
    matrix `LIout' = J(`k', 1, .)
    matrix `LSout' = J(`k', 1, .)
    local fila = 0
    foreach i of local idx_usar {
        local fila = `fila' + 1
        capture matrix `Nout'[`fila', 1] = `Nmat'[1, `i']
        capture matrix `Nsubpout'[`fila', 1] = `Nsubpmat'[1, `i']
        capture matrix `LIout'[`fila', 1] = `Rtable'["ll", `i']
        capture matrix `LSout'[`fila', 1] = `Rtable'["ul", `i']
    }

    * -- Motor Mata: F omnibus + Bonferroni + CLD ----------------------------
    mata: svylet_core("`b_sel'", "`V_sel'", `df_r', `alpha', `k')

    * r() es un solo espacio global -- CUALQUIER comando que corra despues
    * (incluido todo el bloque boot(), que llama svy:/count muchas veces)
    * lo sobreescribe. Guardar TODO lo que hace falta para la salida final
    * en locales/matrices propias ACA, antes de que arranque boot() -- no
    * confiar en que r(...) siga disponible mas adelante en el programa.
    local F_obs     = r(F_omnibus)
    local p_obs     = r(p_omnibus)
    local n_pares   = r(n_pares)
    tempname b_obs
    matrix `b_obs' = r(b)
    forvalues i = 1/`k' {
        local letra_`i' "`r(letra_`i')'"
    }

    * -- boot(): prepivoting del F omnibus (Beran 1988) ----------------------
    * Ver notas del encabezado -- pool de conglomerados de TODOS los grupos,
    * dentro de su estrato original, reasignados al azar en la misma
    * proporcion de tamanos que el dato real (fuerza H0 sin asumir
    * independencia -- Hall & Wilson 1991).
    local p_boot_omni = .
    local B_efectivo = .
    if `boot' > 0 {
        * e(psu)/e(strata) (nombre simple) vs e(psu1)/e(strata1) (marco
        * multietapico) vs e(su)/e(su1) ("sampling unit", el termino que
        * usa el propio svyset en su salida: "Sampling unit 1: cluster5")
        * -- se prueban las tres formas. La hipotesis de psu/psu1 quedo
        * DESCARTADA con datos reales: incluso con conglomerados genuinos
        * de mas de un caso (Number of PSUs=20, df=18), e(psu) y e(psu1)
        * seguian vacios -- el nombre real del macro es otro.
        local psuvar    "`e(psu)'"
        local stratavar "`e(strata)'"
        if "`psuvar'" == "" {
            local psuvar "`e(psu1)'"
        }
        if "`psuvar'" == "" {
            local psuvar "`e(su)'"
        }
        if "`psuvar'" == "" {
            local psuvar "`e(su1)'"
        }
        if "`stratavar'" == "" {
            local stratavar "`e(strata1)'"
        }
        if "`psuvar'" == "" | "`stratavar'" == "" {
            di as err "boot() requiere que svyset tenga psu() y strata() declarados"
            di as err "  e(psu)=" "`e(psu)'" "  e(psu1)=" "`e(psu1)'" ///
                "  e(su)=" "`e(su)'" "  e(su1)=" "`e(su1)'"
            di as err "  e(strata)=" "`e(strata)'" "  e(strata1)=" "`e(strata1)'"
            di as err "  ninguna de las formas probadas (psu/psu1/su/su1)" ///
                " funciono -- correr ereturn list despues del svy: para" ///
                " ver el nombre real del macro en esta version de Stata."
            exit 498
        }

        * Tamanos por grupo medidos en NUMERO DE CONGLOMERADOS distintos
        * (no en numero de filas -- v1.1: la unidad que se reasigna bajo
        * H0 es el conglomerado, asi que las proporciones deben calcularse
        * sobre esa misma unidad).
        tempname n_psu_por_grupo
        matrix `n_psu_por_grupo' = J(`k', 1, .)
        local n_psu_total_orig = 0
        local i = 0
        foreach niv of local niveles_over {
            local i = `i' + 1
            quietly levelsof `psuvar' if `over' == `niv' & `touse', local(__psulet_lv)
            local n_psu_niv : word count `__psulet_lv'
            matrix `n_psu_por_grupo'[`i', 1] = `n_psu_niv'
            local n_psu_total_orig = `n_psu_total_orig' + `n_psu_niv'
        }

        * Acumular los F de cada replica en una MATRIZ de Stata en memoria,
        * no con postfile/tempfile. postfile esta pensado para construir
        * datasets completos, no para juntar 300 numeros sueltos -- y cada
        * capa de archivo temporal (tempfile + postfile, mas los 300
        * preserve/restore del loop) es una oportunidad mas de colision de
        * nombres en el directorio temp de Windows (la causa mas probable
        * del mensaje cosmetico "file ... not found" que aparecia antes,
        * incluso despues de un primer intento de arreglo via frame -- ese
        * intento seguia usando un archivo en disco para leer los
        * resultados, por eso no lo resolvio). Con una matriz en memoria no
        * hay ningun archivo de por medio para esta parte.
        tempname resultados_boot
        matrix `resultados_boot' = J(`boot', 1, .)

        if `bseed' >= 0 set seed `bseed'
        local reps_ok = 0
        forvalues rep = 1/`boot' {
            preserve
            quietly keep if `touse'
            tempvar newpsu psu_first psu_orden rank_psu grupo_asig grupo_pseudo
            capture noisily {
                bsample, cluster(`psuvar') strata(`stratavar') idcluster(`newpsu')

                * v1.1 -- Reasignar el pseudo-grupo a nivel de CONGLOMERADO
                * (`newpsu', el conglomerado ya remuestreado por bsample),
                * NO a nivel de observacion. Un solo sorteo por
                * conglomerado; el conglomerado entero se mueve como
                * bloque indivisible al pseudo-grupo que le toque -- nunca
                * se parte entre dos pseudo-grupos. Esto es lo que
                * preserva la correlacion intra-conglomerado que el propio
                * diseno (y el F observado) asume al calcular varianza
                * (ver referencias en el encabezado).
                quietly bysort `newpsu': gen byte `psu_first' = (_n==1)
                quietly count if `psu_first'
                local n_psu_actual = r(N)

                quietly gen double `psu_orden' = runiform() if `psu_first'
                quietly egen `rank_psu' = rank(`psu_orden') if `psu_first', unique

                quietly gen `grupo_asig' = .
                local acumulado = 0
                local asignado = 0
                local i = 0
                foreach niv of local niveles_over {
                    local i = `i' + 1
                    if `i' == `k' {
                        * el ultimo grupo se lleva el resto exacto, para
                        * garantizar que la suma de todos los ni sea
                        * exactamente n_psu_actual (evita perder o sobrar
                        * conglomerados por redondeo)
                        local ni = `n_psu_actual' - `asignado'
                    }
                    else {
                        local ni_orig = `n_psu_por_grupo'[`i', 1]
                        local ni = round(`n_psu_actual' * `ni_orig' / `n_psu_total_orig')
                    }
                    if `ni' > 0 {
                        quietly replace `grupo_asig' = `niv' if `psu_first' & ///
                            `rank_psu' > `acumulado' & `rank_psu' <= `acumulado' + `ni'
                    }
                    local acumulado = `acumulado' + `ni'
                    local asignado = `asignado' + `ni'
                }
                * Expandir la etiqueta asignada en el "primero de cada
                * conglomerado" a TODAS sus filas -- max() por `newpsu',
                * ya que dentro de cada conglomerado solo esa fila tiene
                * un valor no-missing.
                quietly bysort `newpsu': egen `grupo_pseudo' = max(`grupo_asig')

                quietly svyset `newpsu', strata(`stratavar') singleunit(certainty)
                if "`stat'" == "mean" {
                    quietly svy: mean `varlist', over(`grupo_pseudo')
                }
                else if "`stat'" == "total" {
                    quietly svy: total `varlist', over(`grupo_pseudo')
                }
                else if "`stat'" == "ratio" {
                    quietly svy: ratio `varlist'/`denominator', over(`grupo_pseudo')
                }
                else {
                    quietly svy: proportion `varlist', over(`grupo_pseudo')
                }
                tempname Bb Vb b_sel_b V_sel_b
                matrix `Bb' = e(b)
                matrix `Vb' = e(V)
                _svylet_seleccionar `Bb' `Vb' `k' `level' "`stat'" "`varlist'" `b_sel_b' `V_sel_b'
                mata: st_numscalar("r(Frep)", svylet_fomnibus_only(st_matrix("`b_sel_b'"), st_matrix("`V_sel_b'"), `df_r', `k'))
            }
            if _rc == 0 {
                matrix `resultados_boot'[`rep', 1] = r(Frep)
                local reps_ok = `reps_ok' + 1
            }
            restore
        }

        * Agregacion final directo en Mata sobre la matriz -- sin dataset,
        * sin archivo, sin frame. select() descarta las filas "." (replicas
        * que fallaron); el resto es la misma formula de siempre:
        *   p = (# replicas con F_replica >= F_obs, +1) / (B_efectivo + 1)
        local B_efectivo = 0
        local p_boot_omni = .
        mata: svylet_agregar_boot("`resultados_boot'", `F_obs')
        local B_efectivo = r(B_efectivo)
        if `B_efectivo' > 0 {
            local p_boot_omni = r(p_boot)
        }
    }

    * -- Salida a la ventana de resultados (todavia sin exportar a archivo) --
    di as text _n "{hline 70}"
    local etiqueta_var "`varlist'"
    if "`stat'" == "ratio" local etiqueta_var "`varlist'/`denominator'"
    di as text "svylet -- `stat' de `etiqueta_var', over(`over')"
    di as text "{hline 70}"
    di as text "F de Wald (ANOVA global, ajustado por diseno):"
    if `F_obs' == . {
        di as err "  F no calculable: al menos una categoria de over() tiene" ///
            " varianza no definida (proporcion exactamente 0 o 1 sin" ///
            " variabilidad, u otra causa de SE indefinido) -- ver aviso arriba."
    }
    else {
        di as text "  F(" %3.0f `=`k'-1' ", " %6.0f `=`df_r'-(`k'-1)+1' ") = " as res %9.4f `F_obs' ///
            as text "   Prob > F (analitico) = " as res %6.4f `p_obs'
    }
    if `boot' > 0 {
        if `B_efectivo' > 0 {
            di as text "  Prob > F (bootstrap, prepivotado, B=" %4.0f `B_efectivo' ///
                " de " %4.0f `boot' " replicas validas) = " as res %6.4f `p_boot_omni'
        }
        else {
            di as err "  bootstrap: ninguna replica valida -- revisar diseno/psu/strata"
        }
    }
    di as text _n "Grupos (letras compartidas = NO significativamente distintos," ///
        " Bonferroni analitico sobre " %2.0f `n_pares' " pares -- NO recalibrado por boot();" ///
        " [?] = varianza no definida, NO tratar como sin diferencia):"
    forvalues i = 1/`k' {
        local nom_grupo : word `i' of `niveles_over'
        di as text "  `over' = `nom_grupo'" _col(20) as res %9.4f el(`b_obs', `i', 1) ///
            as text "   grupo: " as res "`letra_`i''"
    }

    return scalar F_omnibus      = `F_obs'
    return scalar p_omnibus      = `p_obs'
    return scalar p_omnibus_boot = `p_boot_omni'
    return scalar B_efectivo     = `B_efectivo'
    return scalar df_num         = `k' - 1
    return scalar df_den          = `df_r' - (`k'-1) + 1
    return scalar k_categorias   = `k'
    * Grados de libertad de diseno SIN el ajuste Korn-Graubard (el que
    * necesita un IC estandar tipo b +/- t(df)*se) -- df_den de arriba
    * es el ajustado para el F omnibus, un numero distinto a proposito
    * (ver help, seccion Korn-Graubard). Devolver los dos por separado
    * para que quien arme una tabla de IC no use por error el ajustado.
    return scalar df_raw          = `df_r'
    * Estimaciones puntuales, varianza, y tamanos de muestra por
    * categoria -- lo que hace falta para armar una tabla de punto
    * (estimacion + error estandar + IC + casos) sin tener que correr
    * svy: total/mean/proportion una SEGUNDA vez por separado. svylet ya
    * corrio esa estimacion internamente -- exponerla evita duplicar el
    * trabajo de estimacion (confirmado con datos reales de ENA: correr
    * la estimacion dos veces en paralelo fue la fuente de varios bugs
    * de alineacion en el pipeline de especies).
    return matrix b               = `b_obs'
    return matrix V               = `V_sel'
    * Swap deliberado -- ver nota junto a "tempname Nmat Nsubpmat Rtable"
    * mas arriba: e(_N) (-> `Nout') resulto ser el tamano SIN ponderar y
    * e(_N_subp) (-> `Nsubpout') el PONDERADO, al reves de la asuncion
    * original.
    return matrix n_ponderado     = `Nsubpout'
    return matrix n_sin_ponderar  = `Nout'
    * Limite inferior/superior tal como los devolvio Stata en r(table)
    * (ver nota mas arriba de por que no se reconstruyen a mano).
    return matrix ci_lower        = `LIout'
    return matrix ci_upper        = `LSout'
    * Letras CLD y nombre de cada categoria de over(), expuestos
    * explicitamente via return local -- antes solo vivian en locales
    * internos (leidos de r(letra_N) puesto por Mata via st_global) que
    * se usaban para el "di" de pantalla, pero nunca se devolvian
    * formalmente. Sin esto, un pipeline que llama a svylet y despues
    * necesita las letras (no solo verlas en pantalla) no tenia forma
    * confiable de leerlas de vuelta.
    forvalues i = 1/`k' {
        local nom_grupo_`i' : word `i' of `niveles_over'
        return local letra_`i' "`letra_`i''"
        return local nombre_categoria_`i' "`nom_grupo_`i''"
    }

    * -- v1.4 -- ref(): p-valores "vs categoria base", crudo y ajustado
    * por Bonferroni (k-1 comparaciones) -- vacios (missing) si ref() no
    * se especifico, y en la posicion de la propia categoria base.
    return scalar ref_idx        = `ref_idx'
    return matrix p_vsref        = `p_vsref'
    return matrix p_vsref_raw    = `p_vsref_raw'
end

* -------------------------------------------------------------------------
* Subrutina compartida: selecciona de una matriz e(b)/e(V) las k filas que
* corresponden a las categorias de over() (para proportion, solo las del
* valor `level' del indicador). Usada tanto por el calculo analitico como
* por cada replica del bootstrap -- un solo lugar donde vive esta logica.
* -------------------------------------------------------------------------
capture program drop _svylet_seleccionar
program define _svylet_seleccionar, rclass
    args Bmat Vmat k level stat varname bout Vout

    local idx_usar ""
    if "`stat'" == "proportion" {
        * Para proportion, e(b) tiene una ECUACION por cada valor de
        * `varname', y DENTRO de cada ecuacion, una columna por categoria
        * de over(). El NOMBRE de la ecuacion depende de si `varname'
        * tiene value label definido -- CONFIRMADO con datos reales las
        * DOS formas (sysuse auto, foreign 0/1):
        *   - CON value label: el texto del label (ej. "Domestic"/"Foreign").
        *   - SIN value label: "_prop_N", donde N es la POSICION ORDINAL
        *     del valor entre los valores distintos de la variable (no el
        *     valor en si) -- ej. con foreign_sinlabel (valores 0 y 1),
        *     la ecuacion del valor 1 se llama "_prop_2" (segundo valor
        *     en orden), NO "_prop_1" ni "1". El intento anterior de usar
        *     el codigo numerico plano ("1") tambien fallo -- confirmado
        *     que esa forma no existe en ningun caso real probado.
        * Se filtra por `coleq', nunca por `colnames' ni por un patron
        * de texto con "@" -- eso resulto ser solo como "matrix list"
        * pinta dos niveles juntos en pantalla, no texto real.
        local vallab : value label `varname'
        local buscado_label ""
        if "`vallab'" != "" {
            local buscado_label : label `vallab' `level'
        }

        * Posicion ordinal del valor `level' entre los valores distintos
        * de `varname' (para armar "_prop_N" en el caso sin value label).
        quietly levelsof `varname', local(valores_x)
        local rank = 0
        local i = 0
        foreach v of local valores_x {
            local i = `i' + 1
            if `v' == `level' {
                local rank = `i'
            }
        }
        local buscado_ordinal = "_prop_`rank'"
        local buscado_numerico = "`level'"

        local ecuaciones : coleq `Bmat'
        local i = 0
        foreach eq of local ecuaciones {
            local i = `i' + 1
            if "`eq'" == "`buscado_label'" | "`eq'" == "`buscado_ordinal'" ///
                | "`eq'" == "`buscado_numerico'" {
                local idx_usar "`idx_usar' `i'"
            }
        }
        if wordcount("`idx_usar'") != `k' {
            local nombres : colnames `Bmat'
            di as err "no se pudieron identificar exactamente " `k' ///
                " filas para el valor " "`level'" " de " "`varname'" " en e(b)."
            di as err "  buscado como texto de label: " "`buscado_label'"
            di as err "  buscado como ordinal sin label: " "`buscado_ordinal'"
            di as err "  buscado como codigo numerico plano: " "`buscado_numerico'"
            di as err "  filas encontradas: " wordcount("`idx_usar'") " de " `k' " esperadas"
            di as err "  ecuacion / nombre de columna reales en e(b):"
            local i = 0
            foreach eq of local ecuaciones {
                local i = `i' + 1
                local nom : word `i' of `nombres'
                di as err "    [" `i' "] eq=" "`eq'" "  col=" "`nom'"
            }
            di as err "  ajustar level() si el valor de exito no es 1."
            exit 498
        }
    }
    else {
        forvalues i = 1/`k' {
            local idx_usar "`idx_usar' `i'"
        }
    }

    matrix `bout' = J(`k', 1, .)
    matrix `Vout' = J(`k', `k', 0)
    local fila = 0
    foreach i of local idx_usar {
        local fila = `fila' + 1
        matrix `bout'[`fila', 1] = `Bmat'[1, `i']
        local col = 0
        foreach j of local idx_usar {
            local col = `col' + 1
            matrix `Vout'[`fila', `col'] = `Vmat'[`i', `j']
        }
    }
    * Exponer los indices de columna usados -- el caller los necesita
    * para extraer en el mismo orden otras cantidades por categoria que
    * no pasan por esta subrutina (ej. e(_N), e(_N_subp)).
    return local idx_usar "`idx_usar'"
end

version 14
mata:
mata clear

// Version liviana: solo el F omnibus, sin Bonferroni/CLD -- es todo lo que
// necesita cada replica del bootstrap (recalcular CLD 1000 veces seria
// trabajo de sobra, ya que boot() solo recalibra el F omnibus, no los pares).
real scalar svylet_fomnibus_only(real matrix b, real matrix V, real scalar df, real scalar k)
{
    real matrix R, RVR
    real scalar i, stat_wald, k_dim, df_ajustado

    k_dim = k - 1
    R = J(k_dim, k, 0)
    for (i=2; i<=k; i++) {
        R[i-1,1] = -1
        R[i-1,i] = 1
    }
    RVR = R*V*R'
    stat_wald = (R*b)' * lusolve(RVR, R*b)
    df_ajustado = df - k_dim + 1
    return((df_ajustado / (k_dim * df)) * stat_wald)
}

// Agregacion final del boot(): sin dataset, sin archivo, sin frame -- opera
// directo sobre la matriz de Stata que acumulo los F de cada replica
// (missing "." en las filas de replicas que fallaron). Mismo criterio de
// siempre: p = (# replicas con F_replica >= F_obs, +1) / (B_efectivo+1),
// corregido por continuidad (Davison & Hinkley 1997).
void svylet_agregar_boot(string scalar matname, real scalar F_obs)
{
    real matrix todas, validas
    real scalar B_efectivo, cruzan

    todas = st_matrix(matname)
    validas = select(todas, todas :!= .)
    B_efectivo = rows(validas)

    st_numscalar("r(B_efectivo)", B_efectivo)
    if (B_efectivo > 0) {
        cruzan = sum(validas :>= F_obs)
        st_numscalar("r(p_boot)", (cruzan + 1) / (B_efectivo + 1))
    }
}

void svylet_core(string scalar bname, string scalar Vname, real scalar df,
                  real scalar alpha, real scalar k)
{
    real matrix b, V, R, RVR, Pmat
    real scalar stat_wald, Fstat, p_omni, i, j, npares, t, se, p_raw, p_adj
    string matrix letras
    real matrix grupos_bin
    real scalar ngrupos, g, tam, ok, a, incluido

    b = st_matrix(bname)
    V = st_matrix(Vname)

    // ---- Deteccion de varianza degenerada/faltante (proporciones en 0
    // o 1 exactos, sin variabilidad -- comun en dominios chicos) -----------
    // Sin esto, un V[i,i] faltante se propaga silenciosamente y Mata trata
    // "." como "mas grande que cualquier numero" en comparaciones, lo que
    // puede terminar agrupando una categoria como "no significativamente
    // distinta" SIN que el test se haya podido calcular en realidad -- eso
    // es un resultado enganoso, no un resultado conservador.
    real colvector var_degenerada
    var_degenerada = J(k,1,0)
    for (i=1; i<=k; i++) {
        if (missing(V[i,i]) | V[i,i] <= 0) var_degenerada[i] = 1
    }
    if (sum(var_degenerada) > 0) {
        printf("{txt}Aviso: varianza no definida (o cero) en %g de %g categorias" +
               " -- probablemente una proporcion exactamente 0 o 1 sin variabilidad.\n",
               sum(var_degenerada), k)
        printf("{txt}  Categorias afectadas (indice dentro de over(), no el valor): ")
        for (i=1; i<=k; i++) {
            if (var_degenerada[i]==1) printf("%g ", i)
        }
        printf("\n{txt}  Sus letras se muestran como [?] -- NO tratar como" +
               " sin diferencia, el test no se pudo calcular ahi.\n")
    }

    // ---- F omnibus: H0 b1=b2=...=bk, contraste contra el grupo 1 ----
    // (matematicamente invariante al grupo de referencia elegido -- ver
    // nota v1.1 en el encabezado del .ado sobre por que refgroup() se
    // elimino en vez de "implementarse": cualquier conjunto de k-1
    // contrastes independientes que genere el mismo subespacio da
    // identico Wald).
    // Ajuste de Korn & Graubard (1990), que es el DEFAULT de Stata para
    // testparm/test despues de un comando svy: -- confirmado contra
    // datos reales (sysuse auto): sin este ajuste, svylet daba F(4,67)
    // mientras que "svy: regress + testparm" daba F(4,64) para el MISMO
    // dato y diseno. La formula (ver "help test", opcion nosvyadjust):
    //   sin ajuste (lo que Stata llama nosvyadjust): W/k ~ F(k,d)
    //   CON ajuste (el default real de Stata):
    //     (d-k+1)/(k*d) * W ~ F(k, d-k+1)
    // con k = dimension del test (aqui, num_grupos-1) y d = e(df_r).
    //
    // v1.6 -- FIX: hasta v1.5, R se armaba con las k categorias
    // COMPLETAS (incluyendo las de varianza degenerada), asi que UNA
    // sola categoria con proporcion exactamente 0 o 1 (sin variabilidad)
    // hacia que TODO el F omnibus quedara en missing -- el "." de esa
    // categoria se propaga por R*V*R' y lusolve() devuelve missing,
    // aunque el resto de las categorias tuvieran variables perfectamente
    // calculables. Confirmado en produccion (tsvy, tabulado real):
    // filas con 4 de 4 anios con "Estimacion" (ninguno realmente
    // faltante) igual daban F en blanco porque 1-3 de esos anios eran
    // exactamente 0%. Igual que ya se hacia para los pares (Pmat, mas
    // abajo: una categoria degenerada solo invalida SUS pares, no todos)
    // ahora el F omnibus se calcula sobre el SUBCONJUNTO de categorias
    // con varianza definida, excluyendo las degeneradas del contraste en
    // vez de dejar que su "." contamine todo. Minimo 2 categorias
    // validas para poder testear una igualdad (con 0 o 1, no hay nada
    // que comparar) -- con 0 degeneradas, idx_validas es 1..k y el
    // resultado es identico al de v1.5.
    real scalar k_dim, df_ajustado, k_valida
    real colvector idx_validas, b_valida
    real matrix V_valida
    idx_validas = select((1::k), var_degenerada :== 0)
    k_valida = rows(idx_validas)
    if (k_valida >= 2) {
        b_valida = b[idx_validas]
        V_valida = V[idx_validas, idx_validas]
        k_dim = k_valida - 1
        R = J(k_dim, k_valida, 0)
        for (i=2; i<=k_valida; i++) {
            R[i-1,1] = -1
            R[i-1,i] = 1
        }
        RVR = R*V_valida*R'
        stat_wald = (R*b_valida)' * lusolve(RVR, R*b_valida)
        df_ajustado = df - k_dim + 1
        Fstat = (df_ajustado / (k_dim * df)) * stat_wald
        p_omni = 1 - F(k_dim, df_ajustado, Fstat)
        if (k_valida < k) {
            printf("{txt}  F omnibus calculado con %g de %g categorias" +
                   " (las de varianza degenerada quedan afuera del" +
                   " contraste, no del reporte -- sus puntos siguen" +
                   " listados, su letra queda en [?]).\n", k_valida, k)
        }
    }
    else {
        k_dim = .
        df_ajustado = .
        Fstat = .
        p_omni = .
        printf("{txt}  F omnibus no calculable: menos de 2 categorias" +
               " con varianza definida (%g de %g).\n", k_valida, k)
    }

    // ---- Pares con correccion Bonferroni ----
    npares = k*(k-1)/2
    Pmat = J(k,k,.)
    for (i=1; i<=k-1; i++) {
        for (j=i+1; j<=k; j++) {
            if (var_degenerada[i]==1 | var_degenerada[j]==1) {
                // No calculable -- se deja explicitamente en missing, NO
                // se fuerza a 1 ni se deja que la aritmetica con "." lo
                // decida por accidente.
                p_adj = .
                Pmat[i,j] = p_adj
                Pmat[j,i] = p_adj
                continue
            }
            se = sqrt(V[i,i] + V[j,j] - 2*V[i,j])
            if (se > 0) {
                t = (b[j]-b[i]) / se
                p_raw = 2*ttail(df, abs(t))   // ttail = cola superior de t(df)
                p_adj = min((p_raw*npares, 1))
            }
            else {
                p_adj = 1
            }
            Pmat[i,j] = p_adj
            Pmat[j,i] = p_adj
        }
    }

    // ---- Compact Letter Display: enumerar subconjuntos maximales
    //      mutuamente no-significativos (misma logica que la version
    //      Python ya validada) ----
    string colvector letra_grupo
    letra_grupo = J(k,1,"")
    real matrix candidatos
    real scalar nc, es_maximal
    string scalar letra_actual
    real scalar letra_idx

    // enumerar todos los subconjuntos no vacios (k<=~8 en la practica, 2^k manejable)
    real scalar total_subsets, s, bitpos
    real colvector miembros
    real matrix subsets_validos
    subsets_validos = J(0, k, .)

    total_subsets = 2^k - 1
    for (s=total_subsets; s>=1; s--) {
        miembros = J(k,1,0)
        for (bitpos=1; bitpos<=k; bitpos++) {
            if ( mod(floor(s / 2^(bitpos-1)), 2) == 1 ) miembros[bitpos] = 1
        }
        if (sum(miembros) < 1) continue
        // Nunca incluir categorias de varianza degenerada en ningun
        // subconjunto candidato -- no pueden ser "iguales" NI "distintas"
        // a nadie si el test no se pudo calcular ahi.
        ok = 1
        for (i=1; i<=k & ok; i++) {
            if (miembros[i]==1 & var_degenerada[i]==1) {
                ok = 0
                break
            }
        }
        if (ok==0) continue
        for (i=1; i<=k & ok; i++) {
            if (miembros[i]==0) continue
            for (j=i+1; j<=k; j++) {
                if (miembros[j]==0) continue
                if (Pmat[i,j] < alpha) {
                    ok = 0
                    break
                }
            }
        }
        if (ok==1) subsets_validos = subsets_validos \ miembros'
    }

    // quedarse solo con los MAXIMALES (no contenidos en otro valido)
    real matrix maximales
    maximales = J(0,k,.)
    real scalar contenido
    for (s=1; s<=rows(subsets_validos); s++) {
        contenido = 0
        for (i=1; i<=rows(subsets_validos); i++) {
            if (i==s) continue
            if (all(subsets_validos[s,] :<= subsets_validos[i,]) &
                sum(subsets_validos[s,]) < sum(subsets_validos[i,])) {
                contenido = 1
                break
            }
        }
        if (contenido==0) maximales = maximales \ subsets_validos[s,]
    }

    letra_idx = 0
    for (g=1; g<=rows(maximales); g++) {
        letra_idx++
        letra_actual = char(96 + letra_idx)   // 97='a'
        for (i=1; i<=k; i++) {
            if (maximales[g,i]==1) letra_grupo[i] = letra_grupo[i] + letra_actual
        }
    }

    // Categorias de varianza degenerada: marcar explicitamente con "?",
    // nunca dejarlas con una letra real ni vacias sin explicacion.
    for (i=1; i<=k; i++) {
        if (var_degenerada[i]==1) letra_grupo[i] = "?"
    }

    st_numscalar("r(F_omnibus)", Fstat)
    st_numscalar("r(p_omnibus)", p_omni)
    st_numscalar("r(n_pares)", npares)
    st_matrix("r(b)", b)
    for (i=1; i<=k; i++) {
        st_global("r(letra_"+strofreal(i)+")", letra_grupo[i])
    }
}

// v1.4 -- ref(): contrastes de Wald de UNA categoria contra una
// categoria BASE fija (ref_idx), k-1 comparaciones -- distinto del F
// omnibus + CLD de arriba (que compara TODOS los pares, k(k-1)/2
// comparaciones). Corrige por Bonferroni sobre k-1 (min(1, p_crudo *
// (k-1))), el mismo criterio que Dunn (1961, JASA) para "comparar
// varios grupos contra un control" -- MAS simple y MAS conservador que
// el metodo de un solo paso de Dunnett (1955, JASA; 1964, Biometrics),
// que usa la distribucion t multivariada para aprovechar que las k-1
// comparaciones comparten la misma categoria base y ganar potencia; no
// implementado aca (ver help para el detalle y las referencias).
void svylet_vsref(string scalar bname, string scalar Vname, real scalar df,
                   real scalar ref_idx, real scalar k,
                   string scalar pname, string scalar prawname)
{
    real matrix b, V
    real colvector p_adj, p_raw
    real scalar i, se, t, ncomp

    b = st_matrix(bname)
    V = st_matrix(Vname)
    p_adj = J(k,1,.)
    p_raw = J(k,1,.)
    ncomp = k - 1

    for (i=1; i<=k; i++) {
        if (i == ref_idx) continue
        if (missing(V[i,i]) | missing(V[ref_idx,ref_idx]) | V[i,i]<=0 | V[ref_idx,ref_idx]<=0) continue
        se = sqrt(V[i,i] + V[ref_idx,ref_idx] - 2*V[i,ref_idx])
        if (se<=0 | missing(se)) continue
        t = (b[i] - b[ref_idx]) / se
        p_raw[i] = 2*ttail(df, abs(t))
        p_adj[i] = min((p_raw[i]*ncomp, 1))
    }
    st_matrix(pname, p_adj)
    st_matrix(prawname, p_raw)
}
end