*! version 1.1.2 11may2011 Ben Jann * based on cloglog.ado, version 1.8.6 04aug2005 program define rrlogit, byable(onecall) prop(/*ml_score*/ swml or) version 9.1 if _by() { local BY `"by `_byvars'`_byrc0':"' } `BY' _vce_parserun rrlogit : `0' if "`s(exit)'" != "" { exit } if replay() { if "`e(cmd)'" != "rrlogit" error 301 if _by() error 190 Display `0' error `e(rc)' exit } `BY' Estimate `0' end program Estimate, eclass byable(recall) /* Parse. */ syntax varlist(numeric ts) [if] [in] [fw iw pw] [, /* */ PWarner(str) PYes(str) p1(str) PNo(str) p2(str) /* */ ASIS FROM(string) /* */ Level(cilevel) or noCONstant Robust CLuster(varname) /* */ OFFset(varname numeric) SCore(string) noLOg noDISPLAY * ] if `"`pwarner'"'=="" local pwarner 1 if `"`p1'"'=="" local p1 `"`pyes'"' if `"`p2'"'=="" local p2 `"`pno'"' if `"`p1'"'=="" local p1 0 if `"`p2'"'=="" local p2 0 if `"`robust'"' != "" | `"`cluster'"' != "" | /* */ `"`weight'"' == "pweight" { local crtype crittype("log pseudolikelihood") } if _by() { _byoptnotallowed score() `"`score'"' } mlopts mlopts, `options' /* Check syntax. */ if `"`score'"'!="" { confirm new variable `score' local nword : word count `score' if `nword' > 1 { di as err "score() must contain the name of only" /* */ " one new variable" exit 198 } tempvar scvar local scopt "score(`scvar')" } if "`constant'"!="" { local nvar : word count `varlist' if `nvar' == 1 { di as err "independent variables required with " /* */ "noconstant option" exit 102 } } /* Mark sample. */ marksample touse, zeroweight if `"`cluster'"'!="" { markout `touse' `cluster', strok local clopt cluster(`cluster') } if "`offset'"!="" { markout `touse' `offset' local offopt "offset(`offset')" } /* Check p's. */ capt assert (`pwarner'>=0) & (`pwarner'<=1) & (`pwarner'!=0.5) if `touse' if _rc { di as err "pwarner() must be in [0,1]; pwarner() must be unequal 0.5" exit 198 } capt assert (`p1'>=0) & (`p2'>=0) & ((`p1'+`p2')<=1) if `touse' if _rc { di as err "pyes() and pno() must be in [0,1]; pyes()+pno() must be < 1" exit 198 } /* Count obs and check values of `y'. */ gettoken y xvars : varlist tsunab yname : `y' loc yname : subinstr local yname "." "_" qui count if `touse' local n `r(N)' if `n' == 0 error 2000 if `n' == 1 error 2001 qui count if `y'==0 & `touse' local n0 `r(N)' if `n0'==0 | `n0'==`n' { di as err "outcome does not vary; remember:" di as err _col(35) "0 = negative outcome," di as err _col(9) /* */ "all other nonmissing values = positive outcome" exit 2000 } if "`log'"!="" | "`display'"!="" { local qui "quietly" } /* If there are negative iweights, -logit- cannot be used. */ local nonegwt 1 if "`weight'"!="" { if "`weight'"=="pweight" { /* pweights create unneeded extra work when getting initial values, etc. */ local wtype "iweight" } else local wtype "`weight'" if "`weight'"=="iweight" { tempname sumw tempvar w qui gen double `w' `exp' if `touse' summarize `w' if `touse', meanonly scalar `sumw' = r(sum) if `sumw' <= 0 { di as err "sum of weights less than " /* */ "or equal to zero" exit 402 } if r(min) < 0 local nonegwt 0 } } /* Remove collinearity. */ _rmcoll `xvars' [`weight'`exp'] if `touse', `constant' local xvars `r(varlist)' /* Run logit to drop any variables and observations. */ if "`asis'"=="" & `nonegwt' { logit `y' `xvars' [`wtype'`exp'] if `touse', /* */ iter(0) `offopt' nocoef nolog `constant' qui replace `touse' = e(sample) _evlist local xvars `s(varlist)' } /* Set rrt probabilities. */ global rrlogit_pw `pwarner' global rrlogit_p1 `p1' global rrlogit_p2 `p2' /* Compute constant-only model. */ if "`constant'"=="" & "`xvars'"!="" { if `"`from'"'=="" { `qui' di as txt _n "Fitting constant-only model:" ml model d2 rrlogit_lf /* */ (`yname': `y'=, `constant' `offopt') /* */ if `touse' [`wtype'`exp'], collinear missing max /* */ nooutput nopreserve wald(0) search(off) `mlopts' /* */ `log' `crtype' nocnsnotes local continu "continue search(off)" `qui' di as txt _n "Fitting full model:" } } /* Fit full model. */ if `"`from'"'!="" local initopt `"init(`from')"' ml model d2 rrlogit_lf (`yname': `y'=`xvars', `constant' `offopt') /* */ if `touse' [`weight'`exp'], collinear missing max nooutput /* */ nopreserve `initopt' `lf0' `mlopts' `log' /* */ `scopt' `robust' `clopt' `continu' /* */ title(Randomized response logistic regression) /* Clear rrt probabilities. */ global rrlogit_pw global rrlogit_p1 global rrlogit_p2 /* Returns. */ eret local cmd if "`score'" != "" { label var `scvar' "Score index for x*b from rrlogit" rename `scvar' `score' eret local scorevars `score' } if "`weight'" == "fweight" { tempvar tmpsum qui gen double `tmpsum' `exp' if `touse' & `y' == 0 summarize `tmpsum' if `touse' & `y' == 0 , meanonly eret scalar N_f = r(sum) } else { qui count if `touse' & `y' == 0 eret scalar N_f = r(N) } eret scalar N_s = e(N) - e(N_f) eret scalar r2_p = 1 - e(ll)/e(ll_0) eret local offset1 eret local offset "`offset'" eret local predict "rrlogit_p" eret local cmd "rrlogit" eret local cmd2 "logit" eret local pwarner `pwarner' eret local pyes `p1' eret local pno `p2' if "`display'"=="" { Display, level(`level') `or' } error `e(rc)' end program define Display syntax [, Level(cilevel) or ] if e(chi2) < 1e5 { local fmt "%9.2f" } else local fmt "%9.2e" local crtype = upper(substr(`"`e(crittype)'"',1,1)) + /* */ substr(`"`e(crittype)'"',2,.) di as txt _n "`e(title)'" _col(49) "Number of obs =" /* */ as res _col(70) %9.0g e(N) _n /* */ as txt _col(49) "Nonzero outcomes =" /* */ as res _col(70) %9.0g e(N_s) _n /* */ as txt `"P(non-negated question) = "' as res e(pwarner) /* */ as txt _col(49) "Zero outcomes =" /* */ as res _col(70) %9.0g e(N_f) _n /* */ as txt `"P(surrogate "yes") = "' as res e(pyes) /* */ as txt _col(49) "`e(chi2type)' chi2(" as res `e(df_m)' /* */ as txt ")" _col(67) "=" as res _col(70) `fmt' e(chi2) _n /* */ as txt `"P(surrogate "no") = "' as res e(pno) /* */ as txt _col(49) "Prob > chi2" _col(67) "=" as res _col(70) /* */ as res %9.4f chiprob(e(df_m),e(chi2)) _n /* */ as txt "`crtype' = " as res %10.0g e(ll) /* */ as txt _col(49) "Pseudo R2" _col(67) "=" as res _col(70) /* */ as res %9.4f e(r2_p) _n ml di, noheader first level(`level') nofootnote `or' _prefix_footnote end