You cannot select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
836 lines
30 KiB
Plaintext
836 lines
30 KiB
Plaintext
10 months ago
|
*! version 1.9.0 2012-07-29 fixes for o. in stata 11+
|
||
|
|
capture program drop fitstat
|
||
|
|
program define fitstat, rclass
|
||
|
|
|
||
|
|
local ecmdversion = c(version) // 190
|
||
|
|
|
||
|
|
version 9.0
|
||
|
|
local cmd `e(cmd)'
|
||
|
|
if ( "`cmd'"!="regress" & "`cmd'"!="logit" & "`cmd'"!="probit" ///
|
||
|
|
& "`cmd'"!="logistic" & "`cmd'"!="cloglog" & "`cmd'"!="ologit" ///
|
||
|
|
& "`cmd'"!="oprobit" & "`cmd'"!="gologit" & "`cmd'"!="mlogit" ///
|
||
|
|
& "`cmd'"!="clogit" & "`cmd'"!="tobit" & "`cmd'"!="cnreg" ///
|
||
|
|
& "`cmd'"!="intreg" & "`cmd'"!="poisson" & "`cmd'"!="nbreg" ///
|
||
|
|
& "`cmd'"!="zip" & "`cmd'"!="zinb" ///
|
||
|
|
& "`cmd'"!="ztp" & "`cmd'"!="ztnb" & "`cmd'"!="slogit" ///
|
||
|
|
& "`cmd'"!="mprobit" & "`cmd'"!="rologit" ///
|
||
|
|
) {
|
||
|
|
di in r "fitstat does not work with the last model estimated."
|
||
|
|
exit
|
||
|
|
}
|
||
|
|
|
||
|
|
tempname n_obs n_obs2 n_parm llF llR df x2 r2 b v_ystar v v_y v_x v_e
|
||
|
|
tempname factor r2mcf r2mcfa r2mz r2ml r2cu r2 r2ef r2adj aic aicn bicx
|
||
|
|
tempname bicp isbin isord isvare isystar islgt ispbt ssr dev lrx2 lrx2p
|
||
|
|
tempname depv2 maxtmp uniqgrp n_obs2 depvgrp ifin v_error bicdif fllR
|
||
|
|
tempname bnocon phat yphat ydev predpr predcat wrong hiprob fbicp fn_obs
|
||
|
|
tempname totcor totone maxrow cntr2 cntr2a chkobs d d2 prevval prevdf
|
||
|
|
tempname results erank Vb noomit
|
||
|
|
|
||
|
|
syntax [if] [in] ///
|
||
|
|
[, Saving(string) Using(string) SAVE Diff FORCE Bic NOSTATA ///
|
||
|
|
Debug(integer 0) ///
|
||
|
|
]
|
||
|
|
|
||
|
|
* 181 mat `b' = get(_b)
|
||
|
|
* 190: remove o. from beta matrix
|
||
|
|
mat `b' = get(_b)
|
||
|
|
mat `Vb' = e(V)
|
||
|
|
local colfn : colfullnames `b'
|
||
|
|
local coln : colnames `b'
|
||
|
|
local coleq : coleq `b'
|
||
|
|
_ms_omit_info `b' // remove columns of zeros from b
|
||
|
|
local cols = colsof(`b')
|
||
|
|
mat `noomit' = J(1,`cols',1) - r(omit)
|
||
|
|
forvalues i = 1/`cols' { // double clear for 0 std error
|
||
|
|
if `Vb'[`i',`i'] == 0 & `noomit'[1,`i']==1 {
|
||
|
|
mat `noomit'[1,`i']==0
|
||
|
|
}
|
||
|
|
}
|
||
|
|
mata: newbeta = ///
|
||
|
|
select(st_matrix(st_local("b")), (st_matrix(st_local("noomit"))))
|
||
|
|
mata: st_matrix(st_local("b"),newbeta)
|
||
|
|
// reassign column names to b
|
||
|
|
local i = 0
|
||
|
|
foreach var of local coln {
|
||
|
|
local ++i
|
||
|
|
local nmvar : word `i' of `coln'
|
||
|
|
local nmeq : word `i' of `coleq'
|
||
|
|
local nmfn : word `i' of `colfn'
|
||
|
|
_ms_parse_parts `var'
|
||
|
|
if "`nmvar'"=="`nmfn'" { // no eq nm
|
||
|
|
local nmeq ""
|
||
|
|
}
|
||
|
|
else { // eq nm
|
||
|
|
local nmeq "[`nmeq']"
|
||
|
|
}
|
||
|
|
if !`r(omit)' & `nmeq'_se[`nmvar']!=0 {
|
||
|
|
local coln2 `coln2' `nmfn'
|
||
|
|
}
|
||
|
|
}
|
||
|
|
version 11: matrix colnames `b' = `coln2'
|
||
|
|
|
||
|
|
local depv `e(depvar)'
|
||
|
|
|
||
|
|
if "`nostata'"!="nostata" { // stata bic and aic
|
||
|
|
tempname stataic statabic stataaic
|
||
|
|
qui estat ic
|
||
|
|
mat `stataic' = r(S)
|
||
|
|
scalar `statabic' = `stataic'[1,6]
|
||
|
|
scalar `stataaic' = `stataic'[1,5]
|
||
|
|
}
|
||
|
|
|
||
|
|
if ("`save'"!="") local saving 0
|
||
|
|
if ("`diff'"!="") local using 0
|
||
|
|
|
||
|
|
*-> get information from last model estimated
|
||
|
|
|
||
|
|
local wtis ""
|
||
|
|
local iswgt "no"
|
||
|
|
if "`e(wtype)'"!="" {
|
||
|
|
local iswgt "yes"
|
||
|
|
local wtis "[`e(wtype)'`e(wexp)']"
|
||
|
|
}
|
||
|
|
|
||
|
|
mat `v_e' = J(1,1,1)
|
||
|
|
sca `n_obs' = e(N)
|
||
|
|
sca `llF' = e(ll)
|
||
|
|
sca `llR' = e(ll_0)
|
||
|
|
sca `df' = e(df_m)
|
||
|
|
if "`cmd'"~="regress" & "`cmd'"~="tobit" & "`cmd'"~="cnreg" {
|
||
|
|
sca `x2' = e(chi2)
|
||
|
|
}
|
||
|
|
if "`cmd'"=="tobit" | "`cmd'"=="cnreg" {
|
||
|
|
sca `x2' = e(F)
|
||
|
|
}
|
||
|
|
if "`cmd'"=="regress" {
|
||
|
|
sca `ssr' = e(rss)
|
||
|
|
sca `r2' = e(r2)
|
||
|
|
sca `r2adj' = e(r2_a)
|
||
|
|
mat `v_e'[1,1] = (e(rmse)*e(rmse)*e(df_r)/e(N))
|
||
|
|
}
|
||
|
|
if "`cmd'"=="clogit" {
|
||
|
|
if "`e(wtype)'"=="iweight" |"`e(wtype)'"=="pweight" ///
|
||
|
|
| "`e(wtype)'"=="aweight" {
|
||
|
|
di in r "cannot use fitstat after clogit with `e(wtype)'"
|
||
|
|
exit
|
||
|
|
}
|
||
|
|
qui egen `uniqgrp' = rank(_n) if e(sample)==1, by(`e(group)')
|
||
|
|
qui replace `uniqgrp' = . if `uniqgrp' ~= 1
|
||
|
|
qui sum `uniqgrp' `wtis', meanonly
|
||
|
|
sca `n_obs' = r(sum_w)
|
||
|
|
}
|
||
|
|
|
||
|
|
*-> classify type of model
|
||
|
|
|
||
|
|
* binary?
|
||
|
|
local isbin "no"
|
||
|
|
if ("`cmd'"=="probit") local isbin "yes"
|
||
|
|
if ("`cmd'"=="logit") local isbin "yes"
|
||
|
|
if ("`cmd'"=="logistic") local isbin "yes"
|
||
|
|
if ("`cmd'"=="cloglog") local isbin "yes"
|
||
|
|
* ordinal?
|
||
|
|
local isord "no"
|
||
|
|
if ("`cmd'"=="oprobit") local isord "yes"
|
||
|
|
if ("`cmd'"=="ologit") local isord "yes"
|
||
|
|
* var(e) assumed?
|
||
|
|
local isvare "no"
|
||
|
|
if ("`cmd'"=="probit") local isvare "yes"
|
||
|
|
if ("`cmd'"=="logit") local isvare "yes"
|
||
|
|
if ("`cmd'"=="logistic") local isvare "yes"
|
||
|
|
if ("`cmd'"=="oprobit") local isvare "yes"
|
||
|
|
if ("`cmd'"=="ologit") local isvare "yes"
|
||
|
|
* is there a ystar?
|
||
|
|
local isystar "`isvare'"
|
||
|
|
if ("`cmd'"=="tobit") local isystar "yes"
|
||
|
|
if ("`cmd'"=="intreg") local isystar "yes"
|
||
|
|
if ("`cmd'"=="cnreg") local isystar "yes"
|
||
|
|
* logit model?
|
||
|
|
local islgt "no"
|
||
|
|
if ("`cmd'"=="logit") local islgt "yes"
|
||
|
|
if ("`cmd'"=="ologit") local islgt "yes"
|
||
|
|
if ("`cmd'"=="logistic") local islgt "yes"
|
||
|
|
* probit model?
|
||
|
|
local ispbt "no"
|
||
|
|
if ("`cmd'"=="probit") local ispbt "yes"
|
||
|
|
if ("`cmd'"=="oprobit") local ispbt "yes"
|
||
|
|
* zip/zinb?
|
||
|
|
local iszi "no"
|
||
|
|
if ("`cmd'"=="zip") local iszi "yes"
|
||
|
|
if ("`cmd'"=="zinb") local iszi "yes"
|
||
|
|
* provide pseudoR2 measures - jf 20Jun2005
|
||
|
|
local pseudo "yes"
|
||
|
|
if ("`cmd'"=="regress") local pseudo "no"
|
||
|
|
if ("`cmd'"=="rologit") local pseudo "no"
|
||
|
|
if ("`cmd'"=="slogit") local pseudo "no"
|
||
|
|
if ("`cmd'"=="mprobit") local pseudo "no"
|
||
|
|
* log-lik intercept only - jf 20Jun2005
|
||
|
|
local llio "yes"
|
||
|
|
if ("`cmd'"=="rologit") local llio "no"
|
||
|
|
if ("`cmd'"=="slogit") local llio "no"
|
||
|
|
if ("`cmd'"=="mprobit") local llio "no"
|
||
|
|
|
||
|
|
|
||
|
|
*-> compute numbers of variables and parameters
|
||
|
|
|
||
|
|
sca `n_parm' = colsof(`b')
|
||
|
|
local n_cat = e(k_cat)
|
||
|
|
* version 1.8.1 2007-09-18 - stata 10 changed returns for some commands
|
||
|
|
if c(stata_version) >= 10 {
|
||
|
|
local n_cat = e(k_out)
|
||
|
|
if ("`cmd'"=="ologit" | "`cmd'"=="oprobit") local n_cat = e(k_cat)
|
||
|
|
}
|
||
|
|
|
||
|
|
if "`n_cat'"=="." {
|
||
|
|
local n_cat = 2
|
||
|
|
}
|
||
|
|
local n_rhs = `n_parm' - `n_cat' + 1
|
||
|
|
if "`cmd'"=="mlogit" {
|
||
|
|
local n_rhs = (`n_parm'/(`n_cat'-1)) - 1
|
||
|
|
}
|
||
|
|
if "`cmd'"=="tobit" | "`cmd'"=="intreg" /// for var(e)
|
||
|
|
| "`cmd'"=="nbreg" | "`cmd'"=="cnreg" /// for alpha
|
||
|
|
| "`cmd'"=="zip" /// for 2nd intercept
|
||
|
|
{
|
||
|
|
local n_rhs = `n_rhs' - 1
|
||
|
|
}
|
||
|
|
if "`cmd'"=="zinb" {
|
||
|
|
local n_rhs = `n_rhs' - 2
|
||
|
|
} // 2nd int & alpha
|
||
|
|
|
||
|
|
local varlist : colnames(`b')
|
||
|
|
parse "`varlist'", parse (" ")
|
||
|
|
local rhsnam ""
|
||
|
|
if "`isord'"=="yes" {
|
||
|
|
local i 1
|
||
|
|
while "``i''"!= "" {
|
||
|
|
local j = substr("``i''",1,4)
|
||
|
|
* pre1.7.0 050216 if "`j'"!="_cut" {
|
||
|
|
* 050216 for Stata 9
|
||
|
|
if "`j'"!="_con" {
|
||
|
|
local rhsnam "`rhsnam' ``i''"
|
||
|
|
}
|
||
|
|
local i = `i' + 1
|
||
|
|
}
|
||
|
|
}
|
||
|
|
if "`isbin'"=="yes" | "`cmd'"=="mlogit" | "`cmd'"=="regress" /*
|
||
|
|
*/ | "`cmd'"=="tobit" | "`cmd'"=="intreg" | "`cmd'"=="cnreg" /*
|
||
|
|
*/ | "`cmd'"=="poisson" | "`cmd'"=="nbreg" {
|
||
|
|
local i 1
|
||
|
|
while "``i''"!= "" {
|
||
|
|
if "``i''"!="_cons" & "``i''"!="_se" {
|
||
|
|
local rhsnam "`rhsnam' ``i''"
|
||
|
|
}
|
||
|
|
local i = `i' + 1
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
*-> compute ll0 for zip and zinb (ll0 is for model with two intercepts)
|
||
|
|
|
||
|
|
if "`iszi'"=="yes" {
|
||
|
|
tempname res zill0
|
||
|
|
tempvar insamp
|
||
|
|
sca `df' = `n_parm' - 2
|
||
|
|
if "`cmd'"=="zinb" {
|
||
|
|
* remove alpha parameter
|
||
|
|
sca `df' = `n_parm' - 3
|
||
|
|
}
|
||
|
|
gen `insamp' = e(sample)
|
||
|
|
local inftyp "`e(inflate)'"
|
||
|
|
if "`inftyp'"=="logit" local inftyp ""
|
||
|
|
local doit "`e(cmd)' `e(depvar)' if `insamp' `wtis',inf(_cons) `inftyp'"
|
||
|
|
_estimates hold `res'
|
||
|
|
quietly `doit'
|
||
|
|
sca `llR' = e(ll) /* sic: not ll_0 */
|
||
|
|
_estimates unhold `res'
|
||
|
|
sca `x2' = -2*(`llR'-`llF')
|
||
|
|
} // iszi model
|
||
|
|
|
||
|
|
if "`bic'" == "" {
|
||
|
|
di _newline in g "Measures of Fit for " in y "`cmd'" in g " of " ///
|
||
|
|
in y "`depv'" _newline
|
||
|
|
}
|
||
|
|
|
||
|
|
if "`isbin'" == "yes" | "`isord'" == "yes" | "`cmd'" == "mlogit" ///
|
||
|
|
| "`cmd'" == "mprobit" {
|
||
|
|
|
||
|
|
* sum does not allow pweights or iweights
|
||
|
|
if "`e(wtype)'"=="pweight" | "`e(wtype)'"=="iweight" {
|
||
|
|
sca `cntr2a' = .
|
||
|
|
sca `cntr2' = .
|
||
|
|
if "`isbin'" == "yes" {
|
||
|
|
sca `r2ef' = .
|
||
|
|
}
|
||
|
|
if "`isbin'"=="yes" {
|
||
|
|
di _n in blu /*
|
||
|
|
*/ "(Efron's R2, Count R2, and Adj Count R2 not calculated if " /*
|
||
|
|
*/ in whi "`e(wtype)'" in blu " used)"
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
di _n in blu /*
|
||
|
|
*/ "(Count R2 and Adj Count R2 not calculated if " /*
|
||
|
|
*/ in whi "`e(wtype)'" in blu " used)"
|
||
|
|
}
|
||
|
|
}
|
||
|
|
* can compute measures if aweight, fweight or no weight
|
||
|
|
else quietly {
|
||
|
|
egen `depvgrp' = group(`depv') if e(sample)==1
|
||
|
|
sca `maxrow' = 0
|
||
|
|
gen `hiprob' = 0 if e(sample)==1
|
||
|
|
gen `predcat' = 0 if e(sample)==1
|
||
|
|
local counter = 1
|
||
|
|
|
||
|
|
* BRM - Efron R2 and Count R2s
|
||
|
|
if "`isbin'" == "yes" {
|
||
|
|
* Efron R2
|
||
|
|
* convert 0/!0 to 0/1 for binary outcomes
|
||
|
|
gen `depv2' = (`depv'!=0) if e(sample)==1
|
||
|
|
predict `phat' if e(sample)==1
|
||
|
|
sum `depv2' `wtis' if e(sample)==1
|
||
|
|
local mn = r(mean)
|
||
|
|
gen `yphat' = (`depv2' - `phat')^2 if e(sample)==1
|
||
|
|
gen `ydev' = (`depv2' - `mn')^2 if e(sample)==1
|
||
|
|
sum `yphat' `wtis' if e(sample)==1
|
||
|
|
local num = r(mean)
|
||
|
|
sum `ydev' `wtis' if e(sample)==1
|
||
|
|
local den = r(mean)
|
||
|
|
scalar `r2ef' = 1 - (`num'/`den')
|
||
|
|
* Count R2 and Adjusted Count R2
|
||
|
|
replace `predcat' = round(`phat', 1) if e(sample)==1
|
||
|
|
gen `wrong' = (`predcat'~=`depv2') if e(sample)==1
|
||
|
|
sum `depv2' `wtis' if `depv2' == 1 & e(sample)==1
|
||
|
|
sca `totone' = r(N)
|
||
|
|
if `totone' > (`n_obs'-`totone') {
|
||
|
|
sca `maxrow' = `totone'
|
||
|
|
}
|
||
|
|
else {
|
||
|
|
sca `maxrow' = `n_obs' - `totone'
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
* ORM - Count R2s
|
||
|
|
if "`isord'" == "yes" {
|
||
|
|
forvalues counter = 1(1)`n_cat' {
|
||
|
|
* ctnam is argument for ologitp/oprobitp
|
||
|
|
tempvar ct`counter'
|
||
|
|
local ctnam "`ctnam' `ct`counter''"
|
||
|
|
* find frequency of the modal category (for adj count R2)
|
||
|
|
sum `depvgrp' if (`depvgrp'==`counter') & (e(sample)==1) `wtis'
|
||
|
|
if r(N) > `maxrow' {
|
||
|
|
sca `maxrow' = r(N)
|
||
|
|
}
|
||
|
|
}
|
||
|
|
predict `ctnam' if e(sample) == 1
|
||
|
|
forvalues counter = 1(1)`n_cat' {
|
||
|
|
replace `predcat' = `counter' ///
|
||
|
|
if (`ct`counter''>`hiprob') & (e(sample)==1)
|
||
|
|
replace `hiprob' = `ct`counter'' ///
|
||
|
|
if (`ct`counter''>`hiprob') & (e(sample)==1)
|
||
|
|
}
|
||
|
|
gen `wrong' = (`predcat'~=`depvgrp') if e(sample)==1
|
||
|
|
}
|
||
|
|
|
||
|
|
* MNLM & and MPROBIT - Count R2s
|
||
|
|
|
||
|
|
if "`cmd'" == "mlogit" ///
|
||
|
|
| "`cmd'" == "mprobit" {
|
||
|
|
while `counter' <= `n_cat' {
|
||
|
|
tempvar ct`counter'
|
||
|
|
sum `depv' `wtis' ///
|
||
|
|
if (`depvgrp'==`counter') & (e(sample)==1), meanonly
|
||
|
|
local outcm = r(mean)
|
||
|
|
predict `ct`counter'', outcome(`outcm')
|
||
|
|
sum `depv' if (`depv'==`outcm') & (e(sample)==1) `wtis'
|
||
|
|
if r(N) > `maxrow' {
|
||
|
|
sca `maxrow' = r(N)
|
||
|
|
}
|
||
|
|
replace `predcat' = `outcm' ///
|
||
|
|
if (`ct`counter''>`hiprob') & (e(sample)==1)
|
||
|
|
replace `hiprob' = `ct`counter'' ///
|
||
|
|
if (`ct`counter''>`hiprob') & (e(sample)==1)
|
||
|
|
local counter = `counter' + 1
|
||
|
|
}
|
||
|
|
gen `wrong' = (`predcat'~=`depv') if e(sample)==1
|
||
|
|
}
|
||
|
|
* compute count R2s
|
||
|
|
sum `wrong' if `wrong' == 0 & e(sample) == 1 `wtis'
|
||
|
|
sca `totcor' = r(N)
|
||
|
|
sca `cntr2' = `totcor'/`n_obs'
|
||
|
|
sca `cntr2a' = (`totcor'-`maxrow')/(`n_obs'-`maxrow')
|
||
|
|
}
|
||
|
|
|
||
|
|
} // is binary, ordinal or nominal
|
||
|
|
|
||
|
|
* CLOGIT - Count R2s
|
||
|
|
if "`cmd'"=="clogit" {
|
||
|
|
quietly {
|
||
|
|
tempname cntr2 clphat clhi clpick clpred clwrong clright /*
|
||
|
|
*/ posval motest1 motest2 mptest1
|
||
|
|
* generate predicted outcomes for each group
|
||
|
|
predict `clphat' if e(sample) == 1
|
||
|
|
egen `clhi' = max(`clphat') if e(sample) == 1, by(`e(group)')
|
||
|
|
gen `clpred' = 1 if `clhi'==`clphat' & e(sample)==1
|
||
|
|
replace `clpred' = 0 if `clhi'!=`clphat' & e(sample)==1
|
||
|
|
* check if multiple predicted 1's per group
|
||
|
|
egen `mptest1' = sum(`clpred') if e(sample), by(`e(group)')
|
||
|
|
* weight of multiple predicted 1's per group split among each 1
|
||
|
|
replace `clpred' = `clpred'/`mptest1' if `clpred'==1 & e(sample)==1
|
||
|
|
* calculate count R2
|
||
|
|
sum `clpred' `wtis' if `clpred'>0 & `clpred'~=. & `depv'==0 & e(sample)==1
|
||
|
|
sca `clwrong' = r(sum)
|
||
|
|
sum `clpred' `wtis' if `clpred'>0 & `clpred'~=. & `depv'!=0 & e(sample)==1
|
||
|
|
sca `clright' = r(sum)
|
||
|
|
sca `cntr2' = `clright' / (`clwrong'+`clright')
|
||
|
|
* warn user if multiple outcomes per group
|
||
|
|
gen `posval'=1 if `depv'!=0 & `depv'~=. & e(sample)==1
|
||
|
|
egen `motest1' = sum(`posval') if `posval'==1, by(`e(group)')
|
||
|
|
sum `motest1'
|
||
|
|
sca `motest2' = r(max)
|
||
|
|
if `motest2' > 1 {
|
||
|
|
noi di in blu "(Multiple outcomes per group for clogit: Count R2 may not be valid)"
|
||
|
|
}
|
||
|
|
} // quietly
|
||
|
|
} // clogit
|
||
|
|
|
||
|
|
*-> compute Var(y*) using cov_x
|
||
|
|
|
||
|
|
if "`isystar'"=="yes" {
|
||
|
|
sca `factor' = 1/(`n_obs'-1)
|
||
|
|
quietly mat accum `v_x' = `depv' `rhsnam' `wtis' if e(sample), ///
|
||
|
|
deviations noconstant
|
||
|
|
sca `v_y' = `factor' * `v_x'[1,1] // sum y/(n-1)
|
||
|
|
mat `v_x' = `factor' * `v_x'[2...,2...] // cov matrix of rhs vars
|
||
|
|
|
||
|
|
if "`e(cmd)'"=="intreg" {
|
||
|
|
mat `v_x' = `v_x'[2...,2...]
|
||
|
|
}
|
||
|
|
mat `bnocon' = `b'[1,1..`n_rhs'] // trim off _con
|
||
|
|
mat `v' = `bnocon' * `v_x' * `bnocon''
|
||
|
|
if "`ispbt'"=="yes" {
|
||
|
|
mat `v_e'[1,1] = 1
|
||
|
|
}
|
||
|
|
if "`islgt'"=="yes" {
|
||
|
|
mat `v_e'[1,1] = (_pi*_pi)/3
|
||
|
|
}
|
||
|
|
if "`e(cmd)'"=="intreg" {
|
||
|
|
mat `v_e'[1,1] = e(sigma)*e(sigma)
|
||
|
|
}
|
||
|
|
if "`e(cmd)'"=="tobit" | "`e(cmd)'"=="cnreg" {
|
||
|
|
mat `v_e'[1,1] = `b'[1,"sigma:_cons"]*`b'[1,"sigma:_cons"]
|
||
|
|
}
|
||
|
|
mat `v' = `v' + `v_e'
|
||
|
|
sca `v_ystar' = `v'[1,1]
|
||
|
|
}
|
||
|
|
|
||
|
|
*-> information measures
|
||
|
|
|
||
|
|
sca `aicn' = ((-2*`llF') + (2*`n_parm'))
|
||
|
|
sca `aic' = `aicn'/`n_obs'
|
||
|
|
sca `dev' = -2*`llF'
|
||
|
|
local devdf = `n_obs'-`n_parm'
|
||
|
|
sca `lrx2' = 2*(`llF'-`llR')
|
||
|
|
local lrx2df = e(df_m)
|
||
|
|
if "`iszi'"=="yes" {
|
||
|
|
local lrx2df = `df'
|
||
|
|
}
|
||
|
|
if "`cmd'"=="mprobit" {
|
||
|
|
scalar `lrx2' = e(chi2)
|
||
|
|
}
|
||
|
|
|
||
|
|
sca `lrx2p' = chiprob(`lrx2df',`lrx2')
|
||
|
|
|
||
|
|
if "`cmd'"=="slogit" {
|
||
|
|
scalar `lrx2' = e(chi2)
|
||
|
|
sca `lrx2p' = e(p)
|
||
|
|
}
|
||
|
|
sca `bicx' = `dev' - (`devdf'*ln(`n_obs'))
|
||
|
|
sca `bicp' = -`lrx2' + (e(df_m)*ln(`n_obs'))
|
||
|
|
capture return scalar dev_df = `devdf'
|
||
|
|
capture return scalar lrx2_df = `lrx2df'
|
||
|
|
|
||
|
|
if "`iszi'"=="yes" {
|
||
|
|
* df = # rhs in both equations
|
||
|
|
sca `bicp' = -`lrx2' + (`df'*ln(e(N)))
|
||
|
|
}
|
||
|
|
|
||
|
|
* these equal usual R2 in regress
|
||
|
|
if "`pseudo'"=="yes" {
|
||
|
|
sca `r2mcf' = 1 - (`llF'/`llR')
|
||
|
|
sca `r2mcfa' = 1 - ((`llF'-`n_parm')/`llR')
|
||
|
|
sca `n_obs2' = 2/`n_obs'
|
||
|
|
sca `r2ml' = 1 - exp(2*(`llR'-`llF')/`n_obs')
|
||
|
|
sca `r2cu' = (`r2ml')/(1-exp(2*`llR'/`n_obs'))
|
||
|
|
}
|
||
|
|
|
||
|
|
if "`isystar'"=="yes" {
|
||
|
|
sca `r2mz' = (`v_ystar' - `v_e'[1,1])/`v_ystar'
|
||
|
|
sca `v_error' = `v_e'[1,1]
|
||
|
|
}
|
||
|
|
|
||
|
|
*-> output
|
||
|
|
|
||
|
|
* define output macros with format where / is parse character:
|
||
|
|
* output_text/tempname_of_scalar/abbrev_for_r()
|
||
|
|
* /df?/width/half(1, 2 or ? for either)
|
||
|
|
local stat2 "Log-Lik Intercept Only/llR/ll_0/no/11/1"
|
||
|
|
local stat3 "Log-Lik Full Model/llF/ll/no/11/2"
|
||
|
|
local stat4 "D/dev/dev/yes/11/1"
|
||
|
|
local stat6 "LR/lrx2/lrx2/yes/11/2"
|
||
|
|
local stat8 "Prob > LR/lrx2p/lrx2_p/no/9/2"
|
||
|
|
local stat9 "McFadden's R2/r2mcf/r2_mf/no/9/1"
|
||
|
|
local stat10 "McFadden's Adj R2/r2mcfa/r2_mfadj/no/9/2"
|
||
|
|
local stat11 "ML (Cox-Snell) R2/r2ml/r2_ml/no/9/1"
|
||
|
|
local stat12 "Cragg-Uhler(Nagelkerke) R2/r2cu/r2_cu/no/7/2"
|
||
|
|
local stat13 "McKelvey & Zavoina's R2/r2mz/r2_mz/no/9/?"
|
||
|
|
local stat14 "Efron's R2/r2ef/r2_ef/no/9/?"
|
||
|
|
local stat15 "Variance of y*/v_ystar/v_ystar/no/11/1"
|
||
|
|
local stat16 "Variance of error/v_error/v_error/no/11/2"
|
||
|
|
local stat17 "Count R2/cntr2/r2_ct/no/9/1"
|
||
|
|
local stat18 "Adj Count R2/cntr2a/r2_ctadj/no/9/2"
|
||
|
|
local stat19 "R2/r2/r2/no/9/1"
|
||
|
|
local stat20 "Adjusted R2/r2adj/r2_adj/no/9/2"
|
||
|
|
local stat21 "AIC/aic/aic/no/11/1"
|
||
|
|
local stat22 "AIC*n/aicn/aic_n/no/11/2"
|
||
|
|
local stat23 "BIC/bicx/bic/no/11/1"
|
||
|
|
local stat24 "BIC'/bicp/bic_p/no/11/2"
|
||
|
|
local stat25 "BIC used by Stata/statabic/statabic/no/11/1"
|
||
|
|
local stat26 "AIC used by Stata/stataaic/stataaic/no/11/2"
|
||
|
|
if "`llio'"=="no" {
|
||
|
|
local stat2 ""
|
||
|
|
local stat3 "Log-Lik Full Model/llF/ll/no/11/1"
|
||
|
|
local stat4 "D/dev/dev/yes/11/2"
|
||
|
|
}
|
||
|
|
if "`cmd'"=="slogit" local stat24 ""
|
||
|
|
if "`cmd'"=="mprobit" | "`cmd'"=="slogit" {
|
||
|
|
local stat6 "Wald X2/lrx2/lrx2/yes/11/1"
|
||
|
|
}
|
||
|
|
if "`cmd'"=="mprobit" | "`cmd'"=="slogit" {
|
||
|
|
local stat8 "Prob > X2/lrx2p/lrx2_p/no/9/2"
|
||
|
|
}
|
||
|
|
|
||
|
|
*-> check models for using and saving
|
||
|
|
|
||
|
|
if "`saving'" != "" {
|
||
|
|
if (length("`saving'")>4) {
|
||
|
|
di in red "saving() name must be < 5 characters long"
|
||
|
|
exit 198
|
||
|
|
}
|
||
|
|
}
|
||
|
|
if "`using'"!= "" {
|
||
|
|
if (length("`using'")>4) {
|
||
|
|
di in red "using() name must be < 5 characters long"
|
||
|
|
exit 198
|
||
|
|
}
|
||
|
|
qui capture mat list fs_`using'
|
||
|
|
if _rc == 111 {
|
||
|
|
di in red "Incorrect using() option: no indices saved as `using'"
|
||
|
|
exit 111
|
||
|
|
}
|
||
|
|
* if the using option has been set, prepare comparison macros
|
||
|
|
sca `fn_obs' = fs_`using'[1,1]
|
||
|
|
if `fn_obs' == -9999 {
|
||
|
|
sca `fn_obs' = .
|
||
|
|
}
|
||
|
|
|
||
|
|
* estimate command saved as row name
|
||
|
|
local fcmd : rownames(fs_`using')
|
||
|
|
local fcmd : word 1 of `fcmd'
|
||
|
|
|
||
|
|
* test if current and using are same type of model
|
||
|
|
local same "no"
|
||
|
|
if ("`fcmd'"=="logit"&"`e(cmd)'"=="logistic") | ///
|
||
|
|
("`fcmd'"=="logistic"&"`e(cmd)'"=="logit") {
|
||
|
|
local same "yes"
|
||
|
|
}
|
||
|
|
if ("`fcmd'"!="`e(cmd)'" & "`same'"=="no") & "`force'"=="" {
|
||
|
|
di in r "Current model estimated by `e(cmd)'" ///
|
||
|
|
in r ", but saved model estimated by `fcmd'."
|
||
|
|
di in r "To make the comparisons anyway, use the " ///
|
||
|
|
"force option."
|
||
|
|
exit
|
||
|
|
}
|
||
|
|
|
||
|
|
if ("`fcmd'"!="`e(cmd)'" & "`same'"=="no") {
|
||
|
|
di in y "Warning: Current model estimated by " /*
|
||
|
|
*/ "`e(cmd)', but saved model estimated by `fcmd'"
|
||
|
|
di
|
||
|
|
}
|
||
|
|
if `n_obs' != `fn_obs' & "`force'"=="force" {
|
||
|
|
di in y "Warning: N's do not match."
|
||
|
|
di
|
||
|
|
}
|
||
|
|
|
||
|
|
di in g _col(30) %9s "Current" _col(48) " Saved" ///
|
||
|
|
_col(65) "Difference"
|
||
|
|
di in g "Model:" _col(30) %9s in y "`e(cmd)'" ///
|
||
|
|
%9s _col(48) in y "`fcmd'"
|
||
|
|
di in g "N:" _col(30) %9.0f in y `n_obs' ///
|
||
|
|
%9.0f _col(48) in y `fn_obs' ///
|
||
|
|
%9.0f _col(66) in y (`n_obs' - `fn_obs')
|
||
|
|
|
||
|
|
* test if N's match
|
||
|
|
if `n_obs' != `fn_obs' & "`force'"=="" {
|
||
|
|
di _newline in r "N's do not match. To make the comparisons, use the " ///
|
||
|
|
in y "force" in r " option."
|
||
|
|
exit
|
||
|
|
}
|
||
|
|
|
||
|
|
* test of ll0 match
|
||
|
|
sca `fllR' = fs_`using'[1,2]
|
||
|
|
if `fllR'==-9999 {
|
||
|
|
sca `fllR' = .
|
||
|
|
}
|
||
|
|
|
||
|
|
if abs(`llR'-`fllR') > .0001 & (`llR'-`fllR') ~= . & "`force'"=="" {
|
||
|
|
di _newline in r /*
|
||
|
|
*/ "Log-likelihoods for intercept model differ for current and saved models."
|
||
|
|
di in r "To make the comparisons, use the " in y "force" in r " option."
|
||
|
|
exit
|
||
|
|
}
|
||
|
|
} /* if "`using'"!="" */
|
||
|
|
|
||
|
|
capture return scalar n_parm = `n_parm'
|
||
|
|
capture return scalar n_rhs = `n_rhs'
|
||
|
|
capture return scalar N = `n_obs'
|
||
|
|
|
||
|
|
*-> create matrices and macros containing results
|
||
|
|
|
||
|
|
local nstats = 24 // number of stats reported in output matrix
|
||
|
|
if "`nostata'"!="stata" local nstats = `nstats' + 2
|
||
|
|
|
||
|
|
mat `results' = J(1, `nstats'+2, .)
|
||
|
|
mat `results'[1,1] = `n_obs'
|
||
|
|
mat `results'[1,`nstats'+1] = `n_parm'
|
||
|
|
mat `results'[1,`nstats'+2] = `n_rhs'
|
||
|
|
local cnames "N" // start bulding column names for output matrix
|
||
|
|
|
||
|
|
forvalues i = 2(1)`nstats' {
|
||
|
|
|
||
|
|
* read individual stats
|
||
|
|
tokenize "`stat`i''", parse(/)
|
||
|
|
* outtxt - name of statistic
|
||
|
|
local outtxt`i' "`1'"
|
||
|
|
* scalnm - name of scalar holding value of statistic
|
||
|
|
local scalnm`i' "`3'"
|
||
|
|
* abnm - name of statistic for return value
|
||
|
|
local abnm`i' "`5'"
|
||
|
|
* df - whether statistic has df
|
||
|
|
local df`i' "`7'"
|
||
|
|
* outwid - width of output - TO DO: try to delete this
|
||
|
|
local outwid`i' "`9'"
|
||
|
|
* outhalf - need to be in first or second half
|
||
|
|
local outhalf`i' "`11'"
|
||
|
|
* cnames - macro of column names for output
|
||
|
|
local cnames "`cnames' `abnm`i''"
|
||
|
|
if "`df`i''"=="yes" {
|
||
|
|
local cnames "`cnames' `abnm`i''_df"
|
||
|
|
}
|
||
|
|
|
||
|
|
* see if scalar exists, if not, no output
|
||
|
|
capture return scalar `abnm`i'' = ``scalnm`i'''
|
||
|
|
if _rc != 0 {
|
||
|
|
if "`outtxt`i''" != "" {
|
||
|
|
mat `results'[1, `i'] = .x
|
||
|
|
if "`df`i''"=="yes" {
|
||
|
|
mat `results'[1,`i'+1] = .x
|
||
|
|
}
|
||
|
|
}
|
||
|
|
}
|
||
|
|
if _rc == 0 {
|
||
|
|
* add to results matrix
|
||
|
|
mat `results'[1, `i'] = ``scalnm`i'''
|
||
|
|
if "`df`i''"=="yes" {
|
||
|
|
mat `results'[1,`i'+1] = ``scalnm`i''df'
|
||
|
|
}
|
||
|
|
|
||
|
|
} // if _rc == 0
|
||
|
|
if "`df`i''"=="yes" {
|
||
|
|
local i = `i' + 1
|
||
|
|
}
|
||
|
|
} // forvalues i
|
||
|
|
|
||
|
|
|
||
|
|
*-> loop through to print results if using option not used
|
||
|
|
|
||
|
|
if "`using'"=="" {
|
||
|
|
|
||
|
|
* arrange numbers in rows and columns for printing
|
||
|
|
|
||
|
|
local trow = 1
|
||
|
|
local tcol = 1
|
||
|
|
forvalues i = 2(1)`nstats' {
|
||
|
|
* figure out if result should be in output
|
||
|
|
if `results'[1,`i']!=.x & "`outtxt`i''" != "" & ///
|
||
|
|
("`bic'"=="" | `i' >= 21) {
|
||
|
|
|
||
|
|
* figure out if need to advance to next row of table
|
||
|
|
if "`outhalf`i''" == "1" & `tcol' == 3 {
|
||
|
|
local ++trow
|
||
|
|
local tcol = 1
|
||
|
|
}
|
||
|
|
if "`outhalf`i''" == "2" & `tcol' == 1 {
|
||
|
|
local tcol = 3
|
||
|
|
}
|
||
|
|
|
||
|
|
local op`trow'_`tcol' = "`outtxt`i''"
|
||
|
|
if "`df`i''" == "yes" {
|
||
|
|
local temp = ``scalnm`i''df'
|
||
|
|
local op`trow'_`tcol' = "`op`trow'_`tcol''(`temp')"
|
||
|
|
}
|
||
|
|
local op`trow'_`tcol' = "`op`trow'_`tcol'':"
|
||
|
|
|
||
|
|
local ++tcol
|
||
|
|
local op`trow'_`tcol' = ``scalnm`i'''
|
||
|
|
|
||
|
|
* take care of output width
|
||
|
|
local fmt`trow'_`tcol' = "%`outwid`i''.3f"
|
||
|
|
|
||
|
|
// move to next column
|
||
|
|
local ++tcol
|
||
|
|
if `tcol' == 5 { // if needed, move to next row
|
||
|
|
local tcol = 1
|
||
|
|
local ++trow
|
||
|
|
}
|
||
|
|
}
|
||
|
|
} // forvalues i
|
||
|
|
|
||
|
|
* print rows and columns
|
||
|
|
|
||
|
|
local numrows = `trow'
|
||
|
|
if `tcol' == 1 {
|
||
|
|
local numrows = `trow' - 1
|
||
|
|
}
|
||
|
|
|
||
|
|
forvalues i = 1(1)`numrows' {
|
||
|
|
di in g "`op`i'_1'" _col(28) in y %11.3f `op`i'_2' ///
|
||
|
|
_col(42) in g "`op`i'_3'" _col(69) in y %11.3f `op`i'_4'
|
||
|
|
}
|
||
|
|
|
||
|
|
} // using option not used
|
||
|
|
|
||
|
|
|
||
|
|
*-> loop through to print results if using option has been used
|
||
|
|
|
||
|
|
if "`using'"!="" {
|
||
|
|
|
||
|
|
local trow = 1
|
||
|
|
forvalues i = 2(1)`nstats' {
|
||
|
|
* figure out if result should be in output
|
||
|
|
if `results'[1,`i']!=.x & ("`bic'"=="" | `i' >= 21) {
|
||
|
|
|
||
|
|
* calculate differences
|
||
|
|
sca `prevval' = fs_`using'[1,`i']
|
||
|
|
tempname difvals
|
||
|
|
sca `difvals' = ``scalnm`i'''-`prevval'
|
||
|
|
|
||
|
|
* for LR prob give LR prob -- computed during LR loop
|
||
|
|
if "`outtxt`i''"=="Prob > LR" {
|
||
|
|
sca `difvals' = `LRchip'
|
||
|
|
}
|
||
|
|
|
||
|
|
* print results
|
||
|
|
if "`df`i''" == "no" {
|
||
|
|
di in g "`outtxt`i''" ///
|
||
|
|
_col(28) in y %11.3f ``scalnm`i''' ///
|
||
|
|
_col(46) in y %11.3f `prevval' ///
|
||
|
|
_col(64) in y %11.3f `difvals'
|
||
|
|
}
|
||
|
|
|
||
|
|
if "`df`i''" == "yes" {
|
||
|
|
local tmp_curdf = ``scalnm`i''df'
|
||
|
|
sca `prevdf' = fs_`using'[1,(`i'+ 1)]
|
||
|
|
local tmp_prevdf = `prevdf'
|
||
|
|
local tmp_difdf = abs(`tmp_curdf' - `tmp_prevdf')
|
||
|
|
|
||
|
|
* for D remove the negative sign from the difference
|
||
|
|
if "`outtxt`i''"=="D" {
|
||
|
|
sca `difvals' = abs(`difvals')
|
||
|
|
}
|
||
|
|
|
||
|
|
* for LR statistic compute LR test prob
|
||
|
|
if "`outtxt`i''"=="LR" | "`outtxt`i''"=="Wald LR" {
|
||
|
|
tempname LRdif LRdf LRchip
|
||
|
|
sca `difvals' = abs(``scalnm`i''' - `prevval')
|
||
|
|
sca `LRdf' = abs(``scalnm`i''df'-`prevdf')
|
||
|
|
sca `LRchip' = chiprob(`LRdf', `difvals')
|
||
|
|
}
|
||
|
|
|
||
|
|
di in g "`outtxt`i''" ///
|
||
|
|
_col(28) in y %11.3f ``scalnm`i''' in g "(`tmp_curdf')" ///
|
||
|
|
_col(46) in y %11.3f `prevval' in g "(`tmp_prevdf')" ///
|
||
|
|
_col(64) in y %11.3f `difvals' in g "(`tmp_difdf')"
|
||
|
|
}
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
} // forvalues i
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
*-> BIC comparisons
|
||
|
|
if "`using'" != "" {
|
||
|
|
if abs(`llR'-`fllR') < .0001 { /* check that LL0s are the same */
|
||
|
|
sca `fbicp' = fs_`using'[1, 24]
|
||
|
|
sca `bicdif' = abs(`bicp'-`fbicp')
|
||
|
|
if `bicdif'~=. {
|
||
|
|
local sup "very strong"
|
||
|
|
if `bicdif'<= .0000000001 {
|
||
|
|
local sup "no"
|
||
|
|
}
|
||
|
|
else if `bicdif' <=2 {
|
||
|
|
local sup "weak"
|
||
|
|
}
|
||
|
|
else if `bicdif' <=6 {
|
||
|
|
local sup "positive"
|
||
|
|
}
|
||
|
|
else if `bicdif' <=10 {
|
||
|
|
local sup "strong"
|
||
|
|
}
|
||
|
|
local modpref "current"
|
||
|
|
if `bicp' > `fbicp' {
|
||
|
|
local modpref "saved"
|
||
|
|
}
|
||
|
|
if `bicdif'< .0000000001 & `bicdif'>-.0000000001 /*
|
||
|
|
*/ {
|
||
|
|
local modpref "either"
|
||
|
|
}
|
||
|
|
di _newline in g "Difference of " %8.3f in y `bicdif' in g /*
|
||
|
|
*/ " in BIC' provides " in y "`sup'" /*
|
||
|
|
*/ in g " support for " in y "`modpref'" in g " model."
|
||
|
|
}
|
||
|
|
} /* if `llR' == `fllR' */
|
||
|
|
} /* if "`using'" != "" */
|
||
|
|
|
||
|
|
if "`bic'"=="" & "`using'" != "" {
|
||
|
|
di _n in g "Note: p-value for difference in LR is only valid if models are nested."
|
||
|
|
}
|
||
|
|
local cnames "`cnames' n_parm n_rhs"
|
||
|
|
mat rownames `results' = `e(cmd)'
|
||
|
|
mat colnames `results' = `cnames'
|
||
|
|
if "`saving'" != "" {
|
||
|
|
capture mat drop fs_`saving'
|
||
|
|
mat fs_`saving' = `results'
|
||
|
|
di _n in blu "(Indices saved in matrix fs_`saving')"
|
||
|
|
}
|
||
|
|
|
||
|
|
quietly `cmd'
|
||
|
|
end
|
||
|
|
|
||
|
|
exit
|
||
|
|
* version 1.8.1 2007-09-18 returns for some commands changed in Stata 10
|
||
|
|
* version 1.8.0 19Jun2006 df correction in lrtest for zinb
|
||
|
|
* version 1.7.9 26Oct2005 fix slogit and mprobit output for LR to Wald X2
|
||
|
|
* version 1.7.8b 20Jun2005 small fix to new output routine
|
||
|
|
* version 1.7.8 20Jun2005 new output routine
|
||
|
|
* version 1.7.7 17Jun2005 streamlined/stata9 improvements to code
|
||
|
|
* version 1.7.6 15Jun2005 added rologit
|
||
|
|
* version 1.7.5 10Jun2005 clogit, tobit, cnreg, Cragg-Uhler (Nagelkerke) R2 for dif, ztp, ztnb
|
||
|
|
* version 1.7.4 10Jun2005 align Variance of y*
|
||
|
|
* version 1.7.2 13Apr2005
|
||
|
|
* version 1.7.1 28Feb2005 mprobit; fixed bic'; added count R2
|
||
|
|
* version 1.7.0 16Feb2005 slogit mprobit ztp, ztnb, format on others
|
||
|
|
* version 1.6.9 12Feb2005 add stata bic and aic
|
||
|
|
* version 1.6.8 07Feb2005 add ztp and ztnb
|
||
|
|
* version 1.6.7 22Jan2005 change LRX2
|