*! version 1.1.3  15dec2022 JM. Domenech, R. Sesma
/*
ciri: Confidence interval for Pearson/Spearman's correlation (immediate command)
If numlist has 2 elements, the CI for the correlation is computed
If numlist has 2 elements and the Rho option is specified, the correlation
			is compared to a hypothetical sample
If numlist has 4 elements, the CI for comparing 2 correlations is computed
*/

program define ciri, rclass
	version 12
	syntax anything(id="argument numlist"), /*
		*/ [PEarson SPearman Rho(numlist max=1 >=-1 <=1) /*
		*/	Level(numlist max=1 >50 <100) nst(string)]

	tempname r1 n1 r0 n0 w se rd t p lo up lo0 up0

	*Tokenize & test parameters
    local len : word count `anything'
	if (`len' != 2 & `len' != 4) print_error "argument numlist must have 2 or 4 elements"

	tokenize `anything'
	scalar `r1' = `1'
	scalar `n1' = `2'
	if (`r1'<-1 | `r1'>1) print_error "invalid data -- coefficient out of range"
	if (`n1'<=3) print_error "invalid data -- sample size out of range"
	local diff = (`len'==4)
	if (`diff') {
		scalar `r0' = `3'
		scalar `n0' = `4'
		if (`r0'<-1 | `r0'>1) print_error "invalid data -- coefficient out of range"
		if (`n0'<=3) print_error "invalid data -- sample size out of range"
	}

	local hypo = ("`rho'"!="")
	if (`len'==4 & `hypo'==1) print_error "too much options -- rho() is incompatible with comparing two correlations"
	if ("`pearson'"!="" & "`spearman'"!="") print_error "pearson and spearman options are incompatible"
	if ("`pearson'"=="" & "`spearman'"=="") local type = 1		//Pearson by default
	if ("`pearson'"!="") local type = 1
	if ("`spearman'"!="") local type = 2

	if ("`level'"=="") local level 95

	if (`type'==1) local title = "Pearson"
	if (`type'==2) local title = "Spearman"
	local ntab = length("`title'")+2

	*Print header
	di
	if (`diff'==0 & `hypo'==0) di as res "CONFIDENCE INTERVAL FOR " upper("`title'") " CORRELATION"
	if (`diff'==0 & `hypo'==1) di as res "COMPARING A CORRELATION WITH A HYPOTETICAL VALUE"
	if (`diff'==1) di as res "CONFIDENCE INTERVAL FOR COMPARING TWO " upper("`title'") " CORRELATIONS"
	if ("`nst'"!="") display as text "{bf:STUDY:} `nst'"

	*Compute & print results
	scalar `w' = ln(sqrt((1+`r1')/(1-`r1')))
	scalar `se'= 1/sqrt(`n1'-3)
	getci `w' `se', level(`level')
	scalar `lo' = r(lo)
	scalar `up' = r(up)

	if (`diff'==0 & `hypo'==0) {
		*Conf. Interval for one correlation
		di as txt "`title' Correlation: r = " as res %7.5f `r1'
        di as txt _col(`ntab') "Sample size: n = " as res `n1'
		di
		di as txt "`level'% CI of r:  " as res %7.5f `lo' as txt "  to  " as res %7.5f `up'

		*For one correlation, use Student's t for the signification
		scalar `t' = `r1'*sqrt(`n1'-2)/sqrt(1-`r1'^2)
		if (`r1'<1) scalar `p' = min(2*ttail(`n1'-2,abs(`t')),1)
		else if (`r1'>=1 & `r1'<.) scalar `p' = 0
		else if (`r1'==.) scalar `p' = .
		di as txt "Significance: t=  " as res %4.2f `t' as txt "    p= " as res %7.5f `p'

		return scalar rho = `r1'
		return scalar lb_rho = `lo'
		return scalar ub_rho = `up'
		return scalar p = `p'
		return scalar N = `n1'
	}
	if (`diff'==0 & `hypo'==1) {
		di as txt "Hypothetical: Rho = " as res %6.4f `rho'
		di as txt "{ralign 19:`title' Corr.: r =} " as res %6.4f `r1' /*
			*/	as txt " (`level'% CI: " as res %6.4f `lo' /*
			*/	as txt " to " as res %6.4f `up' as txt ")"
		di as txt "{ralign 19:Sample size: n =} " as res `n1'

		scalar `rd' = `r1' - `rho'
		scalar `w' = ln(sqrt((1+`r1')/(1-`r1'))) - ln(sqrt((1+`rho')/(1-`rho')))
		getci `w' `se', level(`level')
		di
		di as txt "Difference: r-Rho = " as res %7.5f `rd'
		di as txt " " %4.0g `level' "% CI of Diff.: " as res %7.5f r(lo) /*
			*/	as txt "  to  " as res %7.5f r(up)
		di as txt "Significance (two-sided test):  z=  " as res %4.2f r(z) as txt "  p= " as res %7.5f r(p)

		return scalar rho = `r1'
		return scalar lb_rho = `lo'
		return scalar ub_rho = `up'
		return scalar rd = `rd'
		return scalar lb_rd = r(lo)
		return scalar ub_rd = r(up)
		return scalar p = r(p)
		return scalar N = `n1'
	}
	if (`diff'==1) {
		di as txt "Sample 1: r1 = " as res %6.4f `r1' as txt /*
			*/	" (`level'% CI: " as res %6.4f `lo' as txt " to " /*
			*/ 	as res %6.4f `up' as txt ")"
		di as txt "          n1 = " as res `n1'

		scalar `w' = ln(sqrt((1+`r0')/(1-`r0')))
		scalar `se'= 1/sqrt(`n0'-3)
		getci `w' `se', level(`level')
		scalar `lo0' = r(lo)
		scalar `up0' = r(up)
		di as txt "Sample 0: r0 = " as res %6.4f `r0' as txt /*
			*/	" (`level'% CI: " as res %6.4f `lo0' as txt " to " /*
			*/ 	as res %6.4f `up0' as txt ")"
		di as txt "          n0 = " as res `n0'

		scalar `w' = ln(sqrt((1+`r1')/(1-`r1'))) - ln(sqrt((1+`r0')/(1-`r0')))
		scalar `se'= sqrt(1/(`n1'-3) + 1/(`n0'-3))
		scalar `rd'= `r1' - `r0'
		getci `w' `se', level(`level')
		di
		di as txt "Difference: r1-r0= " as res %7.5f `rd'
		di as txt "  `level'% CI of Diff.: " as res %7.5f r(lo) /*
			*/	as txt "  to  " as res %7.5f r(up)
		di as txt "Significance (two-sided test):  z=  " as res %4.2f r(z) as txt "  p= " as res %7.5f r(p)

		return scalar rho1 = `r1'
		return scalar lb_rho1 = `lo'
		return scalar ub_rho1 = `up'
		return scalar N1 = `n1'
		return scalar rho0 = `r0'
		return scalar lb_rho0 = `lo0'
		return scalar ub_rho0 = `up0'
		return scalar N0 = `n0'
		return scalar rd = `rd'
		return scalar lb_rd = r(lo)
		return scalar ub_rd = r(up)
		return scalar p = r(p)
	}
	return scalar level = `level'
end

program define getci, rclass
	syntax anything, level(numlist max=1 >50 <100)

	tempname w se z_alpha wl wu z

	tokenize `anything'
	scalar `w' = `1'
	scalar `se'= `2'

	scalar `z_alpha' = invnormal((`level'+100)/200)
	scalar `wl' = `w' - `z_alpha'*`se'
	scalar `wu' = `w' + `z_alpha'*`se'
	scalar `z' = `w'/`se'
	return scalar lo = (exp(2*`wl')-1)/(exp(2*`wl')+1)
	return scalar up = (exp(2*`wu')-1)/(exp(2*`wu')+1)
	return scalar z = `z'
	return scalar p = 2*(1-normal(abs(`z')))
end

program define print_error
	args message
	display in red "`message'"
	exit 198
end