// 1.0.1  JRF & AHF  01 FEB 2010
// Calculate Kendall's partial tau for an arbitrary number of confounders.

version 10.0
program define parttau, eclass sortpreserve

syntax varlist(numeric) [if] [in] [, CLuster(string) WCluster ///
			PC ESTimate(string) Level(cilevel) TRANSformation(string) ///
			SHOW Weight(string) RUN2 ]
marksample touse
	local vnum : word count `varlist'
			
	if `"`wcluster'"' == "" {
		local funtype = "bcluster"
	}
	else {
		local funtype = "wcluster"
	}
	if `"`transformation'"' == "" | `"`transformation'"' == "iden" {
		local transformation "identity"
		local transfText "None"
	}
	else if `"`transformation'"' == "sin" | `"`transformation'"' == "rho" {
		local transformation "rho"
		local transfText "Greiner's rho"
	}
	else if `"`transformation'"' == "arcsin" | `"`transformation'"' == "arsin" | `"`transformation'"' == "asin" {
		local transformation "asin"
		local transfText " Daniels' arcsine"
	}
	else if `"`transformation'"' == "z" {
		local transformation "z"
		local transfText "Fisher's z"
	}
	else if `"`transformation'"' == "zrho" {
		local transformation "zrho"
		local transfText "z-transform of Greiner's rho"
	}
	else if `"`transformation'"' == "c" {
		local transformation "c"
		local transfText "Harrell's c"
	}
	else {
		noi di ""
		noi di as error _column(4) "The given transformation does not match any of the offered transformations."
		noi di as error _column(4) "(Note that synonyms or truncations allowed in somersd might not be allowed here.)"
		exit
	}
		
	if `vnum' < 2 {
		noi di ""
		noi di _skip(4) as error "This command requires at least two variables."
		exit
	}
	else if `vnum' == 2 {
		local x : word 1 of `varlist'
		local y : word 2 of `varlist'
		qui replace `touse' = 0 if (`x' == . | `y' == . )
		qui count if `touse' == 1
		local N = r(N)
		
		capture somersd `x' `y' if `touse' != 0, taua level(`level') ///
					transf(`transformation') cluster(`cluster') funtype(`funtype')
		
		if _rc == 0 {			
			noi di ""
			noi di _skip(4) `"With just two variables, we'll run "somersd" with the "taua" option."'
			noi di ""
			
			somersd, level(`level')
		}
		else if `"`run2'"' != "" {
			local run2 = 1
		}
	}
	if `vnum' > 2 | `"`run2'"' != "" {
		unab variablelist : `varlist'
		tokenize `variablelist'

		qui count if `touse' == 1
		local N = r(N)

		local tistring ""
		local thetastring ""
		local numComps = comb(`vnum',2) + `vnum'
		forv i = 1(1)`numComps' {
			tempvar tiplus`i'
			qui gen `tiplus`i'' = .
			local tistring "`tistring'`tiplus`i'' "
			tempvar theta`i'
			qui gen `theta`i'' = .
			local thetastring "`thetastring'`theta`i'' "
		}
		if `"`cluster'"' == "" {
			tempvar clustertemp
			qui egen `clustertemp' = seq() if `touse'
			local nclusters = `N'
		}
		else {
			tempvar clustertemp
			qui egen `clustertemp' = group(`cluster') if `touse'
			qui summ `clustertemp'
			local nclusters = r(max)
		}
		sort `clustertemp'
		if `"`pc'"' == "" {
			local pc = 0
		}
		else {
			local pc = 1
		}
		if `"`estimate'"' == "" {
			tempvar estimate
			qui gen `estimate' = .
		}
		else {
			capture confirm variable `estimate'
			if _rc == 0 {
				noi di `"Note: the variable "`estimate'" already exists and is being replaced."'
			}
			else {
				qui gen `estimate' = .
			}
			qui replace `estimate' = .
			local estimate = "`estimate'"
		}
		if `"`level'"' != "" {
			local hival = 0.5 + `level' / 200
			local cifactor = invnormal(`hival')
		}
		else {
			local cifactor = 1.96
		}
		if `"`weight'"' == "" {
			tempvar weight
			qui gen `weight' = 1
			local weight "`weight'"
		}
		else {
			capture confirm numeric variable `weight'
			if _rc != 0 {
				noi di "{error}The given weight variable does not exist or is not numeric."
				exit
			}
			qui summ `weight'
			if r(min) < 0 {
				noi di "{error}Weights must be nonnegative."
			}
		}
		tempvar viplus
		qui gen `viplus' = .
		
		tempname pointEsts
		
		mata: main(`pc', "`clustertemp'", "`funtype'", `vnum', "`thetastring'", "`variablelist'", ///
					"`tistring'", "`viplus'", "`estimate'", "`weight'", "`transformation'", `cifactor', ///
					"`pointEsts'", "`touse'")

		matrix Sigma = J(`vnum',`vnum',.)
		local counter = 0
		forv i = 1(1)`vnum' {
			forv j = `i'(1)`vnum' {
				local counter = `counter' + 1
				matrix Sigma[`i',`j'] = taus[1,`counter']
				matrix Sigma[`j',`i'] = taus[1,`counter']
			}
		}
		matrix drop taus
			
		tempname primaryvars confounders
		local `primaryvars' = word("`variablelist'",1) + " " + word("`variablelist'",2)
		local `confounders' ""
		forv i=3(1)`vnum' {
			local `confounders' = "``confounders'' " + word("`variablelist'",`i')
		}
		
		noi di ""
		noi di "{text}Primary variables: {result}``primaryvars''"
		noi di "{text}Confounders: {result}``confounders''"
		noi di "{text}Transformation: {result}`transfText'"
		noi di "{text}No. of obs.: {result}`N'"
		if `"`cluster'"' != "" {
			noi di "{text}No. of clusters: {result}`nclusters'"
		}
		noi di ""
		
		if `"`show'"' != "" {
			noi di "{text}Pairwise taus:"
			noi di "{text}{hline 48}"
			noi di "{text}{col 3}#{col 10}Name{col 40}Estimate"
			noi di "{text}{hline 48}"
			forv i = 1(1)`vnum' {
				local name1 = word("`variablelist'",`i')
				local length1 = length("`name1'")
				local dilen1 = min(`length1' + 2,8)
				local name1 = substr("`name1'",1,`dilen1')
				forv j = `i'(1)`vnum' {
					local name2 = word("`variablelist'",`j')
					local length2 = length("`name2'")
					local dilen2 = min(`length2' + 2,8)
					local name2 = substr("`name2'",1,`dilen2')
					noi di "{text}{col 2}`i',`j'{col 10}tau( " ///
							%`dilen1's "`name1', " %`dilen2's "`name2' )" ///
							_column(39) %10.7f as result Sigma[`i',`j']
				}
			}
			noi di "{text}{hline 48}"
			noi di ""
			noi di "{text}Partial tau:"
		}
				
		local rows = rowsof(`pointEsts')
		
		noi di "{text}{hline 13}{c TT}{hline 67}"
		noi di "{col 14}{text}{c |}{text}{col 31}Jackknife"
		noi di "{text}Interval:{col 14}{text}{c |}{col 22}Coef.{col 31}Std. Err.{col 46}z{col 52}P>|z|{col 62}[95% Conf. Interval]"
		noi di "{text}{hline 13}{c +}{hline 67}"
		noi di "{text}   symmetric{col 14}{text}{c |}{col 17}{result}" %10.7f `pointEsts'[1,1] _column(30) %10.7f `pointEsts'[1,2]  ///
				_column(43) %5.2f `pointEsts'[1,1]/`pointEsts'[1,2] ///
				_column(52) %5.3f 2 * (1 - normal(abs(`pointEsts'[1,1]/`pointEsts'[1,2]))) ///
				_column(61) %8.6f `pointEsts'[1,3] "    " %8.6f `pointEsts'[1,4]
		if "`transformation'" != "identity" & "`transformation'" != "rho"  & `"`transformation'"' != "c" {
			noi di "{text}  asymmetric{col 14}{text}{c |}{col 17}{result}" %10.7f `pointEsts'[2,1] as result ///
						_column(61) %8.6f `pointEsts'[2,3] "    " %8.6f `pointEsts'[2,4]
		}
		noi di "{text}{hline 13}{c BT}{hline 67}"
		if `"`transformation'"' != "identity" {
			noi di ""
			if "`transformation'" == "rho" | "`transformation'" == "c" {
				noi di as text "Results for partial tau are in the transformed space."
			}
			else if "`transformation'" == "zrho" {
				noi di as text "1st line is in the transformed, zrho, space."
				noi di as text "2nd line is in the space of rho.
			}
			else {
				noi di as text "1st line is in the transformed space."
				noi di as text "2nd line is in the original space.
			}
		}
		noi di ""
		
		tempname colnames colsofb
		local `colnames' "tau"
		
		matrix V = (`pointEsts'[1,2]^2)
		matrix colnames V = "``colnames''"
		matrix rownames V = "``colnames''"
		matrix colnames b = "``colnames''"
		
		matrix rownames Sigma = `varlist'
		matrix colnames Sigma = `varlist'
		
		ereturn post b V, esample(`touse')
		ereturn matrix taus = Sigma
		
		ereturn scalar SE = `pointEsts'[1,2]
		ereturn scalar N = `N'
		ereturn scalar N_clust = `nclusters'
		
		ereturn local cluster_type = `"`funtype'"'
		ereturn local primaryvars = "``primaryvars''"
		ereturn local confounders = "``confounders''"
		ereturn local vcetype = "Jackknife"
		ereturn local transf = "`transformation'"
		ereturn local cmd = "npt"
	}
end

mata:
	void main(real scalar pc, string scalar clusterv, string scalar funtype, ///
				real scalar nvars, string scalar thetas, string scalar varnames, ///
				string scalar tipluses, string scalar viplusv, string scalar estimatev, ///
				string scalar weightv, string scalar transformation, real scalar cifactor, ///
				string scalar matname, string scalar tousev) 
	{
		st_view(cluster, ., clusterv, tousev)
		clpanel = panelsetup(cluster, 1)
		clstats = panelstats(clpanel)
		
		st_view(V, ., tokens(varnames), tousev)
		st_view(Ti, ., tokens(tipluses), tousev)
		st_view(Th, ., tokens(thetas), tousev)
		st_view(estimate, ., estimatev, tousev)
		st_view(weight, ., weightv, tousev)
		st_view(viplus, ., viplusv, tousev)
		
		numComps = comb(nvars,2) + nvars
		
		taus = J(1,numComps,.)
		
		viplus[.] = makeweight(weight)
		
		allphis = J(clstats[1],2 * (numComps + 1), .)
		
		counter = 0		
		for (i1 = 1; i1 <= nvars; i1++) {
			for (i2 = i1; i2 <= nvars; i2++) {
				counter = counter + 1
				st_subview(x, V, ., i1)
				st_subview(y, V, ., i2)
				st_subview(tiplus, Ti, ., counter)
				st_subview(thetaj, Th, ., counter)
				
				tiplus[.] = maketiplus(x, y)
				phis = makephi(cluster, tiplus, x, y, viplus, weight)
				if (nvars == 2) {
					colOne = 2 * counter - 1
					colTwo = 2 * counter
					allphis[., colOne::colTwo] = phis[.,1::2]
				}
				if (funtype == "bcluster") {
					taus[counter] = (colsum(phis[.,1]) - colsum(phis[.,2])) / (colsum(phis[.,3]) - colsum(phis[.,4]))
				}
				else if (funtype == "wcluster") {
					taus[counter] = colsum(phis[.,2]) / colsum(phis[.,4])
				}
				if (i1 == i2 & pc == 1) {
					for (i3 = 1; i3 <= clstats[1]; i3++) {
						thetaj[i3] = 1
					}
				}
				else {
					thetaj[1..clstats[1]] = makethetaj(phis, funtype)
				}
			}
		}
		st_matrix("taus",taus)
		if (nvars == 2) {
			allphis[.,7::8] = phis[.,3::4]
			influence = makeInfluence(allphis, funtype)
			dispersion = makeDispersion(influence, allphis, clstats[1], funtype, transformation)
			st_matrix("disp",dispersion)
			results = makeResults2by2(dispersion, taus, cifactor, transformation)
		}
		else {
			estimate[.] = makeSigma(Th, nvars)
			overallEst = makeSigma(taus, nvars)
			results = makeResults(overallEst, estimate, cifactor, transformation)
		}
		st_matrix(matname, results)
	}
	
	real colvector makeweight(real colvector weight)
	{
		nobs = rows(weight)

		viplus = J(nobs,1,.)
		for (i1 = 1; i1 <= nobs; i1++) {
			viplus[i1] = 0
			for (i2 = 1; i2 < i1; i2++) {
				viplus[i1] = viplus[i1] + 1
				viplus[i2] = viplus[i2] + 1
			}
		}
		return(viplus)
	}
	
	real colvector maketiplus(real colvector x, real colvector y)
	{
		nobs = rows(x)
		
		tiplus = J(nobs,1,.)
		for (i1 = 1; i1 <= nobs; i1++) {
			tiplus[i1] = 0
			x1 = x[i1]
			y1 = y[i1]
			for (i2 = 1; i2 < i1; i2++) {
				x2 = x[i2]
				y2 = y[i2]
				xsign = (x1 < x2) - (x2 < x1)
				ysign = (y1 < y2) - (y2 < y1)
				tiplus[i1] = tiplus[i1] + xsign * ysign
				tiplus[i2] = tiplus[i2] + xsign * ysign
			}
		}
		return(tiplus)
	}
	
	real matrix makephi(real colvector cluster, real colvector tiplus, ///
						real colvector x, real colvector y, ///
						real colvector viplus, real colvector weight)
	{
		clpanel = panelsetup(cluster, 1)
		clstats = panelstats(clpanel)
				
		phikk = J(clstats[1],1,.)
		phijplus = J(clstats[1],1,.)
		
		phiVkk = J(clstats[1],1,.)
		phiVjplus = J(clstats[1],1,.)
		
		phis = J(clstats[1],4,.)
				
		for (i1 = 1; i1 <= clstats[1]; i1++) {
			cl_tiplus = tiplus[clpanel[i1,1]..clpanel[i1,2]]
			cl_viplus = viplus[clpanel[i1,1]..clpanel[i1,2]]
			phijplus[i1] = colsum(cl_tiplus[.])
			phiVjplus[i1] = colsum(cl_viplus[.])
			
			cl_x = x[clpanel[i1,1]..clpanel[i1,2]]
			cl_y = y[clpanel[i1,1]..clpanel[i1,2]]
			wcl_tiplus = J(rows(cl_x), 1, .)
			wcl_tiplus[.] = maketiplus(cl_x, cl_y)
			
			cl_weight = weight[clpanel[i1,1]..clpanel[i1,2]]
			wcl_viplus = J(rows(cl_weight), 1, .)
			wcl_viplus[.] = makeweight(cl_weight)
			
			phikk[i1] = colsum(wcl_tiplus[.])
			phiVkk[i1] = colsum(wcl_viplus[.])
		}
		
		phis[.,1] = phijplus[.]
		phis[.,2] = phikk[.]
		phis[.,3] = phiVjplus[.]
		phis[.,4] = phiVkk[.]
				
		return(phis)
	}
	
	real colvector makethetaj(real matrix phis, string scalar funtype)
	{
		phikk = phis[.,2]
		phiVkk = phis[.,4]
		
		phikksum = colsum(phikk[.])
		phiVkksum = colsum(phiVkk[.])
		
		thetaj = J(rows(phis),1,.)
		if (funtype == "bcluster") {
			phijplus = phis[.,1]
			phiVjplus = phis[.,3]
			phiplusplus = colsum(phijplus[.])
			phiVplusplus = colsum(phiVjplus[.])
			for (i1 = 1; i1 <= rows(phis); i1++) {
				thetajNumer = (phiplusplus - phikksum) - 2 * (phijplus[i1] - phikk[i1])
				thetajDenom = (phiVplusplus - phiVkksum) - 2 * (phiVjplus[i1] - phiVkk[i1])
				thetaj[i1] = thetajNumer / thetajDenom
			}
		}
		else if (funtype == "wcluster") {
			for (i1 = 1; i1 <= rows(phis); i1++) {
				thetajNumer = (phikksum - phikk[i1])
				thetajDenom = (phiVkksum - phiVkk[i1]) 
				thetaj[i1] = thetajNumer / thetajDenom
			}
		}
		return(thetaj)
	}
	
	real colvector makeSigma(real matrix Th, real scalar nvars)
	{
		Sigma = J(nvars, nvars, .)
		estimate = J(rows(Th),1,.)

		for (i1 = 1; i1 <= rows(Th); i1++) {
			counter = 0
			for (i2 = 1; i2 <= nvars; i2++) {
				for (i3 = i2; i3 <= nvars; i3++) {
					counter = counter + 1
					Sigma[i2,i3] = Th[i1,counter]
					Sigma[i3,i2] = Th[i1,counter]
				}
			}
			if (nvars > 2) {
				Sigma11 = Sigma[1..2,1::2]
				Sigma12 = Sigma[1..2,3::nvars]
				Sigma21 = Sigma[3..nvars,1::2]
				Sigma22 = Sigma[3..nvars,3::nvars]
				invSigma22 = invsym(Sigma22)
				
				Sigma_11dot2 = Sigma11 - Sigma12 * invSigma22 * Sigma21
				estimate[i1] = Sigma_11dot2[1,2] / sqrt(Sigma_11dot2[1,1] * Sigma_11dot2[2,2])
			}
			else {
				estimate[i1] = Sigma[1,2]
			}
		}
		return(estimate)
	}
	
	real matrix makeInfluence(real matrix allPhis, string scalar funtype)
	{
		nrows = rows(allPhis)
		ncols = cols(allPhis)
		
		phiVkk = allPhis[., ncols]
		phiXkk = allPhis[., 2]
		phiYkk = allPhis[., 4]
		
		psiV = J(nrows,1,.)
		psiX = J(nrows,1,.)
		psiY = J(nrows,1,.)
		
		if (funtype == "bcluster") {
			phiVjplus = allPhis[., ncols - 1]
			phiVplusplus = colsum(allPhis[., ncols - 1])
			phiVkksum = colsum(allPhis[., ncols])
			phiXjplus = allPhis[., 1]
			phiXplusplus = colsum(allPhis[., 1])
			phiXkksum = colsum(allPhis[., 2])
			phiYjplus = allPhis[., 3]
			phiYplusplus = colsum(allPhis[., 3])
			phiYkksum = colsum(allPhis[., 4])
			for (i1 = 1; i1 <= nrows; i1++) {
				psiV[i1] = (phiVplusplus - phiVkksum) / (nrows - 1) ///
							- ((phiVplusplus - phiVkksum) - 2 * (phiVjplus[i1] - phiVkk[i1])) / (nrows - 2)
				psiX[i1] = (phiXplusplus - phiXkksum) / (nrows - 1) ///
							- ((phiXplusplus - phiXkksum) - 2 * (phiXjplus[i1] - phiXkk[i1])) / (nrows - 2)
				psiY[i1] = (phiYplusplus - phiYkksum) / (nrows - 1) ///
							- ((phiYplusplus - phiYkksum) - 2 * (phiYjplus[i1] - phiYkk[i1])) / (nrows - 2)
			}
		}
		if (funtype == "wcluster") {
			for (i1 = 1; i1 <= nrows; i1++) {
				psiV[i1] = phiVkk[i1]
				psiX[i1] = phiXkk[i1]
				psiY[i1] = phiYkk[i1]
			}
		}
		
		influence = J(nrows,3,.)
		influence[.,1] = psiV :- (colsum(psiV[.]) / nrows)
		influence[.,2] = psiX :- (colsum(psiX[.]) / nrows)
		influence[.,3] = psiY :- (colsum(psiY[.]) / nrows)
		
		return(influence)
	}
	
	real matrix makeDispersion(real matrix influence, real matrix allPhis, real scalar n_cl, ///
							   string scalar funtype, string scalar transformation)
	{
		if (funtype == "bcluster") {
			Txx = (colsum(allPhis[.,1]) - colsum(allPhis[.,2])) / (n_cl * (n_cl - 1))
			Txy = (colsum(allPhis[.,3]) - colsum(allPhis[.,4])) / (n_cl * (n_cl - 1))
			V = (colsum(allPhis[.,7]) - colsum(allPhis[.,8])) / (n_cl * (n_cl - 1))
		}
		if (funtype == "wcluster") {
			Txx = colsum(allPhis[.,2]) / n_cl
			Txy = colsum(allPhis[.,4]) / n_cl
			V = colsum(allPhis[.,8]) / n_cl
		}
		derivs = (- Txx / (V ^ 2), 1 / V, 0      \  ///
				  - Txy / (V ^ 2), 0    , 1 / V)
		dispersion = derivs * influence' * influence * derivs' / (n_cl * (n_cl - 1))
		
		return(dispersion)
	}
	
	real matrix makeResults2by2(real matrix dispersion, real matrix taus, ///
									real scalar cifactor, string scalar transformation)
	{
		tausApprox = J(1,cols(taus),.)
		for (i = 1; i <= cols(taus); i++) {
			if (taus[i] == 1) {
				tausApprox[i] = 0.999999
			}
			else if (taus[i] == -1) {
				tausApprox[i] = -0.999999
			}
			else {
				tausApprox[i] = taus[i]
			}
		}
		
		if (transformation == "identity") {
			zeta = tausApprox[2]
			lopt = zeta - cifactor * sqrt(dispersion[2,2])
			hipt = zeta + cifactor * sqrt(dispersion[2,2])
			returnMat = (zeta, sqrt(dispersion[2,2]), lopt, hipt)
		}
		else if (transformation == "z") {
			transDerivs = (1 / (1 - tausApprox[1]^2) , 0  \  ///
						   0, 1 / (1 - tausApprox[2]^2))
			transDispersion = transDerivs * dispersion * transDerivs
			zeta = atanh(tausApprox[2])
			slopt = zeta - cifactor * sqrt(transDispersion[2,2])
			shipt = zeta + cifactor * sqrt(transDispersion[2,2])
			alopt = tanh(zeta - cifactor * sqrt(transDispersion[2,2]))
			ahipt = tanh(zeta + cifactor * sqrt(transDispersion[2,2]))
			returnMat = (zeta, sqrt(transDispersion[2,2]), slopt, shipt \
						 taus[2], 						., alopt, ahipt)
		}
		else if (transformation == "asin") {
			transDerivs = (1 / ((1 - tausApprox[1]^2)^(1/2)), 0  \  ///
						   0, 1 / ((1 - tausApprox[2]^2)^(1/2)))
			transDispersion = transDerivs * dispersion * transDerivs
			zeta = asin(tausApprox[2])
			slopt = zeta - cifactor * sqrt(transDispersion[2,2])
			shipt = zeta + cifactor * sqrt(transDispersion[2,2])
			alopt = sin(zeta - cifactor * sqrt(transDispersion[2,2]))
			ahipt = sin(zeta + cifactor * sqrt(transDispersion[2,2]))
			returnMat = (zeta, sqrt(transDispersion[2,2]), slopt, shipt \
						 taus[2], 						., alopt, ahipt)
		}
		else if (transformation == "rho") {
			transDerivs = (pi() / 2 * cos(pi() * tausApprox[1] / 2), 0   \  ///
						   0, pi() / 2 * cos(pi() * tausApprox[2] / 2))
			transDispersion = transDerivs * dispersion * transDerivs
			zeta = sin(pi() / 2 * tausApprox[2])
			slopt = zeta - cifactor * sqrt(transDispersion[2,2])
			shipt = zeta + cifactor * sqrt(transDispersion[2,2])
			returnMat = (zeta, sqrt(transDispersion[2,2]), slopt, shipt)
		}
		else if (transformation == "zrho") {
			transDerivs = ( pi() / 2 * cos(pi() * tausApprox[1] / 2) / (1 - sin(pi() * tausApprox[1] / 2)^2), 0  \  ///
							0, pi() / 2 * cos(pi() * tausApprox[2] / 2) / (1 - sin(pi() * tausApprox[2] / 2 )^2))
			transDispersion = transDerivs * dispersion * transDerivs
			zeta = atanh(sin(pi() / 2 * tausApprox[2]))
			slopt = zeta - cifactor * sqrt(transDispersion[2,2])
			shipt = zeta + cifactor * sqrt(transDispersion[2,2])
			rho = sin(pi() / 2 * tausApprox[2])
			alopt = tanh(zeta - cifactor * sqrt(transDispersion[2,2]))
			ahipt = tanh(zeta + cifactor * sqrt(transDispersion[2,2]))
			returnMat = (zeta, sqrt(transDispersion[2,2]), slopt, shipt \ rho, ., alopt, ahipt)
		}
		else if (transformation == "c") {
			transDerivs = ( 1/2, 0  \  ///
							0, 1/2)
			transDispersion = transDerivs * dispersion * transDerivs
			zeta = (tausApprox[2] + 1) / 2
			slopt = zeta - cifactor * sqrt(transDispersion[2,2])
			shipt = zeta + cifactor * sqrt(transDispersion[2,2])
			returnMat = (zeta, sqrt(transDispersion[2,2]), slopt, shipt)
		}
		st_matrix("b",(zeta))
		return(returnMat)
	}
	
	real matrix makeResults(real scalar overallEst, real colvector estimate, ///
									real scalar cifactor, string scalar transformation)
	{
		if (abs(overallEst) - 1 <= 1e-6) {
			approxTau = 0.999999 * overallEst
		}
		numEsts = rows(estimate) - missing(estimate)
		jkEstimates =  (estimate :* (1 - numEsts)) :+ (numEsts * overallEst)
		jkVar = variance(jkEstimates) / numEsts
		if (transformation == "identity") {
			zeta = overallEst
			zetaVec = jkEstimates
			SE = sqrt(jkVar * 1)
			lopt = zeta - cifactor * SE
			hipt = zeta + cifactor * SE
			returnMat = (zeta, SE, lopt, hipt)
		}
		else if (transformation == "z") {
			zeta = atanh(overallEst)
			zetaVec = atanh(jkEstimates)
			SE = sqrt(jkVar * 1 / (1 - approxTau^2)^2)
			slopt = zeta - cifactor * SE
			shipt = zeta + cifactor * SE
			alopt = tanh(zeta - cifactor * SE)
			ahipt = tanh(zeta + cifactor * SE)
			returnMat = (zeta, SE, slopt, shipt \
						 overallEst, ., alopt, ahipt)
		}
		else if (transformation == "asin") {
			zeta = asin(overallEst)
			zetaVec = asin(jkEstimates)
			SE = sqrt(jkVar * 1 / (1 - approxTau^2))
			slopt = zeta - cifactor * SE
			shipt = zeta + cifactor * SE
			alopt = sin(zeta - cifactor * SE)
			ahipt = sin(zeta + cifactor * SE)
			returnMat = (zeta, SE, slopt, shipt \
						 overallEst, ., alopt, ahipt)
		}
		else if (transformation == "rho") {
			zeta = sin(overallEst * pi() / 2)
			zetaVec = sin(jkEstimates :* pi() / 2)
			SE = sqrt(jkVar * (pi() / 2 * cos(pi() / 2 * approxTau))^2)
			slopt = zeta - cifactor * SE
			shipt = zeta + cifactor * SE
			returnMat = (zeta, SE, slopt, shipt)
		}
		else if (transformation == "zrho") {
			zeta = atanh(sin(overallEst * pi() / 2))
			zetaVec = atanh(sin(jkEstimates :* pi() / 2))
			SE = sqrt(jkVar * (pi() / 2 * cos(pi() / 2 * approxTau) / (1 - (sin(pi() / 2 * approxTau))^2))^2)
			slopt = zeta - cifactor * SE
			shipt = zeta + cifactor * SE
			rho = sin(overallEst * pi() / 2)
			alopt = tanh(zeta - cifactor * SE)
			ahipt = tanh(zeta + cifactor * SE)
			returnMat = (zeta, SE, slopt, shipt \ rho, ., alopt, ahipt)
		}
		else if (transformation == "c") {
			zeta = (overallEst + 1) / 2
			zetaVec = (jkEstimates :+ 1) :/ 2
			SE = sqrt(jkVar / 4)
			slopt = zeta - cifactor * SE
			shipt = zeta + cifactor * SE
			returnMat = (zeta, SE, slopt, shipt)
		}
		st_matrix("b",(zeta))
		return(returnMat)
	}
end