Optimal Visual Testing Confidence Level Identification

VizTest.do

* Choosing the Optimal Confidence Level for Visual Testing
* David A. Armstrong II and William Poirier 
* 
* Requires Stata >= 16 because frames are used to save results.
* Requires the function matsort be installed with:
*     net install matsort.pkg
capture program drop testCI
program testCI
	syntax , [lev1(real .25) lev2(real .99) incr(real .01) a(real .05) easythresh(real .05) inc0 remc usemargins]
	if `lev1' <= 0  | `lev1' >= 1 {
		di "lev1 must be a real value between 0 and 1, exclusive. Please respecify lev1"
		break
	}
	if `lev2' <= 0  | `lev2' >= 1 {
		di "lev1 must be a real value between 0 and 1, exclusive. Please respecify lev2"
		break
	}
	local ldist = `lev2' - `lev1'
	if `incr' >= `ldist' {
		di "incr is intended to be the step increasing moving from lev1 to lev2 and should be much smaller than the distance between lev1 and lev2.  Please respecify incr""
		break
	}
	if `a' > .25 {
		if `a' < 1 { 
			di "The a parameter is the type I error on the individual pairwise tests.  You have specified a type I error rate of `a', are you sure that's right?"
		}
		else{
			di "The a parameter is the type I error on the individual pairwise tests and should be between 0 and 1, exclusive.  Please respecify a." 
			break
		}
	}

	if "`usemargins'" != "" {
		capture local rb = rowsof(r(b))
		if "`rb'" != "" {
			mat tmp = r(b)'
			mat tmpV = r(V)
		} 
		else{
			mat tmp = e(b)'
			mat tmpV = e(V)
		}
	}
	else{
		mat tmp = e(b)'
		mat tmpV = e(V)
	}

	mat bhat = tmp
	mat V = tmpV

	
	local kv = rowsof(V)
	local kv2 = `kv'
	local k = rowsof(bhat)
	local k2 = `k'
	if "`inc0'" != "" {
		local k2 = `k2' + 1
		local kv2 = `kv2' + 1
		mat zmat = J(1,1,0)
		mat rownames zmat = zero
		mat bhat = bhat \ zmat
		mat zrow = J(1, `kv', 0)
		mat rownames zrow = zero
		mat zcol = J(`kv2', 1, 0)
		mat colnames zcol = zero
		mat V = V \ zrow
		mat V = V , zcol
	}

	local nrt = rowsof(bhat)
	
	mata: rn = st_matrixrowstripe("bhat")
	mata: rnms = J(0,1,"")
	forv i = 1/`nrt'{
		mata: tmprnm = invtokens(rn[`i', ], "")
		mata: tmprnm = regexr(tmprnm, "\.", "")
		mata: rnms = rnms \ tmprnm
	}
	mata: rnms = rnms'
	mata: rnmsl = invtokens(rnms[1, ], " ")
	mata: j = rnmsl[1,1]
	mata: st_local("rnms", j)
	mat myV = J(rowsof(V), colsof(V), 0)
	mat myB = J(rowsof(bhat), colsof(bhat), 0)
	mat rownames myB = `rnms'
	mat rownames myV = `rnms'
	mat colnames myV = `rnms'
	local nr = rowsof(bhat)
	forv i = 1/`nr'{ 
		local bb = bhat[`i',1]
		mat myB[`i',1] = `bb'
		forv j = 1/`nr' {
			local vv = V[`i',`j']
			mat myV[`i',`j'] = `vv'
		}
	}
	mat bhat = myB
	mat V = myV
	local rnb: rownames bhat
	if "`remc'" != "" {
		mat f = J(1,1,0)
		local nwb : word count `rnb'
		tokenize `rnb'
		forv s = 1/`nwb'{
			if regexm("``s''", ".*_cons.*") == 0{
				mat f = f \ bhat[`s', ....]
			}
		}
		local rf = rowsof(f)
		mat f = f[2..`rf', ....]
		mat bhat = f
	}
	
	matsort bhat 1 "down"
	local k = rowsof(bhat)
	local rn : rownames(bhat)
	tokenize `rn'
	mat newV = V[....,"`1'"]
	forval i = 2/`k'{
		mat zz = V[....,"``i''"]
		mat newV = newV, zz
	}
	local newrn : colnames(newV)
	tokenize `newrn'
	mat newV2 = newV["`1'",....]
	forval i = 2/`k'{
		mat zzz = newV["``i''",....]
		mat newV2 = newV2 \ zzz
	}
	mat V = newV2
	
	local np = `k'*(`k'-1)/2
	local resdf = _N-rowsof(tmp)
	mat pairs = J(`np', 2, 0)

	local cl = 1
	matrix D = J(`k', `np', 0)
	local upri = `k' -1
	forval i = 1/`upri' {
		local lwrj = `i' + 1
		forval j = `lwrj'/`k' {
			matrix D[`i', `cl'] = 1
			matrix D[`j', `cl'] = -1
			mat pairs[`cl', 1] = `i'
			mat pairs[`cl', 2] = `j'
			local cl = `cl' + 1
		}
	}


	mata: D = st_matrix("D")
	mata: V = st_matrix("V")
	mata: bhat = st_matrix("bhat")
	mata: pairs = st_matrix("pairs")



	mata: diffs = (bhat'*D)'
	mata: diffV = D'*V*D
	mata: diffSE = sqrt(diagonal(diffV))
	local crit = invt(`resdf', 1-`a')
	mata: pwtests = diffs:/diffSE :>`crit'
	mata: ones = J(`np', 1, 1)
	mata: st_numscalar("nsig", pwtests'*ones)
	local qtl = 1-(1-`lev1')/2
	mata: LL = bhat[,1] - invt(`resdf', `qtl'):*sqrt(diagonal(V))[,1]
	mata: UU = bhat[,1] + invt(`resdf', `qtl'):*sqrt(diagonal(V))[,1]
	mata: levs = range(`lev1', `lev2', `incr')
	mata: st_matrix("levs", levs)
	local nlevs = rowsof(levs)
	forval l = 2/`nlevs'  {
		local qtl = 1-(1-levs[`l',1])/2
		mata: LL = LL, bhat[,1] - invt(`resdf', `qtl'):*sqrt(diagonal(V))[,1]
		mata: UU = UU, bhat[,1] + invt(`resdf', `qtl'):*sqrt(diagonal(V))[,1]
		}
	mata: maxU = colmax(UU)
	mata: minL = colmin(LL)
	mata: thresh = (maxU :- minL):*`easythresh'
	mata: st_matrix("thresh", thresh)
	mata: L = LL[pairs[,1], ]
	mata: U = UU[pairs[,2], ]
	mata: citests = L :>= U
	mata: comptests = citests :== pwtests
	mata: ones = J(`np', 1, 1)


	mata: d = (L-U):*pwtests + (U-L):*(1:-pwtests)
	mata: st_matrix("d", d)
	local rd = rowsof(d)
	local cold = colsof(d)
	forv j = 1/`cold'{
		local thrsh = thresh[1,`j']
		forv i = 1/`rd' {
			local dij = d[`i',`j']
			if `dij' > `thrsh' {
				mat d[`i',`j'] = `thrsh'
			}
			if `dij' < -`thrsh'{
				mat d[`i',`j'] = -`thrsh'
			}
		}
	}
	mata: d = st_matrix("d")
	mata: easy = colsum(d)'
	mata: res = levs, (comptests'*ones):/`np', easy
	mata: st_matrix("res", res)
	mata: cm = colmax(res)
	mata: st_numscalar("cm", cm[1,2])
	mata: st_numscalar("maxeasy", cm[1,3])

	mat res_s = J(1,3,0)
	local nlev = rowsof(res)
	forval i = 1/`nlev' {
		local psm = res[`i', 2]
		if `psm' == cm {
			mat res_s = res_s \ res[`i',....]
		}
	}
	local nrs = rowsof(res_s)
	mat res_s = res_s[2..`nrs', ....]
	mata: res_s = st_matrix("res_s")
	mata: cms = colmax(res_s)
	mata: st_numscalar("maxeasy", cms[1,3])
	local nrs = rowsof(res_s)
	forval i = 1/`nrs'{
		if res_s[`i',3] == maxeasy {
			local optlev = res_s[`i',1]
		}
	}
	mata: levnum = range(1, `nlevs', 1)
	mata: levcol = (levs :== `optlev')'*levnum
	mata: alltests = pairs, pwtests, citests[,levcol]
	mata: st_matrix("alltests", alltests)
	mat miss = J(1,4,0)
	local nrm = rowsof(alltests)
	forval i = 1/`nrm' {
		if alltests[`i',3] != alltests[`i',4] { 
			mat miss = miss \ alltests[`i',....]
		}
	}
	local nrm = rowsof(miss)
	if `nrm' > 1 {
	mat miss = miss[2..`nrm', ....]
	local nmt = rowsof(miss)
	mata: miss = st_matrix("miss")
	mata: varnames = st_matrixrowstripe("bhat")[,2]
	mata: sigmat = ("Insignificant" \ "Significant")
	mata: olmat = ("Overlapping" \ "Not overlapping")
	mata: miss_tests = varnames[miss[,1],1], varnames[miss[,2],1], sigmat[miss[,3]:+1, 1], olmat[miss[,4]:+1, 1]
	}
	di " "
	di "Optimal Levels: "
	di " "
	mata: res_s
	di " "
	di "Missed Tests (n=`nmt' of `np')"
	di " "
	if `nrm' > 1{
	mata: miss_tests
	}
	else{
		di "No missed tests!"
	}
	frame create levels level pr_same easiness
	frame change levels 
	mata: st_addobs(rows(res))
	mata: st_store(., "level", res[,1])
	mata: st_store(., "pr_same", res[,2])
	mata: st_store(., "easiness", res[,3])
	local fname = "result-levels-$S_TIME-$S_DATE.dta"
	save "`fname'", replace
	frame change default
	frame drop levels
	if `nr' > 1 {	
		frame create missed 
		frame change missed
		gen str25 bigger = ""
		gen str25 smaller = ""
		gen str25 pw_test = ""
		gen str25 ci_test = ""
		mata: st_addobs(rows(miss_tests))
		mata: st_sstore(., "bigger", miss_tests[,1])
		mata: st_sstore(., "smaller", miss_tests[,2])
		mata: st_sstore(., "pw_test", miss_tests[,3])
		mata: st_sstore(., "ci_test", miss_tests[,4])
		local fname = "result-missed-$S_TIME-$S_DATE.dta"
		save "`fname'", replace
		frame change default
		frame drop missed
	}
	
end