/*
*! version 3.17  16may2022
	handles new metan which returns results on log RR/OR scale
	with thanks to David Fisher 15feb2022
version 3.16  Ian White 02oct2018
version 3.16  Ian White 02oct2018
	corrected error message if called with no observations
	check if logimor() etc are variables, and if so return an error
	[previously their first value was taken. To do: handle them as variables]
version 3.15  3aug2009 - checks for negative imors; lists cases with missing estimate or se
version 3.14  10mar2009 - fixes bug with title("...")
version 3.13  16sep2008 - updated to Stata 10.1: method(MC) is 4-5 times faster
                          as it uses drawnorm and rbeta() instead of my random.ado
version 3.12  27feb2008 - checks metan version before starting. Included in final SJ submission 22aug2008.
version 3.11  19dec2007 - options passed to metan are printed
                          bug fix: method(MC) was wrong
version 3.10  21nov2007 - bug fixes: method(MC), corrlogimor(), nip() work
                          changes: simple dots, drop timer option
                          re-organised output
version 3.9  29oct2007 - bug fix: includes dotter.ado
version 3.8  11oct2007 - state when metan starts (so user can see source of error); small bug in gamblehollis
version 3.7   1oct2007 - logimor defaults to 0 if sdlogimor specified
version 3.6   7sep2007 - new gamblehollis option; report measure used; method(mc) defaults to RR not all 3; if meta not done, explain why; warning if w1 used.
version 3.5   5sep2007 - report weighting scheme; imputes for studies with all reason counts 0 (but exits if all studies are like this)
version 3.4    4sep2007 - new syntax allowing ica0 or ica0(var var) etc.; icab/icaw the right way round; big rewrite of Julian's part
version 3.3   24apr2007 - (Julian) changed adjustments for empty cells
version 3.2   24apr2007 - selects by `measure' not by `rr', `or', `rd'
version 3.1   24apr2007 - graphs fixed
version 3.0    9feb2007 - combined Ian's and Julian's programs
*/

program define metamiss
version 10.1

* Note: this file includes programs expectn, exit498

********************************** CHECK METAN VERSION **********************************

if _N==0 error 2000

tempvar one
gen byte `one'=1
capture {
	metan `one' `one' `one' `one', nointeger nograph
	metan `one' `one', nograph
}
if _rc {
    di as error "Your version of -metan- appears to be too old to be compatible with -metamiss-."
    di as error "Please upgrade to a more recent version (e.g. from http://fmwww.bc.edu/RePEc/bocode/m)."
    di as error "You can do this simply by running " as input "{stata net install metan, replace}"
    exit 498
}
local oldmetan = (`"`r(metan_version)'"'==`""')		// DF 15feb2022: if r(metan_version) not returned, must be v3.04 or older
drop `one'

************************************ COMMON PARSING ************************************

* some disallowed syntaxes
syntax varlist(min=6 max=6) [if] [in], [aca(string) icab(string) icaw(string) *]
if "`aca'"!="" exit498 aca(`aca') not allowed
if "`icab'"!="" exit498 icab(`icab') not allowed
if "`icaw'"!="" exit498 icaw(`icaw') not allowed

* DETERMINE IF OPTIONS ICA0 ETC. HAVE ARGUMENTS (SYNTAX 2) OR NOT (SYNTAX 1)
syntax varlist(min=6 max=6) [if] [in], [aca ica0(string) ica1(string) icapc(string) icape(string) icap(string) icaimor(string) *]
if "`ica0'`ica1'`icapc'`icape'`icap'`icaimor'"!="" {
    if "`aca'"=="aca" exit498 Can't have both ACA and ICA.
    local syntax 2
    local JHoptions ica0(string) ica1(string) icapc(string) icape(string) icap(string) icaimor(string)
}
syntax varlist(min=6 max=6) [if] [in], [ica0 ica1 icapc icape icap icaimor aca icab icaw *]
local words = wordcount("`ica0' `ica1' `icap' `icape' `icapc' `icaimor' `aca' `icab' `icaw'")
if `words' > 1 {
    if "`aca'"=="aca" exit498 Can't have both ACA and ICA.
    else exit498 State only one ICA analysis, or several ICA analyses each with arguments
}
if `words' == 1 {
    if "`syntax'"=="2" exit498 State only one ICA analysis, or several ICA analyses each with arguments
    local syntax 1
    local JHoptions aca ica0 ica1 icapc icape icap icaimor
}

* MAIN PARSING
syntax varlist(min=6 max=6) [if] [in], [                  ///
///Common options:
    or rr rd                                              /// meta-analysis options
    id(varname) log nograph                               /// reporting options
    logimor(string) imor(string)                          /// nonresponse options
///Syntax 3 options:
    SDlogimor(string) CORRlogimor(string)                 /// nonresponse options
    method(string) reps(int -1) nip(int 10)               /// method options
    MISSprior(string) RESPprior(string) noDOTs noreplace  /// full-bayes options
    debug                                                 /// debugging options
///Syntax 1 & 2 options:
    `JHoptions' icab icaw                                 /// main nonresponse options
    w1 w2 w3 w4                                           /// std error options
    listp listx listnum listall                           /// debugging options
///Syntax 4 option:
    GAMBLEhollis                                          ///
*]

* FINISH IDENTIFYING SYNTAX
if "`sdlogimor'"!="" {
    if "`syntax'"!="" exit498 Can't have sdlogimor() with ACA or ICA.
    local syntax 3
    if "`imor'"=="" & "`logimor'"=="" local logimor 0
}

if "`gamblehollis'"!="" {
    if "`syntax'"=="1" | "`syntax'"=="2" exit498 Can't combine gamblehollis with ACA or ICA.
    else if "`syntax'"=="3" exit498 Can't combine gamblehollis with Bayesian analysis.
    else local syntax 4
}

* default to ICA-IMOR if logimor or imor is specified
if "`syntax'"=="" {
    * No method specified
    if "`logimor'`imor'"=="" exit498 Please specify an analysis method.
    local icaimor icaimor
    local syntax 1
}

if `syntax'<=2 & "`sdlogimor'"!="" exit498 Can't have sdlogimor() with ACA or ICA.
if `syntax'<=2 & "`corrlogimor'"!="" exit498 Can't have corrlogimor() with ACA or ICA.


