*! Version 1.1.1  October 28, 1999

*  (C) Copyright 1999 Michael Tomz, Gary King, and Langche Zeng.

*  Current version available at http://gking.harvard.edu

*  RElogitQ: Quantities of Interest from Rare Events Logistic Regression



program define relogit_clear, rclass

   version 6.0

set seed 0210

   if "`e(cmd)'" ~= "relogit" {

      di in r _n "You must run relogit before relogitq."

      exit 198

   }

   ck_xc                                          /* is xc_mrt present?     */

   _strip `"`0'"'                                 /* strip fd wrappers      */

   local fd `r(fd)'                               /* save arg in fd(arg)    */

   local 0 `r(rest)'                              /* rest of command line   */

   syntax [, PR RR(str) CHANGEx(str) LISTX MLE    /*

      */ UNBIased SIMS(numlist integer max=1 >0)  /*

      */ BAYES Level(numlist integer max=1 >0 <100)]



   local analyt : word count `bayes'`unbiased'`mle'  /* analytical methods? */

   if `analyt' > 1 {

      di in r _n "Error: Bayes, Unbiased & MLE options are mutually exclusive"

      exit 198

   }

   local p `e(p)'			             /* proportion of 1's in pop */

   gettoken pLO pHI : p			      /* upper and lower bound    */

   if `analyt' > 0 & "`pHI'" ~= "" {

      di in r _n "Bayes, unbiased and MLE options are not available when " _c

      di in r "the fraction of " _n "1's in the population is not known."

      exit 198

   }

   if "`level'" == "" { local level $S_level }/* level for confid interval */

   local lo = (100-`level')/2                 /* lower bound of confid int */

   local hi = `lo' + `level'                  /* upper bound of confid int */

   if "`fd'" ~= "" {

      if "`changex'" == "" {

         di in r _n "If you request a first difference, you must specify " _c

         di in r "the values of" _n "X by using the changex() option"

         exit 198

      }

      ck_rr `fd'

   }

   else {

      if "`changex'" ~= "" {

         di in r _n "If you include the changex() option, you must request "

         di in r "a first-difference by wrapping one of the options in fd()."

         exit 198

      }

      if "`pr'`rr'" == "" { local pr pr }         /* calculate pr by default*/

   }



   * PARSE CHANGEX AND RELATIVE RISK REQUESTS

   _changex `"`changex'"'                         /* parse changex requests */

   local nx `r(nx)'                               /* # sets of changed x's  */

   local i 1

   while `i' <= `nx' {

      local slist`i' `r(slist`i')'                /* list of changes: start */

      local elist`i' `r(elist`i')'                /* list of changes: end   */

      local i = `i' + 1

   }

   _changex `"`rr'"'                              /* parse relrisk requests */

   local nrr `r(nx)'                              /* # sets of rel risks    */

   local i 1

   while `i' <= `nrr' {

      local srr`i' `r(slist`i')'                  /* list of changes: start */

      local err`i' `r(elist`i')'                  /* list of changes: end   */

      local i = `i' + 1

   }



   * RUN SETX TO PRODUCE CHANGEX AND RR SCENARIOS

   local i 1

   while `i' <= `nx' {                            /* FIRST DIFFERENCES      */

      if "`pHI'" == "" {

         tempname xcS`i' xcE`i'                      /* xcStart, xcEnd vectors */

         setx `slist`i'' $mrt_seto keepmrt           /* set starting x's       */

         matrix `xcS`i'' = r(mrt_xc)                 /* save starting x's      */

         setx `elist`i'' $mrt_seto keepmrt           /* set ending x's         */

         matrix `xcE`i'' = r(mrt_xc)                 /* save ending x's        */

      }

      else {

         tempname xcSL`i' xcSH`i' xcEL`i' xcEH`i'

         setx `slist`i'' $mrt_seto keepmrt           /* set starting x's       */

         matrix `xcSL`i'' = r(mrt_xcL)

         matrix `xcSH`i'' = r(mrt_xcH)

         setx `elist`i'' $mrt_seto keepmrt           /* set ending x's         */

         matrix `xcEL`i'' = r(mrt_xcL)

         matrix `xcEH`i'' = r(mrt_xcH)

      }

      local i = `i' + 1

   }

   local i 1

   while `i' <= `nrr' {                           /* RELATIVE RISKS         */

      if "`pHI'" == "" {

         tempname rrS`i' rrE`i'

         setx `srr`i'' $mrt_seto keepmrt

         matrix `rrS`i'' = r(mrt_xc)

         setx `err`i'' $mrt_seto keepmrt

         matrix `rrE`i'' = r(mrt_xc)

      }

      else {

         tempname rrSL`i' rrSH`i' rrEL`i' rrEH`i'

         setx `srr`i'' $mrt_seto keepmrt

         matrix `rrSL`i'' = r(mrt_xcL)

         matrix `rrSH`i'' = r(mrt_xcH)

         setx `err`i'' $mrt_seto keepmrt

         matrix `rrEL`i'' = r(mrt_xcL)

         matrix `rrEH`i'' = r(mrt_xcH)

      }

      local i = `i' + 1

   }



   * SIMULATE PARAMETERS FROM MULTIVARIATE NORMAL

   tempname b V A row

   matrix `b' = e(b)                                /* parameter estimates  */

   matrix `V' = e(V)                                /* var-cov matrix of est*/

   local nptosim = rowsof(`V')		           /* # params to simulate */

   tempvar keepme			                  /* mark the original ds */

   mark `keepme'					    /* to be restored later */

   if "`sims'" == "" { local sims 1000 }            /* 1000 sims by default */

   if `sims' > _N { qui set obs `sims' }            /* expand ds to fit sims*/

   local i 1

   while `i' <= `nptosim' {

      tempvar c`i'

      qui gen `c`i''=invnorm(uniform()) in 1/`sims'

      local cnames `cnames' `c`i''                  /* collect names of vars*/

      local i = `i' + 1

   }

   _chol `V' `nptosim' `A'                          /* Cholesky decomp of v */

   matrix colnames `A' = `cnames'                   /* rename cols to       */

   local i 1

   while `i' <= `nptosim' {

      matrix `row' = `A'[`i',1...]                  /* get i^th row of A    */

      tempvar b`i'                                  /* temporary variable   */

      local bsims `bsims' `b`i''                    /* names of all b's     */

      matrix score `b`i'' = `row'                   /* c(NxK) * row(1xK)'   */

      qui replace `b`i'' = `b`i'' + `b'[1,`i']      /* add mean             */

      local i = `i' + 1

   }

   qui drop `cnames'



   * CREATE LOW AND HIGH INTERCEPT

   if "`pHI'" ~= "" {

      summarize `e(depvar)' if e(sample), meanonly

      local ybar = r(mean)

      tempvar bHI bLO

      qui gen `bLO'=`b`nptosim''-ln(1-`pLO')+ln(`pLO')-ln(`ybar')+ln(1-`ybar')

      qui gen `bHI'=`b`nptosim''-ln(1-`pHI')+ln(`pHI')-ln(`ybar')+ln(1-`ybar')

      local i 1

      while `i' < `nptosim' {

         local bsimsNC `bsimsNC' `b`i''             /* sims w/o constant   */

         local i = `i' + 1

      }

      local bsimsLO `bsimsNC' `bLO'                 /* sims w/ lo constant */

      local bsimsHI `bsimsNC' `bHI'		    /* sims w/ hi constant */

   }



   if `nptosim' > `e(df_m)' { local k ,1 }         /* append constant to xc? */

   if "`listx'" ~= "" { setx }			   /* list chosen x-values   */



   * PROBABILITIES

   tempname xc xcL xcH AnaXb AnaPr C

   tempvar SimXb SimPr SimXbLO SimXbHI SimPrLO SimPrHI

   if "`pr'" ~= "" {

      if "`pHI'" == "" {

         matrix `xc' = (mrt_xc`k')                  /* xchosen (vector)      */

         matrix colnames `xc' = `bsims'             /* rename cols of xc     */

         qui matrix score `SimXb' = `xc'		    /* simulated xbetas      */

         qui gen `SimPr' = 1/(1+exp(-`SimXb'))      /* simulated Pr(Y=1)     */

         if `analyt' == 0 {

            _pctile `SimPr', p(50)			    /* simulated median      */

            return scalar Pr = r(r1)

         }

         else {

            matrix `AnaXb' = (mrt_xc`k')*`b''         /* Analytical xbeta    */

            scalar `AnaPr' = 1/(1+exp(-`AnaXb'[1,1])) /* Analytical Pr(Y=1)  */

            if "`mle'" == "" {

               mat `C' = (`AnaPr'-`AnaPr'^2)*(.5-`AnaPr')     /* correction  */

               mat `C' = `C'*(mrt_xc`k')*`V'*(mrt_xc`k')'     /* correction  */

               if "`bayes'"~="" {scalar `AnaPr'=`AnaPr'+`C'[1,1]}  /* bayes  */

               else { scalar `AnaPr' = `AnaPr' - `C'[1,1] }        /* unbias */

               _recode `AnaPr'                        /* keep wi [0,1] bound */

            }

            return scalar Pr = `AnaPr'	              /* Analytical mean     */

         }

         qui _pctile `SimPr', p(`lo',`hi')            /* calc percentiles    */

         return scalar PrL = r(r1)	              /* lower bound of ci   */

         return scalar PrU = r(r2)		      /* upper bound of ci   */

      }

      else {

         matrix `xcL' = (mrt_xcL`k')                  /* xchosen|pLO (vector)*/

         matrix colnames `xcL' = `bsimsLO'            /* rename cols         */

         qui matrix score `SimXbLO' = `xcL'	      /* simd xbetas|pLO     */

         qui gen `SimPrLO' = 1/(1+exp(-`SimXbLO'))    /* simd Pr(Y=1|pLO)    */

         matrix `xcH' = (mrt_xcH`k')                  /* xchosen|pHI (vector)*/

         matrix colnames `xcH' = `bsimsHI'	      /* rename cols         */

         qui matrix score `SimXbHI' = `xcH'           /* simd xbetas|pHI     */

         qui gen `SimPrHI' = 1/(1+exp(-`SimXbHI'))    /* simd Pr(Y=1|pHI)    */

         qui _pctile `SimPrLO', p(`lo')               /* calc percentiles    */

         return scalar PrL = r(r1)                    /* lower bound of ci   */

		 scalar PrL = r(r1)		//my additions

         qui _pctile `SimPrHI', p(`hi')               /* calc percentiles    */

         return scalar PrU = r(r1)                    /* upper bound ci      */

		 scalar PrU = r(r1)	//mine

         return scalar Pr = (`return(PrU)'+`return(PrL)')/2 /*mean of bounds */

		 scalar temp = r(Pr)

		 scalar Pr = (`return(PrU)'+`return(PrL)')/2

      }

      di _n in g "Probabilities and confidence intervals based on initial X-values"

      di "   Pr(`e(depvar)'==1)  = " %8.5f `return(Pr)'

      di "   `level'% Confidence Interval = " %8.5f `return(PrL)' " to " `return(PrU)'

   scalar PrU= `return(PrU)'

   scalar PrL = `return(PrL)'

   scalar Pr = `return(Pr)'

   }



   * FIRST DIFFERENCES

   tempname xcS xcE AnaXbS AnaXbE AnaPrS AnaPrE xcSL xcSH xcEL xcEH

   tempname xcdiffL xcplus dPrL dPrH xcdiffH

   tempvar SimXbS SimXbE diff SimXbSL SimXbSH SimXbEL SimXbEH diffL diffH

   tempvar ORLO ARpiLO bsLO nonmon nonmonL nonmonU monL monU ARLLO ARULO

   tempvar ORHI ARpiHI bsHI ARLHI ARUHI

   local i 1

   while `i' <= `nx' {

      if "`pHI'" == "" {

         matrix `xcS' = (`xcS`i''`k')            /* starting values for X   */

         matrix `xcE' = (`xcE`i''`k')		 /* ending values for X     */

         matrix colnames `xcS' = `bsims'         /* rename cols of startvals*/

         qui matrix score `SimXbS' = `xcS'	 /* simulated starting Xb's */

         matrix colnames `xcE' = `bsims'         /* rename cols of endvals  */

         qui matrix score `SimXbE' = `xcE'       /* simulated ending Xb's   */

         qui gen `diff' = 1/(1+exp(-`SimXbE')) - 1/(1+exp(-`SimXbS'))

         if `analyt' == 0 {

            _pctile `diff', p(50)

            return scalar dPr_`i' = r(r1)        /* simulated median 1stdiff*/

         }

         else {

            matrix `AnaXbS' = (`xcS`i''`k')*`b''         /* Analyt XbS      */

            scalar `AnaPrS' = 1/(1+exp(-`AnaXbS'[1,1]))  /* Analyt Pr(Y=1|S)*/

            matrix `AnaXbE' = (`xcE`i''`k')*`b''         /* Analyt XbE      */

            scalar `AnaPrE' = 1/(1+exp(-`AnaXbE'[1,1]))  /* Analyt Pr(Y=1|E)*/

            if "`mle'" == "" {

               mat `C' = (`AnaPrS'-`AnaPrS'^2)*(.5-`AnaPrS')  /* correctn|S */

               mat `C' = `C'*(`xcS`i''`k')*`V'*(`xcS`i''`k')' /* correctn|S */

               if "`bayes'"~=""{scalar `AnaPrS'=`AnaPrS'+`C'[1,1]} /* bayes */

               else { scalar `AnaPrS' = `AnaPrS' - `C'[1,1] }      /* unbias*/

               mat `C'= (`AnaPrE'-`AnaPrE'^2)*(.5-`AnaPrE')   /* correctn|E */

               mat `C' = `C'*(`xcE`i''`k')*`V'*(`xcE`i''`k')' /* correctn|E */

               if "`bayes'"~=""{scalar `AnaPrE'=`AnaPrE'+`C'[1,1]} /* bayes */

               else { scalar `AnaPrE' = `AnaPrE' - `C'[1,1] }      /* unbias*/

               _recode `AnaPrS' `AnaPrE'                  /* keep in [0,1]  */

            }

            return scalar dPr_`i' = `AnaPrE' - `AnaPrS'

         }

         qui _pctile `diff', p(`lo',`hi')        /* calculate percentiles   */

         return scalar dPrL_`i' = r(r1)          /* lower bound of ci       */

         return scalar dPrU_`i' = r(r2)          /* upper bound of ci       */

         drop `SimXbS' `SimXbE' `diff'

      }

      else {

         matrix `xcSL' = (`xcSL`i''`k')        /* starting values for X|pLO*/

         matrix `xcEL' = (`xcEL`i''`k')        /* ending values for X|pLO */

         matrix `xcSH' = (`xcSH`i''`k')        /* starting values for X|pHI*/

         matrix `xcEH' = (`xcEH`i''`k')        /* ending values for X|pHI    */

         * Simulate 1st differences under pLO and pHI

         matrix colnames `xcSL' = `bsimsLO'

         qui matrix score `SimXbSL' = `xcSL'             /* sim XB|S,pLO     */

         matrix colnames `xcSH' = `bsimsHI'

         qui matrix score `SimXbSH' = `xcSH'             /* sim XB|S,pHI     */

         matrix colnames `xcEL' = `bsimsLO'

         qui matrix score `SimXbEL' = `xcEL'             /* sim XB|E,pLO     */

         matrix colnames `xcEH' = `bsimsHI'

         qui matrix score `SimXbEH' = `xcEH'             /* sim XB|E,pHI     */

         qui gen `diffL' = 1/(1+exp(-`SimXbEL')) - 1/(1+exp(-`SimXbSL'))

         qui gen `diffH' = 1/(1+exp(-`SimXbEH')) - 1/(1+exp(-`SimXbSH'))



         * CASE 1: pLO

         * Calculate ARpi assuming pLO

         matrix `xcdiffL' = `xcEL' - `xcSL'

         matrix colnames `xcdiffL' = `bsims'

         qui matrix score `ORLO' = `xcdiffL'

         qui replace `ORLO' = exp(`ORLO')                     /* odds ratio */

         qui gen `ARpiLO' = (sqrt(`ORLO')-1)/(sqrt(`ORLO')+1) /* AR_piLO    */

         * Intercept at which AR is maximized, given pLO

         matrix `xcplus' = -.5*(`xcEL' + `xcSL')         /* neg avg xcE, xcS */

         matrix `xcplus'[1,`nptosim'] = 0                /* omit constant    */

         matrix colnames `xcplus' = `bsims'

         qui matrix score `bsLO' = `xcplus'              /* intercept value  */

         * bounds for AR in monotonic and non-monotinic regions

         qui gen `nonmon'=(`bLO'<=`bsLO')&(`bsLO'<=`bHI')  /* nonmonotonic region*/

         qui egen `nonmonL' = rmin(`diffL' `diffH' `ARpiLO') /* nonmonotonic LO*/

         qui egen `nonmonU' = rmax(`diffL' `diffH' `ARpiLO') /* nonmonotonic Hi*/

         qui egen `monL' = rmin(`diffL' `diffH')           /* monotonic LO   */

         qui egen `monU' = rmax(`diffL' `diffH')           /* monotonic HI   */

         qui gen `ARLLO' = (`nonmon')*`nonmonL' + (1-`nonmon')*`monL'

         qui gen `ARULO' = (`nonmon')*`nonmonU' + (1-`nonmon')*`monU'

         * bounds

         qui _pctile `ARLLO', p(`lo')                  /* calc percentiles    */

         scalar `dPrL' = r(r1)                         /* lower bound of ci   */

         qui _pctile `ARULO', p(`hi')                  /* calc percentiles    */

         scalar `dPrH' = r(r1)                         /* upper bound ci      */

         qui drop `ORLO' `ARpiLO' `bsLO' `nonmon' `nonmonL' `nonmonU' /*

            */ `monL' `monU' `ARLLO' `ARULO'



         * CASE2: piHI

         * Calculate ARpi assuming pHI

         matrix `xcdiffH' = `xcEH' - `xcSH'

         matrix colnames `xcdiffH' = `bsims'

         qui matrix score `ORHI' = `xcdiffH'

         qui replace `ORHI' = exp(`ORHI')                     /* odds ratio */

         qui gen `ARpiHI' = (sqrt(`ORHI')-1)/(sqrt(`ORHI')+1) /* AR_piHI    */

         * Intercept at which AR is maximized, given pHI

         matrix `xcplus' = -.5*(`xcEH' + `xcSH')         /* neg avg xcE, xcS */

         matrix `xcplus'[1,`nptosim'] = 0                /* omit constant    */

         matrix colnames `xcplus' = `bsims'

         qui matrix score `bsHI' = `xcplus'              /* intercept value  */

         * bounds for AR in monotonic and non-monotinic regions

         qui gen `nonmon'=(`bLO'<=`bsHI')&(`bsHI'<=`bHI')  /* nonmonotonic region*/

         qui egen `nonmonL' = rmin(`diffL' `diffH' `ARpiHI') /* nonmonotonic LO*/

         qui egen `nonmonU' = rmax(`diffL' `diffH' `ARpiHI') /* nonmonotonic Hi*/

         qui egen `monL' = rmin(`diffL' `diffH')           /* monotonic LO   */

         qui egen `monU' = rmax(`diffL' `diffH')           /* monotonic HI   */

         qui gen `ARLHI' = (`nonmon')*`nonmonL' + (1-`nonmon')*`monL'

         qui gen `ARUHI' = (`nonmon')*`nonmonU' + (1-`nonmon')*`monU'

         * bounds

         qui _pctile `ARLHI', p(`lo')                  /* calc percentiles    */

         if r(r1) < `dPrL' { scalar `dPrL' = r(r1) }

         qui _pctile `ARUHI', p(`hi')                  /* calc percentiles    */

         if r(r1) > `dPrH' { scalar `dPrH' = r(r1) }

         qui drop `ORHI' `ARpiHI' `bsHI' `nonmon' `nonmonL' `nonmonU' /*

            */ `monL' `monU' `ARLHI' `ARUHI' `SimXbSL' `SimXbSH' /*

            */ `SimXbEL' `SimXbEH' `diffL' `diffH'



         * FINAL BOUNDS

         return scalar dPrL_`i' = `dPrL'              /* lower bound of ci   */

         return scalar dPrU_`i' = `dPrH'               /* upper bound ci      */

         return scalar dPr_`i' = (`return(dPrL_`i')'+`return(dPrU_`i')')/2

      }

      gettoken change changex : changex, parse("&")

      di _n in g "First difference `i': `change'"

      di "   dPr(`e(depvar)'==1)  = " %8.5f `return(dPr_`i')'

      di "   `level'% Confidence Interval = " %8.5f `return(dPrL_`i')' " to " `return(dPrU_`i')'

      gettoken ampers changex : changex, parse("&")

      local i = `i' + 1

   }



   * RELATIVE RISKS

   tempvar SimXbS SimXbE RR

   local i 1

   while `i' <= `nrr' {

      if "`pHI'" == "" {

         matrix `xcS' = (`rrS`i''`k')

         matrix `xcE' = (`rrE`i''`k')

         matrix colnames `xcS' = `bsims'         /* rename cols of startvals*/

         qui matrix score `SimXbS' = `xcS'	 /* simulated starting Xb's */

         matrix colnames `xcE' = `bsims'         /* rename cols of endvals  */

         qui matrix score `SimXbE' = `xcE'       /* simulated ending Xb's   */

         qui gen `RR' = (1/(1+exp(-`SimXbE'))) / (1/(1+exp(-`SimXbS')))

         if `analyt' == 0 {

	     _pctile `RR', p(50)

            return scalar rr_`i' = r(r1)         /* simulated median relrisk*/

         }

         else {

            matrix `AnaXbS' = (`rrS`i''`k')*`b''         /* Analyt XbS      */

            scalar `AnaPrS' = 1/(1+exp(-`AnaXbS'[1,1]))  /* Analyt Pr(Y=1|S)*/

            matrix `AnaXbE' = (`rrE`i''`k')*`b''         /* Analyt XbE      */

            scalar `AnaPrE' = 1/(1+exp(-`AnaXbE'[1,1]))  /* Analyt Pr(Y=1|E)*/

            if "`mle'" == "" {

               mat `C' = (`AnaPrS'-`AnaPrS'^2)*(.5-`AnaPrS')  /* correctn|S */

               mat `C' = `C'*(`rrS`i''`k')*`V'*(`rrS`i''`k')' /* correctn|S */

               if "`bayes'"~=""{scalar `AnaPrS'=`AnaPrS'+`C'[1,1]} /* bayes */

               else { scalar `AnaPrS' = `AnaPrS' - `C'[1,1] }      /* unbias*/

               mat `C'= (`AnaPrE'-`AnaPrE'^2)*(.5-`AnaPrE')   /* correctn|E */

               mat `C' = `C'*(`rrE`i''`k')*`V'*(`rrE`i''`k')' /* correctn|E */

               if "`bayes'"~=""{scalar `AnaPrE'=`AnaPrE'+`C'[1,1]} /* bayes */

               else { scalar `AnaPrE' = `AnaPrE' - `C'[1,1] }      /* unbias*/

               _recode `AnaPrS' `AnaPrE'                  /* keep in [0,1]  */

            }

            return scalar rr_`i' = `AnaPrE' / `AnaPrS'

         }

         qui _pctile `RR', p(`lo',`hi')          /* calculate percentiles   */

         return scalar rrL_`i' = r(r1)           /* lower bound of ci       */

         return scalar rrU_`i' = r(r2)           /* upper bound of ci       */

         drop `SimXbS' `SimXbE' `RR'

      }

      else {

         tempvar SimXbSL SimXbSH SimXbEL SimXbEH rrL rrH RRL RRH

         matrix `xcSL' = (`rrSL`i''`k')

         matrix `xcSH' = (`rrSH`i''`k')

         matrix `xcEL' = (`rrEL`i''`k')

         matrix `xcEH' = (`rrEH`i''`k')

         matrix colnames `xcSL' = `bsimsLO'

         qui matrix score `SimXbSL' = `xcSL'             /* sim XB|S,pLO     */

         matrix colnames `xcSH' = `bsimsHI'

         qui matrix score `SimXbSH' = `xcSH'             /* sim XB|S,pHI     */

         matrix colnames `xcEL' = `bsimsLO'

         qui matrix score `SimXbEL' = `xcEL'             /* sim XB|E,pLO     */

         matrix colnames `xcEH' = `bsimsHI'

         qui matrix score `SimXbEH' = `xcEH'             /* sim XB|E,pHI     */

         * simd rel risks pLO and pHI

         qui gen `rrL' = (1+exp(-`SimXbSL')) / (1+exp(-`SimXbEL'))

         qui gen `rrH' = (1+exp(-`SimXbSH')) / (1+exp(-`SimXbEH'))

         qui egen `RRL' = rmin(`rrL' `rrH')

         qui egen `RRH' = rmax(`rrL' `rrH')

         * final bounds

         qui _pctile `RRL', p(`lo')                   /* calc percentiles    */

         return scalar rrL_`i' = r(r1)                /* lower bound of ci   */

         qui _pctile `RRH', p(`hi')                   /* calc percentiles    */

  	   return scalar rrU_`i' = r(r1)                /* upper bound ci      */

         return scalar rr_`i' = (`return(rrL_`i')'+`return(rrU_`i')')/2

         qui drop `SimXbSL' `SimXbSH' `SimXbEL' `SimXbEH' `rrL' `rrH' /*

            */ `RRL' `RRH'

      }

      gettoken change rr : rr, parse("&")

      di _n in g "Relative Risk `i': `change'"



      di "   rr(`e(depvar)'==1)  = " %8.5f `return(rr_`i')'

      di "   `level'% Confidence Interval = " %8.5f `return(rrL_`i')' " to " `return(rrU_`i')'

      gettoken ampers changex : changex, parse("&")

      local i = `i' + 1

   }



   qui keep if `keepme' == 1



