*! version 1.2 (25 November 2003)		Joao Pedro Wagner de Azevedo
* version 1.1  (19 November 2003)
* version 1.0 (15 July 2003)
* NJC hacks 15 July 2003
* version 0.7 6/6/03 		

program define crtest, rclass

	version 7.0
	tempvar sample tmp new 
	tempname b valnum rcount nxtrow coef output C testn 
		
	*-> confirm mlogit

	if "`e(cmd)'" != "mlogit" {
		di _n in y "crtest" in r " only works after " in y "mlogit"
		exit 498 
	}
	
	*-> get dependent variable from last mlogit

		local depvar = e(depvar)

	*-> get regressors from last mlogit

		matrix `coef' = e(b)
		local cols = colsof(`coef')
		local regs = `cols' / (e(k_cat) - 1)
		matrix `b' = `coef'[1, 1..`regs']
		local rhs : colnames(`b')
		local rhs : subinstr local rhs "_cons" ""
	
	*-> get outcome labels from last mlogit 
	* This section uses the Stata ado file _pecats from J. Scott Long and 
	* Jeremy Freese
	
		_pecats
		local catnms8 "`r(catnms8)'"
				
	*-> get other variables from last mlogit
	
		local ll = e(ll)
		local numcasts = e(k_cat)
		local df = e(df_m) / (e(k_cat) - 1) 

	
	*-> get weight info from last mlogit
	
		if "`e(wtype)'" != "" {
			local wtis "[`e(wtype)'`e(wexp)']"
			local wexp2 "`e(wexp)'"
		}

	*-> check that estimation sample matches n from regression
	
		quietly { 
			generate `sample' = e(sample)
			
			if "`e(wtype)'" == "" | "`e(wtype)'" == "aweight" /* 
				*/ | "`e(wtype)'" == "pweight" {
				count if `sample'
				scalar `testn' = r(N)
			}
			else if "`e(wtype)'" == "fweight" | /* 
				*/ "`e(wtype)'" == "iweight" {
				local wtexp = substr("`e(wexp)'", 3, .)
				gen `tmp' = (`wtexp') * `sample'
				su `tmp', meanonly
				scalar `testn' = round(r(sum),1)
			}
		}	
			
		if e(N) ~= `testn' {
			di  _n in r "data has been altered since " /*
			*/ in y "mlogit" in r " was estimated"
			exit 459 
		}

	*-> check number of categories
	
	if e(k_cat) == 2 {
		di _n in r "Cramer-Ridder test requires at least 3 dependent categories"
		exit 148 
	}

	*-> output column names
	
	mat `output' = (1, 1, 1, 1) 
	mat colnames `output' = "ln L" "ln Lr" LR "P>chi2"

	*-> reference for test
	
	qui tabulate `depvar' if `sample', matrow(`valnum') matcell(`rcount')

	*-> get the highest scalar of the valnum matrix

	local nrows = rowsof(`valnum')
	
	* grab values 
	local c = `valnum'[1,1] 
	forval i = 2 / `nrows' {
		local c = max(`c', `valnum'[`i',1]) 
	}
               
	*->    cycle through all pairs of outcomes with mlogit without weights
	
if "`e(wtype)'" == "" {
	
	gen `new' = 0 
	
	quietly forval count1 = 1 / `numcasts' {
		tempname ccount1
		local c2 = `count1' + 1
			
		forval count2 = `c2' / `numcasts'  {
			tempname ccount2 adj_`count1'_`count2' 
			tempname LR_`count1'_`count2' 
			
			replace `new' = `depvar'
			replace `new' = `c' + 1 if /*
				*/ `depvar'   == `valnum'[`count1', 1] /*
				*/ | `depvar' == `valnum'[`count2', 1]
			mlogit `new' `rhs' if `sample' 
			local llp = e(ll)
			scalar `ccount1' = `rcount'[`count1', 1]
			scalar `ccount2' = `rcount'[`count2', 1]
			scalar `adj_`count1'_`count2'' = /* 
        */ `ccount1' * ln(`ccount1') + `ccount2' * ln(`ccount2') /*
	*/ -(`ccount1' + `ccount2') * ln(`ccount1'+ `ccount2')
			local llr = `llp' + `adj_`count1'_`count2'' 
			scalar `LR_`count1'_`count2'' =  2 * (`ll' - `llr')
			local chi2 = `LR_`count1'_`count2''
			local pval = chi2tail(`df', `LR_`count1'_`count2'')
			local s1`count1'_`count2' : word `count1' of `catnms8'
                	local s2`count1'_`count2' : word `count2' of `catnms8'
			mat `nxtrow'  = (`ll', `llr', `chi2', `pval')
			mat roweq `nxtrow' = "`s1`count1'_`count2''"
			mat rownames `nxtrow' = "`s2`count1'_`count2''"
			mat `output' = `output' \ `nxtrow' 
		}
	}
	*-> print output 

	di _n in g "**** Cramer-Ridder test for combining outcome categories"
	di _n in g /*
	*/ "Ho: Candidates for pooling have the same regressor coefficients "
	di in g "    apart from the intercept"

	matrix `C' = `output'[2..rowsof(`output'),1..4]
	matrix list `C', format(%9.3f) noheader 
	di _n in g /*
	*/ "degrees of freedom for chi-square distribution:  " in w `df'
}



	*->    cycle through all pairs of outcomes with mlogit with weights

