*! ml_sens v1.0.0  metaLong for Stata 14.1
*! Time-Varying Sensitivity Analysis via Longitudinal ITCV
*! Subir Hait, Michigan State University  |  haitsubi@msu.edu
*!
*! Computes the Impact Threshold for a Confounding Variable (ITCV) at each
*! follow-up time using pooled estimates from ml_meta.
*!
*! Reference: Frank, K.A. (2000). Sociological Methods & Research, 29(2), 147-194.
*!
*! Mathematical background
*! -----------------------
*! At time t, let theta_t = pooled effect, sy2_t = weighted variance of yi,
*! se_t = cluster-robust SE, df_t = degrees of freedom.
*!
*! Correlation-scale: r_t = theta_t / sqrt(theta_t^2 + sy2_t)
*! Raw ITCV:          itcv_t = sqrt(|r_t|)
*! Significance-adjusted:
*!   theta_star_t = |theta_t| - crit_t * se_t
*!   if theta_star_t > 0:
*!     r_star_t = theta_star_t / sqrt(theta_star_t^2 + sy2_t)
*!     itcv_alpha_t = sqrt(|r_star_t|)
*!   else: itcv_alpha_t = 0
*!
*! fragile at time t if itcv_alpha_t < delta

program define ml_sens, rclass
    version 14.1

    syntax varlist(min=2 max=2 numeric) [if] [in] , ///
        STUDy(varname)          ///
        TIME(varname numeric)   ///
        METAfile(string asis)   ///  path to saved ml_meta results dataset
        [                       ///
        ALpha(real 0.05)        ///
        DELTA(real 0.15)        ///  fragility threshold
        SAVing(string asis)     ///
        REPLACE                 ///
        ]

    tokenize `varlist'
    local yi `1'
    local vi `2'

    /* ---- mark sample ---- */
    marksample touse
    markout `touse' `time' `yi' `vi'
    quietly replace `touse' = 0 if `study' == ""

    /* ---- load meta results ---- */
    preserve

    quietly use `"`metafile'"', clear

    /* validate it's an ml_meta results file */
    capture confirm variable time theta se df tau2
    if _rc {
        di as error "metafile() does not look like ml_meta output (missing columns)."
        exit 198
    }

    /* extract meta results into matrix */
    local n_times = _N
    tempname meta_mat
    mkmat time k theta se df t_stat p_val ci_lb ci_ub tau2, matrix(`meta_mat')

    restore

    /* ---- result matrix: time theta se df sy r_effect itcv itcv_alpha fragile ---- */
    tempname sens_mat
    matrix `sens_mat' = J(`n_times', 9, .)
    matrix colnames `sens_mat' = time theta se df sy r_effect itcv itcv_alpha fragile

    /* ---- preserve and work on analysis data ---- */
    preserve
    quietly keep if `touse'
    quietly keep `yi' `vi' `study' `time'

    /* ================================================================== */
    /*  Loop over time points from meta results                            */
    /* ================================================================== */
    forvalues row = 1/`n_times' {

        local t_val   = `meta_mat'[`row', 1]
        local theta_t = `meta_mat'[`row', 3]
        local se_t    = `meta_mat'[`row', 4]
        local df_t    = `meta_mat'[`row', 5]
        local tau2_t  = `meta_mat'[`row', 10]

        matrix `sens_mat'[`row', 1] = `t_val'
        matrix `sens_mat'[`row', 2] = `theta_t'
        matrix `sens_mat'[`row', 3] = `se_t'
        matrix `sens_mat'[`row', 4] = `df_t'

        /* Skip if meta estimate is missing */
        if missing(`theta_t') | missing(`se_t') | missing(`df_t') continue

        /* ---- Subset to current time point ---- */
        tempvar in_t
        quietly gen byte `in_t' = (`time' == `t_val')

        /* ---- tau2-adjusted weights ---- */
        local tau2_use = max(0, `tau2_t')
        tempvar wi_RE
        quietly gen double `wi_RE' = 1 / (`vi' + `tau2_use') if `in_t'

        /* ---- Weighted variance of yi (sy2) ---- */
        quietly summ `wi_RE' if `in_t', meanonly
        local sw = r(sum)

        tempvar wdev2
        quietly gen double `wdev2' = `wi_RE' * (`yi' - `theta_t')^2 if `in_t'
        quietly summ `wdev2' if `in_t', meanonly
        local sy2 = r(sum) / `sw'
        local sy  = sqrt(`sy2')

        /* ---- Correlation-scale effect ---- */
        local denom = sqrt(`theta_t'^2 + `sy2')
        if `denom' <= 0 {
            drop `in_t' `wi_RE' `wdev2'
            continue
        }
        local r_effect = `theta_t' / `denom'

        /* ---- Raw ITCV ---- */
        local itcv = sqrt(abs(`r_effect'))

        /* ---- Significance-adjusted ITCV ---- */
        if `df_t' > 1e5 {
            /* z-based */
            local crit = invnormal(1 - `alpha' / 2)
        }
        else {
            local crit = invt(`df_t', 1 - `alpha' / 2)
        }

        local theta_star = abs(`theta_t') - `crit' * `se_t'

        if `theta_star' <= 0 {
            local itcv_alpha = 0
        }
        else {
            local denom_star = sqrt(`theta_star'^2 + `sy2')
            local r_star     = `theta_star' / `denom_star'
            local itcv_alpha = sqrt(abs(`r_star'))
        }

        /* fragile if itcv_alpha < delta */
        local fragile = (`itcv_alpha' < `delta')

        /* Store */
        matrix `sens_mat'[`row', 5] = `sy'
        matrix `sens_mat'[`row', 6] = `r_effect'
        matrix `sens_mat'[`row', 7] = `itcv'
        matrix `sens_mat'[`row', 8] = `itcv_alpha'
        matrix `sens_mat'[`row', 9] = `fragile'

        drop `in_t' `wi_RE' `wdev2'
    }

    /* ================================================================== */
    /*  Build results dataset                                              */
    /* ================================================================== */
    quietly {
        clear
        svmat double `sens_mat', names(col)

        label var time       "Follow-up time"
        label var theta      "Pooled effect"
        label var se         "Cluster-robust SE"
        label var df         "Degrees of freedom"
        label var sy         "Weighted SD of effects (sqrt of sy2)"
        label var r_effect   "Effect on correlation scale (r)"
        label var itcv       "Raw ITCV: threshold to nullify estimate"
        label var itcv_alpha "Adj. ITCV: threshold to lose significance"
        label var fragile    "1 if itcv_alpha < delta (`delta')"

        char _dta[ml_type]   "ml_sens"
        char _dta[ml_alpha]  "`alpha'"
        char _dta[ml_delta]  "`delta'"
    }

    /* ================================================================== */
    /*  Compute summary statistics                                         */
    /* ================================================================== */
    quietly summ itcv_alpha, meanonly
    local itcv_min  = r(min)
    local itcv_mean = r(mean)

    quietly count if fragile == 1
    local n_fragile = r(N)
    quietly count if !missing(itcv_alpha)
    local n_valid = r(N)
    local frag_prop = cond(`n_valid' > 0, `n_fragile' / `n_valid', .)

    /* ================================================================== */
    /*  Display                                                            */
    /* ================================================================== */
    di _newline
    di as txt "  {hline 72}"
    di as txt "  metaLong: Longitudinal Sensitivity Analysis (ITCV)"
    di as txt "  Fragility threshold (delta) = `delta'"
    di as txt "  {hline 72}"
    di as txt _col(4) "Time" _col(12) "theta" _col(21) "sy" ///
               _col(30) "r_effect" _col(40) "ITCV" _col(50) "ITCV_adj" ///
               _col(60) "Fragile?"
    di as txt "  {hline 72}"

    forvalues row = 1/`n_times' {
        local t_r   = `sens_mat'[`row', 1]
        local th_r  = `sens_mat'[`row', 2]
        local sy_r  = `sens_mat'[`row', 5]
        local re_r  = `sens_mat'[`row', 6]
        local iv_r  = `sens_mat'[`row', 7]
        local ia_r  = `sens_mat'[`row', 8]
        local fg_r  = `sens_mat'[`row', 9]

        if missing(`th_r') {
            di as txt %7.1f `t_r' "    (missing estimate)"
            continue
        }

        local frag_lbl = cond(`fg_r' == 1, "YES", "no")
        local frag_col = cond(`fg_r' == 1, "as err", "as res")

        di as txt  %7.1f `t_r' ///
           as res  %9.4f `th_r' %9.4f `sy_r' %9.4f `re_r' ///
                   %9.4f `iv_r' %9.4f `ia_r' " " ///
           `frag_col' "  `frag_lbl'"
    }
    di as txt "  {hline 72}"
    di as txt "  ITCV_adj min  = " as res %6.4f `itcv_min'
    di as txt "  ITCV_adj mean = " as res %6.4f `itcv_mean'
    di as txt "  Fragile prop  = " as res %6.4f `frag_prop' ///
               as txt "  (`n_fragile' / `n_valid' time points)"

    /* ================================================================== */
    /*  Save                                                               */
    /* ================================================================== */
    if `"`saving'"' != "" {
        if "`replace'" != "" quietly save `saving', replace
        else                  quietly save `saving'
        di as txt _newline "  Results saved to: " as res `"`saving'"'
    }

    restore

    /* ================================================================== */
    /*  Return                                                             */
    /* ================================================================== */
    return matrix sens      = `sens_mat'
    return scalar itcv_min  = `itcv_min'
    return scalar itcv_mean = `itcv_mean'
    return scalar frag_prop = `frag_prop'
    return scalar delta     = `delta'
    return scalar alpha     = `alpha'

end