di as text "*******************************************************************"
di as text "******** METAMISS: meta-analysis allowing for missing data ********"
if "`aca'"=="aca" di as text "********               Available cases analysis            ********"
else if `syntax'==1 di as text "********                  Simple imputation                ********"
else if `syntax'==2 di as text "********             Imputation using reasons              ********"
else if `syntax'==3 di as text "********           Bayesian analysis using priors          ********"
else if `syntax'==4 di as text "********               Gamble-Hollis analysis              ********"
else exit498 Syntax not identified
di as text "*******************************************************************"

* END OF MAIN PARSING AND IDENTIFYING SYNTAX: NOW SORT OUT SOME OPTIONS

if `syntax'<=2 & "`icaimor'"=="" {
    if "`imor'"!="" di as error "Warning: icaimor not specified, imor(`imor') ignored."
    if "`logimor'"!="" di as error "Warning: icaimor not specified, logimor(`logimor') ignored."
}

* SORT OUT METHOD
if `syntax'!=3 {
    if "`method'"!="" {
       di as error "Warning: method(`method') ignored."
       local method
    }
    if `nip'!=10 di as error  "Warning: nip(`nip') ignored."
    if `reps'!=-1 di as error  "Warning: reps(`reps') ignored."
    if "`missprior'"!="" di as error  "Warning: missprior(`missprior') ignored."
    if "`respprior'"!="" di as error  "Warning: respprior(`respprior') ignored."
}
local method = lower("`method'")

* SORT OUT MEASURES
if "`method'"=="mc" {
    if "`rd'"=="rd" local measure RD
    if "`rr'"=="rr" local measure `measure' logRR
    if "`or'"=="or" local measure `measure' logOR
    if "`measure'"=="" local measure logRR
}
else {
   local measure = upper("`or'`rr'`rd'")
   if length("`measure'")>2 {
      di as error "Please specify only one of or, rr, rd."
      exit 498
   }
   if length("`measure'")==0 local measure RR
   if "`rd'"=="rd" & "`log'"=="log" {
      di as error "Warning: log option ignored with RD."
      local log
   }
}
local nmeasures = wordcount("`measure'")
di as text "Measure: `measure'" _c
if "`log'"=="log" di as text " (log scale)" _c
di "."

* parse imor() to imorE, imorC, etc.
local parselist logimor imor
if `syntax'==3 local parselist `parselist' sdlogimor
foreach thing in `parselist' {
   if "``thing''"!="" {
       tokenize "``thing''"
       local `thing'E `1'
       if "`2'"=="" local `thing'C `1'
       else local `thing'C `2'
       if "`3'"!="" {
           di in red "Too many arguments in `thing'"
           exit 498
       }
   }
}
if `syntax'==2 {
   foreach thing in ica0 ica1 icap icapc icape icaimor {
      if "``thing''"!="" {
          tokenize "``thing''"
          local `thing'E `1'
          local `thing'C `2'
          if "`2'"=="" exit498 Too few arguments in `thing'(``thing'')
          if "`3'"!="" exit498 Too many arguments in `thing'(``thing'')
      }
   }
}

if "`imor'"!="" & "`logimor'"!="" {
    di as error "Can't specify both logimor() and imor()"
    exit 498
}
foreach arm in E C {
   if "`imor'"!="" {
      local logimor`arm'=log(`imor`arm'')
      if `imor`arm''==0 local logimor`arm' -50
      if `imor`arm''<0 {
         di as error "Negative IMOR in arm `arm'"
         exit 498
      }
   }
   if "`logimor'"!="" {
      if `logimor`arm''>50 & `logimor`arm''!=. local logimor`arm' 50
      if `logimor`arm''<-50 local logimor`arm' -50
      local imor`arm'=exp(`logimor`arm'')
      if `logimor`arm''<-50 local imor`arm' 0
   }
}

* check if logimor() etc are variables, and if so return an error
tempvar thing
gen `thing' = 0
foreach arg in logimorE logimorC imorE imorC sdlogimorE sdlogimorC corrlogimor {
	if mi("``arg''") continue
	qui replace `thing' = ``arg''
	summ `thing', meanonly
	local arg = subinstr("`arg'","C","",1)
	local arg = subinstr("`arg'","E","",1)
	if r(max)>r(min) {
		di as error "Sorry, metamiss cannot at present handle `arg'() varying across individuals"
		exit 498
	}
}

********************************* END OF COMMON PARSING *********************************

tempvar eststar sestar

********************************* IAN'S PARSING *********************************

if `syntax'==3 {

* SORT OUT METHOD
if "`method'"=="" {
   if `reps'~=-1 local method mc
   else local method gh
}
if `reps'==-1 local reps 100
if "`method'" != "mc" & "`method'" != "taylor" & "`method'" != "gh" {
   di as error "Please specify method(MC|Taylor|GH)."
}

* PARSE PRIORS FOR INFORMATIVENESS PARAMETERS logimor
if "`corrlogimor'"=="" local corrlogimor 0
assert `corrlogimor'>=-1 & `corrlogimor'<=1

********************************* END OF IAN'S PARSING *********************************

********************************* IAN'S CALCULATION *********************************

if "`debug'"=="debug" di as text "Starting Ian's calculation"

