set more off
* This file includes all the code in the paper's text plus
* additional code, e.g. it shows how to exactly replicate
* other published analyses and then convert the results into oglm
* equivalents.
version 10
********************************************************************
* Example 1: Allison replication

capture program drop palogit
* Allison's original program with a typo corrected.
* Has been renamed palogit from glogit because a
* program named glogit already exists in Stata.
* This is referred to in the text but not shown.
program define palogit
	version 6
	args lnf theta delta
	quietly replace `lnf' = /// 
		$ML_y1*`theta'*(1+`delta') - ///
		ln(1+exp(`theta'*(1+`delta')))
end
use "http://www.indiana.edu/~jslsoc/stata/spex_data/tenure01.dta", clear
* Allison restricted the sample to the first 10 years as an Assistant Prof
keep if year <= 10
* Allison's Table 1 - men only
quietly logit tenure female year yearsq select articles prestige if female==0
quietly estimates store male
* Allison's Table 1 - females only
quietly logit tenure female year yearsq select articles prestige if female==1
quietly estimates store female
* oglm replication of Allison's delta models from his Table 2
quietly oglm tenure year yearsq select articles prestige female, het(female) store(oglm1)
* Compute Allison's delta 
display (1 - exp(.3022305))/ exp(.3022305) 
quietly oglm tenure year yearsq select articles prestige female f_articles, het(female) store(oglm2)
* Compute Allison's delta 
display (1 - exp(.1774193))/ exp(.1774193) 
esttab male female oglm1 oglm2, stats(N ll) mtitle

* Supplemental -- Using Allison's own code...
* Pooled model with delta added
quietly ml model lf palogit (tenure = year yearsq select articles prestige female ) ///
	(delta: female,nocons), maximize
quietly est store delta1
* pooled model with delta and interaction term
quietly ml model lf palogit (tenure = year yearsq select articles prestige female  f_articles) ///
	(delta: female,nocons), maximize
quietly est store delta2
esttab male female delta1 delta2, stats(N ll) mtitle
esttab delta1 oglm1 delta2 oglm2, stats(N ll) mtitle
********************************************************************

* Example 2: Hauser & Andrew Replication
* Make sure data are in current directory, or else change the use command.
* Data are also available at http://www.nd.edu/~rwilliam/oglm/lrpc.dta
use lrpc, clear
* Following is modified from pplogit in Hauser & Andrew's Appendix.
* Double precision is used.  Hauser and Andrew's program
* estimated lambda - 1, this program estimates lambda directly.
* Output has been made more readable.
* In the lrppc model, the equation now labeled "intercepts" would also
* include any other variables whose effects are not multiplied by lambda.
capture program drop lrpc01
program define lrpc01
	tempvar theta
	version 8
	args lnf betas intercepts lambdas
	gen double `theta' = `intercepts' + (`lambdas' * `betas')
	quietly replace `lnf' = ln(exp(`theta')/(1+exp(`theta'))) if $ML_y1==1
	quietly replace `lnf' = ln(1/(1+exp(`theta'))) if $ML_y1==0
end

* The next line constrains lamda 1 to be 1
constraint 1 [lambdas]trans1 = 1

* H & A Table 5, Model 4 - lrpc
ml model lf lrpc01 ///
	(betas: dunc sibsttl9 ln_inc_trunc edhifaom edhimoom broken farm16 south, nocons) ///
	(intercepts: outcome = trans1 trans2 trans3 trans4 trans5 trans6, nocons) ///
	(lambdas: trans1 trans2 trans3 trans4 trans5 trans6, nocons)  ///
	 , maximize constraint(1) nolog
ml display
quietly est store lrpc

* oglm equivalent
oglm outcome dunc sibsttl9 ln_inc_trunc edhifaom edhimoom broken farm16 south trans2 trans3 trans4 trans5 trans6 , ///
	het(trans2 trans3 trans4 trans5 trans6) store(olrpc)