if "`e(wtype)'" != "" {
	
	quietly mat colnames `output' = "ln pseudo-L" "ln Lr" LR "P>chi2"
	
	gen `new' = 0 
	
	if _caller() < 8 {
	
	quietly forval count1 = 1 / `numcasts' {
		tempname ccount1
		local c2 = `count1' + 1
		forval count2 = `c2' / `numcasts'  {
			tempname ccount2 adj_`count1'_`count2' 
			tempname LR_`count1'_`count2' 
			replace `new' = `depvar'
			replace `new' = `c' + 1 if /*
				*/ `depvar'   == `valnum'[`count1', 1] /*
				*/ | `depvar' == `valnum'[`count2', 1]
			mlogit `new' `rhs' `wtis' if `sample' 			
			local llp = e(ll)
			svyprop `depvar' if `depvar'==`count1' [pweight=`wexp2']  
			scalar `ccount1' =r(N_pop)
			svyprop `depvar' if `depvar'==`count2' [pweight=`wexp2'] 
			scalar `ccount2' =r(N_pop)
			scalar `adj_`count1'_`count2'' = /* 
        */ `ccount1' * ln(`ccount1') + `ccount2' * ln(`ccount2') /*
	*/ -(`ccount1' + `ccount2') * ln(`ccount1'+ `ccount2')
			local llr = `llp' + `adj_`count1'_`count2'' 
			scalar `LR_`count1'_`count2'' =  2 * (`ll' - `llr')
			local chi2 = `LR_`count1'_`count2''
			local pval = chi2tail(`df', `LR_`count1'_`count2'')
			local s1`count1'_`count2' : word `count1' of `catnms8'
                	local s2`count1'_`count2' : word `count2' of `catnms8'
			mat `nxtrow'  = (`ll', `llr', `chi2', `pval')
			mat roweq `nxtrow' = "`s1`count1'_`count2''"
			mat rownames `nxtrow' = "`s2`count1'_`count2''"
			mat `output' = `output' \ `nxtrow' 
			}
		}
	}

quietly else {
  	
  	version 8
  	svyset `wtis'
		
		forval count1 = 1 / `numcasts' {
			tempname ccount1
			local c2 = `count1' + 1
			forval count2 = `c2' / `numcasts'  {
				tempname ccount2 adj_`count1'_`count2' 
				tempname LR_`count1'_`count2' 
				replace `new' = `depvar'
				replace `new' = `c' + 1 if /*
					*/ `depvar'   == `valnum'[`count1', 1] /*
					*/ | `depvar' == `valnum'[`count2', 1]
				mlogit `new' `rhs' `wtis' if `sample' 			
				local llp = e(ll)
				svyprop `depvar' if `depvar'==`count1' 
				scalar `ccount1' =r(N_pop)
				svyprop `depvar' if `depvar'==`count2'  
				scalar `ccount2' =r(N_pop)
				scalar `adj_`count1'_`count2'' = /* 
	        */ `ccount1' * ln(`ccount1') + `ccount2' * ln(`ccount2') /*
		*/ -(`ccount1' + `ccount2') * ln(`ccount1'+ `ccount2')
				local llr = `llp' + `adj_`count1'_`count2'' 
				scalar `LR_`count1'_`count2'' =  2 * (`ll' - `llr')
				local chi2 = `LR_`count1'_`count2''
				local pval = chi2tail(`df', `LR_`count1'_`count2'')
				local s1`count1'_`count2' : word `count1' of `catnms8'
	                	local s2`count1'_`count2' : word `count2' of `catnms8'
				mat `nxtrow'  = (`ll', `llr', `chi2', `pval')
				mat roweq `nxtrow' = "`s1`count1'_`count2''"
				mat rownames `nxtrow' = "`s2`count1'_`count2''"
				mat `output' = `output' \ `nxtrow' 
			}
		}
	}
	
	

	*-> print output 
	
        di _n in g "**** Cramer-Ridder test for combining outcome categories"
        di _n in g /*
        */ "Ho: Candidates for pooling have the same regressor coefficients "
        di in g "    apart from the intercept"
	
	matrix `C' = `output'[2..rowsof(`output'),1..4]
	matrix list `C', format(%9.3f) noheader 
	di _n in g /*
	*/ "degrees of freedom for chi-square distribution:  " in w `df'
	di _n in y "crtest uses the " in r e(crittype) in y " for estimations with pweight"
	
}




end