*! version 1.12 jacs/ss May 1997  STB-38 sbe16

program define meta
     version 5.0
     local varlist "req ex min(2) max(2)"
     local if "opt"
     local in "opt"
     local options "EForm ID(string) Level(integer $S_level) PRint EBayes"
     local options "`options' GRaph(string) YLABel(string) XLOG"
     local options "`options' Symbol(string) FMult(real 1)  BOXYsca(real 1)"
     local options "`options' BOXSHad(integer 0) CLine YTick GAP"
     local options "`options' LTRunc(string) RTRunc(string) *"

     parse "`*'"

     if "`id'"~=""{
          confirm variable `id'
     }

     if "`graph'"~="" {
       capture assert "`graph'"=="r" | "`graph'"=="f" | "`graph'"=="e"
       if _rc~=0 {
         di in re "Must specify graph(f) or graph(r) or graph(e)"
         exit 198
       }
       if "`graph'"=="e" {
         local ebayes "ebayes"
         local print "print"
       }
     }

     if "`ylabel'"~="" {
       di in re "ylabel option not permitted"
       exit 198
     }

     if "`xlog'"~="" {
       di in re "xlog option not permitted (use eform option)"
       exit 198
     }

     if "`symbol'"~="" {
       di in re "symbol option not permitted"
       exit 198
     }

     if "`ytick'"~="" {
       di in re "ytick option not permitted"
       exit 198
     }

     if "`gap'"~="" {
       di in re "gap option not permitted"
       exit 198
     }

     if "`fmult'"~="" {
       capture assert `fmult'>0
       if _rc~=0 {
         di in re "Label font scaling factor must be >0"
         exit 198
       }
       local fmult "fmult(`fmult')"
     }

     if "`boxysca'"~="" {
       capture assert `boxysca'<=1 & `boxysca'>0
       if _rc~=0 {
         di in re "Y scaling factor for box must be between 0 and 1"
         exit 198
       }
       local boxysca "boxysca(`boxysca')"
     }

     if "`boxshad'"~="" {
       capture assert `boxshad'<=4 & `boxshad'>=0
       if _rc~=0 {
         di in re "Box shading must be between 0 and 4"
         exit 198
       }
       local boxshad "boxshad(`boxshad')"
     }

     tempvar touse w wpsi qi w2 wstar wstarps v
     tempname sumwpsi sumw sumw2 cappsi rcappsi lc uc rlc ruc Q k tausq vareb seeb
     tempname swstar swstarp  rt ft slab

     parse "`varlist'", parse(" ")
     local psi `1'
     local se `2'
     gen `v'=`se'^2

     preserve

     quietly {
       mark `touse' `if' `in'
       markout `touse' `psi' `v'
       keep if `touse'
     }


     * FIXED EFFECTS
     gen `w'=1/`v'

     gen `wpsi'=`w'*`psi'

     qui summ `wpsi'
     scalar def `sumwpsi'=_result(3)*_result(1)
     local k=_result(1)
     qui summ `w'
     scalar def `sumw'=_result(3)*_result(1)

     scalar def `cappsi'=`sumwpsi'/`sumw'
     scalar def `lc'=`cappsi'-invnorm(`level'*0.005 + 0.5)*(`sumw'^(-0.5))
     scalar def `uc'=`cappsi'+invnorm(`level'*0.005 + 0.5)*(`sumw'^(-0.5))
     scalar def `ft'=`cappsi'/(`sumw'^(-0.5))     

     * RANDOM EFFECTS
     gen `qi'=`w'*((`psi'-`cappsi')^2)
     qui summ `qi'
     scalar def `Q'=_result(3)*_result(1)
          gen `w2'=`w'^2
     qui summ `w2'
     scalar def `sumw2'=_result(3)*_result(1)
     scalar def `tausq'=max(0,(`Q'-(`k'-1))/(`sumw'-(`sumw2'/`sumw')))
     gen `wstar'=(`v'+`tausq')^(-1)
     
     gen `wstarps'=`wstar'*`psi'
     qui summ `wstarps'
     scalar def `swstarp'=_result(3)*_result(1)
     qui summ `wstar'
     scalar def `swstar'=_result(3)*_result(1)

     scalar def `rcappsi'=`swstarp'/`swstar'
     scalar def `rlc'=`rcappsi'-invnorm(`level'*0.005 + 0.5)*((`swstar')^(-0.5))
     scalar def `ruc'=`rcappsi'+invnorm(`level'*0.005 + 0.5)*((`swstar')^(-0.5))
     scalar def `rt'=`rcappsi'/((`swstar')^(-0.5))
     if "`eform'"~="" {
           scalar define `cappsi'=exp(`cappsi')
           scalar define `lc'=exp(`lc')
           scalar define `uc'=exp(`uc')
           scalar define `rcappsi'=exp(`rcappsi')
           scalar define `rlc'=exp(`rlc')
           scalar define `ruc'=exp(`ruc')
           local ef=" (exponential form)"
           local efu="-------------------"
     }


gen Est=`cappsi'
gen Lower=`lc'
gen Upper=`uc'
gen z_value=`ft'
gen p_value=2*min((1-normprob(`ft')),normprob(`ft'))

format Est Lower Upper z_value p_value %6.3f


     di in gr _n "Meta-analysis of `k' studies `ef'"
     di in ye "----------------------------`efu'"
     di in gr _n "Fixed and random effects pooled estimates, "
     di in gr "lower and upper `level'% confidence limits, and"
     di in gr "asymptotic z-test for null hypothesis that true effect=0"_n 
     
     di in ye _n "Fixed effects estimation"

list Est Lower Upper z_value p_value if _n==_N, noobs nodisplay


     di in bl _n "Test for heterogeneity: Q= " in ye %6.3f /*
     */ `Q' in bl " on " in ye (`k'-1) in bl " degrees of freedom (p=" /*
     */ in ye %6.3f chiprob(`k'-1,`Q') in bl ")"
     di in bl  "Der Simonian and Laird estimate of between studies variance = " in ye %6.3f /*
     */ `tausq'

     di in ye _n "Random effects estimation"

qui {
replace Est=`rcappsi'
replace Lower=`rlc'
replace Upper=`ruc'
replace z_value=`rt'
replace p_value=2*min((1-normprob(`rt')),normprob(`rt'))
}

format Est Lower Upper z_value p_value %6.3f

list Est Lower Upper z_value p_value if _n==_N, noobs nodisplay

     if "`print'"~="" {
     di in gr _n(2)"Weights given to each study in fixed and random effects estimation, "
     di in gr "estimates of effect in each study, "
     di in gr "and lower and upper `level'% confidence limits"_n 
     

     if "`ebayes'"~="" & `tausq'~=0 {
     di in bl _n"Note: estimates and confidence limits are empirical Bayes"
     }
     if "`ebayes'"~="" & `tausq'==0 {
     di in bl _n"Note: between studies variance is 0, so empirical Bayes estimates "
     di in bl "cannot be calculated - estimates and confidence limits reported are "
     di in bl "from the data"
     }
     if "`ebayes'"=="" {
     di in bl _n"Note: estimates and confidence limits are from the data"
     }
     
     
     gen Fixed=`w'
     gen Rand=`wstar'

    if "`ebayes'"~="" & `tausq'~=0{
       if "`eform'"~="" {
         scalar define `rcappsi'=log(`rcappsi')
       }  
       qui replace Est=(`psi'/(`v') + `rcappsi'/`tausq') / (1/(`v') + 1/`tausq')
       qui gen `vareb'=(`tausq'*`v')/(`tausq'+`v')+((`v'/(`tausq'+`v'))^2)/`swstar'
       qui gen `seeb'=sqrt(`vareb')
       qui replace Lower=Est-invnorm(`level'*0.005 + 0.5)*sqrt(`vareb')
       qui replace Upper=Est+invnorm(`level'*0.005 + 0.5)*sqrt(`vareb')
       if "`eform'"~="" {
         scalar define `rcappsi'=exp(`rcappsi')
       }  
     }

     if "`ebayes'"=="" | ("`ebayes'"~="" & `tausq'==0) {
     qui replace Est=`psi'
     qui replace Lower=`psi'- invnorm(`level'*0.005 + 0.5)*(`v'^(0.5))
     qui replace Upper=`psi'+ invnorm(`level'*0.005 + 0.5)*(`v'^(0.5))
     }
   
     if "`eform'"~="" {
     qui replace Est=exp(Est)
     qui replace Lower=exp(Lower)
     qui replace Upper=exp(Upper)
     }
     else {}
     
     format Fixed Rand Est Lower Upper %6.2f
 
     
     if "`id'"~="" {
     gen str12 Study=`id'
     }
     else {
     gen Study=_n
     }
        
     list Study Fixed Rand Est Lower Upper, noobs nodisplay
     }

    
     if "`id'"~="" {
       local id "id(`id')"
     }

     if "`graph'"=="r" | "`graph'"=="e" {
       local cappsi=`rcappsi'
       local lc=`rlc'
       local uc=`ruc'
     }

     if "`graph'"~="" {
       if "`ltrunc'"~="" {
         local ltrunc "ltrunc(`ltrunc')"
       }
       if "`rtrunc'"~="" {
         local rtrunc "rtrunc(`rtrunc')"
       }
       if "`graph'"~="e" {
         metagrph `psi' `v', `id' cappsi(`cappsi') poollc(`lc') pooluc(`uc') /*
         */ `eform' `fmult' `boxysca' `boxshad' `cline' `ltrunc' `rtrunc' /*
         */ `options'
       }

       if "`graph'"=="e" {
        if `tausq'==0 {
         di in re "Note: between studies variance is 0, so can't plot empirical Bayes estimates "
        }
        else {
          if "`eform'"~="" {
            qui replace Est=log(Est)
          }
          local xvarlab : variable label `psi'
          if "`xvarlab'"~="" {
            label var Est "`xvarlab' (Empirical Bayes)"
          }
          else {
            label var Est "Empirical Bayes estimate"
          } 
          metagrph Est `vareb', `id' cappsi(`cappsi') poollc(`lc') /*
          */ pooluc(`uc') `eform' `fmult' `boxysca' `boxshad' `cline' `ltrunc' /*
          */ `rtrunc' `options'
         }
       }
     }

     global S_1 = `cappsi'
     global S_2 = (`sumw'^(-0.5))
     global S_3 = `lc'
     global S_4 = `uc'
     global S_5 = `ft'
     global S_6 = tprob(`k'-1,`ft')
     global S_7 = `rcappsi'
     global S_8 = (`swstar'^(-0.5))
     global S_9 = `rlc'
     global S_0 = `ruc'
     global S_11 = `rt'
     global S_12 = tprob(`k'-1,`rt')
     global S_13 = `tausq'

     if "`ebayes'"~="" & `tausq'~=0 {
       restore
       cap drop ebest
       cap drop ebse

       if "`eform'"~="" {
         scalar define `rcappsi'=log(`rcappsi')
       }  
       gen ebest = (`psi'/(`v') + `rcappsi'/`tausq') / (1/(`v') + 1/`tausq') 
       gen ebse=sqrt((`tausq'*`v')/(`tausq'+`v')+((`v'/(`tausq'+`v'))^2)/`swstar')

       if "`eform'"~="" {
         qui replace ebest=exp(ebest)
         preserve
       }
       else preserve
     }    

    restore

end

program define metagrph
  version 5.0

  local varlist "req ex min(2) max(2)"

  local options "ID(string) CAPPSI(string) POOLLC(string)"
  local options "`options' POOLUC(string) SAving(string) EFORM"
  local options "`options' LEVEL(integer $S_level) FMult(real 1)"
  local options "`options' BOXYsca(real 1) BOXSHad(integer 0) CLine"
  local options "`options' LTRunc(string) RTRunc(string) *"
  parse "`*'"

  tempvar se obsno lci uci idlen
  tempname obslab k

  parse "`varlist'", parse(" ")
  local psi `1'
  local v `2'
  
  gen `obsno'=_n
  gsort -`obsno'

  if _N>20 {
    local fdiv1=20/_N
  }
  else {
    local fdiv1=1
  }

  local k=_N
  quietly {
    gen `se'=sqrt(`v')
    gen `lci'=`psi'-invt(`k'-1,`level'*0.005 + 0.5)*`se'
    gen `uci'=`psi'+invt(`k'-1,`level'*0.005 + 0.5)*`se'

    if "`eform'"~="" {
      replace `psi'=exp(`psi')
      replace `lci'=exp(`lci')
      replace `uci'=exp(`uci')
      local xlog "xlog"
    }

    if "`ltrunc'"~="" {
      quietly summ `psi'
      if `ltrunc'>_result(5) {
        di in re "Left truncation must be less than all effect estimates"
        exit 198
      }
      quietly replace `lci'=`ltrunc' if `lci'<`ltrunc'
    }
    if "`rtrunc'"~="" {
      quietly summ `psi'
      if `rtrunc'<_result(6) {
        di in re "Right truncation must be greater than all effect estimates"
        exit 198
      }
      quietly replace `uci'=`rtrunc' if `uci'>`rtrunc'
    }

    if "`saving'"~="" {
      local saving "saving(`saving')"
    }
  
  quietly replace `obsno'=_n+1
  local i 1
  local ylab="0"
  local ytick="2"
  while `i'<=_N {
    local i=`i'+1

    if `i'>2 {
      local ytick "`ytick',`i'"
    }
  }
  local i=`i'+1

  if _N<26 {
    local ytick "ytick(`ytick')"
  }
  else {
    local ytick ""
  }

  local nobs1=_N+1
  local nobs2=_N+2
  quietly {
    set obs `nobs2'
    summ `lci'
    replace `psi'=_result(5) in `nobs1'
    summ `uci'
    replace `psi'=_result(6) in `nobs2'

    replace `obsno'=0 in `nobs1'
    replace `obsno'=`nobs2' in `nobs2'

    label define `obslab' 0 "Combined" `nobs2' " ", add
    label values `obsno' `obslab'
  }

  sort `obsno'

  graph `obsno' `psi', ylab(`ylab') `ytick' s(i) gap(10) `xlog' `options'

  parse "$S_G1", parse(",")
* noi display "* `*'"
  local leftgph `3'
  local dr `9'
  local dc `11'

  parse "$S_G2", parse(",")
  local leftdat `3'
  local rgttext=(`leftdat'-`leftgph')*.75

  local imax=_N
  local i=2
  tempvar boxwid

  quietly gen `boxwid'=1/`se'

  if "`id'"~="" {
    gen `idlen'=length(`id')
    quietly summarize `idlen'
    local idleng=_result(6)
    
    if `idleng'>8 {
      local fdiv2=8/`idleng'
    }
    else {
      local fdiv2=1
    }
  
    local fdiv=min(`fdiv1',`fdiv2')
*   noi display "`fdiv'  `fdiv1'  `fdiv2'"

    local dr=`dr'*.7*`fdiv'*`fmult'
    local dc=`dc'*.7*`fdiv'*`fmult'

/*
    if `fdiv'<1 {
      noi display ""
      noi display "Font reduction"
      noi display "`fdiv' `dr' `dc'"
    }
*/
  }

  gph open, `saving'
  graph
  local ay=_result(5)
  local by=_result(6)
  local ax=_result(7)
  local bx=_result(8)

  while `i'<`imax' {

* row value
    local st=`obsno'[`i']
    local row=(`st'*`ay') + `by'

* label y axis
    if "`id'"~="" {
      gph font `dr' `dc'
      local textrow=`row'+(`dc'/2)
      local chari=`id'[`i']

      gph text `textrow' `rgttext' 0 1 `chari'
    }

* draw box with area proportional to inverse variance

    local mu=`psi'[`i']

    if "`eform'"~="" {
      local col=(log(`mu')*`ax') + `bx'
    }
    else {
      local col=(`mu'*`ax') + `bx'
    }

* derive maximum box size (yrange)
    quietly summ `obsno'
    local ymin=_result(5)
    local ymax=_result(6)
    local ynum=_result(1)
    local rmin=(`ymin'*`ay') + `by'
    local rmax=(`ymax'*`ay') + `by'
    local yrange=abs(`rmax'-`rmin')/(2*`ynum')

* draw box
    quietly summ `boxwid'
    local widmax=_result(6)
    local width=`boxwid'[`i']
    local mult=`yrange'*`width'/`widmax'
    local xscale=320/233
    local boxlr=`row'-(`mult'*`boxysca')
    local boxur=`row'+(`mult'*`boxysca')
    local boxlc=`col'-(`mult'/`xscale')
    local boxuc=`col'+(`mult'/`xscale')
    gph box `boxlr' `boxlc' `boxur' `boxuc' `boxshad '

