*!cart.ado WvP ,  20/5/2002
*!first version 6/4/98
*!id() option added to save identifier in _CART.dta
program define cart
version 7.0
local varlist "req ex min(1) "
local if "opt"
local in "opt"
#delimit ;
local options "Pval(real 0.05) PNOMinal MINSize(int 10) MINFail(int 10) Time(string)
               Fail(string)  NAme(string) SUmby(string) TABby(string)
         Strata(string) Adjust(string) AT(string) ID(string) ";
#delimit cr
local cmd cart `*'
parse "`*'"
quietly {
  preserve
  if "`if'`in'"~="" { keep `if'`in' }
  local if
  local in
  if "`name'"=="" {
    local name  _cart
    local save  _cart
    local post _cartx
  }
  else {
    local save `name'
    local post `name'x
  }
  global name `name'
  if "`fail'"=="" {
     noisily display "fail(var) must be specified!"
     restore
     exit
  }
  unabbrev "`fail'"
  local fail "$S_1"
  drop if `fail'==.
  replace `fail'=`fail'>0.001    /* explicit definition of failures */
  local ffail "f(`fail')"

  if "`time'"=="" {
     noisily display "time(var) must be specified!"
     restore
     exit
  }
  unabbrev "`time'"
  local time "$S_1"
  drop if `time'==.|`time'<0
  replace `time'=max(.001,`time')   /* no zero times */
  local ttime "t(`time')"
  ttitle `time'
  local title1 $S_1
  keep `id' `adjust' `strata' `time' `fail' `varlist' `sumby' `tabby'
  *drop records with missing values
  parse "`varlist'",parse(" ")
  local vars `varlist'
  local labdef label define varnr
  local i=1
  while "`1'"~="" {
    drop if `1'==.    /* listwise deletion */
    global VAR`i' `1'
    local varlab:  var label `1'
    if "`varlab'"=="" {local varlab `1' }
    local labdef  "`labdef' `i'"`varlab'""
    local vlab`i' : value label `1'
    local vlab `vlab' `vlab`i''
    local i=`i'+1
    macro shift
  }
  if "`vlab'"~="" { label save `vlab' using vlab ,replace }
  if "`strata'"~="" {
    unabbrev `strata'
    local strata $S_1
    parse "`strata'",parse(" ")
    while "`1'"~="" {
      drop if `1'==.
      macro shift
    }
    local titstr `strata'
    local strata strata(`strata')
  }
  if "`adjust'"~="" {
    unabbrev `adjust'
    local adjust $S_1
    parse "`adjust'",parse(" ")
    while "`1'"~="" {
      drop if `1'==.
      macro shift
    }
    local titadj `adjust'
    local adjust adjust(`adjust')
  }

  /* Two new variables will be generated with fixed name:
     _h_     will contain the exp(b'x) for each subject for adjusted (stratified) Cox model.
             If no adjustment is done b=0 and _h_=1.
     _tmin_  used for reverse order of time. required for the calculation of the expected
                           number of events per subject within a CART group
  */
  gen _tmin_=-`time'
  if "`adjust'"!="" {
    stset `time' `fail'
* nois di "`titadj'"
* pause
    qui  stcox `titadj' ,`strata'
    predict _h_ , hr
  }
  else {    gen byte _h_=1   }
  global NVAR =`i'-1
  global post `post'
  global RHR1=1
  postfile $post group level from order varnr value minval maxval dec chi2 pnom padj n f e rhr hr  using `post',replace
  local minsize =max(1,`minsize')
  local minfail =max(1,`minfail')
  gen byte _CART=1
  label var _CART "CART group number"
  gen  _Corder=10
  label var _Corder "CART order of group"

  tempvar nf atr n1  dx
  global CARTi 1      /* index for current CART group */
  global CARTn 1      /* current number of CART groups */
  global CARTl 0      /* to be used as level of split */
  global RHR1=1
  while $CARTi<=$CARTn {
    nois display "CARTi=$CARTi  , CARTn=$CARTn "
    pause
    noisily carts `vars' , `ffail' `ttime' pval(`pval') `pnominal' minsize(`minsize') minfail(`minfail') `strata' `adjust'
    pause
  }
  save `save',replace
  postclose $post

  if $CARTn>1 {
    uc $post
    sort order
    count if order<=10
    local node1= _result(1)+1
    char _dta[node1]  `node1'   /* stores location of first split node */
    char _dta[cmd]    `cmd'     /* stores command and characteristics */
    char _dta[time]   `title1'
    char _dta[vars]   `vars'
    char _dta[strata] `titstr'
    char _dta[adjust] `titadj'

    replace order=_n + (order>10)
    sort group level
    quietly by group:gen byte last=(_n==_N)  /* end groups if last==1 */
    format chi2 hr rhr %6.2f
    label var group "Group"
    label var level "Level of split"
    label var order "Order branch/node"
    label var last  "End group indicator"
    label var from  "Split from group:"
    `labdef'
    change varnr . varnr  Split on variable
    label var value "Split cutpoint "
    label var minval "Minimum of range after split"
    label var maxval "Maximum of range after split"
    label var dec    "Number of decimals for value"
    label var chi2   "Chisquare value"
    label var pnom   "Nominal P-value"
    label var padj   "Adjusted P-value min Miller & Siegmund, Worsley"
    label var n "Total Number in group"
    label var f "Number of failures"
    label var e  "Expected # failures | upper level"
    label var rhr "Relative hazard rate [overall]"
    label var hr "Hazard Ratio | upper level"
    change minval min
    change maxval max
    change value val

    gen str20 mima =""  /* gen string variable that contains description of the split*/
    if "`vlab'"~="" { do vlab }
    local i=1
    while `i'<=_N {
      local j=varnr[`i']
      local min=min[`i']
      local max=max[`i']
      local dec= dec[`i']
      if "`vlab`j''"~="" {
        local mi:label `vlab`j'' `min'
        if `max'>`min' {local ma:label `vlab`j'' `max' }
      }
      else {
        fns mi `min' ,dec(`dec')
        local mi $mi
        if `max'>`min' {
          fns ma  `max',dec(`dec')
          local ma $ma
        }
      }
      if `max' >`min' { replace mima="`mi'-`ma'" in `i' }
      else { replace mima="`mi'" in `i' }
      local i=`i'+1
    }
    format rhr hr %4.2f
    format e  %4.1f
    eformat _all,vn val lab
    save,replace
*    relog  `name' ,240
    sort level group
    nois dis _n "CART history:"
    nois  qlist group level from order varnr val mima dec chi2 pnom padj n f e hr rhr  last,noo nodis
    nois dis _n "CART end groups:"
    nois  qlist group level from order varnr val mima dec chi2 pnom padj n f e hr rhr  if last==1,noo nodis

    uc `save'
*    keep `vars' `time' `fail' _EXP _CART
*    save,replace
    nois dis _n "Ranges of variables:" _n
    if $CARTn <=15 {
      nois sumby `vars' ,by(_CART) hea wi(2)  st(n min max)
    }
    else {
      nois sumby `vars'  if _CART<=17 ,by(_CART) hea wi(2)   st(n min max) nobytot
      if $CARTn> 17 {
        nois sumby `vars'  if _CART>17&_CART<=34 ,by(_CART) hea wi(2)   st(n min max) nobytot
      }
      if $CARTn> 34 {
        nois sumby `vars'  if _CART>34 ,by(_CART) hea wi(2)   st(n min max) nobytot
      }
    }
    nois dis _n
    if "`at'"!="" {
      capture nois sby _CART ,time(`time') dead(`fail')   st(nr f med se at) at(`at') hea
    }

    if "`sumby'"~="" {
      nois dis _n "Summary of some variables:" _n
      if $CARTn <=15 { nois sumby `sumby'  ,by(_CART) hea wi(2) }
      else {
        nois sumby `sumby' if _CART<=17  ,by(_CART) hea wi(2) nobytot
        if $CARTn> 17 {
          nois sumby `sumby' if _CART>17&_CART<=34 ,by(_CART) hea wi(2) nobytot
        }
        if $CARTn> 34 {
          nois sumby `sumby' if _CART>34 ,by(_CART) hea wi(2)   nobytot
        }
      }
    }

    if "`tabby'"~="" {
      nois dis _n "Tabulation of some variables:" _n
      if $CARTn <=15 { nois tabby `tabby'  ,by(_CART) hea wi(2) nozero }
      else {
        nois tabby `tabby' if _CART<=17  ,by(_CART) hea wi(2) nobytot  nozero
        if $CARTn> 17 {
          nois tabby `tabby' if _CART>17&_CART<=34 ,by(_CART) hea wi(2) nozero nobytot
        }
        if $CARTn> 34 {
          nois tabby `tabby' if _CART>34 ,by(_CART) hea wi(2)   nobytot  nozero
        }
      }
    }
*    cllog
    tree file($post) save(`name',replace) pval(`pval') `pnominal'
  }
  else {noisily display "No split possible"}
}

end

*!carts.ado Subroutine called by CART
*!CART analysis performs the searching within a group for optimal cutpoint and best covariate
/*
Syntax
  carts varlist , fail() time() pval() minsize() minfail() pnominal strata() adjust()

  Option adjust() for adjustment for other covariates is in fact not (yet) used!
  It is available for option to recalculate _h_ by refitting
  the Cox model within each CART group - if such an adjustment is needed.
*/
program define carts
version 7.0
local varlist "req ex min(1) "
local options "Pval(real 0.05) PNOMinal MINSize(int 10) MINFail(int 10)  Time(string) Fail(string) STRata(string) ADJust(string) "
parse "`*'"
quietly {
  tempfile temp1 temp2
  tempvar nf atr cum  sf se soe c tij n1 n2 totf
  * all calculations restricted to current group: _CART= $CARTi
  save "`temp1'",replace
  keep if _CART==$CARTi
  nois dis "Group $CARTi,Level $CARTl ,Number " _N

  if "`strata'"!="" {
    local qbystrata qby `strata' :
    local bystrata by(`strata')
  }

  egen `totf'=sum(`fail') , `bystrata'

  if _N<2*`minsize'|`totf'<`minfail' {
    global CARTi=$CARTi +1
    global CARTl ${L$CARTi}
    use "`temp1'",clear
    exit
  }

/* May 20, 2002 I come to the conclusion that the approach followed below is too simplistic:
   to determine the expected values by multiplication by a factor.
   Another and better way has been implemented.
  *Adapt the _EXP values such that they sum to `totf'
  su _EXP
  local factor=`totf'/(_result(1)*_result(3))
  gen `cum'=_EXP*`factor'
*/
  tempvar H H1 exp
  sort `strata' _tmin_
  `qbystrata' gen `H'=sum(_h_)
  qby `strata' _tmin_:replace `H'=`H'[_N]
  sort `strata' `time'
  `qbystrata' gen `H1'=sum(`fail'/`H')
  qby `strata' `time' :replace `H1'=`H1'[_N]
  gen `exp'= _h_*`H1'  /* expected number of failures within this CART group for this subject */

  local cc=0
  local pp=1
  local i=1
  while   `i'<=$NVAR {
    local x ${VAR`i'}
    sort `x' `time' `fail'
    count if `x'>`x'[`minsize']
    if _result(1)>=`minsize' {  /*else not splittable on var i */

      if "`strata'"=="" {
        gen `sf'=sum(`fail')
        gen `se'=sum(`exp')
        quietly by `x':replace `sf'=. if _n<_N
        quietly by `x':replace `se'=. if _n<_N
        gen `soe'=`sf'-`se'
        gen `c'=(`soe'^2)/(`se'*(`totf'-`se')/`totf') if `totf'>0 &`se'>0& `totf'-`se'>0
        drop `sf' `se' `soe'
      }
      else {
        sort `strata' `x' `time' `fail'
        `qbystrata' gen `se'=sum(`exp')
        quietly by `strata' `x':replace `se'=`se'[_N]
        *first we calculate the variances within strata and the differences for successive cutpoints
        tempvar Var dVar
        gen `Var'=(`se'*(`totf'-`se')/`totf')             if `totf'>0 &`se'>0& `totf'-`se'>0
        replace `Var'=0 if `Var'==.
        `qbystrata' gen  `dVar'=`Var'-max(0,`Var'[_n-1])   /* we need the differences! */
        sort `x' `time'  `strata' `fail'
        gen `soe'=sum(`fail' - `exp')
        qby `x': replace `soe'=. if _n<_N
        replace `Var'=sum(`dVar')                    /* now we have the variance for a cutpoint */
        gen   `c'=(`soe'^2)/`Var'     if `Var'>0       /* now turned into Chi2 value */
        drop  `se' `soe' `Var' `dVar'
      }

      su `c' if `x'>=`x'[`minsize']&`x'<`x'[_N-`minsize']
      local ncutp=r(N)
      local max=_result(6)
      local pnom=chiprob(1,`max')  /* unadjusted P-value */

* nois display "Max chi square= `max', P nominal= `pnom'"
* nois display "pnominal `pnominal'"
*pause on
*nois su `sf' `se' `soe' `c'
*pause
*pause off
      /* now we have to adjust for the multiple testing. The adjustment is optional via the otion
         "pnominal". The adjustment will be done in the same way as in SRD.ADO.
         Adjustment will not be done if the nominal P-value is already larger than required `pval'
         or larger than the smaller p-value on another factor. It will only be done if there is
         really multiple testing, i.e. at least 2 potential cutpoints

         Note, that even with the option pnominal we calculate the adjusted P-value. Only it will not be
         used for the decisions.
      */
      if `pnom'<min(`pp',`pval')  {
        if `ncutp'>1 {
          *  adjustment Formula Altman Miller&Siegmund
          * Calculate fractions eps1 and eps2 for Miller&Siegmund correction
          count if `x'<=`x'[`minsize']
          local eps1=r(N)/_N
          count if `x'>=`x'[_N-`minsize']
          local eps2=r(N)/_N
          local pcor1=normd(`max'^.5)*(4+(`max'-1)*log((1-`eps1')*(1-`eps2')/(`eps1'*`eps2')))/`max'^.5

          *   adjustment improved Bonferroni (Worsley)
          local totn=_N
          qby `x': gen int `n1'=_N if _n==1
          gen int `n2'=sum(`n1')
          qby `x': replace `n2'=`n2'[_N]
          replace `n2'=. if `n2'<`minsize'|`n2'>_N-`minsize'
          su `n2'
          replace `n2'=. if `n2'>=r(max)-.0001
          gen `tij'= (`totn'*`n1'[_n+1]/(`n2'[_n+1]*(`totn'-`n2')))^0.5
          replace `tij'= (`tij'-(`max'/4-1)*`tij'^3/6) if `tij'!=.
          replace `tij'=sum(`tij')
          local pcor2=chiprob(1,`max')+(1/_pi)*exp(-`max'/2)*`tij'[_N]
          ** The result is sometimes negative! This implies an error somewhere.
          ** For the time being pcor2 will be ignored if pcor2<pnom
          if `pcor2'<`pnom' {local pcor2 . }
          local pcor=min(`pcor1',`pcor2')
          drop `n1' `n2' `tij'
        }
        else { local pcor `pnom' }
* nois displ "pcor `pcor', pp `pp', pval `pval', `pnominal'"
        if (`pcor'<`pp'&`pcor'<`pval') | "`pnominal'"!="" {
          su `x' if abs(`c'-`max')<.000001
          local value =_result(5)  /* if more cutpoints with the same max Chi square value the smallest is chosen */
          if `value'<0 {local minval0=0.99999999*`value' }
          else {local minval0=1.00000001*`value'}
          local var    `i'
          local varnam `x'
          if `ncutp'>1 {local ppadj `pcor'}
          else         {local ppadj .     }
          local cc `max'
          local ppnom `pnom'
          local pp =min(`pnom',`ppadj')
          if "`pnominal'"!="" { local pp `pnom' }
        }
      }
      drop `c'
    }
    local i=`i'+1
  }

  if `pp'<`pval' {
    global CARTn   =$CARTn +1
    global CARTl   =$CARTl +1
    global L$CARTn = $CARTl
    local x `varnam'

    su `exp' if `x'<=`minval0'
    local e=_result(1)*_result(3)
    count if `x'<=`minval0' &`fail'==1
    local f=_result(1)
    count if  `x'<=`minval0'
    local nn=_result(1)
*    nois displ "x1 $x1 , nn `nn', N " _N
*    pause
    capture drop `totf'
    count if `fail'==1
    local totf=r(N)
    local hr=`f'*(`totf'-`e')/(`e'*(`totf'-`f'))

    su `x' if `x'<=`minval0'
    local minval=_result(5)
    su `x' if `x'>`minval0'
    local minval2=_result(5)
    local maxval=_result(6)
    local format:format `x'
    local type: type `x'
    if "`type'"!="float" &"`type'"!="double"  {
      local minval=int(`minval')
      local minval2=int(`minval2')
      local maxval=int(`maxval')
    }
    local dec =substr("`format'",index("`format'",".")+1,1)
    local order=_Corder[1]

    if `hr'<1 {
      global RHR$CARTn= ${RHR$CARTi}*(`totf'-`f')/(`totf'-`e')
      global RHR$CARTi= ${RHR$CARTi}*`f'/`e'
      replace _CART=$CARTn if `x'>`minval0'
      replace _Corder=_Corder+.5^$CARTl if `x'>`minval0'
      replace _Corder=_Corder-.5^$CARTl if `x'<=`minval0'
      post $post ($CARTi) ($CARTl) ($CARTi) (`order'-.5^$CARTl) (`var') (`value') (`minval')  (`value')  (`dec') (`cc') (`ppnom') (`ppadj') (`nn') (`f') (`e')  (${RHR$CARTi}) (`hr')
      post $post ($CARTn) ($CARTl) ($CARTi) (`order'+.5^$CARTl) (`var') (`value') (`minval2') (`maxval') (`dec') (`cc') (`ppnom') (`ppadj') (_N-`nn') (`totf'-`f') (`totf'-`e') (${RHR$CARTn}) (1/`hr')
    }
    else      {
      global RHR$CARTn= ${RHR$CARTi}*`f'/`e'
      global RHR$CARTi= ${RHR$CARTi}*(`totf'-`f')/(`totf'-`e')
      replace _CART=$CARTn if `x'<=`minval0'
      replace _Corder=_Corder-.5^$CARTl if `x'>`minval0'
      replace _Corder=_Corder+.5^$CARTl if `x'<=`minval0'
      post $post ($CARTi) ($CARTl) ($CARTi) (`order'-.5^$CARTl) (`var') (`value') (`minval2') (`maxval') (`dec') (`cc') (`ppnom') (`ppadj') (_N-`nn') (`totf'-`f') (`totf'-`e') (${RHR$CARTi}) (1/`hr')
      post $post ($CARTn) ($CARTl) ($CARTi) (`order'+.5^$CARTl) (`var') (`value') (`minval')  (`value')  (`dec') (`cc') (`ppnom') (`ppadj') (`nn') (`f') (`e') (${RHR$CARTn}) (`hr')
    }
    save "`temp2'",replace
    use "`temp1'",clear

    replace _CART=0  if  _CART==$CARTi
    append using "`temp2'"
    drop if  _CART==0
    save $name,replace
  }
  else { /* no split of this group */
    global CARTi=$CARTi +1
    global CARTl ${L$CARTi}
    use "`temp1'",clear
    save $name,replace
  }
}
end



*!tree.ado   Subroutine for CART.
*!10/4/98 -> 5/8/98
/*
Generates the tree and necessary variables starting with the data in the CART history file.

Syntax :
  tree file(filename of posted CART history file) save(string) pnominal

filename is the name of the CART history file

A tree consists of split-nodes and branches.

The split-nodes are also numberd in order of generation starting with 1 as
	top-node, i.e. the first split. Per node there are in general
	two branches. However, we will try to take into account future gene-
	ralisations with more than two branches per node. For example at the
	start of the tree, if one forces a split in three categories on a
	specified variable.
	A node with n>2 branches is thought of as cocnsiting of n-1 subnodes
	at the same level.

*/

program define tree
capture preserve
local options "File(string) Save(string) PNOMinal Pval(real 0.05) "
parse ",`*'"
uc `file'

quietly {
 *Create help variables for location and position of nodes and branches */

 capture confirm new var next
 if _rc==0 {  /* otherwise all help variables are already defined */
  sort group level
  by group:gen next=varnr[_n+1]  /* Var nr of next split variable, if not end-node */
 }

 if "`pnominal'"!="" { local t1note Split if nominal P<`pval' }
 else                { local t1note Split if (adjusted) P<`pval' }

 su level
 local nlev=_result(6)
 local nbr =_N+1
 cvar mima
 local lmima:char mima[width]
 local dy=int(20000/`nbr')

 local fh1=int(min(`dy'/2,700))  /* font group descriptor text */
* local fh1=int(min(max(`dy'/2,400),700))  /* font group descriptor text */
 local fw1=int(`fh1'/2)
 local fh2=int(1.1*`fh1')   /* font split variable text (standard font) */
 local fw2=int(1.1*`fw1')

* local fh2=int(max(1.1*`fh1',300))   /* font split variable text (standard font) */
* local fw2=int(`fh2'/2)
 local fh1x=int(`fh2'/1.1)
 local fw1x=int(`fh1x'/2)

 local dx=int(max(10,`lmima')*1.7*`fw1')
 if (`nlev'+.2)*`dx'+14*1.7*`fw2' >30000 {
	local f=30000/((`nlev'+.2)*`dx'+14*1.7*`fw2')
*	local f=90000/((`nlev'+.2)*`dx'+14*1.7*`fw2')
	local fh1=int(`fh1'*`f')
	local fw1=int(`fw1'*`f')
	local dx=int(`dx'*`f')
 }
 local font1 `fh1' `fw1'
* local font1 `fh1x' `fw1x'
 local font2 `fh2' `fw2'
 local ax=500
 local ay=2000

 local nc1=`ax'+`dx'*(`nlev'+.2)+3*1.6*`fw2'
 local nc2=`ax'+`dx'*(`nlev'+.2)+7*1.6*`fw2'
 local nc3=`ax'+`dx'*(`nlev'+.2)+12*1.6*`fw2'
 local rtop=`ay'+`dy'-2*`fh2'
 local node1 :char _dta[node1]
 local r1=`ay'+`node1'*`dy'
 local c1=`ax'
 sort group level
 local varn=varnr[1]
 local text1:label varnr `varn'

*All elements for a tree are available. Let us make one.
/*
nois display "r1 `r1'"
nois display "c1 `c1'"
nois display "nc1 `nc1'"
nois display "nc2 `nc2'"
nois display "nc3 `nc3'"
nois display "dx `dx'"
nois display "dy `dy'"
nois display "ax `ax'"
nois display "ay `ay'"
nois display "font1 `font1'"
nois display "font2 `font2'"
nois display "text1 `text1'"
pause
*/
  gph open  ,saving(`save')
  gph pen 1
  gph font `font2'
  gph text `rtop' `nc1'  0 1 N
  gph text `rtop' `nc2'  0 1 F
  gph text `rtop' `nc3'  0 1 RHR

  gph point  `r1' `c1' `fw1' 1
  local c1=`c1'+30
  gph pen 2
  gph font `font2'
  gph text `r1' `c1'  0 -1 `text1'

  sort from level order
  gph pen 1
  local i=1
  while `i'<=_N {
    if from[`i']==from[`i'-1]&level[`i']==level[`i'-1] {
	local x=`ax'+`dx'*(level[`i']-1)
	local y1=`ay'+`dy'*(order[`i'-1])
	local y2=`ay'+`dy'*(order[`i'])
	gph line `y1' `x' `y2' `x'
    }
    local i=`i'+1
  }

  sort order
  local i=1
  while `i'<=_N {
   	local last=last[`i']
	local next=next[`i']
	local text1=mima[`i']

	local x1=`ax'+`dx'*(level[`i']-1)
	local x2=`ax'+`dx'*(level[`i'])
	local x3=`x2'+30
	local y1 =`ay'+`dy'*order[`i']
	local y2=`y1'-130
        gph pen 1
	gph line `y1' `x1' `y1' `x2'
	gph pen 1
	gph font `font1'
	gph text `y2' `x1'  0 -1 `text1'
	if "`next'"~="." {
		local text2:label varnr `next'
		gph point `y1' `x2' `fw1' 1
		gph pen 2
		gph font `font2'
		gph text `y1' `x3'  0 -1 `text2'
	}
	else {
		gph point `y1' `x2' `fw1' 2
		local no=n[`i']
		local nf=f[`i']
		local rhr=rhr[`i']
		fns rhr `rhr' ,dec(2)
		gph font `font2'
		gph text `y1' `nc1'  0 1 `no'
		gph text `y1' `nc2'  0 1 `nf'
		gph text `y1' `nc3'  0 1 $rhr
	}
	local i=`i'+1
  }
  local t1 :char _dta[time]
  local t2 :char _dta[vars]
  local t3 :char _dta[adjust]
  local t4 :char _dta[strata]
  local nx=("`t3'"~="")+("`t4'"~="")
  local yt=1.2*(`fh2'+(1+`nx')*`fh1')
  local f=min(1,`ay'/`yt')
  local y1=1.2*`f'*`fh2'
  local y2=`y1'+1.2*`fh1'
  local y3=`y2'+1.2*`fh1'
  local y4=`y3'+1.2*`fh1'*(`nx'==2)
  local fh1=`f'*`fh1'
  local fw1=`f'*`fw1'
  local fh2=`f'*`fh2'
  local fw2=`f'*`fw2'
  gph pen 2
  gph font `fh2' `fw2'
  gph text `y1' `ax' 0 -1 CART analysis `t1' - `t1note'
  gph pen 1
  gph font `fh1' `fw1'
  gph text `y2' `ax' 0 -1 With variables: `t2'
  if  "`t3'"~="" {  gph text `y3' `ax' 0 -1 Adjusted for: `t3'}
  if  "`t4'"~="" {  gph text `y4' `ax' 0 -1 Stratified by: `t4'}
  gph close


end













