*! 1.0.2 MLB 01 Aug 2012 *! 1.0.1 MLB 25 May 2010 *! 1.0.0 MLB 26 Apr 2008 * Fit a zero one inflated beta /*------------------------------------------------ playback request */ program zoib, eclass byable(onecall) if c(stata_version) >= 11 { version 11 global fv "fv" } else if c(stata_version) >= 9 { version 9 } else { version 8.2 } if replay() { if "`e(cmd)'" != "zoib" { di as err "results for zoib 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 Estimate, eclass byable(recall) syntax varlist(numeric $fv) [if] [in] [fw pw aw] [, /// noONE noZERO ZEROInflate(varlist $fv) ONEInflate(varlist $fv) /// PHIvar(varlist $fv) Cluster(varname) Level(integer $S_level) noLOG Robust * ] if "`one'" != "" & "`oneinflate'"!= "" { di as err "cannot specify both noone and oneinflate" exit 198 } if "`zero'" != "" & "`zeroinflate'"!= "" { di as err "cannot specify both nozero and zeroinflate" exit 198 } gettoken y x: varlist if c(stata_version) >= 11 { _fv_check_depvar `y' } marksample touse markout `touse' `varlist' `zeroinflate' `oneinflate' `cluster' qui count if (`y' < 0 | `y' > 1) & `touse' if r(N) { noi di " " noi di as txt "warning: {res:`y'} has `r(N)' values < 0 or > 1;" _c noi di as txt " not used in calculations" } qui replace `touse' = 0 if `y' < 0 | `y' > 1 // No inflate equation when there are no zeros or ones in the dependent variable qui count if `y' == 0 & `touse' if r(N) == 0 { if "`zeroinflate'" != "" { di as err "the zeroinflate() option can only be specified if the dependent variable contains zeros" exit 198 } local zero "zero" } qui count if `y' == 1 & `touse' if r(N) == 0 { if "`oneinflate'" != "" { di as err "the oneinflate() option can only be specified if the dependent variable contains ones" exit 198 } local one "one" } local param = cond("`one'" != "", /* */ cond("`zero'" != "", "beta", "zib"), /* */ cond("`zero'" != "","oib","zoib")) if "`one'" != "" { qui count if `y' == 1 if r(N) > 0 { noi di as txt "warning: {res:`y'} has `r(N)' values == 1;" _c noi di as txt " not used in calculations" } qui replace `touse' = 0 if `y' == 1 } if "`zero'" != "" { qui count if `y' == 0 if r(N) > 0 { noi di as txt "warning: {res:`y'} has `r(N)' values == 0;" _c noi di as txt " not used in calculations" } qui replace `touse' = 0 if `y' == 0 } qui count if `touse' if r(N) == 0 error 2000 local title "ML fit of `param'" local wtype `weight' local wtexp `"`exp'"' if "`weight'" != "" local wgt `"[`weight'`exp']"' if "`cluster'" != "" { local clopt "cluster(`cluster')" } if "`level'" != "" local level "level(`level')" local log = cond("`log'" == "", "noisily", "quietly") mlopts mlopts, `options' Init `y' if `touse' `wgt', param("`param'") tempname init matrix `init' = r(init) if "`param'" == "zib" { `log' ml model lf zib_lf (proportion: `y' = `x') /// (zeroinflate: `zeroinflate') /// (ln_phi: `phivar') /// `wgt' if `touse' , maximize /// collinear title(`title') /// search(off) `clopt' `level' /// `mlopts' `stdopts' `modopts' /// init(`init') `robust' } if "`param'" == "oib" { `log' ml model lf oib_lf (proportion: `y' = `x') /// (oneinflate: `oneinflate') /// (ln_phi: `phivar') /// `wgt' if `touse' , maximize /// collinear title(`title') /// search(off) `clopt' `level' /// `mlopts' `stdopts' `modopts' /// init(`init') `robust' } if "`param'" == "zoib" { `log' ml model lf zoib_lf (proportion: `y' = `x') /// (oneinflate: `oneinflate') /// (zeroinflate: `zeroinflate') /// (ln_phi: `phivar') /// `wgt' if `touse' , maximize /// collinear title(`title') /// search(off) `clopt' `level' /// `mlopts' `stdopts' `modopts' /// init(`init') `robust' } if "`param'" == "beta" { `log' ml model lf zoib_beta_lf (proportion: `y' = `x') /// (ln_phi: `phivar') /// `wgt' if `touse' , maximize /// collinear title(`title') /// search(off) `clopt' `level' /// `mlopts' `stdopts' `modopts' /// init(`init') `robust' } eret local cmd "zoib" eret local depvar "`y'" ereturn local predict "zoib_p" Display, `level' `diopts' end program Init, rclass syntax varname [if] [aw pw fw], param(string) tempname init mu phi m var marksample touse if "`weight'" == "pweight" { local wgt "[aweight`exp']" } else { local wgt "[`weight'`exp']" } qui sum `varlist' `wgt' if `touse' & `varlist' > 0 & `varlist' < 1 scalar `m' = r(mean) scalar `var' = r(Var) if "`param'" == "beta" { scalar `mu' = ln(`m'/(1-`m')) scalar `phi' = ln((`m'*(1-`m'))/`var' - 1) matrix `init' = `mu', `phi' matrix colnames `init' = proportion:_cons ln_phi:_cons } else if "`param'" == "zib" { scalar `mu' = ln(`m'/(1-`m')) scalar `phi' = ln((`m'*(1-`m'))/`var' - 1) tempvar zero tempname zi qui gen byte `zero' = `varlist' == 0 if `touse' sum `zero', meanonly scalar `zi' = ln(`r(mean)'/(1-`r(mean)')) matrix `init' = `mu', `zi', `phi' matrix colnames `init' = proportion:_cons zeroinflate:_cons ln_phi:_cons } else if "`param'" == "oib" { scalar `mu' = ln(`m'/(1-`m')) scalar `phi' = ln((`m'*(1-`m'))/`var' - 1) tempvar one tempname oi qui gen byte `one' = `varlist' == 1 if `touse' sum `one', meanonly scalar `oi' = ln(`r(mean)'/(1-`r(mean)')) matrix `init' = `mu', `oi', `phi' matrix colnames `init' = proportion:_cons oneinflate:_cons ln_phi:_cons } else if "`param'" == "zoib" { scalar `mu' = ln(`m'/(1-`m')) scalar `phi' = ln((`m'*(1-`m'))/`var' - 1) tempvar zero tempname zi qui gen byte `zero' = `varlist' == 0 if `touse' sum `zero', meanonly scalar `zi' = ln(`r(mean)'/(1-`r(mean)')) tempvar one tempname oi qui gen byte `one' = `varlist' == 1 if `touse' sum `one', meanonly scalar `oi' = ln(`r(mean)'/(1-`r(mean)')) matrix `init' = `mu', `zi', `oi', `phi' matrix colnames `init' = proportion:_cons zeroinflate:_cons oneinflate:_cons ln_phi:_cons } return matrix init = `init' end program Display syntax [, Level(int $S_level) *] local diopts "`options'" if `level' < 10 | `level' > 99 local level = 95 ml display, level(`level') `or' `diopts' end