* confidence interval
    local lc=`lci'[`i']
    local uc=`uci'[`i']

    if "`eform'"~="" {
      local cleft=(log(`lc')*`ax') + `bx'
      local cright=(log(`uc')*`ax') + `bx'
    }
    else {
      local cleft=(`lc'*`ax') + `bx'
      local cright=(`uc'*`ax') + `bx'
    }
    gph line  `row' `cleft' `row' `cright'

    local i=`i'+1
  }

* diamond for overall estimate

* row value
  local st=`obsno'[1]
  local row=(`st'*`ay') + `by'
  local rowup=`row'+(`yrange'/3)
  local rowdn=`row'-(`yrange'/3)

  if "`eform'"~="" {
    local cmiddle=(log(`cappsi')*`ax') + `bx'
    local cleft=(log(`poollc')*`ax') + `bx'
    local cright=(log(`pooluc')*`ax') + `bx'
  }
  else {
    local cmiddle=(`cappsi'*`ax') + `bx'
    local cleft=(`poollc'*`ax') + `bx'
    local cright=(`pooluc'*`ax') + `bx'
  }
  gph line  `rowup' `cmiddle' `row' `cright'
  gph line  `row' `cright' `rowdn' `cmiddle'
  gph line  `rowup' `cmiddle' `row' `cleft'
  gph line  `row' `cleft' `rowdn' `cmiddle'

* dotted line at the combined estimate
  if "`cline'"~="" {
    local top=`obsno'[_N-1] + 0.5
    local rowup=(`top'*`ay') + `by'
    local incr=(`rowup'-`rowdn')/100
    local j 0
    while `j'<50 {
      local i=2*`j'
      local lower=`rowdn'+(`i'*`incr')
      local upper=`rowdn'+((`i'+1)*`incr')
      gph line `lower' `cmiddle' `upper' `cmiddle'
      local j=`j'+1
    }
  }
  gph close

  }
end