*! version 1.0.2  TJS 17mar2000  (STB-58: sbe39.1)
program define metatrim
version 6.0
*
*  syntax: 
*  metatrim varlist [if] [in] [ , Vari CI REffect PRint Funnel EForm
*           LEvel(real 95) ESTimat(R|L|Q) GRaph IDvar(varname) 
*           SAve(filename [, REPLACE] * ]

if "`*'" == "" {
  di "Syntax is:"
  di in wh "metatrim " in gr "{ theta { se | var } | " _c
  di in gr "exp(theta) | ll ul [cl] } [ " in wh "if " _c
  di in gr "exp] [" in wh "in " in gr "range]"
  di in gr "              [ " in wh ", "  _c
  di in gr "{ " in wh "v" in gr "ar | " in wh "ci" in gr " } " _c
  di in wh "re" in gr "ffect " in wh "pr" in gr "int " in wh "est" _c
  di in gr "imator" in wh "(" in gr " { r | L | q } " in wh ")"
  di in gr "                         lev" in gr "el " in wh "f" _c 
  di in gr "unnel" in wh "(" in gr "#" in wh ") id" in gr "var" _c
  di in wh "(" in gr "var" in wh ") gr" in gr "aph graph_options ]"
  di _n in gr "           where { a | b |...} means choose" _c
  di in gr " one and only one of {a, b, ...}"
  exit
}

* Setup

syntax varlist(numeric min=2 max=4) [if] [in]  /*
  */ [ , Vari CI REffect Funnel LEvel(real 95) /*
  */     ESTimat(str) GRaph IDvar(varname) SAve(string) * ]
marksample touse
tokenize `varlist'

tempvar theta setheta var zz effect sr delta ceffect order
tempvar sgnrank tn minsrnk diff sdiff fill fillse id

local theta `1'
if "`3'" == "" { local setheta `2' }
  else {
    tempvar ll ul cl
    local ll `2'
    local ul `3'
    local cl `4'
  }

* input error traps

if "`ci'" != "" & "`vari'" != "" {
  di _n in re "Error: options 'ci' and 'var' cannot " _c
  di in re "be specified together."
  exit
}
if "`ci'" == "ci" & "`ul'" != "" {
  di _n in bl "Note: option 'ci' specified."
}
if "`ci'" == "ci" & "`ul'" == "" {
  di _n in re "Error: option 'ci' specified but varlist " _c
  di in re "has only 2 variables."
  exit
}
if "`ci'" != "ci" & "`vari'" != "vari" & "`ul'" != "" {
  di _n in bl "Warning: varlist has 3 variables but option " _c
  di in bl "'ci' not specified; 'ci' assumed."
  local ci "ci"
  local vari ""
}
if "`vari'" == "vari" & "`ul'" != "" {
  di _n in re "Error: option 'var' specified but varlist " _c
  di in re "has more than 2 variables."
  exit
}
if "`vari'" == "vari" & "`ul'" == "" {
  di _n in bl "Note: option 'var' specified."
}
if "`vari'" != "vari" & "`ul'" == "" {
  di _n in bl "Note: default data input format (theta, " _c
  di in bl "se_theta) assumed."
}
if "`estimat'" == "" { local est = "L" }
if "`estimat'" != "" {
  local est = substr(upper(trim("`estimat'")),1,1)
  if index("RLQ", "`est'") == 0 { 
    di _n in bl "Warning: invalid parameter for estimat(); " _c
    di    in bl "default linear estimator used."
    local est = "L"
  }
}

* Select data to analyze

if "`ul'" == "" { markout `touse' `theta' `setheta' }
           else { markout `touse' `theta' `ll' `ul' }

preserve                /* Note: data preserved here */
qui drop if !`touse'
if _N == 0 { error 2000 }
sort `theta'
qui gen `order' = _n

* Do calculations

if "`vari'" == "vari" { qui replace `setheta' = sqrt(`setheta')}

if "`ci'" == "ci" {
  capture confirm variable `cl'
  if _rc~=0 { qui gen `zz'  = invnorm(.975) }
  else { qui replace `cl' = `cl' * 100 if `cl' < 1
    qui gen `zz' = -1 * invnorm((1- `cl' / 100) / 2 )
    qui replace `zz' = invnorm(.025) if `zz'==.
  }
  qui gen   `setheta' = ( ln(`ul') - ln(`ll')) / 2 / `zz'
  qui replace `theta' = ln(`theta')
}

local estname = "Linear"
if "`est'" == "R" { local estname = "Run" }
if "`est'" == "Q" { local estname = "Quadratic" }

local matname = "Fixed-effects model"
if "`reffect'" == "reffect" { local matname = "Random-effects model" }

if "`idvar'" != "" {
  local sf: format `idvar'
  local sf = substr("`sf'", 2, length("`sf'")-2)
  local sn = max(int(real("`sf'")),5)
  qui gen str`sf' `id'=`idvar'
}
  else { 
    local len = int(log10(_n)) + 8
    qui gen str`len' `id' = "study " + string(_n)
  }

local          n = _N
local       trim = 0
local       flag = 1
local          i = 0
qui gen     `sr' = 0
qui gen `effect' = `theta'

* Looped Calculations 

while `flag' {
  `q' meta `effect' `setheta' `if' 
  if "`reffect'" == "reffect" { local d = $S_7 }
                         else { local d = $S_1 }
  qui {
    gen       `delta' = `d'
    gen     `ceffect' = `theta' - `delta'
    egen    `sgnrank' = rank(abs(`ceffect'))
    replace `sgnrank' = `sgnrank' * sign(`ceffect') 
    egen    `minsrnk' = min(`sgnrank')
    egen    `tn'      = sum(`sgnrank') if `sgnrank' > 0
    sort `order' 
    local Tn = `tn'[_N]
    if "`est'" == "R" { local trim = `n' - abs(`minsrnk') - 1.5 }
      else if "`est'" == "Q" { 
        local trim = `n' -0.5 - sqrt(2 * `n'^2 - 4 * `Tn' + 0.25) }
      else { 
        local trim = (4 * `Tn' - `n' * (`n' + 1)) / (2 * `n' - 1) }
    local trim = int(max(0, `trim' + 0.5))
    gen     `diff'  = `sr' - `sgnrank'
    replace `sr'    = `sgnrank'
    egen    `sdiff' = sum(abs(`diff'))
    if `sdiff' == 0 | `i' > 9 { local flag = 0 }
    local i = `i' + 1
    local q "qui"
    sort `order' 
    qui replace `effect' = . if _n > _N - `trim' 
  }
  if `i' == 1 {
    di
    di in gr "Trimming estimator: " in ye "`estname'"
    di in gr "Meta-analysis type: " in ye "`matname'"  
    di 
    di in gr "iteration |  estimate    Tn    # to trim     diff"
    di in gr "----------+--------------------------------------"
  }
  di in ye %5.0f `i' in gr "     |" in ye %9.3f `d' %7.0f `Tn' _c
  di in ye %10.0f `trim' %12.0f `sdiff'[_N]

/* 
  di 
  di "         theta      delta    ceffect    sgnrank" _c
  di "         tn       diff" _c
  list `theta' `delta' `ceffect' `sgnrank' `tn' `diff'
*/

  cap drop `delta' `ceffect' `sgnrank' `tn' `diff' `sdiff' `minsrnk'
}

* Filled analysis setup

if `trim' == 0 { di _n in bl "Note: no trimming performed; data unchanged" }
if `i' > 9     { di _n in bl "Warning: iterative algorithm did not converge" }
local new = _N + `trim'
sort `order'
qui {
  if `trim' > 0 { set obs `new' }
  gen       filled = `theta'
  gen     `fillse' = `setheta'
  gen       `fill' = .
  gen      `delta' = `d'
  replace   filled = 2 * `delta' - filled[_n - `trim'] if _n > _N - `trim'
  replace `fillse' =             `fillse'[_n - `trim'] if _n > _N - `trim'
  replace   `fill' =                            filled if _n > _N - `trim'
  replace     `id' =     "fill " + string(_N - _n + 1) if _n > _N - `trim'
}
sort filled `order'
label var filled "theta, filled"
if "`graph'" != "" {
  if "`reffect'" != "" { local gr = "gr(r)" }
                  else { local gr = "gr(f)" }
}
di _n in gr "Filled " _c

* Do filled meta-analysis

meta filled `fillse', id(`id') `gr' `options'

* Graph a funnel plot

if "`funnel'" == "funnel" {
  if "`reffect'" != "" { global Xtheta = $S_7 }
                  else { global Xtheta = $S_1 }
 
  fungph filled `fillse' `fill', level(`level') `options'
}

* Save filled data?

 if "`save'" != "" {
   local c = index("`save'",",")
   if `c' != 0 {
     local save = substr("`save'",1,`c'-1) + substr("`save'",`c'+1, .)
   }
   local save1 : word 1 of `save'
   local replace : word 2 of `save'
   if "`replace'" == "replace" { capture erase `save1'.dta }
   capture confirm new file `save1'.dta
   if _rc == 0 {
     qui keep `id' filled `fillse'
     rename `id' id
     rename `fillse' fillse
     lab var id "study ID"
     lab var fillse "S.E. of theta, filled"
     lab data "trimed and filled $S_FN"
     di 
	save `save1'
   }
   else {
     local rc = _rc
     di _n in re "File `save1' exists. Use 'replace' option: save(filename, replace)."
     exit `rc'
   }
 }

exit
end

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

program define fungph
version 6.0

* creates a filled funnel plot graph

* Setup

syntax varlist(min=3 max=3) [if] [in] [, CI L1title(str)  /*
  */ L2title(str) Connect(str) Symbol(str) SOrt Pen(str)  /*
  */ T2title(str) B2title(str) XLAbel YLAbel YLAbel(str)  /*
  */ XLAbel(str) LEVel(integer $S_level) GAp(str) GWeight /*
  */ EForm /* ...the following are trapped but never used...
  */ GRaph(str) PRint CLine XLIne(str) FMult(str) EBayes  /*
  */ BOXYsca(str) BOXSHad(str) LTRunc(str) RTRunc(str)    /*
  */  * ]
tokenize `varlist'

tempvar theta setheta

local theta    `1'
local setheta  `2'
local fill     `3'

preserve        /* Note: data preserved here */

* Graph options

if "`connect'" == "" { local connect "lll.." }
  else {
    local lll = length("`connect'")
    if      `lll' == 1 { local connect "`connect'll.." }
    else if `lll' == 2 { local connect "`connect'l.." }
    else if `lll' == 3 { local connect "`connect'.." }
    else if `lll' == 4 { local connect "`connect'." }
  }
local connect "co(`connect')"

if "`symbol'" == "" { local symbol "iiioS" }
  else {
    local lll = length("`symbol'")
    if      `lll' == 1 { local symbol "`symbol'iioS" }
    else if `lll' == 2 { local symbol "`symbol'ioS" }
    else if `lll' == 3 { local symbol "`symbol'oS" }
    else if `lll' == 4 { local symbol "`symbol'S" }
  }
local symbol "sy(`symbol')"

if "`pen'" == "" { local pen "35522" }
  else {
    local lll = length("`pen'")
    if      `lll' == 1 { local pen "`pen'5522" }
    else if `lll' == 2 {
    local pen = "`pen'" + substr("`pen'",2,1) + "22"
  }
  else if `lll' == 3 { local pen "`pen'22" }
  else if `lll' == 4 {
    local pen = "`pen'" + substr("`pen'",4,1)
  }
}
local pen "pen(`pen')"

if "`l2title'" == "" {
  local l2title : variable label `theta'
  if "`l2title'" == "" { local l2title "`theta'" }
}
local titlel2 = "`l2title'"
if "`ci'" == "" { local l2title "l2(`l2title')" }
           else { local l2title "l2(log[`l2title'])" }

if "`l1title'" == "" { local l1title "" "" }
local l1title "l1(`l1title')"

if "`b2title'" == "" { local b2title "`titlel2'" }
if "`ci'" == "" { local b2title = "b2(s.e. of: `b2title')" }
           else { local b2title = "b2(s.e. of: log[`b2title'])" }

if `"`t2title'"' == `""' {
  local t2title = "Filled funnel plot with pseudo"
  local t2title = "`t2title' `level'% confidence limits"
  local t2title "t2(`"`t2title'"')"
}
  else if `"`t2title'"' == `"."' { local t2title }  
  else { local t2title "t2(`"`t2title'"')" }

if "`xlabel'" == "" { local xlabel "xla" }
  else { local xlabel "xlabel(`xlabel')" }

if "`ylabel'" == "" { local ylabel "yla" }
  else { local ylabel "ylabel(`ylabel')" }

if "`sort'" == "" { local sort "sort" }

if "`gap'" == "" { local gap "gap(3)" }
  else { local gap "gap(`gap')" }

* calculations for plot

tempvar ll2 ul2 RRm 

qui {
  if `level' < 1 { local `level' =`level' * 100 }
  local z  = -1 * invnorm((1 - `level' / 100) / 2)
  local obs1=_N+1
  set obs `obs1'
  replace `setheta'=0 in `obs1'
  replace `theta' = . in `obs1'
  if "`eform'" != "" { global Xtheta = log($Xtheta) } 
  gen `RRm' = $Xtheta
  gen `ll2' = `RRm' - `z' * `setheta'
  gen `ul2' = `RRm' + `z' * `setheta'
}

if "`gweight'" != "" {
  tempvar ww
  qui gen `ww' = `w'
  qui replace `ww' = 1 if `setheta' == 0 
  local gw "[w=`ww']"
}
  else { local gw "" }

* do plot

#delimit ;
graph `RRm' `ll2' `ul2' `theta' `fill' `setheta' `gw',
  `connect' `symbol' `t2title' `l2title' `b2title'
  `l1title' `xlabel' `ylabel' `sort' `pen' `gap' `options';
#delimit cr

exit
end