*! version 3.0 10jul2019, Benjamin R. Shear and sean f. reardon

* requires hetop_lf.ado

 capture program drop hetop
 program define hetop , eclass

	version 13.1

	syntax namelist(min=2 max=2) , 		/// arg1 is the group id variable; arg2 is the stem for the category counts variable
		NUMCats(integer)				/// technically redundant but a good check that there are the expected number of categories
		[ 											///
			MODtype(string) 						/// "hetop" or "homop"
			IDentify(string)						/// "sums" or "refgroup" or "cuts"
			SETRef(integer 0) 						/// group to use as reference group, if identify == refgroup
			PHOP(string)							/// 0/1 indicator; constrain groups with `phop'==1 to have a common SD, can optionally specify "mean"
			SETCuts(numlist ascending min=1 max=9)	/// set known cut scores
			PKVALS(varname numeric max=1) 			/// name of variable containing group proportions in population; default is to use observed group proportions
			STARTFRom(namelist max=1)				/// matrix with desired starting values
			KAPPA(integer 1)					 	/// use a specified value of kappa in calculations
			SAVE(string)							/// if specified, save estimates as new variables in sheet with names like mstarstring; optional to include se's
			ADDCONStraints(string)					/// use pre-defined constraints
			MINSIZE(integer 1)						/// drop groups smaller than this
			INITVALS 								/// compute initial values
			GAPS 									/// return matrix of gaps and SEs of gaps
			NOIsily		 							/// display -oglm- output during estimation
			HOMOP									/// easier way to specify a homop model
			CSD(real -1)							/// constrain common SD to this value for HOMOP models; default will change to 1 if left as -1
			*										/// ML options
		]

	preserve

	tokenize `namelist'
	local grpid `1'
	local catname `2'

	mlopts mlopts , `options'

	// set default options:

	if "`modtype'" == "" {
		if "`homop'" == "homop"	local modtype "homop"
		if "`homop'" == ""		local modtype "hetop"
	}
	if "`modtype'" == "hetop" & "`homop'" == "homop" {
		noi di in red "WARNING: modtype(hetop) and homop options are not " ///
		"consistent. homop model will be fit."
		local modtype "homop"
	}
	if `setref' > 0 & "`identify'" == ""	local identify "refgroup"
	if "`identify'" == ""					local identify "refgroup"
	local ref_group "n"		// default, only gets reset if identify=refgroup

	// evaluate and parse syntax

	* if setcuts is specified, ignore any other identification options
	if "`setcuts'" != "" {

		local identify "setcuts"

		local num_setcuts : word count `setcuts'
		if `=`numcats'-1' != `num_setcuts' {
			noi di as error "number of cutscores != K-1"
			error 499
		}

		forv i = 1/`num_setcuts' {
			local set_c`i' : word `i' of `setcuts'
		}

		* maybe add if `num_setcuts' == 1 {local modtype "homop"}

	}

	* verify that id variable integer > 0 and uniquely identifies groups
	qui capture assert `grpid' == int(`grpid')
	qui sum `grpid'
	if r(min) < 1 | _rc != 0 {
		noi di in red "ID variable must be positive integer values only"
		error 499
	}

	cap bys `grpid' : assert _N==1
	if _rc != 0 {
		noi di in red 	"ID variable `grpid' does not uniquely "///
						"identify observations"
		error 499
	}

	* can't do HETOP with fewer than 3 categories
	if "`modtype'" == "hetop" & `numcats' < 3 {
		noi di in red "ERROR: must have > 2 categories to fit HETOP model."
		error 499
	}

	* verify that variable save names not taken
	if "`save'" != "" {
		parse_save_opt `save'
		local newvarnames "`s(allvarnames)'"
		foreach v in `newvarnames' {
			cap confirm variable `v' , exact
			if _rc == 0 {
				noi di in red "variable `v' exists but save option " ///
					"specified to use this name. please rename `v' or alter " ///
					"save variable name options."
				error 499
			}
		}
	}

	* confirm identification
	if inlist("`identify'", "sums", "refgroup", "cuts", "setcuts") == 0 {
		noi di in red "ERROR: identify() option must be 'sums', 'refgroup'," ///
		" or 'cuts'."
		error 499
	}

	* confirm modeltype selection
	if inlist("`modtype'", "hetop", "homop") == 0 {
		noi di in red "ERROR: modtype() must be 'hetop' or 'homop'."
		error 499
	}

	* confirm phopvar is only 0 and 1 and not specified with homop option
	* ignore phop if fitting model with a single group
	if _N == 1 {
		if "`phop'" != "" {
			noi di in yellow "NOTE: only 1 group; ignoring PHOP variable."
		}
		local phop ""
	}
	if "`phop'" != "" & "`modtype'" == "homop" {
		noi di in red "WARNING: phop option with homop model is redundant. " ///
			"Ignoring phop variable."
		local phop ""
	}
	if "`phop'" != "" {
		parse_phop_opt `phop'
		local phop `s(phopvar)'
		local phop_mean `s(phop_mean)'
		qui count if (`phop'!=0 & `phop'!=1)
		if r(N) > 0 {
			noi di in red "ERROR: phop variable can only contain 0 or 1."
			error 499
		}
	}

	* can't specify STARTFROM and INITVALS
	if "`startfrom'" != "" {
		if "`initvals'" != "" {
			noi di in red "ERROR: cannot specify both startfrom() and initvals."
			error 499
		}
		cap confirm matrix `startfrom'
		if _rc != 0 {
			noi di in red "ERROR: `startfrom' must be a matrix."
			exit _rc
		}
	}

	// check and clean frequency count data

	* replace any missing frequencies with 0s
	qui foreach var of varlist `catname'* {
		qui count if `var' == .
		if r(N) > 0 noi di in yellow ///
			"NOTE: replaced " r(N) " missing frequencies in `var' with 0."
		replace `var' = 0 if `var' == .
	}

	* verify category name & number of category variables match
	unab vars : `catname'*
	local numcatsfound = `: word count `vars''
	if `numcatsfound' != `numcats' {
		noi di in red 	"ERROR: Number of categories found does not match " ///
						"number supplied to numcats."
		error 499
	}
	local numcuts = `numcats' - 1

	// restrict sample and issue warnings about existence

	tempvar Nk
	egen `Nk' = rowtotal(`catname'*)

	* drop groups with no observations or that are below `minsize'
	qui count if `Nk' < `minsize'	// `minsize'=1 by default
	if r(N) > 0 noi di in red _n "WARNING: dropping " r(N) ///
		" groups with N<`minsize'"
	qui drop if `Nk' < `minsize'

	* check for categories with no observations across entire sample
	foreach var of varlist `catname'* {
		qui sum `var'
		if (r(sum) <= 0 | r(sum) == .) & "`identify'" != "setcuts" {
			noi di in red "ERROR: `var' contains no observations (possibly" ///
			" after removing groups below the minsize() threshold."		///
			" All categories must have at least 1 observation"	///
			" in the total sample."
			error 499
		}
		if (r(sum) <= 0 | r(sum) == .) & "`identify'" == "setcuts" {
			noi di in red "WARNING: at least one category has no " ///
			"observations, but estimation with fixed cut scores will proceed."
		}
	}

	* issue warning about groups with data in 2 cells for HETOP
	tempvar num0
	egen `num0' = anycount(`catname'*) , val(0)
	qui count if `num0' == (`numcats'-2)
	if "`modtype'" == "hetop" & r(N) > 0 {
		di in red "WARNING: some groups have data in only 2 " ///
		"categories and ML estimates may not be properly defined."
	}

	* now issue warnings/errors for groups with data in only 1 category
	qui count if `num0' == (`numcats'-1)

	* error for 1-count groups with only 2 categories and HOMOP
	if "`modtype'" == "homop" & r(N) > 0 & `numcats' == 2 {
		noi di in red "ERROR: cannot fit HOMOP model when some groups have " ///
			"data in only 1 category and numcats=2. There are " r(N) ///
			" such groups."
		error 499
	}
	else if r(N) > 0 {
		// this scenario possible with PHOP
		noi di in red "WARNING: some groups have data in only 1 " ///
		"category and ML estimates may not be properly defined."
	}


					// !! sample is now finalized !!  //


	// create additional variables and subset to analytic sample

	* sort by grpid, restrict variables
	sort `grpid'
	qui sum `grpid'
	local K = r(N)		// number of groups
	keep `grpid' `catname'* `pkvals' `Nk' `phop'
	sort `grpid'

	* create pk values if none are supplied
	if "`pkvals'" == "" {
		tempvar pk
		qui sum `Nk'
		gen double `pk' = `Nk' / r(sum)
		local pkvals "`pk'"
		noi di in yellow "NOTE: no pk values specified; " ///
			"using nk/N as proportions."
	}

	sort `grpid'

	* these will be used below

	tempname P nk nkHI
	mkmat `pkvals' , matrix(`P')
	matrix `P' = `P''
	mkmat `Nk' , matrix(`nk')
	matrix `nk' = `nk''
	matrix `nkHI' = J(1 , `K' , 0)
	forv i = 1/`K' {
		matrix `nkHI'[1,`i'] = 1/`nk'[1,`i']
	}

	// set up identifying constraints

	if "`identify'" == "refgroup" {

		* check specified ref group or automatically generate one
		if `setref' == 0 {
			noi di in yellow "NOTE: no reference group (or 0) specified; " ///
				"automatic group will be specified."
			qui get_ref , refrank(1) catname(`catname') grpid(`grpid') ///
							numcats(`numcats')
			local ref_group = r(refid)
		}

		else {
			qui count if `grpid' == `setref'
			if r(N) != 1 {
				noi di "specified reference group `setref' not found"
				error 499
			}
			local ref_group = `setref'
			noi di in yellow "NOTE: user specified `setref' as reference group."
		}
	}
	else if "`identify'" == "cuts" {

		noi di in yellow "NOTE: setting cut scores for model identification."

		* if only 1 cut score, or HOMOP fix to 0
		if `numcats' == 2 & "`modtype'" == "homop" {
			constraint free
			constraint `=r(free)' [cut1]_cons = 0
			local cutconstraints "`=r(free)'"
		}
		else if `numcats' > 2 & "`modtype'" == "homop" {
			constraint free
			constraint `=r(free)' [cut2]_cons = 0
			local cutconstraints "`=r(free)'"
		}
		else {
			* if 2 or more and HETOP, fix to -1 and 0
			constraint free
			constraint `=r(free)' [cut1]_cons = -1
			local cutconstraints "`=r(free)'"
			constraint free
			constraint `=r(free)' [cut2]_cons = 0
			local cutconstraints "`cutconstraints' `=r(free)'"
		}

	}
	else if "`identify'" == "sums" {

		noi di in yellow "NOTE: using constraint on sum of estimates " ///
			"to identify model."

		qui levelsof `grpid' , local(grpconstraints)

		local w1 `: word 1 of `grpconstraints''
		local mnconstraints "`=`P'[1,1]' * [mean]:`w1'.`grpid'"
		local sdconstraints "`=`P'[1,1]' * [lnsigma]:`w1'.`grpid'"

		forv i = 2/`K' {

			local w1 `: word `i' of `grpconstraints''
			local mnconstraints "`mnconstraints' + `=`P'[1,`i']' * [mean]:`w1'.`grpid'"
			local sdconstraints "`sdconstraints' + `=`P'[1,`i']' * [lnsigma]:`w1'.`grpid'"

		}

		local mnconstraints "`mnconstraints' = 0"
		local sdconstraints "`sdconstraints' = 0"

		constraint free
		local pmc = r(free)
		constraint `pmc' `mnconstraints'

		if "`modtype'" == "hetop" {
			constraint free
			local smc = r(free)
			constraint `smc' `sdconstraints'
		}

		local sumconstraints "`pmc' `smc'"
	}
	else if "`identify'" == "setcuts" {

		local cutconstraints ""
		forv i = 1/`=`numcats'-1' {
			constraint free
			constraint `=r(free)' [cut`i']_cons = `set_c`i''
			local cutconstraints "`cutconstraints' `=r(free)'"
		}

	}
	***

	* add phop constraints if phop indicator supplied

	if "`phop'" != "" & "`modtype'" == "hetop" {

		qui count if `phop' == 1
		if r(N) == _N {
			noi di in red "WARNING: Constant `phop' means this is a HOMOP model " ///
				"and should probably be respecified."
		}
		qui count if `phop' == 0
		if r(N) == _N {
			noi di in red "WARNING: `phop'=0 for all observations. Model is " ///
				"equivalent to HETOP and will ignore `phop'."
			local phop ""
			local phop_mean = 0
		}

		if "`phop'" != "" {

			qui levelsof `grpid' if `phop' == 1 , local(grpconstraints)
			qui count if `phop' == 1
			local numconstraints = r(N) - 1
			local grpconstraints "`grpconstraints'"

			forv i = 1/`numconstraints' {

				local w1 `: word `i' of `grpconstraints''
				local w2 `: word `=`i'+1' of `grpconstraints''
				capture constraint free

					if _rc != 0 {
						noi di in red "ERROR: attempting to set too many " ///
							"constraints with the phop command. " ///
							"May need to use an alternative model."
						error _rc
						exit
					}

				local smc = r(free)
				constraint `smc' [lnsigma]:`w1'.`grpid' = [lnsigma]:`w2'.`grpid'
				local phopconstraints "`phopconstraints' `smc'"

			}
		}

		if "`phop_mean'" == "1" {

			qui count if `phop' == 1

			if r(N) < _N {

				// define extra row of Cns matrix for the mean constraint

				tempname cnsaddrow
				tempvar phopfree_cns
				qui g double `phopfree_cns' = 0
				qui count if `phop' == 0
				scalar sumwt = 1/`r(N)'
				qui replace `phopfree_cns' = -1*scalar(sumwt) if `phop'==0
				qui sum `grpid' if `phop' == 1
				qui replace `phopfree_cns' = 1 if `grpid' == r(min)
				qui mkmat `phopfree_cns' , mat(`cnsaddrow')
				mat `cnsaddrow' = `cnsaddrow''

				tempname cnsaddrow_l cnsaddrow_r
				mat `cnsaddrow_l' = J(1,`K',0)
				mat `cnsaddrow_r' = J(1,`=`numcuts'+1',0)
				mat `cnsaddrow' = `cnsaddrow_l', `cnsaddrow', `cnsaddrow_r'

				local cnsaddrow_true = 1

			}
		}
	}
	*****


	****************************************************************************
	* set up equations and globals for hetop_lf, run -ml-

	local eq1label mean
	local eq2label lnsigma

	unab vlist : `catname'*
	local vlist : subinstr local vlist "`catname'" "" , all
	numlist "`vlist'" , sort
	local j=1
	foreach i in `r(numlist)' {
		global FW_`j' "`catname'`i'"
		local j=`j'+1
	}

	* equation for means

	local eqs (`eq1label':ib`ref_group'.`grpid', noconstant)
	ereturn local mean_eq "(`eq1label':ib`ref_group'.`grpid', noconstant)"

	* equation for ln(SD)
	if "`modtype'" == "homop" {
		if "`setcuts'" == "" {
			global fixedsd = 1
			if `csd' < 0 local csd = 1
			global clsd = `=ln(`csd')'
			local num_b_parms = `=`K'+`numcuts''
		}
		else {
			if `num_setcuts' == 1 {
				global fixedsd = 1
				if `csd' < 0 local csd = 1
				global clsd = `=ln(`csd')'
				local num_b_parms = `=`K'+`numcuts''
			}
			else {
				if `csd' < 0 {
					local csd = 1
					global fixedsd = 0
					local eqs 				 "`eqs' (`eq2label': )"
					ereturn local lnsigma_eq "(`eq2label': )"
					local num_b_parms = `=`K'+1+`numcuts''
				}
				else {
					global fixedsd = 1
					global clsd = `=ln(`csd')'
					local num_b_parms = `=`K'+`numcuts''
				}
			}
		}
	}
	else if "`modtype'" == "hetop" {
		local csd = 1
		global fixedsd = 0

		local eqs 		   "`eqs' (`eq2label':ib`ref_group'.`grpid', noconstant)"
		ereturn local lnsigma_eq "(`eq2label':ib`ref_group'.`grpid', noconstant)"

		local num_b_parms = `=2*`K'+`numcuts''
	}

	* cut score equations
	forval i = 1/`numcuts' {
		local eqs "`eqs' (cut`i':)"
	}

	local allconstraints ///
		`refgroup_constraints' `sumconstraints' `phopconstraints' ///
		`cutconstraints' `homopconstraints' `addconstraints'

	if "`allconstraints'" != "" | "`cnsaddrow_true'" != "" {

		// set up Cns constraints matrix

		tempname cnsb cnsV cnsmat

		local cnsb_m_names
		local cnsb_s_names
		local cnsb_c_names

		qui levelsof `grpid' , local(grpid_levels)

		// note: names must reflect the ib`refgroup' equation in order to
		// correctly set constraints in makecns
		foreach lev in `grpid_levels' {
			if "`lev'" != "`ref_group'" {
				local cnsb_m_names `cnsb_m_names' `eq1label':`lev'.`grpid'
			}
			else {
				local cnsb_m_names `cnsb_m_names' `eq1label':`lev'b.`grpid'
			}
		}
		foreach lev in `grpid_levels' {
			if "`lev'" != "`ref_group'" {
				local cnsb_s_names `cnsb_s_names' `eq2label':`lev'.`grpid'
			}
			else {
				local cnsb_s_names `cnsb_s_names' `eq2label':`lev'b.`grpid'
			}
		}
		forv i = 1/`numcuts' {
			local cnsb_c_names `cnsb_c_names' cut`i':_cons
		}

		// dummy solutions
		if "`modtype'" == "hetop" {
			mat `cnsb' = J(1,`num_b_parms',0)
			mat colnames `cnsb' = `cnsb_m_names' `cnsb_s_names' `cnsb_c_names'
		}
		else if "`modtype'" == "homop" {
			mat `cnsb' = J(1,`num_b_parms',0)
			if `num_b_parms' == `=`K'+`numcuts'' {
				mat colnames `cnsb' = `cnsb_m_names' `cnsb_c_names'
			}
			else if `num_b_parms' == `=`K'+1+`numcuts'' {
				mat colnames `cnsb' = `cnsb_m_names' `eq2label':_cons `cnsb_c_names'
			}
		}

		mat `cnsV' = `cnsb''*`cnsb'
		ereturn post `cnsb' `cnsV'

		makecns `allconstraints'
		mat `cnsmat' = e(Cns)

		// add extra phop,mean row if necessary
		if "`cnsaddrow_true'" == "1" {
			mat `cnsmat' = `cnsmat' \ `cnsaddrow'
		}

		local constraintcommand "constraints(`cnsmat')"

	}

	****************************************************************************
	* NOTES to user

	if "`modtype'" == "hetop" {

		if "`phop'" != "" di in yellow "NOTE: fitting heteroskedastic " ///
			"ordered probit model with constraints in `phop'."

		if "`phop'" == "" di in yellow "NOTE: fitting heteroskedastic " ///
			"ordered probit model."

		*local scale_equation "hetero(ib`ref_group'.`grpid')"

	}

	else if "`modtype'" == "homop" {

		di in yellow "NOTE: fitting homoskedastic ordered probit model " ///
			"with common ln(SD) set to `=ln(`csd')'"

	}

	di in yellow "NOTE: user specified --`grpid'-- as group id variable."
	di in yellow "NOTE: user specified --`catname'*-- as category variables."
	di in yellow "NOTE: user specified `numcats' categories."
	if "`initvals'" != "" {
		di in yellow "NOTE: custom starting values will be used."
	}

	****************************************************************************
	* starting values

	tempname initvalsmat
	if "`startfrom'" != "" {
		matrix `initvalsmat' = `startfrom'
	}
	else {
		di in yellow "calculating initial values..."
		get_initial_values , grpid(`grpid') refgroup(`ref_group') ///
					catname(`catname') numgrps(`K') modtype(`modtype') ///
					numcuts(`numcuts') weighttype(iweight) ///
					identify(`identify') ///
					setcuts(`setcuts') `initvals'
		mat `initvalsmat' = r(initvalsmat)
	}
	local startvalcall "init(`initvalsmat', skip)"

	****************************************************************************
	* estimate model

	* get total sample size
	qui sum `Nk'
	local totaln = r(sum)

	* workaround to estimate fixed cut score model with single group
	* warning: single-group estimation not recommended
	if _N == 1 {
		noi di in red "WARNING: only 1 group used for estimation."
		qui expand 2
		qui foreach var of varlist `catname'* {
			replace `var' = `var'/2
		}
		local docollapse 1
	}

	di in yellow "...estimating model..."

	capture `noisily' ml model lf hetop_lf `eqs' , ///
		maximize ///
		`constraintcommand' ///
		`startvalcall' ///
		`mlopts' ///
		obs(`totaln') ///
		wald(0)

	if "`docollapse'" == "1" {
		qui collapse (sum) `catname'*  (mean) `Nk' `pk' , by(`grpid')
	}

	local numcuts2 `numcuts'
	local numcuts3 `numcuts'
	ereturn scalar k_aux = `numcuts2'
	ereturn scalar k_exp = `numcuts3'
	ereturn scalar k_eq_model = e(k_eq) - e(k_aux)
	ereturn local equation "`eqs'"

	****************************************************************************

	// process and return results

	if _rc != 0 {

		noi di in red "WARNING: estimation error code _rc = " _rc
		constraint drop `allconstraints'
		macro drop fixedsd FW_* clsd
		error 499

	}

	else if _rc == 0 {

		di in yellow "...estimation done."

		if e(converged) == 0 noi di in red "WARNING: model failed to converge."

		sort `grpid'

		** SE calculations and dereferencing

		tempvar ninv1 ninv

		tempname Vfull Bfull P2 One PI n ntilde Q a ///
			Mprime Gprime Sprime Cprime ///
			Mprime_se Sprime_se Gprime_se Cprime_se ///
			Vprime Omegaprime Lambdaprime Wprime Zprime ///
			sigmaw sigmab sigmaprime ///
			Mstar Sstar Mstar_se Sstar_se Cstar Cstar_se ///
			icchatratio icchat R T varsigprime Vstar Wstar Zstar icchatvar ///
			S G Gvar1 ///
			Braw Bprime Aprime Dprime U Bstar Astar Dstar ///
			Mraw Sraw Mraw_se Sraw_se Craw Craw_se ///
			Zraw Wraw Vraw

		matrix `Vfull' = e(V)
		matrix `Bfull' = e(b)'

		matrix `P2' = hadamard(`P', `P')
		matrix `One' = J(1,`K',1)
		matrix `PI' = I(`K') - `One''*`P'

		gen double `ninv1' = 1 / (`Nk' - 1)
		mkmat `ninv1' , mat(`ninv1')
		drop `ninv1'

		mkmat `Nk' , mat(`n')
		matrix `n' = `n''

		* adjust the "omega-hat-bar-G" (ohbg) term for model type

		qui if "`modtype'" == "hetop" {
			matrix `ntilde' = ( (1/`K') * (`One' * `ninv1') )
			matrix `ntilde' = invsym(`ntilde')
			local ohbg = 1 / (2*`ntilde'[1,1])
			mat drop `ntilde'
		}

		qui if "`modtype'" == "homop" {
			qui sum `Nk'
			if $fixedsd == 0 local ohbg = 1 / (2 * (r(sum)-`K'))
			if $fixedsd == 1 local ohbg = 0
		}

		qui if "`modtype'" == "hetop" & "`phop'" != "" {
			tempvar nt
			g double `nt' = `Nk' - 1
			qui sum `nt' if `phop' == 1
			replace `nt' = r(sum) if `phop' == 1
			replace `nt' = 1/(2*`nt')
			qui sum `nt'
			local ohbg = (1/`K')*r(sum)
		}

		g double `ninv' = 1/`Nk'
		mkmat `ninv' , mat(`ninv')
		drop `ninv'

		matrix `Q' = (1/(1+(2*`ohbg'))) * ///
			hadamard( hadamard( `ninv'' , (`P' + `n' - `One') ) , `P' )

		/*

		create the following vectors/matrices.

		Mprime: estimated means in estimation metric.
		Vprime: estimated sampling variance/covariance of the Mprime vector
		Gprime: estimated ln(SD) in estimation metric.
		Omegaprime: estimated sampling variance/covariance of Gprime vector
		Lambdaprime: sampling covariances of Mprime and Gprime
		Cprime: estimated thresholds
		Bprime: sampling variance/covariance of thresholds
		Sprime: exp(Gprime)
		Wprime: estimated sampling variance/covariance of Sprime
		Zprime: estimated covariance of Mprime and Sprime
		Aprime: sampling var/cov of means X cuts
		Dprime: sampling var/cov of gammas X cuts

		*/

		qui levelsof `grpid' , local(nameid)
		local mnames
		local snames
		forv i = 1/`K' {
			local j : word `i' of `nameid'
			local mnames	"`mnames' `j'.`grpid'"
			local snames 	"`snames' `j'.`grpid'"
		}

		matrix `Mprime' = `Bfull'[1 .. `K', 1]
		matrix `Vprime' = `Vfull'[1..`K',1..`K']

		if "`modtype'" == "hetop" {

			matrix `Gprime' = `Bfull'[`=`K'+1' .. `=`K'*2', 1]
			matrix `Omegaprime' = `Vfull'[`=`K'+1'..`=2*`K'',`=`K'+1'..`=2*`K'']
			matrix `Lambdaprime' = `Vfull'[1..`K',`=`K'+1'..`=2*`K'']
			matrix `Cprime' = `Bfull'[`=`K'*2+1' .. `=`K'*2+`numcuts'' , 1]
			matrix `Bprime' = `Vfull'[`=`K'*2+1' .. `=`K'*2+`numcuts'' , ///
				`=`K'*2+1' .. `=`K'*2+`numcuts'']
			matrix `Aprime' = `Vfull'[1..`K',`=`K'*2+1' .. `=`K'*2+`numcuts'']
			matrix `Dprime' = `Vfull'[`=`K'+1'..`=2*`K'', ///
				`=`K'*2+1' .. `=`K'*2+`numcuts'']

		}

		else {
			// else model is homop
			if $fixedsd == 1 {
				matrix `Gprime' = J(`K',1,`=${clsd}')
				matrix `Omegaprime' = J(`K',`K',0)
				matrix `Lambdaprime' = J(`K',`K',0)
				matrix `Cprime' = `Bfull'[`=`K'+1' .. `=`K'+`numcuts'' , 1]
				matrix `Bprime' = `Vfull'[`=`K'+1' .. `=`K'+`numcuts'' , ///
				                             `=`K'+1' .. `=`K'+`numcuts'']
				matrix `Aprime' = `Vfull'[1..`K',`=`K'+1' .. `=`K'+`numcuts'']
				matrix `Dprime' = J(`K',`numcuts',0)
			}
			else if $fixedsd == 0 {
				// only here if "`identify'"=="setcuts"
				local sdest = `Bfull'[`=`K'+1',1]
				matrix `Gprime' = J(`K',1,`sdest')
				matrix `Omegaprime' = J(`K',`K',`Vfull'[`=`K'+1',`=`K'+1'])
				matrix `Lambdaprime' = `Vfull'[`=`K'+1',1..`K']
				matrix `Lambdaprime' = `Lambdaprime''*`One'
				matrix `Cprime' = `Bfull'[`=`K'+2' .. `=`K'+`numcuts'+1' , 1]
				matrix `Bprime' = `Vfull'[`=`K'+2' .. `=`K'+`numcuts'+1' , ///
				                          `=`K'+2' .. `=`K'+`numcuts'+1']
				matrix `Aprime' = J(`K',`numcuts',0)
				matrix `Dprime' = J(`K',`numcuts',0)
			}

		}

		* store the means and SDs in the estimation "raw" metric

		matrix `Mraw' = `Mprime'
		matrix `Vraw' = `Vprime'
		matrix `Mraw_se' = vecdiag(`Vraw')'

		matrix roweq	`Mraw' = mean:
		matrix rownames `Mraw' = `mnames'
		matrix roweq	`Mraw_se' = mean:
		matrix rownames `Mraw_se' = `mnames'

		matrix roweq	`Vraw' = mean:
		matrix coleq	`Vraw' = mean:
		matrix rownames `Vraw' = `mnames'
		matrix colnames `Vraw' = `mnames'

		matrix `Sraw' = `Gprime'
		forv k = 1/`K' {
			matrix `Sraw'[`k',1] = exp(`Sraw'[`k',1])
			matrix `Mraw_se'[`k',1] = `Mraw_se'[`k',1]^.5
		}

		matrix `Wraw' = diag(`Sraw') * `Omegaprime' * diag(`Sraw')
		matrix `Sraw_se' = vecdiag(`Wraw')'
		forv k = 1/`K' {
			matrix `Sraw_se'[`k',1] = `Sraw_se'[`k',1]^.5
		}

		matrix roweq	`Wraw' = sigma:
		matrix coleq	`Wraw' = sigma:
		matrix rownames `Wraw' = `snames'
		matrix colnames `Wraw' = `snames'

		matrix `Zraw' = `Lambdaprime'*diag(`Sraw')

		matrix roweq	`Zraw' = mean:
		matrix coleq	`Zraw' = sigma:
		matrix rownames `Zraw' = `mnames'
		matrix colnames `Zraw' = `snames'

		matrix roweq	`Sraw' = sigma:
		matrix rownames `Sraw' = `snames'
		matrix roweq	`Sraw_se' = sigma:
		matrix rownames `Sraw_se' = `snames'

		matrix `Braw' = `Bprime'
		matrix `Craw' = `Cprime'
		matrix `Craw_se' = vecdiag(`Bprime')'
		forv i = 1/`=`numcats'-1' {
			mat `Craw_se'[`i',1] = `Craw_se'[`i',1]^.5
		}

		if "`identify'" != "sums" {

			// convert from double-prime/raw to single-prime metric if

			tempname A OneK1 MprimeS
			matrix `A' = `P' * `Gprime'
			local a = `A'[1,1]

			local kap = 1
			if `kappa' == 2 {
				matrix `varPG' = `P' * `Omegaprime' * `P''
				local kap = (1 + (1/2) * `varPG'[1,1])
			}

			local delta = exp(`a')
			matrix `OneK1' = J(1,`numcuts',1)
			matrix `MprimeS' = (exp(`a')^-1) * ( `PI' * `Mprime' )	// single prime Mprime vector

			matrix `Cprime' = (`delta'^(-1))*(`Cprime' - `OneK1''*(`P'*`Mprime'))	// convert to single prime Cprime

			matrix `Bprime' = ///
				exp(-2*`a')* ///
				(`kap'^2 *(`Bprime'-(`Aprime''*`P'')*`OneK1'-`OneK1''*(`P'*`Aprime')+`OneK1''*(`P'*`Vprime'*`P'')*`OneK1') ///
				-(`delta'/`kap')*hadamard(`Cprime'*`OneK1',(`OneK1''*`P'*`Dprime'-`OneK1''*(`P'*`Lambdaprime'*`P'')*`OneK1')) ///
				-(`delta'/`kap')*hadamard(`OneK1''*`Cprime'',(`Dprime''*`P''*`OneK1'-`OneK1''*(`P'*`Lambdaprime'*`P'')*`OneK1')) ///
				+((`delta'^2)/(`kap'^4))*hadamard(hadamard((`Cprime'*`OneK1'),(`OneK1''*`Cprime'')),(`OneK1''*(`P'*`Omegaprime'*`P'')*`OneK1')))

			matrix `Aprime' = (1/`delta'^2) * (`kap'^2) * ///
				(`PI'*(`Aprime'-`Vprime'*`P''*`OneK1')) ///
				- (`delta'*`kap')^(-1)*hadamard(`MprimeS'*`OneK1',`One''*`P'*`Dprime'-`One''*(`P'*`Lambdaprime'*`P'')*`OneK1') ///
				- (`delta'*`kap')^(-1)*hadamard(`One''*`Cprime'', `PI'*`Lambdaprime'*`P''*`OneK1') ///
				+ (`kap')^(-4)*hadamard(hadamard(`MprimeS'*`OneK1', `One''*`Cprime''), `One''*(`P'*`Omegaprime'*`P'')*`OneK1')

			matrix `Dprime' = (`kap'/`delta') * ///
				(`PI'*(`Dprime'-`Lambdaprime''*`P''*`OneK1')) - ///
				`kap'*hadamard(`One''*`Cprime'', `PI'*(`Omegaprime''*`P''*`OneK1'))

			matrix `Vprime' = (exp(`a')^-2) * ( ///
											( `kap'^2 * `PI' * `Vprime' * `PI'' ) ///
											- (`kap'^(-1)) * (`PI' * `Mprime' * `P' * `Lambdaprime'' * `PI'' + `PI' * `Lambdaprime' * `P'' * `Mprime'' * `PI'' ) ///
											+ (`kap'^(-4)) * (`PI' * `Mprime' * `Mprime'' * `PI'' * ( `P' * `Omegaprime' * `P'') ) ///
											)

			matrix `Lambdaprime' = (exp(`a')^-1) * ( ( `kap' * `PI' * `Lambdaprime' * `PI'') - (`kap'^(-2)) * `PI' * `Mprime' * `P' * `Omegaprime' * `PI'' )

			matrix `Mprime' = (exp(`a')^-1) * ( `PI' * `Mprime' )
			matrix `Gprime' = `PI' * `Gprime'
			matrix `Omegaprime' = `PI' * `Omegaprime' * `PI''

		}

		* create Sprime, Wprime and Zprime

		matrix `Sprime' = `Gprime'
		forv k = 1/`K' {
			matrix `Sprime'[`k',1] = exp(`Sprime'[`k',1])
		}

		matrix roweq	`Mprime' = mean:
		matrix rownames `Mprime' = `mnames'
		matrix roweq	`Vprime' = mean:
		matrix coleq	`Vprime' = mean:
		matrix rownames `Vprime' = `mnames'
		matrix colnames `Vprime' = `mnames'

		matrix roweq	`Aprime' = mean:
		matrix rownames	`Aprime' = `mnames'

		matrix roweq	`Sprime' = sigma:
		matrix rownames	`Sprime' = `snames'

		matrix `Wprime' = diag(`Sprime') * `Omegaprime' * diag(`Sprime')
		matrix `Zprime' = `Lambdaprime' * diag(`Sprime')

		if `K' > 1 {	// standardizing, etc. not appropriate when number of groups = 1

			* --------------------------------- *
			* Mstar and Sstar
			* --------------------------------- *

			matrix `sigmaw' = (1/(1+(2*`ohbg'))) * (`P' * hadamard( `Sprime' , `Sprime' ))
			matrix `sigmab' = (`P' * hadamard( `Mprime' , `Mprime' )) + ( (1/(1+(2*`ohbg'))) * ( hadamard( `ninv'' , (`P2' - `P') ) * hadamard( `Sprime' , `Sprime' ) ) )

			matrix `sigmaprime' = cholesky( `P' * hadamard(`Mprime',`Mprime') + `Q' * hadamard(`Sprime', `Sprime') )

			matrix `Mstar' = invsym(`sigmaprime')*`Mprime'
			matrix `Sstar' = invsym(`sigmaprime')*`Sprime'
			matrix `Cstar' = invsym(`sigmaprime')*`Cprime'

			matrix `icchatratio' = `sigmab' * invsym(`sigmaw' + `sigmab')
			matrix `icchat' = 1 - ((1/(1+(2*`ohbg'))) * `P' * hadamard(`Sstar', `Sstar'))

			* --------------------------------- *
			* Vstar and Wstar and Zstar
			* --------------------------------- *

			capture {

				matrix `R' = `P' * diag(`Mstar') * `Vprime' + `Q' * diag(`Sstar') * `Zprime''
				matrix `T' = `P' * diag(`Mstar') * `Zprime' + `Q' * diag(`Sstar') * `Wprime'
				matrix `varsigprime' = invsym(hadamard(`sigmaprime', `sigmaprime')) * (`P' * diag(`Mprime') * `Vprime' * diag(`Mprime') * `P'' + `Q' * diag(`Sprime') * `Wprime' * diag(`Sprime') * `Q'' + 2 * `P' * diag(`Mprime') * `Zprime' * diag(`Sprime') * `Q'')
				matrix `Vstar' = invsym(hadamard(`sigmaprime', `sigmaprime')) * (`Vprime' - (`Mstar' * `R' + `R'' * `Mstar'') + `Mstar' * `Mstar'' * `varsigprime')
				matrix `Wstar' = invsym(hadamard(`sigmaprime', `sigmaprime')) * (`Wprime' - (`Sstar' * `T' + `T'' * `Sstar'') + `Sstar' * `Sstar'' * `varsigprime')
				matrix `Zstar' = invsym(hadamard(`sigmaprime', `sigmaprime')) * (`Zprime' - (`Mstar' * `T' + `R'' * `Sstar'') + `Mstar' * `Sstar'' * `varsigprime')

				matrix `icchatvar' = 4*(1/(1+(2*`ohbg')))^2 * `P'*(diag(`Sstar')*`Wstar'*diag(`Sstar'))*`P''

				matrix `U' = `P'*diag(`Mstar')*`Aprime'+`Q'*diag(`Sstar')*`Dprime'

				matrix `Bstar' = invsym(hadamard(`sigmaprime', `sigmaprime'))*(`Bprime'-(`Cstar'*`U'+`U''*`Cstar'')+(`Cstar'*`Cstar'')*`varsigprime')
				matrix `Astar' = invsym(hadamard(`sigmaprime', `sigmaprime'))*(`Aprime'-(`Mstar'*`U'+`R''*`Cstar'')+`Mstar'*`Cstar''*`varsigprime')
				matrix `Dstar' = invsym(hadamard(`sigmaprime', `sigmaprime'))*(`Dprime'-(`Sstar'*`U'+`T''*`Cstar'')+`Sstar'*`Cstar''*`varsigprime')

			}
			*****

			local varmatsRC = _rc		// record errors

			* standard errors
			local mseRC = 0
			cap matrix `Mstar_se' = vecdiag(cholesky(diag(vecdiag(`Vstar'))))'
			if _rc != 0 {
				matrix `Mstar_se' = J(`K' , 1 , .)
				local mseRC = _rc
			}

			local cseRC = 0
			cap matrix `Cstar_se' = vecdiag(cholesky(diag(vecdiag(`Bstar'))))'
			if _rc != 0 {
				matrix `Cstar_se' = J(`numcuts' , 1 , .)
				local cseRC = _rc
			}

			local sseRC = 0
			cap matrix `Sstar_se' = vecdiag(cholesky(diag(vecdiag(`Wstar'))))'
			if _rc != 0 {
				matrix `Sstar_se' = J(`K' , 1 , .)
				local sseRC = _rc
			}

			* --------------------------------- *
			* Gaps matrix G if requested
			* --------------------------------- *

			if "`gaps'" != "" {

				// gap matrix G
				matrix `S' = (1/2)*(hadamard(`Sstar'*`One', `Sstar'*`One')+hadamard(`One''*`Sstar'', `One''*`Sstar''))
				matrix `G' = J(`K',`K',0)
				forv i = 1/`K' {
					forv j = 1/`K' {

						mat `S'[`i',`j'] = `S'[`i',`j']^.5
						mat `G'[`i',`j'] = 1 / `S'[`i',`j']

					}
				}

				matrix `G' = hadamard( `Mstar' * `One' - `One'' * `Mstar'' , `G' )

				* sampling variances, appendix formula B3

				forv i = 1/1 {
					matrix `Gvar`i'' = J(`K', `K', 0)
				}


				forv g = 1/`K' {
					forv h = 1/`K' {
						if `g' != `h' {

							local deltagh = (`Vstar'[`g',`g'] + `Vstar'[`h',`h'] - 2*`Vstar'[`g',`h'])
							local etagh = 1/(4*`S'[`g',`h']^2)*(`Sstar'[`g',1]^2*`Wstar'[`g',`g']+ `Sstar'[`h',1]^2*`Wstar'[`h',`h'] + 2*`Sstar'[`g',1]*`Sstar'[`h',1]*`Wstar'[`g',`h'])

							matrix `Gvar1'[`g',`h'] = `deltagh'/`S'[`g',`h']^2 * (1 + `G'[`g',`h']^2 * (`etagh'/`deltagh') - (`etagh'/`S'[`g',`h']^2))

						}
					}
				}
			}
		}
	}  // conclude calculations done
	*****

	// returns

	* matrices

	ereturn matrix pk = `P'
	ereturn matrix PI = `PI'

	* return matrix with order of "best" reference groups
	tempname rr
	qui get_ref , refrank(1) catname(`catname') grpid(`grpid') numcats(`numcats')
	mat `rr' = r(refrank)
	ereturn matrix refrank = `rr'

	* return initial values if used
	ereturn matrix initvalsmat = `initvalsmat'
	if "`initvals'" == "" & "`startfrom'" == "" {
		ereturn scalar initvals = 0
	}
	else {
		ereturn scalar initvals = 1
	}

	* return prime and star estimates

	if "`identify'" != "setcuts" {
		matrix `Cprime_se' = vecdiag(cholesky(diag(vecdiag(`Bprime'))))'
	}
	else if "`identify'" == "setcuts" {
		matrix `Cprime_se' = vecdiag(cholesky(diag(vecdiag(`Bprime'))))'
	}

	if `K' > 1 matrix `Mprime_se' = vecdiag(cholesky(diag(vecdiag(`Vprime'))))'

	if "`modtype'" == "hetop" & `K' > 1 {
		cap matrix `Gprime_se' = vecdiag(cholesky(diag(vecdiag(`Omegaprime'))))'
		if _rc != 0 matrix `Gprime_se' = J(`K',1,.)
		cap matrix `Sprime_se' = vecdiag(cholesky(diag(vecdiag(`Wprime'))))'
		if _rc matrix `Sprime_se' = J(`K',1,.)
	}
	else {
		matrix `Gprime_se' = J(`K',1,0)
		matrix `Sprime_se' = J(`K',1,0)
		matrix roweq `Gprime_se' = lnsigma:
		matrix roweq `Sprime_se' = sigma:
		matrix rownames `Sprime_se' = `snames'
		matrix rownames `Gprime_se' = `snames'
	}

	ereturn matrix Braw = `Braw'
	ereturn matrix Zraw = `Zraw'
	ereturn matrix Wraw = `Wraw'
	ereturn matrix Vraw = `Vraw'
	ereturn matrix craw_se = `Craw_se'
	ereturn matrix craw = `Craw'
	ereturn matrix sraw_se = `Sraw_se'
	ereturn matrix sraw = `Sraw'
	ereturn matrix mraw_se = `Mraw_se'
	ereturn matrix mraw = `Mraw'

	if `K' > 1 {

	* return gap estimates
	if "`gaps'" != "" {
		ereturn matrix G = `G'
		ereturn matrix Gvar1 = `Gvar1'	// formula version 1
	}

	ereturn matrix Bprime = `Bprime'
	ereturn matrix Aprime = `Aprime'
	ereturn matrix Zprime = `Zprime'
	ereturn matrix Wprime = `Wprime'
	ereturn matrix Vprime = `Vprime'

	ereturn matrix cprime_se = `Cprime_se'
	ereturn matrix cprime = `Cprime'

	ereturn matrix gprime_se = `Gprime_se'
	ereturn matrix gprime = `Gprime'

	ereturn matrix sprime_se = `Sprime_se'
	ereturn matrix sprime = `Sprime'

	ereturn matrix mprime_se = `Mprime_se'
	ereturn matrix mprime = `Mprime'

	ereturn matrix Bstar = `Bstar'
	ereturn matrix Astar = `Astar'
	ereturn matrix Dstar = `Dstar'
	ereturn matrix Zstar = `Zstar'
	ereturn matrix Wstar = `Wstar'
	ereturn matrix Vstar = `Vstar'

	ereturn matrix cstar_se = `Cstar_se'
	ereturn matrix cstar = `Cstar'

	ereturn matrix sstar_se = `Sstar_se'
	ereturn matrix sstar = `Sstar'

	ereturn matrix mstar_se = `Mstar_se'
	ereturn matrix mstar = `Mstar'

	* scalars and strings
	ereturn scalar mseRC = `mseRC'
	ereturn scalar sseRC = `sseRC'
	ereturn scalar cseRC = `cseRC'
	ereturn scalar varmatsRC = `varmatsRC'

	if `mseRC' != 0 | `sseRC' != 0 | `cseRC' != 0 | `varmatsRC' != 0 {
		noi di in red "Warning: problem with de-referenced SEs. " ///
		"See mseRC, sseRC and/or varmatsRC."
	}

	* estimated icc and SE
	ereturn scalar icchat = `icchat'[1,1]
	ereturn scalar icchatratio = `icchatratio'[1,1]		// alternate formula: sigmab / (sigmab + sigmaw)
	ereturn scalar icchat_var  = `icchatvar'[1,1]

	* sigmaprime estimate
	ereturn scalar sigmaprime = `sigmaprime'[1,1]
	ereturn scalar sigmaw = `sigmaw'[1,1]
	ereturn scalar sigmab = `sigmab'[1,1]
	ereturn scalar varsigprime = `varsigprime'[1,1]

	}

	ereturn scalar numgrps = `K'

	if "`identify'" != "refgroup" {
		ereturn scalar refgrp = 0
	}
	else {
		ereturn scalar refgrp = `ref_group'
	}

	ereturn local identify 	= "`identify'"			// "sums" or "refgroup"
	ereturn local modtype 	= "`modtype'"			// "hetop" or "homop"
	if "`phop'" == "" local phop "."
	ereturn local phop 		= "`phop'"				// phop constraint variable, blank if not used
	ereturn local cmd 		= "hetop"
	if "`modtype'" == "homop" ereturn local csd=`csd'
	if "`modtype'" != "homop" ereturn local csd=.
	ereturn local levelvars = "`vars'"

	// cleanup constraints and macros

	constraint drop `allconstraints'
	macro drop FW_* fixedsd clsd


	*************************************************
	***** save estimates if requested *****

	if "`save'" != "" {

		sort `grpid'

		parse_save_opt `save'

		foreach v in `s(stypes)' {

			do_save `v' , `s(se)' `s(cov)' savename(`s(stubname)')

		}

		keep `grpid' `newvarnames'
		tempfile results
		qui save "`results'"
	}

	************************************************
	* restore and return original data file

	restore

	qui if "`save'" != "" {
		tempvar curorder
		g `curorder' = _n
		merge 1:1 `grpid' using "`results'" , nogen
		sort `curorder'
	}

	ereturn local grpid 	= "`grpid'"

end

*************************************************


*************************************************
* parse PHOP options
cap program drop parse_phop_opt
program define parse_phop_opt , sclass
	version 13.1
	syntax varname [ , mean weight ]
	tokenize `varlist'
	local phop `1'
	sreturn local phopvar `phop'
	if "`mean'" != "" sreturn local phop_mean = 1
end
*************************************************


*************************************************
* getsubmat, 21jun2017
* helper function to get matrix subset
* based somewhat on -matselrc- from nick cox
* programmed here to avoid need for external package

cap program drop getsubmat
program define getsubmat
	version 13.1
	syntax namelist(min=2 max=2) , [ Row(numlist) Col(numlist) ]
	tokenize `namelist'
	local m1 `1'
	local m2 `2'

	tempname a1 a2

	if "`row'" == "" {
		local nr = rowsof(`m1')
		qui numlist "1(1)`nr'"
		local row "`r(numlist)'"
	}
	if "`col'" == "" {
		local nc = colsof(`m1')
		qui numlist "1(1)`nc'"
		local col "`r(numlist)'"
	}

	local i = 1
	foreach num in `row' {
		if `i' == 1 matrix `a1' = `m1'[`num',1...]
		if `i' != 1 matrix `a1' = `a1' \ `m1'[`num',1...]
		local i = `i'+1
	}
	local i = 1
	foreach num in `col' {
		if `i' == 1 matrix `a2' = `a1'[1...,`num']
		if `i' != 1 matrix `a2' = `a2' , `a1'[1...,`num']
		local i = `i'+1
	}

	matrix `m2' = `a2'

end

*************************************************


*************************************************
* parse save syntax

cap program drop parse_save_opt
program define parse_save_opt , sclass
	version 13.1

	syntax [ anything ] [ , STAR PRIME RAW SE COV ]

	tokenize `anything'
	local stubname `1'

	local allvarnames

	foreach v in `star' `prime' `raw' {
		local allvarnames "`allvarnames' m`v'`stubname' s`v'`stubname'"
		if "`se'" != "" local allvarnames "`allvarnames' m`v'`stubname'_se s`v'`stubname'_se"
		if "`cov'" != "" local allvarnames "`allvarnames' z`v'`stubname'_cov"
	}

	sreturn local stubname `stubname'
	sreturn local se `se'
	sreturn local cov `cov'
	sreturn local stypes "`star' `prime' `raw'"
	sreturn local allvarnames "`allvarnames'"

end

*************************************************


*************************************************
* run save command

cap program drop do_save
program define do_save
	version 13.1
	syntax anything [ , SE COV SAVENAME(string) ]

	tokenize `anything'
	local savetype `1'

	foreach v in m s {
		tempname m
		mat `m' = e(`v'`savetype')
		svmat `m'
		rename `m'1 `v'`savetype'`savename'
	}

	if "`se'" != "" {
		foreach v in m s {
			tempname m
			mat `m' = e(`v'`savetype'_se)
			svmat `m'
			rename `m'1 `v'`savetype'`savename'_se
		}
	}

	if "`cov'" != "" {
		tempname z
		matrix `z' = e(Z`savetype')
		matrix `z' = vecdiag(`z')'
		svmat `z'
		rename `z'1 z`savetype'`savename'_cov
	}

end

*************************************************


*************************************************
* determine best reference groups

cap program drop get_ref
program define get_ref , rclass
	version 13.1
	syntax , ///
		[ ///
		REFRANK(integer 1) CATNAME(string) GRPID(string) NUMCATS(integer 0) ///
		]

	* limit to those with at least 3 non-zero cells
	* identify those above median size (of those with at least 3 cells)
	* then sort on proportion distance metric (binned)
	* then sort on sample size

	tempname refmat

	preserve

		qui reshape long `catname', i(`grpid') j(y)
		qui levelsof y , local(ylevels)
		reshape wide

		egen Nk = rowtotal(`catname'*)

		qui su Nk , d
		local Ntot = r(sum)

		egen num0s = anycount(`catname'*) , val(0)
		gen has3 = 0
		replace has3 = 1 if (`numcats' - num0s) > 2

		qui sum Nk if has3 == 1 , detail
		local medN = r(p50)
		gen abvmed = 0
		replace abvmed = 1 if (Nk >= `medN')

		foreach i in `ylevels' {
			qui su `catname'`i'
			gen ptmp`i' = (( `catname'`i' / Nk ) / ( r(sum) / `Ntot' )) - 1
			gen abstmp`i' = abs(ptmp`i')
		}

		egen psum = rowtotal(abstmp*)
		gen psumrd = round(psum, .25)

		gen Nkrank = 1/Nk

		gsort -has3 -abvmed +psumrd +Nkrank +`grpid'
		g tmp_refgrp_rank = _n

		mkmat `grpid' , mat(`refmat')

		qui sum `grpid' if tmp_refgrp_rank == `refrank'
		local useid = r(mean)

	restore

	return matrix refrank = `refmat'
	return scalar refid = `useid'