end



****************************** SUB-PROCEDURES ********************************



program define ck_xc

   * Checks for presence of xc matrix, simulated parameters, depvar

   version 6.0

   tokenize "`e(p)'"

   if "`2'" == "" {

      capture mat list mrt_xc

      if _rc == 0  {

         local xccols : colnames(mrt_xc)

         if "`xccols'" == "`e(rhsvars)'" { local xcfound yes }

      }

   }

   else {

      capture mat list mrt_xcL                         /* does xcL mat exist? */

      local rc = _rc

      capture mat list mrt_xcH			       /* does xcH mat exist? */

      local rc = `rc' + _rc

      if `rc' == 0 {                                   /* if xc does exist... */

         local xccolsL : colnames(mrt_xcL)             /* names of vars in xcL*/

         local xccolsH : colnames(mrt_xcH)             /* names of vars in xcH*/

         if "`xccolsL'"=="`e(rhsvars)'"  &             /* if names match ivars

            */ & "`xccolsH'"=="`e(rhsvars)'"           /*

            */ { local xcfound yes }                   /* then xc matrix is OK*/

      }

   }

   if "`xcfound'" ~= "yes" {

      di in r _n "X-values missing or invalid.  Please rerun -setx-."

      exit 198

   }



end



program define _strip, rclass

   * Strips fd and changex from command line

   version 6.0

   args cmdline

   tokenize `"`cmdline'"', parse("@")

   if "`1'" == "@" | "`2'" == "@" {

      di in r _n "The -relogitq- command cannot contain the @ character."

      exit 198

   }

   _parsep fd `cmdline'                   /* strip-off changex(arg)   */

   return local fd `r(Arg)'

   return local rest `r(RestOpt)'

end



program define _parsep, rclass

   * Parsing with embedded parentheses

   * Based on Jeroen Weesie, "Parsing Options with Embedded Parentheses,"

   * Stata Technical Bulletin 40 (Nov 1997), insert ip22.

   * Limitation: input can't contain any @'s

   version 6.0

   if "`*'" == "" {

      return local Arg        /* pass empty argument  */

      return local RestOpt    /* pass empty remainder */

      exit

   }

   local optname `1'                 /* pluck-off option name */

   mac shift

   while "`1'" ~= "" {               /* separate all tokens with @ */

      local optamp `optamp'`1'@

      mac shift

   }

   local H 0

   local Mode0 None

   local ProcArg 0

   tokenize "`optamp'", parse("@()[]")

   while "`1'" ~= "" {

      if "`1'" == "`optname'" & `H' == 0 {

         if "`2'" == "(" {

            local H = `H' + 1

            local Mode`H' p

            local ProcArg 1

            mac shift

         }

      }

      else {

         if "`1'" == "(" {

            local H = `H' + 1

            local Mode`H' p

         }

         else if "`1'" == "[" {

            local H = `H' + 1

            local Mode`H' b

         }

         else if "`1'" == ")" {

            if "`Mode`H''" ~= "p" {

               di in r "too many ')' or ']'"

               exit 132

            }

            local H = `H' - 1

         }

         else if "`1'" == "]" {

            if "`Mode`H''" ~= "b" {

               di in r "too many ')' or ']'"

               exit 132

            }

            local H = `H' - 1

         }

         else if "`1'" == "@" {

            local 1 " "

         }

         if `ProcArg' == 1 {

            if `H' > 0 { local Arg "`Arg'`1'" }

         }

         else { local RestOpt "`RestOpt'`1'" }

         if `H' == 0 {

            local ProcArg 0

            local Arg "`Arg' "

         }

      }

      mac shift

   }

   return local Arg `Arg'

   return local RestOpt `RestOpt'

