*! version 1.1.5 03may2018 by alexis dot dinno at pdx dot edu *! perform the Conover-Iman test of multiple comparisons using rank sums ******************************************************************************** * Syntax: conovertest varname [if exp] [in range] , by(groupvar) ma(method) * [ nokwallis nolabel wrap list rmc level(#) altp ] program define conovertest if ( int(_caller())<8 ) { di in r "conovertest- does not support this version of Stata." _newline di as txt "Requests for a version compatible with versions of STATA earlier than v8 are " di as txt "untenable since I do not have access to the software." _newline di as txt "All requests are welcome and will be considered." exit } else if ( int(_caller())<14 ) { conovertest8 `0' } else { conovertest14 `0' } end program define conovertest8, rclass sort version 8 syntax varname [if] [in], BY(varname) [ma(string) NOKWALLIS NOLABEL WRAP LIST /* */ RMC level(cilevel) ALTP] preserve * Validate nolabel and labels if needed if "`nolabel'" == "" { local labelname : value label `by' if "`labelname'" == "" { capture confirm numeric variable `by' if (!_rc) { di as res _newline "Warning: by() values are unlabeled, option nolabel implicit" _newline local nolabel = "nolabel" } else { encode `by', generate(_TempGroup) local groupvarname = "`by'" drop `by' rename _TempGroup `groupvarname' local labelname : value label `by' local nolabel = "" } } } * Validate ma if (lower("`ma'") == "" ) { local ma = "none" } if (lower("`ma'") != "none" & lower("`ma'") != "bonferroni" & lower("`ma'") != "sidak" & lower("`ma'") != "holm" & lower("`ma'") != "hs" & lower("`ma'") != "hochberg" & lower("`ma'") != "bh" & lower("`ma'") != "by") { noi: di as err "option ma() must be one of none, bonferroni, sidak, hs, hochberg, bh, or by" exit 198 } if (lower("`ma'")=="none") { local Name = "No adjustment" } if (lower("`ma'")=="bonferroni") { local Name = "Bonferroni" } if (lower("`ma'")=="sidak") { local Name = "Sid{c a'}k" } if (lower("`ma'")=="holm") { local Name = "Holm" } if (lower("`ma'")=="hs") { local Name = "Holm-Sid{c a'}k" } if (lower("`ma'")=="hochberg") { local Name = "Hochberg" } if (lower("`ma'")=="bh") { local Name = "Benjamini-Hochberg" } if (lower("`ma'")=="by") { local Name = "Benjamini-Yekuteili" } if (lower("`nokwallis'")=="") { kwallis `varlist' `if' `in', by(`by') local chi2 = r(chi2) local df = r(df) local chi2_adj = r(chi2_adj) } if (lower("`nokwallis'")!="") { qui: kwallis `varlist' `if' `in', by(`by') local chi2 = r(chi2) local df = r(df) local chi2_adj = r(chi2_adj) } qui: drop if `varlist' == . capture confirm numeric variable `by' if (!_rc) { qui: drop if `by' == . } else { qui: drop if `by' == "" encode `by', generate(_TempGroup) local groupvarname = "`by'" drop `by' rename _TempGroup `groupvarname' } if "`if'" != "" { qui: keep `if' } if "`in'" != "" { qui: keep `in' } generate __order = _n egen ranks = rank(`varlist') by `by', sort: egen __meanrank = sum(ranks) qui: tab __meanrank egen __group = group(`by') su __group, meanonly local k = r(max) qui: drop __group local m = `k'*(`k'-1)/2 if (c(matsize) < `m') { set matsize `m' } local N = r(N) egen __group = group(`by') su __group, meanonly forvalues i = 1/`r(max)' { qui: tabstat ranks if __group == `i', s(mean) save local meanrank`i' = el(r(StatTot),1,1) qui: count if __group == `i' local n`i' = r(N) } local var_mA = 1/(`N'-1) local var_mB = 0 forvalues i = 1(1)`N' { local var_mB = `var_mB' + (ranks[`i']^2 ) } local var_mB = `var_mB' - (`N'* ( ((`N'+1)^2)/4) ) local S2 = `var_mA'*`var_mB' * Pause to validate nolabel if "`nolabel'" == "" { qui: egen __blocks = group(`by') `if' `in' qui: tab __block `if' `in', nofreq local b = r(r) local labelname : value label `by' qui: levelsof `by' local groupvalues = r(levels) forvalues i = 1(1)`k' { * get the ith value label local grouplabel = "`:label `labelname' `:word `i' of `groupvalues'''" if "`grouplabel'" == "" | real("`grouplabel'") !=. { di as res _newline "Warning: by() values are incompletely labeled, option nolabel implicit" local nolabel = "nolabel" break } } } qui: sort __order di _newline as txt %~76s "Conover-Iman Pairwise Comparison of `varlist' by `by'" di as txt %~76s "(`Name')" local colhead = "" local kminusone = `k'-1 capture confirm numeric variable `by' if (!_rc) { local groupisnum = 1 } else { local groupisnum = 0 } qui: levelsof `by' local groupvalues = r(levels) matrix T = J(1,`m',0) local nu = `N' - `k' matrix P = J(1,`m',0) matrix Reject = J(1,`m',0) forvalues i = 2/`k' { local iminusone = `i'-1 forvalues j = 1/`iminusone' { if "`rmc'" == "" { local t = (`meanrank`j''-`meanrank`i'')/(sqrt(`S2'*((`N'-1-`chi2_adj')/(`N'-`k')))*sqrt(1/`n`j'' + 1/`n`i'')) } if "`rmc'" != "" { local t = (`meanrank`i''-`meanrank`j'')/(sqrt(`S2'*((`N'-1-`chi2_adj')/(`N'-`k')))*sqrt(1/`n`j'' + 1/`n`i'')) } * If p = P(T <= |t|) if ("`altp'"=="") { local p = ttail(`nu',abs(`t')) } * Otherwise, if p = P(|T| <= |t|) else { local p = 2*ttail(`nu',abs(`t')) } * Static multiple comparisons adjustments if (lower("`ma'") == "bonferroni") { local p = min(1,`p'*`m') } if (lower("`ma'") == "sidak") { local p = min(1,1 - (1-`p')^`m') } local index = ((`i'-2)*(`i'-1)/2) + `j' matrix T[1,`index'] = `t' matrix P[1,`index'] = `p' local alpha = (100 - `level')/100 * If p = P(T <= |t|) if ("`altp'"=="") { if el(matrix(P),1,`index') <= `alpha'/2 { matrix Reject[1,`index'] = 1 } } * Otherwise, if p = P(|T| <= |t|) else { if el(matrix(P),1,`index') <= `alpha' { matrix Reject[1,`index'] = 1 } } } } * Sequential multiple comparrisons adjustments noi: di as inp "Here?" quietly { local alpha = (100 - `level')/100 if (lower("`ma'")=="holm") { matrix Reject = J(1,`m',0) mata: Ps = st_matrix("P")\range(1,`m',1)'\rangen(0,0,`m')' mata: Psort = sort(Ps',1)' forvalues i = 1/`m' { local adjust = `m'+1-`i' mata: Psort[1,`i'] = min((1,Psort[1,`i'] :*`adjust')) * If p = P(T <= |t|) if ("`altp'"=="") { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha'/2 } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha'/2) & Psort[3,`i'-1] != 0 } } * Otherwise, if p = P(|T| <= |t|) else { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha' } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha') & Psort[3,`i'-1] != 0 } } } mata: Psort = sort(Psort',2)' mata: st_matrix("P",Psort[1,]) mata: st_matrix("Reject",Psort[3,]) } if (lower("`ma'")=="hs") { matrix Reject = J(1,`m',0) mata: Ps = st_matrix("P")\range(1,`m',1)'\rangen(0,0,`m')' mata: Psort = sort(Ps',1)' forvalues i = 1/`m' { local adjust = `m'+1-`i' mata: Psort[1,`i'] = min((1,(1 - ((1 - Psort[1,`i']) ^`adjust')))) * If p = P(T <= |t|) if ("`altp'"=="") { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha'/2 } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha'/2) & Psort[3,`i'-1] != 0 } } * Otherwise, if p = P(|T| <= |t|) else { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha' } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha') & Psort[3,`i'-1] != 0 } } } mata: Psort = sort(Psort',2)' mata: st_matrix("P",Psort[1,]) mata: st_matrix("Reject",Psort[3,]) } if (lower("`ma'")=="hochberg") { matrix Reject = J(1,`m',0) mata: Ps = st_matrix("P")\range(1,`m',1)'\rangen(0,0,`m')' mata: Psort = sort(Ps',-1)' forvalues i = 1/`m' { local adjust = `i' mata: Psort[1,`i'] = min((1,Psort[1,`i'] :*`adjust')) * If p = P(T <= |t|) if ("`altp'"=="") { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha'/2 } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha'/2) | Psort[3,`i'-1] == 1 } } * Otherwise, if p = P(|T| <= |t|) else { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha' } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha') | Psort[3,`i'-1] == 1 } } } mata: Psort = sort(Psort',2)' mata: st_matrix("P",Psort[1,]) mata: st_matrix("Reject",Psort[3,]) } if (lower("`ma'")=="bh") { matrix Reject = J(1,`m',0) mata: Ps = st_matrix("P")\range(1,`m',1)'\rangen(0,0,`m')' mata: Psort = sort(Ps',-1)' forvalues i = 1/`m' { local adjust = (`m'/(`m'+1-`i')) mata: Psort[1,`i'] = min((1,Psort[1,`i'] :*`adjust')) * If p = P(T <= |t|) if ("`altp'"=="") { if (`i' == 1) { noi: mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha'/2 } if (`i' > 1) { noi: mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha'/2) | Psort[3,(`i'-1)] == 1 } } * Otherwise, if p = P(|T| <= |t|) else { if (`i' == 1) { noi: mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha' } if (`i' > 1) { noi: mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha') | Psort[3,(`i'-1)] == 1 } } } mata: Psort = sort(Psort',2)' mata: st_matrix("P",Psort[1,]) mata: st_matrix("Reject",Psort[3,]) } if (lower("`ma'")=="by") { matrix Reject = J(1,`m',0) mata: Ps = st_matrix("P")\range(1,`m',1)'\rangen(0,0,`m')' mata: Psort = sort(Ps',-1)' local adjustB = 0 forvalues i=1/`m' { local C = `C' + 1/`i' } forvalues i = 1/`m' { local adjust = (`m'/(`m'+1-`i')) * `C' mata: Psort[1,`i'] = min((1,Psort[1,`i'] :* `adjust' )) * If p = P(T <= |t|) if ("`altp'"=="") { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha'/2 } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha'/2) | Psort[3,`i'-1] == 1 } } * Otherwise, if p = P(|T| <= |t|) else { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha' } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha') | Psort[3,`i'-1] == 1 } } } mata: Psort = sort(Psort',2)' mata: st_matrix("P",Psort[1,]) mata: st_matrix("Reject",Psort[3,]) } } ******************************************************************************** *** OUTPUT ******************************************************************************** * Need to determine how many tables (reps) to output local reps = int((`k'-1)/6) local laststart = `reps'*6 + 1 local kminusone = `k' - 1 * Replication loop for >7 groups, no wrap if ("`wrap'"=="") { if (`k' > 7) { forvalues rep = 1/`reps' { local colstart = (6*`rep')-5 local colstop = 6*`rep' _conovertestheader `by' `colstart' `colstop' "`nolabel'" `groupisnum' "`rmc'" * Table body local start = `colstart'+1 forvalues i = `start'/`k' { local colstop = min(`i'-1,6*`rep') _conovertestttable `by' `i' T `colstart' `colstop' "`nolabel'" `groupisnum' if (`i' < `k') { _conovertestptable `i' P `colstart' `colstop' Reject 0 } else { _conovertestptable `i' P `colstart' `colstop' Reject 1 } } } * End of table if (`laststart' < `k') { _conovertestheader `by' `laststart' `kminusone' "`nolabel'" `groupisnum' "`rmc'" * Table body local start = `laststart'+1 forvalues i = `start'/`k' { local colstop = `i' - 1 _conovertestttable `by' `i' T `laststart' `colstop' "`nolabel'" `groupisnum' if (`i' < `k') { _conovertestptable `i' P `laststart' `colstop' Reject 0 } else { _conovertestptable `i' P `laststart' `colstop' Reject 1 } } } } * Replication loop for <=7 groups if (`k' <= 7) { _conovertestheader `by' 1 `kminusone' "`nolabel'" `groupisnum' "`rmc'" * Table body forvalues i = 2/`k' { local colstop = `i'-1 _conovertestttable `by' `i' T 1 `colstop' "`nolabel'" `groupisnum' if (`i' < `k') { _conovertestptable `i' P 1 `colstop' Reject 0 } else { _conovertestptable `i' P 1 `colstop' Reject 1 } } } } * Replication loop for >7 groups, with wrap if ("`wrap'"=="wrap") { _conovertestheader `by' 1 `kminusone' "`nolabel'" `groupisnum' "`rmc'" * Table body forvalues i = 2/`k' { local colstop = `i'-1 _conovertestttable `by' `i' T 1 `colstop' "`nolabel'" `groupisnum' if (`i' < `k') { _conovertestptable `i' P 1 `colstop' Reject 0 } else { _conovertestptable `i' P 1 `colstop' Reject 1 } } } * Output pairwise comparisons as a list if requested. if ("`list'"=="list") { _conovertestlist `by' T P Reject `k' "`ma'" "`rmc'" "`nolabel'" } * Output level of significance local alpha = (100 - `level')/100 if length("`level'") > 2 { local precision = length("`level'") - 1 } else { local precision = 2 } * If p = P(T <= |t|) if ("`altp'"=="") { if ("`ma'" == "none") { di _newline as text "alpha = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(T <= |t|) <= alpha/2" } else if ("`ma'" == "bonferroni" | "`ma'" == "sidak") { di _newline as text "Family-wise Error Rate = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(T <= |t|) <= FWER/2" } else if ("`ma'" == "holm" | "`ma'" == "hs") { di _newline as text "Family-wise Error Rate = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(T <= |t|) <= FWER/2 with stopping rule" } else if ("`ma'" == "hochberg" | "`ma'" == "bh" | "`ma'" == "by") { di _newline as text "False Discovery Rate = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(T <= |t|) <= FDR/2 with stopping rule" } } * Otherwise, if p = P(|T| <= |t|) else { if ("`ma'" == "none") { di _newline as text "alpha = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(|T| <= |t|) <= alpha" } else if ("`ma'" == "bonferroni" | "`ma'" == "sidak") { di _newline as text "Family-wise Error Rate = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(|T| <= |t|) <= FWER" } else if ("`ma'" == "holm" | "`ma'" == "hs") { di _newline as text "Family-wise Error Rate = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(|T| <= |t|) <= FWER with stopping rule" } else if ("`ma'" == "hochberg" | "`ma'" == "bh" | "`ma'" == "by") { di _newline as text "False Discovery Rate = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(|T| <= |t|) <= FDR with stopping rule" } } restore * Saves * If p = P(T <= |t|) if ("`altp'"=="") { ret matrix P = P } * Otherwise, if p = P(|T| <= |t|) else { ret matrix altP = P } ret matrix T = T ret scalar df = `df' ret scalar chi2_adj = `chi2_adj' end program define conovertest14, rclass sort version 14 syntax varname [if] [in], BY(varname) [ma(string) NOKWALLIS NOLABEL WRAP LIST /* */ RMC level(cilevel) ALTP ] preserve * Validate nolabel and labels if needed if "`nolabel'" == "" { local labelname : value label `by' if "`labelname'" == "" { capture confirm numeric variable `by' if (!_rc) { di as res _newline "Warning: by() values are unlabeled, option nolabel implicit" _newline local nolabel = "nolabel" } else { encode `by', generate(_TempGroup) local groupvarname = "`by'" drop `by' rename _TempGroup `groupvarname' local labelname : value label `by' local nolabel = "" } } } * Validate ma if (lower("`ma'") == "" ) { local ma = "none" } if (lower("`ma'") != "none" & lower("`ma'") != "bonferroni" & lower("`ma'") != "sidak" & lower("`ma'") != "holm" & lower("`ma'") != "hs" & lower("`ma'") != "hochberg" & lower("`ma'") != "bh" & lower("`ma'") != "by") { noi: di as err "option ma() must be one of none, bonferroni, sidak, hs, hochberg, bh, or by" exit 198 } if (lower("`ma'")=="none") { local Name = "No adjustment" } if (lower("`ma'")=="bonferroni") { local Name = "Bonferroni" } if (lower("`ma'")=="sidak") { local Name = "Sid{c a'}k" } if (lower("`ma'")=="holm") { local Name = "Holm" } if (lower("`ma'")=="hs") { local Name = "Holm-Sid{c a'}k" } if (lower("`ma'")=="hochberg") { local Name = "Hochberg" } if (lower("`ma'")=="bh") { local Name = "Benjamini-Hochberg" } if (lower("`ma'")=="by") { local Name = "Benjamini-Yekuteili" } if (lower("`nokwallis'")=="") { kwallis `varlist' `if' `in', by(`by') local chi2 = r(chi2) local df = r(df) local chi2_adj = r(chi2_adj) } if (lower("`nokwallis'")!="") { qui: kwallis `varlist' `if' `in', by(`by') local chi2 = r(chi2) local df = r(df) local chi2_adj = r(chi2_adj) } qui: drop if `varlist' == . capture confirm numeric variable `by' if (!_rc) { qui: drop if `by' == . } else { qui: drop if `by' == "" encode `by', generate(_TempGroup) local groupvarname = "`by'" drop `by' rename _TempGroup `groupvarname' } if "`if'" != "" { qui: keep `if' } if "`in'" != "" { qui: keep `in' } generate __order = _n egen ranks = rank(`varlist') by `by', sort: egen __meanrank = sum(ranks) qui: tab __meanrank egen __group = group(`by') su __group, meanonly local k = r(max) qui: drop __group local m = `k'*(`k'-1)/2 if (c(matsize) < `m') { set matsize `m' } local N = r(N) egen __group = group(`by') su __group, meanonly forvalues i = 1/`r(max)' { qui: tabstat ranks if __group == `i', s(mean) save local meanrank`i' = el(r(StatTotal),1,1) qui: count if __group == `i' local n`i' = r(N) } local var_mA = 1/(`N'-1) local var_mB = 0 forvalues i = 1(1)`N' { local var_mB = `var_mB' + (ranks[`i']^2 ) } local var_mB = `var_mB' - (`N'* ( ((`N'+1)^2)/4) ) local S2 = `var_mA'*`var_mB' * Pause to validate nolabel if "`nolabel'" == "" { qui: egen __blocks = group(`by') `if' `in' qui: tab __block `if' `in', nofreq local b = r(r) local labelname : value label `by' qui: levelsof `by' local groupvalues = r(levels) forvalues i = 1(1)`k' { * get the ith value label local grouplabel = "`:label `labelname' `:word `i' of `groupvalues'''" if "`grouplabel'" == "" | real("`grouplabel'") !=. { di as res _newline "Warning: by() values are incompletely labeled, option nolabel implicit" local nolabel = "nolabel" break } } } qui: sort __order di _newline as txt %~76s "Conover-Iman Pairwise Comparison of `varlist' by `by'" di as txt %~76s "(`Name')" local colhead = "" local kminusone = `k'-1 capture confirm numeric variable `by' if (!_rc) { local groupisnum = 1 } else { local groupisnum = 0 } qui: levelsof `by' local groupvalues = r(levels) matrix T = J(1,`m',0) local nu = `N' - `k' matrix P = J(1,`m',0) matrix Reject = J(1,`m',0) forvalues i = 2/`k' { local iminusone = `i'-1 forvalues j = 1/`iminusone' { if "`rmc'" == "" { local t = (`meanrank`j''-`meanrank`i'')/(sqrt(`S2'*((`N'-1-`chi2_adj')/(`N'-`k')))*sqrt(1/`n`j'' + 1/`n`i'')) } if "`rmc'" != "" { local t = (`meanrank`i''-`meanrank`j'')/(sqrt(`S2'*((`N'-1-`chi2_adj')/(`N'-`k')))*sqrt(1/`n`j'' + 1/`n`i'')) } * If p = P(T <= |t|) if ("`altp'"=="") { local p = ttail(`nu',abs(`t')) } * Otherwise, if p = P(|T| <= |t|) else { local p = 2*ttail(`nu',abs(`t')) } * Static multiple comparisons adjustments if (lower("`ma'") == "bonferroni") { local p = min(1,`p'*`m') } if (lower("`ma'") == "sidak") { local p = min(1,1 - (1-`p')^`m') } local index = ((`i'-2)*(`i'-1)/2) + `j' matrix T[1,`index'] = `t' matrix P[1,`index'] = `p' local alpha = (100 - `level')/100 * If p = P(T <= |t|) if ("`altp'"=="") { if el(matrix(P),1,`index') <= `alpha'/2 { matrix Reject[1,`index'] = 1 } } * Otherwise, if p = P(|T| <= |t|) else { if el(matrix(P),1,`index') <= `alpha' { matrix Reject[1,`index'] = 1 } } } } * Sequential multiple comparrisons adjustments quietly { local alpha = (100 - `level')/100 if (lower("`ma'")=="holm") { matrix Reject = J(1,`m',0) mata: Ps = st_matrix("P")\range(1,`m',1)'\rangen(0,0,`m')' mata: Psort = sort(Ps',1)' forvalues i = 1/`m' { local adjust = `m'+1-`i' mata: Psort[1,`i'] = min((1,Psort[1,`i'] :*`adjust')) * If p = P(T <= |t|) if ("`altp'"=="") { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha'/2 } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha'/2) & Psort[3,`i'-1] != 0 } } * Otherwise, if p = P(|T| <= |t|) else { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha' } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha') & Psort[3,`i'-1] != 0 } } } mata: Psort = sort(Psort',2)' mata: st_matrix("P",Psort[1,]) mata: st_matrix("Reject",Psort[3,]) } if (lower("`ma'")=="hs") { matrix Reject = J(1,`m',0) mata: Ps = st_matrix("P")\range(1,`m',1)'\rangen(0,0,`m')' mata: Psort = sort(Ps',1)' forvalues i = 1/`m' { local adjust = `m'+1-`i' mata: Psort[1,`i'] = min((1,(1 - ((1 - Psort[1,`i']) ^`adjust')))) * If p = P(T <= |t|) if ("`altp'"=="") { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha'/2 } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha'/2) & Psort[3,`i'-1] != 0 } } * Otherwise, if p = P(|T| <= |t|) else { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha' } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha') & Psort[3,`i'-1] != 0 } } } mata: Psort = sort(Psort',2)' mata: st_matrix("P",Psort[1,]) mata: st_matrix("Reject",Psort[3,]) } if (lower("`ma'")=="hochberg") { matrix Reject = J(1,`m',0) mata: Ps = st_matrix("P")\range(1,`m',1)'\rangen(0,0,`m')' mata: Psort = sort(Ps',-1)' forvalues i = 1/`m' { local adjust = `i' mata: Psort[1,`i'] = min((1,Psort[1,`i'] :*`adjust')) * If p = P(T <= |t|) if ("`altp'"=="") { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha'/2 } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha'/2) | Psort[3,`i'-1] == 1 } } * Otherwise, if p = P(|T| <= |t|) else { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha' } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha') | Psort[3,`i'-1] == 1 } } } mata: Psort = sort(Psort',2)' mata: st_matrix("P",Psort[1,]) mata: st_matrix("Reject",Psort[3,]) } if (lower("`ma'")=="bh") { matrix Reject = J(1,`m',0) mata: Ps = st_matrix("P")\range(1,`m',1)'\rangen(0,0,`m')' mata: Psort = sort(Ps',-1)' forvalues i = 1/`m' { local adjust = (`m'/(`m'+1-`i')) mata: Psort[1,`i'] = min((1,Psort[1,`i'] :*`adjust')) * If p = P(T <= |t|) if ("`altp'"=="") { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha'/2 } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha'/2) | Psort[3,`i'-1] == 1 } } * Otherwise, if p = P(|T| <= |t|) else { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha' } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha') | Psort[3,`i'-1] == 1 } } } mata: Psort = sort(Psort',2)' mata: st_matrix("P",Psort[1,]) mata: st_matrix("Reject",Psort[3,]) } if (lower("`ma'")=="by") { matrix Reject = J(1,`m',0) mata: Ps = st_matrix("P")\range(1,`m',1)'\rangen(0,0,`m')' mata: Psort = sort(Ps',-1)' local adjustB = 0 forvalues i=1/`m' { local C = `C' + 1/`i' } forvalues i = 1/`m' { local adjust = (`m'/(`m'+1-`i')) * `C' mata: Psort[1,`i'] = min((1,Psort[1,`i'] :* `adjust' )) * If p = P(T <= |t|) if ("`altp'"=="") { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha'/2 } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha'/2) | Psort[3,`i'-1] == 1 } } * Otherwise, if p = P(|T| <= |t|) else { if (`i' == 1) { mata: Psort[3,`i'] = Psort[1,`i'] <= `alpha' } if (`i' > 1) { mata: Psort[3,`i'] = (Psort[1,`i'] <= `alpha') | Psort[3,`i'-1] == 1 } } } mata: Psort = sort(Psort',2)' mata: st_matrix("P",Psort[1,]) mata: st_matrix("Reject",Psort[3,]) } } ******************************************************************************** *** OUTPUT ******************************************************************************** * Need to determine how many tables (reps) to output local reps = int((`k'-1)/6) local laststart = `reps'*6 + 1 local kminusone = `k' - 1 * Replication loop for >7 groups, no wrap if ("`wrap'"=="") { if (`k' > 7) { forvalues rep = 1/`reps' { local colstart = (6*`rep')-5 local colstop = 6*`rep' _conovertestheader `by' `colstart' `colstop' "`nolabel'" `groupisnum' "`rmc'" * Table body local start = `colstart'+1 forvalues i = `start'/`k' { local colstop = min(`i'-1,6*`rep') _conovertestttable `by' `i' T `colstart' `colstop' "`nolabel'" `groupisnum' if (`i' < `k') { _conovertestptable `i' P `colstart' `colstop' Reject 0 } else { _conovertestptable `i' P `colstart' `colstop' Reject 1 } } } * End of table if (`laststart' < `k') { _conovertestheader `by' `laststart' `kminusone' "`nolabel'" `groupisnum' "`rmc'" * Table body local start = `laststart'+1 forvalues i = `start'/`k' { local colstop = `i' - 1 _conovertestttable `by' `i' T `laststart' `colstop' "`nolabel'" `groupisnum' if (`i' < `k') { _conovertestptable `i' P `laststart' `colstop' Reject 0 } else { _conovertestptable `i' P `laststart' `colstop' Reject 1 } } } } * Replication loop for <=7 groups if (`k' <= 7) { _conovertestheader `by' 1 `kminusone' "`nolabel'" `groupisnum' "`rmc'" * Table body forvalues i = 2/`k' { local colstop = `i'-1 _conovertestttable `by' `i' T 1 `colstop' "`nolabel'" `groupisnum' if (`i' < `k') { _conovertestptable `i' P 1 `colstop' Reject 0 } else { _conovertestptable `i' P 1 `colstop' Reject 1 } } } } * Replication loop for >7 groups, with wrap if ("`wrap'"=="wrap") { _conovertestheader `by' 1 `kminusone' "`nolabel'" `groupisnum' "`rmc'" * Table body forvalues i = 2/`k' { local colstop = `i'-1 _conovertestttable `by' `i' T 1 `colstop' "`nolabel'" `groupisnum' if (`i' < `k') { _conovertestptable `i' P 1 `colstop' Reject 0 } else { _conovertestptable `i' P 1 `colstop' Reject 1 } } } * Output pairwise comparisons as a list if requested. if ("`list'"=="list") { _conovertestlist `by' T P Reject `k' "`ma'" "`rmc'" "`nolabel'" } * Output level of significance local alpha = (100 - `level')/100 if ustrlen("`level'") > 2 { local precision = ustrlen("`level'") - 1 } else { local precision = 2 } * If p = P(T <= |t|) if ("`altp'"=="") { if ("`ma'" == "none") { di _newline as text "alpha = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(T <= |t|) <= alpha/2" } else if ("`ma'" == "bonferroni" | "`ma'" == "sidak") { di _newline as text "Family-wise Error Rate = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(T <= |t|) <= FWER/2" } else if ("`ma'" == "holm" | "`ma'" == "hs") { di _newline as text "Family-wise Error Rate = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(T <= |t|) <= FWER/2 with stopping rule" } else if ("`ma'" == "hochberg" | "`ma'" == "bh" | "`ma'" == "by") { di _newline as text "False Discovery Rate = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(T <= |t|) <= FDR/2 with stopping rule" } } * Otherwise, if p = P(|T| <= |t|) else { if ("`ma'" == "none") { di _newline as text "alpha = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(|T| <= |t|) <= alpha" } else if ("`ma'" == "bonferroni" | "`ma'" == "sidak") { di _newline as text "Family-wise Error Rate = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(|T| <= |t|) <= FWER" } else if ("`ma'" == "holm" | "`ma'" == "hs") { di _newline as text "Family-wise Error Rate = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(|T| <= |t|) <= FWER with stopping rule" } else if ("`ma'" == "hochberg" | "`ma'" == "bh" | "`ma'" == "by") { di _newline as text "False Discovery Rate = " as res %6.`precision'f `alpha' di as txt "Reject Ho if p = P(|T| <= |t|) <= FDR with stopping rule" } } restore * Saves * If p = P(T <= |t|) if ("`altp'"=="") { ret matrix P = P } * Otherwise, if p = P(|T| <= |t|) else { ret matrix altp = P } ret matrix T = T ret scalar df = `df' ret scalar chi2_adj = `chi2_adj' end program define _conovertestheader /* This program displays Conover-Iman test table headers. Syntax: _conovertestheader groupvar colstart colstop nolabel groupisnum rmc */ args groupvar colstart colstop nolabel groupisnum rmc if "`rmc'" == "" { di as txt "Col Mean-{c |}" di as txt "Row Mean {c |}" _continue } else { di as txt "Row Mean-{c |}" di as txt "Col Mean {c |}" _continue } qui: levelsof `groupvar' local groupvalues = r(levels) forvalues col = `colstart'/`colstop' { if (lower("`nolabel'")=="") { local labelname : value label `groupvar' local vallab = substr("`:label `labelname' `:word `col' of `groupvalues'''",1,8) local pad = 8-length("`vallab'") local colhead = "{dup `pad': }`vallab'" if ("`vallab'" == "") { local colhead = "" local nolabel = "nolabel" } } if (lower("`nolabel'")=="nolabel") { if (`groupisnum' == 0) { local pad = max(0,8-length(substr(`groupvar'[`col'],1,8))) local colhead = " {dup `pad': }"+ substr(`groupvar'[`col'],1,8) } if (`groupisnum' == 1) { local pad = max(0,8-length(substr(word("`groupvalues'",`col'),1,8))) local colhead = " {dup `pad': }"+ substr(word("`groupvalues'",`col'),1,8) } } di as txt " " %8s "`colhead'" _continue } di local separatorlength = 10 + 11*(`colstop'-`colstart')+1 di as txt "{hline 9}{c +}{hline `separatorlength'}" end program define _conovertestttable /* This program displays Conover-Iman test t values. Syntax: _conovertestttable groupvar index T colstart colstop nolabel groupisnum */ args groupvar index T colstart colstop nolabel groupisnum qui: levelsof `groupvar' local groupvalues = r(levels) * Validate labels if (lower("`nolabel'")=="") { local testlab: label (`groupvar') `colstart' 8, strict if ("`testlab'" == "") { local nolabel = "nolabel" } } * Row headers if (lower("`nolabel'")=="") { local labelname : value label `groupvar' local vallab = substr("`:label `labelname' `:word `index' of `groupvalues'''",1,8) di as txt %8s "`vallab'" " {c |}" _continue } if (lower("`nolabel'")=="nolabel") { if (`groupisnum' == 0) { di as txt %8s substr(`groupvar'[`index'],1,8) " {c |}" _continue } if (`groupisnum' == 1) { di as txt %8s substr(word("`groupvalues'",`index'),1,8) " {c |}" _continue } } * Table t entries forvalues i = `colstart'/`colstop' { di as res " " %9.6f el(matrix(`T'),1,((`index'-2)*(`index'-1))/2 + `i') _continue } di end program define _conovertestptable /* This program displays Conover-Iman test p values. Syntax: _conovertestptable index P colstart colstop Reject last */ args index P colstart colstop Reject last * Blank row header di as txt " {c |}" _continue * Table p entries forvalues i = `colstart'/`colstop' { if ( el(matrix(`Reject'),1,((`index'-2)*(`index'-1))/2 + `i') == 0) { di as res " " %6.4f el(matrix(`P'),1,((`index'-2)*(`index'-1))/2 + `i') _continue } else { di as res " {ul on}" %6.4f el(matrix(`P'),1,((`index'-2)*(`index'-1))/2 + `i') "{ul off}" _continue } } di * Close out with another blank row header if (`last' == 0) { di as txt " {c |}" } end program define _conovertestlist /* This program displays Conover-Iman test t and p values as a list. Syntax: _conovertestlist groupvar T P Reject k ma rmc nolabel */ args groupvar T P Reject k ma rmc nolabel qui: levelsof `groupvar' local groupvalues = r(levels) * Output list header, depending on whether the ma option is used if (lower("`ma'")=="none" || "`ma'"=="") { noi: di as text _newline "List of pairwise comparisons: t statistic (p-value)" } else { noi: di as text _newline "List of pairwise comparisons: t statistic (adjusted p-value)" } * flag whether any rejections occur local maxreject = 0 forvalues i = 1(1)`k' { if el(matrix(`Reject'),1,`i') == 1 { local maxreject = 1 } } * Validate labels if (lower("`nolabel'")=="") { forvalues i = 1(1)`k' { local testlab: label (`groupvar') `i', strict if ("`testlab'" == "") { local nolabel = "nolabel" } } } * List headers if (lower("`nolabel'")=="") { * get the labelname for the group labels local labelname : value label `groupvar' if ("`labelname'") != "" { * get the length of the largest group label (whether explicitly labeled or not) local firstlongest = 1 forvalues i = 1(1)`k' { * get the length of the ith value label local currentlength = ustrlen("`:label `labelname' `:word `i' of `groupvalues'''") if `currentlength' > `firstlongest' { local firstlongest = `currentlength' } } * get the secondlongest group label length (whether explicitly labeled or not) local secondlongest = 1 forvalues i = 1(1)`k' { * get the length of the ith value label local currentlength = ustrlen("`:label `labelname' `:word `i' of `groupvalues'''") if `currentlength' < `firstlongest' & `currentlength' > `secondlongest' { local secondlongest = `currentlength' } } * replace secondlongest if the value of firstlongest appears twice local twice = 0 forvalues i = 1(1)`k' { local currentlength = ustrlen("`:label `labelname' `:word `i' of `groupvalues'''") if `currentlength' == `firstlongest' { local twice = `twice' + 1 } if `twice' == 2 { local secondlongest = `firstlongest' } } * stringlength will be the sum of the two largest values in lengths plus 6 local stringlength = `firstlongest' + `secondlongest' + 6 * output underline length, dependent upon group labels local underline = `stringlength' - 2 noi: di as text "{hline `underline'}{c TT}{hline 19}{hline `maxreject'}" * Start listing output! local index = 0 forvalues i = 2(1)`k' { local jmax = `i'-1 forvalues j = 1(1)`jmax' { local index = `index' + 1 local labeli = "`:label `labelname' `:word `i' of `groupvalues'''" local labelj = "`:label `labelname' `:word `j' of `groupvalues'''" local bufferlength = max(`stringlength' - (ustrlen("`labeli'") + ustrlen("`labelj'") + 4) - 2,0) if "`rmc'" == "" { noi: di as text "`labelj' - `labeli'{dup `bufferlength': } {c |} " as res %9.6f el(matrix(`T'),1,`index') as text " (" as res %6.4f el(matrix(`P'),1,`index') as text ")" _continue } if "`rmc'" == "rmc" { noi: di as text "`labeli' - `labelj'{dup `bufferlength': } {c |} " as res %9.6f el(matrix(`T'),1,`index') as text " (" as res %6.4f el(matrix(`P'),1,`index') as text ")" _continue } if el(matrix(`Reject'),1,`index') == 1 { noi: di as res "*" _continue } noi: di } } } else { local nolabel = "nolabel" } } if (lower("`nolabel'")=="nolabel") { * get the lengths of the two largest group numbers local firstlongest = 1 local secondlongest = 1 forvalues i = 1(1)`k' { * get the length of the ith value label local currentlength = floor(log10(real(word("`groupvalues'",`i')))) + 1 if `currentlength' > `firstlongest' { local firstlongest = `currentlength' } if `currentlength' < `firstlongest' & `currentlength' > `secondlongest' { local secondlongest = `currentlength' } } * replace secondlongest if the value of firstlongest appears twice local twice = 0 forvalues i = 1(1)`k' { local currentlength = floor(log10(real(word("`groupvalues'",`i')))) + 1 if `currentlength' == `firstlongest' { local twice = `twice' + 1 } if `twice' == 2 { local secondlongest = `firstlongest' } } * stringlength will be the sum of the two largest values in lengths plus 6 local stringlength = `firstlongest' + `secondlongest' + 6 * output underline length, dependent upon group labels local underline = `stringlength' - 2 noi: di as text "{hline `underline'}{c TT}{hline 19}{hline `maxreject'}" * Start listing output! local index = 0 forvalues i = 2(1)`k' { local jmax = `i'-1 forvalues j = 1(1)`jmax' { local index = `index' + 1 local labeli = word("`groupvalues'",`i') local labelj = word("`groupvalues'",`j') local bufferlength = max(`stringlength' - (floor(log10(`labeli')) + floor(log10(`labelj')) + 6) - 2,0) if "`rmc'" == "" { noi: di as text "`labelj' - `labeli'{dup `bufferlength': } {c |} " as res %9.6f el(matrix(`T'),1,`index') as text " (" as res %6.4f el(matrix(`P'),1,`index') as text ")" _continue } if "`rmc'" == "rmc" { noi: di as text "`labeli' - `labelj'{dup `bufferlength': } {c |} " as res %9.6f el(matrix(`T'),1,`index') as text " (" as res %6.4f el(matrix(`P'),1,`index') as text ")" _continue } if el(matrix(`Reject'),1,`index') == 1 { noi: di as res "*" _continue } noi: di } } } end