* SORT OUT DATA INCLUDING ZEROES
tokenize "`varlist'"
tempvar anyzero rE fE mE rC fC mC
qui {
   gen `anyzero' = min(`1', `2', `4', `5')==0
   gen `rE' = `1'+`anyzero'/2
   gen `fE' = `2'+`anyzero'/2
   gen `mE' = `3'
   gen `rC' = `4'+`anyzero'/2
   gen `fC' = `5'+`anyzero'/2
   gen `mC' = `6'
   count if `anyzero'==1
}
if r(N)>0 di as text "Zero cells detected: adding 1/2 to " as result r(N) as text " studies."

* PRINT DETAILS
di as text "Priors used:  Group 1: N(" as result "`logimorE'" as text "," as result "`sdlogimorE'" as text "^2). Group 2: N(" as result "`logimorC'" as text "," as result "`sdlogimorC'" as text "^2). Correlation: " as result "`corrlogimor'" as text "."

***************** TAYLOR AND GH METHODS **********************
qui {
   * BASICS AND VAROBS
   if "`method'"=="taylor" | "`method'"=="gh" {
      tempvar varobs varimor
      foreach R in E C {
         * "star" indicates values logimorE, logimorC
         tempvar imor`R' alpha`R' piobs`R' pimiss`R' pistar`R' dpistar`R'_dpi`R' dpistar`R'_dalpha`R' varpiobs`R' varalpha`R' varpistar`R' K`R'
         gen `imor`R'' = exp(`logimor`R'')
         replace `imor`R'' = 0 if `logimor`R'' < -50
         replace `imor`R'' = 9999 if `logimor`R'' > 50
         label var `imor`R'' "imor`R'"
         gen `alpha`R'' = `m`R''/(`r`R''+`f`R''+`m`R'')
         gen `piobs`R'' = `r`R''/(`r`R''+`f`R'')
         gen `pimiss`R'' = `imor`R''*`piobs`R''/(`imor`R''*`piobs`R''+1-`piobs`R'')
         gen `pistar`R'' = (1-`alpha`R'')*`piobs`R'' + `alpha`R''*`pimiss`R''
         gen `dpistar`R'_dpi`R'' = 1-`alpha`R''+`alpha`R''*`imor`R''/(`imor`R''*`piobs`R''+1-`piobs`R'')^2
         gen `dpistar`R'_dalpha`R'' = (`imor`R''-1)*`piobs`R''*(1-`piobs`R'')/(`imor`R''*`piobs`R''+1-`piobs`R'')
         gen `varpiobs`R'' = `piobs`R''*(1-`piobs`R'')/(`r`R''+`f`R'')
         gen `varalpha`R'' = `alpha`R''*(1-`alpha`R'')/(`r`R''+`f`R''+`m`R'')
         gen `varpistar`R'' = (`dpistar`R'_dpi`R'')^2 * `varpiobs`R'' + (`dpistar`R'_dalpha`R'')^2 * `varalpha`R''
         gen `K`R'' = `alpha`R''*`imor`R''/(`imor`R''*`piobs`R''+1-`piobs`R'')^2
      }

      if "`measure'"=="RD" gen `varobs' = `varpistarE' + `varpistarC'
      if "`measure'"=="RR" gen `varobs' = `varpistarE'/`pistarE'^2 + `varpistarC'/`pistarC'^2
      if "`measure'"=="OR" gen `varobs' = `varpistarE'/(`pistarE'*(1-`pistarE'))^2 + `varpistarC'/(`pistarC'*(1-`pistarC'))^2

      * TAYLOR METHOD
      if "`method'"=="taylor" {
            noi di in green "Method: Taylor series approximation."
            if "`measure'"=="RD" {
              gen `eststar' = `pistarE' - `pistarC'
              gen `varimor' = (`KE'*`piobsE'*(1-`piobsE')*(`sdlogimorE'))^2 + (`KC'*`piobsC'*(1-`piobsC')*(`sdlogimorC'))^2 - 2*`KE'*`KC'*`piobsE'*(1-`piobsE')*`piobsC'*(1-`piobsC')*(`sdlogimorE')*(`sdlogimorC')*(`corrlogimor')
            }
            if "`measure'"=="RR" {
              gen `eststar' = log(`pistarE' / `pistarC')
              gen `varimor' = (`KE'*(1-`piobsE')*(`sdlogimorE'))^2 + (`KC'*(1-`piobsC')*(`sdlogimorC'))^2 - 2*`KE'*`KC'*(1-`piobsE')*(1-`piobsC')*(`sdlogimorE')*(`sdlogimorC')*(`corrlogimor')
            }
            if "`measure'"=="OR" {
              gen `eststar' = log(`pistarE'/(1-`pistarE')) - log(`pistarC'/(1-`pistarC'))
              gen `varimor' = (`KE'*(`sdlogimorE'))^2 + (`KC'*(`sdlogimorC'))^2 - 2*`KE'*`KC'*(`sdlogimorE')*(`sdlogimorC')*(`corrlogimor')
            }
      }

      * GAUSS-HERMITE METHOD
      if "`method'"=="gh" {
         noi di as text "Method: Gauss-Hermite quadrature (" as result `nip' as text " integration points)."
         if "`debug'"=="debug" noi di as text "Performing GH quadrature"
         tempvar a1 a2
         gen `a1'=sqrt((1+(`corrlogimor'))/2)
         gen `a2'=sqrt((1-(`corrlogimor'))/2)
         local imorE exp(`logimorE'+(`sdlogimorE')*(`a1'*Z1+`a2'*Z2))
         local imorC exp(`logimorC'+(`sdlogimorC')*(`a1'*Z1-`a2'*Z2))
         foreach R in E C {
            local pimiss`R' (`piobs`R''*`imor`R'' / (1-`piobs`R''+`piobs`R''*`imor`R''))
            local pistar`R' ((1-`alpha`R'')*`piobs`R'' + `alpha`R''*`pimiss`R'')
         }
         if "`measure'"=="RD" local formula `pistarE'-`pistarC'
         if "`measure'"=="RR" local formula log(`pistarE') - log(`pistarC')
         if "`measure'"=="OR" local formula log(`pistarE'/(1-`pistarE')) - log(`pistarC'/(1-`pistarC'))
         qui expectn `formula', gen(`eststar') genvar(`varimor') token(Z1 Z2) nip(`nip') `debug'
      }
    gen `sestar' = sqrt(`varobs' + `varimor')
    }
}
***************** MONTE CARLO METHOD **********************
qui {
if "`method'"=="mc" {
   noi di in green "Method: Monte Carlo (`reps' draws)."
   foreach meas in `measure' {
      foreach adj in RAW STAR {
         foreach out in EST VAR {
            if "`replace'"=="noreplace" confirm new var `out'`adj'_`meas'
            else cap drop `out'`adj'_`meas'
         }
      }
   }

   if `reps'==-1 local reps 1000

   * PARSE PRIORS FOR P(MISSING)
   * note beta(c,d) with mean c/(c+d) is like c 1s and d 0s
   tokenize "`missprior'"
   if "`1'"=="" {
      local aE1 1
      local aE0 1
   }
   else {
      local aE1 `1'
      local aE0 `2'
   }
   if "`3'"=="" {
      local aC1 `aE1'
      local aC0 `aE0'
   }
   else {
      local aC1 `3'
      local aC0 `4'
   }

   * PARSE PRIORS FOR P(OUTCOME|NON-MISSING)
   tokenize "`respprior'"
   if "`1'"=="" {
      local bE1 1
      local bE0 1
   }
   else {
      local bE1 `1'
      local bE0 `2'
   }
   if "`3'"=="" {
      local bC1 `bE1'
      local bC0 `bE0'
   }
   else {
      local bC1 `3'
      local bC0 `4'
   }

   * COMPUTE DRAWS FROM POSTERIOR
   tempvar newlogimorE newlogimorC x0 x1 alphaE alphaC piobsE piobsC pistarE pistarC RAW_RD STAR_RD RAW_logRR STAR_logRR RAW_logOR STAR_logOR

   foreach meas in `measure' {
      foreach adj in RAW STAR {
         tempvar sum`adj'_`meas' sum`adj'_`meas'2
         gen `sum`adj'_`meas''  = 0
         gen `sum`adj'_`meas'2' = 0
      }
   }
   forvalues i=1/`reps' {
      if "`dots'"!="nodots" noi di "." _c
      if "`debug'"=="debug" {
         noi di "Drawing logimorE " _c
      }
      drawnorm `newlogimorE' `newlogimorC' `if' `in', means(`logimorE', `logimorC') sds(`sdlogimorE', `sdlogimorC') corr(1,`corrlogimor' \ `corrlogimor',1)
*      random `newlogimorE' `if' `in', dist(normal)
      if "`debug'"=="debug" {
         noi di " `logimorC' " _c
      }
*      random `newlogimorC' `if' `in', dist(normal) mean((`corrlogimor')*`newlogimorE') var(1-(`corrlogimor')^2)
      foreach R in E C {
*         replace `newlogimor`R'' = `logimor`R'' + `newlogimor`R''*(`sdlogimor`R'')

         if "`debug'"=="debug" {
            noi di " alpha`R' " _c
         }
*         random `x1' `if' `in', dist(gamma) mean(`a`R'1'+`m`R'')        var(`a`R'1'+`m`R'')
*         random `x0' `if' `in', dist(gamma) mean(`a`R'0'+`f`R''+`r`R'') var(`a`R'0'+`f`R''+`r`R'')
*         gen `alpha`R''=`x1'/(`x0'+`x1')
*         drop `x0' `x1'
         gen `alpha`R''=rbeta(`a`R'1'+`m`R'', `a`R'0'+`f`R''+`r`R'')


         if "`debug'"=="debug" {
            noi di " piobs`R' " _c
         }
*         random `x1' `if' `in', dist(gamma) mean(`b`R'1'+`r`R'') var(`b`R'1'+`r`R'')
*         random `x0' `if' `in', dist(gamma) mean(`b`R'0'+`f`R'') var(`b`R'0'+`f`R'')
*         gen `piobs`R''=`x1'/(`x0'+`x1')
*         drop `x0' `x1'
         gen `piobs`R''=rbeta(`b`R'1'+`r`R'', `b`R'0'+`f`R'')

         gen `pistar`R''=(1-`alpha`R'')*`piobs`R'' + `alpha`R''*exp(`newlogimor`R'')*`piobs`R'' / (1-`piobs`R''+exp(`newlogimor`R'')*`piobs`R'')
         drop `newlogimor`R'' `alpha`R''
      }

      gen `RAW_RD'=`piobsE'-`piobsC'
      gen `STAR_RD'=`pistarE'-`pistarC'

      gen `RAW_logRR'=log(`piobsE'/`piobsC')
      gen `STAR_logRR'=log(`pistarE'/`pistarC')

      gen `RAW_logOR' = log(`piobsE'/(1-`piobsE')) - log(`piobsC'/(1-`piobsC'))
      gen `STAR_logOR'= log(`pistarE'/(1-`pistarE')) - log(`pistarC'/(1-`pistarC'))

      drop `piobsE' `piobsC' `pistarE' `pistarC'

      foreach meas in `measure' {
         foreach adj in RAW STAR {
            replace `sum`adj'_`meas''  = `sum`adj'_`meas'' + ``adj'_`meas''
            replace `sum`adj'_`meas'2' = `sum`adj'_`meas'2' + ``adj'_`meas''^2
         }
      }

      drop `RAW_RD' `STAR_RD' `RAW_logRR' `STAR_logRR' `RAW_logOR' `STAR_logOR'

      if "`debug'"=="debug" {
         noi di " done." _newline
      }
   }
   * SUMMARISE POSTERIOR
   foreach meas in `measure' {
      foreach adj in RAW STAR {
         gen EST`adj'_`meas' = `sum`adj'_`meas''/`reps'
         gen VAR`adj'_`meas' = `sum`adj'_`meas'2'/`reps' - EST`adj'_`meas'^2
      }
   }
   if `nmeasures'==1 {
      gen `eststar' = ESTSTAR_`measure'
      gen `sestar' = sqrt(VARSTAR_`measure')
   }

}

} /* end of quietly */

