📄 rocfit_8.ado
字号:
*! version 7.0.15 17mar2005
prog def rocfit_8, eclass
version 6, missing
local vv : display "version " string(_caller()) ", missing:"
if replay() {
if `"`e(cmd)'"' != "rocfit" { error 301 }
syntax [, Level(cilevel) *]
}
else {
syntax varlist(numeric min=2 max=2) [if] [in] /*
*/ [ fweight] [, Level(cilevel) /*
*/ FROM(string) MLMethod(string) /*
*/ MLOpt(string) noLOG CONTinuous(string) Generate(string) * ]
marksample touse
tempvar sample
gen int `sample'=`touse'
tokenize `varlist'
local D = `"`1'"'
local C = `"`2'"'
local ind `"`2'"'
local DC = `"`2'"'
if "`from'"!="" { local iniopt init(`from') }
if "`mlmethod'" == "" { local mlmetho "lf" }
if "`offset'" !="" { local offopt offset(`offset') }
mlopts options, `options'
cap assert `D'==0 | `D'==1 if `touse'
if _rc~=0 {
noi di in red "true status variable `D' must be 0 or 1"
exit 198
}
qui summ `D' if `touse', meanonly
if r(min) == r(max) {
di in red "outcome does not vary"
exit 198
}
if "`log'"!="" {
local nlog="*"
}
if "`continuous'"~="" & "`continuous'"~="." {
if `continuous'>_N {
di in red /*
*/ "continuous() must be less than" /*
*/ " or equal to the number of" /*
*/ " observations"
error 198
}
if `continuous' < 3 {
di in red /*
*/ "continuous() must be greater than" /*
*/ " or equal to 3"
error 198
}
tempvar pct
qui sum `C', meanonly
local cmax=r(max)
local cmin=r(min)
local A=`continuous'/(`cmax'+1-`cmin')
local I=1-`A'*`cmin'
qui gen long `pct'=int(`A'*`C'+`I')
tempvar C
qui gen `C'=`pct'
qui drop `pct'
}
if "`generate'"~="" {
if "`continuous'"=="" {
noi di in red /*
*/ "generate() cannot be specified without "/*
*/ "continuous() option"
error 198
}
if "`continuous'"=="." {
noi di in red /*
*/ "generate() cannot be specified with " /*
*/ "continuous(.) option; "
noi di in red "continuous(.) requests that "/*
*/ "class var not be grouped"
error 198
}
qui gen long `generate'=`C'
qui compress `generate'
}
/* BEGIN-create freq weighted summary data */
if `"`weight'"' =="" {
tempvar wv
qui gen int `wv' = 1 if `touse'
local nowt 1
local weight="fweight"
}
else {
tempvar wv
qui gen double `wv' `exp'
local nowt 0
}
local mymat: set matsize
qui inspect `C' if `touse'
if r(N_unique)>20 & "`continuous'"=="" {
di in red /*
*/ "variable " abbrev("`C'",12) /*
*/ " has over 20 unique values; " /*
*/ "use continuous(#) to group "
di in red /*
*/ "into # categories " /*
*/ "or continuous(.) if variable is already categorical"
exit 198
}
if r(N_unique)<. & r(N_unique)>`mymat'-3 {
if "`continuous'"=="" {
di in red /*
*/ "matsize too small; type -help matsize-"
di in red /*
*/ "you may want to specify continuous() option"
exit 908
}
else {
di in red /*
*/ "matsize needs to be at least: " /*
*/ `r(N_unique)'+3
error 908
}
}
if r(N_unique)>=. {
tempvar cnt
sort `touse' `C'
qui by `touse' `C': gen byte `cnt' = 1 /*
*/ if _n==1 & `touse'
qui summ `cnt', meanonly
if r(sum)>800 {
di in red /*
*/ "`C' takes on more than 800 unique values"
di in red "unable to estimate model;" /*
*/ " use roctab command"
exit 134
}
drop `cnt'
}
tempvar newwt
sort `touse' `D' `C'
qui by `touse' `D' `C': gen double `newwt'=sum(`wv') if `touse'
qui by `touse' `D' `C': replace `newwt'=. if _n~=_N
qui replace `newwt'=. if `newwt'==0
qui replace `touse'=0 if `newwt'>=.
qui replace `wv'=`newwt'
local wt=`"[fweight=`wv']"'
tempvar MN
sort `touse' `D' `C'
qui by `touse' `D' : gen long `MN' = sum(`wv')
qui by `touse' `D' : replace `MN' = `MN'[_N]
qui replace `MN'=. if `touse'==0
drop `newwt'
*preserve
*qui keep if `touse'
/* END -create freq weighted summary data */
/* BEGIN-get initial values */
tempvar p
qui logistic `D' `C' `wt' if `touse', asis
tempvar b
mat `b'=e(b)
local bc=colsof(`b')
if `bc'>1 {
local coef=_b[`C']
}
else {
di in green "unable to estimate model;" /*
*/ " Use roctab command"
est clear
exit
}
qui predict double `p', p
tempvar prob sens spec
qui sum `p', meanonly
cap assert reldif(`p', r(min))<1e-12 if `p'<.
if _rc==0 {
_rocsen if `touse', class(`C') genprob(`prob') /*
*/ gensens(`sens') genspec(`spec')
}
else {
qui `vv' lsens if `touse', nograph genprob(`prob') /*
*/ gensens(`sens') genspec(`spec')
}
qui replace `sens'=. if `prob'==0 | `prob'==1
qui replace `spec'=. if `prob'==0 | `prob'==1
qui replace `prob'=. if `sens'>=. | `spec'>=.
sort `prob'
qui replace `prob'=. if `sens'==1 & `spec'==0 & _n==1
tempvar g
qui egen `g'=group(`C') if `touse'
qui sum `g' , meanonly
local max=r(max)
local i= 1
while `i'<`max' {
local zv "`zv' /z`i'"
local nzv "`nzv' z`i'"
local zname "`zname' Z:_cut`i'"
local i=`i'+1
}
*noi di in red "`zv'"
tempvar zfpf ztpf ztnf
qui gen double `ztpf'=invnorm(`sens') if `touse'
qui gen double `zfpf'=invnorm(1-`spec') if `touse'
qui gen double `ztnf'=invnorm(`spec') if `touse'
if `max'>2 {
cap qui regress `ztpf' `zfpf'
if _rc==0 {
local a0 = _b[_cons]
local b0 = _b[`zfpf']
}
else {
local a0 = 1
local b0 = 1
}
}
else {
local a0 = 1
local b0 = 1
}
drop `ztpf' `zfpf' `sens' `spec'
tempname bc0
qui mat `bc0' = ( `a0', `b0')
qui mat colnames `bc0' = a:_cons b:_cons
tempname ZV
sort `prob'
local i= 1
while `i'<`max' {
tempname Z
if `ztnf'[`i'] <. {
cap mat `Z'=(`ztnf'[`i']+`i'*1e-6)
}
cap mat colnames `Z' = z`i':_cons
if `i'==1 {
cap mat `ZV'= `Z'
}
else {
cap mat `ZV'=`ZV' , `Z'
}
if _rc~=0 {
local i=`max'+1
local this 1
}
else { local i=`i'+1 }
}
drop `ztnf'
if "`this'"~="" {
sort `prob'
local i= 1
while `i'<`max' {
if `prob'[`i']>=. {
qui replace `prob'=0.9 if _n==`i'
}
tempname Z
qui mat `Z'=(invnorm(`prob'[`i']))
qui mat colnames `Z' = z`i':_cons
qui mat `bc0'=`bc0' , `Z'
local i=`i'+1
}
}
else {
qui mat `bc0'=`bc0' , `ZV'
}
if `coef'<0 {
local tel=colsof(`bc0')
tempname telm
qui mat `telm'=(-`bc0'[1,1], `bc0'[1,2])
qui mat colnames `telm' = a:_cons b:_cons
local j 1
local i=`tel'
while `i' >2 {
tempname telc
mat `telc' = -`bc0'[1,`i']
qui mat colnames `telc' = z`j':_cons
local i=`i'-1
local j=`j'+1
mat `telm'= `telm', `telc'
}
mat `bc0'= `telm'
}
/* noi mat l `bc0' */
`nlog' di _n in gr "Fitting binormal model:"
global MD `D'
global MC `g'
global ZV `nzv'
tempvar p
qui gen double `p'=.
global MP `p'
sort `D' `C'
ml model `mlmethod' rocf_lf /*
*/ (a: `D' `C' = ) /*
*/ /b `zv' `wt' if `touse', /*
*/ `robust' `cluopt' `scopt' init(`bc0') `mlopt' /*
*/ missing collin nopreserve /*
*/ max search(off) `log' `options'
global MD
global MC
global ZV
tempname B
qui mat `B'=e(b)
local a = `B'[1,1]
local b = `B'[1,2]
mat colnames `B' = :intercept :slope `zname'
tempname x2 df
tempvar num
qui gen double `num'= /*
*/ (`MN'/$MP)*(($MP - `wv'/`MN')^2) if `touse'
qui replace `num'=sum(`num')
scalar `x2'=`num'[_N]
scalar `df'=`max'-3
est scalar chi2_gf=`x2'
est scalar df_gf=`df'
est scalar p_gf=chiprob(`df',`x2')
est local chi2type="GOF"
est repost b=`B' , rename esample(`sample')
if `nowt'==1 {
est local wtype
est local wexp
}
else {
if "`e(wtype)'" != "" {
est local wexp `"`exp'"'
}
}
ROCInd `level'
est local title = "Binormal model of "+abbrev("`D'",12)+ /*
*/ " on " + abbrev("`ind'",12)
est local predict _no_predict
est local depvar "`D' `DC'"
est local cmd rocfit
est scalar area =`s(area)'
est scalar se_area =`s(se_area)'
est scalar deltam =`s(deltam)'
est scalar se_delm =`s(se_delm)'
est scalar de =`s(de)'
est scalar se_de =`s(se_de)'
est scalar da =`s(da)'
est scalar se_da =`s(se_da)'
est scalar neg=`MN'
est scalar pos=e(N)-`MN'
global MP
}
di _n in gr `"`e(title)'"' _col(51) in gr "Number of obs" /*
*/ _col(67) "= " in ye %10.0g e(N)
di in gr "Goodness-of-fit chi2(" in ye "`e(df_gf)'" /*
*/ in gr ") = " in ye %10.2f e(chi2_gf)
di in gr "Prob > chi2 = " in ye %10.4f e(p_gf)
di in gr "Log likelihood = " in ye %10.0g e(ll)
DISTab1 `level'
DISTab2 `level'
end
program define ROCInd, sclass
/* program to calculated indices */
args level
tempname B V Z a b va vb vab
mat `B'=e(b)
mat `V'=e(V)
scalar `a' = `B'[1,1]
scalar `b' = `B'[1,2]
scalar `va'=`V'[1,1]
scalar `vb'=`V'[2,2]
scalar `vab'=`V'[1,2]
scalar `Z'=`a'/sqrt(1+`b'^2)
local alpha=invnorm(.5+`level'/200)
/* AREA */
tempname ga gb ase area delta de da
scalar `area'=normprob(`Z')
scalar `ga'=normd(`Z')*(1/sqrt(1+`b'^2))
scalar `gb'=normd(`Z')*(-`a'*`b')/((1+`b'^2)^(3/2))
scalar `ase'=(`ga'^2)*(`va')+(`gb'^2)*(`vb') + 2*`ga'*`gb'*`vab'
scalar `ase'=sqrt(`ase')
sret local se_area = `ase'
sret local area = `area'
/* DELTA M */
scalar `delta' = `a'/`b'
scalar `ga'=1/`b'
scalar `gb'=-`a'/(`b'^2)
scalar `ase'=(`ga'^2)*(`va')+(`gb'^2)*(`vb') + 2*`ga'*`gb'*`vab'
scalar `ase'=sqrt(`ase')
sret local se_delm = `ase'
sret local deltam = `delta'
/* DE Line */
scalar `de' = 2*`a'/(`b'+1)
scalar `ga'=2/(`b'+1)
scalar `gb'=-2*`a'/((`b'+1)^2)
scalar `ase'=(`ga'^2)*(`va')+(`gb'^2)*(`vb') + 2*`ga'*`gb'*`vab'
scalar `ase'=sqrt(`ase')
sret local se_de = `ase'
sret local de = `de'
/* DA Line */
scalar `da' = sqrt(2)*invnorm(`area')
scalar `ga'=sqrt(2)/sqrt(1+`b'^2)
scalar `gb'=-sqrt(2)*`a'*`b'/sqrt((1+`b'^2)^3)
scalar `ase'=(`ga'^2)*(`va')+(`gb'^2)*(`vb') + 2*`ga'*`gb'*`vab'
scalar `ase'=sqrt(`ase')
sret local se_da = `ase'
sret local da = `da'
end
prog def DISTab1
args level max
tempname B V
qui mat `B'=e(b)
qui mat `V'=e(V)
local max=colsof(`B')
local alpha=invnorm(.5+`level'/200)
di
DivLine top
di in smcl in gr " {c |} Coef. Std. Err." /*
*/ `" z P>|z| [`=strsubdp("`level'")'% Conf. Interval]"'
DivLine mid
ELine "intercept" `B'[1,1] sqrt(`V'[1,1]) `alpha' 0
ELine "slope (*)" `B'[1,2] sqrt(`V'[2,2]) `alpha' 1
DivLine mid
local i 3
while `i'<=`max' {
local j = `i'-2
ELine "/cut`j'" `B'[1,`i'] sqrt(`V'[`i',`i']) `alpha' 0
local i=`i'+1
}
DivLine bot
DivLine top
end
program define DivLine
args where
if "`where'" == "top" {
di in smcl in gr "{hline 13}{c TT}{hline 64}"
}
else if "`where'"=="mid" {
di in smcl in gr "{hline 13}{c +}{hline 64}"
}
else di in smcl in gr "{hline 13}{c BT}{hline 64}"
end
program define ELine
args name coef se alpha zero
local z = (`coef'-`zero')/`se'
local p = 2*norm(-abs(`z'))
local lb = `coef' - `alpha'*`se'
local ub = `coef' + `alpha'*`se'
di in smcl in gr %12s "`name'" " {c |}" /*
*/ in ye %11.6f `coef' /*
*/ " " %10.6f `se' /*
*/ %9.2f `z' /*
*/ " " %6.3f `p' /*
*/ " " %10.6f `lb' /*
*/ " " %10.6f `ub'
end
program define ILine
args name coef se alpha
local lb = `coef' - `alpha'*`se'
local ub = `coef' + `alpha'*`se'
di in smcl in gr %12s "`name'" " {c |}" /*
*/ in ye %11.6f `coef' /*
*/ " " %10.6f `se' /*
*/ _col(57) /*
*/ %10.6f `lb' /*
*/ " " %10.6f `ub'
end
prog def DISTab2
args level
local alpha=invnorm(.5+`level'/200)
local areau=`e(area)'+ `alpha'*`e(se_area)'
local areal=`e(area)'- `alpha'*`e(se_area)'
local deltau=`e(deltam)'+ `alpha'*`e(se_delm)'
local deltal=`e(deltam)'- `alpha'*`e(se_delm)'
local deu=`e(de)'+ `alpha'*`e(se_de)'
local del=`e(de)'- `alpha'*`e(se_de)'
local dau=`e(da)'+ `alpha'*`e(se_da)'
local dal=`e(da)'- `alpha'*`e(se_da)'
local cil `=string(`level')'
local cil `=length("`cil'")'
if `cil' == 2 {
local spaces " "
}
else if `cil' == 4 {
local spaces " "
}
else {
local spaces " "
}
di in smcl in gr _col(14) /*
*/ "{c |} Indices from binormal fit"
di in smcl in gr " Index {c |} Estimate Std. Err." /*
*/ `" `spaces'[`=strsubdp("`level'")'% Conf. Interval]"'
DivLine mid
ILine "ROC area" `e(area)' `e(se_area)' `alpha'
ILine "delta(m)" `e(deltam)' `e(se_delm)' `alpha'
ILine "d(e)" `e(de)' `e(se_de)' `alpha'
ILine "d(a)" `e(da)' `e(se_da)' `alpha'
DivLine bot
di in gr "(*) z test for slope==1"
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -