// qv: quasi-variances (Firth 2003)
	// verson 1.03 27JUN2014
	// Aspen Chen aspenchensoc@gmail.com

capture program drop qv
program qv, eclass
version 11
if !replay()	{
	syntax varlist(fv) [, Ref(string) LEVel(cilevel) plot saving(passthru)]

	// retrieve command line & df
	if "`e(mi)'"=="mi"	{
		loca cmdline "`e(cmdline_mi)'"
		loca df=`e(df_r_mi)'
	}
	else	{
		loca cmdline "`e(cmdline)'"
		loca df=e(df_m)
	}

	// confirm
	loca fvops = "`s(fvops)'" == "true"
	if `fvops' {
		// variables in last estimation
		if strpos("`cmdline'", "`varlist'")==0	{
			di as error `""`varlist'" not included in previous estimation"'
			exit
		}

		// one variable specified
		if wordcount("`varlist'")>1	{
			di as error "only one variable allowed when factor variables specified"
			exit
		}
		}
	else	{
		// variables in last estimation
		foreach var in `varlist'	{
			if strpos("`cmdline'", "`var'")==0	{
				di as error "`var' not included in previous estimation"
				exit
			}
		}

		// dichotomous variables
		foreach var in `varlist'	{
			qui su `var'
			if	`r(min)'!=0 | `r(max)'!=1 {
				di as error "`var' contains values other than 0 or 1"
				exit
			}
			qui tab `var'
			if `r(r)'!=2	{
				di as error "`var' contains values other than 0 or 1"
				exit
			}
		}
	}

	// parse options
		// ref group
		if "`ref'"==""	{
			loca rf="ref_grp"
		}
		else if wordcount("`ref'")!=1	{
			di as error "only one word allowed in option ref()"
			exit
		}
		else loca rf="`ref'"

		// ci level
		if	"`level'"==""	{
			loca level=95
		}
			loca cilev=(1-(1-(`level'/100))/2)

	// matrices
	tempname b V FLG BE QV _V1 _B1 _V2 _B2 _V1 _B2
		mat `b'=e(b)
		mat `V'=e(V)

	if `fvops'	{
		fvexpand `varlist'
		loca n2=wordcount("`r(varlist)'")
		loca flist "`r(varlist)'"
		loca lv_min=word("`flist'",1)						// first level
		loca lv_max=word("`flist'",-1)						// last level
		loca mx_min=colnumb(`b',"`lv_min'")					// order in matrix
		loca mx_max=colnumb(`b',"`lv_max'")					// order in matrix
		mat `BE'=`b'[1,`mx_min'..`mx_max']					// submatrix for coefficients
		mat `QV'=`V'[`mx_min'..`mx_max',`mx_min'..`mx_max']	// submatrix for var-cov
	}
	else {
		loca n1=wordcount("`varlist'")			// # of all vars
		loca n2=wordcount("`varlist'")+1		// # of vars in contrast list
		loca cols=colsof(`b')
		mat `FLG'=J(1,`cols', 0)
		foreach var in `varlist'	{
			loca i=colnumb(`b',"`var'")
			mat `FLG'[1,`i']=1
		}
		mata: _V1 = select(st_matrix(st_local("V")),(st_matrix(st_local("FLG"))))
		mata: _V1 = select(_V1, (st_matrix(st_local("FLG")))')
		mata: st_matrix("`_V1'",_V1)
		mata: _B1 = select(st_matrix(st_local("b")),(st_matrix(st_local("FLG"))))
		mata: st_matrix("`_B1'",_B1)

		loca elist: colname `b'
		foreach var in `elist'	{
			if regexm("`varlist'","`var'")==0	{
				loca elist=subinword("`elist'","`var'","",.) // same order as last estimate
			}
		}
		loca elist=itrim("`elist'")
		mat colname `_V1'=`elist'
		mat rowname `_V1'=`elist'

		mat `_V2'=J(`n1',`n1',0)
		mat `_B2'=J(1,`n1',0)

		loca i=1
		foreach v1 in `varlist'	{
			loca j=1
				foreach v2 in `varlist'	{
					loca a=colnumb(`_V1',"`v1'")
					loca b=colnumb(`_V1',"`v2'")
					mat `_V2'[`i',`j']=`_V1'[`a',`b']
					loca ++j
				}
			loca c=colnumb(`_V1',"`v1'")
			mat `_B2'[1,`i']=`_B1'[1,`c']
			loca ++i
		}

		// extra row & column
		mat `QV'=J(`n2', `n2', 0)
		mat `QV'[2,2]=`_V2'
		mat `BE'=J(1,`n2',0)
		mat `BE'[1,2]=`_B2'

		// names
		loca flist "`rf' `varlist'"
		matrix colnames `QV'=`flist'
		matrix rownames `QV'=`flist'
	}

	// contrasts
	loca cslist ""
	forvalues x=1/`n2'	{
	loca cslist "`cslist' cs`x'"
	}

	tempvar vij `cslist'
	loca num_cs=`n2'*(`n2'-1)/2

	qui forvalues x=1/`n2'	{
		gen `cs`x''=0 if _n<=`num_cs'
	}

	qui gen `vij'=.	// log(var(Beta_j-Beta_k))
	loca counter=1
	qui forvalues x=1/`n2'	{
		forvalues y=`x'/`n2'	{
			if `x'!=`y'	{
				replace `cs`x''=1 if _n==`counter'
				replace `cs`y''=1 if _n==`counter'
				replace `vij'=ln(`QV'[`x',`x']+`QV'[`y',`y']-2*`QV'[`x',`y']) if _n==`counter'
				loca counter=`counter'+1
			}
		}
	}

	// Estimate q_j+q_k
		capture qui glm `vij' `cs1'-`cs`n2'', f(gau) l(glm_exp) irls nocon	// irls
		if _rc	{
			qui glm `vij' `cs1'-`cs`n2'', f(gau) l(glm_exp) ml nocon search	// ml
		}

	// Set output matrix
	mat QV=e(b)'
	matrix roweq QV=""
	matrix rownames QV=`flist'
	matrix colnames QV=qv
	mat R=J(`n2', 6, 0)
	qui forvalues x=1/`n2'	{
		mat R[`x',1]=`BE'[1,`x']
		mat R[`x',2]=sqrt(`QV'[`x',`x'])
		mat R[`x',3]=sqrt(QV[`x',1])
		mat R[`x',4]=(QV[`x',1])
		mat R[`x',5]=R[`x',1]-invttail(`df',1-`cilev')*R[`x',3]
		mat R[`x',6]=R[`x',1]+invttail(`df',1-`cilev')*R[`x',3]
	}
		mat qvlb=R[1..`n2',5]
		mat qvub=R[1..`n2',6]
		mat rowname qvlb=`flist'
		mat colname qvlb=lb
		mat rowname qvub=`flist'
		mat colname qvub=ub

	// print results
	di _newline(2)
	di "----------------------------------------------------------------------------"
	di _col(2) "Category" _col(20) "Coef." _col(34) "SE" _col(41) "Quasi-SE" _col(55) "lb(QV)" _col(67) "ub(QV)"
	di "----------------------------------------------------------------------------"
	loca counter=1
	qui foreach ctg in `flist'	{
		noi di _col(2) `"`ctg'"' _col(16) %10.4f R[`counter',1] _col(28) %10.4f R[`counter',2] _col(39) %10.4f R[`counter',3] /*
		*/ _col(51) %10.4f R[`counter',5] _col(64) %10.4f R[`counter',6]
		loca ++counter
	}
	di "----------------------------------------------------------------------------"
	di "*lb and ub calculated at the `level'% level"

	// plot
	if "`plot'"!=""		{
		tempvar group beta ub lb
		tempname lb

		qui gen `group'=_n in 1/`n2'
		loca grplbl ""

		qui forval i=1/`n2'	{
			loca lbl=word("`flist'", `i')
			loca grplbl="`grplbl'"+"`i'"+" "+"`lbl'"+" "
		}
		label define `lb' `grplbl'
		label value `group' `lb'

		qui gen double `beta'=.
		qui gen double `ub'=.
		qui gen double `lb'=.
		qui forval i=1/`n2'	{
			replace `beta'=R[`i',1] if _n==`i'
			replace `lb'=R[`i',5] if _n==`i'
			replace `ub'=R[`i',6] if _n==`i'
		}

		graph twoway scatter `group' `beta' || rspike `ub' `lb' `group', hor `saving'		///
		legend(off) ytitle("") ylabel(1(1)`n2', valuelab alter angle(horizontal))
	}

	// return result
	qui `cmdline' 			// reproduce original estimate
	ereturn mat qvub=qvub	// upper bounds
	ereturn mat qvlb=qvlb	// lower bounds
	ereturn mat qv=QV		// quasi variance matrix
}

else {
	if "`e(cmd)'"!="qv" error 301
}
end