********************************** END OF MONTE CARLO PART **********************************
}
********************************** END OF IAN'S CALCULATION **********************************

else if `syntax'==1 | `syntax'==2  {

********************************* JULIAN'S PARSING *********************************

if `syntax'==1 {
if "`aca'"=="aca" local method "Available cases analysis"
else if "`ica0'"! = ""  local method "ICA-0 (impute zeros)"
else if "`ica1'"! = ""  local method "ICA-1 (impute ones)"
else if "`icape'"!= ""  local method "ICA-pE (impute experimental group risk)"
else if "`icapc'"!= ""  local method "ICA-pC (impute control group risk)"
else if "`icap'"! = ""  local method "ICA-p (impute group-specific risk)"
else if "`icab'"! = ""  local method "ICA-b ('best'-case [0 in exp, 1 in control])"
else if "`icaw'"! = ""  local method "ICA-w ('worst'-case [1 in exp, 0 in control])"
else if "`icaimor'"!="" local method "ICA-IMOR (impute using IMORs `imorE' `imorC')"
}
else if `syntax'==2 {
if "`ica0'"! = ""  local list `list' ICA-0
if "`ica1'"! = ""  local list `list' ICA-1
if "`icape'"!= ""  local list `list' ICA-pE
if "`icapc'"!= ""  local list `list' ICA-pC
if "`icap'"! = ""  local list `list' ICA-p
if "`icaimor'"!="" local list `list' ICA-IMOR
local method "ICA-r combining `list'"
}
if "`aca'"!="aca" di as text "Method: `method'."

