#delim ;
prog def robit, eclass byable(recall) prop(ml_score swml svyb svyj svyr mi);
version 16.0;
*
 Robit wrapper for glm.
 with added estimation results.
*!Authors: Roger Newson, Milena Falcaro
*!Date: 26 October 2022
*;


*
 Check that all necessary Stata packages are present
*;
qui _whichlist robit1 robit2 robit3 robit4 robit5 robit6 robit7 robit8 robit9 robit10,
  pack(xlink);
local incpackages `"`r(incomplete)'"';
if "`incpackages'"!="" {;
  disp as error "robit requires the following SSC packages to be installed:"
    _n as result `"`incpackages'"';
  disp as error "To install these packages, type:";
  foreach IP in `incpackages' {;
    disp as result "ssc install `IP'";
  };
  error 498;
};


if(replay()){;
*
 Beginning of replay section (not indented)
*;

if "`e(cmd)'"!="glm" error 301;
local cmdline `"`e(cmdline)'"';
local cmd2: word 1 of `cmdline';
if `"`cmd2'"'!="robit" error 301;
if _by() error 190;
syntax [, Level(cilevel) noHEADer notable COEFLegend
  ];
_robitplay, level(`level') 
  `table' `header' `coeflegend'
  ;

*
 End of replay section (not indented)
*;
};
else {;
*
 Beginning of non-replay section (not indented)
*;

syntax varlist(numeric fv ts) [if] [in] [fweight pweight iweight aweight],
  DFreedom(numlist integer >=1 <=8 min=1 max=1) [
  Level(cilevel)
  noCONStant
  ASIS
  OFFset(passthru)
  CONSTraints(passthru)
  VCE(passthru)
  COLlinear
  COEFLegend
  DIFficult
  FROM(passthru)
  noHEADer
  noTABLE
  ];
/*
 dfreedom() specifies the robit degrees of freedom.
 Other options are passed to glm.
*/


* Identify sample to use *;
marksample touse;


* Check that y-variable is binary in sample to use *;
local yvar: word 1 of `varlist';
cap assert inlist(`yvar',0,1) if `touse';
if _rc {;
  disp as error "dependent variable must be binary with values 0 or 1";
  error 498;
};
summ `yvar' if `touse', meanonly;
if missing(r(min)) {;
  error 2000;
};
if r(min)==r(max) {;
  disp as error "dependent variable is binary but does not vary"
    _n "all values are " r(min);
  error 498;
};


* Run glm *;
glm `varlist' if `touse' [`weight'`exp'] , level(`level') `eform'
  family(bernoulli) link(robit`dfreedom')
  `constant'
  `asis'
  `offset'
  `constraints'
  `vce'
  `collinear'
  `coeflegend'
  `difficult'
  `from'
  noheader
  notable
  ;


* Change e(cmdline) (but not e(cmd)) to refer to robit *;
ereturn local cmdline `"robit `0'"';


* Add e(depvarsum) *;
local yvar "`e(depvar)'";
if "`e(wtype)'"=="fweight" {;
  qui summ `yvar' [`e(wtype)'`e(wexp)'] if e(sample), meanonly;
};
else {;
  qui summ `yvar' if e(sample), meanonly;
};
ereturn scalar depvarsum=r(sum);


*
 Display results
*;
_robitplay, level(`level') `header' `table' `coeflegend';


*
 End of non-replay section (not indented)
*;
};


end;


#delim ;
prog def _whichlist, rclass;
version 16.0;
/*
 Input a list of which input items
 and optionally a package list
 and output lists of present and absent items
 and complete and incomplete packages.
*/


syntax anything(name=itemlist) [ , Packagelist(namelist) NOIsily ];
*
 packagelist() specifies a list of packages
   for the items to belong to.
 noisily specifies that whichlist will have the output generated by which
   for each item in the item list.
*;
local Nitem: word count `itemlist';

*
 Extend packagelist if required
*;
if "`packagelist'"!="" {;
  local Npackage: word count `packagelist';
  if `Npackage' < `Nitem' {;
    local lastpackage: word `Npackage' of `packagelist';
    forv i1=`=`Npackage'+1'(1)`Nitem' {;
      local packagelist "`packagelist' `lastpackage'";
    };
  };
};


*
 Create present, absent, complete, and incomplete lists
*;
if "`packagelist'"=="" {;
  * Create present and absent lists only *;
  forv i1=1(1)`Nitem' {;
    local itemcur: word `i1' of `itemlist';
    cap `noisily' which `itemcur';
    if _rc local absent `"`absent' `itemcur'"';
    else local present `"`present' `itemcur'"';
  };
};
else {;
  * Create present, absent, complete, and incomplete lists *;
  forv i1=1(1)`Nitem' {;
    local itemcur: word `i1' of `itemlist';
    local packagecur: word `i1' of `packagelist';
    cap `noisily' which `itemcur';
    if _rc {;
      local absent `"`absent' `itemcur'"';
      local incomplete "`incomplete' `packagecur'";
    };
    else {;
      local present `"`present' `itemcur'"';
      local complete "`complete' `packagecur'";
    };
  };
  local incomplete: list uniq incomplete;
  local incomplete: list sort incomplete;
  local complete: list uniq complete;
  local complete: list complete - incomplete;
  local complete: list sort complete;
};
local present: list uniq present;
local present: list sort present;
local absent: list uniq absent;
local absent: list sort absent;


*
 List results
*;
if `"`present'"'!="" {;
  disp as text "Present items:";
  disp as result `"`present'"';
};
if `"`absent'"'!="" {;
  disp as text "Absent items:";
  disp as result `"`absent'"';
};
if "`complete'"!="" {;
  disp as text "Complete packages:";
  disp as result `"`complete'"';
};
if "`incomplete'"!="" {;
  disp as text "Incomplete packages:";
  disp as result `"`incomplete'"';
};


*
 Return results
*;
foreach R in incomplete complete absent present {;
  return local `R' `"``R''"';
};


end;


program define _robitplay;
version 16.0;
*
 Display output
*;

syntax [, Level(cilevel) noHEADer notable COEFLegend ];
/*
 level() specifies a confidence level.
 noheader specifies no header.
 notable specified no table.
*/


* Check that confidence level is within range *;
if((`level'<10.00)|(`level'>99.99)){;
  disp as error "Level must be between 10.00 and 99.99 inclusive";
  exit 198;
};

* Display estimates *;
disp as text "Model: `e(linkt)'" _n;
if "`header'"!="noheader" {;
  disp as text "Number of obs: " as result e(N);
  disp as text "`e(chi2type)' chi2(`=e(df_m)'): " as result e(chi2);
  disp as text "Prob > chi2: " as result e(p);
  if inlist("`e(vce)'","oim","opg") {;
    disp as text "Log likelihood: " as result e(ll);
  };
  else {;
    disp as text "Log pseudolikelihood: " as result e(ll);
  };
  disp;
};
if "`table'"!="notable" ereturn display, level(`level') `coeflegend';

end;