*! version 1.0.0 Ariel Linden 02sep2014
*! code partially based on signrank version 2.2.8  28sep2004

program define alignedpairs, rclass byable(recall)
	version 13
	syntax varname [=/exp] [if] [in] , [ Level(real 95) ]

	preserve
	tempname tp tn v unv z adj0 effect nN Calpha lb_1 ub_1 alpha
	tempvar touse diff ranks t zavg

quietly {
	mark `touse' `if' `in'
	keep if `touse' 	//keeping the touse data eliminates problems with the permutations

		gen double `diff' = `varlist'-(`exp') // if `touse'
		egen double `ranks' = rank(abs(`diff')) //if `touse'

/* We do want to OMIT the ranks corresponding to `diff'==0 in the sums.  */

		gen double `t' = sum(cond(`diff'>0,`ranks',0))
		scalar `tp' = `t'[_N]

		replace `t' = sum(cond(`diff'<0,`ranks',0))
		scalar `tn' = `t'[_N]

		replace `t' = sum(cond(`diff'~=0,`ranks'*`ranks',0))
		scalar `v' = `t'[_N]/4
		scalar `z' = (`tp'-`tn')/(2*sqrt(`v'))

		count //if `touse'
		local n = r(N)
		scalar `unv' = `n'*(`n'+1)*(2*`n'+1)/24

		count if `diff' == 0 //& `touse'
		local n0 = r(N)
		scalar `adj0' = -`n0'*(`n0'+1)*(2*`n0'+1)/24

		count if `diff' > 0 //& `touse'
		local np = r(N)
		local nn = `n' - `np' - `n0'

		* HL permutations of x = y
		local N = `n' * (`n' + 1) / 2
		set obs `N'
		gen `zavg' = .
		local k = 1
		local J = 1
		forval i = 1/`n' {
		forval j = `J'/`n' {
		replace `zavg' = (`diff'[`i'] + `diff'[`j'])/2 in `k++'
		}
		local ++J
		}

		sort `zavg'
		centile `zavg'
		scalar `effect' = r(c_1)
		scalar `nN' = r(N)	//get n(n+1)/2 for ub_1

		scalar `alpha'=(100-`level')/100
		scalar `Calpha' = round(`n'*(`n'+1)/4 - invnorm(1-`alpha'/2)*(`n'*(`n'+1)*(2*`n'+1)/24)^.5,1)

	* set Calfa to 1 for cases when it is 0 (small samples)
		if `Calpha' <1 {
		scalar `Calpha' = 1
		}
		else {
		scalar `Calpha' = `Calpha'
		}

		scalar `lb_1' = round(`zavg'[`Calpha'],.001)
		scalar `ub_1' = round(`zavg'[`nN'-`Calpha'+1],.001)
}

	di _n in gr `"Wilcoxon signed-rank test"' _n
	di in smcl in gr `"        sign {c |}      obs   sum ranks    expected"'
	di in smcl in gr "{hline 13}{c +}{hline 33}"
	ditablin positive `np' `tp' (`tp'+`tn')/2
	ditablin negative `nn' `tn' (`tp'+`tn')/2
	ditablin zero     `n0' `n0'*(`n0'+1)/2 `n0'*(`n0'+1)/2
	di in smcl in gr "{hline 13}{c +}{hline 33}"
	ditablin all `n' `n'*(`n'+1)/2 `n'*(`n'+1)/2

	if `unv' < 1e7  local vfmt `"%10.2f"'
	else            local vfmt `"%10.0g"'

	di in smcl in gr _n `"unadjusted variance"' _col(22) ///
	   in ye `vfmt' `unv' _n							///
	   in gr `"adjustment for ties"' _col(22)			///
	   in ye `vfmt' `v'-`unv'-`adj0' _n					///
	   in gr `"adjustment for zeros"' _col(22)			///
	   in ye `vfmt' `adj0' _n							///
	   in gr _col(22) "{hline 10}" _n					///
	   in gr `"adjusted variance"' _col(22)				///
	   in ye `vfmt' `v' _n(2)							///
	   in gr `"Ho: `varlist' = `exp'"' _n				///
	   in gr _col(14) `"z = "'							///
	   in ye %7.3f `z' _n								///
	   in gr _col(5) `"Prob > |z| = "'					///
	   in ye %8.4f 2*normprob(-abs(`z'))

	di _n in gr `"Aligned signed-rank (Hodges-Lehmann) estimate"' _n
	di in smcl in gr `"             {c |}  Estimate   [`level'% Conf. Interval]"'
	di in smcl in gr "{hline 13}{c +}{hline 33}"
	ditablin1 effect `effect' `lb_1' `ub_1'
	di in smcl in gr "{hline 13}{c BT}{hline 33}"

	/* return scalars for signrank */
	ret scalar sum_pos = `tp'
	ret scalar sum_neg = `tn'
	ret scalar z = `z'
	ret scalar Var_a = `v'
	ret scalar N_pos = `np'
	ret scalar N_neg = `nn'
	ret scalar N_tie = `n0'

	/* return scalars for H-L estimates */
	ret scalar estimate = `effect'
	ret scalar lb_1 = `lb_1'
	ret scalar ub_1 = `ub_1'
	ret scalar hl_obs = `nN'
	restore

end

program define ditablin
        di in smcl in gr %12s `"`1'"' `" {c |}"' in ye	///
                _col(17) %7.0g `2'						///
                _col(26) %10.0g `3'						///
                _col(38) %10.0g `4'
end

program define ditablin1
        di in smcl in gr %12s `"`1'"' `" {c |}"' in ye	///
                _col(17) %7.3g `2'						///
                _col(26) %10.3g `3'						///
                _col(38) %10.3g `4'
end