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.

298 lines
9.6 KiB
Plaintext

*! version 1.6.0 3/29/01
capture program drop brant
program define brant, rclass
version 6
tempvar touse
tempname bout d pvals ivchi ivout step1 step2 ologit
tempname XpWmmX iXpWmmX XpWmlX XpWllX iXpWllX DB DBp iDvBDp
syntax [, detail]
if "`e(cmd)'"!="ologit" {
di in r "brant can only be used after ologit"
exit
}
*to make output stata 6 or stata 7 compatible
cap version 7
if _rc!=0 {
local vers7 "no"
local smcl ""
local dash "-"
local vline "|"
local plussgn "+"
local topt "-"
local bottomt "-"
}
else { local vers7 "yes"
local smcl "in smcl "
local dash "{c -}"
local vline "{c |}"
local plussgn "{c +}"
local topt "{c TT}"
local bottomt "{c BT}"
}
version 6.0
local ocmd "`e(cmd)'"
if "`ocmd'"=="ologit" { local bcmd "logit" }
local depvar "`e(depvar)'"
gen `touse' = e(sample)
local wtis ""
if "`e(wtype)'"!="" {
di in r "-brant- does not work with ologit models with weights"
error 999
}
_perhs
local rhsnms "`r(rhsnms)'"
local nrhs "`r(nrhs)'"
_pecats
local numcats = `r(numcats)'
local catvals "`r(catvals)'"
local catnms "`r(catnms)'"
local catnms8 "`r(catnms8)'"
estimates hold `ologit'
*** estimate series of binary logits
local i = 1
while `i' <= `numcats'-1 {
local splitat : word `i' of `catvals'
tempvar dummy
quietly gen `dummy' = 0 if `depvar' <= `splitat' & `touse'==1
quietly replace `dummy' = 1 if `depvar' > `splitat' & `touse'==1
quietly `bcmd' `dummy' `rhsnms' `wtis' if `touse' == 1
_perhs
local binnrhs = "`r(nrhs)'"
if `nrhs' != `binnrhs' {
di in r "not all independent variables can be retained in all binary logits"
di in r "brant test cannot be computed"
exit 999
}
tempvar prob`i'
quietly predict `prob`i''
tempname b`i' V`i' bc`i'
mat `b`i'' = e(b)
mat `b`i'' = `b`i''[1, 1..`nrhs']
mat `V`i'' = e(V)
mat `V`i'' = `V`i''[1..`nrhs', 1..`nrhs']
mat `bc`i'' = e(b) /* with constant--for detail output only */
mat `bc`i'' = `bc`i'''
local outname "y>`splitat'"
local outname = substr("`outname'", 1, 8)
mat colnames `bc`i'' = "`outname'"
mat `bout' = nullmat(`bout'), `bc`i''
local i = `i' + 1
}
*** make variables for W(ml) matrices
local i = 1
while `i' <= `numcats'-1 {
local i2 = `i'
while `i2' <= `numcats'- 1 {
tempvar w`i'_`i2'
quietly gen `w`i'_`i2'' = `prob`i2'' - (`prob`i''*`prob`i2'')
local i2 = `i2' + 1
}
local i = `i' + 1
}
*** calculate variance Bm, Bl
local i = 1
while `i' <= `numcats'-1 {
local i2 = `i'
while `i2' <= `numcats'- 1 {
quietly {
* inverse(X'W(mm)X)
matrix accum `XpWmmX' = `rhsnms' [iw=`w`i'_`i''] if `touse'==1
matrix `iXpWmmX' = inv(`XpWmmX')
* X'W(ml)X
matrix accum `XpWmlX' = `rhsnms' [iw=`w`i'_`i2''] if `touse'==1
* inverse(X'W(ll)X)
matrix accum `XpWllX' = `rhsnms' [iw=`w`i2'_`i2''] if `touse'==1
matrix `iXpWllX' = inv(`XpWllX')
* product of three matrices
matrix `step1' = `iXpWmmX' * `XpWmlX'
tempname vb`i'_`i2'
matrix `vb`i'_`i2'' = `step1' * `iXpWllX'
}
mat `vb`i'_`i2''= `vb`i'_`i2''[1..`nrhs',1..`nrhs']
local i2 = `i2' + 1
}
local i = `i' + 1
}
* define var(B) matrix
local i = 1
while `i' <= `numcats'-1 {
tempname row`i'
local i2 = 1
while `i2' <= `numcats'- 1 {
quietly {
if `i'==`i2' { mat `row`i'' = nullmat(`row`i''), `V`i'' }
if `i'<`i2' { mat `row`i'' = nullmat(`row`i'') , `vb`i'_`i2'' }
if `i'>`i2' { mat `row`i'' = nullmat(`row`i'') , `vb`i2'_`i''' }
}
local i2 = `i2' + 1
}
local i = `i' + 1
}
* combine matrices
tempname varb
local i = 1
while `i' <= `numcats'-1 {
mat `varb' = nullmat(`varb') \ `row`i''
local i = `i' + 1
}
* make beta vector
tempname bstar
local i = 1
while `i' <= `numcats'-1 {
mat `bstar' = nullmat(`bstar') , `b`i''
local i = `i' + 1
}
mat `bstar' = `bstar''
* create design matrix for wald test; make I, -I, and 0 matrices
tempname id negid zero
local dim = `nrhs'
mat `id' = I(`dim')
mat rownames `id' = `rhsnms'
mat colnames `id' = `rhsnms'
mat `negid' = -1*`id'
mat rownames `negid' = `rhsnms'
mat colnames `negid' = `rhsnms'
mat `zero' = J(`dim', `dim', 0)
mat rownames `zero' = `rhsnms'
mat colnames `zero' = `rhsnms'
* dummy matrix
local i = 1
while `i' <= `numcats'-2 {
tempname drow`i'
local i2 = 1
while `i2' <= `numcats'- 1 {
quietly {
tempname feed
if `i2'==1 { mat `feed' = `id' }
else if `i2'-`i'==1 { mat `feed' = `negid' }
else { mat `feed' = `zero' }
mat `drow`i'' = nullmat(`drow`i'') , `feed'
}
local i2 = `i2' + 1
}
local i = `i' + 1
}
* combine matrices
local i = 1
while `i' <= `numcats'-2 {
mat `d' = nullmat(`d') \ `drow`i''
local i = `i' + 1
}
* terms of wald test
mat `DB' = `d' * `bstar'
mat `DBp' = `DB''
mat `step1' = `d'*`varb'
mat `step2' = `step1' * (`d'')
mat `iDvBDp' = inv(`step2')
*** calculate wald stat
tempname step1 wald waldout pout dfout
mat `step1' = `DBp' * `iDvBDp'
mat `wald' = `step1' * `DB'
sca `waldout' = `wald'[1,1]
sca `dfout' = `nrhs'*(`numcats'-2)
sca `pout' = chiprob(`dfout', `waldout')
tempname dtemp vbtemp bstemp
local i = 1
while `i' <= `nrhs' {
tempname d`i' vb`i' bstar`i'
local i2 = 1
while `i2' <= `numcats'-1 {
local row = ((`nrhs')*(`i2'-1)) + (`i')
tempname drow vbrow
local i3 = 1
while `i3' <= `numcats'-1 {
local column = ((`nrhs')*(`i3'-1)) + (`i')
if (`i2'<`numcats'-1) {
mat `dtemp' = `d'[`row',`column']
mat `drow' = nullmat(`drow') , `dtemp'
}
mat `vbtemp' = `varb'[`row',`column']
mat `vbrow' = nullmat(`vbrow') , `vbtemp'
local i3 = `i3' + 1
}
if (`i2'<`numcats'-1) { mat `d`i'' = nullmat(`d`i'') \ `drow' }
mat `vb`i'' = nullmat(`vb`i'') \ `vbrow'
mat `bstemp' = `bstar'[`row', 1]
mat `bstar`i'' = nullmat(`bstar`i'') \ `bstemp'
local i2 = `i2' + 1
}
local i = `i' + 1
}
*** wald test for each independent variable
tempname waldiv
local i = 1
while `i' <= `nrhs' {
tempname DB DBp iDvBDp step1 step2
mat `DB' = `d`i'' * `bstar`i''
mat `DBp' = `DB''
mat `step1' = `d`i''*`vb`i''
mat `step2' = `step1' * (`d`i''')
mat `iDvBDp' = inv(`step2')
tempname step1 wald`i'
mat `step1' = `DBp' * `iDvBDp'
mat `wald`i'' = `step1' * `DB'
mat `waldiv' = nullmat(`waldiv') \ `wald`i''
local i = `i' + 1
}
if "`detail'"!="" {
di _n in gr "Estimated coefficients from j-1 binary regressions"
mat list `bout', noheader
}
di _n in g "Brant Test of Parallel Regression Assumption"
di _n `smcl' in g " Variable `vline' chi2 p>chi2 df"
di `smcl' _dup(13) in g "`dash'" "`plussgn'" _dup(26) in g "`dash'"
di `smcl' in g " All `vline'" in y /*
*/ %10.2f `waldout' %9.3f `pout' %6.0f `dfout'
di `smcl' _dup(13) in g "`dash'" "`plussgn'" _dup(26) in g "`dash'"
* calculate p for individual wald tests
mat `pvals' = J(`nrhs', 1, 0)
local i = 1
local df = `numcats'-2
while `i' <= `nrhs' {
sca `ivchi' = `waldiv'[`i',1]
if `ivchi' >= 0 {
mat `pvals'[`i',1] = chiprob(`df',`ivchi')
}
if `ivchi' < 0 {
mat `pvals'[`i',1] = -999
}
local vnm : word `i' of `rhsnms'
*added for stata 7 compatibility
local printnm "`vnm'"
if "`vers7'"=="yes" { local printnm = abbrev("`printnm'", 12) }
di `smcl' in g %12s "`printnm'" _col(14) "`vline'" in y /*
*/ %10.2f `ivchi' %9.3f `pvals'[`i',1] %6.0f `df'
local i = `i' + 1
}
di `smcl' _dup(13) in g "`dash'" "`bottomt'" _dup(26) in g "`dash'"
di _n in g /*
*/ "A significant test statistic provides evidence that the parallel"
di in g "regression assumption has been violated."
mat `ivout' = `waldiv', `pvals'
mat rownames `ivout' = `rhsnms'
mat colnames `ivout' = chi2 p>chi2
estimates unhold `ologit'
return scalar chi2 = `waldout'
return scalar p = `pout'
return scalar df = `dfout'
return matrix ivtests `ivout'
end