end



program define ck_fd

   tokenize "`changex'", parse(",")

   if "`1'" == "," | "`2'" == "," {

      di in r _n "The changex option cannot contain a comma."

      exit 198

   }

end



program define _changex, rclass

   * parses the changex commands

   version 6.0

   args changex

   local nrhsvar : word count `e(rhsvars)'                /* # rhs vars  */

   local nx 0                                             /* # sets of x */

   while "`changex'" ~= "" {

      local nx = `nx' + 1                                 /* next set    */

      gettoken changes changex : changex, parse("&")      /* changes in  */

      while "`changes'" ~= "" {                           /*   this set  */

         gettoken vars changes : changes, match(paren)    /* variable(s) */

         gettoken fun1 changes : changes, match(paren1)   /* function 1  */

         gettoken fun2 changes : changes, match(paren2)   /* function 2  */

         while "`vars'" ~= "" {

            gettoken var vars : vars                      /* get varname */

            if "`paren1'" == "" {local slist`nx' `slist`nx'' `var' `fun1'}

            else {local slist`nx' `slist`nx'' `var' (`fun1')}

            if "`paren2'" == "" {local elist`nx' `elist`nx'' `var' `fun2'}

            else {local elist`nx' `elist`nx'' `var' (`fun2')}

         }

      }

      return local slist`nx' `slist`nx''    /* return starting list #nx */

      return local elist`nx' `elist`nx''    /* return ending list #nx   */

      gettoken amper changex : changex, parse("&")       /* strip-off & */

   }

   return local nx `nx'

end



program define ck_pr, rclass

   version 6.0

   args pr prval

   tempname catdi allcats

   if "`pr'" == "no" {                           /* pr was not specified */

      matrix `catdi' = J(1,2,0)

      while "`prval'" ~= "" {

         gettoken val prval : prval

         if `val' == 0 { matrix `catdi'[1,1] = 1 }

         else if `val' == 1 { matrix `catdi'[1,2] = 1 }

         else {

            di in r _n "Error: `val' is not a valid outcome for `e(depvar)'"

            exit 198

         }

      }

   }

   else {                                           /* pr was specified */

      if "`prval'" ~= "" {

         di in g _n "Note: You specified both the pr and the prval() " /*

           */ "options." _n "pr takes precedence over prval(), so -relogitq-"/*

           */ " will display " _n "Pr(Y=0) and Pr(Y=1), rather than the " /*

           */ "values listed in prval()."

      }

      matrix `catdi' = J(1,2,1)

   }

   return matrix catdi `catdi'

end



program define _recode

   version 6.0

   while "`0'" ~= "" {

      gettoken p 0 : 0

      if `p' < 0 { scalar `p' = 0 }

      else if `p' > 1 { scalar `p' = 1 }

   }

end





*! version 1.3  April 24, 1999  Michael Tomz

* Cholesky decomposition of an arbitrary matrix

* Input: V, the original matrix

*        k, the dimension of the matrix

*        A, the name of a new matrix to be created

* Output: A, the lower-triangle cholesky of V

program define _chol

   version 6.0

   args V k A

   capt matrix `A' = cholesky(`V')             /* square root of VC matrix  */

   if _rc == 506 {                             /* If VC ~ pos definite, ... */

      tempname eye transf varianc vcd tmp

      mat `eye' = I(`k')                       /* identity matrix           */

      mat `transf' = J(1,`k',0)                /* initialize transf matrix  */

      local i 1

      while `i' <= `k' {

         scalar `varianc' = `V'[`i',`i']       /* variance of parameter `i' */

         if `varianc' ~= 0 {                   /* if has variance, add row  */

            mat `tmp' = `eye'[`i',1...]        /* of `eye' to transf matrix */

            mat `transf' = `transf' \ `tmp'

         }

         local i = `i' + 1

      }

      mat `transf' = `transf'[2...,1...]       /* drop 1st row of transf mat*/

      mat `vcd' = `transf'*`V'*`transf''       /* decomposed VC (no 0 rows) */

      mat `A' = cholesky(`vcd')                /* square root of decomp VC  */

      mat `A' = `transf''*`A'*`transf'         /* rebuild full sq root mat  */

   }

   else if _rc { matrix `A' = cholesky(`V') }  /* redisplay error message   */

end



program define ck_rr

   version 6.0

   args fd

   local 0 ,`fd'

   syntax [, PR RR(string) ]

   if "`rr'" ~= "" {

      di in r _n "Error: relogitq does not allow first differences of " /*

         */ "relative risks."

      exit 198

   }

end