* quit if imor reasons indicated but IMOR not given
if ("`icaimor'" != "" & "`imorE'" == "") {
    di as error "icaimor specified but no IMOR provided"
    exit 498
}
* now if user has not entered imors, replace with ones to avoid errors in later code (coeff will always be zero)
foreach imorR in imorE imorC {
    if "``imorR''"=="" local `imorR' = 1
}

* check only one weight specified
local weight `w1'`w2'`w3'`w4'
if length("`weight'")>2 {
    di as error "Please specify only one of w1/w2/w3/w4."
    exit 498
}
if length("`weight'")==0 {
    if "`aca'"=="aca" {
       local w2 w2
       local weight w2
    }
    else {
       local w4 w4
       local weight w4
    }
}
if "`aca'"!="aca" {
   di as text "Weighting scheme: `weight'" _c
   if "`weight'"=="w1" di as error " (warning: standard errors may be too small)" _c
   di as text "."
}
************************ END OF JULIAN'S PARSING ****************************

************************ JULIAN'S CALCULATION ****************************
if "`debug'"=="debug" di as text "Starting Julian's calculation"

* SORT OUT DATA NOT INCLUDING ZEROES
tokenize "`varlist'"
tempvar rE fE mE rC fC mC
qui {
   gen `rE' = `1'
   gen `fE' = `2'
   gen `mE' = `3'
   gen `rC' = `4'
   gen `fC' = `5'
   gen `mC' = `6'
   foreach arm in E C {
      foreach reason in 0 1 pe pc p imor {
         tempvar num`reason'`arm'
         if "`aca'"!="" gen `num`reason'`arm'' = 0
         else if `syntax'==1 {
            if "`ica`reason''"=="ica`reason'" qui gen `num`reason'`arm'' = `m`arm''
            else gen `num`reason'`arm'' = 0
         }
         else if `syntax'==2 {
            if "`ica`reason'`arm''"!="" qui gen `num`reason'`arm'' = `ica`reason'`arm''
            else gen `num`reason'`arm'' = 0
         }
         char `num`reason'`arm''[varname] num`reason'`arm'
         local numlist`arm' `numlist`arm'' `num`reason'`arm''
      }
   }
   * IMPLEMENT ICA-B, ICA-W, ACCA
   if "`icab'"!="" replace `num1E' = `mE'
   if "`icab'"!="" replace `num0C' = `mC'
   if "`icaw'"!="" replace `num0E' = `mE'
   if "`icaw'"!="" replace `num1C' = `mC'
   if "`aca'"!="" replace `mE' = 0
   if "`aca'"!="" replace `mC' = 0
}

if "`listnum'"=="listnum" {
    di _newline "Missing values by reason, from syntax:" _c
    foreach arm in E C {
        l `id' `numlist`arm'' `if' `in', subvarname
    }
}

* for studies with missing data, but no reasons, assign number for each reason to overall proportion (across studies) with that reason
foreach arm in E C {
   tempvar numall`arm'
   qui egen `numall`arm'' = rsum(`num0`arm'' `num1`arm'' `numpe`arm'' `numpc`arm'' `nump`arm'' `numimor`arm'')
   summ `numall`arm'' `if' `in', meanonly
   local denom = r(sum)
   qui count if `m`arm''>0
   if `denom'==0 & r(N)>0 {
      di as error "Reason-specific counts sum to zero: can't assign missing data to reasons in " _c
      if "`arm'"=="C" di as error "control arm."
      if "`arm'"=="E" di as error "intervention arm."
      exit 498
   }
   foreach reason in 0 1 pe pc p imor {
      summ `num`reason'`arm'' `if' `in', meanonly
      qui replace `num`reason'`arm'' = r(sum) / `denom' if `numall`arm''==0
      qui replace `num`reason'`arm'' = 0 if `m`arm''==0
   }
}

* calculate basic sample sizes and probabilities
tempvar nE nC aE aC pE pC numallE numallC

* calculate numbers that indicate whether there are zero cells after imputations of 0s and 1s

tempvar temprE tempfE temprC tempfC addpointfive

quiet{
    gen `addpointfive' = 0
    gen `temprE' = `rE' + `num1E'*`mE'
    gen `tempfE' = `fE' + `num0E'*`mE'
    gen `temprC' = `rC' + `num1C'*`mC'
    gen `tempfC' = `fC' + `num0C'*`mC'
    replace `addpointfive' = 1 if (`temprE'==0 & `tempfE'~=0 & `temprC'~=0 & `tempfC'~=0 | `temprE'~=0 & `tempfE'==0 & `temprC'~=0 & `tempfC'~=0 | `temprC'==0 & `tempfC'~=0 & `temprE'~=0 & `tempfE'~=0  | `temprC'~=0 & `tempfC'==0 & `temprE'~=0 & `tempfE'~=0 )
    count if `addpointfive'==1
}
if r(N)>0 di as text "Zero cells detected: adding 1/2 to " as result r(N) as text " studies."

* create new data set with 0.5s added where necessary

quiet{
foreach arm in E C {
    replace `r`arm'' = `r`arm''+`addpointfive'/2
    replace `f`arm'' = `f`arm''+`addpointfive'/2
    gen `n`arm'' = `r`arm'' + `f`arm'' + `m`arm''
    gen `p`arm'' = `r`arm'' / (`r`arm'' + `f`arm'')
    gen `a`arm'' = `m`arm'' / `n`arm''
    foreach start in r f n p a n r f m num0 num1 numpe numpc nump numimor {
       char ``start'`arm''[varname] `start'`arm'
    }
}
}
if "`listx'"=="listx" {
    l `id' `rE' `fE' `nE' `pE' `aE' `rC' `fC' `nC' `pC' `aC' `if' `in', subvarname
}

* calculate probability of imputing each specific reason
* asup*E's sum to aE, asup*C's sum to aC
quiet {
   foreach arm in E C {
      gen `numall`arm'' = `num0`arm''+`num1`arm''+`numpe`arm''+`numpc`arm''+`nump`arm''+`numimor`arm''
      foreach reason in 0 1 pe pc p imor {
         tempvar asup`reason'`arm'
         gen `asup`reason'`arm'' = `a`arm''* `num`reason'`arm'' / `numall`arm'' if `numall`arm''~=0
         replace `asup`reason'`arm'' = 0 if `numall`arm''==0
         replace `num`reason'`arm'' = `n`arm''*`asup`reason'`arm''
      }
   }
}

if "`listall'"=="listall" {
    di _newline "Missing values by reason, scaled and imputed:" _c
    foreach arm in E C {
        l `id' `n`arm'' `r`arm'' `f`arm'' `m`arm'' `num0`arm'' `num1`arm'' `numpe`arm'' `numpc`arm'' `nump`arm'' `numimor`arm'' `if' `in', subvarname
    }
}

