*! version 2.3.20 SRH 7 Sept 2011 program define gllamm, eclass version 7.0 timer on 1 if replay() { if "`e(cmd)'" ~= "gllamm" { di in red "gllamm was not the last command" error 301 } Replay `0' } else { Estimate `0' } timer off 1 end program define procstr, eclass version 6.0 tempname bc b Vc Vr V ll esamp noi disp "processing constraints" scalar `ll' = e(ll) local df = e(df_m) local dof if "`df'"~="."{ local dof "dof(`df')" } local k = e(k) capture matrix `Vr' = e(Vr) capture robclus "`e(robclus)'" matrix `bc' = e(b) matrix `Vc' = e(V) local y = "`e(depvar)'" matrix `b' = `bc'*M_T' + M_a matrix colnames `b' = $HG_coln matrix coleq `b' = $HG_cole * matrix list `b' * matrix list `Vc' matrix `V' = M_T*`Vc'*M_T' * disp "computed V" * matrix list `V' estimates post `b' `V' M_C, $HG_obs `dof' est local ll =`ll' est local k = `k' est local depvar "`y'" capture est matrix Vr `Vr' capture est local robclus "`robclus'" * disp "posted results" end program define Replay, eclass version 6.0 syntax [, Level(real $S_level) EFORM ALLC ROBUST CLUSTER(varname) DOTS SCorefil(passthru) noROb noDISplay EVAL] tempname M_nffc M_nu Vs if "`eval'"==""{ if "`e(pweight)'"~="" { local robust robust } if "`robust'"~=""|"`cluster'"~=""{ if "`cluster'"~=""{ local cluster cluster(`cluster') } gllarob, `cluster' `dots' `scorefil' `rob' } if ("`robust'"==""|"`rob'"~="")&"`cluster'"==""{ * disp "reposting model-based standard errors" matrix `Vs' = e(Vs) estimates repost V =`Vs' } } if "`display'"==""{ local const = e(const) local tplv = e(tplv) matrix `M_nffc' = e(nffc) capture matrix `M_nu' = e(nu) * capture matrix list `M_nu' if _rc == 0{ disp " " local j = 1 while `j' <= `tplv'{ disp in gr "number of level `j' units = " in ye `M_nu'[1,`j'] local j = `j' + 1 } disp " " } local cn = e(cn) if `cn'>0{ disp in gr "Condition Number = " in ye `cn' } else{ disp in gr "Condition Number could not be computed" } disp " " * matrix list `M_nffc' local numeq = colsof(`M_nffc') if `M_nffc'[1,1]==0{ local numeq = `numeq' -1 } if `numeq' > 1{ local first neq(`numeq') } else{ local first first } local primess if `e(prior)' == 1 { local primess "with prior " } if e(ll_0)==.|`M_nffc'[1,1]==0{ local nohead "noheader" if `const'==0{ disp in gr "gllamm model `primess'" } else{ disp in gre "gllamm model `primess'with constraints" *matrix dispCns } disp " " if `e(prior)' == 1 { disp in gr "log posterior = " in ye e(ll) disp in gr "log likelihood = " in ye e(a_ll) } else{ disp in gr "log likelihood = " in ye e(ll) } } if "`eform'"~=""{ local eform "eform(exp(b))" } disp " " if "`cluster'"~=""|("`robust'"~=""&"`rob'"==""){ if "`cluster'"~=""{ disp "Robust standard errors for clustered data: `cluster'" } else{ disp "Robust standard errors" } } if `M_nffc'[1,1]>0|`numeq'>0 { quietly q born if $S_1>15502{ /* version 8 or later */ _coef_table, level(`level') `first' `eform' } else{ /* older versions */ if `const' == 0{ noi ml display, level(`level') `nohead' `first' `eform' } else{ matrix dispCns noi estimates display, level(`level') `first' `eform' } } } else{ disp in gr "No fixed effects" } noi disprand if "`allc'"~=""{ if `const' == 0{ noi ml display, level(`level') `nohead' } else{ noi estimates display, level(`level') } } } end program define Estimate, eclass version 6.0 syntax varlist(min=1) [if] [in] , I(string) [NRf(numlist integer min=1 >=1) /* */ Eqs(string) GEqs(string) PEqs(string) /* */ noCORrel noCOnstant BMATrix(string) INTER(string) FRLoad(numlist integer min=1 >=1) /* */ Family(string) DEnom(varname numeric min=1) Link(string) EXpanded(string) /* */ Offset(varname numeric) Basecategory(integer 999) /* */ THresh(string) ETHresh(string) COMPosite(varlist numeric min=3) * ] *disp in re "`opt'" *disp in re "1: `varlist'" local origif `"`if'"' local origin `"`in'"' local origvar `"`varlist'"' local opt "`options'" local 0 "`origvar' `origif' `origin', `opt'" syntax varlist(min=1) [if] [in] [, /* */ Weightf(string) LV(varname numeric min=1) FV(varname numeric min=1) S(string) /* */ NATS IP(string) NIp(numlist integer min=1 >=1) ADapt Constraints(numlist) /* */ FRom(string) LONG Gateaux(passthru) LF0(passthru) SEarch(passthru) /* */ ROBust CLuster(varname) PWeight(string) ITerate(passthru)/* */ DOts noLOg TRace noDISplay noESt EVal Level(real $S_level) INit noDIFficult /* */ EFORM ALLC ADOONLY SKIP COPY PRior(passthru) ] /* does not allow extra options! */ *disp in re "2: `varlist'" *disp in re "i: `i'" *disp in re "frload: `frload'" *disp in re "display: `display'" *disp in re "iterate: `iterate'" tempname mat mnip mnbrf global HG_error=0 /* deal with adoonly */ global HG_noC = 0 global HG_noC1 = 0 global HG_noC2 = 0 if "`adoonly'"=="" { qui q born if $S_1 < 16350 { global HG_noC2 = 1 if $S_1 < 15722 { noi di noi di in gr "You must have the Stata 8 executable updated on or after " _c noi di in ye %d 15722 noi di in gr "in order to use all internal routines" global HG_noC1 = 1 if $S_1 < 15274 { global HG_noC = 1 noi di in gr "Option adoonly assumed." noi di } } } } if "`adoonly'"!="" { global HG_noC = 1 /* link and family */ global HG_noC1 = 1 /* lnupdate */ global HG_noC2 = 1 /* don't want to set HG_noC to 0 */ } /* deal with trace */ if "`trace'"!="" { local noi "noisily" } /* deal with dots */ global HG_dots = 0 if "`dots'"!="" { global HG_dots = 1 } /* deal with init */ global HG_init=0 if "`init'"~="" { global HG_init=1 } /* deal with if and in */ marksample touse qui count if `touse' if _result(1) <= 1 { di in red "insufficient observations" exit 2001 } /* deal with varlist */ tokenize `varlist' local y "`1'" macro shift /* `*' is list of dependent variables */ local indep "`*'" local num: word count `indep' /* number of independent variables */ markout `touse' `y' `indep' /* deal with Link and Family */ global HG_lev1=0 global HG_famil global HG_linko global HG_link matrix M_olog=(0) capture matrix drop M_oth global HG_mlog=0 global HG_nolog = 0 global HG_lv global HG_fv global HG_smlog=0 global HG_oth = 0 local l: word count `family' if `l'>1 { `noi' qui disp "more than one family" if "`fv'"==""{ disp in re "need fv option" exit 198 } else{ confirm variable `fv' global HG_fv `fv' } parse "`family'", parse(" ") local n=1 while "`1'"~=""{ qui count if `fv'==`n' if _result(1)==0{ disp "family `1' not used" } fm "`1'" if "$S_2"=="gauss"{ if $HG_lev1==0{ global HG_lev1=1 } else if $HG_lev1==2{ global HG_lev1=3 } } else if "$S_2"=="gamma"{ if $HG_lev1==0{ global HG_lev1=2 } else if $HG_lev1==1{ global HG_lev1=3 } } global HG_famil "$HG_famil $S_2" local n = `n'+1 mac shift } } local k: word count `link' local mll = 0 if `k'>1{ `noi' qui disp "more than one link" if "`lv'"==""{ disp in re "need lv option" exit 198 } else{ confirm variable `lv' global HG_lv `lv' } parse "`link'", parse(" ") local n=1 while "`1'"~=""{ qui count if $HG_lv==`n' if _result(1)==0{ disp "link `1' not used" } lnk "`1'" if "$S_1"=="sprobit"|"$S_1"=="soprobit"{ if $HG_lev1 == 2{ global HG_lev1 = 3 } else{ global HG_lev1 = 1 } } if "$S_1"=="ll"{ global HG_noC = $HG_noC2 } /* nominal */ if "$S_1"=="mlogit"|"$S_1"=="smlogit"{ if $HG_mlog>0{ disp in re "can only have one mlogit link" exit 198 } global HG_mlog=`n' if "$S_1"=="smlogit"{ if $HG_lev1 == 2{ global HG_lev1 = 3 } else{ global HG_lev1 = 1 } } tempvar first sort `touse' $HG_lv `y' qui by `touse' $HG_lv `y': gen byte `first' = cond(_n==1,1,0) mkmat `y' if `first' == 1 & `touse' & $HG_lv == `n', mat(M_respm) if "$S_1"=="smlogit"{global HG_smlog=1} } /* ordinal */ else if "$S_1"=="ologit"|"$S_1"=="oprobit"|"$S_1"=="ocll"|"$S_1"=="soprobit"{ global HG_linko "$HG_linko $S_1" if $HG_nolog>0{ * disp "more than one ordinal response" matrix M_olog = M_olog,`n' } else{ capture matrix drop M_nresp matrix M_olog[1,1] = `n' tempvar first sort `touse' $HG_lv `y' qui by `touse' $HG_lv `y': gen byte `first' = cond(_n==1,1,0) } mkmat `y' if `first' == 1 & `touse' & $HG_lv == `n', mat(`mat') local ll = rowsof(`mat') * matrix list `mat' * disp "adding `ll' to M_nresp" matrix M_nresp = nullmat(M_nresp),`ll' if `ll'>`mll'{ local mll = `ll' } global HG_nolog = $HG_nolog + 1 } /* other */ else { global HG_link "$HG_link $S_1" matrix M_oth = nullmat(M_oth),`n' global HG_oth=1 } local n = `n'+1 mac shift } if $HG_nolog>0{ tempname junk global HG_lvolo "`junk'" qui gen byte $HG_lvolo = 0 matrix M_resp = J(`mll',$HG_nolog,0) local no = 1 local totresp = 0 while `no'<=$HG_nolog{ local olog = M_olog[1,`no'] qui replace $HG_lvolo = 1 if $HG_lv == `olog' mkmat `y' if `first' == 1 & `touse' & $HG_lv == `olog', mat(`mat') local ii = 1 while `ii'<= M_nresp[1,`no']{ * disp "M_resp[`ii',`no'] = mat[`ii',1]" matrix M_resp[`ii',`no'] = `mat'[`ii',1] local ii = `ii' + 1 } local totresp = `totresp' + M_nresp[1,`no'] local no = `no' + 1 } } if $HG_mlog>0{ if $HG_nolog==0{ tempname junk global HG_lvolo "`junk'" qui gen byte $HG_lvolo = 0 } qui replace $HG_lvolo = 1 if $HG_lv == $HG_mlog } } else if `k'<=1&`l'<=1{ /* no more than one link and one family given */ lnkfm "`link'" "`family'" global HG_link = "$S_1" global HG_famil = "$S_2" if "$HG_link"=="ologit"|"$HG_link"=="oprobit"|"$HG_link"=="ocll"|"$HG_link"=="soprobit"{ global HG_linko = "$HG_link" global HG_nolog = 1 matrix M_olog[1,1] = 1 } if "$HG_link"=="smlogit"|"$HG_link"=="mlogit"{global HG_mlog=1} if "$HG_famil"=="gauss"{global HG_lev1=1} if "$HG_famil"=="gamma"{global HG_lev1=2} if "$HG_link"=="sprobit"{global HG_lev1=1} if "$HG_link"=="soprobit"{global HG_lev1=1} if "$HG_link"=="smlogit"{global HG_lev1=1} if "$HG_link"=="ll"{global HG_noC = $HG_noC2 } if $HG_mlog==0&$HG_nolog==0{global HG_oth = 1} } else if `k'==1{ lnk "`lnk'" global HG_link = "$S_1" if "$HG_link"=="ologit"|"$HG_link"=="oprobit"|"$HG_link"=="ocll"|"$HG_link"=="soprobit"{ global HG_nolog = 1 matrix M_olog[1,1] = 1 global HG_linko = "$HG_link" } if "$HG_link"=="smlogit"|"$HG_link"=="mlogit"{global HG_mlog=1} if "$HG_link"=="sprobit"{global HG_lev1=1} if "$HG_link"=="smlogit"{global HG_lev1=1} if "$HG_link"=="soprobit"{global HG_lev1=1} if "$HG_link"=="ll"{global HG_noC = $HG_noC2 } if $HG_mlog==0&$HG_nolog==0{global HG_oth = 1} } if `l'==1{ fm "`family'" global HG_famil = "$S_2" if "$HG_famil"=="gauss"{global HG_lev1=1} if "$HG_famil"=="gamma"{global HG_lev1=2} if $HG_mlog==0&$HG_nolog==0{global HG_oth = 1} } if ((`k'>1&`l'==0)|(`l'>1&`k'==0))&$HG_oth==1{ disp in re /* */ "both link() and fam() required for multiple links or families" exit 198 } markout `touse' $HG_lv $HG_fv /* deal with noCORrel */ global HG_cor = 1 if "`correl'"~=""{ global HG_cor = 0 } /* deal with DEnom */ global HG_denom local f=0 parse "$HG_famil", parse(" ") while "`1'"~=""&`f'==0{ if "`1'"=="binom"{ local f=1 } mac shift } if `f'==1{ if "`denom'"~=""{ confirm variable `denom' global HG_denom "`denom'" } else{ tempvar den quietly gen byte `den'=1 global HG_denom "`den'" } } else{ if "`denom'"~=""{ disp in blue/* */"option denom(`denom') given but binomial family not used" } } markout `touse' `denom' /* deal with offset */ global HG_off if "`offset'"~=""{ global HG_off "`offset'" local offset "offset(`offset')" } markout `touse' $HG_off /* deal with ip */ global HG_gauss = 1 global HG_free = 0 global HG_cip = 1 global HG_mult = 0 if "`ip'"~=""{ if "`ip'"=="g"{ global HG_gauss = 1 } else if "`ip'"=="l"{ global HG_gauss = 0 } else if "`ip'"=="f"{ global HG_free = 1 } else if "`ip'"=="m"{ global HG_mult = 1 global HG_gauss = 0 } else if "`ip'"=="fn"{ global HG_free = 1 global HG_cip = 0 } else { disp in re "ip option `ip' not valid" exit 198 } } global HG_lzpr lzprobg global HG_zip zipg * disp in re "HG_mult = " $HG_mult if $HG_mult{ global HG_lzpr lzprobm } else if $HG_free{ global HG_lzpr lzprobf global HG_zip zipf } /* deal with expanded */ global HG_ind global HG_exp = 0 global HG_expf = 0 if "`expanded'"~=""{ global HG_exp = 1 if $HG_mlog==0{ disp in re "expanded option only valid with mlogit link" exit 198 } local k: word count `expanded' if `k'~=3{ disp in re "expanded option must have three arguments" } local exp: word 1 of `expanded' confirm variable `exp' global HG_mlg `exp' local k: word 2 of `expanded' * 11/11/06 confirm variable `k' global HG_ind `k' local k: word 3 of `expanded' if "`k'"=="o"{ global HG_expf=1 } else{ if "$HG_link"~="mlogit"&"$HG_link"~="smlogit"{ disp in re "must use o in expanded option when combining mlogit with other links" exit 198 } } * 11/11/06 markout `touse' $HG_mlg $HG_ind } else{ if $HG_mlog>0&"$HG_link"~="mlogit"&"$HG_link"~="smlogit"{ disp in re "must use expanded option when combining mlogit with other links" exit 198 } tempvar ind gen byte `ind' = 1 global HG_ind `ind' global HG_exp = 0 } /* deal with composite */ global HG_comp = 0 global HG_coall if "`composite'"~=""{ local k: word count `composite' global HG_comp = `k' - 2 local k: word 1 of `composite' global HG_mlg `k' local k: word 2 of `composite' global HG_ind `k' local kk = 1 while `kk'<= $HG_comp { local jj = `kk' + 2 local k: word `jj' of `composite' global HG_co`kk' `k' global HG_coall $HG_coall `k' local kk = `kk' + 1 } global HG_noC = $HG_noC2 /* use ado-code for link if not latest Stata 8 */ } /* deal with I (turn list around)*/ if ("`i'"==""){ disp in red "i() required" global HG_error=1 exit 198 } local tplv: word count `i' global HG_tplv = `tplv'+1 global HG_clus local k = `tplv' while `k'>=1{ local clus: word `k' of `i' confirm numeric variable `clus' markout `touse' `clus', strok local k=`k'-1 global HG_clus "$HG_clus `clus'" } if "`expanded'"==""&"`composite'"==""{ tempvar id gen long `id'=_n global HG_clus "$HG_clus `id'" } else{ global HG_clus "$HG_clus $HG_mlg" } /* deal with weightf */ tempvar wt quietly gen double `wt'=1 local j = 1 if "`weightf'"==""{ while (`j'<=$HG_tplv){ tempname junk global HG_wt`j' "`junk'" gen double ${HG_wt`j'}=1 local j = `j' + 1 } global HG_weigh } else{ global HG_weigh "`weightf'" local wtvars local found = 0 while (`j'<=$HG_tplv){ capture confirm variable `weightf'`j' /* frequency weight */ if _rc ~= 0 { tempname junk global HG_wt`j' "`junk'" gen double ${HG_wt`j'} = 1 /* can become non-integeger later */ } else{ tempname junk global HG_wt`j' "`junk'" gen double ${HG_wt`j'}=`weightf'`j' local wtvars `wtvars' `weightf'`j' quietly replace `wt'=`wt'*${HG_wt`j'} local found = `found' + 1 } local j = `j' + 1 } if `found' == 0 { disp in red "weight variables `weightf' not found" global HG_error=1 exit 111 } markout `touse' `wtvars' } /* deal with probability weights */ local pw local wtvars if "`pweight'"~=""{ local robust robust tempname wtp local pw pweight global HG_pwt "`pweight'" quietly gen double `wtp' = 1 local j = 1 local found = 0 while (`j'<=$HG_tplv){ capture confirm variable `pweight'`j' /* probability weight */ if _rc == 0 { quietly replace `wtp'=`wtp'*`pweight'`j' local wtvars `wtvars' `pweight'`j' local found = `found' + 1 } local j = `j' + 1 } if `found' == 0 { disp in red "probability weight variables not found" global HG_error=1 exit 111 } markout `touse' `wtvars' } /* deal with cluster */ if "`cluster'"~=""{ * check that top-level cluster does not vary within cluster? local k: word count $HG_clus local top: word 1 of $HG_clus qui sort `top' `cluster' capture qui by `top' : assert `cluster'[1]==`cluster'[_N] if _rc>0{ disp in re "`cluster' varies within `top'" exit 198 } markout `touse' `cluster' local cluster cluster(`cluster') } /* deal with categorical response variables */ if "$HG_link" == "mlogit"|"$HG_link" == "smlogit"{ sort `touse' `y' tempvar first qui by `touse' `y': gen byte `first' = cond(_n==1,1,0) mkmat `y' if `first' == 1 & `touse', mat(M_respm) } else if /* */ "$HG_link" == "ologit"|"$HG_link" == "ocll"|"$HG_link" == "oprobit"|"$HG_link"=="soprobit"{ sort `touse' `y' tempvar first qui by `touse' `y': gen byte `first' = cond(_n==1,1,0) mkmat `y' if `first' == 1 & `touse', mat(M_resp) local totresp = rowsof(M_resp) matrix M_nresp = (`totresp') } /* deal with base-category */ if `basecategory'~=999{ if "$HG_link" ~= "mlogit"&"$HG_link" ~= "smlogit"{ disp in red "basecategory ignored because response not nominal" } } if $HG_mlog>0&$HG_expf==0{ tempname bas if `basecategory'==999{ scalar `bas' = M_respm[1,1] matrix `bas' = (`bas') local basecat = M_respm[1,1] * disp in re "`basecat'" } else{ matrix `bas' = (`basecategory') local basecat = `basecategory' } local n = rowsof(M_respm) local j = 1 local found = 0 while `j'<=`n'{ local el = M_respm[`j',1] if `el'==`basecat'{ local found = 1 } else{ matrix `bas' = `bas'\ `el' } local j = `j' + 1 } if `found' == 0 { disp in re "basecategory = `basecat' not one of the categories" exit 198 } matrix M_respm = `bas' local el = M_respm[1,1] local basecat basecat(`el') } /* deal with noCOns */ if "`constant'"~=""{ if $HG_nolog>0{ disp in re "noconstant option not allowed for ordinal links:" _n "This is the default because all thresholds estimated" exit 198 } local cns } else{ if $HG_cip == 0{ disp in re "are you sure you need a constant with ip(fn) option?" } local num=`num'+1 local cns "_cons" } matrix M_nffc=(`num') if `num'>0 { global HG_fixe (`y': `y'=`indep', `constant') local dep } else{ global HG_fixe local dep "`y'=" } /* fixed effects matrix */ global HG_ethr = 0 if `num' > 0 { matrix M_initf=J(1,`num',0) matrix coleq M_initf=`y' matrix colnames M_initf=`indep' `cns' } else{ cap matrix drop M_initf } if $HG_nolog==0{ if "`thresh'"~=""{ disp in re "thresh option ignored" } if "`ethresh'"~=""{ disp in re "ethresh option ignored" } } else if $HG_nolog>0{ if "`ethresh'"~=""{ if "`thresh'"~=""{ disp in re "cannot use both ethresh() and thresh() options! exit 198 } global HG_ethr = 1 local thresh `ethresh' } if "`thresh'"~=""{ local k: word count `thresh' if `k'~=$HG_nolog{ disp in re "number of threshold equations should be " $HG_nolog exit 198 } } global HG_fixe local n = rowsof(M_resp) matrix M_nffc[1,1] = `num'-1 if `num'>1{ global HG_fixe (`y': `y'=`indep', nocons) matrix `mat' = M_initf[1,1..`num'-1] local ce: coleq(`mat') local cn `indep' matrix M_initf=J(1,`num'-1,0) } else{ capture matrix drop M_initf } local el = M_nffc[1,1] local ii = 1 local nxt = M_nffc[1,1] + 1 local ntr = 1 local vars local rhs "_cons" while `ii'<= $HG_nolog{ local j = 1 if "`thresh'"~=""{ local eqnam: word `ii' of `thresh' eq ? "`eqnam'" local vars "$S_1" markout `touse' `vars' local ntr: word count `vars' local ntr = `ntr' + 1 local rhs "`vars' _cons" } while `j'< M_nresp[1,`ii']{ * disp "`ii'th ordinal response, level `j'" local el = `el' + `ntr' matrix M_nffc = M_nffc, `el' matrix `mat'=J(1,`ntr',0) matrix coleq `mat' = _cut`ii'`j' local cee: coleq(`mat') local ce `ce' `cee' local cn `cn' `rhs' global HG_fixe $HG_fixe (_cut`ii'`j':`vars') if `j' == 1 & `ii'==1 & `num' == 1{ global HG_fixe (_cut`ii'`j':`y'= `vars') } local j = `j' + 1 if $HG_ethr{ matrix `mat'[1,`ntr'] = - 0.5 } else{ matrix `mat'[1,`ntr'] = `j' - (M_nresp[1,`ii']+1)/2 } matrix M_initf = nullmat(M_initf), `mat' local nxt = `nxt' + 1 } local ii = `ii' + 1 } matrix colnames M_initf=`cn' matrix coleq M_initf=`ce' * matrix list M_initf } if ($HG_mlog>0)&$HG_expf==0{ global HG_fixe local n = rowsof(M_respm) matrix `bas'=M_initf matrix drop M_initf matrix drop M_nffc local j = 2 while `j'<=`n'{ local el = M_respm[`j',1] matrix coleq `bas' = c`el' matrix M_initf = nullmat(M_initf), `bas' matrix M_nffc = nullmat(M_nffc), (`j'-1)*`num' if `j' == 2{ global HG_fixe $HG_fixe ( c`el':`y' = `indep', `constant') } else{ global HG_fixe $HG_fixe ( c`el':`indep', `constant') } local j = `j' + 1 } local num = `num'*(`n' - 1) } * matrix list M_nffc * matrix list M_initf /* display information */ quietly `noi'{ disp " " disp in gr "General model information" *disp in smcl in gr "{hline 78}" _n di in gr _dup(78) "-" _n disp in gr "dependent variable:" in ye " `y'" if $HG_oth{ disp in gr "family:" in ye " $HG_famil" disp in gr "link:" in ye " $HG_link" } if "$HG_linko"~=""{ disp in gr "ordinal responses:" in ye " $HG_linko" } if $HG_mlog>0{ if $HG_smlog==1 { disp in gr "nominal responses:" in ye " smlogit" } else{ disp in gr "nominal responses:" in ye " mlogit" } } if "$HG_denom"~=""{ if "`denom'"~=""{ disp in gr "denominator:" in ye " `denom'" } else{ disp in gr "denominator:" in ye " 1" } } if "`offset'"~=""{ disp in gr "offset:" in ye " $HG_off" } local m = colsof(M_nffc) if `m'==1&M_nffc[1,1]>0{ local cuts: colnames(M_initf) disp in gr "equation for fixed effects " in ye " `cuts'" } else if `m'==1{ disp in gr "equation for fixed effects " in ye " none" } else{ disp in gr "equations for fixed effects" local j = 1 local nxt = 1 local prev = 0 while `j'<=`m'{ local n = M_nffc[1,`j'] - `prev' if `n'>0{ local prev = M_nffc[1,`j'] matrix `mat' = M_initf[1,`nxt'..`nxt'+`n'-1] local nxt = `nxt' + `n' local ce: coleq(`mat') local ce: word 1 of `ce' local cn: colnames(`mat') disp in gr " `ce': " in ye " `cn'" } local j = `j' + 1 } disp " " } } /* end qui `noi' */ /* deal with inter */ global HG_inter = 0 if "`inter'"~=""{ global HG_inter=1 local j: word count `inter' if `j'~=2{ disp in red "inter should have two arguments" exit 198 } local j: word 1 of `inter' capture confirm number `j' if _rc>0{ disp in red "arguments of inter must be numbers" exit 198 } global HG_l = `j' local j: word 2 of `inter' capture confirm number `j' if _rc>0{ disp in red "arguments of inter must be numbers" exit 198 } global HG_r = `j' } /* initialise macros */ quietly `noi' initmacs "`nrf'" "`nip'" "`eqs'" "`geqs'" "`peqs'" "`s'" "`nats'" "`bmatrix'" "`touse'" "`dep'" "`frload'" qui count if `touse' if _result(1) <= 1 { di in red "insufficient observations" exit 2001 } /* deal with noESt */ if "`est'"~=""{ exit 0 } /* only use observations satisfying if and in and having nonmissing values */ preserve quietly keep if `touse' /* work out number of units at each level */ qui summ `wt' if `touse', meanonly local lobs = r(sum) tempvar cw f qui gen double `cw' = `wt' qui gen byte `f' = 1 matrix M_nu=J(1,$HG_tplv,0) matrix M_nu[1,1]=`lobs' local sortlst $HG_clus local j = 1 local k = $HG_tplv quietly `noi' disp in gr "number of level 1 units = " in ye `lobs' while `j'<$HG_tplv{ *disp "sort `sortlst'" sort `sortlst' tokenize "`sortlst'" local `k' " " local sortlst "`*'" *disp "replace cw = cw/wt`j'" qui replace `cw' = `cw'/${HG_wt`j'} *disp "by `sortlst': replace f=_n==1" qui by `sortlst': replace `f' = _n==1 qui summ `cw' if `f' ==1, meanonly local lobs = r(sum) quietly `noi' disp in gr "number of level " `j'+1 " units = " in ye `lobs' matrix M_nu[1,`j'+1] = `lobs' local j = `j' + 1 local k = `k' - 1 } /* use pweights */ if "`pw'"~=""{ local j = 1 while (`j'<=$HG_tplv){ capture confirm variable `pweight'`j' /* probability weight */ if _rc == 0 { quietly replace ${HG_wt`j'}=${HG_wt`j'}*`pweight'`j' } local j = `j' + 1 } } /* check if weights are integer */ qui cap summ `y' if `touse' [fweight=`wt'], meanonly if _rc>0 { global HG_obs local ftype pweight disp in re "weights are non-integer" } else { local lobs = M_nu[1,1] global HG_obs obs(`lobs') local ftype fweight } if "`pweight'"~=""{ quietly replace `wt' = `wt'*`wtp' local ftype pweight } /* deal with from */ if "`from'"~=""{ capture qui matrix list `from' /* do not comment out this line! */ local rc=_rc if `rc'>1{ disp in red "`from' not a matrix" exit 111 } tempname amat matrix `amat' = `from' local from "`amat'" } /* deal with eval */ if "`eval'"~=""{ if "`from'"==""{ disp in re "eval not valid without from option" exit 198 } } /* deal with constraints (and from long) */ global HG_const = 0 if "`constra'"~=""{ global HG_const = 1 if $HG_init==0{ tempname b V matrix `b' = nullmat(M_initf), nullmat(M_initr) * noi matrix list `b' if "`from'"~=""& "`long'"~=""{ local nb = colsof(`b') local nf = colsof("`from'") * disp "nb = " `nb' * disp "nf = " `nf' if "`gateaux'"~=""{ local tprf=M_nrfc[2,$HG_tplv]-M_nrfc[2,$HG_tplv-1] local nnf = `nf' + `tprf' + 1 if `nnf'~=`nb'{ disp in re "from matrix has `nf' columns and should have " `nb'-`tprf'-1 exit 198 } matrix `from' = `from',`b'[1,`nf'+1...] ** * local lfix = `M_nffc'[1,colsof(`M_nffc')] * matrix `from' = `from'[1,1..`lfix'],`b'[1,`nf'+1..`nf'+`tprf'+1],`from'[1,`lfix'+1...] * noi matrix list `from' } else{ capture ml model d0 gllam_ll $HG_fixe $HG_eqs, /* */ noscvars waldtest(0) nopreserve missing collinear * disp "ml init from, `skip' `copy'" ml init `from', `skip' `copy' matrix `from' = ML_b } matrix `b' = `from' * matrix list `b' } global HG_coln: colnames(`b') global HG_cole: coleq(`b') * matrix list `b' matrix `V' = `b''*`b' estimates post `b' `V' matrix `b' = e(b) matrix makeCns `constra' qui `noi' disp in gr "Constraints:" qui `noi' matrix dispCns qui `noi' disp " " matcproc M_T M_a M_C matrix M_inshrt = `b'*M_T local n = colsof(M_inshrt) qui `noi' disp "estimating `n' parameters" local i = 1 local lst "`y'" gen byte __0066 = 1 while `i'< `n'{ local lst `lst' "eq`i'" local i = `i' + 1 } global HG_eqs matrix coleq M_inshrt = `lst' matrix colnames M_inshrt = __0066 *matrix list M_inshrt *matrix `b' = M_inshrt*M_T' + M_a *matrix list `b' if "`gateaux'"~=""{ local nf = `nf' - (`nb' - `n') matrix `from' = M_inshrt[1,1..`nf'] } else if "`from'"~=""&"`long'"~=""{ matrix `from' = M_inshrt } } } if "`from'"~=""{ local from "from(`from')" } /* deal with prior */ global HP_prior = 0 if "`prior'"!=""{ if "`noi'"!=""{ local loud "loud" } *disp in re "prior loud is `prior' `loud'" init_prior, `prior' `loud' global HP_sprd = 0 } * disp in re "HP_prior: " $HP_prior /* deal with adapt */ global HG_adapt=0 if "`adapt'"~=""{ if $HG_free==1|($HG_gauss==0&$HG_mult==0){ disp in re "adapt can only be used with ip(g) or ip(m) option" exit 198 } global HG_adapt = 1 } /* Nothing to estimate if init option and no fixed effects or level-1 variance */ if $HG_init==1&M_nffc[1,$HG_tpff]==0&$HG_lev1==0{ disp in re "Nothing to estimate" exit 498 } local fit = 0 if (M_nffc[1,$HG_tpff]>0|$HG_lev1>0)&$HG_init==0&"`from'"==""&$HG_comp==0{ /* initial values for fixed effects using Stata's commands */ local lnk $HG_link if "$HG_link"=="recip"{ local lnk pow -1 } qui `noi' disp " " qui `noi' disp in gr "Initial values for fixed effects" if $HG_const&$HG_init==0 { qui `noi' disp in gr "(Not applying constraints at this point)" } qui `noi' disp " " if ("$HG_famil"=="gauss")&("$HG_link"=="ident")& "`s'"==""{ tempvar yn if "`offset'"~=""{ quietly gen double `yn' = `y' - $HG_off } else{ gen `yn' = `y' } quietly `noi' reg `yn' `indep' [`ftype'=`wt'], `constant' matrix M_initr[1,1]=ln(_result(9)) local fit = 1 qui drop `yn' } else if ($HG_nolog+$HG_oth+$HG_mlog==1)&("$HG_famil"=="binom"|$HG_nolog==1|/* */ $HG_mlog==1)&$HG_exp==0{ local fit = 1 local mnd = 1 if "$HG_denom"~=""{ qui summ $HG_denom, meanonly local mnd = r(mean) } if `mnd'>1 { if $HG_mlog>0 { disp in re "can't have denominator > 1 for mlogit" exit 198 } if ($HG_nolog>0) { disp in re "can't have denominator > 1 for ordinal response" exit 198 } qui `noi' glm `y' `indep' [`ftype'=`wt'], link(`lnk') /* */ fam(binom `denom') `constant' `offset' `log' } else{ if "$HG_link"=="logit"{ qui `noi' logit `y' `indep' [`ftype'=`wt'], `constant' `offset' `log' } else if "$HG_link"=="probit"{ qui `noi' probit `y' `indep' [`ftype'=`wt'], `constant' `offset' `log' } else if "$HG_link"=="cll"{ qui `noi' cloglog `y' `indep' [`ftype'=`wt'], `constant' `offset' `log' } else if "$HG_link"=="ll"{ tempvar yn qui gen `yn' = 1-`y' qui `noi' cloglog `yn' `indep' [`ftype'=`wt'], `constant' `offset' `log' } else if $HG_mlog==1{ qui `noi' mlogit `y' `indep' [`ftype'=`wt'] if $HG_ind==1, `constant' `basecat' `log' } else if "$HG_linko"=="ologit"&"`thresh'"==""{ qui `noi' ologit `y' `indep' [`ftype'=`wt'], `offset' `log' } else if "$HG_linko"=="oprobit"&"`thresh'"==""{ qui `noi' oprobit `y' `indep' [`ftype'=`wt'], `offset' `log' } else if "$HG_linko"=="ocll"|"$HG_link"=="sprobit"|"$HG_linko"=="soprobit"|$HG_nolog>1|"`thresh'"~=""{ local fit = 0 } else{ qui `noi' glm `y' `indep' [`ftype'=`wt'], link(`lnk') /* */ fam(binom `denom') `constant' `offset' `log' local fit = 1 } } } else if ("$HG_famil"=="poiss")&("$HG_link"=="log"){ qui `noi' poisson `y' `indep' [`ftype'=`wt'], `constant' `offset' `log' local fit = 1 } else if ("$HG_famil"=="gamma"& M_nbrf[1,1]==1){ qui `noi' glm `y' `indep' [`ftype'=`wt'], link(`lnk')/* */ fam(gamma) `constant' `offset' `log' matrix M_initr[1,1]= -ln($S_E_dc) local fit = 1 } } if `fit' == 0 &("`from'"==""|$HG_init==1)& (M_nffc[1,$HG_tpff]>0|$HG_lev1>0) { /* fit level 1 model using gllamm */ /* preserve macros */ qui `noi' disp in green "(using gllamm for inital values)" local eqs "$HG_eqs" local tprf = $HG_tprf local tplv = $HG_tplv local tpi = $HG_tpi local const = $HG_const local link $HG_link local linko $HG_linko local lev1 = $HG_lev1 local ngeqs = $HG_ngeqs local sprior = $HP_prior tempvar keep quietly gen int `keep' = $HG_wt1 quietly replace $HG_wt1 = `wt' matrix `mnip' = M_nip matrix `mnbrf' = M_nbrf local adapt = $HG_adapt /* change global macros */ local frm if $HG_init==0{ global HG_const = 0 } global HG_ngeqs = 0 if $HG_init==0{ /* scale of probit usually not identified without random effects */ global HG_link local k: word count `link' local kk = 1 while `kk' <= `k'{ local ll: word `kk' of `link' if "`ll'" == "sprobit"{ global HG_link "$HG_link probit" } else{ global HG_link $HG_link `ll' } local kk = `kk' + 1 } *disp in re "links were `link' and are $HG_link" global HG_linko local k: word count `linko' local kk = 1 while `kk' <= `k'{ local ll: word `kk' of `linko' if "`ll'" == "soprobit"{ global HG_linko "$HG_linko oprobit" } else{ global HG_linko "$HG_linko `ll'" } local kk = `kk' + 1 } *disp in re "ordinal links were `linko' and are $HG_linko" if $HG_lev1>0{ /* check if any of the families is the normal or gamma*/ local found = 0 local k: word count $HG_famil local kk = 1 while `kk'<= `k'{ local ll: word `kk' of $HG_famil if "`ll'"=="gauss"| "`ll'"=="gamma"{ local found = 1 } local kk = `kk' + 1 } if `found'==0{ global HG_lev1 = 0 matrix M_nbrf = (0) } } } matrix M_nip=(1,1\1,1) if $HG_lev1>0{ global HG_eqs $HG_s1 global HG_tprf=1 global HG_tpi=1 global HP_prior=0 } else{ global HG_eqs global HG_tprf=0 global HG_tpi=1 } if "`from'"~=""{ local frm `from' } global HG_adapt = 0 /* fit model for initial values */ global HG_tplv=1 /* Level 1 model */ local opt if $HG_init{ /* apply constraints if init option is used */ if $HG_const==1{ tempname b V matrix `b' = M_initf, M_initr capture ml model d0 gllam_ll $HG_fixe $HG_eqs, /* */ noscvars waldtest(0) nopreserve missing collinear *matrix list `b' ml init `b', skip /* uses parameters needed for model above */ matrix `b' = ML_b /* unconstrained param., no random effects */ global HG_coln: colnames(`b') global HG_cole: coleq(`b') if "`from'"~=""& "`long'"~=""{ /* assumes copy option */ local nb = colsof(`b') local nf = colsof(`amat') if `nb'~=`nf'{ disp in re "from matrix must be dimension " `nb' exit 189 } disp in gr "(copy option assumed for from matrix)" matrix `b' = `amat' matrix coleq `b' = $HG_cole matrix coln `b' = $HG_coln } matrix `V' = `b''*`b' estimates post `b' `V' matrix `b' = e(b) matrix makeCns `constra' qui `noi' disp in gr "Constraints:" qui `noi' matrix dispCns qui `noi' disp " " matcproc M_T M_a M_C matrix M_inshrt = `b'*M_T /* constrained parameters, no r. effects */ local n = colsof(M_inshrt) qui `noi' disp "estimating `n' parameters" if "`from'"~=""&"`long'"==""{ local nf = colsof(`amat') if `n'~=`nf'{ disp in re "from matrix must be dimension " `n' exit 189 } matrix M_inshrt = `amat' } local i = 1 local lst "`y'" gen byte __0066 = 1 while `i'< `n'{ local lst `lst' "eq`i'" local i = `i' + 1 } global HG_eqs matrix coleq M_inshrt = `lst' matrix colnames M_inshrt = __0066 local frm "from(M_inshrt)" matrix coleq M_inshrt = `lst' matrix colnames M_inshrt = __0066 local n = colsof(M_inshrt) global HG_fixe (`y': `y' =__0066, nocons) local i = 1 while `i'< `n'{ global HG_fixe $HG_fixe (eq`i': __0066, nocons) local i = `i' + 1 } } local opt `options' } * noi qui `noi' hglm_ml `y', /* */ $HG_obs `log' title("fixed effects model") /* */ `frm' `trace' skip `difficult' `opt' `copy' `iterate' `eval' local fit = 1 if $HG_init==0 { quietly `noi' ml display, level(`level') nohead } if $HG_init{ if $HG_error==0{ if "`eval'"~=""{ local robust local cluster local pw } noi prepare, `robust' `cluster' `pw' `dots' `noi' * disp in re "running delmacs" delmacs restore estimates repost, esample(`touse') estimate local cmd "gllamm" estimate local predict "gllapred" * disp in re "running replay" noi Replay, level(`level') `eform' `allc' `robust' `cluster' `display' `eval' exit 0 } } if $HG_lev1>0{ local num=M_nbrf[1,1] matrix `mat'=e(b) if $HG_nats{ matrix `mat'=`mat'[1,"s1:"] } else{ matrix `mat'=`mat'[1,"lns1:"] } local i=1 while `i'<=`num'{ matrix M_initr[1,`i']=`mat'[1,`i'] local i=`i'+1 } } if $HG_init{ delmacs restore estimates repost, esample(`touse') estimate local cmd "gllamm" estimate local predict "gllapred" exit 0 } /* restore global macros */ global HG_tplv=`tplv' global HG_eqs "`eqs'" global HG_tprf=`tprf' global HP_prior=`sprior' global HG_tpi=`tpi' global HG_link "`link'" global HG_linko "`linko'" global HG_ngeqs = `ngeqs' quietly replace $HG_wt1=`keep' matrix M_nip=`mnip' matrix M_nbrf = `mnbrf' global HG_const = `const' global HG_lev1 = `lev1' global HG_adapt = `adapt' } if `fit'{ /* put estimates in `b' */ local cn: colnames(M_initf) local ce: coleq(M_initf) matrix M_initf=e(b) capture matrix colnames M_initf = `cn' capture matrix coleq M_initf = `ce' local num=M_nffc[1,$HG_tpff] if `num'>0 { local nn=colsof(M_initf) if `nn'<`num'{ disp in re "variables have been dropped, can't continue" exit 198 } matrix M_initf=M_initf[1,1..`num'] * matrix list M_initf } if $HG_const==1{ matrix `b' = nullmat(M_initf), nullmat(M_initr) matrix M_inshrt = `b'*M_T } if $HG_error==1{ exit } } /* estimation */ *qui `noi' disp in smcl in gr "{hline 78}" _n qui `noi' di in gr _dup(78) "-" _n qui `noi' dis qui `noi' dis "start running on $S_DATE at $S_TIME" * check this: *local skip if $HG_const==1{ matrix coleq M_inshrt = `lst' matrix colnames M_inshrt = __0066 local n = colsof(M_inshrt) global HG_fixe (`y': `y' =__0066, nocons) local i = 1 while `i'< `n'{ global HG_fixe $HG_fixe (eq`i': __0066, nocons) local i = `i' + 1 } } * disp "`trace' `options' " * disp "$HG_obs `log' `from'" * disp "`search' `lf0' `gateaux' `skip' `difficult' `eval' " /* means and sds for adapt */ if $HG_adapt{ local rf = 1 while `rf'<=$HG_tprf{ tempname junk global HG_MU`rf' "`junk'" tempname junk global HG_SD`rf' "`junk'" gen double ${HG_MU`rf'}=0 gen double ${HG_SD`rf'}=1 local rf2 = `rf' + 1 while `rf2' < $HG_tprf { tempname junk global HG_C`rf2'`rf' "`junk'" gen double ${HG_C`rf2'`rf'}=0 local rf2 = `rf2' + 1 } local rf = `rf' + 1 } } * noi capture noi capture noi hglm_ml `y', `trace' `options' /* */ $HG_obs `log' title("gllamm model") `from' `iterate'/* */ `search' `lf0' `gateaux' `skip' `copy' `difficult' `eval' if _rc>0{ global HG_error=1 } qui `noi' dis "finish running on $S_DATE at $S_TIME" qui `noi' dis " " if $HG_error==0{ if "`eval'"~=""{ /* do not want to compute robust standard errors */ local robust local cluster local pw } noi prepare, `robust' `cluster' `pw' `dots' `noi' * disp "running delmacs" delmacs * disp "restore" restore estimates repost, esample(`touse') estimate local cmd "gllamm" estimate local predict "gllapred" * disp "running replay" noi Replay, level(`level') `eform' `allc' `robust' `cluster' `display' `eval' } end program define prepare syntax [, ROBUST CLUSTER(varname) PWEIGHT DOTS NOISILY] * disp "options are: `robust' `cluster' `pweight' `dots' `noisily'" tempname b v X U matrix `b' = e(b) local n = colsof(`b') matrix M_Vs = e(V) capture matrix `v' = inv(M_Vs) if _rc==0{ matrix symeigen `X' `U' = `v' global HG_cn = sqrt(`U'[1,1]/`U'[1,`n']) } else{ global HG_cn = -1 } if $HG_const { matrix M_Vs = M_T*M_Vs*M_T' } /* deal with robust */ if "`robust'"~=""|"`cluster'"~=""|"`pweight'"~=""{ if "`cluster'"~=""{ global HG_rbcls "`cluster'" *disp "HG_rbcls is $HG_rbcls" local cluster cluster(`cluster') } * disp "calling gllarob" qui `noisily' gllarob, first macs `cluster' `dots' } * disp "HG_const = " $HG_const * disp "running remcor" qui remcor `b' if $HG_const { * disp "running procstr" qui procstr } end program define hglm_ml version 6.0 syntax varlist(min=1)[, TITLE(passthru) LF0(numlist) noLOg TRace /* */ OBS(passthru) FROM(string) SEarch(integer 0) Gateaux(numlist min=3 max=3) skip copy/* */ noDIFficult EVal ITerate(numlist) *] *disp in re "running hglm_ml" *disp in re "iterate: `iterate'" if "`log'"=="" { local log "noisily" } else{ local log } if "`trace'"~="" { local noi "noisily" } parse "`varlist'", parse(" ") local y "`1'" tempvar mysamp tempname b f V M_init M_initr a lnf mlnf ip deriv local adapt = $HG_adapt if "`from'"~=""{ matrix `M_init'=`from' if "`eval'"~=""|`adapt'==1{ *noi disp "$HG_fixe $HG_eqs" ml model d0 gllam_ll $HG_fixe $HG_eqs, /* */ noscvars waldtest(0) nopreserve missing collinear *disp "ml init M_init, `skip' `copy'" ml init `M_init', `skip' `copy' matrix ML_g=ML_b if `adapt'==0{ qui `noi' capture ml maximize, search(off) /* */ iterate(0) novce nooutput nowarn /* `options' */ disp in gr "log-likelihood = " in ye e(ll) exit 0 } else{ matrix `M_init' = ML_b global ML_y1 `y' noi gllam_ll 1 "`M_init'" "`lnf'" "junk" "junk" 1 } * matrix list `M_init' } if "`gateaux'"~=""&$HG_free==0{ disp in re "option gateaux not allowed (ignored) for fixed integration points" } else if "`gateaux'"~=""&$HG_free==1{ qui `noi' disp in gr "Gateaux derivative" if $HG_tplv>2{ disp "searching for additional point at level " $HG_tplv } local ll=$HG_tplv-1 local tprf=M_nrfc[2,$HG_tplv]-M_nrfc[2,`ll'] /* number of locations */ capture local mf = colsof(M_initf) if _rc>0 {local mf = 0} capture local mr = colsof(M_initr) if _rc>0 {local mr = 0} if $HG_const{ local nreq = colsof(M_inshrt) - `tprf' - 1 local cn: colnames(M_inshrt) local ce: coleq(M_inshrt) } else{ local nreq = `mf'+`mr'-`tprf'-1 } if `nreq'~=colsof(`M_init'){ disp in re "initial value vector should have length `nreq'" matrix list `from' global HG_error=1 exit 198 } *set trace on local l = `mr' - `tprf'-1 /* length of previous M_initr */ local locp = $HG_befB - `mf' local beg = `locp'-`tprf' /* position of first new mass-point parameter */ * matrix list M_initr matrix `a' = M_initr[1,`beg'..`locp'] /* new mass-point parameters */ * noi disp "new mass-point par" * matrix list `a' if $HG_cip==0{ /* new point is one before last since last has no prob. par. */ local locp = `locp' - `tprf' } local befB = $HG_befB - `tprf' - 1 /* befB for previous parameter vector */ global HG_befB if `nreq' > `befB' { matrix `M_init' = `M_init'[1,1..`befB'],`a',`M_init'[1,`befB'+1...] } else { matrix `M_init' = `M_init'[1,1..`befB'],`a' } * matrix list `M_init' local locp = `mf' + `locp' *disp "position of new p-parameter = " `locp' local nreq = `locp' - `tprf' -1 *disp "position of first new location parameter = " `nreq'+1 if $HG_cip==0{ * new point must be one before last since new probability tiny local jl = 1 while `jl'<=`tprf'{ matrix `M_init'[1,`locp'+`jl']=`M_init'[1,`nreq'+`jl'] local jl = `jl' + 1 } } *matrix list `M_init' tokenize "`gateaux'" local min = `1' local max = `2' local num = `3' local stp = (`max'-`min')/(`num'-1) matrix `M_init'[1,`locp']=-6 /* mass of new masspoint */ scalar `mlnf'=0 matrix `ip'=M_ip matrix `ip'[1,1]=1 *recursive loop matrix `ip'[1,`tprf']=1 local k = `nreq' + `tprf' matrix `M_init'[1,`k']=`min' local nxtrf = `tprf'+1 matrix `ip'[1,`nxtrf']=`num' local rf = `tprf' while `rf' <= `tprf'{ *reset ip up to random effect `rf' while (`rf'>1) { local rf = `rf'-1 matrix `ip'[1,`rf'] = 1 local k = `nreq' + `rf' matrix `M_init'[1,`k']=`min' } * update lowest digit local rf = 1 while `ip'[1,`rf'] <= `num'{ local k = `nreq' + `rf' matrix `M_init'[1,`k'] = `min' + (`ip'[1,`rf']-1)*`stp' * matrix list `M_init' global ML_y1 `y' gllam_ll 0 "`M_init'" "`lnf'" noi di in gr "." _c * noisily disp "likelihood=" `lnf' if (`lnf'>`mlnf'|`mlnf'==0)&`lnf'~=.{ scalar `mlnf'=`lnf' matrix M_initr=`M_init' } matrix `ip'[1,`rf'] = `ip'[1,`rf'] + 1 } matrix `ip'[1,`rf'] = `num' /* lowest digit has reached the top */ while `ip'[1,`rf']==`num'&`rf'<=`tprf'{ local rf = `rf' + 1 } * rf is first r.eff that is not complete or rf>nrf if `rf'<=`tprf'{ matrix `ip'[1,`rf'] = `ip'[1,`rf'] + 1 local k = `nreq' + `rf' matrix `M_init'[1,`k'] = `min' + (`ip'[1,`rf']-1)*`stp' } } if "`lf0'"~=""{ local junk: word 2 of `lf0' * disp in re "junk = " `junk' * disp in re "mlnf - lf0 is " `mlnf' " - " `junk' scalar `deriv' = `mlnf'-`junk' disp " " disp in ye "maximum gateaux derivative is " `deriv' * matrix list `M_initr' if `deriv'<0.00001{ disp in re "maximum gateaux derivative less than 0.00001" global HG_error=1 exit } } else{ disp in ye "no gateaux derivarives could be calculated without lf0() option" matrix list `M_initr' } matrix `M_init' = M_initr * starting log odds for new location matrix `M_init'[1,`locp']=-3 if $HG_const{ matrix colnames `M_init' = `cn' matrix coleq `M_init' = `ce' } * matrix list `M_init' } /* end if gateaux */ } /* end if from */ else{ /* no from() */ if "`gateaux'"~=""{ disp in red "gateaux can't be used without option from()" exit 198 } if "`eval'"~=""{ disp in red "eval option only allowed with from()" exit 198 } capture matrix `M_init'=M_initf if $HG_tprf|$HG_lev1>1{ matrix `M_initr'=M_initr local max=3 local min=0 scalar `mlnf' = 0 local f1= M_nbrf[1,1]+1 local l=colsof(M_initr) local m=1 if `search'>1{ if $HG_const==1{ disp in re "search option does not work yet with constraints" exit 198 } else{ qui `noi' disp in gr /* */ "searching for initial values for random effects" qui `noi' disp "likelihood:" } } while `m'<=`search'{ /* begin search */ * matrix list M_initr matrix `a'=`M_init',M_initr *matrix list `a' global ML_y1 `y' noisily gllam_ll 0 "`a'" "`lnf'" qui `noi' disp in ye %10.0g `lnf' " " _c if mod(`m',6)==0 {qui `noi' disp } *qui `noi' disp "likelihood=" `lnf' if (`lnf'>`mlnf'|`m'==1)&`lnf'~=. { scalar `mlnf'=`lnf' matrix `M_initr'=M_initr } local k=`f1' while `k'<=`l'{ matrix M_initr[1,`k']=`min' + (`max'-`min')*uniform() local k=`k'+1 } local m = `m' + 1 } /* end search */ matrix `M_init' = nullmat(`M_init'),`M_initr' } if $HG_const{ matrix `M_init' = M_inshrt } } if "`difficult'"~=""{ local difficu /* erase macro */ } else{ local difficu "difficult" /* default */ } * disp "$HG_fixe $HG_eqs, init(`M_init',`skip') " * disp "`lf0' `obs' `trace' `difficu' `options'" *matrix list `M_init' if "`lf0'"~="" { local lf0 "lf0(`lf0')" } * matrix list `M_init' if `adapt'{ tempname pa1 tempname pa2 tempname ad1 tempname ad2 ad0 * global HG_adapt=0 global ML_y1 `y' if "`from'"==""{ noi gllam_ll 1 "`M_init'" "`lnf'" "junk" "junk" 1 } scalar `ad2' = `lnf' scalar `ad1' = 0 local i = 2 qui `log' di in gr _n "Running adaptive quadrature" qui `noi' di in gr _dup(78) "-" _n "Iteration 0 of adaptive quadrature:" *qui `noi' di in smcl in gr "{hline 78}" _n "Iteration 0 of adaptive quadrature:" qui `noi' di in gr "Initial parameters:" _n qui `noi' mat list `M_init', noheader noblank format(%9.0g) /* first calculation of adaptive quadrature points ==> ad2 */ global HG_adapt=1 qui `noi' disp in gr _n "Updated log likelihood:" qui `noi' disp in ye %10.0g `ad2' " " _c while abs((`ad1'-`ad2')/`ad2')>1e-8&`i'<120&`ad2'~=.{ scalar `ad1' = `ad2' noi gllam_ll 1 "`M_init'" "`ad2'" "junk" "junk" 1 qui `noi' disp in ye %10.0g `ad2' " " _c if mod(`i',6)==0 {qui `noi' disp } local i = `i' + 1 } if `i'>=120{ disp in re "Convergence not achieved: try with more quadrature points" global HG_error=1 exit } if `ad2'==.{ disp in re "Log-likelihood cannot be computed" global HG_error=1 exit } qui `noi' di in gr _n _col(52) "log likelihood = " in ye %10.0g /* */ scalar(`ad2') if "`noi'"==""{ qui `log' di in gr "Iteration 0: log likelihood = " in ye %10.0g /* */ scalar(`ad2') } if "`eval'"~=""{ qui `noi' di in gr _dup(78) "-" *qui `noi' di in smcl in gr "{hline 78}" di in gr "log likelihood = " in ye %10.0g /* */ scalar(`ad2') _n qui `noi' capture ml maximize, search(off) /* */ iterate(0) novce nooutput nowarn /* `options' */ *delmacs *exit 1 exit 0 } if "`iterate'"~=""{ local iter=`iterate' } else{ local iter=150 } global HG_adapt=1 /* loop, update parameters then adaptive quadrature */ capture `log' ml model d0 gllam_ll $HG_fixe $HG_eqs, /* */ noscvars `lf0' `obs' `title' /* */ waldtest(0) nopreserve missing collinear if `iter'==0{ qui `noi' di in gr _dup(78) "-" disp in green "Adaptive quadrature has not converged" ml init `M_init', `skip' `copy' capture ml maximize, search(off) `difficu' /* */ `trace' iterate(0) nolog /* */ nooutput noclear /* `options' */ exit 0 } local it = 0 scalar `ad0' = `ad2'*1.2 scalar `pa2' = 0 local only1 iteronly1 local fst = 0 while abs((`ad0'-`ad2')/`ad0')>1e-6&`it'<`iter'{ scalar `ad0' = `ad2' local it = `it' + 1 qui `noi' di in gr _dup(78) "-" _n "Iteration `it' of adaptive quadrature:" *qui `noi' di in smcl in gr "{hline 78}" _n "Iteration `it' of adaptive quadrature:" qui `noi' disp in gr "Updated parameters:" _n /* update parameters ==> pa2 */ ml init `M_init', `skip' `copy' *ml report * noi capture noi capture ml maximize, search(off) `difficu' /* */ `trace' iterate(1) nolog /* */ nooutput noclear `only1' hessian gradient /* `options' */ * local rc = _rc if `rc' == 198&`fst'==0{ local only1 *noi capture noi capture ml maximize, search(off) `difficu' /* */ `trace' iterate(1) nolog /* */ nooutput noclear `only1' /* `options' technique(bfgs) */ local rc = _rc local fst = 1 } if `rc'>1 { di in red "(error occurred in ML computation)" di in red "(use trace option and check correctness " /* */ "of initial model)" global HG_error=1 exit `rc' } qui `noi' mat list $ML_b, noheader noblank format(%9.0g) qui `noi' di /* blank line */ scalar `pa1' = `pa2' scalar `pa2' = e(ll) matrix `M_init'=e(b) /* update adaptive quadrature ==> ad2 */ local j = 2 scalar `ad1' = `pa2' qui `noi' disp in gr "Updated log likelihood: " qui `noi' disp in ye %10.0g `ad1' " " _c while (abs((`ad1'-`ad2')/`ad2')>1e-8)&`j'<120&`ad2'~=.{ global ML_y1 `y' scalar `ad1' = `ad2' noi gllam_ll 1 "`M_init'" "`ad2'" "junk" "junk" 1 qui `noi' disp in ye %10.0g `ad2' " " _c if mod(`j',6)==0 { qui `noi' disp } local j = `j' + 1 } if `ad2'==.{ disp in re "Log-likelihood cannot be computed" global HG_error=1 exit } qui `noi' di in gr _n _col(52) "log likelihood = " in ye %10.0g /* */ scalar(`ad2') if "`noi'"==""{ qui `log' di in gr "Iteration " `it' ": log likelihood = " in ye %10.0g /* */ scalar(`ad2') } *qui `noi' disp in gr "log-likelihood is " in ye `ad2' in gre /* */ " and was " in ye `ad0' in gre ", relative change: " /* */ in ye abs((`ad2'-`ad0')/`ad0') } qui `noi' di in gr _dup(78) "-" _n *qui `noi' di in smcl in gr "{hline 78}" _n *tempname v *capture matrix `v'=e(V) *capture matrix `v' = inv(`v') *if _rc>0{ if `it'==`iter'&((`ad0'-`ad2')/`ad0')>1e-6 { disp in green "Adaptive quadrature did not converge in `iter' iterations" exit 0 } qui `log' disp in gr _n _n "Adaptive quadrature has converged, running Newton-Raphson" _c ml init `M_init', skip copy capture `log' ml maximize, search(off) `difficu' /* */ `trace' nooutput iterate(`iter') /* `options' */ local rc = _rc if `rc'>1 { di in red "(error occurred in ML computation)" di in red "(use trace option and check correctness " /* */ "of initial model)" global HG_error=1 exit `rc' } *} } else{ /* not adaptive */ if "`iterate'"~=""{ local iter "iterate(`iterate')" } else{ local iter } capture `log' ml model d0 gllam_ll $HG_fixe $HG_eqs, /* */ maximize search(off) /* */ init(`M_init', `skip' `copy') noscvars `lf0' `obs' `title' `trace' /* */ waldtest(0) nopreserve missing `difficu' collinear `iter'/* `options' */ * technique(bfgs) gtol(1e-4) } local rc = _rc if `rc'>1 { di in red "(error occurred in ML computation)" di in red "(use trace option and check correctness " /* */ "of initial model)" global HG_error=1 exit `rc' } if `rc'==1 { di in red /* */ "(Maximization aborted)" delmacs global HG_error=1 exit 1 } else if $HG_error==1{ disp in red "some error has occurred" exit } end program define lnkfm version 6.0 args link fam global S_1 /* link */ global S_2 /* family */ lnk "`1'" fm "`2'" if "$S_1" == "" { if "$S_2" == "gauss" { global S_1 "ident" } if "$S_2" == "poiss" { global S_1 "log" } if "$S_2" == "binom" { global S_1 "logit" } if "$S_2" == "gamma" { global S_1 "recip" } } /* if ("$S_1"=="mlogit"|"$S_1"=="smlogit")&"$S_2"~="binom"{ disp in red "mlogit link must be combined with binomial probability" exit 198 } */ if ("$S_1"=="mlogit"|"$S_1"=="smlogit"|"$S_1"=="ologit"|"$S_1"=="oprobit"|"$S_1"=="soprobit"|"$S_1"=="ocll"){ global S_2 } end program define fm version 6.0 args fam local f = lower(trim("`fam'")) local l = length("`f'") if "`f'" == substr("gaussian",1,max(`l',3)) { global S_2 "gauss" } else if "`f'" == substr("normal",1,max(`l',3)) { global S_2 "gauss" } else if "`f'" == substr("poisson",1,max(`l',3)) { global S_2 "poiss" } else if "`f'" == substr("binomial",1,max(`l',3)) { global S_2 "binom" } else if "`f'" == substr("gamma",1,max(`l',3)) { global S_2 "gamma" } else if "`f'" != "" { noi di in red "unknown family() `fam'" exit 198 } if "$S_2" == "" { global S_2 "gauss" } end program define lnk version 6.0 args link local f = lower(trim("`link'")) local l = length("`f'") if "`f'" == substr("identity",1,max(`l',2)) { global S_1 "ident" } else if "`f'" == substr("log",1,max(`l',3)) { global S_1 "log" } else if "`f'" == substr("logit",1,max(`l',4)) { global S_1 "logit" } else if "`f'" == substr("mlogit",1,max(`l',3)) { global S_1 "mlogit" } else if "`f'" == substr("smlogit",1,max(`l',3)) { global S_1 "smlogit" } else if "`f'" == substr("ologit",1,max(`l',3)) { global S_1 "ologit" } else if "`f'" == substr("oprobit",1,max(`l',3)) { global S_1 "oprobit" } else if "`f'" == substr("probit",1,max(`l',3)) { global S_1 "probit"} else if "`f'" == substr("ocll",1,max(`l',3)) { global S_1 "ocll"} else if "`f'" == substr("cll",1,max(`l',3)) { global S_1 "cll"} else if "`f'" == substr("ll",1,max(`l',2)) { global S_1 "ll"} else if "`f'" == substr("sprobit",1,max(`l',3)) { global S_1 "sprobit"} else if "`f'" == substr("soprobit",1,max(`l',3)) { global S_1 "soprobit"} else if "`f'"==substr("reciprocal",1,max(`l',3)) { global S_1 "recip" } else if "`f'" != "" { noi di in red "unknown link() `link'" exit 198 } end program define delmacs, eclass version 6.0 /* deletes all global macros and matrices and store some results in e()*/ tempname var if "$HG_tplv"==""{ * macros already gone exit } local nrfold = M_nrfc[2,1] local lev = 2 while (`lev'<=$HG_tplv){ local i2 = M_nrfc[2,`lev'] local i1 = `nrfold'+1 local i = `i1' local nrfold = M_nrfc[2,`lev'] local n = M_nrfc[1,`lev'] local n = M_nip[2,`n'] capture est matrix zps`n' M_zps`n' while `i' <= `i2'{ local n = M_nip[2,`i'] capture est matrix zlc`n' M_zlc`n' capture est matrix zps`n' M_zps`n' local i = `i' + 1 } local lev = `lev' + 1 } if $HG_free==0&$HG_init==0{ est matrix chol CHmat } if $HG_free{ est matrix mnp M_np } est matrix nrfc M_nrfc est matrix nffc M_nffc est matrix nbrf M_nbrf est matrix nu M_nu capture est matrix Vs M_Vs capture est matrix mresp M_resp capture est matrix mrespm M_respm capture est matrix frld M_frld if $HG_ngeqs>0{ est matrix mngeqs M_ngeqs } matrix drop M_ip est matrix nip M_nip capture est matrix mb M_b matrix drop M_znow capture matrix drop M_initf capture matrix drop M_initr capture matrix drop M_chol capture est matrix mb M_b est matrix olog M_olog capture est matrix moth M_oth if $HG_const == 1{ capture drop __0066 est matrix a M_a * est matrix C M_C est matrix T M_T est local coln $HG_coln est local cole $HG_cole global HG_coln global HG_cole } /* globals defined in gllam_ll */ local i=1 while (`i'<=$HG_tpff){ global HG_xb`i' local i= `i'+1 } local i = 1 while (`i'<=$HG_tprf){ global HG_s`i' local i= `i'+1 } local i = 1 while (`i'<=$HG_tplv){ global HG_wt`i' local i = `i' + 1 } if $HG_adapt{ macro drop HG_zuoff macro drop HG_SD* macro drop HG_MU* macro drop HG_E* macro drop HG_V* macro drop HG_C* } macro drop HG_zip macro drop HG_lzpr est local nats=$HG_nats global HG_nats est local noC=$HG_noC global HG_noC est local noC1=$HG_noC1 global HG_noC1 global HG_noC2 est local adapt=$HG_adapt global HG_adapt est local const = $HG_const global HG_const global HG_fixe est local inter = $HG_inter global HG_inter global HG_dots est local ngeqs = $HG_ngeqs global HG_ngeqs est local nolog = $HG_nolog if $HG_nolog>0{ est local k_eform = 1 /* controls how many eqns are exponentiated */ } else{ est local k_eform = $HG_tpff } global HG_nolog est local ethr = $HG_ethr global HG_ethr est local mlog = $HG_mlog global HG_mlog est local smlog = $HG_smlog global HG_smlog global HG_lvolo est local oth = $HG_oth global HG_oth est local lev1 = $HG_lev1 global HG_lev1 est local bmat = $HG_bmat global HG_bmat est local tplv = $HG_tplv global HG_tplv est local tprf = $HG_tprf global HG_tprf est local prior = $HP_prior if $HP_prior == 1{ local a_ll = `e(ll)' - $HP_res est local a_ll = `a_ll' est local lpr = $HP_res global HP_res if $HP_invga==1{ est local invga $HP_invga est local shape $shape est local rate $rate global HP_invga global shape global rate } if $HP_invwi==1{ est local invwi $HP_invwi est local df $df est matrix scale scale global HP_invwi global df } if $HP_foldt==1{ est local foldt $HP_foldt est local df $df est local scale $scale est local location $location global HP_foldt global df global scale global location } if $HP_logno==1{ est local logno $HP_logno est local meanlog $meanlong est local sdlog $sdlog global HP_logno global meanlog global sdlog } if $HP_gamma==1{ est local gamma $HP_gamma est local scale $HP_scale est local var $HP_var est local shape $HP_shape global HP_gamma global HP_scale global HP_var global HP_shape } if $HP_corre==1{ est local corre $HP_corre est local alpha $alpha est local beta $beta global HP_corre global alpha global beta } if $HP_boxco==1{ est local boxdo $HP_boxco est local scale $scale est local lambda $lambda global HP_boxco global scale global lambda } if $HP_spect==1{ est local spect $HP_spect est local alpha $alpha est local beta $beta global HP_spect global alpha global beta } if $HP_wisha==1{ est local wisha $HP_wisha est local df $df est matrix scale scale global HP_wisha global df } } global HP_prior est local tpi = $HG_tpi global HG_tpi est local tpff = $HG_tpff global HG_tpff est local clus "$HG_clus" global HG_clus est local weight "$HG_weigh" global HG_weigh est local pweight "$HG_pwt" global HG_pwt global which global HG_gauss est local free = $HG_free global HG_free est local mult = $HG_mult global HG_mult est local cip = $HG_cip est local famil "$HG_famil" global HG_famil est local link "$HG_link" global HG_link est local linko "$HG_linko" global HG_linko capture est local exp $HG_exp global HG_exp capture est local expf $HG_expf global HG_expf est local lv "$HG_lv" global HG_lv est local fv "$HG_fv" global HG_fv global HG_nump global HG_eqs global HG_obs est local offset "$HG_off" global HG_off est local denom "$HG_denom" global HG_denom est local cor = $HG_cor global HG_cor est local s1 "$HG_s1" global HG_s1 capture est local init $HG_init global HG_init capture est local ind "$HG_ind" global HG_ind capture est local comp $HG_comp global HG_comp capture est local coall "$HG_coall" global HG_coall capture est local cn = $HG_cn global HG_cn capture est local robclus "$HG_rbcls" global HG_rbcls global HG_befB global HG_cip end program define initmacs version 6.0 /* defines all global macros */ args nrf nip eqs geqs peqs s nats bmatrix touse dep frload tempname mat disp " " disp in gr "Random effects information for" in ye " $HG_tplv" in gr " level model" *disp in smcl in gr "{hline 78}" _n di in gr _dup(78) "-" _n /* deal with nrf */ matrix M_nrfc=J(2,$HG_tplv,1) if "`nrf'"==""|$HG_free|$HG_mult{ local k=1 while (`k'<=$HG_tplv){ matrix M_nrfc[1,`k']=`k' matrix M_nrfc[2,`k']=`k' local k=`k'+1 } } if "`nrf'"~=""{ local greater = 0 local k: word count `nrf' if `k'~=$HG_tplv-1 { if $HG_tplv==1{ disp in red "option nrf is meaningless for 1-level model" } else{ disp in red "option nrf() does not contain " $HG_tplv-1 " argument(s)" } exit 198 } parse "`nrf'", parse(" ") local k=2 while (`k'<=$HG_tplv){ matrix M_nrfc[2,`k']=`1' if `1'>1{ local greater=1 } local k=`k'+1 mac shift } /* make cumulative */ local k=2 while (`k'<=$HG_tplv){ matrix M_nrfc[2,`k']=M_nrfc[2,`k'-1]+M_nrfc[2,`k'] if $HG_free==0&$HG_mult==0{matrix M_nrfc[1,`k']=M_nrfc[2,`k']} local k=`k'+1 } if `greater'>0{ if "`eqs'"==""{ disp in re "eqs() option required" exit 198 } } } * matrix list M_nrfc global HG_tprf=M_nrfc[2,$HG_tplv] /* number of random effects */ global HG_tpi=M_nrfc[1,$HG_tplv] /* number of integration loops + 1 */ if $HG_tplv==$HG_tprf{ if $HG_cor==0{ disp "option nocorrel ignored because no multiple r. effects per level" } } /* deal with nip and set up zloc and zps */ if "`nip'"~=""{ local k: word count `nip' if `k'==1{ if `nip' == 1{ matrix M_nip=J(2,$HG_tprf,1) if $HG_free { if $HG_cip==1 { global HG_init=1 } } else{ disp in re "(co)variances of latent variables will not be identified with nip(1)" ghquad 1 } } } if `k'>1 | $HG_tpi == 2 { if `k'~=$HG_tpi-1{ disp in red "option nip() has `k' arguments, need 1 or " $HG_tpi-1 exit 198 } matrix M_nip=J(2,$HG_tprf,1) local i=1 while `i'<$HG_tpi{ local k: word `i' of `nip' local l = `i' + 1 matrix M_nip[1,`l']= `k' if $HG_mult{ * disp in re "k = `k'" local n = M_nrfc[2,`l'] - M_nrfc[2,`l'-1] if `n'>1 { mint `n' `k' matrix M_nip[1,`l'] = $S_1 } else{ local kk = (`k' + 1)/2 if int(`kk')~=`kk'{ disp in re "arguments in nip() option must be odd with ip(m)" exit 198 } ghquad `kk' matrix M_nip[1,`l']= `kk' } } else if $HG_free==0{ if $HG_gauss{ ghquad `k' } else{ lebesque `k' } } local i = `i' + 1 } /* end i */ } else{ /* one argument given for more than one loop */ matrix M_nip=J(2,$HG_tprf,`nip') matrix M_nip[1,1] = 1 if $HG_mult{ local l = 2 while `l' <= $HG_tpi{ local n = M_nrfc[2,`l'] - M_nrfc[2,`l'-1] if `n'>1{ mint `n' `nip' matrix M_nip[1,`l'] = $S_1 } else{ local kk = (`nip' + 1)/2 if int(`kk')~=`kk'{ disp in re "arguments in nip() option must be odd with ip(m)" exit 198 } ghquad `kk' matrix M_nip[1,`l']= `kk' } local l = `l' + 1 } } else if $HG_free==0{ if $HG_gauss{ ghquad `nip' } else{ lebesque `nip' } } } } else{ /* no nip argument given */ matrix M_nip=J(2,$HG_tprf,8) matrix M_nip[1,1] = 1 if $HG_mult{ local l = 2 while `l' <= $HG_tpi{ local n = M_nrfc[2,`l'] - M_nrfc[2,`l'-1] if `n'>1{ sphern `n' matrix M_nip[1,`l'] = $S_1 } else{ matrix M_nip[1,`l'] = 8 ghquad 8 } local l = `l' + 1 } } else if $HG_free==0{ if $HG_gauss{ ghquad 8 } else{ lebesque 8 } } } local i = M_nrfc[2,1]+1 while `i'<= $HG_tprf{ if $HG_free{ matrix M_nip[2,`i'] = `i' } else{ matrix M_nip[2,`i'] = M_nip[1,`i'] } local i = `i' + 1 } * noi matrix list M_nip capture matrix drop M_initr *matrix M_zlc8 = J(1,8,1) *matrix M_zps8 = J(1,8,ln(1/8)) /* deal with Eqs */ matrix M_frld = J(1,$HG_tprf,0) local depv `dep' matrix M_nbrf=(0) global HG_eqs global HG_nats = 0 local lns1 lns1 local log log if $HG_lev1>0{ if "`nats'"~=""{ global HG_nats = 1 local lns1 s1 local log } disp in gr "***level 1 equation:" if "`s'"~=""{ eq ? "`s'" local vars "$S_1" markout `touse' `vars' global HG_eqs "$HG_eqs (`lns1': `depv' `vars',nocons)" global HG_s1 "(`lns1': `depv' `vars',nocons)" } else{ local vars "_cons" global HG_eqs "$HG_eqs (`lns1': `depv')" global HG_s1 "(`lns1': `depv')" } local depv disp " " if $HG_lev1==1{disp in gr " `log' standard deviation"} else if $HG_lev1==2{disp in gr " `log' coefficient of variation"} else if $HG_lev1==3{disp in gr " `log'(sqrt(phi))"} disp in ye " `lns1': `vars'" local num: word count `vars' matrix M_nbrf=(`num') matrix `mat'=J(1,`num',-1) matrix colnames `mat'=`vars' matrix coleq `mat'=`lns1' matrix M_initr=nullmat(M_initr),`mat' } else{ matrix M_nbrf=(0) if "`s'"~=""{ disp in re "s() option ignored because families do not include dispersion parameters" } if "`nats'"~=""{ disp in re "nats option ignored because families do not include dispersion parameters" } } matrix M_np = J(1,$HG_tpi,0) if "`peqs'" ~= ""{ local k: word count `peqs' if $HG_free == 0 { disp in re "peqs() argument meaningless for quadrature estimation, ignored" } else if `k'~=$HG_tpi-1{ disp in red `k' " equations specified for prior probabilities: `peqs', need" $HG_tpi - 1 exit 198 } global HG_noC1 = $HG_noC2 } if "`eqs'"~=""{ local k: word count `eqs' if `k'~=$HG_tprf-1{ disp in red `k' " equations specified: `eqs', need " $HG_tprf-1 exit 198 } * check that they are equations and find number of variables in each: nbrf local lev=2 local l=1 local ic=0 while (`lev'<=$HG_tplv){ disp " " local m=$HG_tplv-`lev'+1 local clusnam: word `m' of $HG_clus disp " " disp in gr "***level `lev' (" in ye "`clusnam'" in gr ") equation(s):" if $HG_init==1{ di in gr "(init option: parameters will not be estimated)" di " " } local clusnam=substr("`clusnam'",1,3) local lv = `lev'-1 local clusnam "`clusnam'`lv'_" local i1=M_nrfc[2, `lev'-1] local j1=M_nrfc[2, `lev'] local nrf=`j1'-`i1' disp " (`nrf' random effect(s))" disp " " local rfl = 1 /* MASS POINTS */ if $HG_free { tempname pmat if $HG_cor==0{ disp "option nocorrel irrelevant for free masses" } local k = 1 local nloc = M_nip[1, `lev'] if $HG_cip{ local nloc = `nloc' - 1} * new local pnum = 0 local pvars _cons matrix `pmat' = J(1,1,0) if "`peqs'"~=""{ local levm = `lev' - 1 local peqn: word `levm' of `peqs' eq ? "`peqn'" local pvars "$S_1" markout `touse' `pvars' local pnum: word count `pvars' matrix M_np[1,`lev'] = `pnum' + 1 matrix `pmat' = J(1, `pnum', 0) } matrix colnames `pmat' = `pvars' * end new while `k' <= `nloc'{ disp " " disp in gre "class `k'" local j = `i1' while `j'< `j1'{ local eqnam: word `j' of `eqs' eq ? "`eqnam'" local vars "$S_1" markout `touse' `vars' local num: word count `vars' local ex = 0 if "`frload'"~=""{ local nex: word count `frload' local iex = 1 while `iex'<=`nex'{ local aex: word `iex' of `frload' if `aex' == `j'{ local ex = 1 } local iex = `iex' + 1 } if `ex'{ matrix M_frld[1,`j'+1] = 1 local num = `num' + 1 } } matrix `mat'=(`num') matrix M_nbrf=M_nbrf,`mat' if (`num'>1){ parse "`vars'", parse(" ") local vars1 "`1'" if `k'==1{ if `ex'==0{ mac shift } local vars2 "`*'" local eqnaml "`clusnam'`rfl'l" eq "`eqnaml': `vars2'" eq ? "`eqnaml'" disp " " disp in gr " lambdas for random effect " in ye `j' disp in ye " `eqnaml': `vars2'" global HG_eqs "$HG_eqs (`eqnaml': `depv' `vars2', nocons)" local depv local num=`num'-1 * initial loading on masspoints local lod = 1.1 + (`j'-1)*(-1)^(`j')/5 matrix `mat'=J(1,`num',`lod') matrix colnames `mat'= `vars2' matrix coleq `mat'=`eqnaml' matrix M_initr = nullmat(M_initr), `mat' } } else{local vars1 `vars'} disp " " disp in gr " location for random effect " in ye `j' local eqnam "z`lev'_`j'_`k'" if `nrf'==1{ local eqnam "z`lev'_1_`k'" } eq "`eqnam'":`vars1' eq ? "`eqnam'" disp in ye " `eqnam': `vars1'" global HG_eqs "$HG_eqs (`eqnam': `depv' `vars1', nocons)" local depv markout `touse' `vars1' * initial locations of mass points *local val = int((`k'+1)/2)*(-1)^`k'/10 local val = int((`k'+1)/2)*(-1)^`k' matrix `mat'=(`val') matrix colnames `mat'=`vars1' matrix coleq `mat'=`eqnam' matrix M_initr=nullmat(M_initr),`mat' local j = `j' + 1 local rfl = `rfl' + 1 } if `k'< M_nip[1, `lev']{ local eqnam "p`lev'_`k'" eq "`eqnam'": eq ? "`eqnam'" disp " " disp in gr " log odds for level " in ye `lev' if `pnum'>0{ disp in ye " `eqnam': `pvars' _cons" global HG_eqs "$HG_eqs (`eqnam':`pvars')" matrix coleq `pmat'=`eqnam' matrix M_initr=nullmat(M_initr),`pmat' } else{ disp in ye " `eqnam': _cons" global HG_eqs "$HG_eqs (`eqnam': `depv')" local depv * initial log odds for masspoints } * set constant matrix `mat'=(-.4) matrix colnames `mat'=_cons matrix coleq `mat'=`eqnam' matrix M_initr=nullmat(M_initr),`mat' } local k = `k' + 1 } * matrix list M_initr * disp "$HG_eqs" } /* STD DEVS */ else{ local j = `i1' while (`j'<`j1'){ local eqnam: word `l' of `eqs' eq ? "`eqnam'" local vars "$S_1" local num: word count `vars' local ex = 0 if "`frload'"~=""{ local nex: word count `frload' local iex = 1 while `iex'<=`nex'{ local aex: word `iex' of `frload' if `aex' == `j'{ local ex = 1 } local iex = `iex' + 1 } if `ex'{ matrix M_frld[1,`j'+1] = 1 local num = `num' + 1 } } matrix `mat'=(`num') matrix M_nbrf=M_nbrf,`mat' markout `touse' `vars' if "`vars'"==""{ local vars "_cons"} if `num'>1{ * vars1 is variable of first loading (fix at one) parse "`vars'", parse(" ") local vars1 "`1'" if `ex'==0{ mac shift } local vars "`*'" local eqnaml "`clusnam'`rfl'l" eq "`eqnaml'": `vars' eq ? "`eqnaml'" disp " " disp in gr " lambdas for random effect " in ye `j' disp in ye " `eqnaml': `vars'" global HG_eqs "$HG_eqs (`eqnaml': `depv' `vars', nocons)" local depv * initial values of loadings local lod = 1.1 + (`j'-1)*(-1)^(`j')/5 /*different loading for diff r.eff*/ matrix `mat'=J(1,`num'-1,`lod') matrix colnames `mat'=`vars' matrix coleq `mat'=`eqnaml' matrix M_initr=nullmat(M_initr),`mat' } else{ local vars1 `vars' } * variance local eqnam "`clusnam'`rfl'" eq "`eqnam'": `vars1' if `nrf'==1|$HG_cor==0{ disp in gr " standard deviation for random effect " in ye `j' } else{ disp " " disp in gr /* */" diagonal element of cholesky decomp. of covariance matrix" } disp in ye " `eqnam' : `vars1'" global HG_eqs "$HG_eqs (`eqnam': `depv' `vars1', nocons)" local depv * initial value of standard deviation matrix `mat' = (0.5) matrix colnames `mat' = `vars1' matrix coleq `mat' = `eqnam' matrix M_initr=nullmat(M_initr),`mat' local l=`l'+1 local j=`j'+1 local rfl = `rfl' + 1 } if `nrf' > 1&$HG_cor==1{ /* generate equations for covariance parameters */ disp " " disp in gr " off-diagonal elements" local ii=2 *local num = $HG_tplv-`lev'+1 *local eqnam: word `num' of $HG_clus *local eqnam = substr("`eqnam'",1,4) while (`ii'<=`nrf'){ local jj=1 while (`jj'<`ii'){ local eqnaml "`clusnam'`ii'_`jj'" eq "`eqnaml'": eq ? "`eqnaml'" disp in ye " `eqnaml': _cons" global HG_eqs "$HG_eqs (`eqnaml':)" matrix `mat'=(0) matrix colnames `mat'=_cons matrix coleq `mat'=`eqnaml' matrix M_initr=nullmat(M_initr),`mat' local jj = `jj' + 1 } local ii=`ii'+1 } } } /* end else $HG_free */ local lev=`lev'+1 } /* lev loop */ } /* endif equ given */ else{ /* random intercepts */ if M_nrfc[1,$HG_tplv]~=$HG_tplv{ "must specify equations for random effects" exit 198 } local k=$HG_tprf-1 matrix `mat'=J(1,`k',1) matrix M_nbrf=M_nbrf,`mat' local lev=2 disp " " while (`lev'<=$HG_tplv){ local l=$HG_tplv-`lev'+1 local clusnam: word `l' of $HG_clus disp " " disp in gr "***level `lev' (" in ye "`clusnam'" in gr ") equation(s):" if $HG_init==1{ di in gr "(init option: parameters will not be estimated)" di " " } local clusnam = substr("`clusnam'",1,4) local lv = `lev' - 1 local clusnam "`clusnam'`lv'" /*MASS POINTS */ if ($HG_free){ local k = 1 local nloc = M_nip[1, `lev'] if $HG_cip{ local nloc = `nloc' - 1} * new local pnum = 0 local pvars _cons tempname pmat matrix `pmat' = J(1,1,0) if "`peqs'"~=""{ local levm = `lev' - 1 local peqn: word `levm' of `peqs' eq ? "`peqn'" local pvars "$S_1" markout `touse' `pvars' local pnum: word count `pvars' matrix M_np[1,`lev'] = `pnum' + 1 matrix `pmat' = J(1, `pnum', 0) } matrix colnames `pmat' = `pvars' * end new while `k' <= `nloc'{ disp " " disp in gre "class `k'" local j = 1 local eqnam "z`lev'_1_`k'" disp in gr " location for random effect" disp in ye " `eqnam': _cons" global HG_eqs "$HG_eqs (`eqnam': `depv')" local depv * initial locations of mass points *local val = int((`k'+1)/2)*(-1)^`k'/10 local val = int((`k'+1)/2)*(-1)^`k' matrix `mat'=(`val') matrix colnames `mat'=_cons matrix coleq `mat'=`eqnam' matrix M_initr=nullmat(M_initr),`mat' if `k'0{ disp in ye " `eqnam': `pvars' _cons" global HG_eqs "$HG_eqs (`eqnam':`pvars')" matrix coleq `pmat'=`eqnam' matrix M_initr=nullmat(M_initr),`pmat' } else{ disp in ye " `eqnam': _cons" global HG_eqs "$HG_eqs (`eqnam':)" * initial log odds for masspoints } * set constant matrix `mat'=(-.4) matrix colnames `mat'=_cons matrix coleq `mat'=`eqnam' matrix M_initr=nullmat(M_initr),`mat' } local k = `k' + 1 } } /* ST. DEVS */ else{ local eqnam "`clusnam'" disp " " disp in gr " standard deviation of random effect" disp in ye " `eqnam': _cons" global HG_eqs "$HG_eqs (`eqnam':`depv')" local depv * initial value for sd matrix `mat'=(0.5) matrix colnames `mat'=_cons matrix coleq `mat'=`eqnam' matrix M_initr=nullmat(M_initr),`mat' local cons `cons'1 } local lev=`lev'+1 } } disp " " matrix `mat' = nullmat(M_initr), nullmat(M_initf) /* M_initr may not exist */ global HG_befB = colsof(`mat') /* deal with Bmatrix */ global HG_bmat = 0 if "`bmatrix'"~=""{ if $HG_tprf<2{ disp in re "bmatrix can only be used for more than 1 random effect" exit 198 } * capture matrix list `bmatrix' if _rc>0{ disp in re "bmatrix is not a matrix" exit 198 } local bn = colsof(`bmatrix') if rowsof(`bmatrix')~=`bn'{ disp in re "bmatrix must be square" exit 198 } if `bn'~=$HG_tprf-1{ disp in re "number of rows and columns of B matrix must be " $HG_tprf-1 exit 198 } matrix M_b=`bmatrix' global HG_bmat = 1 disp in gr "B-matrix:" local i = 1 while `i' <= `bn'{ local j = 1 while `j'<= `bn'{ if M_b[`i',`j']>0{ local eqnam b`i'_`j' disp " " disp in ye " `eqnam': _cons" global HG_eqs "$HG_eqs (`eqnam':)" * initial value for sd matrix `mat'=(0.5) matrix colnames `mat'=_cons matrix coleq `mat'=`eqnam' matrix M_initr=nullmat(M_initr),`mat' local cons `cons'1 } local j = `j' + 1 } local i = `i' + 1 } disp " " } /* total number of fixed linear predictors */ global HG_tpff = colsof(M_nffc) /* deal with geqs */ global HG_ngeqs = 0 if "`geqs'"~=""{ * M_ngeqs: first row says which random effect, second how many terms local num: word count `geqs' global HG_ngeqs = `num' matrix M_ngeqs=J(3,`num',0) local nxt = M_nffc[1,$HG_tpff] + colsof(M_initr) + 1 disp in gr "Regressions of random effects on covariates:" tokenize `geqs' local i = 1 while "`1'"~="" { local k = substr("`1'",2,1) local k = `k' + 1 if `k'>$HG_tprf { disp in red "eq `1' refers to a random effects that does not exist" exit 198 } local j = 1 while `j'<=`i'{ if M_ngeqs[1,`j']==`k' { disp in red "more than one geq given for random effect" `k'-1 exit 198 } local j = `j' + 1 } eq ? "`1'" local vars "$S_1" local num: word count `vars' matrix `mat'=J(1,`num',0) matrix colnames `mat'=`vars' matrix coleq `mat'=`1' matrix M_initr=nullmat(M_initr),`mat' markout `touse' `vars' disp in gr " equation for random effect " in ye `k'-1 disp in ye " `1': `vars'" global HG_eqs "$HG_eqs (`1': `vars', nocons)" matrix M_ngeqs[1,`i']=`k' matrix M_ngeqs[2,`i']=`num' matrix M_ngeqs[3,`i']=`nxt' local nxt = `nxt' + `num' local i = `i' + 1 mac shift } disp " " } global which = 15 /* the "clock" ip and znow*/ local k = $HG_tprf+2 matrix M_ip = J(1,`k',1) local k = $HG_tprf - 1 matrix M_znow =J(1,`k',1) end program define mint args n k if `k' == 5 { multn `n' } else if `k' == 7 { sphern `n' } else if `n'==2&(`k' == 9|`k' == 11|`k' == 15){ int2k `k' } else{ disp in re "nip must be 5 or 7 if dim >2 and 5, 7, 9, 11 or 15 if dim=2" exit 198 } end program define sphern version 6.0 /* degree 7 rule for n>=3 used by Naylor and Smith */ args n if `n' == 2 { spher2 } else if `n' == 3{ spher3 } else if `n' == 6{ spher6 } /* else if `n' == 4{ spher4 } */ else{ tempname r s t B C D r1 r2 A1 A2 x X y a norm im local num =2*( 2^`n' + 2*`n'^2) * generate U_n: 7-1 page 295 scalar `r' = sqrt(2) scalar `s' = sqrt(2/`n') scalar `t' = 1 scalar `B' = ln(8-`n') - ln(`n') - ln(`n'+2) - ln(`n'+4) + ln(2) - lngamma(`n'/2) scalar `C' = -`n'*ln(2) + 3*ln(`n') - ln(`n') - ln(`n'+2) - ln(`n'+4) + ln(2) - lngamma(`n'/2) scalar `D' = ln(4) - ln(`n') - ln(`n'+2) - ln(`n'+4) + ln(2) - lngamma(`n'/2) scalar `r2' = sqrt( (`n' + 2 - sqrt(2*(`n'+2)))/2 ) scalar `r1' = sqrt( (`n' + 2 + sqrt(2*(`n'+2)))/2 ) scalar `A2' = ln( ( `n' + 2 + sqrt(2*(`n'+2)) )/(4*(`n'+2)) )+ lngamma(`n'/2) scalar `A1' = ln( ( `n' + 2 - sqrt(2*(`n'+2)) )/(4*(`n'+2)) )+ lngamma(`n'/2) * Do r and B matrix `x' = J(1,`n',0) local ir = 1 while `ir' <= `n'{ local j = 1 while `j' <= 2 { matrix `x'[1,`ir'] = (-1)^`j'*`r' matrix `X' = nullmat(`X')\ `x' matrix `y' = nullmat(`y'),`B' matrix `x'[1,`ir'] = 0 local j = `j' + 1 } local ir = `ir' + 1 } * Do t and D matrix `x' = J(1,`n',0) local it1 = 1 while `it1' <= `n'{ local j1 = 1 while `j1' <= 2{ matrix `x'[1,`it1'] = (-1)^`j1'*`t' local it2 = `it1' + 1 while `it2' <= `n'{ local j2 = 1 while `j2' <= 2{ * disp "it1 = `it1' and it2 = `it2'" matrix `x'[1,`it2'] = (-1)^`j2'*`t' matrix `X' = `X'\ `x' matrix `y' = `y', `D' matrix `x'[1,`it2'] = 0 local j2 = `j2' + 1 } local it2 = `it2' + 1 } matrix `x'[1,`it1'] = 0 local j1 = `j1' + 1 } local it1 = `it1' + 1 } * Do s and C matrix `x' = J(1,`n',`s') matrix `im' = J(1,`n'+1,0) matrix `im'[1,1] = 1 local pm = `n'+1 while `pm' <= `n'+1{ /* set previous digits to 0 */ while `pm'>1{ matrix `x'[1,`pm'-1] = `s' matrix `im'[1,`pm'] = 0 local pm = `pm' - 1 } matrix `X' = `X'\ `x' matrix `y' = `y', `C' * matrix list `im' local pm = 1 while `im'[1,`pm'] == 1{ local pm = `pm' + 1 } * disp "pm = " `pm' /* pm is first incomplete digit */ if `pm'<= `n' + 1{ matrix `im'[1,`pm'] = 1 matrix `x'[1,`pm'-1] = -`s' local pm = `pm' - 1 } } *matrix list `X' *matrix list `y' matrix M_zlc`num' = `X'*`r1'\ `X'*`r2' matrix M_zlc`num' = M_zlc`num'' matrix `a'= J(1,`num'/2,`A1'),J(1,`num'/2,`A2') matrix M_zps`num' = `y',`y' matrix M_zps`num' = M_zps`num'+`a' global S_1 = `num' } /* end else */ end program define spher2 version 6.0 *p. 324, 7-1 tempname r s t A B C x capture matrix drop M_zlc12 capture matrix drop M_zps12 scalar `r' = sqrt(6) scalar `s' = sqrt((9-3*sqrt(5))/4) scalar `t' = sqrt((9+3*sqrt(5))/4) scalar `A' = ln(1/36) scalar `B' = ln(5 + 2*sqrt(5)) - ln(45) scalar `C' = ln(5 - 2*sqrt(5)) - ln(45) *r matrix `x' = J(1,2,0) local ir = 1 while `ir' <= 2{ local j = 1 while `j' <= 2{ matrix `x'[1,`ir'] = (-1)^`j'*`r' matrix M_zlc12 = nullmat(M_zlc12)\ `x' matrix M_zps12 = nullmat(M_zps12), `A' matrix `x'[1,`ir'] = 0 local j = `j' + 1 } local ir = `ir' + 1 } *t matrix `x' = J(1,2,`t') matrix M_zlc12 = M_zlc12\ `x' matrix M_zps12 = M_zps12, `C' local m = 1 /* number of - */ while `m' <= 2{ local i = 1 /* pos of first - */ while `i'<= 2{ matrix `x'[1,`i'] = -`t' local j = `i' + 1 if `m'>1{ while `j' <= 2 { matrix `x'[1,`j'] = -`t' matrix M_zlc12 = M_zlc12\ `x' matrix M_zps12 = M_zps12, `C' matrix `x'[1,`j'] = `t' local j = `j' + 1 } } else{ matrix M_zlc12 = M_zlc12\ `x' matrix M_zps12 = M_zps12, `C' } matrix `x'[1, `i'] = `t' local i = `i' + 1 } local m = `m' + 1 } * s matrix `x' = J(1,2,`s') matrix M_zps12 = M_zps12, `B' matrix M_zlc12 = M_zlc12\ `x' local m = 1 /* number of - */ while `m' <= 2{ local i = 1 /* pos of first - */ while `i'<= 2{ matrix `x'[1,`i'] = -`s' local j = `i' + 1 if `m'>1{ while `j' <= 2 { matrix `x'[1,`j'] = -`s' matrix M_zlc12 = M_zlc12\ `x' matrix M_zps12 = M_zps12, `B' matrix `x'[1,`j'] = `s' local j = `j' + 1 } } else{ matrix M_zlc12 = M_zlc12\ `x' matrix M_zps12 = M_zps12, `B' } matrix `x'[1, `i'] = `s' local i = `i' + 1 } local m = `m' + 1 } matrix M_zlc12 = M_zlc12' global S_1 = 12 end program define spher3 version 6.0 * Stroud (1971) p. 327, E_3^{r^2} 7-1 args n tempname r s t A B C D x im local num = 27 capture matrix drop M_zlc27 capture matrix drop M_zps27 scalar `r' = sqrt((15 + sqrt(15))/2) scalar `s' = sqrt(6 - sqrt(15)) scalar `t' = sqrt(9 + 2*sqrt(15)) scalar `A' = ln(720 + 8*sqrt(15)) - ln(2205) scalar `B' = ln(270 - 46*sqrt(15)) - ln(15435) scalar `C' = ln(162 + 41*sqrt(15)) - ln(6174) scalar `D' = ln(783 - 202*sqrt(15)) - ln(24696) local n = 3 * D matrix `x' = J(1,`n',`t') matrix `im' = J(1,`n'+1,0) matrix `im'[1,1] = 1 local pm = `n'+1 while `pm' <= `n'+1{ /* set previous digits to 0 */ while `pm'>1{ matrix `x'[1,`pm'-1] = `t' matrix `im'[1,`pm'] = 0 local pm = `pm' - 1 } matrix M_zlc`num' = nullmat(M_zlc`num')\ `x' matrix M_zps`num' = nullmat(M_zps`num'), `D' * matrix list `im' local pm = 1 while `im'[1,`pm'] == 1{ local pm = `pm' + 1 } * disp "pm = " `pm' /* pm is first incomplete digit */ if `pm'<= `n' + 1{ matrix `im'[1,`pm'] = 1 matrix `x'[1,`pm'-1] = -`t' local pm = `pm' - 1 } } * B matrix `x' = J(1,`n',0) local ir = 1 while `ir' <= `n'{ local j = 1 while `j' <= 2{ matrix `x'[1,`ir'] = (-1)^`j'*`r' matrix M_zlc`num' = M_zlc`num'\ `x' matrix M_zps`num' = M_zps`num',`B' matrix `x'[1,`ir'] = 0 local j = `j' + 1 } local ir = `ir' + 1 } * C matrix `x' = J(1,`n',0) local is1 = 1 while `is1' <= `n'{ local j1 = 1 while `j1' <= 2{ matrix `x'[1,`is1'] = (-1)^`j1'*`s' local is2 = `is1' + 1 while `is2' <= `n'{ local j2 = 1 while `j2' <= 2{ * disp "is1 = `is1' and is2 = `is2'" matrix `x'[1,`is2'] = (-1)^`j2'*`s' matrix M_zlc`num' = M_zlc`num'\ `x' matrix M_zps`num' = M_zps`num',`C' matrix `x'[1,`is2'] = 0 local j2 = `j2' + 1 } local is2 = `is2' + 1 } matrix `x'[1,`is1'] = 0 local j1 = `j1' + 1 } local is1 = `is1' + 1 } * A matrix M_zps`num' = M_zps`num',`A' matrix `x' = J(1,`n',0) matrix M_zlc`num' = M_zlc`num'\ `x' matrix M_zlc`num' = M_zlc`num'' global S_1 = `num' end program define spher4 version 6.0 * Stroud (1971) p. 329, E_4^{r^2} 7-1 args n tempname r s t A B C D x im local num = 49 capture matrix drop M_zlc`num' capture matrix drop M_zps`num' scalar `s' = sqrt(3 - sqrt(3)) scalar `r' = 2*`s' scalar `t' = sqrt(6 + 2*sqrt(3)) scalar `A' = -ln(4) scalar `B' = ln(9 + 5*sqrt(3)) - ln(576) scalar `C' = ln(9 + -5*sqrt(3)) - ln(576) local n = 4 * C and t matrix `x' = J(1,`n',0) local is1 = 1 while `is1' <= `n'{ local j1 = 1 while `j1' <= 2{ matrix `x'[1,`is1'] = (-1)^`j1'*`t' local is2 = `is1' + 1 while `is2' <= `n'{ local j2 = 1 while `j2' <= 2{ * disp "is1 = `is1' and is2 = `is2'" matrix `x'[1,`is2'] = (-1)^`j2'*`t' matrix M_zlc`num' = nullmat(M_zlc`num')\ `x' matrix M_zps`num' = nullmat(M_zps`num'),`C' matrix `x'[1,`is2'] = 0 local j2 = `j2' + 1 } local is2 = `is2' + 1 } matrix `x'[1,`is1'] = 0 local j1 = `j1' + 1 } local is1 = `is1' + 1 } * B and s matrix `x' = J(1,`n',`s') matrix `im' = J(1,`n'+1,0) matrix `im'[1,1] = 1 local pm = `n'+1 while `pm' <= `n'+1{ /* set previous digits to 0 */ while `pm'>1{ matrix `x'[1,`pm'-1] = `s' matrix `im'[1,`pm'] = 0 local pm = `pm' - 1 } matrix M_zlc`num' = nullmat(M_zlc`num')\ `x' matrix M_zps`num' = nullmat(M_zps`num'), `B' * matrix list `im' local pm = 1 while `im'[1,`pm'] == 1{ local pm = `pm' + 1 } * disp "pm = " `pm' /* pm is first incomplete digit */ if `pm'<= `n' + 1{ matrix `im'[1,`pm'] = 1 matrix `x'[1,`pm'-1] = -`s' local pm = `pm' - 1 } } * B and r matrix `x' = J(1,`n',0) local ir = 1 while `ir' <= `n'{ local j = 1 while `j' <= 2{ matrix `x'[1,`ir'] = (-1)^`j'*`r' matrix M_zlc`num' = M_zlc`num'\ `x' matrix M_zps`num' = M_zps`num',`B' matrix `x'[1,`ir'] = 0 local j = `j' + 1 } local ir = `ir' + 1 } * A matrix M_zps`num' = M_zps`num',`A' matrix `x' = J(1,`n',0) matrix M_zlc`num' = M_zlc`num'\ `x' matrix M_zlc`num' = M_zlc`num'' global S_1 = `num' end capture program drop spher6 program define spher6 * Stroud 7-1 (p. 318) tempname r s t B C D A x X y a norm im local n = 6 local num = 2^`n' + 2*`n'^2 + 1 capture matrix drop M_zlc`num' capture matrix drop M_zps`num' /* scalar `r' = sqrt( (3*(8-`n')-(`n'-2)*sqrt(3*(8-`n')))/(2*(5-`n')) ) scalar `s' = sqrt( (3*`n' - 2*sqrt(3*(8-`n')))/(2*(3*`n' - 8)) ) scalar `t' = sqrt( (6+sqrt(3*(8-`n')))/2 ) scalar `B' = (8-`n')/( 8*`r'^6) *2/exp(lngamma(`n'/2)) scalar `C' = 1 /( 2^(`n'+3)*`s'^6) *2/exp(lngamma(`n'/2)) scalar `D' = 1/(16*`t'^6) *2/exp(lngamma(`n'/2)) scalar `A' = 2/exp(lngamma(`n'/2)) - 2*`n'*`B' - 2^`n'*`C' - 2*`n'*(`n'-1)*`D' disp in re "r = " `r' " s = " `s' " t = " `t' disp in re "A = " `A' " B = " `B' " C = " `C' " D = " `D' */ * - - + scalar `r' = sqrt( (3*(8-`n')-(`n'-2)*sqrt(3*(8-`n')))/(5-`n') ) scalar `s' = sqrt( (3*`n' - 2*sqrt(3*(8-`n')))/(3*`n' - 8) ) scalar `t' = sqrt( 6+sqrt(3*(8-`n')) ) scalar `B' = ln(8-`n') - ln(8) - 6*( ln(`r') - ln(2)/2 ) + ln(2) - lngamma(`n'/2) scalar `C' = -(`n'+ 3)*ln(2)-6*( ln(`s') - ln(2)/2 ) + ln(2) - lngamma(`n'/2) scalar `D' = -ln(16) - 6*( ln(`t') - ln(2)/2 ) + ln(2) - lngamma(`n'/2) scalar `A' = ln( 2/exp(lngamma(`n'/2)) - 2*`n'*exp(`B') - 2^`n'*exp(`C') - 2*`n'*(`n'-1)*exp(`D') ) *disp in re "r = " `r' " s = " `s' " t = " `t' *disp in re "A = " `A' " B = " `B' " C = " `C' " D = " `D' * A matrix `X' = J(1,`n',0) matrix `y' = (`A') * Do r and B matrix `x' = J(1,`n',0) local ir = 1 while `ir' <= `n'{ local j = 1 while `j' <= 2 { matrix `x'[1,`ir'] = (-1)^`j'*`r' matrix `X' = nullmat(`X')\ `x' matrix `y' = nullmat(`y'),`B' matrix `x'[1,`ir'] = 0 local j = `j' + 1 } local ir = `ir' + 1 } * Do t and D matrix `x' = J(1,`n',0) local it1 = 1 while `it1' <= `n'{ local j1 = 1 while `j1' <= 2{ matrix `x'[1,`it1'] = (-1)^`j1'*`t' local it2 = `it1' + 1 while `it2' <= `n'{ local j2 = 1 while `j2' <= 2{ * disp "it1 = `it1' and it2 = `it2'" matrix `x'[1,`it2'] = (-1)^`j2'*`t' matrix `X' = `X'\ `x' matrix `y' = `y', `D' matrix `x'[1,`it2'] = 0 local j2 = `j2' + 1 } local it2 = `it2' + 1 } matrix `x'[1,`it1'] = 0 local j1 = `j1' + 1 } local it1 = `it1' + 1 } * Do s and C matrix `x' = J(1,`n',`s') matrix `im' = J(1,`n'+1,0) matrix `im'[1,1] = 1 local pm = `n'+1 while `pm' <= `n'+1{ /* set previous digits to 0 */ while `pm'>1{ matrix `x'[1,`pm'-1] = `s' matrix `im'[1,`pm'] = 0 local pm = `pm' - 1 } matrix `X' = `X'\ `x' matrix `y' = `y', `C' * matrix list `im' local pm = 1 while `im'[1,`pm'] == 1{ local pm = `pm' + 1 } * disp "pm = " `pm' /* pm is first incomplete digit */ if `pm'<= `n' + 1{ matrix `im'[1,`pm'] = 1 matrix `x'[1,`pm'-1] = -`s' local pm = `pm' - 1 } } matrix M_zlc`num' = `X' matrix M_zps`num' = `y' matrix M_zlc`num' = M_zlc`num'' global S_1 = `num' end program define fsrs version 6.0 args r s x n if abs(`s')<1e-10&abs(`r')<1e-10{ matrix `x' = (0,0) scalar `n' = 1 } else if abs(`s')<1e-10{ matrix `x' = (`r',0\ -`r',0\ 0,`r'\ 0,-`r') scalar `n' = 4 } else if abs(`r'-`s')<1e-10{ matrix `x' = (`r',`s'\ -`r',`s'\ `r',-`s'\ -`r',-`s') scalar `n' = 4 } else{ matrix `x' = (`r',`s'\ -`r',`s'\ `r',-`s'\ -`r',-`s'\ `s',`r'\ -`s',`r'\ `s',-`r'\ -`s',-`r') scalar `n' = 8 } end program define int2k version 6.0 args k if `k' == 9{ *E_2^{r^2}: 9-1 local l = 4 local num = 20 tempname r1 r2 r3 r4 s1 s2 s3 s4 B1 B2 B3 B4 /* scalar `r1' = 1.538189001320852*sqrt(2) scalar `r2' = 1.224744871391589*sqrt(2) scalar `r3' = 0.4817165220011443*sqrt(2) scalar `r4' = 2.607349811958554*sqrt(2) scalar `s4' = 0.9663217712794149*sqrt(2) scalar `s1' = 0 scalar `s2' = `r2' scalar `s3' = `r3' scalar `B1' = ln(0.1237222328857347) - ln(_pi) scalar `B2' = ln(0.06544984694978697) - ln(_pi) scalar `B3' = ln(0.5935280476180875) - ln(_pi) scalar `B4' = ln(0.001349017971918148) - ln(_pi) */ scalar `r3' = 1.538189001320852*sqrt(2) scalar `r2' = 1.224744871391589*sqrt(2) scalar `r4' = 0.4817165220011443*sqrt(2) scalar `r1' = 2.607349811958554*sqrt(2) scalar `s1' = 0.9663217712794149*sqrt(2) scalar `s3' = 0 scalar `s2' = `r2' scalar `s4' = `r4' scalar `B3' = ln(0.1237222328857347) - ln(_pi) scalar `B2' = ln(0.06544984694978697) - ln(_pi) scalar `B4' = ln(0.5935280476180875) - ln(_pi) scalar `B1' = ln(0.001349017971918148) - ln(_pi) } if `k' == 11{ *E_2^{r^2}: 11-1 local l = 5 local num = 28 tempname r1 r2 r3 r4 r5 s1 s2 s3 s4 s5 B1 B2 B3 B4 B5 /* scalar `r1' = 2.757816396257008*sqrt(2) scalar `r2' = 1.732050807568877*sqrt(2) scalar `r3' = 0.6280515301597559*sqrt(2) scalar `r4' = 1.224744871391589*sqrt(2) scalar `r5' = 0.7071067811865475*sqrt(2) scalar `s4' = 2.121320343559643*sqrt(2) scalar `s5' = 1.224744871391589*sqrt(2) scalar `s1' = 0 scalar `s2' = 0 scalar `s3' = 0 scalar `B1' = ln(0.0008176645817675417) - ln(_pi) scalar `B2' = ln(0.04363323129985824) - ln(_pi) scalar `B3' = ln(0.5373255214498174) - ln(_pi) scalar `B4' = ln(0.003636102608321520) - ln(_pi) scalar `B5' = ln(0.09817477042468103) - ln(_pi) */ scalar `r1' = 2.757816396257008*sqrt(2) scalar `r3' = 1.732050807568877*sqrt(2) scalar `r5' = 0.6280515301597559*sqrt(2) scalar `r2' = 1.224744871391589*sqrt(2) scalar `r4' = 0.7071067811865475*sqrt(2) scalar `s2' = 2.121320343559643*sqrt(2) scalar `s4' = 1.224744871391589*sqrt(2) scalar `s1' = 0 scalar `s3' = 0 scalar `s5' = 0 scalar `B1' = ln(0.0008176645817675417) - ln(_pi) scalar `B3' = ln(0.04363323129985824) - ln(_pi) scalar `B5' = ln(0.5373255214498174) - ln(_pi) scalar `B2' = ln(0.003636102608321520) - ln(_pi) scalar `B4' = ln(0.09817477042468103) - ln(_pi) } if `k' == 15{ *E_2^{r^2}: 15-1 local l = 9 local num = 44 tempname r1 r2 r3 r4 r5 r6 r7 r8 r9 s1 s2 s3 s4 s5 s6 s7 s8 s9 B1 B2 B3 B4 B5 B6 B7 B8 B9 /* scalar `r1' = 3.538388728121807*sqrt(2) scalar `r2' = 2.359676416877929*sqrt(2) scalar `r3' = 1.312801844620926*sqrt(2) scalar `r4' = 0.5389559482114205*sqrt(2) scalar `r5' = 2.300279949805658*sqrt(2) scalar `r6' = 1.581138830084189*sqrt(2) scalar `r7' = 0.8418504335819279*sqrt(2) scalar `r8' = 2.685533581755341*sqrt(2) scalar `r9' = 1.740847514397403*sqrt(2) scalar `s8' = 1.112384431771456*sqrt(2) scalar `s9' = 0.7210826504868960*sqrt(2) scalar `s1' = 0 scalar `s2' = 0 scalar `s3' = 0 scalar `s4' = 0 scalar `s5' = `r5' scalar `s6' = `r6' scalar `s7' = `r7' scalar `B1' = ln(0.000008006483569659628) - ln(_pi) scalar `B2' = ln(0.003604577420838264) - ln(_pi) scalar `B3' = ln(0.1187609330759137) - ln(_pi) scalar `B4' = ln(0.4372488543791402) - ln(_pi) scalar `B5' = ln(0.00003671735075832989) - ln(_pi) scalar `B6' = ln(0.005654866776461627) - ln(_pi) scalar `B7' = ln(0.1777774268424240) - ln(_pi) scalar `B8' = ln(0.0002735449647853290) - ln(_pi) scalar `B9' = ln(0.02087984556938594) - ln(_pi) */ scalar `r1' = 3.538388728121807*sqrt(2) scalar `r4' = 2.359676416877929*sqrt(2) scalar `r7' = 1.312801844620926*sqrt(2) scalar `r9' = 0.5389559482114205*sqrt(2) scalar `r2' = 2.300279949805658*sqrt(2) scalar `r5' = 1.581138830084189*sqrt(2) scalar `r8' = 0.8418504335819279*sqrt(2) scalar `r3' = 2.685533581755341*sqrt(2) scalar `r6' = 1.740847514397403*sqrt(2) scalar `s3' = 1.112384431771456*sqrt(2) scalar `s6' = 0.7210826504868960*sqrt(2) scalar `s1' = 0 scalar `s4' = 0 scalar `s7' = 0 scalar `s9' = 0 scalar `s2' = `r2' scalar `s5' = `r5' scalar `s8' = `r8' scalar `B1' = ln(0.000008006483569659628) - ln(_pi) scalar `B4' = ln(0.003604577420838264) - ln(_pi) scalar `B7' = ln(0.1187609330759137) - ln(_pi) scalar `B9' = ln(0.4372488543791402) - ln(_pi) scalar `B2' = ln(0.00003671735075832989) - ln(_pi) scalar `B5' = ln(0.005654866776461627) - ln(_pi) scalar `B8' = ln(0.1777774268424240) - ln(_pi) scalar `B3' = ln(0.0002735449647853290) - ln(_pi) scalar `B6' = ln(0.02087984556938594) - ln(_pi) } tempname x n capture matrix drop M_zlc`num' local i = 1 while `i'<=`l'{ fsrs `r`i'' `s`i'' `x' `n' * matrix list `x' matrix M_zlc`num' = nullmat(M_zlc`num')\ `x' matrix `x' = J(1,`n',`B`i'') matrix M_zps`num' = nullmat(M_zps`num'), `x' local i = `i' + 1 } matrix M_zlc`num' = M_zlc`num'' global S_1 = `num' end program define multn version 6.0 * Stroud (1971) p. 317, 5-2 args n tempname r s A B C x local num = 2*`n'^2 + 1 capture matrix drop M_zlc`num' capture matrix drop M_zps`num' scalar `r' = sqrt((`n'+2)) scalar `s' = sqrt((`n'+2)/2) scalar `A' = ln(2) - ln(`n'+2) scalar `B' = ln(4-`n')-ln(2) -2*ln(`n'+2) scalar `C' = -2*ln(`n'+2) matrix M_zps`num' = (`A') matrix `x' = J(1,`n',0) matrix M_zlc`num' = `x' local ir = 1 while `ir' <= `n'{ local j = 1 while `j' <= 2{ matrix `x'[1,`ir'] = (-1)^`j'*`r' matrix M_zlc`num' = M_zlc`num'\ `x' matrix M_zps`num' = M_zps`num',`B' matrix `x'[1,`ir'] = 0 local j = `j' + 1 } local ir = `ir' + 1 } local is1 = 1 while `is1' <= `n'{ local j1 = 1 while `j1' <= 2{ matrix `x'[1,`is1'] = (-1)^`j1'*`s' local is2 = `is1' + 1 while `is2' <= `n'{ local j2 = 1 while `j2' <= 2{ * disp "is1 = `is1' and is2 = `is2'" matrix `x'[1,`is2'] = (-1)^`j2'*`s' matrix M_zlc`num' = M_zlc`num'\ `x' matrix M_zps`num' = M_zps`num',`C' matrix `x'[1,`is2'] = 0 local j2 = `j2' + 1 } local is2 = `is2' + 1 } matrix `x'[1,`is1'] = 0 local j1 = `j1' + 1 } local is1 = `is1' + 1 } matrix M_zlc`num' = M_zlc`num'' global S_1 = `num' end program define ghquad * stolen from rfprobit (Bill Sribney) version 4.0 local n `1' tempname xx ww a b local i 1 local m = int((`n' + 1)/2) matrix M_zlc`n' = J(1,`m',0) matrix M_zps`n' = M_zlc`n' while `i' <= `m' { if `i' == 1 { scalar `xx' = sqrt(2*`n'+1)-1.85575*(2*`n'+1)^(-1/6) } else if `i' == 2 { scalar `xx' = `xx'-1.14*`n'^0.426/`xx' } else if `i' == 3 { scalar `xx' = 1.86*`xx'-0.86*M_zlc`n'[1,1] } else if `i' == 4 { scalar `xx' = 1.91*`xx'-0.91*M_zlc`n'[1,2] } else { local im2 = `i' -2 scalar `xx' = 2*`xx'-M_zlc`n'[1,`im2'] } hermite `n' `xx' `ww' matrix M_zlc`n'[1,`i'] = `xx' matrix M_zps`n'[1,`i'] = ln(`ww') - ln(_pi)/2 local i = `i' + 1 } if mod(`n', 2) == 1 { matrix M_zlc`n'[1,`m'] = 0} /* start in tails */ matrix `b' = (1,1) matrix M_zps`n' = M_zps`n'#`b' matrix M_zps`n' = M_zps`n'[1,1..`n'] matrix `b' = (1,-1) matrix M_zlc`n' = M_zlc`n'#`b' matrix M_zlc`n' = M_zlc`n'[1,1..`n'] /* other alternative (start in centre) */ /* matrix `b' = J(1,`n',0) local i = 1 while ( `i'<=`n'){ matrix `b'[1, `i'] = M_zlc`n'[1, `n'+1-`i'] local i = `i' + 1 } matrix M_zlc`n' = `b' local i = 1 while ( `i'<=`n'){ matrix `b'[1, `i'] = M_zps`n'[1, `n'+1-`i'] local i = `i' + 1 } matrix M_zps`n' = `b' */ /* end other alternative */ matrix M_zlc`n' = M_zlc`n'* sqrt(2) end program define hermite /* integer n, scalar x, scalar w */ * stolen from rfprobit (Bill Sribney) version 4.0 local n "`1'" local x "`2'" local w "`3'" local last = `n' + 2 tempname i p matrix `p' = J(1,`last',0) scalar `i' = 1 while `i' <= 10 { matrix `p'[1,1]=0 matrix `p'[1,2] = _pi^(-0.25) local k = 3 while `k'<=`last'{ matrix `p'[1,`k'] = `x'*sqrt(2/(`k'-2))*`p'[1,`k'-1] /* */ - sqrt((`k'-3)/(`k'-2))*`p'[1,`k'-2] local k = `k' + 1 } scalar `w' = sqrt(2*`n')*`p'[1,`last'-1] scalar `x' = `x' - `p'[1,`last']/`w' if abs(`p'[1,`last']/`w') < 3e-14 { scalar `w' = 2/(`w'*`w') exit } scalar `i' = `i' + 1 } di in red "hermite did not converge" exit 499 end program define lebesque version 5.0 local n `1' tempname pt a b scalar `a' = 1/`n' matrix M_zps`n' = J(1,`n',`a') local i = 1 local m = int((`n' + 1)/2) matrix M_zlc`n' = J(1,`m',0) while(`i'<=`m'){ scalar `pt' = `i'/`n' -1/(2*`n') matrix M_zlc`n'[1,`i']=invnorm(`pt') local i = `i' + 1 } /* start in tails */ matrix `b' = (1,-1) matrix M_zlc`n' = M_zlc`n'#`b' matrix M_zlc`n' = M_zlc`n'[1,1..`n'] /* other alternative: left to right */ /* while ( `i'<=`n'){ matrix M_zlc`n'[1, `i'] = -M_zlc`n'[1, `n'+1-`i'] local i = `i' + 1 } */ end program define disprand version 6.0 * displays additional information about random effects * disp "running disprand " disp " " if "e(tplv)" == ""{ * estimates not found exit } tempname var b se cor mn0 mm0 matrix `b' = e(b) local names: colnames(`b') tempname M_nrfc M_nip M_nbrf M_nffc M_b V M_frld matrix `V' = e(V) matrix `M_nrfc' = e(nrfc) matrix `M_nip' = e(nip) matrix `M_nbrf' = e(nbrf) matrix `M_nffc' = e(nffc) matrix `M_frld' = e(frld) local ngeqs = e(ngeqs) local bmat = e(bmat) if `bmat' ==1{matrix `M_b' = e(mb)} local bmat = e(bmat) local iscor = e(cor) local nxt = `M_nffc'[1,colsof(`M_nffc')]+1 local free = e(free) local tplv = e(tplv) local lev1 = e(lev1) local tprf = e(tprf) local cip = e(cip) local nats = e(nats) if `free'{ tempname M_np matrix `M_np' = e(mnp) } local nrfold = `M_nrfc'[2,1] if `M_nbrf'[1,1]>0{ if `lev1' == 1 {disp in gr "Variance at level 1"} else if `lev1' == 2 {disp in gr "Squared Coefficient of Variation"} else if `lev1' == 3 {disp in gr "Dispersion at level 1"} *disp in smcl in gr "{hline 78}" _n di in gr _dup(78) "-" _n if `M_nbrf'[1,1]==1{ if `nats'{ scalar `var' = `b'[1, `nxt']^2 scalar `se' = 2*sqrt(`var'*`V'[`nxt',`nxt']) disp in gr " " in ye `var' " (" `se' ")" } else{ scalar `var' = exp(2*`b'[1, `nxt']) scalar `se' = 2*`var'*sqrt(`V'[`nxt',`nxt']) disp in gr " " in ye `var' " (" `se' ")" } local nxt = `nxt' + 1 } else{ local log "log " if `nats'{ local log } disp " " if `lev1'==1{disp in gr " equation for `log'standard deviation: "} else if `lev1'==2{disp in gr " equation for `log'coefficient of variation"} else if `lev1'==3{disp in gr " equation for `log'(sqrt(phi))"} disp " " local i = 1 while `i' <= `M_nbrf'[1,1]{ scalar `var' = `b'[1,`nxt'] scalar `se' = sqrt(`V'[`nxt',`nxt']) local nna: word `nxt' of `names' disp in gr " `nna': " in ye `var' " (" `se' ")" local i = `i' + 1 local nxt = `nxt' + 1 } } } if `tplv' > 1{ local lev = 2 if `free' == 1{ disp " " disp in gr "Probabilities and locations of random effects" } else{ disp " " disp in gr "Variances and covariances of random effects" } *disp in smcl in gr "{hline 78}" _n di in gr _dup(78) "-" _n *disp in gr "-----------------------------------------------------------------------------" while (`lev' <= `tplv'){ local nip = `M_nip'[1,`lev'] local sof = `M_nrfc'[2,`lev'-1] /* `M_nrfc'[1,`lev'-1] */ disp " " local cl = `tplv' - `lev' + 1 local cl: word `cl' of `e(clus)' disp in gr "***level `lev' (" in ye "`cl'" in gr ")" if `free' == 1{ tempname M_zps`lev' matrix `M_zps`lev'' = e(zps`lev') } local i2 = `M_nrfc'[2,`lev'] local i1 = `nrfold'+1 local num = `i2' -`i1' + 1 /* number of random effects */ if `free'==0{ * get standard errors of variances from those of cholesky decomp. *disp "sechol `lev' `num' `nxt'" qui sechol `lev' `num' `nxt' } local k = 1 local i = `i1' local ii = 1 local nrfold = `M_nrfc'[2,`lev'] while `i'<= `i2'{ local n=`M_nip'[2,`i'] if `free'==1{ tempname M_zlc`n' matrix `M_zlc`n'' = e(zlc`n') local j = 2 local zz=string(`M_zlc`n''[1,1],"%7.0gc") local mm "`zz'" scalar `mn0' = `M_zlc`n''[1,1]*exp(`M_zps`lev''[1,1]) while `j'<=`nip'{ scalar `mn0' = `mn0' + `M_zlc`n''[1,`j']*exp(`M_zps`lev''[1,`j']) local zz=string(`M_zlc`n''[1,`j'],"%7.0gc") local mm "`mm'" ", " "`zz'" local j = `j' + 1 } disp " " disp in gr " loc`ii': " in ye "`mm'" } local j = `i1' local jj = 1 while (`j'<=`i'){ if `free'==1{ local m = `M_nip'[2,`j'] capture tempname M_zlc`m' matrix `M_zlc`m'' = e(zlc`m') scalar `mm0'=0 local mm = 1 while `mm'<=`nip'{ scalar `mm0' = `mm0' + `M_zlc`m''[1,`mm']*exp(`M_zps`lev''[1,`mm']) local mm = `mm' + 1 } local l = 1 scalar `var' = 0 while `l'<=`nip'{ scalar `var' = `var' + /* */ (`M_zlc`n''[1,`l']-`mn0')*(`M_zlc`m''[1,`l']-`mm0')*exp(`M_zps`lev''[1,`l']) local l = `l' + 1 } if `i' == `j'{ disp in gr " var(`ii'): " in ye `var' local nb = `M_nbrf'[1,`ii'+`sof'] if `nb'>1{ disp " " disp in gr " loadings for random effect " `ii' local load = `nxt' + `nb' - 1 if `M_frld'[1,`j']==0{ local nna: word `load' of `names' disp in gr " `nna': " in ye 1 " (fixed)" } *disp in gr " coefficient of" local load = 1 while `load'<=`nb'-1{ local nna: word `nxt' of `names' scalar `var'=`b'[1,`nxt'] scalar `se' = sqrt(`V'[`nxt',`nxt']) disp in gr " `nna': " in ye `var' " (" `se' ")" local nxt = `nxt' + 1 local load = `load' + 1 } disp " " } local nxt = `nxt' + 1 /* skip location parameter */ * disp "increased nxt to" `nxt' if `i'==`i2'{ local l = 2 local zz=string(exp(`M_zps`lev''[1,1]),"%6.0gc") if `nip'>1{ local mm "0`zz'" } else{ local mm "1" } while `l'<=`nip'{ local zz=string(exp(`M_zps`lev''[1,`l']),"%6.0gc") local mm "`mm'" ", " "0`zz'" local l = `l' + 1 } disp in gr " prob: " in ye "`mm'" *new *** display log odds with standard errors local npar = `M_np'[1,`lev'] if `npar' == 0 { local npar = 1} if `npar'>1{ disp " " disp in gr " log odds parameters" } local l = 1 while `l'< `nip'{ if `npar'>1{ if `l'>1 { disp " " } disp in gr " class `l'" } local load = 1 local ncur = `nxt' + (`l' - 1)*`num' + (`l' - 1)*`npar' - 1 while `load'<=`npar'{ local ncur = `ncur' + 1 * disp "nxt = `nxt', nb = `nb', num = `num' and ncur = `ncur'" local nna: word `ncur' of `names' scalar `var'=`b'[1,`ncur'] scalar `se' = sqrt(`V'[`ncur',`ncur']) if `npar'>1{ disp in gr " `nna': " in ye `var' " (" `se' ")" } local load = `load' + 1 } local l = `l' + 1 } local nxt = `ncur' + 1 } * end new } else{ disp in gr "cov(`ii',`jj'): " in ye `var' } } else{/* free=0 */ * disp "k= " `k' ", i= " `i' ", j= " `j' ", ii= " `ii' ", jj= " `jj' scalar `var' = M_cov[`ii', `jj'] scalar `se' = sqrt(M_se[`k', `k']) if `i' == `j'{ disp " " disp in gr " var(`ii'): " in ye `var' " (" `se' ")" local nb = `M_nbrf'[1,`ii'+`sof'] if `nb'>1{ disp " " disp in gr " loadings for random effect " `ii' local load = `nxt' + `nb' -1 if `M_frld'[1,`j']==0{ local nna: word `load' of `names' disp in gr " `nna': " in ye 1 " (fixed)" } * disp in gr " coefficient of" local load = 1 while `load'<=`nb'-1{ local nna: word `nxt' of `names' * disp "nxt = " `nxt' scalar `var'=`b'[1,`nxt'] scalar `se' = sqrt(`V'[`nxt',`nxt']) disp in gr " `nna': " in ye `var' " (" `se' ")" local nxt = `nxt' + 1 local load = `load' + 1 } disp " " } * skip variance parameter local nxt = `nxt' + 1 } else{ if `iscor'==0{ disp in gr " cov(`ii',`jj'): " in ye "fixed at 0" } else{ scalar `cor' = `var'/sqrt(M_cov[`ii',`ii']*M_cov[`jj',`jj']) disp in gr " cov(`ii',`jj'): " in ye `var' " (" `se' ")" /* */ " cor(`ii',`jj'): " `cor' *local nxt = `nxt' + 1 } } } local j = `j' + 1 local jj = `jj' + 1 local k = `k' + 1 } local i = `i' + 1 local ii = `ii' + 1 } local lev = `lev' + 1 /* skip off-diagonal cholesky parameters */ if `iscor'~=0&`free'==0{local nxt = `nxt' + `num'*(`num'-1)/2} /* -1? */ *disp "next nxt is " `nxt' /* if `free'{ local nxt = `nxt'+(`nip'-1)*(`num'+1) if `cip'==0{ local nxt = `nxt'+`num' } local nxt = `nxt' - 1 } */ *disp "next nxt is " `nxt' } if `tprf'>1&`bmat'==1{ disp " " disp in gr "B-matrix:" *disp in smcl in gr "{hline 78}" _n di in gr _dup(78) "-" _n disp " " disp " " * disp "nxt = " `nxt' local i = 1 while `i'<`tprf'{ local j = 1 while `j' < `tprf'{ if `M_b'[`i',`j']>0{ scalar `var' =`b'[1,`nxt'] scalar `se' = sqrt(`V'[`nxt',`nxt']) disp in gr " B(`i',`j'): " in ye `var' " (" `se' ")" local nxt = `nxt' + 1 } local j = `j' + 1 } local i = `i' + 1 } } if `ngeqs'>0{ disp " " disp in gr "Regressions of latent variables on covariates" *disp in smcl in gr "{hline 78}" _n di in gr _dup(78) "-" _n disp " " tempname mngeqs matrix `mngeqs' = e(mngeqs) local i = 1 while `i'<=`ngeqs'{ local k = `mngeqs'[1,`i'] local n = `mngeqs'[2,`i'] disp in gr " random effect " in ye `k'-1 in gr " has " in ye `n' in gr " covariates:" local nxt2 = `nxt'+`n'-1 local j = 1 while `j' <= `n'{ local nna: word `nxt' of `names' scalar `var'=`b'[1,`nxt'] scalar `se' = sqrt(`V'[`nxt',`nxt']) disp in gr " `nna': " in ye `var' " (" `se' ")" local nxt = `nxt' + 1 local j = `j' + 1 } local i = `i' + 1 } } *disp in smcl in gr "{hline 78}" _n di in gr _dup(78) "-" _n disp " " } /* endif toplv >1 */ end program define sechol version 6.0 args lev num nxt * num is number of random effects local l = `num'*(`num' + 1)/2 *disp "lev = `lev' num = `num' nxt = `nxt' l= `l'" tempname b V C L zero a H M_nbrf M_nrfc ind matrix `M_nbrf' = e(nbrf) matrix `M_nrfc' = e(nrfc) local iscor = e(cor) matrix `b' = e(b) matrix `V' = e(V) local sof = `M_nrfc'[2,`lev'-1] /* was `M_nrfc'[1,`lev'-1] */ local i = 1 local k = 1 matrix `C' = J(`l',`l',0) matrix `L' = J(`num',`num',0) matrix `ind' = `L' * get L matrix while `i' <= `num'{ * skip loading parameters local nb = `M_nbrf'[1,`i'+`sof'] local nxt = `nxt' + `nb' -1 disp "nxt = " `nxt' matrix `L'[`i',`i'] = `b'[1, `nxt'] * matrix list `L' matrix `ind'[`i',`i'] = `nxt' local nxt = `nxt' + 1 local i = `i' + 1 } local i = 2 while `i' <= `num'&`iscor'==1{ local j = 1 while `j' < `i'{ matrix `L'[`i',`j'] = `b'[1, `nxt'] matrix `ind'[`i',`j'] = `nxt' local nxt = `nxt' + 1 local j = `j' + 1 } local i = `i' + 1 } * disp "L and ind" * matrix list `L' * matrix list `ind' * get C matrix local ll1 = 1 local i = 1 while `i' <= `num'{ local j = 1 while `j' <= `i'{ local nxt1 = `ind'[`i', `j'] local ll2 = 1 local ii = 1 while `ii' <= `num'{ local jj = 1 while `jj' <= `ii'{ local nxt2 = `ind'[`ii', `jj'] disp "ll1 = " `ll1' " ll2 = " `ll2' " nxt1 = " `nxt1' " nxt2 = " `nxt2' if `iscor' == 1{ matrix `C'[`ll1', `ll2'] = `V'[`nxt1',`nxt2'] matrix `C'[`ll2', `ll1'] = `C'[`ll1', `ll2'] } else if `i'==`j'&`ii'==`jj'{ matrix `C'[`ll1', `ll2'] = `V'[`nxt1',`nxt2'] matrix `C'[`ll2', `ll1'] = `C'[`ll1', `ll2'] } local ll2 = `ll2' + 1 local jj = `jj' + 1 } local ii = `ii' + 1 } local ll1 = `ll1' + 1 local j = `j' + 1 } local i = `i' + 1 } * disp "C" * matrix list `C' matrix `zero' = J(`num', `num', 0) local k = 1 local i = 1 local n = `num' * (`num' + 1)/2 matrix `H' = J(`n',`n',0) while `i' <= `num' { local j = 1 while `j' <= `i'{ * derivative of LL' with respect to i,j th element of L mat `a' = `zero' mat `a'[`i',`j'] = 1 mat `a' = `a'*(`L')' mat `a' = `a' + (`a')' disp "a" * matrix list `a' local ii = 1 local kk = 1 while `ii'<=`num'{ local jj = 1 while `jj' <= `ii'{ matrix `H'[`kk',`k'] = `a'[`ii',`jj'] local jj = `jj' + 1 local kk = `kk' + 1 } local ii= `ii' + 1 } local j = `j' + 1 local k = `k' + 1 } local i = `i' + 1 } * disp "H" * matrix list `H' matrix M_se = `H'*`C'*(`H')' matrix M_cov = `L'*(`L')' * matrix list M_se * matrix list M_cov end program define timer version 6.0 end