*! version 2.0 dgc/mh Sep 1997  STB-41 ssa10.1

program define staalen
    version 5.0
    st_is
    local if "opt"
    local in "opt"
    local options "BY(string) FRom(real -1.e20) TO(real 1.e+20) FCodes(string)"
    local options "`options' RISK RMult(string) ORigin(string) SCale(real 1)"
    local options "`options' L1(string) B2(string) ylog * "
    parse "`*'"
    * Extract st data

    local tout : char _dta[st_t]
    local ttype double
    local tin : char _dta[st_t0]
    local fail : char _dta[st_d]
    local wv : char _dta[st_wv]
    local wt : char _dta[st_wt]
    local aw : char _dta[st_w]

    * Checks and error messages

    if "`aw'"!="" {
        di in re "weighted analysis not implemented yet"
        exit
    }
    if "`by'"!="" {
        local nby : word count `by' 
        if `nby' != 1 {
            di in re "Only one variable allowed in by() option"
            exit
        }
        confirm var `by'
    }
    if "`origin'"!="" {
       capture confirm var `origin'
       if _rc==0 {
          local vor "`origin'"
       }
       else {
          capture confirm number `origin'
          if _rc!=0 {
             di in re "illegal origin option"
             exit
          }
       }
    }
    * Show st data characteristics 

    di in gr "Nelson-Aaalen estimate"
    st_show
    if "`fcodes'"!="" {
       di
       di in gr "   failure codes:  " in ye "`fcodes'"
    }
    if "`origin'"!="" {
       di in gr "     time origin:  " in ye "`origin'"
    }
    if `scale' !=1 {
       di in gr "    time scaling:  " in ye "`scale'"
    }
    if "`rmult'"!="" {
       di in gr  " rt multipliers:  " in ye "`rmult'"
    }
    else {
       tempvar rmult
       gen `rmult' = 1
    }
   
    * x-axis label

    if "`b2'"=="" {
        local b2 "Time"
    }

    * y-axis label

    if "`l1'"=="" {
        if `rmult'==1 {
	    if "`risk'"=="" {local l1 "Cumulative rate"}
	    else {local l1 "Cumulative risk"}
	}
	else {
	    if "`risk'"=="" {local l1 "Cumulative baseline rate"}
	    else {local l1 "Cumulative baseline risk"}
	}
    }


    preserve
    tempvar touse timin timout
    gen double `timout' = `tout'
    if "`tin'"=="" {
       gen double `timin' = 0.0
    }
    else {
       gen double `timin' = `tin'
    }
    mark `touse' `aw' `if' `in'
    markout `touse' `vor' `timin' `timout' `fail' `by' `rmult'
    quietly {
        keep if `touse'
        keep `vor' `timin' `timout' `fail' `by' `rmult'
    }
    if "`wv'"=="" {
        tempvar wv
        gen `wv' = 1
    }

    quietly {
        if "`origin'"!="" {
            replace `timin' = `timin' - (`origin')
            replace `timout' = `timout' - (`origin')
        }
        if `scale'!=1 {
            replace `timin' = `timin'/`scale'
            replace `timout' = `timout'/`scale'
        }
	tempvar id t cum tcf grp ntry
        if "`fcodes'" == "" {
            replace `fail' = (`fail'!=0)
        }
        else {
            recode `fail' `fcodes' = 1  * = 0
        }
	gen long `id'=_n
	drop if `timin'>`to'|`timout'<`from'
	replace `timin'=`from' if `timin'<`from'
	replace `timout'=`to' if `timout'>`to'
	replace `fail'=0 if `timout'==`to'
	expand 2
	sort `id'
	by `id': gen `ntry' = (_n==1)
	gen `ttype' `t'= cond(`ntry', `timin', `timout')
	/* In case of tied times, order is failures<censorings<entries */
	gen int `tcf' = cond(`ntry', 2, 1-`fail')
	replace `fail'=0 if `ntry'
	replace `rmult' = -`rmult' if !`ntry'
	sort `by' `tcf' `t'
	by `by' `tcf' `t': replace `rmult' = `rmult'+`rmult'[_n-1] if _n>1 
	by `by' `tcf' `t': replace `fail' = `fail'+`fail'[_n-1] if _n>1
	by `by' `tcf' `t': drop if _n!=_N
	sort `by' `t' `tcf'
	if "`by'"=="" {
	    replace `rmult' = sum(`rmult')
	    gen `cum' = cond(_n==1, 0, `fail'/`rmult'[_n-1])
	    keep if `cum'>0|_n==1|_n==_N
	    replace `cum' = sum(`cum')
	    if "`risk'"!="" {replace `cum' = 1 - exp(-`cum')}
	    if "`ylog'"!="" {replace `cum' = . if `cum' <= 0}
	    local vlist "`cum'"
	    local syms "i"
	    local connect "J"
	    local pen "2"
	}
	else {
	    by `by': replace `rmult' = sum(`rmult')
	    by `by': gen `cum' = cond(_n==1, 0, `fail'/`rmult'[_n-1])
	    by `by': keep if `cum'>0|_n==1|_n==_N
	    by `by': replace `cum' = sum(`cum')
	    by `by': gen long `grp' = (_n==1)
	    replace `grp' = sum(`grp')
	    local group 1
	    while `group' <= `grp'[_N] {
		tempvar gsc
		gen `gsc' = `cum' if `group'==`grp'
		if "`risk'"!="" {replace `gsc' = 1 - exp(-`gsc')}
	        if "`ylog'"!="" {replace `gsc' = . if `gsc' <= 0}
		local vlist "`vlist' `gsc'"
		local syms "`syms'i"
		local connect "`connect'J"
		local group = `group'+1
		local pen "`pen'`group'"
	    }
	    tempvar gsc
	    by `by': gen `gsc' = `cum' if _n==_N
	    if "`risk'"!="" {replace `gsc' = 1 - exp(-`gsc')}
	    if "`ylog'"!="" {replace `gsc' = . if `gsc' <= 0}
	    local vlist "`vlist' `gsc'"
	    local syms "`syms'[`by']"
	    local connect "`connect'."
	    local pen "`pen'1"
	}
    }
    graph `vlist' `t', symbol(`syms') connect(`connect') /*
	*/ l1(`l1') b2(`b2') pen (`pen') `ylog' `options'
    end