* calculate risks among missing participants using:
*       p`arm' - AC success fraction
*       a`arm' - fraction of missing data
*       asup`reason'`arm' - fraction of missing data by reason

foreach arm in E C {
    tempvar pfromimor`arm' p`arm'star
    if "`icaimor'"!="" qui gen `pfromimor`arm'' = (`p`arm''*`imor`arm'') / (`p`arm''*`imor`arm'' + 1 - `p`arm'')
    else gen `pfromimor`arm'' = 0
    qui gen `p`arm'star' = `p`arm''*(1-`a`arm'') + `asup1`arm'' + `asuppc`arm''*`pC' + `asuppe`arm''*`pE' + `asupp`arm''*`p`arm'' + `asupimor`arm''*`pfromimor`arm''
    char `pfromimor`arm''[varname] pfromimor`arm'
    char `p`arm'star'[varname] p`arm'star
}


* calculate new 'effective 2x2 table based on new overall risks, incorporating missing data imputations, with original sample size

tempvar erE efE enE erC efC enC

quiet {
    gen `erE' = `nE' * `pEstar'
    gen `efE' = `nE' * (1-`pEstar')
    gen `enE' = `erE' + `efE'
    gen `erC' = `nC' * `pCstar'
    gen `efC' = `nC' * (1-`pCstar')
    gen `enC' = `erC' + `efC'
}