end
*************************************************


*************************************************
* generate starting values

cap program drop get_initial_values
program define get_initial_values , rclass
	version 13.1
	syntax , ///
	[ ///
		grpid(string) refgroup(string) catname(string) csd(integer 1) ///
		numgrps(integer 0) modtype(string) numcuts(integer 0) ///
		weighttype(string) identify(string) ///
		setcuts(numlist ascending min=1 max=9) ///
		initvals ///
	]

	local K = `numgrps'
	if "`refgroup'" == "n" local refgroup = -1

	tempname bc IV

	mat `IV' = J(1,`=`K'*2',0)

	preserve
		keep `catname'* `grpid'
		qui reshape long `catname' , i(`grpid') j(y)
		qui replace `catname' = 0 if `catname' == .

		* cut scores in standard normal
		qui proportion y [`weighttype' = `catname']
		mat `bc' = e(b)
		mat `bc' = `bc'[1,1..`numcuts']
		forval j = 2/`numcuts' {
			local i = `j' - 1
			mat `bc'[1,`j'] = `bc'[1,`j'] + `bc'[1,`i']
		}
		forv j = 1/`numcuts' {
			if `bc'[1,`j'] == 0 mat `bc'[1,`j'] = 1/10000
			if `bc'[1,`j'] == 1 mat `bc'[1,`j'] = 1-(1/10000)
			mat `bc'[1,`j'] = invnormal(`bc'[1,`j'])
		}

		* means and SDs in standard normal
		tempvar N M S dev
		bys `grpid' : egen `N' = sum(`catname')
		bys `grpid' : egen `M' = sum(y*`catname')
		qui replace `M' = (1/`N')*`M'
		qui generate `dev' = (y-`M')^2
		bys `grpid' : egen `S' = sum(`dev'*`catname')
		qui replace `S' = ((1/`N')*`S')^.5
		qui sum y [`weighttype'=`catname']
		qui replace `M' = (`M'-`r(mean)')/`r(sd)'
		qui replace `S' = `S'/`r(sd)'

		qui collapse `M' `S' `N' , by(`grpid')
		qui replace `M' = 0 if `M' == .
		qui replace `S' = 1 if `S' == 0 | `S' == .

		* re-scale based on identification
		if "`identify'" == "refgroup" {
			qui sum `M' if `grpid' == `refgroup'
			local rm = `r(mean)'
			qui sum `S' if `grpid' == `refgroup'
			local rs = `r(mean)'
			qui replace `M' = (`M'-`rm')/`rs'
			qui replace `S' = `S'/`rs'
			tempname rmatm rmats
			matrix `rmatm' = J(1,`numcuts',`=-1*`rm'')
			matrix `bc' = (`bc'+`rmatm')*(1/`rs')
		}
		else if "`identify'" == "cuts" {
			if `numcuts' == 1 {
				qui sum `S'
				local rs = `r(mean)'/`csd'
				local rm = `bc'[1,1]
				qui replace `M' = `M'/`rs'-`rm' // (`M'-`rm')/`rs'
				matrix `bc'[1,1]=0
			}
			else {
				local c1 = `bc'[1,1]
				local c2 = `bc'[1,2]
				local rslp = 1/(`c2'-`c1')
				local rint = -1*(`rslp'*`c2')
				qui replace `M' = `rslp'*`M'+`rint'
				qui replace `S' = `rslp'*`S'
				forv i = 1/`numcuts' {
					matrix `bc'[1,`i'] = `bc'[1,`i']*`rslp'+`rint'
				}
			}
		}
		else if "`identify'" == "setcuts"{
			if `numcuts' == 1 {
				qui sum `S'
				local rs = `r(mean)'/`csd'
				local rm = `setcuts'-`bc'[1,1]
				qui replace `M' = `M'/`rs'+`rm'
				matrix `bc'[1,1]=`setcuts'
			}
			else {
				tempname X Y B
				qui mata: st_matrix("`Y'", strtoreal(tokens(st_local("setcuts"))))
				matrix `Y' = `Y''
				matrix `X' = J(`numcuts',1,1) , `bc''
				matrix `B' = invsym(`X''*`X')*`X''*`Y'
				qui replace `M' = `B'[2,1]*`M'+`B'[1,1]
				qui replace `S' = `B'[2,1]*`S'
				forv i = 1/`numcuts' {
					local c : word `i' of `setcuts'
					matrix `bc'[1,`i'] = `c'
				}
			}
		}

		qui replace `S' = ln(`S')

		if "`identify'" == "sums" {
			qui sum `M' [`weighttype'=`N']
			local gm = r(mean)
			qui sum `S' [`weighttype'=`N']
			qui replace `S' = `S'-`r(mean)'
			qui replace `M' = (`M'-`gm')/exp(`r(mean)')
		}

		* create initial value matrix and column names
		sort `grpid'
		qui levelsof `grpid' , local(grplevels)
		qui forval k = 1/`K' {
			local j : word `k' of `grplevels'
			local Minitcolumns `Minitcolumns' mean:`j'.`grpid'
			local Sinitcolumns `Sinitcolumns' lnsigma:`j'.`grpid'
		}
		if "`modtype'" == "hetop" local numparms = `numgrps'*2
		if "`modtype'" == "homop" & $fixedsd == 0 {
			local numparms = `numgrps'+1
			local Sinitcolumns "lnsigma:_cons"
		}
		if "`modtype'" == "homop" & $fixedsd == 1 {
			local numparms = `numgrps'
			local Sinitcolumns ""
		}
		matrix `IV' = J(1,`numparms',0)

		* add appropriate initial values
		if "`initvals'" != "" {
			if "`modtype'" == "hetop" {
				qui levelsof `grpid' , local(grplevels)
				qui forval k = 1/`K' {
					local j : word `k' of `grplevels'
					qui sum `M' if `grpid' == `j'
					matrix `IV'[1,`k'] = `r(mean)'
					qui sum `S' if `grpid' == `j'
					matrix `IV'[1,`=`K'+`k''] = `r(mean)'
				}

			}
			else if "`modtype'" == "homop" & $fixedsd == 0 {
				qui levelsof `grpid' , local(grplevels)
				qui forval k = 1/`K' {
					local j : word `k' of `grplevels'
					qui sum `M' if `grpid' == `j'
					matrix `IV'[1,`k'] = `r(mean)'
				}
				qui sum `S' [`weighttype'=`N']
				matrix `IV'[1,`=`numgrps'+1'] = `r(mean)'
			}
			else {
				// "`modtype'" == "homop" & $fixedsd == 1
				qui levelsof `grpid' , local(grplevels)
				qui forval k = 1/`K' {
					local j : word `k' of `grplevels'
					qui sum `M' if `grpid' == `j'
					matrix `IV'[1,`k'] = `r(mean)'
				}
			}
		}

		* label matrix columns
		forval i = 1/`numcuts' {
			local bcnames `bcnames' cut`i':_cons
		}
		matrix colnames `bc' = `bcnames'
		matrix colnames `IV' = `Minitcolumns' `Sinitcolumns'
		matrix `IV' = `IV' , `bc'

	restore

	return matrix initvalsmat = `IV'

end

*************************************************