*! version 1.0.4 28jan2009 Ben Jann prog logitcprplot, sortpreserve version 9.2 syntax varlist [if] [in] [ , SAMPle(numlist max==1 <=100 >=1) /// Generate(name) /// Replace /// NORLine RLine(str asis) /// Lowess Lowess2(str asis) /// LPoly LPoly2(str asis) /// Stata 10 required FPfit FPfit2(str asis) /// RCspline RCspline2(str asis) /// PSpline PSpline2(str asis) /// -pspline- required RSpline RSpline2(str asis) /// -uvrs- required RUNning RUNning2(str asis) /// -running- required AUTOsmoo AUTOsmoo2(str asis) /// -autosmoo- required noGRaph /// addplot(str asis) * ] if `"`generate'"'!="" & "`replace'"=="" { confirm new var `generate' } if `"`lpoly2'"'!="" local lpoly lpoly _parse_lpoly, `lpoly2' if `"`fpfit2'"'!="" local fpfit fpfit _parse_fpfit, `fpfit2' if `"`lowess2'"'!="" local lowess lowess if !inlist(`"`e(cmd)'"',"logit","logistic") { di as err "last estimates not found" exit 301 } if "`sample'"!="" { capt which gsample local rc = _rc if _rc { di as error "-gsample- user command required; see {stata ssc describe gsample}" } capt findfile lmoremata.mlib if _rc { di as error "-moremata- user package required; see {stata ssc describe moremata}" } if `rc' | _rc exit 499 } // compute predictions and partial residuals marksample touse, novar qui replace `touse' = 0 if e(sample)!=1 local wtype `"`e(wtype)'"' local wexp `"`e(wexp)'"' tempvar pr r xb rplus local depvar = e(depvar) qui predict `pr' if `touse', pr qui gen `r' = (`depvar' - `pr') / (`pr'*(1-`pr')) local terms foreach var of local varlist { local terms `"`terms' + _b[`var']*`var'"' } qui gen `xb' = 0 `terms' if `touse' qui gen `rplus' = `r' + `xb' gettoken firstvar : varlist // compute rcspline smooth if `"`rcspline'`rcspline2'"'!="" { local rcspline rcspline tempvar rcvar rcci_l rcci_u _rcspline `rplus' `firstvar' if `touse' [`wtype'`wexp'], /// gen(`rcvar' `rcci_l' `rcci_u') `rcspline2' } // compute pspline smooth if `"`pspline'`pspline2'"'!="" { capt which pspline if _rc { di as error "-pspline- user command required; see {stata ssc describe pspline}" error 499 } local pspline pspline tempvar psvar _pspline `rplus' `firstvar' if `touse' [`wtype'`wexp'], /// gen(`psvar') `pspline2' } // compute rspline smooth if `"`rspline'`rspline2'"'!="" { capt which uvrs if _rc { di as error "-uvrs- user command required; see {stata net sj 7-1 st0120}" error 499 } local rspline rspline tempvar rsvar rsci_l rsci_u _rspline `rplus' `firstvar' if `touse' [`wtype'`wexp'], /// gen(`rsvar' `rsci_l' `rsci_u') `rspline2' } // compute running smooth if `"`running'`running2'"'!="" { capt which running if _rc { di as error "-running- user command required; see {stata net sj 5-2 sed9_2}" error 499 } local running running tempvar runvar runci_l runci_u _running `rplus' `firstvar' if `touse' [`wtype'`wexp'], /// gen(`runvar' `runci_l' `runci_u') `running2' } // compute autosmoo smooth if `"`autosmoo'`autosmoo2'"'!="" { capt which autosmoo if _rc { di as error "-autosmoo- user command required; see {stata net stb 41 gr27}" error 499 } local autosmoo autosmoo tempvar autovar _autosmoo `rplus' `firstvar' if `touse' [`wtype'`wexp'], /// gen(`autovar') `autosmoo2' } // select obs to be plotted; also affects lowess and lpoly (rcspline will always use all obs) if "`sample'"!="" { tempname touse2 gsample `sample' if `touse', wor generate(`touse2') percent } else { local touse2 `touse' } tempname touse3 // unique X values qui bys `touse' `firstvar': gen byte `touse3' = _n==1 & `touse' // compile graph local i = 0 local p = 1 + ("`norline'"=="") local lord if `"`lpci'"'!="" { local lpcigraph (lpolyci `rplus' `firstvar' if `touse2' [`wtype'`wexp'], /// clstyle(p`++p') yvarlabel("lpoly smooth") `lpoly2') local ++i local lord "`lord' `++i'" } if `"`fpci'"'!="" { local fpcigraph (fpfitci `rplus' `firstvar' if `touse' [`wtype'`wexp'], /// clstyle(p`++p') yvarlabel("fracpoly fit") `fpfit2') local ++i local lord "`lord' `++i'" } if `"`rcspline'"'!="" { if `"`rcci'"'!="" { local rccigraph (rarea `rcci_l' `rcci_u' `firstvar' if `touse3', /// sort yvarlabel("rcspline CI") `rcci2') local ++i } } if `"`rspline'"'!="" { if `"`rsci'"'!="" { local rscigraph (rarea `rsci_l' `rsci_u' `firstvar' if `touse3', /// sort yvarlabel("rspline CI") `rsci2') local ++i } } if `"`running'"'!="" { if `"`runningci'"'!="" { local runciraph (rarea `runci_l' `runci_u' `firstvar' if `touse3', /// sort yvarlabel("running CI") `runningci2') local ++i } } local ++i if "`norline'"=="" { local rlgraph (line `xb' `firstvar' if `touse3', /// sort pstyle(p2) yvarlabel("model fit") `rline') local lord "`++i' `lord'" } if `"`rcspline'"'!="" { local rcgraph (line `rcvar' `firstvar' if `touse3', /// sort pstyle(p`++p') yvarlabel("rcspline smooth") `rcspline2') local lord "`lord' `++i'" } if `"`pspline'"'!="" { local psgraph (line `psvar' `firstvar' if `touse3', /// sort pstyle(p`++p') yvarlabel("pspline smooth") `pspline2') local lord "`lord' `++i'" } if `"`rspline'"'!="" { local rsgraph (line `rsvar' `firstvar' if `touse3', /// sort pstyle(p`++p') yvarlabel("rspline smooth") `rspline2') local lord "`lord' `++i'" } if `"`running'"'!="" { local rungraph (line `runvar' `firstvar' if `touse3', /// sort pstyle(p`++p') yvarlabel("running smooth") `running2') local lord "`lord' `++i'" } if `"`autosmoo'"'!="" { local autograph (line `autovar' `firstvar' if `touse3', /// sort pstyle(p`++p') yvarlabel("auto smooth") `autosmoo2') local lord "`lord' `++i'" } if `"`lpoly'"'!="" & `"`lpci'"'=="" { local lpgraph (lpoly `rplus' `firstvar' if `touse2' [`wtype'`wexp'], /// pstyle(p`++p') yvarlabel("lpoly smooth") `lpoly2') local lord "`lord' `++i'" } if `"`fpfit'"'!="" & `"`fpci'"'=="" { local fpgraph (fpfit `rplus' `firstvar' if `touse' [`wtype'`wexp'], /// pstyle(p`++p') pstyle(p`++p') yvarlabel("fracpoly fit") `fpfit2') local lord "`lord' `++i'" } if `"`lowess'"'!="" { local lsgraph (lowess `rplus' `firstvar' if `touse2' [`wtype'`wexp'], /// pstyle(p`++p') yvarlabel("lowess smooth") `lowess2') local lord "`lord' `++i'" } if `"`graph'"'=="" { twoway `lpcigraph' `fpcigraph' `rccigraph' `rscigraph' `runciraph' /// (scatter `rplus' `firstvar' if `touse2', /// pstyle(p1) yti("partial logit residual") xtit(`firstvar') /// legend(order(`lord') cols(3)) `options') /// `rlgraph' `rcgraph' `psgraph' `rsgraph' `rungraph' `autograph' /// `lpgraph' `fpgraph' `lsgraph' /// || `addplot' } // return results if `"`generate'"'!="" { if "`replace'"!="" { capt confirm var `generate', exact if _rc==0 drop `generate' } rename `rplus' `generate' } end program _parse_lpoly syntax [, ci CI2(str asis) * ] if `"`ci2'"'!="" { local ci ci } if "`ci'"!="" { local ci2 pstyle(ci) `ci2' } c_local lpci `ci' c_local lpoly2 `ci2' `options' end program _parse_fpfit syntax [, ci CI2(str asis) * ] if `"`ci2'"'!="" { local ci ci } if "`ci'"!="" { local ci2 pstyle(ci) `ci2' } c_local fpci `ci' c_local fpfit2 `ci2' `options' end program _rcspline // loosely following -rcspline- by Nick Cox syntax varlist(numeric min=2 max=2) [if] [in] [fw aw iw pw], gen(str) [ /// NKnots(passthru) Knots(passthru) DIsplayknots regressopts(str) /// CI CI2(str asis) Level(cilevel) * ] if `"`ci2'"'!="" { local ci ci } if "`ci'"!="" { local ci2 pstyle(ci) `ci2' } c_local rcci `ci' c_local rcci2 `ci2' c_local rcspline2 `options' gettoken y x : varlist gettoken pred gen : gen gettoken ci_l gen : gen gettoken ci_u gen : gen tempname hcurrent _est hold `hcurrent', restore marksample touse tempname stub * mkspline `stub' = `x' if `touse' [`weight'`exp'], cubic `nknots' `knots' `displayknots' cubicspline `stub' = `x' if `touse' [`weight'`exp'], `nknots' `knots' `displayknots' qui regress `y' `stub'* if `touse' [`weight' `exp'] , `regressopts' qui predict `pred' if e(sample) if "`ci'"!="" { tempvar se qui predict `se' if e(sample), stdp local z = invttail(e(df_r), (100 - `level') / 200) qui gen `ci_l' = `pred' - `z' * `se' qui gen `ci_u' = `pred' + `z' * `se' } end program _pspline syntax varlist(numeric min=2 max=2) [if] [in] [fw aw iw pw], gen(str) [ /// Degree(passthru) NKnots(passthru) Knots(passthru) ESTOPts(passthru) noPENalty /// * ] c_local pspline2 `options' tempname hcurrent _est hold `hcurrent', restore qui pspline `varlist' `if' `in' [`weight' `exp'], /// nograph gen(`gen') /// `degree' `nknots' `knots' `estopts' `penalty' end program _rspline syntax varlist(numeric min=2 max=2) [if] [in] [fw aw iw pw], gen(str) [ /// ALpha(passthru) DEGree(passthru) df(passthru) KNots(passthru) regressopts(str) /// CI CI2(str asis) Level(cilevel) * ] if `"`ci2'"'!="" { local ci ci } if "`ci'"!="" { local ci2 pstyle(ci) `ci2' } c_local rsci `ci' c_local rsci2 `ci2' c_local rspline2 `options' gettoken y x : varlist gettoken pred gen : gen gettoken ci_l gen : gen gettoken ci_u gen : gen tempname hcurrent _est hold `hcurrent', restore tempname X qui gen `X' = `x' // so that variables added by -uvrs- are dropped qui uvrs regress `y' `X' `if' `in' [`weight' `exp'], /// `alpha' `degree' `df' `knots' `regressopts' qui fracpred `pred' if "`ci'"!="" { tempvar se qui fracpred `se', stdp local z = invttail(e(df_r), (100 - `level') / 200) qui gen `ci_l' = `pred' - `z' * `se' qui gen `ci_u' = `pred' + `z' * `se' } end program _running syntax varlist(numeric min=2 max=2) [if] [in] [fw aw iw pw], gen(str) [ /// Double Knn(passthru) Mean Repeat(passthru) SPan(passthru) TWice /// CI CI2(str asis) Level(cilevel) * ] if `"`ci2'"'!="" { local ci ci } if "`ci'"!="" { local ci2 pstyle(ci) `ci2' } c_local runningci `ci' c_local runningci2 `ci2' c_local running2 `options' gettoken pred gen : gen gettoken ci_l gen : gen gettoken ci_u gen : gen if "`ci'"!="" { tempvar se local gense gense(`se') } qui running `varlist' `if' `in' [`weight' `exp'], /// nograph generate(`pred') `gense' /// `double' `knn' `mean' `repeat' `span' `twice' if "`ci'"!="" { local z = invnormal(1 - (100 - `level') / 200) qui gen `ci_l' = `pred' - `z' * `se' qui gen `ci_u' = `pred' + `z' * `se' } end program _autosmoo syntax varlist(numeric min=2 max=2) [if] [in] [fw aw iw pw], gen(str) [ /// kmin(passthru) kmax(passthru) Repeat(passthru) /// * ] c_local autosmoo2 `options' qui autosmoo `varlist' `if' `in' [`weight' `exp'], /// nograph gen(`gen') /// `kmin' `kmax' `repeat' end * based on code from mkspline.ado by StataCorp, version 1.2.3, 14may2007 * adapted to allow pweights * - without weights: -altdef- formula => same as mkspline * - with fweight: standard formula => possibly different than mkspline (?) * - with pweight: standard formula => not supported in mkspline program define cubicspline, rclass sortpreserve version 6, missing gettoken name 0 : 0, parse(" =") gettoken eqsign : 0, parse(" =") gettoken eqsign 0 : 0, parse(" =") syntax varname [if] [in] [fweight pweight] [, /// NKnots(numlist max=1) Knots(numlist) DIsplayknots ] marksample touse if "`nknots'"!="" { local nk `nknots' } else { local nk=5 } if "`knots'"!="" { local nc 0 tokenize "`knots'" while "`1'" != "" { local nc = `nc' + 1 local t`nc' "`1'" mac shift } } if "`nknots'"!="" & "`knots'"=="" { local nc `nk' } if "`nknots'"=="" & "`knots'"!="" { local nk `nc' } if "`nknots'"!="" & "`knots'"!="" { if `nc' != `nk' { display as error /// "Count in nknots must be the same as the number of knots specified." error 498 } } sort `varlist' if "`knots'"!="" { if `nc' < 2 { display as error /// "Restricted cubic splines must have at least 2 knots." error 498 } local prevt=`t1' local j=2 while `j'<=`nk' { if `t`j'' <= `prevt' { disp as error /// "Knots must be in increasing order." error 498 } local prevt=`t`j'' local j=`j'+1 } } else { local altdef altdef if "`weight'"!="" { local altdef if "`weight'"=="pweight" { local suwgt "aweight" } else { local suwgt "`weight'" } } if `nk' == 3 { _pctile `varlist' if `touse' [`weight' `exp'], /// percentiles(10 50 90) `altdef' } else if `nk'== 4 { _pctile `varlist' if `touse' [`weight' `exp'], /// percentiles(5 35 65 95) `altdef' } else if `nk'== 5 { _pctile `varlist' if `touse' [`weight' `exp'], /// percentiles(5 27.5 50 72.5 95) `altdef' } else if `nk'== 6 { _pctile `varlist' if `touse' [`weight' `exp'], /// percentiles(5 23 41 59 77 95) `altdef' } else if `nk'== 7 { _pctile `varlist' if `touse' [`weight' `exp'], /// percentiles(2.5 18.33 34.17 50 65.83 81.67 97.5) `altdef' } else { display as error /// "Restricted cubic splines with `nk' knots at default values not implemented." display as error /// "Number of knots specified in nknots() must be between 3 and 7." error 498 } forvalues i=1 / `nk' { local t`i' = r(r`i') } qui sum `varlist' [`suwgt' `exp'], meanonly local cenmin = 5/(r(N)+1)*100 local cenmax = (r(N)-4)/(r(N)+1)*100 _pctile `varlist' if `touse' [`weight' `exp'], /// percentiles(`cenmin' `cenmax') `altdef' local min = r(r1) local max = r(r2) if `t1' < `min' { local t1 = `min' } if `t`nk'' > `max' { local t`nk' = `max' } } local j = 1 matrix loc = J(1,`nk',0) local col "" while `j' <= `nk' { mat loc[1,`j'] = `t`j'' local col "`col' knot`j'" local j = `j'+1 } mat colnames loc = `col' matrix rownames loc = "`varlist'" local km1 = `nk' - 1 if `t1' >= `t2' | `t`nk'' <= `t`km1'' { display as error /// "Sample size too small for this many knots." error 198 } qui gen `name'1=`varlist' if `touse' local j = 1 while `j' <= `nk' { qui gen _Xm`j' = `varlist' - `t`j'' uplus _Xm`j' _Xm`j'p local j = `j'+1 } local j = 1 while `j' <= `nk' -2 { local jp1 = `j' + 1 qui gen `name'`jp1' = (_Xm`j'p^3 - (_Xm`km1'p^3)* /// (`t`nk'' - `t`j'')/(`t`nk'' - `t`km1'') /// + (_Xm`nk'p^3 )*(`t`km1'' - `t`j'')/ /// (`t`nk'' - `t`km1'')) / (`t`nk'' - `t1')^2 /// if `touse' local j = `j' + 1 } drop _Xm* _Xm*p if "`displayknots'"!="" { matlist loc } return scalar N_knots = `nk' return matrix knots loc end program define uplus version 7 args u uplus qui gen `uplus' = `u' quietly replace `uplus' = 0 if `u' <= 0 end