*! version 1.1.0 15dec2022 JM. Domenech, R. Sesma /* Kruskal-Wallis equality-of-populations rank test and non-parametric pairwise comparisons across the levels of factor variables */ program define pwkwallis, rclass version 12 syntax varname(numeric) [if] [in], by(varname numeric) [BONferroni] tempname sumR2 N H S2 SE0 P tempvar R R2 if "`varlist'" == "`by'" { display in red "outcome variable can't be the group defining variable" exit 198 } local adj = ("`bonferroni'"!="") marksample touse //ifin marksample *Kruskal-Wallis kwallis `varlist' if `touse', by(`by') scalar `H' = r(chi2_adj) quietly { egen `R' = rank(`varlist') if `touse' //Sort answers generate `R2' = `R'^2 if `touse' summarize `R2' //R2 sum, N (number of cases) scalar `sumR2' = r(sum) scalar `N' = r(N) } preserve quietly collapse (sum) Ri = `R' (count) ni = `R' if `touse', by(`by') //Ri, ni with collapse quietly generate Rm = Ri/ni //Rm: mean Ri local k = _N //k: number of groups scalar `S2' = (1/(`N'-1))*(`sumR2' - (`N'*(`N'+1)^2)/4) scalar `SE0' = `S2'*(`N'-1-`H')/(`N'-`k') *Compute contrast, std err, t, p no adj pairwise quietly levelsof `by', local(levels) local ncomb = (`k'*(`k'-1))/2 //number of combinations of n elements chosing r: n!/(r!(n-r)!); if r=2, ncomb = (n*(n-1))/2 matrix `P' = J(`ncomb',4,.) local c 1 local rowlbls "" forvalues i = 1/`k' { local t = `i'+1 forvalues j = `t'/`k' { matrix `P'[`c',1] = Rm[`i'] - Rm[`j'] matrix `P'[`c',2] = sqrt(`SE0'*(1/ni[`i'] + 1/ni[`j'])) matrix `P'[`c',3] = abs(`P'[`c',1])/`P'[`c',2] matrix `P'[`c',4] = 2*ttail(`N'-`k',`P'[`c',3]) tempname lbl`c' local vali : word `i' of `levels' local valj : word `j' of `levels' scalar `lbl`c'' = "`vali' vs `valj'" local rowlbls = "`rowlbls' `vali'vs`valj'" local c = `c'+1 } } restore if (`adj') { matrix `P' = `P' , J(`ncomb',1,.) //add column for manual bonferroni correction, NOT using _mtest command forvalues i = 1/`ncomb' { matrix `P'[`i',5] = 1 - (1-`P'[`i',4])^(`k'*(`k'-1)/2) } *for Holm and Sidak correction, use _mtest command _mtest adjust `P', mtest(holm) pindex(4) append _mtest adjust r(result), mtest(sidak) pindex(4) append matrix `P' = r(result) } *Print results local len = cond(`adj',51,27) di di as res "Pairwise comparisons" di as txt "{hline 13}{c TT}{hline `len'}" di as txt _col(14) "{c |}" _col(29) "{ralign 13:Non corrected}" _c if (`adj') di as txt "{ralign 8:Bonferr}{ralign 8:Holm}{ralign 8:Sidak}" _c di as txt _n _col(14) "{c |} {ralign 9:Contrast} {ralign 6:t}{ralign 8:P>|t|}" _c if (`adj') di as txt "{ralign 8:P>|t|}{ralign 8:P>|t|}{ralign 8:P>|t|}" _c di as txt _n "{hline 13}{c +}{hline `len'}" local c = abbrev("`by'",12) di as txt "{ralign 12:`c'} {c |}" _c local nrows = rowsof(`P') local ncols = colsof(`P') forvalues i=1/`nrows' { local c = `lbl`i'' di as txt _newline "{ralign 12:`c'} {c |}" _c forvalues j=1/`ncols' { local fmt = cond(`j'<=2,"%9.0g",cond(`j'==3,"%6.0g","%6.4f")) if (`j'!=2) di as res " " `fmt' `P'[`i',`j'] _c } } di as txt _n "{hline 13}{c BT}{hline `len'}" *Stored results matrix `P' = (`P'[1..`nrows',1],`P'[1..`nrows',3..`ncols']) // delete StdErr column matrix rownames `P' = `rowlbls' if (!`adj') matrix colnames `P' = Contrast t NonCorr else matrix colnames `P' = Contrast t NonCorr Bonferr Holm Sidak return matrix results = `P' end