*! version 1.2.0 Stephen P. Jenkins, December 2005 *! corrected bug in SEs displayed for Prob(Types), *! fixed help file reference to starting values for Prob(Types) *! version 1.1.1 Stephen P. Jenkins, September 2004 *! Discrete time proportional hazard model with Heckman-Singer mass point heterogeneity /*------------------------------------------------ playback request */ program define hshaz, eclass byable(onecall) sortpreserve version 8.2 if replay() { if "`e(cmd)'" != "hshaz" { noi di as error "results for hshaz not found" exit 301 } if _by() { error 190 } Display `0' exit `rc' } if _by() { by `_byvars'`_byrc0': Estimate `0' } else Estimate `0' end /*------------------------------------------------ estimation */ program define Estimate, eclass byable(recall) syntax varlist(numeric) [if] [in] [fw iw] [, /// Id(varname) Dead(varname numeric) Seq(varname numeric) /// Nmp(integer 2) M2(real 1) P2(real .3) M3(real -1) P3(real .3) /// M4(real .1) P4(real .1) M5(real -.1) P5(real .1) /// NOCONS noBeta0 /// Level(integer $S_level) EForm /// noLOG * ] local title "Discrete time PH model, with discrete mixture" local option0 `options' local wtype `weight' local wtexp `"`exp'"' if "`weight'" != "" { local wgt `"[`weight'`exp']"' } if "`level'" != "" { local level "level(`level')" } local log = cond("`log'" == "", "noisily", "quietly") local log2 = cond("`beta0'" == "", "noisily", "quietly") marksample touse markout `touse' `varlist' qui count if `touse' if r(N) == 0 { error 2000 } mlopts mlopts, `options' set more off if "`id'" == "" { di as error "Variable identifying each person must " di as error " be specified in id(idvar) " exit 198 } if "`dead'" == "" { di as error "Per-interval censoring indicator variable must " di as error "be specified in dead(deadvar):" di "see -hshaz- help file" exit 198 } capture assert `dead'==1 | `dead'==0 if _rc~=0 { di as error "Per-interval censoring indicator variable must " di as error "equal one or zero: see -hshaz- help file" exit 198 } if "`seq'" == "" { di as error "Integer-valued variable identifying " di as error "spell time interval be specified " _c di as error "in seq(seqvar) " exit 198 } capture assert (`nmp' == 2 | `nmp' == 3 | `nmp' == 4 | `nmp' == 5) if _rc ~= 0 { di as error "Number of mass points must be 2 (default), 3, 4, or 5" exit 198 } capture assert `m2' ~= 0 if _rc ~= 0 { di as error "Choose non-zero starting value for m2" exit 198 } capture assert (`p2' > 0 & `p2' < 1) if _rc ~= 0 { di as error "p2 must be between 0 and 1" exit 198 } if `nmp' == 2 { local m3 local p3 local m4 local p4 local m5 local p5 } if `nmp' == 3 { local m4 local p4 local m5 local p5 } if `nmp' == 4 { local m5 local p5 } if `nmp' > 2 { capture assert `m3' ~= 0 if _rc ~= 0 { di as error "Choose non-zero starting value for m3" exit 198 } capture assert (`p3' > 0 & `p3' < 1) if _rc ~= 0 { di as error "p3 must be between 0 and 1" exit 198 } } if `nmp' > 3 { capture assert `m4' ~= 0 if _rc ~= 0 { di as error "Choose non-zero starting value for m4" exit 198 } capture assert (`p4' > 0 & `p4' < 1) if _rc ~= 0 { di as error "p4 must be between 0 and 1" exit 198 } } if `nmp' > 4 { capture assert `m5' ~= 0 if _rc ~= 0 { di as error "Choose non-zero starting value for m5" exit 198 } capture assert (`p5' > 0 & `p5' < 1) if _rc ~= 0 { di as error "p5 must be between 0 and 1" exit 198 } } if `nmp' == 3 { capture assert (`p2' + `p3' < 1) if _rc ~= 0 { di as error "Must have p2 + p3 < 1" exit 198 } } if `nmp' == 4 { capture assert (`p2' + `p3' + `p4' < 1) if _rc ~= 0 { di as error "Must have p2 + p3 + p4 < 1" exit 198 } } if `nmp' == 5 { capture assert (`p2' + `p3' + `p4' + `p5' < 1) if _rc ~= 0 { di as error "Must have p2 + p3 + p4 + p5 < 1" exit 198 } } tempvar mysamp tempname b b0 b1 b2 lnf V mass2 logitp20 logitp2 /// mass3 logitp30 logitp3 mass4 logitp40 logitp4 /// mass5 logitp50 logitp5 scalar `logitp20' = logit(`p2') if `nmp' >= 3 { scalar `logitp30' = logit(`p3') } if `nmp' >= 4 { scalar `logitp40' = logit(`p4') } if `nmp' >= 5 { scalar `logitp50' = logit(`p5') } set more off _rmcoll `varlist' `wgt' if `touse', `nocons' local varlist `r(varlist)' /* get estimates of no-het model */ if "`beta0'" == "" { di as txt "Discrete time PH model without frailty" } `log2' glm `dead' `varlist' `wgt' if `touse', f(b) l(c) /// `level' `eform' `nocons' nolog quietly { local LL1 = e(ll) matrix `b0' = e(b) matrix coleq `b0' = hazard // now pack out -b0- with starting // value for mass point and logitp matrix `mass2' = (`m2') matrix colnames `mass2' = m2:_cons matrix `logitp2' = (`logitp20') matrix colnames `logitp2' = logitp2:_cons if `nmp' >= 3 { matrix `mass3' = (`m3') matrix colnames `mass3' = m3:_cons matrix `logitp3' = (`logitp30') matrix colnames `logitp3' = logitp3:_cons } if `nmp' >= 4 { matrix `mass4' = (`m4') matrix colnames `mass4' = m4:_cons matrix `logitp4' = (`logitp40') matrix colnames `logitp4' = logitp4:_cons } if `nmp' >= 5 { matrix `mass5' = (`m5') matrix colnames `mass5' = m5:_cons matrix `logitp5' = (`logitp50') matrix colnames `logitp5' = logitp5:_cons } if `nmp' == 2 { matrix `b1' = `b0',`mass2',`logitp2' local eqq "(m2: ) (logitp2: ) " } if `nmp' == 3 { matrix `b1' = `b0',`mass2',`mass3',`logitp2',`logitp3' local eqq "(m2: ) (m3: ) (logitp2: ) (logitp3: ) " } if `nmp' == 4 { matrix `b1' = `b0',`mass2',`mass3',`mass4',`logitp2',`logitp3',`logitp4' local eqq "(m2: ) (m3: ) (m4: ) (logitp2: ) (logitp3: ) (logitp4: ) " } if `nmp' == 5 { matrix `b1' = `b0',`mass2',`mass3',`mass4', `mass5',`logitp2',`logitp3',`logitp4',`logitp5' local eqq "(m2: ) (m3: ) (m4: ) (m5: ) (logitp2: ) (logitp3: ) (logitp4: ) (logitp5: ) " } /* estimate full model */ sort `id' `touse' `seq' global S_E_id "`id'" global S_E_dd "`dead'" global S_E_nmp "`nmp'" `log' ml model d0 hshaz_ll (hazard: `dead' = `varlist', `nocons') `eqq' /// `wgt' if `touse', maximize init(`b1') waldtest(0) search(off) /// collinear title(`title') `level' `mlopts' `stdopts' `modopts' tempname tag egen `tag' = tag(`id') if e(sample) count if `tag' & e(sample) eret scalar N_spell = r(N) eret local cmd "hshaz" eret local depvar "`dead'" eret local idvar "`id'" eret local seqvar "`seq'" eret scalar ll_nofr = `LL1' // prob(Type 2) = invlogit(logitp) // s.e.(prob) = s.e.(logitp) * invlogit(logitp) * invlogit(-logitp) eret scalar nmp = `nmp' eret scalar m1 = 0 eret scalar m2 = [m2]_cons eret scalar se_m2 = [m2]_se[_cons] eret scalar pr2 = invlogit( [logitp2]_cons ) eret scalar se_pr2 = [logitp2]_se[_cons] * invlogit( [logitp2]_cons) * invlogit( -[logitp2]_cons ) if `nmp' > 2 { eret scalar m3 = [m3]_cons eret scalar se_m3 = [m3]_se[_cons] eret scalar pr3 = invlogit( [logitp3]_cons ) eret scalar se_pr3 = [logitp3]_se[_cons] * invlogit( [logitp3]_cons) * invlogit( -[logitp3]_cons ) } if `nmp' > 3 { eret scalar m4 = [m4]_cons eret scalar se_m4 = [m4]_se[_cons] eret scalar pr4 = invlogit( [logitp4]_cons ) eret scalar se_pr4 = [logitp4]_se[_cons] * invlogit( [logitp4]_cons) * invlogit( -[logitp4]_cons ) } if `nmp' > 4 { eret scalar m5 = [m5]_cons eret scalar se_m5 = [m5]_se[_cons] eret scalar pr5 = invlogit( [logitp5]_cons ) eret scalar se_pr5 = [logitp5]_se[_cons] * invlogit( [logitp5]_cons) * invlogit( -[logitp5]_cons ) } if `nmp' == 2 { eret scalar pr1 = 1 - `e(pr2)' } if `nmp' == 3 { eret scalar pr1 = 1 - `e(pr2)' - `e(pr3)' } if `nmp' == 4 { eret scalar pr1 = 1 - `e(pr2)' - `e(pr3)' - `e(pr4)' } if `nmp' == 5 { eret scalar pr1 = 1 - `e(pr2)' - `e(pr3)' - `e(pr4)' - `e(pr5)' } } Display, `level' `eform' end program define Display , eclass syntax [,Level(int $S_level) EForm *] local diopts "`options'" ml display, level(`level') `eform' plus if `level' < 10 | `level' > 99 { local level = 95 } if `e(nmp)' == 2 { _diparm logitp2, level(`level') func(1-invlogit(@)) /// deriv( -invlogit(@)*(1-invlogit(@)) ) label("Prob. Type 1") prob ci(logit) eret scalar se_pr1 = r(se) _diparm logitp2, invlogit level(`level') label("Prob. Type 2") prob ci(logit) } if `e(nmp)' == 3 { _diparm logitp2 logitp3, level(`level') func(1-invlogit(@1)-invlogit(@2)) /// deriv( -invlogit(@1)*(1-invlogit(@1)) /// -invlogit(@2)*(1-invlogit(@2)) ) /// label("Prob. Type 1") prob ci(logit) eret scalar se_pr1 = r(se) _diparm logitp2, invlogit level(`level') label("Prob. Type 2") prob ci(logit) _diparm logitp3, invlogit level(`level') label("Prob. Type 3") prob ci(logit) } if `e(nmp)' == 4 { _diparm logitp2 logitp3 logitp4, level(`level') /// func(1-invlogit(@1)-invlogit(@2)-invlogit(@3)) /// deriv( -invlogit(@1)*(1-invlogit(@1)) /// -invlogit(@2)*(1-invlogit(@2)) /// -invlogit(@3)*(1-invlogit(@3)) ) /// label("Prob. Type 1") prob ci(logit) eret scalar se_pr1 = r(se) _diparm logitp2, invlogit level(`level') label("Prob. Type 2") prob ci(logit) _diparm logitp3, invlogit level(`level') label("Prob. Type 3") prob ci(logit) _diparm logitp4, invlogit level(`level') label("Prob. Type 4") prob ci(logit) } if `e(nmp)' == 5 { _diparm logitp2 logitp3 logitp4 logitp5, level(`level') /// func(1-invlogit(@1)-invlogit(@2)-invlogit(@3)-invlogit(@4)) /// deriv( -invlogit(@1)*(1-invlogit(@1)) /// -invlogit(@2)*(1-invlogit(@2)) /// -invlogit(@3)*(1-invlogit(@3)) /// -invlogit(@4)*(1-invlogit(@4)) ) /// label("Prob. Type 1") prob ci(logit) eret scalar se_pr1 = r(se) _diparm logitp2, invlogit level(`level') label("Prob. Type 2") prob ci(logit) _diparm logitp3, invlogit level(`level') label("Prob. Type 3") prob ci(logit) _diparm logitp4, invlogit level(`level') label("Prob. Type 4") prob ci(logit) _diparm logitp5, invlogit level(`level') label("Prob. Type 5") prob ci(logit) } di in smcl in gr "{hline 13}{c BT}{hline 64}" di in gr "Note: m1 = " as res "0" end