*! version 2.0.0  17jul2024
program pwmc

    version 12.1

    if ( replay() ) {

        pwmc_display `0'
        exit

    }

    pwmc_estimate_and_display `0'

end


/*  _________________________________________________________________________
                                                                 estimate  */

program pwmc_estimate_and_display

    syntax varname(numeric) [ if ] [ in ] [ fweight ] ///
        , Over(varname numeric)                       /// required
    [                                                 ///
        MCOMPare(passthru)                            ///
        PROCedure(passthru)       /// synonym mcompare(); no longer documented
        noADJust       /// synonym for mcompare(noadjust)
        HC3                                           ///
        Welch                                         ///
        Level(cilevel)                                ///
        CIeffects                                     ///
        PVEffects                                     ///
        PValues       /// retained synonym for pveffects; no longer documented
        EFFects                                       ///
        VARLabels                                     ///
        VALLabels                                     ///
        CFORMAT(passthru)                             ///
        PFORMAT(passthru)                             ///
        SFORMAT(passthru)                             ///
        noTABle                                       ///
    ]

    return_stats `varlist' `if' `in' [`weight' `exp'] , over(`over')

    mata : pwmc_estimate()

    pwmc_display ,                                    ///
        `mcompare' `procedure' `adjust'               ///
        `hc3' `welch' level(`level')                  ///
        `cieffects' `pveffects' `pvalues' `effects'   ///
        `varlabels' `vallabels'                       ///
        `cformat' `pformat' `sformat'                 ///
        `table'

end


/*  _________________________________________________________________________
                                                                  display  */

program pwmc_display

    if ( !inlist("`r(cmd)'","pwmc","pwmci") ) {

        display as err "last command not pwmc"
        exit 301

    }

    if ("`r(cmd)'" == "pwmci") ///
        mata : pwmc_estimate()

    if ("`r(cmd2)'" != "pwmci") {

        local VARLabels VARlabels
        local VALLabels VALLabels

    }

    syntax                   ///
    [ ,                      ///
        MCOMPare(passthru)   ///
        PROCedure(passthru)  ///
        noADJust             ///
        HC3                  ///
        Welch                ///
        Level(cilevel)       ///
        CIeffects            ///
        PVEffects            ///
        PValues              ///
        EFFects              ///
        `VARLabels'          ///
        `VALLabels'          ///
        CFORMAT(string asis) ///
        PFORMAT(string asis) ///
        SFORMAT(string asis) ///
        noTABle              ///
    ]

    parse_mcompare , `mcompare' `procedure' `adjust'

    confirm_variable , `varlabels' `vallabels'

    mata : pwmc_mctables(`level',"`hc3'","`welch'","`mcompare'")

    if ( ("`table'"=="notable") | (!c(noisily)) ) ///
        exit

    return_effects , `pveffects' `pvalues' `cieffects' `effects'

    return_fmt c %9.7g `cformat'
    return_fmt p %5.3f `pformat'
    return_fmt s %8.2f `sformat'

    return_table_layout , `hc3' `adjust' `varlabels' `vallabels'

    tempname table

    pwmc_display_title

    pwmc_display_table_header , `adjust'

    if ("`adjust'" == "noadjust") {

        matrix `table' = r(table_vs)
        matrix `table' = `table'[1..6,1...]'

        pwmc_display_table `table' , noadjust

        if ("`r(mcmethod_vs)'" == "") ///
            exit

        display

        return_table_layout , `hc3' `varlabels' `vallabels'

        pwmc_display_table_header

    }

    matrix `table' = r(table_mc_d)

    pwmc_display_table `table'

end


program pwmc_display_title

    display

    if (r(k)<3) ///
        display as txt "Two-sample t test with unequal variances"
    else ///
        display as txt "Pairwise comparisons of means with unequal variances"

end


program pwmc_display_table_header

    syntax [ , noADJust ]

    line TT , `adjust'
    display as txt _col(`r(T_col_1_width)') ///
        "{c |}" `r(T_se_title2)' `r(T_p_title2)' `r(T_ci_title2)'
    display as txt %`=`r(T_col_1_width)'-2's `"`r(T_depvar_label)'"' _continue
    display as txt _col(`r(T_col_1_width)') ///
        "{c |}" `r(T_diff_title)' `r(T_se_title)' `r(T_t_title)' `r(T_p_title)' `r(T_ci_title)'

    line + , `adjust'

end

program pwmc_display_table

    syntax name(name = table) [ , noADJust ]

    if ("`adjust'" == "noadjust") ///
        local method noadjust
    else ///
        local method = strupper("`r(mcmethod_vs)'")

    local nmethods : word count `method'

    local kstar = r(ks)

    local row = 1

    if (`"`r(T_over_label)'"' != "") {

        display as txt %`=`r(T_col_1_width)'-2's ///
            `"`r(T_over_label)'"' _continue
        display as txt _col(`r(T_col_1_width)') "{c |}"

    }

    forvalues i = 1/`kstar' {

        display as txt %`=`r(T_col_1_width)'-3's `"`r(T_vs_label`i')'"' _continue
        display as txt _col(`r(T_col_1_width)') "{c |}" _continue

        *display as txt _col(`r(T_col_1_width)') "{c |}" _continue

        display as res `r(T_diff_fmt)' `table'[`row',1] _continue // Diff.
        display as res `r(T_se_fmt)'   `table'[`row',2] _continue // Std. Err.

        if ("`r(pveffects)'" == "pveffects") ///
            display as res `r(T_t_fmt)' `table'[`row',3] _continue // t

        foreach m of local method {

            display as res `BAR' _continue

            if ("`r(pveffects)'" == "pveffects") ///
                display as res `r(T_p_fmt)' `table'[`row',4] _continue // P>|t|

            if ("`r(cieffects)'" == "cieffects") {

                display as res ///
                    `r(T_ci_ll_fmt)' `table'[`row',5] /// CI ll
                    `r(T_ci_ul_fmt)' `table'[`row',6] /// CI ul
                _continue

            }

            if (`nmethods' > 1) ///
                display as txt _col(`=`r(T_col_last)'+2') %4s "(`m')" _continue

            display // end of line

            local ++row

            local BAR _col(`r(T_col_1_width)') "{c |}"

        }

        local BAR // void

    }

    line BT , `adjust'

    if (`nmethods' == 1) ///
        exit

    local C  Dunnett C
    local GH Games and Howell
    local T2 Tamhane T2

    display as txt "Key:" _continue
    foreach m of local method {
        display as txt _col(`=8-strlen("`m'")') "`m': ``m''"
    }

end


/*  _________________________________________________________________________
                                                                utilities  */

program return_stats , rclass

    syntax varlist [ if ] [ in ] [ fweight ] , OVER(string)

    marksample touse
    quietly count if `touse'
    if ( !r(N) ) ///
        error 2000

    capture assert `over' == trunc(`over') if `touse' , fast
    if ( _rc ) {

        display as err "`over' may not contain noninteger values"
        exit 498

    }

    tempname L
    quietly tabulate `over' if `touse' , matrow(`L')
    local k = r(r)
    if (`k' < 2) {

        display as err "`over' has too few levels"
        exit 498

    }

    tempname stats
    matrix `stats' = J(3,`k',.)
    matrix rownames `stats' = mean sd n
    forvalues i = 1/`k' {

        quietly summarize `varlist' ///
            if `touse' & (`over' == `L'[`i',1]) [`weight' `exp']

        matrix `stats'[1,`i'] = r(mean)
        matrix `stats'[2,`i'] = r(sd)
        matrix `stats'[3,`i'] = r(N)

    }

    return visible scalar ks          = `k'*(`k'-1)/2
    return visible scalar k           = `k'
    return visible local  wtype         "`weight'"
    return visible local  wexp          `"`exp'"'
    return visible local  over          `over'
    return visible local  depvar        `varlist'
    return hidden  matrix levels_over = `L'
    return hidden  matrix stats       = `stats'

end


program parse_mcompare

    syntax                     ///
    [ ,                        ///
        MCOMPare(string asis)  ///
        PROCedure(string asis) ///
        noADJust               ///
    ]

    if (`"`mcompare'`procedure'"' != "") {

        parse_mcompare_procedure , mcompare(`mcompare') procedure(`procedure')

        local posof_noadjust : list posof "noadjust" in mcompare
        if ( `posof_noadjust' ) {

            local mcompare : subinstr local mcompare "noadjust" ""
            local mcompare `mcompare' // strip whitespaces
            local adjust noadjust

        }

        if ("`procedure'" != "") ///
            local option_name procedure()
        else ///
            local option_name mcompare()

    }

    if ( ("`mcompare'"=="") & ("`adjust'"=="") ) ///
        local mcompare "c gh t2" // default

    if (r(k) < 3) {

        if ("`option_name'" != "")  ///
            display as txt "too few levels; option `option_name' ignored"

        local mcompare // void
        local adjust noadjust

    }

    c_local adjust   : copy local adjust
    c_local mcompare : copy local mcompare

end


program parse_mcompare_procedure

    syntax                     ///
    [ ,                        ///
        MCOMPare(string asis)  ///
        PROCedure(string asis) ///
    ]

    if (`"`procedure'"' != "") {

        if ( !inlist(`"`mcompare'"',"",`"`procedure'"') ) {

            display as err "invalid option procedure()"
            exit 198

        }

        local mcompare : copy local procedure
        local method procedure

    }

    foreach m of local mcompare {

        if ( !inlist(strlower(`"`m'"'),"c","gh","t2") ) {

            local 0 , `m'
            capture syntax , NOADJust
            if ( _rc ) {

                if (`"`procedure'"' == "") {

                    display as err "option mcompare() invalid"
                    local method method

                }

                display as err "unknown `method' `m'"
                exit 198

            }

            local m noadjust

        }

        local mcmethod_vs `mcmethod_vs' `= strlower("`m'")'

    }

    c_local mcompare : list uniq mcmethod_vs