if "`listp'"=="listp" {
    list `id' `pE' `pEstar' `pC' `pCstar' `if' `in', subvarname
}

tempvar RR lnRR W1selnRR W2selnRR OR lnOR W1selnOR W2selnOR VE VC W4varlnRR W4selnRR


if "`debug'"=="debug" di "Calculating estimates ..."
qui {
   if "`measure'"=="RD" gen `eststar' = `pEstar' - `pCstar'
   else if "`measure'"=="RR" gen `eststar' = log(`pEstar'/`pCstar')
   else if "`measure'"=="OR" gen `eststar' = log((`pEstar'/(1-`pEstar'))/(`pCstar'/(1-`pCstar')))
   else exit498 Error in local measure = `measure'
}

* calculate std errors
if "`debug'"=="debug" di "Calculating standard errors ..."
if "`weight'"=="w1" {
   if "`measure'"=="RD" qui gen `sestar' = sqrt(`erE'*(`enE'-`erE')/(`enE')^3 + `erC'*(`enC'-`erC')/(`enC')^3 )
   if "`measure'"=="RR" qui gen `sestar' = sqrt(1/`erE' - 1/`enE' + 1/`erC' - 1/`enC')
   if "`measure'"=="OR" qui gen `sestar' = sqrt(1/`erE' + 1/(`enE'-`erE') + 1/`erC' + 1/(`enC'-`erC'))
}
else if "`weight'"=="w2" {
   * ACA fixed-effect analysis for standard errors
   qui metan `rE' `fE' `rC' `fC', fixedi `or' `rr' `rd' nointeger nograph
   if "`measure'"=="RD" | !`oldmetan' qui gen `sestar' = _seES		// DF 15feb2022: if metan v4.0+, use _seES rather than _selogES
   else qui gen `sestar' = _selogES
}
else if "`weight'"=="w3" {
   * apply pstar probabilities to original ACA sample sizes (after zero cell adjustment)
   tempvar W3rE W3rC W3fE W3fC W3selnRR W3selnOR
   quiet {
      gen `W3rE' = `pEstar' * (`rE'+`fE')
      gen `W3rC' = `pCstar' * (`rC'+`fC')
      gen `W3fE' = (1-`pEstar') * (`rE'+`fE')
      gen `W3fC' = (1-`pCstar') * (`rC'+`fC')
      if "`measure'"=="RD" gen `sestar' = sqrt(`W3rE'*`W3fE'/(`rE'+`fE')^3 + `W3rC'*`W3fC'/(`rC'+`fC')^3)
      if "`measure'"=="RR" gen `sestar' = sqrt(1/`W3rE' - 1/(`rE'+`fE') + 1/`W3rC' - 1/(`rC'+`fC'))
      if "`measure'"=="OR" gen `sestar' = sqrt(1/`W3rE' + 1/`W3fE' + 1/`W3rC' + 1/`W3fC')
   }
}
else if "`weight'"=="w4" {
   * from formula derived from Taylor series
   * var of pstar adds over multiple reasons for missingness - here 0, 1, E, C, P and imor
   * formulae use 0-cell adjusted PSTARs, Ns, Ps
   tempvar OR_CE OR_EC varpEstar varpCstar firstsumE lastsumE firstsumC lastsumC

   quiet {
      gen `OR_CE' = (`pC' / (1-`pC')) / (`pE' / (1-`pE'))
      gen `OR_EC' = (`pE' / (1-`pE')) / (`pC' / (1-`pC'))

      replace `OR_CE' = 9999999 if (`pC'==1 | `pE'==0)
      replace `OR_EC' = 9999999 if (`pE'==1 | `pC'==0)

      gen `firstsumE' = - `asup0E' - `asup1E' + `asuppcE' * (`OR_CE' / ((`pE'*`OR_CE' + 1 - `pE' )^2) - 1) + `asupimorE' * (`imorE' / ((`pE'*`imorE' + 1 - `pE' )^2) - 1)
      gen `firstsumC' = - `asup0C' - `asup1C' + `asuppeC' * (`OR_EC' / ((`pC'*`OR_EC' + 1 - `pC' )^2) - 1) + `asupimorC' * (`imorC' / ((`pC'*`imorC' + 1 - `pC' )^2) - 1)

      gen `lastsumE' = `asup0E' * (`pEstar'^2 - (`pE'-`pEstar')^2) + `asup1E' * ((1-`pEstar')^2 - (`pE'-`pEstar')^2) + `asuppcE' * ((`pC'-`pEstar')^2 - (`pE'-`pEstar')^2) + `asupimorE' * ((`pfromimorE'-`pEstar')^2 - (`pE'-`pEstar')^2)
      gen `lastsumC' = `asup0C' * (`pCstar'^2 - (`pC'-`pCstar')^2) + `asup1C' * ((1-`pCstar')^2 - (`pC'-`pCstar')^2) + `asuppeC' * ((`pE'-`pCstar')^2 - (`pC'-`pCstar')^2) + `asupimorC' * ((`pfromimorC'-`pCstar')^2 - (`pC'-`pCstar')^2)

      gen `varpEstar' = (`pE' * (1-`pE') / (`nE' - `mE') ) * (1 + `firstsumE')^2 + (1/`nE') * ((`pE'-`pEstar')^2 + `lastsumE')
      gen `varpCstar' = (`pC' * (1-`pC') / (`nC' - `mC') ) * (1 + `firstsumC')^2 + (1/`nC') * ((`pC'-`pCstar')^2 + `lastsumC')

      if "`measure'"=="RD" gen `sestar' = sqrt( cond( `pEstar'==0, 0, `varpEstar') +  cond(`pCstar'==0, 0, `varpCstar'))

      if "`measure'"=="RR"  gen `sestar' = sqrt( cond( `pEstar'==0, 0, `varpEstar'/(`pEstar'^2)) +  cond(`pCstar'==0, 0, `varpCstar'/(`pCstar'^2)))

      if "`measure'"=="OR" gen `sestar' = sqrt( cond( `pEstar'==0, 0, `varpEstar'/((`pEstar'*(1-`pEstar'))^2) ) + cond(`pCstar'==0, 0, `varpCstar'/((`pCstar'*(1-`pCstar'))^2) ) )
   }
}

************************ END OF JULIAN'S CALCULATION ****************************
}

else if `syntax'==4 {
************************* GAMBLE-HOLLIS CALCULATION *****************************
qui {
   noi di "Method: Gamble-Hollis."
   * add halves if needed in ACA (to force their addition in ICA-b,ICA-w)
   tokenize "`varlist'"
   tempvar anyzero rE fE mE rC fC mC rmE fmE rmC fmC lower upper
   gen `anyzero' = min(`1', `2', `4', `5')==0
   gen `rE' = `1'+`anyzero'/2
   gen `fE' = `2'+`anyzero'/2
   gen `mE' = `3'
   gen `rC' = `4'+`anyzero'/2
   gen `fC' = `5'+`anyzero'/2
   gen `mC' = `6'
   gen `rmE'=`rE'+`mE'
   gen `fmE'=`fE'+`mE'
   gen `rmC'=`rC'+`mC'
   gen `fmC'=`fC'+`mC'
   count if `anyzero'==1
   if r(N)>0 noi di as text "Zero cells detected: adding 1/2 to " as result r(N) as text " studies."

   if "`rd'"!="rd" local logopt log
   metan `rE' `fE' `rC' `fC' `if' `in', `or' `rr' `rd' `options' `logopt' nointeger nograph
   gen `eststar' = _ES
   metan `rmE' `fE' `rC' `fmC' `if' `in', `or' `rr' `rd' `options' `logopt' nointeger nograph
   gen `upper' = _UCI
   metan `rE' `fmE' `rmC' `fC' `if' `in', `or' `rr' `rd' `options' `logopt' nointeger nograph
   gen `lower' = _LCI
   gen `sestar' = (`upper'-`lower')/(2*invnorm(1/2 + $S_level/200))
}

********************* END OF GAMBLE-HOLLIS CALCULATION **************************
}
* RUN METAN
if `nmeasures'==1 {
    cap drop _SS
    cap drop _SSmiss
    qui gen _SS = `rE'+`fE'+`rC'+`fC'
    qui gen _SSmiss = `mE'+`mC'
    label var _SS "Sample size (observed)"
    label var _SSmiss "Sample size (missing)"
    if "`measure'"=="RD" local semeasurename `measure'
    else local semeasurename log`measure'
    if "`measure'"!="RD" & "`log'"~="log" local eform eform
    if "`measure'"!="RD" & "`log'"=="log" local measurename log`measure'
    else local measurename `measure'
    local idopt = cond("`id'"~="","label(namevar=`id')","")
    di as text _newline "(Calling metan " _c
    if `"`idopt'`options'`eform'`graph'"'!="" di as text `"with options: `idopt' `options' `eform' `graph'"' _c
    di as text " ...)"
    if "`debug'"=="debug" di `"metan `eststar' `sestar' `if' `in', `idopt' `options' `eform' `graph'"'

	qui count if mi(`eststar', !`sestar')
	if r(N)>0 {
		di as error "Missing values of estimate and/or standard error found"
		foreach name in eststar sestar pEstar pCstar {
			char ``name''[varname] `name'
		}
		l `id' `eststar' `sestar' `pEstar' `pCstar' if mi(`eststar',! `sestar'), subvarname
		exit 498
	}

    metan `eststar' `sestar' `if' `in', `idopt' `options' `eform' `graph' effect(`log' `measure')
	if `oldmetan' {		// added DF 15feb2022
		qui gen _ES = `eststar'
		if "`eform'"=="eform" qui replace _ES=exp(_ES)
		if "`measure'"!="RD" & "`log'"!="log" local selog log
		qui gen _se`selog'ES = `sestar'
		label var _ES        "`measurename'"
		label var _LCI       "Lower CI (`measurename')"
		label var _UCI       "Upper CI (`measurename')"
		label var _se`selog'ES "se(`semeasurename')"
		move _ES _LCI
		move _se`selog'ES _LCI
	}
}
else {
    di _newline as text "More than one measure specified: meta-analysis not performed."
    di as text "You can run the meta-analyses using ESTSTAR* and VARSTAR*."
}

end

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

prog def expectn
syntax anything, gen(name) [genvar(name) nip(int 10) maxit(int 1000) tol(int 10) trace token(string) list debug]
confirm new var `gen'
if "`genvar'"!="" confirm new var `genvar'
if "`token'"=="" local token @
di as text "Finding expectation of `anything', assuming `token'~N(0,1) ..."
di as text "(using Gauss-Hermite quadrature with `nip' integration points)"
* GENERATE INTEGRATION POINTS AND WEIGHTS
local pim4 .751125544464925
local m=(`nip'+1)/2
forvalues i=1/`m' {
   if "`trace'"=="trace" di "Starting i=`i'"
   if `i'==1 local z=sqrt(float(2*`nip'+1))-1.85575*(2*`nip'+1)^(-0.16667)
   else if `i'==2 local z=`z'-1.14*`nip'^0.426/`z'
   else if `i'==3 local z=1.86*`z'-0.86*`zlag2'
   else if `i'==4 local z=1.91*`z'-0.91*`zlag2'
   else local z=2*`z'-`zlag2'

   local stop 0
   local iter 0
   while `stop'<1 {
      local ++iter
      local p1=`pim4'
      local p2=0
      forvalues j=1/`nip' {
         local p3=`p2'
         local p2=`p1'
         local p1=`z'*sqrt(2/`j')*`p2'-sqrt((`j'-1)/`j')*`p3'
      }
      local pp=sqrt(2*`nip')*`p2'
      local z1=`z'
      local z=`z1'-`p1'/`pp'
      local stop = abs(`z'-`z1')<10^(-`tol')
      if "`trace'"=="trace" di "Iteration `iter', change = " `z'-`z1'
      if `iter'>`maxit' {
         di as error "Too many iterations for i=`i'"
         exit 498
      }
   }
   if "`trace'"=="trace" di "i=`i' required `iter' iterations"

   if `i'==`m' & abs(`z')<1E-10 local z 0

   local j=`nip'+1-`i'
   local w`i' = 2/(`pp'*`pp'*sqrt(_pi))
   local x`i' = -`z'*sqrt(2)
   local w`j' = `w`i''
   local x`j' = -`x`i''
   local zlag2 `zlag1'
   local zlag1 `z'

}

if "`list'"=="list" {
   di as text _col(10) "Integration point" _col(30) "Weight"
   forvalues i=1/`nip' {
      di as result `i' _col(10) `x`i'' _col(30) `w`i''
   }
}

