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.
1091 lines
34 KiB
Plaintext
1091 lines
34 KiB
Plaintext
*! v.2.0.2 Confirmatory factor analysis, by Stas Kolenikov, skolenik at gmail dot com, 08 Sep 2009
|
|
program define confa, eclass properties( svyr svyb svyj ) sortpreserve
|
|
version 10.0
|
|
|
|
if replay() {
|
|
if ("`e(cmd)'" != "confa") error 301
|
|
Replay `0'
|
|
}
|
|
|
|
else {
|
|
Estimate `0'
|
|
}
|
|
|
|
end
|
|
|
|
program define Estimate, eclass properties( svyr svyb svyj )
|
|
|
|
syntax anything [if] [in] [aw pw iw/ ], [ ///
|
|
UNITvar(str) /// provides the list of factors where unit variance identification is used
|
|
FREE /// estimate all parameters as free -- the user provides identification through constraints
|
|
CONSTRaint(numlist) /// see the previous one
|
|
FROM(str) /// starting values, compliant with -ml init- syntax
|
|
LEVel(int $S_level) ROBust VCE(string) CLUster(passthru) /// standard errors and inference
|
|
LOGLEVel(int 1) /// logging level
|
|
CORRelated(string) /// correlated measurement errors
|
|
SUBtractone /// subtract one from the sample size in places
|
|
USENames ///
|
|
MISSing /// allow for special treatment of missing data
|
|
SVY * ]
|
|
|
|
* preliminary work
|
|
global CONFA_loglevel = `loglevel'
|
|
|
|
cap bothlist a \ b, global( CONFA_t )
|
|
if _rc==199 {
|
|
* listutil not installed
|
|
di as err "listutil not found, trying to install from SSC..."
|
|
ssc install listutil
|
|
}
|
|
|
|
if "`subtractone'"~="" local subtractone -1
|
|
|
|
*** MISSING
|
|
tempvar touse
|
|
marksample touse, zeroweight
|
|
global CONFA_touse `touse'
|
|
if $CONFA_loglevel > 2 tab $CONFA_touse
|
|
|
|
* weights?
|
|
if "`weight'" ~= "" {
|
|
global CONFA_wgt [`weight'=`exp']
|
|
global CONFA_wgti [iw=`exp']
|
|
}
|
|
|
|
* initial values?
|
|
if "`from'"~="" {
|
|
if "`from'" == "iv" | "`from'" == "IV" | "`from'" == "ivreg" | "`from'" == "2SLS" {
|
|
if "`unitvar'" ~= "" {
|
|
di as err "cannot specify from(`from') and unitvar at the same time"
|
|
CleanIt
|
|
exit 198
|
|
}
|
|
else global CONFA_init IV
|
|
}
|
|
else if "`from'" == "ones" global CONFA_init ones
|
|
else if "`from'" == "smart" global CONFA_init smart
|
|
else {
|
|
gettoken isitmat isitnot : from , parse(",")
|
|
cap confirm matrix `isitmat'
|
|
if _rc {
|
|
di "{err}Warning: matrix `isitmat' not found"
|
|
* do something sensible instead
|
|
global CONFA_init smart
|
|
}
|
|
else {
|
|
* do nothing --- let the hell break loose
|
|
global CONFA_init
|
|
}
|
|
}
|
|
}
|
|
|
|
* vce?
|
|
if "`vce'" ~= "" {
|
|
gettoken vce1 rest : vce, parse(" ,")
|
|
CheckVCE , `vce1'
|
|
local lvce = length("`vce'")
|
|
if `"`vce'"' != substr("sbentler", 1, max(2, `lvce')) ///
|
|
& `"`vce'"' != substr("satorrabentler", 1, max(3, `lvce')) {
|
|
local vceopt vce(`vce')
|
|
}
|
|
}
|
|
|
|
if $CONFA_loglevel > 2 di as text "Parsing..."
|
|
cap noi Parse `anything'
|
|
if _rc {
|
|
CleanIt
|
|
exit 198
|
|
}
|
|
|
|
* copy everything down -- -ivreg- cleans -sreturn-
|
|
local obsvar=s(obsvar)
|
|
global CONFA_obsvar `obsvar'
|
|
|
|
local nobsvar : word count `obsvar'
|
|
local nfactors = s(nfactors)
|
|
forvalues k=1/`nfactors' {
|
|
local indicators`k' = s(indicators`k')
|
|
local name`k' = s(name`k')
|
|
local factorlist `factorlist' `name`k''
|
|
}
|
|
|
|
* begin collecting the equations, starting values, bounds, and model structure
|
|
if $CONFA_loglevel > 2 di as text "Setting the structure up..."
|
|
Structure, unitvar(`unitvar') correlated(`correlated') `usenames'
|
|
local nicecorr $CONFA_t
|
|
* produces a bunch of globals
|
|
|
|
* did we need all this starting values business at all?
|
|
gettoken isitmat isitnot : from , parse(",")
|
|
cap confirm matrix `isitmat'
|
|
if ~_rc | strpos("`from'",",") {
|
|
* the user has provided the starting values
|
|
global CONFA_start `from'
|
|
}
|
|
else {
|
|
* not a matrix, no comma: use our ugly computations
|
|
global CONFA_start $CONFA_start, copy
|
|
}
|
|
* if "$CONFA_bounds" ~= "" global CONFA_bounds bounds($CONFA_bounds)
|
|
|
|
if $CONFA_loglevel > 3 di `" ml model lf confa_lf $CONFA_toML $CONFA_wgt, constraint($CONFA_constr `constraint') `svy' `robust' `cluster' init($CONFA_start) bounds($CONFA_bounds) `options' maximize"'
|
|
|
|
tempvar misspat touse1
|
|
global CONFA_miss `misspat'
|
|
qui gen byte `touse1' = 1-$CONFA_touse
|
|
if "`missing'" != "" {
|
|
if $CONFA_loglevel > 1 di "{txt}Working on missing values..."
|
|
* cycle over the observed variables, create missing indicators
|
|
forvalues k=1/`nobsvar' {
|
|
local thisvar : word `k' of `obsvar'
|
|
tempvar miss`k'
|
|
qui gen byte `miss`k'' = mi( `thisvar' ) if $CONFA_touse
|
|
local misslist `misslist' `miss`k''
|
|
}
|
|
|
|
* sort by pattern: relevant observations first
|
|
* when $CONFA_touse==0 `misslist' will be missing
|
|
qui {
|
|
bysort `touse1' `misslist' : gen long $CONFA_miss = (_n==1)
|
|
replace $CONFA_miss = sum( $CONFA_miss )
|
|
replace $CONFA_miss = . if mi( $CONFA_touse )
|
|
}
|
|
cap assert $CONFA_miss == 1 if $CONFA_touse
|
|
local anymissing = _rc
|
|
if !`anymissing' {
|
|
di "{txt}Option missing specified, but no missing data found"
|
|
}
|
|
else {
|
|
qui tab $CONFA_miss
|
|
di _n "{txt}Note: {res}" r(r) "{txt} patterns of missing data found"
|
|
}
|
|
if $CONFA_loglevel > 3 li $CONFA_miss `misslist'
|
|
}
|
|
* if "`anymissing'"=="0" | "`missing'" == "" {
|
|
else {
|
|
* -missing- option is omitted
|
|
if $CONFA_loglevel > 1 di "{err}NOT {txt}working on missing values"
|
|
markout $CONFA_touse `obsvar'
|
|
qui gen byte $CONFA_miss = 1 if $CONFA_touse
|
|
if $CONFA_loglevel > 2 {
|
|
sum `obsvar' if $CONFA_touse
|
|
tab $CONFA_miss, missing
|
|
}
|
|
}
|
|
|
|
|
|
cap noi ml model lf confa_lfm $CONFA_toML $CONFA_wgt if $CONFA_touse, ///
|
|
constraint($CONFA_constr `constraint') `svy' `robust' `cluster' `vceopt' ///
|
|
init($CONFA_start) bounds($CONFA_bounds) `options' `missing' ///
|
|
maximize
|
|
|
|
local mlrc = _rc
|
|
if `mlrc' {
|
|
CleanIt
|
|
error `mlrc'
|
|
}
|
|
|
|
* parametric matrices
|
|
tempname bb
|
|
mat `bb' = e(b)
|
|
global CONFA_loglevel -1
|
|
* to indicate to CONFA_StrucToSigma() that the matrices should be posted to Stata
|
|
qui mata : CONFA_StrucToSigma(st_matrix("`bb'"))
|
|
global CONFA_loglevel `loglevel'
|
|
* now, post all those matrices to ereturn
|
|
|
|
mat rownames CONFA_Sigma = `obsvar'
|
|
mat colnames CONFA_Sigma = `obsvar'
|
|
mat rownames CONFA_Lambda = `obsvar'
|
|
mat colnames CONFA_Lambda = `factorlist'
|
|
mat rownames CONFA_Phi = `factorlist'
|
|
mat colnames CONFA_Phi = `factorlist'
|
|
mat colnames CONFA_Theta = `obsvar'
|
|
mat rownames CONFA_Theta = `obsvar'
|
|
ereturn matrix Sigma = CONFA_Sigma, copy
|
|
ereturn matrix Lambda = CONFA_Lambda, copy
|
|
ereturn matrix Phi = CONFA_Phi, copy
|
|
ereturn matrix Theta = CONFA_Theta, copy
|
|
|
|
if "`missing'"!= "" ereturn local missing missing
|
|
|
|
eret local observed `obsvar'
|
|
eret local factors `factorlist'
|
|
if "`unitvar'" ~= "" {
|
|
* need to unwrap the contents of `unitvar'...
|
|
* or change its defintion from passthru to string
|
|
if "`unitvar'" == "_all" eret local unitvar `factorlist'
|
|
else eret local unitvar `unitvar'
|
|
}
|
|
forvalues k=1/`nfactors' {
|
|
eret local factor`k' `name`k'' : `indicators`k''
|
|
}
|
|
if "`correlated'"!="" eret local correlated `nicecorr'
|
|
|
|
if "`svy'`cluster'`exp'`robust'" == "" & "`vce'"!="robust" & substr("`vce'",1,2)!="cl" & "`missing'"=="" {
|
|
|
|
* if the data are not i.i.d., LRT is not applicable
|
|
* don't know what to do with missing data
|
|
|
|
tempname S Sindep trind
|
|
qui mat accum `S' = `obsvar' $CONFA_wgti if $CONFA_touse, dev nocons
|
|
mat `S' = `S' / ( e(N) `subtractone' )
|
|
mat `Sindep' = diag(vecdiag(`S'))
|
|
|
|
* degrees of freedom
|
|
local nconstr = `: word count $CONFA_constr' + `: word count `constraint''
|
|
local pstar = `nobsvar' * (`nobsvar' + 1) / 2
|
|
local df_m = rowsof(CONFA_Struc) - `nobsvar' - `nconstr'
|
|
ereturn scalar pstar = `pstar'
|
|
|
|
* test against independence
|
|
mat `trind' = trace( syminv(`Sindep') * `S' )
|
|
local trind = `trind'[1,1]
|
|
ereturn scalar ll_indep = -0.5 * `nobsvar' * e(N) * ln(2*_pi) - 0.5 * e(N) * ln(det(`Sindep')) - 0.5 * e(N) * `trind'
|
|
ereturn scalar lr_indep = 2*(e(ll)-e(ll_indep))
|
|
ereturn scalar df_indep = `pstar' - `nobsvar'
|
|
ereturn scalar p_indep = chi2tail(e(df_indep),e(lr_indep))
|
|
|
|
* goodness of fit test
|
|
ereturn scalar ll_0 = -0.5 * `nobsvar' * e(N) * ln(2*_pi) - 0.5 * e(N) * ln(det(`S')) - 0.5 * `nobsvar' * e(N)
|
|
ereturn scalar df_u = `pstar' - `df_m'
|
|
ereturn scalar lr_u = cond(e(df_u)==0,0,-2*(e(ll)-e(ll_0)))
|
|
ereturn scalar p_u = chi2tail(e(df_u),e(lr_u))
|
|
|
|
* make the g.o.f. test the default test
|
|
ereturn scalar df_m = `df_m'
|
|
ereturn local chi2type LR
|
|
ereturn scalar chi2 = e(lr_u)
|
|
ereturn scalar p = e(p_u)
|
|
|
|
* other crap
|
|
ereturn matrix S = `S'
|
|
|
|
if `"`vce'"'==substr("satorrabentler",1,max(3, length("`vce'"))) ///
|
|
| "`vce'" ==substr("sbentler",1,max(4, length("`vce'"))) {
|
|
* repost Satorra-Bentler covariance matrix
|
|
* not defined for complex survey data,
|
|
|
|
cap noi SatorraBentler, constraint(`constraint') `missing'
|
|
|
|
if _rc {
|
|
di as err "Satorra-Bentler standard errors are not supported; revert to vce(oim)"
|
|
}
|
|
else {
|
|
|
|
tempname SBVar SBV Delta Gamma VV U trUG2 Tdf
|
|
mat `SBVar' = r(SBVar)
|
|
mat `Delta' = r(Delta)
|
|
mat `Gamma' = r(Gamma)
|
|
mat `SBV' = r(SBV)
|
|
mat `VV' = e(V)
|
|
mat `SBVar' = ( `VV'[1..`nobsvar',1..`nobsvar'], `VV'[1..`nobsvar',`nobsvar'+1 ...] ///
|
|
\ `VV'[`nobsvar'+1...,1..`nobsvar'], `SBVar'[`nobsvar'+1...,`nobsvar'+1...] )
|
|
ereturn repost V = `SBVar'
|
|
ereturn matrix SBGamma = `Gamma', copy
|
|
ereturn matrix SBDelta = `Delta', copy
|
|
ereturn matrix SBV = `SBV', copy
|
|
ereturn local vce SatorraBentler
|
|
ereturn local vcetype "Satorra-Bentler"
|
|
|
|
* compute the corrected tests, too
|
|
* only takes care of the covariance structure
|
|
* Satorra-Bentler 1994
|
|
mat `U' = `SBV' - `SBV'*`Delta'*syminv(`Delta''*`SBV'*`Delta')*`Delta''*`SBV'
|
|
ereturn matrix SBU = `U'
|
|
mat `U' = trace( e(SBU)*e(SBGamma) )
|
|
ereturn scalar SBc = `U'[1,1]/e(df_u)
|
|
ereturn scalar Tsc = e(lr_u)/e(SBc) * (e(N) `subtractone' ) / e(N)
|
|
ereturn scalar p_Tsc = chi2tail( e(df_u), e(Tsc) )
|
|
|
|
mat `trUG2' = trace( e(SBU)*`Gamma'*e(SBU)*`Gamma')
|
|
ereturn scalar SBd = `U'[1,1]*`U'[1,1]/`trUG2'[1,1]
|
|
ereturn scalar Tadj = ( e(SBd)/`U'[1,1]) * e(lr_u) * (e(N) `subtractone' ) / e(N)
|
|
ereturn scalar p_Tadj = chi2tail( e(SBd), e(Tadj) )
|
|
|
|
* saddlepoint approximation comes here!!!
|
|
|
|
* Yuan-Bentler 1997
|
|
ereturn scalar T2 = e(lr_u)/(1+e(lr_u)/e(N) )
|
|
ereturn scalar p_T2 = chi2tail( e(df_u), e(T2) )
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
* are we done yet?
|
|
ereturn matrix CONFA_Struc = CONFA_Struc
|
|
ereturn local predict confa_p
|
|
ereturn local estat_cmd confa_estat
|
|
ereturn local cmd confa
|
|
|
|
Replay
|
|
|
|
CleanIt
|
|
|
|
end
|
|
|
|
program define CleanIt
|
|
* just in case
|
|
return clear
|
|
|
|
* release the constraints
|
|
constr drop $CONFA_constr
|
|
|
|
* clear the globals
|
|
if $CONFA_loglevel < 3 {
|
|
global CONFA_constr
|
|
global CONFA_init
|
|
global CONFA_loglevel
|
|
global CONFA_toML
|
|
global CONFA_start
|
|
global CONFA_bounds
|
|
global CONFA_args
|
|
global CONFA_constr
|
|
global CONFA_obsvar
|
|
global CONFA_t
|
|
global CONFA_wgt
|
|
global CONFA_wgti
|
|
}
|
|
end
|
|
|
|
program define Parse, sclass
|
|
|
|
* number of factors?
|
|
local input `0'
|
|
|
|
mata: st_local("nfactors",strofreal(CONFA_NF(`"`input'"')))
|
|
|
|
if `nfactors' == 0 {
|
|
* something terrible happened
|
|
di as err "incorrect factor specification"
|
|
exit 198
|
|
}
|
|
|
|
sreturn local nfactors = `nfactors'
|
|
|
|
tokenize `input', parse("()")
|
|
|
|
local k = 0
|
|
while "`1'"~="" {
|
|
* right now, `1' should contain an opening bracket
|
|
if "`1'"~="(" {
|
|
* the first character is not a "("
|
|
di as err "incorrect factor specification"
|
|
exit 198
|
|
}
|
|
else {
|
|
* the first character IS a "("
|
|
mac shift
|
|
* right now, `1' should contain a factor-type statement
|
|
local ++k
|
|
local factor`k' `1'
|
|
mac shift
|
|
* right now, `1' should contain a closing bracket
|
|
if "`1'"~=")" {
|
|
* the first character is not a ")"
|
|
di as err "incorrect factor specification"
|
|
exit 198
|
|
}
|
|
else mac shift
|
|
* it may contain a space, I guess
|
|
* then -mac shift- it again
|
|
if trim("`1'")=="" mac shift
|
|
}
|
|
}
|
|
forvalues k=1/`nfactors' {
|
|
* now, parse each factor statement
|
|
tokenize `factor`k'', parse(":")
|
|
sreturn local name`k' `1'
|
|
* `2' is the colon
|
|
unab indicators : `3'
|
|
sreturn local indicators`k' `indicators'
|
|
local obsvar `obsvar' `indicators'
|
|
}
|
|
|
|
cap uniqlist `obsvar'
|
|
if _rc == 199 {
|
|
* uniqlist not found
|
|
di as err "uniqlist not found, trying to install from SSC..."
|
|
ssc install listutil
|
|
uniqlist `obsvar'
|
|
}
|
|
local obsvar = r(list)
|
|
|
|
* mata: st_local("obsvar",CONFA_UL(`"`obsvar'"'))
|
|
|
|
sreturn local obsvar `obsvar'
|
|
sreturn local nobsvar = `: word count `obsvar''
|
|
|
|
end
|
|
|
|
program define Structure
|
|
|
|
syntax , [unitvar(str) correlated(str) usenames]
|
|
|
|
* implement usenames:
|
|
* the parameters go along with the factor and variable names
|
|
* rather than matrix indices
|
|
|
|
* utilize all the sreturn results
|
|
if $CONFA_loglevel > 3 sreturn list
|
|
|
|
* copy everything down -- -ivreg- cleans -sreturn-
|
|
local obsvar=s(obsvar)
|
|
local nobsvar : word count `obsvar'
|
|
local nfactors = s(nfactors)
|
|
|
|
forvalues k=1/`nfactors' {
|
|
local indicators`k' = s(indicators`k')
|
|
local name`k' = s(name`k')
|
|
local factorlist `factorlist' `name`k''
|
|
}
|
|
|
|
if "`unitvar'" == "_all" {
|
|
local unitvar `factorlist'
|
|
}
|
|
|
|
* set up the labeling system
|
|
if "`usenames'" != "" {
|
|
* give the parameters varname labels
|
|
forvalues k=1/`nobsvar' {
|
|
local o`k' : word `k' of `obsvar'
|
|
}
|
|
forvalues k=1/`nfactors' {
|
|
local f`k' `name`k''
|
|
}
|
|
}
|
|
else {
|
|
* give the parameters numberic lables
|
|
forvalues k=1/`nobsvar' {
|
|
local o`k' `k'
|
|
}
|
|
forvalues k=1/`nfactors' {
|
|
local f`k' `k'
|
|
}
|
|
}
|
|
|
|
* returns:
|
|
* - ML equations
|
|
* - ML bounds
|
|
* - the structure matrix
|
|
* - ML statement for the likelihood evaluator
|
|
|
|
* initialize everything
|
|
local eqno = 0
|
|
global CONFA_toML
|
|
global CONFA_start
|
|
global CONFA_args
|
|
global CONFA_constr
|
|
global CONFA_bounds
|
|
mata : CONFA_Struc = J(0,4,.)
|
|
|
|
* process the means first
|
|
tokenize `obsvar'
|
|
forvalues j=1/`nobsvar' {
|
|
* 1. equations to ML
|
|
local ++eqno
|
|
global CONFA_toML $CONFA_toML (mean_`o`j'':)
|
|
* 2. starting values
|
|
sum ``j'', mean
|
|
global CONFA_start $CONFA_start `r(mean)'
|
|
* 3. confa_lf arguments
|
|
global CONFA_args $CONFA_args mean_`o`j''
|
|
* 4. CONFA structure
|
|
mata : CONFA_Struc = CONFA_Struc \ (1, `eqno', `j', 0)
|
|
}
|
|
|
|
* next, process lambda's
|
|
forvalues k=1/`nfactors' {
|
|
* determine if unitvar is needed here
|
|
bothlist `name`k'' \ `unitvar', global(CONFA_t)
|
|
if "$CONFA_t" ~= "" {
|
|
* identification by unit variance, no scaling variables
|
|
local scalevar
|
|
}
|
|
else {
|
|
* identification by the scaling variable: the 1st one on the list
|
|
local scalevar : word 1 of `indicators`k''
|
|
}
|
|
forvalues j=1/`nobsvar' {
|
|
* determine whether `k'-th factor loads on `j'-th variable
|
|
if strpos( "`indicators`k''", "``j''") {
|
|
* 1. equations to ML
|
|
local ++eqno
|
|
global CONFA_toML $CONFA_toML (lambda_`o`j''_`f`k'':)
|
|
* 2. starting values
|
|
local r2_`j' = 0.5
|
|
if "``j''" == "`scalevar'" {
|
|
* the current one is the scaling variable
|
|
* set up the constraints, initialize to 1
|
|
global CONFA_start $CONFA_start 1
|
|
constraint free
|
|
local nconstr = r(free)
|
|
constraint `nconstr' [lambda_`o`j''_`f`k'']_cons = 1
|
|
global CONFA_constr $CONFA_constr `nconstr'
|
|
}
|
|
else if "$CONFA_init" == "IV" {
|
|
* initialize by a simple version of instrumental variables
|
|
* use the remaining indicators of this factor as instruments
|
|
dellist `indicators`k'', delete(`scalevar' ``j'')
|
|
local ivlist = r(list)
|
|
if "`ivlist'" == "." {
|
|
di as err "Warning: no instruments available for ``j''"
|
|
local ivl = 1
|
|
}
|
|
else {
|
|
qui ivreg ``j'' (`scalevar' = `ivlist')
|
|
local ivl = _b[`scalevar']
|
|
}
|
|
global CONFA_start $CONFA_start `ivl'
|
|
if !mi(e(r2)) local r2_`j' = e(r2)
|
|
}
|
|
else if "$CONFA_init" == "ones" {
|
|
global CONFA_start $CONFA_start 1
|
|
}
|
|
else {
|
|
* no init options
|
|
global CONFA_start $CONFA_start 0
|
|
}
|
|
global CONFA_bounds $CONFA_bounds /lambda_`o`j''_`f`k'' -100 100
|
|
* 3. confa_lf arguments
|
|
global CONFA_args $CONFA_args lambda_`o`j''_`f`k''
|
|
* 4. CONFA structure
|
|
mata : CONFA_Struc = CONFA_Struc \ (2, `eqno', `j', `k')
|
|
}
|
|
}
|
|
}
|
|
|
|
* next, process Phi matrix
|
|
forvalues k=1/`nfactors' {
|
|
local scalevar1 : word 1 of `indicators`k''
|
|
|
|
foreach kk of numlist `k'/1 {
|
|
* 1. equations to ML
|
|
local ++eqno
|
|
global CONFA_toML $CONFA_toML (phi_`f`kk''_`f`k'':)
|
|
* 2. starting values
|
|
if `k' == `kk' {
|
|
* diagonal entry
|
|
bothlist `name`k'' \ `unitvar', global(CONFA_t)
|
|
if "$CONFA_t" ~= "" {
|
|
* identification by unit variance
|
|
constraint free
|
|
local nconstr = r(free)
|
|
constraint `nconstr' [phi_`f`k''_`f`k'']_cons = 1
|
|
global CONFA_constr $CONFA_constr `nconstr'
|
|
local v`k' = 1
|
|
}
|
|
else {
|
|
* identification by the scaling variable
|
|
if "$CONFA_init" == "smart" | "$CONFA_init" == "IV" {
|
|
qui sum `scalevar1'
|
|
local v`k' = r(Var)*0.5
|
|
}
|
|
else if "$CONFA_init" == "ones" local v`k' = 1
|
|
else local v`k' = 0
|
|
}
|
|
global CONFA_start $CONFA_start `v`k''
|
|
global CONFA_bounds $CONFA_bounds /phi_`f`k''_`f`kk'' 0 1000
|
|
}
|
|
else {
|
|
* off-diagonal entry
|
|
if "$CONFA_init" == "smart" | "$CONFA_init" == "IV" {
|
|
local scalevar2 : word 1 of `indicators`kk''
|
|
qui corr `scalevar1' `scalevar2'
|
|
local v = 0.5*r(rho)*sqrt(`v`k''*`v`kk'')
|
|
}
|
|
else if "$CONFA_init" == "ones" local v = 0.5
|
|
else local v = 0
|
|
local vv = 1.5*abs(`v') + 0.01
|
|
global CONFA_start $CONFA_start `v'
|
|
global CONFA_bounds $CONFA_bounds /phi_`f`kk''_`f`k'' -`vv' `vv'
|
|
}
|
|
* 3. confa_lf arguments
|
|
global CONFA_args $CONFA_args phi_`f`kk''_`f`k''
|
|
* 4. CONFA structure
|
|
mata : CONFA_Struc = CONFA_Struc \ (3, `eqno', `kk', `k')
|
|
}
|
|
|
|
}
|
|
|
|
* residual variances
|
|
forvalues j=1/`nobsvar' {
|
|
* 1. equations to ML
|
|
local ++eqno
|
|
global CONFA_toML $CONFA_toML (theta_`o`j'':)
|
|
* 2. starting values
|
|
if "$CONFA_init" == "ones" {
|
|
local v_`j' = 1
|
|
}
|
|
else if "$CONFA_init" == "IV" | "$CONFA_init" == "smart" {
|
|
qui sum ``j''
|
|
local v_`j' = r(Var)*(1-`r2_`j'')
|
|
}
|
|
else local v_`j' = 0.01
|
|
global CONFA_start $CONFA_start `v_`j''
|
|
global CONFA_bounds $CONFA_bounds /theta_`o`j'' 0 1000
|
|
* 3. confa_lf arguments
|
|
global CONFA_args $CONFA_args theta_`o`j''
|
|
* 4. CONFA structure
|
|
mata : CONFA_Struc = CONFA_Struc \ (4, `eqno', `j', 0)
|
|
}
|
|
|
|
* the error correlations
|
|
while "`correlated'" != "" {
|
|
gettoken corrpair correlated : correlated , match(m)
|
|
gettoken corr1 corrpair : corrpair, parse(":")
|
|
unab corr1 : `corr1'
|
|
gettoken sc corr2 : corrpair, parse(":")
|
|
unab corr2 : `corr2'
|
|
|
|
* make sure both are present in the list of observed variables
|
|
|
|
poslist `obsvar' \ `corr1', global(CONFA_t)
|
|
local k1 = $CONFA_t
|
|
if `k1' == 0 {
|
|
di as err "`corr1' is not among the observed variables"
|
|
CleanIt
|
|
exit 198
|
|
}
|
|
poslist `obsvar' \ `corr2', global(CONFA_t)
|
|
local k2 = $CONFA_t
|
|
if `k2' == 0 {
|
|
di as err "`corr2' is not among the observed variables"
|
|
CleanIt
|
|
exit 198
|
|
}
|
|
|
|
* will be empty @ the first call
|
|
local nicecorr `nicecorr' (`corr1':`corr2')
|
|
|
|
* 1. equations to ML
|
|
local ++eqno
|
|
global CONFA_toML $CONFA_toML (theta_`o`k1''_`o`k2'':)
|
|
* 2. starting values
|
|
global CONFA_start $CONFA_start 0
|
|
local vv = sqrt(`v_`k1''*`v_`k2'')
|
|
global CONFA_bounds $CONFA_bounds /theta_`o`k1''_`o`k2'' -`vv' `vv'
|
|
* 3. confa_lf arguments
|
|
global CONFA_args $CONFA_args theta_`o`k1''_`o`k2''
|
|
* 4. CONFA structure
|
|
mata : CONFA_Struc = CONFA_Struc \ (5, `eqno', `k1', `k2')
|
|
}
|
|
if "`nicecorr'"!="" global CONFA_t `nicecorr'
|
|
|
|
if $CONFA_loglevel > 3 {
|
|
di as text "ML input (" as res `: word count $CONFA_toML' as text "): " as res "$CONFA_toML"
|
|
di as text "Starting values (" as res `: word count $CONFA_start' as text "): " as res "$CONFA_start"
|
|
di as text "Likelihood evaluator (" as res `: word count $CONFA_args' as text"): " as res "$CONFA_args"
|
|
di as text "Constraints (" as res `nfactors' as text "): " as res "$CONFA_constr"
|
|
di as text "Correlated errors: " as res "`nicecorr'"
|
|
constraint dir $CONFA_constr
|
|
mata : CONFA_Struc
|
|
}
|
|
mata : st_matrix("CONFA_Struc",CONFA_Struc)
|
|
|
|
|
|
end
|
|
|
|
program define Replay
|
|
|
|
syntax, [posvar llu(str) level(passthru)]
|
|
|
|
* get the implied matrix
|
|
tempname bb Sigma
|
|
mat `bb' = e(b)
|
|
* mata : st_matrix("Sigma",CONFA_StrucToSigma(st_matrix("`bb'")))
|
|
mat `Sigma' = e(Sigma)
|
|
mat CONFA_Struc = e(CONFA_Struc)
|
|
|
|
* determine what kind of labeling has been used
|
|
* RATHER FRAGILE: checking for mean_1 rather than trying to find
|
|
* whether option usenames was specified
|
|
cap local whatis = [mean_1]_cons
|
|
if _rc {
|
|
* mean_1 not found => labeling by names
|
|
forvalues k=1/`: word count `e(observed)' ' {
|
|
local o`k' : word `k' of `e(observed)'
|
|
}
|
|
forvalues k=1/`: word count `e(factors)' ' {
|
|
local f`k' : word `k' of `e(factors)'
|
|
}
|
|
}
|
|
else {
|
|
* mean_1 was found => labeling by numbers
|
|
forvalues k=1/`: word count `e(observed)' ' {
|
|
local o`k' `k'
|
|
}
|
|
forvalues k=1/`: word count `e(factors)' ' {
|
|
local f`k' `k'
|
|
}
|
|
}
|
|
|
|
* header
|
|
di _n as text "`e(crittype)' = " as res e(ll) _col(59) as text "Number of obs = " as res e(N)
|
|
di as text "{hline 13}{c TT}{hline 64}"
|
|
if "`e(vcetype)'" ~= "" {
|
|
di as text " {c |} {center 15:`e(vcetype)'}"
|
|
}
|
|
di as text " {c |} Coef. Std. Err. z P>|z| [$S_level% Conf. Interval]"
|
|
di as text "{hline 13}{c +}{hline 64}"
|
|
|
|
tokenize `e(observed)'
|
|
local nobsvar : word count `e(observed)'
|
|
|
|
* let's go equation by equation
|
|
local eqno = 0
|
|
|
|
* Means
|
|
_diparm __lab__, label("Means") eqlabel
|
|
forvalues j = 1/`nobsvar' {
|
|
local ++eqno
|
|
_diparm mean_`o`j'' , label("``j''") prob `level'
|
|
}
|
|
|
|
* Loadings
|
|
_diparm __lab__, label("Loadings") eqlabel
|
|
local ++eqno // to point to the next line
|
|
forvalues k=1/`: word count `e(factors)' ' {
|
|
_diparm __lab__ , label("`: word `k' of `e(factors)' '")
|
|
while CONFA_Struc[`eqno',1]<=2 & CONFA_Struc[`eqno',4]==`k' {
|
|
local j = CONFA_Struc[`eqno',3]
|
|
_diparm lambda_`o`j''_`f`k'', label("``j''") prob `level'
|
|
local ++eqno
|
|
}
|
|
}
|
|
|
|
* Factor covariance
|
|
_diparm __lab__, label("Factor cov.") eqlabel
|
|
forvalues k=1/`: word count `e(factors)'' {
|
|
foreach kk of numlist `k'/1 {
|
|
_diparm phi_`f`kk''_`f`k'', label("`: word `kk' of `e(factors)''-`: word `k' of `e(factors)''") prob `level'
|
|
local ++eqno
|
|
}
|
|
}
|
|
|
|
* Error variances
|
|
_diparm __lab__, label("Var[error]") eqlabel
|
|
forvalues j= 1/`nobsvar' {
|
|
_diparm theta_`o`j'' , label("``j''") prob `level'
|
|
local ++eqno
|
|
}
|
|
|
|
* Error correlations
|
|
if `eqno' <= rowsof(CONFA_Struc) & CONFA_Struc[`eqno',1] == 5 {
|
|
_diparm __lab__, label("Cov[error]") eqlabel
|
|
while (`eqno' <= rowsof(CONFA_Struc) & CONFA_Struc[`eqno',1] == 5) {
|
|
local k1 = CONFA_Struc[`eqno',3]
|
|
local k2 = CONFA_Struc[`eqno',4]
|
|
_diparm theta_`o`k1''_`o`k2'', label("``k1''-``k2''") prob `level'
|
|
* range check: what am I supposed to check here? Hm...
|
|
local ++eqno
|
|
}
|
|
}
|
|
|
|
if "`e(vcetype)'"~="Robust" & "`e(missing)'"=="" {
|
|
di as text "{hline 13}{c +}{hline 64}"
|
|
di as text "R2{col 14}{c |}"
|
|
forvalues j = 1/`nobsvar' {
|
|
qui sum ``j'' if e(sample)
|
|
local r2 = (`Sigma'[`j',`j']-_b[theta_`o`j'':_cons])/r(Var)
|
|
di as text %12s "``j''" "{col 14}{c |}{col 20}" as res %6.4f `r2'
|
|
}
|
|
}
|
|
|
|
|
|
di as text "{hline 13}{c BT}{hline 64}"
|
|
|
|
if e(df_u)>0 {
|
|
di as text _n "Goodness of fit test: LR = " as res %6.3f e(lr_u) ///
|
|
as text _col(40) "; Prob[chi2(" as res %2.0f e(df_u) as text ") > LR] = " as res %6.4f e(p_u)
|
|
}
|
|
else {
|
|
di as text "No degrees of freedom to perform the goodness of fit test"
|
|
}
|
|
di as text "Test vs independence: LR = " as res %6.3f e(lr_indep) ///
|
|
as text _col(40) "; Prob[chi2(" as res %2.0f e(df_indep) as text ") > LR] = " as res %6.4f e(p_indep)
|
|
|
|
if "`e(vce)'" == "SatorraBentler" & e(df_u)>0 {
|
|
* need to report all those corrected statistics
|
|
|
|
di as text _n "Satorra-Bentler Tsc" _col(26) "= " as res %6.3f e(Tsc) ///
|
|
as text _col(40) "; Prob[chi2(" as res %2.0f e(df_u) as text ") > Tsc ] = " as res %6.4f e(p_Tsc)
|
|
|
|
di as text "Satorra-Bentler Tadj" _col(26) "= " as res %6.3f e(Tadj) ///
|
|
as text _col(40) "; Prob[chi2(" as res %4.1f e(SBd) as text ") > Tadj] = " as res %6.4f e(p_Tadj)
|
|
|
|
di as text "Yuan-Bentler T2" _col(26) "= " as res %6.3f e(T2) ///
|
|
as text _col(40) "; Prob[chi2(" as res %2.0f e(df_u) as text ") > T2 ] = " as res %6.4f e(p_T2)
|
|
}
|
|
|
|
if "`e(vce)'" == "BollenStine" {
|
|
* need to report Bollen-Stine measures
|
|
di as text _n "Bollen-Stine simulated Prob[ LR > " as res %6.4f e(lr_u) as text " ] = " as res %6.4f e(p_u_BS) ///
|
|
as text _n "Based on " as res e(B_BS) as text " replications. " ///
|
|
as text "The bootstrap 90% interval: (" as res %6.3f e(T_BS_05) as text "," ///
|
|
as res %6.3f e(T_BS_95) as text ")"
|
|
|
|
|
|
}
|
|
|
|
mat drop CONFA_Struc
|
|
|
|
end
|
|
|
|
**************************** Satorra-Bentler covariance matrix code
|
|
|
|
program SatorraBentler, rclass
|
|
syntax [, noisily constraint(numlist) missing]
|
|
|
|
if "`missing'"!="" {
|
|
di "{err}cannot specify Satorra-Bentler standard errors with missing data"
|
|
exit 198
|
|
}
|
|
|
|
* assume the maximization completed, the results are in memory as -ereturn data-
|
|
* we shall just return the resulting matrix
|
|
|
|
* assume sample is restricted to e(sample)
|
|
* preserve
|
|
* keep if e(sample)
|
|
|
|
* get the variable names
|
|
tempname VV bb
|
|
mat `bb' = e(b)
|
|
mat `VV' = e(V)
|
|
local p : word count $CONFA_obsvar
|
|
qui count if $CONFA_touse
|
|
local NN = r(N)
|
|
|
|
* compute the implied covariance matrix
|
|
tempname Lambda Theta Phi Sigma
|
|
mata : st_matrix("`Sigma'",CONFA_StrucToSigma(st_matrix("`bb'")))
|
|
|
|
* compute the empirical cov matrix
|
|
tempname SampleCov
|
|
qui mat accum `SampleCov' = $CONFA_obsvar $CONFA_wgti if $CONFA_touse , nocons dev
|
|
* divide by sum of weights instead???
|
|
mat `SampleCov' = `SampleCov' / (`NN'-1)
|
|
|
|
* compute the matrix Gamma (fourth moments)
|
|
if $CONFA_loglevel > 4 {
|
|
di as text "Computing the Gamma matrix of fourth moments..."
|
|
}
|
|
tempname Gamma
|
|
SBGamma $CONFA_obsvar if $CONFA_touse
|
|
mat `Gamma' = r(Gamma)
|
|
return add
|
|
|
|
* compute the V matrix, the normal theory weight
|
|
if $CONFA_loglevel > 4 {
|
|
di as text "Computing the V matrix..."
|
|
}
|
|
SBV `SampleCov' `noisily'
|
|
if !mi(r(needmatsize)) {
|
|
di as err "matsize too small; need at least " r(needmatsize)
|
|
exit 908
|
|
}
|
|
tempname V
|
|
mat `V' = r(SBV)
|
|
return add
|
|
|
|
* compute the Delta matrix
|
|
if $CONFA_loglevel > 4 {
|
|
di as text "Computing the Delta matrix..."
|
|
}
|
|
|
|
tempname Delta DeltaId
|
|
noi mata : SBStrucToDelta("`Delta'")
|
|
|
|
*** put the pieces together now
|
|
|
|
* enact the constraints!
|
|
SBconstr `bb', constraint(`constraint')
|
|
|
|
* zero out the rows of Delta that correspond to fixed parameters
|
|
mat `DeltaId' = `Delta' * diag( r(Fixed) )
|
|
|
|
local dcnames : colfullnames `bb'
|
|
local drnames : rownames `Gamma'
|
|
mat colnames `DeltaId' = `dcnames'
|
|
mat rownames `DeltaId' = `drnames'
|
|
return matrix Delta = `DeltaId', copy
|
|
|
|
tempname VVV
|
|
mat `VVV' = ( `DeltaId'' * `V' * `DeltaId' )
|
|
mat `VVV' = syminv(`VVV')
|
|
mat `VVV' = `VVV' * ( `DeltaId'' * `V' * `Gamma' * `V' * `DeltaId' ) * `VVV'/`NN'
|
|
|
|
* add the covariance matrix for the means, which is just Sigma/_N
|
|
* weights!
|
|
* third moments!
|
|
return matrix SBVar = `VVV'
|
|
|
|
end
|
|
* of satorrabentler
|
|
|
|
* Compute Gamma: the fourth moments matrix -- check!
|
|
program define SBGamma, rclass
|
|
syntax varlist [if] [in]
|
|
unab varlist : `varlist'
|
|
tokenize `varlist'
|
|
|
|
marksample touse
|
|
|
|
local p: word count `varlist'
|
|
|
|
forvalues k=1/`p' {
|
|
* make up the deviations; weights are used in a weird way
|
|
*** MISSING: change r(mean) to _b[whatever] ?
|
|
qui sum ``k'' $CONFA_wgti if `touse', meanonly
|
|
tempvar d`k'
|
|
qui g double `d`k'' = ``k'' - r(mean) if `touse'
|
|
local dlist `dlist' `d`k''
|
|
}
|
|
|
|
local pstar = `p'*(`p'+1)/2
|
|
forvalues k=1/`pstar' {
|
|
tempvar b`k'
|
|
qui g double `b`k'' = .
|
|
local blist `blist' `b`k''
|
|
}
|
|
|
|
|
|
* convert into vech (z_i-bar z)(z_i-bar z)'
|
|
mata : SBvechZZtoB("`dlist'","`blist'")
|
|
|
|
* blist now should contain the moments around the sample means
|
|
* we need to get their covariance matrix
|
|
|
|
tempname Gamma
|
|
qui mat accum `Gamma' = `blist' $CONFA_wgti if `touse', dev nocons
|
|
mat `Gamma' = `Gamma'/(r(N)-1)
|
|
mata : Gamma = st_matrix( "`Gamma'" )
|
|
|
|
* make nice row and column names
|
|
forvalues i=1/`p' {
|
|
forvalues j=`i'/`p' {
|
|
local namelist `namelist' ``i''_X_``j''
|
|
}
|
|
}
|
|
mat colnames `Gamma' = `namelist'
|
|
mat rownames `Gamma' = `namelist'
|
|
|
|
return matrix Gamma = `Gamma'
|
|
|
|
end
|
|
* of computing Gamma
|
|
|
|
* compute V = 1/2 D' (Sigma \otimes Sigma) D
|
|
* normal theory weight matrix, see Satorra (1992), eq (24) -- check!
|
|
program define SBV, rclass
|
|
args A noisily
|
|
tempname D Ainv V
|
|
local p = rowsof(`A')
|
|
if $CONFA_loglevel > 3 di as text "Computing the duplication matrix..."
|
|
mata : Dupl(`p',"`D'")
|
|
mat `Ainv' = syminv(`A')
|
|
cap mat `V' = .5*`D''* (`Ainv' # `Ainv') * `D'
|
|
if _rc == 908 {
|
|
* need a larger matrix
|
|
return scalar needmatsize = rowsof(`A')*rowsof(`A')
|
|
}
|
|
else {
|
|
return matrix SBV = `V'
|
|
}
|
|
end
|
|
* of computing V
|
|
|
|
program define SBconstr, rclass
|
|
* need to figure out whether a constraint has the form [parameter]_cons = value,
|
|
* and to nullify the corresponding column
|
|
syntax anything, [constraint(numlist)]
|
|
local bb `anything'
|
|
* that's the name of the parameter vector, a copy of e(b)
|
|
tempname Iq
|
|
mat `Iq' = J(1,colsof(`bb'),1)
|
|
tokenize $CONFA_constr `constraint'
|
|
while "`1'" ~= "" {
|
|
constraint get `1'
|
|
local constr `r(contents)'
|
|
gettoken param value : constr, parse("=")
|
|
* is the RHS indeed a number?
|
|
local value = substr("`value'",2,.)
|
|
confirm number `value'
|
|
* parse the square brackets and turn them into colon
|
|
* replace the opening brackets with nothing, and closing brackets, with colon
|
|
* that way, we will get "parameter:_cons", which is the format of e(b) labels
|
|
local param = subinstr("`param'","["," ",1)
|
|
local param = subinstr("`param'","]",":",1)
|
|
local param = trim("`param'")
|
|
local coln = colnumb(`bb',"`param'" )
|
|
mat `Iq'[1,`coln']=0
|
|
|
|
mac shift
|
|
}
|
|
return matrix Fixed = `Iq'
|
|
end
|
|
|
|
program define CheckVCE
|
|
syntax [anything] , [ROBust CLuster oim opg SBentler SATorrabentler BOOTstrap JACKknife]
|
|
if "`bootstrap'" ~= "" {
|
|
di "{err}vce(bootstrap) not allowed, but you can run {inp}bootstrap ... : confa ... {err}instead."
|
|
CleanIt
|
|
exit 198
|
|
}
|
|
if "`jackknife'" ~= "" {
|
|
di "{err}vce(jackknife) not allowed, but you can run {inp}jackknife ... : confa ... {err}instead."
|
|
CleanIt
|
|
exit 198
|
|
}
|
|
end
|
|
|
|
exit
|
|
|
|
if "$SOCST" == "c:\-socialstat" {
|
|
// at home, run the Mata file
|
|
do C:\-Mizzou\CONFA\confa.mata
|
|
}
|
|
else {
|
|
// for public release, add Mata code
|
|
mata : mata mlib index
|
|
}
|
|
|
|
Globals used:
|
|
CONFA_init -- initialization type
|
|
CONFA_loglevel -- detail level
|
|
CONFA_toML -- model statement for -ml model-
|
|
CONFA_start -- default starting values
|
|
CONFA_bounds -- ml search bounds
|
|
CONFA_args -- the list of parameters, to appear in -confa_lf-
|
|
CONFA_constr -- the list of constraints
|
|
CONFA_obsvar -- the list of observed variables
|
|
CONFA_wgt -- weight specification
|
|
CONFA_wgti -- iweight
|
|
CONFA_t -- temporary global for -listutil-
|
|
|
|
Structure matrix:
|
|
CONFA_Struc -- the model structure: (parameter type, equation number, index1, index2)
|
|
|
|
History:
|
|
v.1.0 -- Jan 09, 2008
|
|
-- basic formulation without -cluster-, -robust-, -weights-, -svy-
|
|
v.1.1 -- Mar 21, 2008
|
|
-- Satorra-Bentler?
|
|
v.1.2 -- Sep 16, 2008
|
|
-- Ken Higbee comments
|
|
v.1.5 -- usenames
|
|
-- Mata moved to lconfa.mlib
|
|
-- survey-compatible
|
|
v.1.6 -- listwise deletion for missing data
|
|
-- what kind of idiot should Stas be to not pay attention to this???
|
|
-- informative message about matsize in Satorra-Bentler calculations
|
|
v.2.0 -- FIML missing data
|
|
-- prepared for revision in SJ
|
|
v.2.0.1 -- fixed -if- in Satorra-Bentler calculations
|
|
v.2.0.2 -- fixed reporting of correlations with -unitvar-: confa.ado, confa_estat.ado
|
|
v.2.1 -- someday?
|
|
-- Bartlett correction: (N - 1 - (2p+4m+5)/6)
|
|
-- F-statistic in place of chi-square, both normal theory and S-B
|