*! version 2.3.0 30aug2016 Richard Williams, rwilliam@nd.edu * 1.0.0 - Initial release * 1.0.1 - BETA experimental log link added * 1.0.2 - BETA hetero/scale options added * 1.0.3 - BETA Bug fixes * 1.0.4 - BETA flip option added * 1.0.5 - BETA Cauchit (Inverse Cauchy) link added * 1.1.0 - 2nd official release * 1.1.1 - Minor release - showeqns option added to ml display; version changed to 9.2 * 1.1.2 - Fixed problem with pweights resulting in incorrect pseudo R^2. * pweights are now rescaled. (Feb 1, 2007) * 1.1.3 - Added svyb and svyj as supported properties (May 6, 2007) * 1.1.4 - Fixed esoteric bug with value labels. Help revised to indicate * better citations & references. * 2.0.0 - Updated for Stata 11. Supports factor variables and the margins * command. Those with Stata 9 or 10 should use oglm9 instead. * 2.0.1 - Predict option tweaked to default to outcome #1 if * no outcome specified for pr or score. This helps the routine to * work with the spost13 commands. * 2.1.0 - Added support for more Stata display options. Added check program * for when svy: is used. svy option dropped; use svy: prefix instead. * 2.1.1 - Fixed bug in or option. Undocumented Collinear option added. * Added by request. Not recommended, use at your own risk. * Help file updated. * 2.1.2 - Fixed a problem of perfect collinearity with rmcoll when using subsamples * Fixed problem with ml options not working correctly * renamed mlopts and diopts to be consistent with other Stata commands * 2.2.0 - ereturns marginsdefault. In Stata 14+, this will cause the margins command * to default to using all the outcomes rather than just the first. * 2.3.0 - Added modifications provided by Ben Shear, Stanford, to specify * initial values via the STARTVALS option program oglm, eclass byable(recall) sortpreserve properties(swml svyr svyb svyj or irr rrr hr eform mi fvaddcons) version 11.2 macro drop oglmby oglmx oglmh Link dv_* // Replay results if that is all that is requested. if replay() { if "`e(cmd)'" != "oglm" { display as error "oglm was not the last estimation command" exit 301 } if _by() { display as error "You cannot use the by command " /// "when replaying oglm results" exit 190 } Replay `0' exit } else { // This not a replay, so estimate the model. // The bysample code identifies the cases currently selected by by. // Cases not in the by sample will later be marked out. tempvar bysample quietly generate byte `bysample' = 1 quietly if _by() replace `bysample' = . if `_byindex'!=_byindex() global oglmby `bysample' Estimate `0' } // Drop global vars macro drop oglmby oglmx oglmh Link dv_* end ************************ program Estimate, eclass version 11.1 syntax [varlist(default=none fv)] [if] [in] /// Standard Stata [pweight fweight iweight] [, /// ROBust CLuster(varname) /// Constraints(string) /// SCore(passthru) /// Level(cilevel) /// display options or rrr irr hr EForm LOG /// LRForce STOre(name) /// special oglm options LInk(string) /// Link fncs, e.g. logit, probit HETero(varlist fv) SCALE(varlist fv) eq2(varlist fv) /// Synonyms, use one only flip /// Reverse location & scale equations HC LS /// Affects equation labeling TRUstme FORCE /// Override some otherwise fatal errors COLlinear /// Don't check for collinearity. Affects both equations STARTVALS(name) /// Name of matrix containing desired initial values; ben shear * /// ml options ] // Syntax checks, special oglm options // nolog is the default unless log is specified if "`log'"=="" local nolog "nolog" // make sure only one eform specified. If so, it will // be local macro eform local eform `or' `irr' `rrr' `hr' `eform' if `:list sizeof eform' > 1 { opts_exclusive "`eform'" } // Only one of hc or ls may be specified. If so, it will // be local macro eqlabel local eqlabel `hc' `ls' if `:list sizeof eqlabel' > 1 { opts_exclusive "`eqlabel'" } // if trustme/force is specified, some otherwise fatal exit commands become warnings only if "`trustme'"!="" | "`force'"!="" { local comment "* " local warning "Warning! " } // Mark the estimation sample. More marking to come. marksample touse // Syntax checks, normal Stata options // display and maximization options _get_diopts diopts options, `options' local diopts `diopts' `eform' level(`level') mlopts mlopts, `options' // Routine checks for cluster, weight if "`cluster'" != "" { local clopt cluster(`cluster') } if "`weight'" != "" { local wgt "[`weight'`exp']" } // Stata commands are inconsistent in accepting colons as an alternative to /. // Comma as an alternative to space sometimes works, sometimes does not. // We therefore change colons to / and commas to spaces, which always work. local constraints: subinstr local constraints ":" "/", all local constraints: subinstr local constraints "," " ", all // Markout other missing values from the estimation sample // Make sure we have cases! markout `touse' `offset' `exposure' markout `touse' `cluster', strok // limit to bysample markout `touse' $oglmby quietly count if `touse' !=0 if r(N) == 0 { display as error "There are no observations left!" display as error "Double-check your sample selection." exit 2000 } // Additional checks on the varlist specfications // Get Y and X variables gettoken y x: varlist // get rid of collinear explanatory variables _rmcoll `x' if `touse', expand `collinear' local x `"`r(varlist)'"' _fv_check_depvar `y' // Link functions // logit is the default link if none has been specified. This // could be changed if someone preferred a different default. if "`link'"=="" local link "logit" // logit link if "`link'"=="logit" | "`link'"=="l" { local link "logit" local link_title "Ordered Logistic Regression" } // probit is also supported else if "`link'"=="probit" | "`link'"=="p" | { local link "probit" local link_title "Ordered Probit Regression" } // cloglog - SPSS calls this nloglog else if "`link'"=="cloglog" | "`link'"=="c" | { local link "cloglog" local link_title "Ordered Cloglog Regression" } // loglog - SPSS calls this cloglog else if "`link'"=="loglog" | "`link'"=="ll" { local link "loglog" local link_title "Ordered Loglog Regression" } // cauchit else if "`link'"=="cauchit" | "`link'"=="ca" | "`link'"=="cau" { local link "cauchit" local link_title "Ordered Cauchit Regression" } // log - This link is experimental, may not work right else if "`link'"=="log" { local link "log" local link_title "Ordered Log Regression" } // non-supported link else { display "" display as error /// "{yellow}`link'{red} is not a legal link function" display "" exit 198 } global Link `link' // heteroskedasticity/ scale/ eq2 option // Can only specify one of hetero, scale, eq2 local het_options = 0 foreach opt in hetero scale eq2 { if "``opt''"!="" local het_options = `het_options' + 1 } if `het_options' > 1 { display as error "hetero, scale & eq2 are synonyms - use only one of them" exit 198 } local hetero `scale' `hetero' `eq2' if "`hetero'" !="" { _rmcoll `hetero' if `touse', expand `collinear' local hetero `r(varlist)' } // Flip the 2 equations if so requested. This is mainly useful // if using the sw or nestreg prefixes. if "`flip'"!="" { local vhetero `x' local vx `hetero' local hetero `vhetero' local x `vx' } // oglmx = 1 if there are X vars, 0 otherwise local Numx: word count `x' global oglmx = "`x'" != "" // oglmh = 1 if there are hetero vars, 0 otherwise local Numh: word count `hetero' global oglmh = "`hetero'" != "" // labeling of equations & model if $oglmh { local link_title Heteroskedastic `link_title' markout `touse' `hetero' } if "`eqlabel'"=="" { local eq1label `y' local eq2label lnsigma } else if "`eqlabel'"=="hc" { local eq1label choice local eq2label variance } else if "`eqlabel'"=="ls" { local eq1label location local eq2label scale } // Generate all the new things we will need for the models // Get Y values from tab tempname Y_Values /* Vector containing the values of the DV */ quietly tab `y' if `touse', matrow(`Y_Values') local M = r(r) if `M' > 20 { di di as error "`warning'`y' has `M' categories - a maximum of 20 is normally allowed" di as error "Make sure you are using an ordinal dependent variable" di `comment'exit 149 } matrix `Y_Values' = `Y_Values'' /* Transpose the matrix */ local numcuts = `M' - 1 local lastcat = `Y_Values'[1, `M'] /* last category will be base category */ // ll program will use these values to determine 1rst, 2nd, 3rd, etc. values of Y macro drop dv_* forval i = 1/`M' { global dv_`i' = `Y_Values'[1, `i'] } // Build equations for full model to be estimated. There are 4 possible // types of models, depending on whether or not Xs and hetero are specified. // First do for models with X's & hetero if $oglmx & $oglmh { local eqs (`eq1label': `y'=`x', noconstant) (`eq2label': `hetero', noconstant) forval i = 1/`numcuts' { local eqs "`eqs' (cut`i':)" } } // Next do Xs but no hetero else if $oglmx & !$oglmh { local eqs (`eq1label': `y'=`x', noconstant) forval i = 1/`numcuts' { local eqs "`eqs' (cut`i':)" } } // Next do hetero only, no Xs else if !$oglmx & $oglmh { local eqs (`eq2label': `y' = `hetero', noconstant) forval i = 1/`numcuts' { local eqs "`eqs' (cut`i':)" } } // Finally, no Xs & no hetero, i.e. cutpoints only else if !$oglmx & !$oglmh { local eqs (cut1: `y'=) forval i = 2/`numcuts' { local eqs "`eqs' (cut`i':)" } } // Start values for each type of link will be in matrix b0 // The Start_Values subroutine will generate the start values tempname b0 // Estimate the constant-only model Start_Values `y', /// wgt(`wgt') touse(`touse') `robust' clopt(`clopt') /// numcuts(`numcuts') b0(`b0') /// link(`link') eqsx(`eqsx') /// constraints(`constraints') /// mlopts(`mlopts') local initopt `s(initopt)' if "`s(LL0)'"!="" local LL0 = `s(LL0)' if "`startvals'" != "" mat `b0' = `startvals' // ben shear added this // Estimate the final model; and add some stats ml model lf oglm_ll `eqs' `wgt' if `touse', /// constraints(`constraints') `robust' `clopt' /// waldtest(-`numcuts') `initopt' title(`link_title') /// collinear missing maximize nocnsnotes /// `score' `nolog' `mlopts' ereturn repost, buildfvinfo ADDCONS // e(k_aux) tells Stata how many cutpoints there are. If you don't // have this, everything gets printed out as separate equations. local numcuts2 `numcuts' ereturn scalar k_aux = `numcuts' ereturn scalar k_exp = `numcuts2' if "`LL0'"!="" ereturn scalar ll_0 = `LL0' // Compute McFadden's Pseudo R^2 if necessary info exists, i.e. it won't // be possible after using svy because svy does not return it. // Also, Stata ml does not rescale pweights, but oglm will if (e(ll) < . & e(ll_0) < .) { if "`weight'" == "pweight" { local wgtvar: word 2 of `=e(wexp)' quietly sum `wgtvar' if e(sample) ereturn scalar ll = e(ll)/ r(mean) } ereturn scalar r2_p = 1 - (e(ll) / e(ll_0) ) } // Stata reports a Wald rather than LR test whenever constraints have been imposed. // The lrforce parameter makes Stata report a LR test instead, provided // crittype is log likelihood. Use this at your own risk local crittype = e(crittype) local chi2type = e(chi2type) if "`lrforce'"!="" & "`crittype'"=="log likelihood" & "`chi2type'"!="LR" { ereturn local chi2type "LR" ereturn scalar chi2 = -2 * (e(ll_0) - e(ll)) ereturn scalar p = chi2tail(e(df_m),e(chi2)) } // Adapted from ologit.ado Stata 14 -- additional margins support for V 14 // I use some different names than ologit because I created equivalents // of them earlier, e.g. M instead of ncat forval i = 1/`M' { local j = `Y_Values'[1,`i'] local mdflt `mdflt' predict(pr outcome(`j')) } ereturn local marginsdefault `"`mdflt'"' // Return assorted values ereturn local xvars `x' ereturn local hetero `hetero' ereturn scalar k_cat = `M' ereturn scalar basecat = `lastcat' ereturn scalar ibasecat = `M' ereturn matrix cat = `Y_Values' ereturn scalar k_eform = 1 ereturn scalar k_eq_model = e(k_eq) - e(k_aux) ereturn local link `link' // Final cleanup macro drop oglmx oglmh Link dv_* ereturn local predict "oglm_p" ereturn local cmd "oglm" ereturn local marginsprop addcons ereturn local marginsok pr xb ereturn local marginsnotok sigma stdp // display the results Replay, `diopts' store(`store') end ************************ program Start_Values, sclass * Gets the start values (i.e. the values of the cutpoints-only model) version 11.1 syntax varname , /// [robust clopt(string) /// wgt(string) touse(varname) b0(name) numcuts(integer -1) /// link(string) eqsx(string) /// constraints(string) /// mlopts(string) ] local numcuts = `numcuts' local M = `numcuts' + 1 local y `varlist' // Matrix b0 will contain the cumulative probabilities for each category // LL0 will contain the log likelihood 0 for non-svy models quietly proportion `y' `wgt' if `touse' , `clopt' nolabel mat `b0' = e(b) forval j = 1/`M' { local LL0 = `LL0' + `e(N)' * `b0'[1, `j'] * ln(`b0'[1,`j']) } mat `b0' = `b0'[1, 1..`numcuts'] // Now we convert the b's to cumulative percentages forval j = 2/`numcuts' { local i = `j' - 1 mat `b0'[1, `j'] = `b0'[1, `j'] + `b0'[1, `i'] } // Convert as needed for the link function used. Remembers signs are // the opposite of gologit2 if "`link'"=="logit" { forval j = 1/`numcuts' { mat `b0'[1,`j'] = logit(`b0'[1,`j']) } } else if "`link'"=="probit" { forval j = 1/`numcuts' { mat `b0'[1,`j'] = invnormal(`b0'[1,`j']) } } else if "`link'"=="cloglog" { forval j = 1/`numcuts' { mat `b0'[1,`j'] = -cloglog(1-`b0'[1,`j']) } } else if "`link'"=="loglog" { forval j = 1/`numcuts' { mat `b0'[1,`j'] = cloglog(`b0'[1,`j']) } } else if "`link'"=="cauchit" { forval j = 1/`numcuts' { mat `b0'[1,`j'] = tan(_pi * (`b0'[1,`j'] - .5)) } } else if "`link'"=="log" { forval j = 1/`numcuts' { mat `b0'[1,`j'] = -ln(1-`b0'[1,`j']) } } forval i = 1/`numcuts' { local columnames `columnames' cut`i':_cons } matrix colnames `b0' = `columnames' // Below, Ben Shear added the " , skip" option which has no impact on default functioning of -oglm- local initopt init(`b0' , skip) lf0(`numcuts' `LL0') sreturn local initopt `initopt' if "`LL0'"!="" sreturn local LL0 `LL0' end ************************ program Replay, eclass version 11.1 syntax [, /// Level(cilevel) /// or irr rrr hr EForm /// STOre(name) /// * ] // make sure only one eform specified. If so, it will // be local macro eform local eform `or' `irr' `rrr' `hr' `eform' if `:list sizeof eform' > 1 { opts_exclusive "`eform'" } // Make sure equation names show if hetero equation is being used if "`e(hetero)'" !="" local showeqns showeqns _get_diopts diopts , `options' local diopts `diopts' `eform' level(`level') `showeqns' ml display , `diopts' // Store results if requested if "`store'"!="" estimates store `store' end