*! version 0.1 HS, Feb 20, 2004 *! version 0.2 HS, Jun 15, 2006 /* wtdpred Computes predicted density for waiting time distribution. As default the latest estimated wtdml-model is used, i.e. when no options are given. If, however, just one option is given, they must all be set. If you supply options manually, be aware that the program does little with respect to checking consistency of the supplied options - use at your own risk! */ program define wtdpred version 8.2 syntax varlist(min = 2 max = 2) [if] [in] [, prevd(string) hdens(string) ddens(string) cens(string) coefvek(string) tailar(real 0)] tokenize `varlist' local time `1' local pred `2' tempvar tmptime tmppred stdtime hinc S0mort S1mort tempfile orgdat tempname pval D1 D0 H phi p10 p11 curar gparm /* If there are no options, we use last estimated wtdml model */ if "`options'" == "" { foreach parm in prevd hdens ddens cens { local `parm' = e(`parm') } local tailar = e(tailar) tempname coefvek matrix `coefvek' = e(b) } marksample touse, novar qui { gen double `stdtime' = `time' if `touse' scalar `pval' = exp(`coefvek'[1, 1]) / (1 + exp(`coefvek'[1, 1])) tempvar gpred gen `gpred' = . if "`prevd'" == "exp" { local coeflg = colsof(`coefvek') scalar `gparm' = exp(`coefvek'[1, `coeflg']) replace `gpred' = `gparm' * exp( - (`gparm' * `stdtime')) } if "`prevd'" == "wei" { local coeflg = colsof(`coefvek') local coeflg1 = `coeflg' - 1 matrix `gparm' = `coefvek'[1, `coeflg1'..`coeflg'] frwei `stdtime' `gpred', gparm("`gparm'") } if "`prevd'" == "lnorm" { local coeflg = colsof(`coefvek') local coeflg1 = `coeflg' - 1 matrix `gparm' = `coefvek'[1, `coeflg1'..`coeflg'] frlnorm `stdtime' `gpred', gparm("`gparm'") } if "`hdens'" == "unif" { scalar `H' = exp(`coefvek'[1, 2]) / (1 + exp(`coefvek'[1, 2])) gen `hinc' = 1 } if "`hdens'" == "exp" { scalar `H' = exp(- exp(`coefvek'[1, 2])) gen `hinc' = . tempname Hrate scalar `Hrate' = exp(`coefvek'[1, 2]) cf_exp `stdtime' `hinc', beta(`Hrate') } if "`cens'" == "none" { replace `pred' = (`pval' * `gpred' + (1 - `pval') * (1 - `H') * `hinc') /* */ / (`pval' * (1 - `tailar') + (1 - `pval') * (1 - `H')) } if index("`cens'", "dep") ~= 0 { if "`ddens'" == "" { di in r "ddens must be specified when cens() is set to `cens'" exit } sa `orgdat' if "`prevd'" == "exp" { drop _all set obs 1000 gen double `tmptime' = (_n - .5) / _N gen double `tmppred' = . replace `tmppred' = `gparm' * exp( - (`gparm' * `tmptime')) } if "`prevd'" == "wei" { drop _all set obs 1000 gen double `tmptime' = (_n - .5) / _N gen double `tmppred' = . frwei `tmptime' `tmppred', gparm("`gparm'") } if "`prevd'" == "lnorm" { drop _all set obs 1000 gen double `tmptime' = (_n - .5) / _N gen double `tmppred' = . frlnorm `tmptime' `tmppred', gparm("`gparm'") } if "`cens'" == "depphi" { local philoc = colnumb(`coefvek', "lnphi:_cons") if `philoc' == . { di in r "Required parameter -lnphi- is not supplied" exit } else { scalar `phi' = exp(`coefvek'[1, `philoc']) } } else { local philoc = colnumb(`coefvek', "lnphi:_cons") if `philoc' ~= . { di in r "The supplied coefficent vector contains a lnphi parameter" di in r "but you did not use option -cens(depphi)-" exit } scalar `phi' = 1 } if "`hdens'" == "unif" { gen `hinc' = 1 } if "`hdens'" == "exp" { gen `hinc' = . cf_exp `tmptime' `hinc', beta(`Hrate') } if "`ddens'" == "unif" { if index("`cens'", "dep") == 1 { scalar `D1' = exp(`coefvek'[1, 3]) / (1 + exp(`coefvek'[1, 3])) scalar `D0' = exp(`coefvek'[1, 4]) / (1 + exp(`coefvek'[1, 4])) } if index("`cens'", "indep") == 1 { scalar `D1' = exp(`coefvek'[1, 3]) / (1 + exp(`coefvek'[1, 3])) scalar `D0' = `D1' } gen `S1mort' = 1 - `tmptime' gen `S0mort' = 1 - `tmptime' } if "$hmodel" == "exp" { if index("`cens'", "dep") == 1 { scalar `D1' = exp(- exp(`coefvek'[1, 3])) scalar `D0' = exp(- exp(`coefvek'[1, 4])) } if index("`cens'", "indep") == 1 { scalar `D1' = exp(- exp(`coefvek'[1, 3])) scalar `D0' = `D1' } gen `S1mort' = . gen `S0mort' = . forv d = 0/1 { tempname Drate`d' scalar `Drate`d'' = - log(`D`d'') cs_exp `tmptime' `S`d'mort', beta(`Drate`d'') } } scalar `p11' = (1 - `H') * (1 - `D0') * `phi' scalar `p10' = (1 - `H') - `p11' replace `tmppred' = `pval' * `tmppred' * (`D1' + /* */ (1 - `D1') * `S1mort') + /* */ (1 - `pval') * `hinc' * (`p10' + `p11' * `S0mort') su `tmppred', mean scalar `curar' = r(mean) scalar curar = `curar' use `orgdat', clear if "`ddens'" == "unif" { gen `S1mort' = 1 - `stdtime' gen `S0mort' = 1 - `stdtime' } if "`ddens'" == "exp" { gen `S1mort' = . gen `S0mort' = . forv d = 0/1 { cs_exp `stdtime' `S`d'mort', beta(`Drate`d'') } } replace `pred' = `pval' * `gpred' * (`D1' + /* */ (1 - `D1') * `S1mort') + /* */ (1 - `pval') * `hinc' * (`p10' + `p11' * `S0mort') replace `pred' = `pred' / `curar' } } end