end


program confirm_variable

    syntax        ///
    [ ,           ///
        VARLabels ///
        VALLabels ///
    ]

    if ("`varlabels'" == "varlabels") {

        novarabbrev confirm variable `r(depvar)'
        novarabbrev confirm variable `r(over)'

    }
    else if ("`vallabels'" == "vallabels") ///
        novarabbrev confirm variable `r(over)'

end


program return_effects

    syntax        ///
    [ ,           ///
        PVEFFECTS ///
        PVALUES   ///
        CIEFFECTS ///
        EFFECTS   ///
    ]

    if ("`effects'" == "effects") {

        local cieffects cieffects
        local pveffects pveffects

    }
    else {

        if ("`pvalues'" == "pvalues") ///
            local pveffects pveffects
        else if ("`pveffects'" == "") ///
            local cieffects cieffects

    }

    mata : st_global("r(pveffects)",st_local("pveffects"),"hidden")
    mata : st_global("r(cieffects)",st_local("cieffects"),"hidden")

end


program return_fmt

    args letter fmt fmt2

    if ("`c(`letter'format)'" != "") ///
        local fmt `c(`letter'format)'

    if ("`fmt2'" != "") {

        if (fmtwidth(`"`fmt2'"') <= fmtwidth("`fmt'")) ///
            local fmt `fmt2'
        else ///
            display as txt "note: invalid `letter'ormat(), using default"

    }

    mata : st_global("r(`letter'fmt)","`fmt'","hidden")