* DO THE INTEGRATION
tempvar argument sumarg2 sumarg sumw
gen `sumarg2' = 0
gen `sumarg' = 0
gen `sumw' = 0
qui gen `argument' = .
local wordcount=wordcount("`token'")
forvalues j=1/`wordcount' {
   local token`j'=word("`token'",`j')
}
if `wordcount'==1 {
   forvalues i1=1/`nip' {
      local arg : subinstr local anything "`token1'" "(`x`i1'')", all
      qui replace `argument' = `arg'
      qui count if missing(`argument')
      if r(N)>0 {
          di /*as error*/ "expectn: missing expression (`argument') in " r(N) " cases at `token1'=`x`i1''"
          l if missing(`argument')
*          exit 498
      }
      qui replace `sumarg2' = `sumarg2' + `w`i1''*(`argument')^2
      qui replace `sumarg' = `sumarg' + `w`i1''*(`argument')
      qui replace `sumw' = `sumw' + `w`i1''
if "`debug'"=="debug" noi di `arg', `sumw', `sumarg', `sumarg2'
   }
}
else if `wordcount'==2 {
   forvalues i1=1/`nip' {
      forvalues i2=1/`nip' {
         local arg : subinstr local anything "`token1'" "(`x`i1'')", all
         local arg : subinstr local arg "`token2'" "(`x`i2'')", all
*        note: the following syntax fails as it truncates strings
*        local arg = subinstr("`anything'","`token1'","(`x`i1'')",.)
         qui replace `argument' = `arg'
         qui count if missing(`argument')
         if r(N)>0 {
             di /*as error*/ "expectn: missing expression `argument' in " r(N) " cases at `token1'=`x`i1'', `token2'=`x`i2''"
             l if missing(`argument')
*             exit 498
         }
         local w=`w`i1''*`w`i2''
         qui replace `sumarg2' = `sumarg2' + `w'*(`argument')^2
         qui replace `sumarg' = `sumarg' + `w'*`argument'
         qui replace `sumw' = `sumw' + `w'
if "`debug'"=="debug" noi di `arg', `sumw', `sumarg', `sumarg2'
      }
   }
}
else if `wordcount'==3 {
   forvalues i1=1/`nip' {
      forvalues i2=1/`nip' {
         forvalues i3=1/`nip' {
            local arg : subinstr local anything "`token1'" "(`x`i1'')", all
            local arg : subinstr local arg "`token2'" "(`x`i2'')", all
            local arg : subinstr local arg "`token3'" "(`x`i3'')", all
            qui replace `argument' = `arg'
            qui count if missing(`argument')
            if r(N)>0 {
                di /*as error*/ "expectn: missing expression `argument' in " r(N) " cases at `token1'=`x`i1'', `token2'=`x`i2'', `token3'=`x`i3''"
                l if missing(`argument')
*                exit 498
            }
            local w=`w`i1''*`w`i2''*`w`i3''
            qui replace `sumarg2' = `sumarg2' + `w'*(`argument')^2
            qui replace `sumarg' = `sumarg' + `w'*(`argument')
            qui replace `sumw' = `sumw' + `w'
if "`debug'"=="debug" noi di `arg', `sumw', `sumarg', `sumarg2'
         }
      }
   }
}
else {
    di as error "This program hasn't yet been adapted to handle >3 Normal deviates"
    di as error "Please fiddle with the code"
    exit 498
}
cap assert abs(`sumw'-1)<1E-6
if _rc {
    di as error "Weights failed to sum to 1 in following cases:"
    l `sumw' if abs(`sumw'-1)>=1E-6
}
gen `gen' = `sumarg'/`sumw'
label var `gen' "Expectation of `anything'"
if "`genvar'"!="" {
    gen `genvar' = `sumarg2'/`sumw'-`gen'^2
    label var `genvar' "Variance of `anything'"
}
di as text "Variable `gen' created"
end

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

prog def exit498
di as error `"`0'"'
exit 498
end

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