* Convert oglm's choice equations trans1-trans6 into H & A's intercepts
* Do trans1 separately - it is just the negative of oglm's cutpoint
forval i = 1/6 {
	display
	display as yellow "trans`i':"
	if "`i'"=="1" {
		display "oglm's cutpoint coefficient = [#3]_cons = " %6.3f [#3]_cons
		display "Hauser and Andrew's intercept ="
		nlcom -[#3]_cons
	}
	else {
		display "oglm's choice coefficient = [#1]trans`i' = " %6.3f [#1]trans`i'
		display "Hauser and Andrew's intercept ="
		nlcom ([#1]trans`i' - [#3]_cons) * 1/exp([lnsigma]trans`i') 
	}
}

* Convert oglm's lnsigma parameters into H & A's lambdas
* Remember that lambda 1 is fixed at 1
forval i = 2/6 {
	display
	display as yellow "trans`i':"
	display "oglm's sigma = exp([lnsigma]trans`i') = exp(" %6.3f `=[lnsigma]trans`i '' ")= " %9.3f exp([lnsigma]trans`i')
	display "Hauser and Andrew's lambda = 1/sigma ="
	nlcom 1/exp([lnsigma]trans`i')
}
*** H & A Model Table 5 Model 6: An intercept for each transition, proportional effects of 
* socioeconomic variables, interactions of broken, farm, and south with transition.
* This is the second hetero choice model (equivalent to H & A’s LRPPC).
oglm outcome trans2 trans3 trans4 trans5 trans6 broken farm16 south /// 
trans2Xbroken trans2Xfarm16 trans2Xsouth trans3Xbroken trans3Xfarm16 trans3Xsouth /// 
trans4Xbroken trans4Xfarm16 trans4Xsouth trans5Xbroken trans5Xfarm16 trans5Xsouth /// 
trans6Xbroken trans6Xfarm16 trans6Xsouth dunc sibsttl9 /// 
ln_inc_trunc edhifaom edhimoom, het(trans2 trans3 trans4 trans5 trans6) store(m6)


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

* Example 2 - Supplemental: Biochemist data and the LRPC
* This is referred to in the footnote 8 but not shown.

use "http://www.indiana.edu/~jslsoc/stata/spex_data/tenure01.dta", clear
* Allison restricted the sample to the first 10 years as an Assistant Prof
keep if pdasample

capture program drop palogit
* Allison's original program with a typo corrected.
* Has been renamed palogit from glogit because a
* program named glogit already exists in Stata.
program define palogit
	version 6
	args lnf theta delta
	quietly replace `lnf' = /// 
		$ML_y1*`theta'*(1+`delta') - ///
		ln(1+exp(`theta'*(1+`delta')))
end
* Allison's original model, using his program & biochemist data
ml model lf palogit (tenure = year yearsq select articles prestige female ) ///
	(delta: female,nocons), maximize nolog
ml display

capture program drop lrpc02
* Hauser & Andrew's original LRPC program
* Code has been made more efficient and readable,
* but results are the same.  Note that it
* actually estimates and reports 
* lambda - 1 rather than lamba.
program define lrpc02
	tempvar theta
	version 8
	args lnf  intercepts lambdaminus1 betas
	gen double `theta' = `intercepts' + `betas' + (`lambdaminus1' * `betas')
	quietly replace `lnf' = ln(exp(`theta')/(1+exp(`theta'))) if $ML_y1==1
	quietly replace `lnf' = ln(1/(1+exp(`theta'))) if $ML_y1==0
end

* Hauser & Andrews original LRPC parameterization used with Allison's data
ml model lf lrpc02 ///
	(intercepts: tenure = male female, nocons) ///
	(lambdaminus1: female, nocons) ///
	(betas: year yearsq select articles prestige, nocons), max nolog
ml display

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

* Example 3: gologit2 versus oglm

use http://www.indiana.edu/~jslsoc/stata/spex_data/ordwarm2.dta, clear
ologit  warm yr89 male white age ed prst, nolog
gologit2  warm yr89 male white age ed prst, autofit lrf store(gologit2)
oglm  warm yr89 male white age ed prst, het(yr89 male) store(oglm)
lrtest gologit2 oglm, stats force

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

* Example 4: A seemingly trivial change with non-trivial consequences

use "http://www.indiana.edu/~jslsoc/stata/spex_data/tenure01.dta", clear
* Allison restricted the sample to the first 10 years as an Assistant Prof
keep if pdasample

* As before, use female in the equations
quietly oglm  tenure year yearsq select articles prestige female , het(female) store(oglm_f)
* Now use male instead
quietly oglm  tenure year yearsq select articles prestige male , het(male) store(oglm_m)
* Do females only logit model again, using oglm
quietly oglm  tenure year yearsq select articles prestige if female, store(females)
* Do males only logit model again, using oglm
quietly oglm  tenure year yearsq select articles prestige if male, store(males)
esttab oglm_f oglm_m males females, stats(N ll chi2 df_m) mtitle

* Test prestige under the male versus female models
* Female is in the model:
quietly oglm  tenure (year yearsq select articles female), het(female) store(f1)
quietly oglm  tenure (year yearsq select articles female prestige), het(female) store(f2)
lrtest f1 f2, stats
* Male is in the model:
quietly oglm  tenure (year yearsq select articles male), het(male) store(m1)
quietly oglm  tenure (year yearsq select articles male prestige), het(male) store(m2)
lrtest m1 m2, stats

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

* Example 5: Stepwise regressions & Model diagnostics

use "http://www.indiana.edu/~jslsoc/stata/spex_data/tenure01.dta", clear
* Allison restricted the sample to the first 10 years as an Assistant Prof
keep if pdasample

sw, pe(.01) lr: oglm  tenure female year yearsq select  articles prestige, ///
	eq2(female  year yearsq select  articles prestige ) flip store(sw1)
gen articles2 = articles^2
oglm  tenure female year yearsq select articles articles2 prestige, het(articles) store(sw2)
lrtest sw1 sw2, stats

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