end


program return_table_layout

    syntax        ///
    [ ,           ///
        HC3       ///
        noADJust  ///
        VARLabels ///
        VALLabels ///
    ]

    if ("`adjust'" == "noadjust") ///
        local method noadjust
    else ///
        local method `r(mcmethod_vs)'

    local nmethods : word count `method'

    return_labels , `varlabels' `vallabels'

    local diff_title _col(`=`r(T_col_1_width)'+8') "Diff."
    local diff_fmt   _col(`=`r(T_col_1_width)'+3') `r(cfmt)'

    local se_title _col(`=`r(T_col_1_width)'+15') "Std. Err."
    local se_fmt   _col(`=`r(T_col_1_width)'+14') `r(cfmt)'

    if ("`hc3'" == "hc3") ///
        local se_title2 _col(`=`r(T_col_1_width)'+14') "Robust HC3"
    else ///
        local se_title2 _col(`=`r(T_col_1_width)'+16') "Robust"

    if ("`r(pveffects)'" == "pveffects") {

        local t_title _col(`=`r(T_col_1_width)'+30') "t"
        local t_fmt   _col(`=`r(T_col_1_width)'+24') `r(sfmt)'

        local p_title _col(`=`r(T_col_1_width)'+35') "P>|t|"
        local p_fmt   _col(`=`r(T_col_1_width)'+35') `r(pfmt)'

        if (`nmethods' == 1) {

            if ("`method'" == "noadjust") ///
                local p_title2 _col(`=`r(T_col_1_width)'+30') "Unadjusted"
            else if ("`method'" == "c") ///
                local p_title2 _col(`=`r(T_col_1_width)'+29') "Dunnett's C"
            else if ("`method'" == "gh") ///
                local p_title2 _col(`=`r(T_col_1_width)'+26') "Games / Howell"
            else if ("`method'" == "t2") ///
                local p_title2 _col(`=`r(T_col_1_width)'+28') "Tamhane's T2"

        }

        local col_last = `r(T_col_1_width)' + 35 + fmtwidth("`r(pfmt)'")

    }

    if ("`r(cieffects)'" == "cieffects") {

        if ("`r(pveffects)'" == "pveffects") {

            if (`nmethods' == 1) {

                if ("`method'" == "noadjust") ///
                    local p_title2 // void
                else ///
                    local p_title2 _col(`=`r(T_col_1_width)'+32') "Adjusted"

            }
            else local p_title2 // void

            local col_ci_title = `r(T_col_1_width)' + 45 - (strlen("`r(level)'")-2)
            local ci_ll_fmt _col(`=`r(T_col_1_width)'+44') `r(cfmt)'
            local ci_ul_fmt _col(`=`r(T_col_1_width)'+56') `r(cfmt)'

            local col_last = `r(T_col_1_width)' + 56 + fmtwidth("`r(cfmt)'")

        }
        else {

            local col_ci_title = `r(T_col_1_width)' + 28 - (strlen("`r(level)'")-2)
            local ci_ll_fmt _col(`=`r(T_col_1_width)'+27') `r(cfmt)'
            local ci_ul_fmt _col(`=`r(T_col_1_width)'+39') `r(cfmt)'

            local col_last = `r(T_col_1_width)' + 39 + fmtwidth("`r(cfmt)'")

        }

        local ci_title  _col(`col_ci_title') "[`r(level)'% Conf. Interval]"

        if (`nmethods' == 1) {

            if ("`method'" == "noadjust") ///
                local ci_title2 _col(`=`col_ci_title'+6') "Unadjusted"
            else if ("`method'" == "c") ///
                local ci_title2 _col(`=`col_ci_title'+5') "Dunnett's C"
            else if ("`method'" == "gh") ///
                local ci_title2 _col(`=`col_ci_title'+2') "Games and Howell"
            else if ("`method'" == "t2") ///
                local ci_title2 _col(`=`col_ci_title'+4') "Tamhane's T2"

        }

    }

    mata : st_global("r(T_diff_title)", st_local("diff_title"),   "hidden")
    mata : st_global("r(T_se_title)",   st_local("se_title"),     "hidden")
    mata : st_global("r(T_t_title)",    st_local("t_title"),      "hidden")
    mata : st_global("r(T_p_title)",    st_local("p_title"),      "hidden")
    mata : st_global("r(T_ci_title)",   st_local("ci_title"),     "hidden")

    mata : st_global("r(T_se_title2)",  st_local("se_title2"),    "hidden")
    mata : st_global("r(T_p_title2)",   st_local("p_title2"),     "hidden")
    mata : st_global("r(T_ci_title2)",  st_local("ci_title2"),    "hidden")

    mata : st_global("r(T_diff_fmt)",   st_local("diff_fmt"),     "hidden")
    mata : st_global("r(T_se_fmt)",     st_local("se_fmt"),       "hidden")
    mata : st_global("r(T_t_fmt)",      st_local("t_fmt"),        "hidden")
    mata : st_global("r(T_p_fmt)",      st_local("p_fmt"),        "hidden")
    mata : st_global("r(T_ci_ll_fmt)",  st_local("ci_ll_fmt"),    "hidden")
    mata : st_global("r(T_ci_ul_fmt)",  st_local("ci_ul_fmt"),    "hidden")

    mata : st_global("r(T_col_last)",   st_local("col_last"),     "hidden")

end


program return_labels

    syntax        ///
    [ ,           ///
        VARLABELs ///
        VALLABELs ///
    ]

    if (c(stata_version) >= 14) ///
        local u u

    local col_1_width 14 // default

    foreach var in depvar over {

        if ("`varlabels'" == "varlabels") ///
            local `var'_label : variable label `r(`var')'

        if (`"``var'_label'"' == "") ///
            local `var'_label `r(`var')'

        mata : st_global("r(T_`var'_label)",st_local("`var'_label"),"hidden")

        local col_1_width = max(`col_1_width',`u'strlen(`"``var'_label'"')+2)
        mata : st_global("r(T_col_1_width)","`col_1_width'","hidden")

    }

    local vs_names `r(vs_names)'
    local i 0

    foreach vs of local vs_names {

        assert regexm("`vs'","^([0-9]+)vs([0-9]+)$")

        local vs_1 = regexs(1)
        local vs_2 = regexs(2)

        if ("`r(levels_over)'" == "matrix") {

            local vs_1 = el(r(levels_over),`vs_1',1)
            local vs_2 = el(r(levels_over),`vs_2',1)

            if ("`vallabels'" == "vallabels") {

                local vs_1 : label (`r(over)') `vs_1'
                local vs_2 : label (`r(over)') `vs_2'

            }

        }

        local label `"`vs_1' vs `vs_2'"'

        mata : st_global("r(T_vs_label`++i')",st_local("label"),"hidden")

        local col_1_width = max(`col_1_width',`u'strlen(`"`label'"')+3)
        mata : st_global("r(T_col_1_width)","`col_1_width'","hidden")

    }

end


program line

    syntax anything(name = separator) [ , noADJust ]
    args separator

    local w_1 = `r(T_col_1_width)' - 1
    local w_2 = `r(T_col_last)' - `r(T_col_1_width)' - 1

    if ("`adjust'" != "noadjust") ///
        if (`: word count `r(mcmethod_vs)'' > 1) ///
            local w_2 = `w_2' + 6

    display as txt "{hline `w_1'}{c `separator'}{hline `w_2'}"

end


/*  _________________________________________________________________________
                                                                     Mata  */

version 12.1


mata :


mata set matastrict   on
mata set mataoptimize on


class pwmc_ado
{
    public :

        void get_st_r()
        void estimate()
        void set_st_r()
        void st_r_mat()

    protected :

        void vs_names()

        real   scalar    k
        real   scalar    kstar
        real   rowvector xbar
        real   rowvector sd
        real   rowvector n

        real   rowvector b_vs
        real   rowvector se_vs
        real   rowvector df_vs
        real   rowvector se_hc3_vs
        real   rowvector df_welch_vs

        string matrix    vs_names
}


void pwmc_ado::get_st_r()
{
    k     = st_numscalar("r(k)")
    kstar = st_numscalar("r(ks)")
    xbar  = st_matrix("r(stats)")[1,]
    sd    = st_matrix("r(stats)")[2,]
    n     = st_matrix("r(stats)")[3,]
}


void pwmc_ado::estimate()
{
    real rowvector nm1, np1, sd2n, sd2nm1
    real scalar    i, j
    real rowvector jj


    nm1    = n:-1
    np1    = n:+1
    sd2n   = (sd:^2):/n   // equivalent to vce(hc2)
    sd2nm1 = (sd:^2):/nm1 // equivalent to vce(hc3)

    b_vs        = J(1,kstar,.z)
    se_vs       = J(1,kstar,.z)
    df_vs       = J(1,kstar,.z)
    se_hc3_vs   = J(1,kstar,.z)
    df_welch_vs = J(1,kstar,.z)

    for (i=j=1; i<k; i++) {

        jj = (j..(k-i)+(j-1))

        b_vs[jj] = xbar[|i+1\ k|] :- xbar[i]

        se_vs[jj] = sqrt(sd2n[|i+1\ k|] :+ sd2n[i])

        df_vs[jj] =
            (se_vs[jj]:^4) :/
            (
                (sd[|i+1\ k|]:^4 :/ (n[|i+1\ k|]:^2:*nm1[|i+1\ k|]))
                :+ (sd[i]:^4:/(n[i]:^2*nm1[i]))
            )

        se_hc3_vs[jj] = sqrt(sd2nm1[|i+1\ k|] :+ sd2nm1[i])

        df_welch_vs[jj] =
            (se_vs[jj]:^4) :/
            (
                (sd[|i+1\ k|]:^4 :/ (n[|i+1\ k|]:^2:*np1[|i+1\ k|]))
                :+ (sd[i]:^4:/(n[i]:^2*np1[i]))
            ) :- 2

        j = j + (k-i)

    }
}


void pwmc_ado::set_st_r()
{
    vs_names()

    st_global("r(cmd2)", st_global("r(cmd)"))
    st_global("r(cmd)", "pwmc")

    st_r_mat("b_vs",        b_vs,        "", vs_names, "hidden")
    st_r_mat("se_vs",       se_vs,       "", vs_names, "hidden")
    st_r_mat("df_vs",       df_vs,       "", vs_names, "hidden")
    st_r_mat("se_hc3_vs",   se_hc3_vs,   "", vs_names, "hidden")
    st_r_mat("df_welch_vs", df_welch_vs, "", vs_names, "hidden")
}


void pwmc_ado::vs_names()
{
    string colvector vs_level
    real   scalar    i, j
    real   rowvector jj


    vs_names = J(kstar,1,"")
    vs_level = strofreal((1::k))

    for (i=j=1; i<k; i++) {

        jj = (j::(k-i)+(j-1))
        vs_names[jj] = vs_level[i+1::k]:+"vs":+vs_level[i]:+vs_names[jj]
        j  = j + (k-i)

    }

    if (st_global("r(over)") != "")
        vs_names = vs_names:+("."+st_global("r(over)"))

    vs_names = (J(kstar,1,""),vs_names)
}


void pwmc_ado::st_r_mat(

    string scalar matname,
    real   matrix Matrix,
  | string matrix rownames,
    string matrix colnames,
    string scalar hcat

    )
{
    string scalar rmatname


    rmatname = sprintf("r(%s)",matname)

    if (args() < 5)
        st_matrix(rmatname,Matrix)
    else
        st_matrix(rmatname,Matrix,hcat)

    if ( (args()<3) | (length(Matrix)==0) )
        return

    if (rownames != "")
        st_matrixrowstripe(rmatname,rownames)

    if (colnames != "")
        st_matrixcolstripe(rmatname,colnames)
}


    /*  _________________________________  tables  */

class pwmc_tables extends pwmc_ado
{
    public :

        void ci_level()
        void table_vs()
        void mcompare()
        void set_st_r()

    private :

        real   matrix dunnett_c()

        real   scalar alpha

        real   matrix table_vs
        real   matrix table_mc   // by methods
        real   matrix table_mc_d // by differences

        real   matrix pvalue_adj
        real   matrix crit_adj
        real   matrix ll_adj
        real   matrix ul_adj

        string matrix mc_names
}


void pwmc_tables::ci_level(real scalar level) alpha = level/100


void pwmc_tables::table_vs(

    real scalar hc3,
    real scalar welch

    )
{
    table_vs = J(9,kstar,0)

    table_vs[1,] = st_matrix("r(b_vs)")
    table_vs[2,] = st_matrix((hc3   ? "r(se_hc3_vs)"   : "r(se_vs)"))
    table_vs[7,] = st_matrix((welch ? "r(df_welch_vs)" : "r(df_vs)"))

    table_vs[3,] = table_vs[1,]:/table_vs[2,]
    table_vs[4,] = 2*ttail(table_vs[7,],abs(table_vs[3,]))
    table_vs[8,] = invttail(table_vs[7,],(1-alpha)/2)
    table_vs[5,] = table_vs[1,]:-table_vs[8,]:*table_vs[2,]
    table_vs[6,] = table_vs[1,]:+table_vs[8,]:*table_vs[2,]
}


void pwmc_tables::mcompare(string rowvector mcompare)
{
    real rowvector df
    real colvector row


    if ( !cols(mcompare) )
        return

    pvalue_adj = crit_adj = J(cols(mcompare),kstar,.z)

    df = table_vs[7,]

    if ( any(row=select((1..cols(mcompare)), (mcompare:=="c"))) ) {
       crit_adj[row,] = dunnett_c(invtukeyprob(k,(n:-1),alpha) :* (sd:^2):/n)
    }

    if ( any(row=select((1..cols(mcompare)), (mcompare:=="gh"))) ) {
        pvalue_adj[row,] = 1:-tukeyprob(k,df,abs(table_vs[3,])*sqrt(2))
        pvalue_adj[row,] = colmin(J(1, kstar, 1)\ pvalue_adj[row,])
        crit_adj[row,] = invtukeyprob(k,df,alpha)/sqrt(2)
    }

    if ( any(row=select((1..cols(mcompare)), (mcompare:=="t2"))) ) {
        pvalue_adj[row,] = 1 :- (1:-table_vs[4,]):^kstar
        pvalue_adj[row,] = colmin(J(1, kstar, 1)\ pvalue_adj[row,])
        crit_adj[row,] = invttail(df,(1-alpha^(1/kstar))/2)
    }

    ll_adj = table_vs[1,]:-crit_adj:*table_vs[2,]
    ul_adj = table_vs[1,]:+crit_adj:*table_vs[2,]

    table_mc = J(cols(mcompare),1,table_vs[(1::3),]')
    table_mc = (table_mc,vec(pvalue_adj'),vec(ll_adj'),vec(ul_adj'))

    table_mc_d = colshape(J(1,cols(mcompare),table_vs[(1::3),]'),3)
    table_mc_d = (table_mc_d,vec(pvalue_adj),vec(ll_adj),vec(ul_adj))
}


real matrix pwmc_tables::dunnett_c(real rowvector SR)
{
    real rowvector V, C
    real rowvector jj
    real scalar    i, j


    V = table_vs[2,]:^2
    C = J(1, kstar, .z)

    for (i=j=1; i<k; i++) {

        jj = (j..(k-i)+(j-1))
        C[jj] = ((SR[|i+1\ k|]:+SR[i]):/V[jj])/sqrt(2)
        j  = j + (k-i)

    }

    return(C)
}


void pwmc_tables::set_st_r(string rowvector mcompare)
{
    real scalar i


    vs_names()
    mc_names = (J(cols(mcompare),1,""), mcompare')

    st_numscalar("r(level)", alpha*100)

    st_global("r(mcmethod_vs)", invtokens(mc_names[,2]'))

    st_r_mat("table_vs", table_vs,
        (J(9,1,""), ("b"\"se"\"t"\"pvalue"\"ll"\"ul"\"df"\"crit"\"eform")),
        vs_names
    )

    st_r_mat("table_mc", table_mc,
        (
         colshape(J(1,kstar,mc_names[,2]),1),
         J(cols(mcompare),1,vs_names[,2])
        ),
        (J(6,1,""),("b"\"se"\"t"\"pvalue_adj"\"ll_adj"\"ul_adj")),
        "hidden"
    )

    st_r_mat("table_mc_d", table_mc_d,
        (
         colshape(J(1,cols(mcompare),vs_names[,2]),1),
         J(kstar,1,mc_names[,2])
        ),
        (J(6,1,""),("b"\"se"\"t"\"pvalue_adj"\"ll_adj"\"ul_adj")),
        "hidden"
    )

    st_global("r(procedure)", st_global("r(mcmethod_vs)"), "hidden")

    st_r_mat("t",     table_vs[3,],    "", vs_names, "hidden")
    st_r_mat("nuhat", table_vs[7,],    "", vs_names, "hidden")
    st_r_mat("Var",   table_vs[2,]:^2, "", vs_names, "hidden")
    st_r_mat("diff",  table_vs[1,],    "", vs_names, "hidden")

    st_r_mat("A",     crit_adj',   vs_names, mc_names, "hidden")
    st_r_mat("p_adj", pvalue_adj', vs_names, mc_names, "hidden")

    st_r_mat("ci", (vec(ll_adj'), vec(ul_adj')),
        (
         colshape(J(1,kstar,mc_names[,2]), 1),
         J(cols(mcompare),1,vs_names[,2])
        ),
        ((""\""), ("ll"\"ul")),
        "hidden"
    )

    for (i=1; i<=rows(mc_names); i++) {

        st_r_mat("A_"+mc_names[i, 2],  crit_adj[i,], "", vs_names, "hidden")
        st_r_mat("ll_"+mc_names[i, 2], ll_adj[i,],   "", vs_names, "hidden")
        st_r_mat("ul_"+mc_names[i, 2], ul_adj[i,],   "", vs_names, "hidden")

    }

    if (st_global("r(over)") != "")
        st_global(
            "r(vs_names)",
            invtokens(substr(vs_names[,2],1,strpos(vs_names[,2],"."):-1)'),
            "hidden"
            )
    else
        st_global("r(vs_names)",invtokens(vs_names[,2]'),"hidden")

}


/*  _________________________________________________________________________
                                                          entry point ado  */

void pwmc_estimate()
{
    class pwmc_ado scalar P


    P.get_st_r()

    P.estimate()

    P.set_st_r()
}


void pwmc_mctables(

    real   scalar level,
    string scalar hc3,
    string scalar welch,
    string scalar mcompare

    )
{
    class pwmc_tables scalar T


    T.get_st_r()

    T.ci_level(level)

    T.table_vs((hc3=="hc3"),(welch=="welch"))

    T.mcompare(tokens(mcompare))

    T.set_st_r(tokens(mcompare))
}


end


exit


/*  _________________________________________________________________________
                                                              version history

2.0.0   17jul2024   complete rewrite
                    new options -noadjust-, -hc3-, and -welch-
                    fweights allowed (not documented)
                    no longer requires external Mata function
                        (bug fix under Linux)
                    new output (again)
                    changed returned results; old results hidden
1.1.0   07jan2014   new external Mata function mPwmc.mo
                    calculate adjusted p-values
                    changed returned results (old results hidden)
                    new output (new code)
                    -replay()- results
                    new option -pvalues-
                    new options -pformat()-, -sformat()-, -notable-
                    option -mcompare- as synonym for -procedure-
1.0.0   28jan2013   first release on SSC