*! version 2.0, sean f. reardon, 30dec2005 /*------------------------------------------------- The hlmrun program runs an existing hlm command file. estimates from the model are returned as e-class -------------------------------------------------*/ capture program drop hlmrun program define hlmrun, eclass ereturn clear version 8.2 syntax anything (id="hlm cmdfile" name=filename) /// [, hlm2 hlm3 noView MDM(string) STore(string) ROBUST] gettoken path ext: filename, p(".") if "`ext'"=="" local filename "`filename'.hlm" else if "`ext'" ~= ".hlm" { di in re "file must be an .hlm file" error 603 } confirm file `"`filename'"' if "`hlm2'"=="hlm2" & "`hlm3'"=="hlm3" error 197 else if "`hlm2'"=="" & "`hlm3'"=="" local hlm2 "hlm2" if "`mdm'"=="" local mdm "$mdmfile" if "`mdm'"=="" { di in re "no mdm file specified" exit } *----------------------------------------------------------------- *NOTE if $varsinorder == "", then a problem reading in vc matrices *from models with > 60 fixed effects estimates *but not easy to get varsinorder from .hlm commmand file *----------------------------------------------------------------- *extract modeltype and variables from hlm command file *-------------------------------------------------- tempname cmd file open `cmd' using "`filename'", read file read `cmd' line while r(eof) == 0 { if `"`=lower(substr("`line'",1,7))'"' == "nonlin:" { loc modeltype `"`=trim(lower(subinstr("`line'","nonlin:","",1)))'"' tokenize "`modeltype'", parse(",") loc modeltype "`1'" if "`modeltype'" == "n" loc modeltype "linear" } else if `"`=lower(substr("`line'",1,7))'"' == "level1:" { loc out `"`=trim(lower(subinstr("`line'","level1:","",1)))'"' tokenize "`out'", parse("=") loc out "`1'" } else if `"`=lower(substr("`line'",1,8))'"' == "laplace:" { loc laplace `"`=lower(substr("`line'",9,1))'"' } else if `"`=lower(substr("`line'",1,9))'"' == "laplace6:" { loc l6 `"`=lower(substr("`line'",10,1))'"' } else if `"`=lower(substr("`line'",1,11))'"' == "emlaplace5:" { loc eml5 `"`=lower(substr("`line'",12,1))'"' } else if `"`=lower(substr("`line'",1,9))'"' == "deviance:" { tokenize "`line'", parse(":") loc dev = `3' } else if `"`=lower(substr("`line'",1,3))'"' == "df:" { tokenize "`line'", parse(":") loc df = `3' } file read `cmd' line } file close `cmd' *determine which results are stored in e(b) and e(V) *--------------------------------------------------- if "`modeltype'" == "linear" { if "`store'"~="lin" & "`store'"~="" { di in re "store option invalid for linear model" exit 197 } if "`store'" == "" loc store = "lin" if "`robust'"=="robust" loc store = "r" } else if "`laplace'`l6'`eml5'" == "" { if "`store'"~="us" & "`store'"~="pa" & "`store'"~="" { di in re "store option invalid" exit 197 } if "`store'" == "" loc store = "us" if "`robust'" == "robust" loc store = "`store'r" } else if "`laplace'`l6'" ~= "" { if "`store'"~="us" & "`store'"~="pa" & "`store'"~="eml" /// & "`store'"~="l6" & "`store'"~="" { di in re "store option invalid for `modeltype' model" exit 197 if "`robust'" == "robust" & ("`store'" == "eml" | "`store'" == "l6" /// di in gr "note: robust standard errors not computed for laplace estimates" } if "`laplace'" == "y" & "`store'" == "" loc store = "eml" else if "`l6'" == "y" & "`store'" == "" loc store "l6" if "`robust'" == "robust" & "`store'" ~= "eml" & "`store'" ~= "l6" /// loc store = "`store'r" } *temporarily rename existing tau & vc files *------------------------------------------ local d = 0 loc rc = 0 foreach f in gamvc gamvcr gamvcus gamvcusr gamvcpa gamvcpar /// gamvceml gamvcl tauvc tauvceml tauvcl6 { capture confirm new file "`f'.dat" loc rc = `rc' + _rc } while `rc' != 0 { loc rc = 0 local d = `d' + 1 foreach f in gamvc gamvcr gamvcus gamvcusr gamvcpa gamvcpar /// gamvceml gamvcl tauvc tauvceml tauvcl6 { capture confirm new file "`f'_temp`d'.dat" loc rc = `rc' + _rc } } if `d' != 0 { foreach f in gamvc gamvcr gamvcus gamvcusr gamvcpa gamvcpar /// gamvceml gamvcl tauvc tauvceml tauvcl6 { capture confirm file "`f'.dat" if _rc==0 { copy "`f'.dat" "`f'_temp`d'.dat" erase "`f'.dat" } } } global vctemp = `d' *run hlm command file to fit HLM model *------------------------------------- !`hlm2'`hlm3' "`mdm'" "`filename'" if "`view'"=="" view `path'.txt *store the estimates of gammas and their VC matrices *------------------------------------------------------------------ loc nvars: word count $varsinorder loc matlines=int(`=`nvars'+59'/60) foreach f in null r us usr pa par eml l { if "`f'" == "null" loc f "" capture confirm file "gamvc`f'.dat" if _rc == 0 { tempname bvc b_`f' V_`f' file open `bvc' using "gamvc`f'.dat", read loc gammas "" loc l = 1 while `l' <= `matlines' { file read `bvc' line loc gammas "`gammas' `line'" loc l = `l' + 1 } loc vcmat "" loc r = 1 while `r' <= `nvars' { loc l = 1 while `l' <= `matlines' { file read `bvc' line loc vcmat "`vcmat' `line'" loc l = `l' + 1 } if `r' < `nvars' loc vcmat "`vcmat' \" loc r = `r' + 1 } file close `bvc' capture erase "gamvc`f'.dat" if "`f'" == "" loc f "lin" if "`f'" == "l" loc f "l6" matrix input b_`f' = (`gammas') matrix input V_`f' = (`vcmat') matrix rownames b_`f' = "y1" matrix colnames b_`f' = $varsinorder matrix rownames V_`f' = $varsinorder matrix colnames V_`f' = $varsinorder matrix se_`f'=vecdiag(cholesky(diag(vecdiag(V_`f')))) } } capture confirm matrix b_`store' if _rc == 0 { ereturn post b_`store' V_`store', depname(`out') matrix b_`store' = e(b) matrix V_`store' = e(V) loc disp "y" } else { di in re /// `"note: requested results ("`store'") not produced, e(V) and e(b) not stored"' loc disp "n" } foreach f in lin r us usr pa par eml l6 { capture confirm matrix b_`f' if _rc == 0 { if "`f'" != "r" & "`f'" != "usr" & "`f'" != "par" /// ereturn matrix b_`f' = b_`f' else matrix drop b_`f' ereturn matrix V_`f' = V_`f' ereturn matrix se_`f' = se_`f' } } *store the tau matrices *---------------------- capture confirm file "tauvc.dat" if _rc == 0 & "`hlm2'" == "hlm2" { tempname tau file open `tau' using "tauvc.dat", read file read `tau' line loc rtau: word count `line' loc taumat "" loc r = 1 while `r' < `rtau' { loc taumat "`taumat'`line' \ " file read `tau' line loc r = `r' + 1 } loc taumat "`taumat'`line'" matrix input tau2 = (`taumat') ereturn matrix tau2 = tau2 * - - - - - - - - - - - - - - - - - - - - - - - - - *need to check dimension of tau information matrix *before activating the following lines: *tau information matrix not included in HLM2 tauvc.dat file? * loc rdim = `=`rtau'*(`rtau'+1)/2' * loc tauimat "" * loc r = 1 * file read `tau' line * while `r' < `rdim' { * loc tauimat "`tauimat'`line' \ " * file read `tau' line * loc r = `r' + 1 * } * loc tauimat "`tauimat'`line'" * matrix input invtaui = (`tauimat') * matrix se_tau = vecdiag(cholesky(diag(vecdiag(invtaui)))) * ereturn matrix invtaui = invtaui * ereturn matrix se_tau = se_tau * - - - - - - - - - - - - - - - - - - - - - - - - - file read `tau' line loc sig2chk: word count `line' if `sig2chk' == 1 ereturn scalar sig2 = `line' file close `tau' capture erase "tauvc.dat" } else if _rc == 0 & "`hlm3'" == "hlm3" { tempname tau file open `tau' using "tauvc.dat", read file read `tau' line loc rtau2: word count `line' loc tau2mat "" loc r = 1 while `r' < `rtau2' { loc tau2mat "`tau2mat'`line' \ " file read `tau' line loc r = `r' + 1 } loc tau2mat "`tau2mat'`line'" file read `tau' line loc rtau3: word count `line' loc tau3mat "" loc r = 1 while `r' < `rtau3' { loc tau3mat "`tau3mat'`line' \ " file read `tau' line loc r = `r' + 1 } loc tau3mat "`tau3mat'`line'" matrix input tau2 = (`tau2mat') matrix input tau3 = (`tau3mat') ereturn matrix tau2 = tau2 ereturn matrix tau3 = tau3 * - - - - - - - - - - - - - - - - - - - - - - - - - *need to check dimension of tau information matrix *before activating the following lines: loc rdim = `=(`rtau2'*(`rtau2'+1)/2)+(`rtau3'*(`rtau3'+1)/2)' loc tauimat "" loc r = 1 file read `tau' line while `r' < `rdim' { loc tauimat "`tauimat'`line' \ " file read `tau' line loc r = `r' + 1 } loc tauimat "`tauimat'`line'" matrix input invtaui = (`tauimat') matrix se_tau = vecdiag(cholesky(diag(vecdiag(invtaui)))) ereturn matrix invtaui = invtaui ereturn matrix se_tau = se_tau * - - - - - - - - - - - - - - - - - - - - - - - - - file read `tau' line loc sig2chk: word count `line' if `sig2chk' == 1 ereturn scalar sig2 = `line' file close `tau' capture erase "tauvc.dat" } capture confirm file "tauvceml.dat" if _rc == 0 { tempname taueml file open `taueml' using "tauvceml.dat", read file read `taueml' line loc rtau: word count `line' loc taumat "" loc r = 1 while `r' < `rtau' { loc taumat "`taumat'`line' \ " file read `taueml' line loc r = `r' + 1 } loc taumat "`taumat'`line'" file close `taueml' matrix input tau2_eml = (`taumat') ereturn matrix tau2_eml = tau2_eml capture erase "tauvceml.dat" } capture confirm file "tauvcl6.dat" if _rc == 0 { tempname taul6 file open `taul6' using "tauvcl6.dat", read file read `taul6' line loc rtau: word count `line' loc taumat "" loc r = 1 while `r' < `rtau' { loc taumat "`taumat'`line' \ " file read `taul6' line loc r = `r' + 1 } loc taumat "`taumat'`line'" matrix input tau2_l6 = (`taumat') ereturn matrix tau2_l6 = tau2_l6 * - - - - - - - - - - - - - - - - - - - - - - - - - *need to check dimension of tau information matrix *before activating the following lines: * loc rdim = `=`rtau'*(`rtau'+1)/2' * loc tauimat "" * loc r = 1 * file read `taul6' line * while `r' < `rdim' { * loc tauimat "`tauimat'`line' \ " * file read `taul6' line * loc r = `r' + 1 * } * loc tauimat "`tauimat'`line'" * matrix input invtaui_l6 = (`tauimat') * matrix se_tau_l6 = vecdiag(cholesky(diag(vecdiag(invtaui_l6)))) * ereturn matrix invtaui_l6 = invtaui_l6 * ereturn matrix se_tau_l6 = se_tau_l6 * - - - - - - - - - - - - - - - - - - - - - - - - - file close `taul6' capture erase "tauvcl6.dat" } *store the deviance statistics *----------------------------- mac drop df dev dev_us dev_pa dev_eml dev_l6 tempname outfile file open `outfile' using `"`path'.txt"', read file read `outfile' line loc line: subinstr loc line `"'"' `"""', all *set trace on while r(eof) == 0 { tokenize `"`line'"' if `"`=lower(substr(trim(`"`line'"'),1,22))'"' == "the maximum number of " { if "`5'" == "level-1" ereturn scalar N_1 = `8' if "`5'" == "level-2" ereturn scalar N_2 = `8' if "`5'" == "level-3" ereturn scalar N_3 = `8' } if "`=lower(substr(`"`1'"',2,8))'" == "results" { loc w = 1 while "``w''" ~= "" { if "`=lower(`"``w''"')'" == "unit-specific" loc results = "_us" if "`=lower(`"``w''"')'" == "population" loc results = "_pa" if "`=lower(`"``w''"')'" == "laplace-6" loc results = "_l6" if "`=lower(`"``w''"')'" == "laplace-2" loc results = "_eml" loc w = `w' + 1 } file read `outfile' line loc line: subinstr loc line `"'"' `"""', all tokenize `"`line'"' loc w = 1 while "``w''" ~= "" { if "`=lower(`"``w''"')'" == "unit-specific" loc results = "_us" if "`=lower(`"``w''"')'" == "population" loc results = "_pa" if "`=lower(`"``w''"')'" == "laplace-6" loc results = "_l6" if "`=lower(`"``w''"')'" == "laplace-2" loc results = "_eml" loc w = `w' + 1 } } if "`1'" == "Deviance" { global dev`results' = `3' ereturn scalar dev`results' = ${dev`results'} } if `"`=lower(substr(trim(`"`line'"'),1,30))'"' /// == "number of estimated parameters" global df = `6' file read `outfile' line loc line: subinstr loc line `"'"' `"""', all } file close `outfile' ereturn local cmd "`hlm2'`hlm3'" ereturn local model "`modeltype'" ereturn local depvar "`out'" ereturn local cmdfile `"`filename'"' local setype "non-robust" if "`robust'" === "robust" & "`store'" ~= "eml" local setype "robust" ereturn local setype "`setype'" local esttype "" if "`store'" == "us" | "`store'" == "usr" local esttype "unit-specific" if "`store'" == "pa" | "`store'" == "par" local esttype "population-average" if "`store'" == "eml" local esttype "em laplace" if "`store'" == "l6" local esttype "laplace-6" ereturn local esttype "`esttype'" if "`esttype'" ~= "" local esttype " `esttype'" if "$df" ~= "" ereturn scalar df = $df if "`disp'" == "y" { di in gr "dependent variable: " in ye "`out'" di in gr "model type: " in ye "`modeltype'" di in gr "results: " in ye "estimated`esttype' gammas with `setype' standard errors" ereturn display } if "`dev'" ~= "" & "`e(df)'" ~= "" { loc results: subinstr loc store "r" "", all loc results "_`results'" if "`results'" == "_lin" local results "" di if (`e(dev`results')'-`dev')*(`e(df)'-`df')>0 /// di in re "Caution: Models are not nested" di in gr "{col 25}Deviance{col 40}df" di in gr "--------------------------------------------- di in gr "Current model:{col 25}" e(dev`results') " {col 40}" e(df) di in gr "Comparison model:{col 25}" `dev' " {col 40}" `df' di in gr "Difference:{col 25}" in ye abs(`e(dev`results')'-`dev') /// " {col 40}" abs(`e(df)'-`df') di di in ye "P-value:{col 24}" /// 1-chi2(abs(`e(df)'-`df'),abs(`e(dev`results')'-`dev')) di in gr "--------------------------------------------- } *drop new tau and vc files and rename old ones *--------------------------------------------- if $vctemp != 0 { foreach f in gamvc gamvcr gamvcus gamvcusr gamvcpa gamvcpar /// gamvceml gamvcl tauvc tauvceml tauvcl6 { capture confirm file "`f'_temp${vctemp}.dat" if _rc==0 copy "`f'_temp${vctemp}.dat" "`f'.dat" if _rc==0 erase "`f'_temp${vctemp}.dat" } } mac drop vctemp end