*! version 2.2.16 SRH 18 July 2006
program define gllam_ll
version 6.0
args todo bo lnf junk1 junk2 what res
* what = 1: update posterior means and standard deviations
* what = 2: posterior probabilities
* what = 3: cluster likelihood contributions*cluster frequency weights
* what = 4: cluster likelihood contributions
* what = 5-7: posterior mean of statistic

/*
noi disp " "
noi disp $HG_SD1[$which]
noi disp $HG_SD2[$which]
noi disp $HG_C21[$which]
noi disp " "
*/

if "`what'"==""{
    local what = 0
}
if `what' == 1 {
    global HG_zip zipg1
    if $HG_free{
        global HG_zip zipf1
    }
}
else if $HG_adapt{
    global HG_zip zipga
}

*disp in re "what=" `what'

* matrix list `bo'
if $HG_dots {
    noi disp in gr "." _c
}
/* ----------------------------------------------------------------------------- */
/* set up variables and macros needed */

tempname b mzlc

local toplev = $HG_tplv
local topi = $HG_tpi
local clus $HG_clus

sort `clus'

* reset the the clock and reset znow
local i = 1
matrix M_ip[1,1] = 1 /* in case topi=0 */
while (`i' <= `topi'){
    matrix M_ip[1,`i'] = 1
    local i = `i' + 1
}

/* -------------------------------------------------------------------------------- */
/* set up lprod1 ... lint1 */

local i=1
tempvar extra
quietly gen double `extra'=0
while `i' <= `toplev'{
    tempvar  lint`i'
    gen double `lint`i'' = 0.0   /* used to integrate, therefore must be zero */
    tempvar lfac`i'
    quietly gen double `lfac`i'' = 0
    if (`i'>1){
        tempvar  lprod`i'
        quietly gen double `lprod`i'' = .
        tempvar lfac`i'
        quietly gen double `lfac`i'' = 0
    }
    local i = `i' + 1
}

/* set up names for HG_xb`i' */
local i = 1
while (`i' <= $HG_tpff){
    tempname junk
    global HG_xb`i' "`junk'"
    local i = `i' + 1
}
/* set up names for HG_s`i' */
local i = 1
while (`i'<=$HG_tprf){
    tempname junk
    global HG_s`i' "`junk'"
    local i = `i'+1
}


if $HG_free{
/* set up names for HG_p`lev'`k' */
    local lev = 2
    while `lev'<=$HG_tplv{
        local npar = M_np[1,`lev']
        if `npar'>0{
            local k = 1
            while `k'<=M_nip[1, `lev']{
                * disp in re "creating HG_p`lev'`k'"
                tempname junk
                global HG_p`lev'`k' "`junk'"
                local k = `k' + 1
            }
        }
        local lev = `lev' + 1
    }
}


/* set up names for HG_E`rf'`lev' and  HG_V`rf'`lev' */
if `what'==1{
    local lev = 1
    while `lev'<=$HG_tplv{
        if `lev'<$HG_tplv{
            local maxrf =M_nrfc[2,`lev'+1]
        }
        else{
            local maxrf =M_nrfc[2,$HG_tplv]
        }

        local rf = 1
        while `rf'<`maxrf'{
            * disp in re "creating HG_V`rf'`rf'`lev'"
            tempname junk
            global HG_E`rf'`lev' "`junk'"
            gen double ${HG_E`rf'`lev'}=0
            tempname junk
            global HG_V`rf'`rf'`lev' "`junk'"
            gen double ${HG_V`rf'`rf'`lev'}=0
            local rf2 = `rf' + 1
            while `rf2'<$HG_tprf{ /* lower diagonal elements */
                * disp in re "creating HG_V`rf2'`rf'`lev'"
                tempname junk
                global HG_V`rf2'`rf'`lev' "`junk'"
                gen double ${HG_V`rf2'`rf'`lev'}=0
                local rf2 = `rf2' + 1
            }
            local rf = `rf' + 1
        }

        local lev = `lev' + 1
    }
}
else if $HG_adapt{
    local llev = 1
    while `llev'<$HG_tplv{
        local lrf = M_nrfc[2,`llev']
        while `lrf'<M_nrfc[2,`llev'+1]{
            tempname junk
            global HG_E`lrf'`llev' "`junk'"
            gen double ${HG_E`lrf'`llev'}=0
            local lrf = `lrf' + 1
        }
        local llev = `llev' + 1
    }
}

/* set up names for posterior mean of statistic */
if `what'>=5{
    tempvar lpred
    local lev = 1
    while `lev'<=$HG_tplv{
        tempname junk
        global HG_stat`lev' "`junk'"
        gen double ${HG_stat`lev'}=0
        local lev =`lev' + 1
    }
}


qui remcor "`bo'"

*noi matrix list `bo'
/*
if `npar'>0{
disp in re "after remcor: HG_p21[$which] = "  $HG_p21[$which]
disp in re "after remcor: HG_p22[$which] = "  $HG_p22[$which]
}
*/

if $HG_error==1{
    disp in re "error in remcor"
    scalar `lnf' = .
    exit
}

* disp "HG_xb1 after remcor: HG_xb1[$which] = " $HG_xb1[$which]
* disp "HG_s2 after remcor: HG_s2[$which] = " $HG_s2[$which]
* disp "HG_s1 after remcor: HG_s1[$which] = $HG_s1 = " $HG_s1[$which]

if $HG_adapt{
    *noi disp " "
    local i = $HG_tprf - 1
    while `i' >= 1 {
        local ip = `i' + 1
        local k = `i' - 1
        while `k' >= 1 {
            local kp = `k' + 1
            * disp in re "replace HG_s`ip' = HG_s`ip' + HG_C`i'`k'*HG_s`kp'"
            qui replace ${HG_s`ip'} = ${HG_s`ip'} + ${HG_C`i'`k'}*${HG_s`kp'}
            local k = `k' - 1
        }
        local i = `i' - 1
    }

    local i = 2
    tempname junk
    global HG_zuoff "`junk'"
    qui gen double $HG_zuoff = 0
    while `i'<=$HG_tprf{
        local im = `i' - 1
        * disp in re "replace HG_zuoff = HG_zuoff + HG_s`i'*HG_MU`im'"
        qui replace $HG_zuoff = $HG_zuoff + ${HG_s`i'}*${HG_MU`im'}

        * disp in re "replace HG_s`i' = HG_SD`im'*HG_s`i'"
        qui replace ${HG_s`i'} = ${HG_SD`im'}*${HG_s`i'}

        local i = `i' + 1
    }
    *dis in re "HG_zuoff[$which] = " $HG_zuoff[$which]
}


local i = `toplev'
while `i' > 1 {
    *noi disp "$HG_zip `i'$
    quietly $HG_zip `i'
    local i = `i' - 1
}


timer on 2
* disp "STARTING LOOP"
/* --------------------------------------------------------------------------------- */
/* recursive loop through r.effs (levels)    */
/* topi nested loops: irf from 1 ro nip(rf) */
/* ip is "clock": (irf) stages of loops      */

local levno = `toplev'
local rf = `topi'
while (`rf' <= `topi') { /* for each r.eff */

/* ----------------------------------------------------------------------------------*/
/* reset ip to 1st point for all lower r.effs... */
/* update znow                                   */
    * disp "reset ip up to random effect " `rf'
    while (`rf' > 1) {
        local rf = `rf' - 1
        * disp `rf'
        matrix M_ip[1,`rf'] = 1
    }
    while (`levno' > 1){
    /* update znow for all new ips    */
        *noi disp "$HG_zip `levno'$
        $HG_zip `levno'
        local levno = `levno' - 1
    }
/* --------------------------------------------------------------------------------- */
/* set lint1 to lpyz for new znow  */

    local rf = 1
    local levno = 1
    local sortlst `clus' /* cluster variables to aggregate over */

    *matrix list M_ip
    timer on 3
    qui lpyz `lint1'
    timer off 3

    if `what'>=5{
        if $HG_tplv>1{
            matrix score double `lpred' = M_znow
        }
        else{
            qui gen double `lpred' = 0
        }
        if $HG_adapt{
            qui replace `lpred' = $HG_xb1 + $HG_zuoff + `lpred' - $HG_offs
        }
        else{
            qui replace `lpred' = $HG_xb1 + `lpred' - $HG_offs
        }
    timer on 5
        gllas_yu $HG_stat1 `lpred' `res' `what'
        sort `sortlst' /* added Oct 2004 */
    timer off 5
        qui drop `lpred'
    }

    *noi disp in re " after lpyz lint1 = " `lint1'[$which]
    *noi summ `lint1' if `lint1'==-100
    quietly count if `lint1'==.& $HG_ind==1
    if r(N) > 0{
        * overflow problem
        * disp in re "overflow at level 1 ( " r(N) " missing values)"
        * list $HG_clus if `lint1'==.
        *matrix list `bo'
        *lpyz `lint1'
        if `what'!=3{ /* what = 3 used by robust & need valid ll contributions */
            scalar `lnf' = .
            exit
        }
    }
    quietly replace `lint1' = `lint1' * $HG_wt1

/* --------------------------------------------------------------------------------- */
/* update lint for all completed levels up to */
/* highest completed level, reset lower lints */
/* to zero (for models including a random effect) */

    while (M_ip[1,`rf'] == M_nip[1,`rf'] & `rf' <= `topi'){
    * digit equals its max => increment next digit
        if (`rf' == M_nrfc[1,`levno'] & `levno' < `toplev'){
        * done last r.eff of current level
            * disp "********** level " `levno' " complete ************"

            * next level
            local lprev = `levno'
            local levno = `levno' + 1

            /* change sortlst */
            local prvsort `sortlst'
            local l = `toplev' - `levno' + 2
            tokenize "`sortlst'"
            * take away var. of level to sum over
            local `l' " "
            local sortlst "`*'"
/*------------------------------------------------------------------------------------ */
/* change lprod`levno' and  */
/* update lint`levno'       */

            if "`f`levno''"==""{ /* first term for lint`levno' */
                local f`levno'=1
                * disp "first term for level `levno'"
            }
            else{
                local f`levno'=0
                * disp "next term for level `levno'"
            }
timer on 4
if $HG_noC1 {
    lnupdate `levno' `lint`levno'' `lint`lprev'' `lprod`levno'' /*
        */ `lfac`levno'' `lfac`lprev'' `extra' /*
        */ `lnf' "`prvsort'" "`sortlst'" `f`levno'' `what'
}
else {
    _gllamm_lnu, levno(`levno') lintlv(`lint`levno'') /*
        */ lintprv(`lint`lprev'') lprodlv(`lprod`levno'') /*
    */ lfaclv(`lfac`levno'') lfacprv(`lfac`lprev'') extra(`extra') /*
        */ lnf(`lnf')  prvsort(`prvsort') sortlst(`sortlst') /*
    */ first(`f`levno'') what(`what')
}
timer off 4

            local f`lprev'

*args levno lintlv lintprv lprodlv lfaclv lfacprv extra lnf prevsort sortlst first
/* ------------------------------------------------------------------------------------------- */
        } /* next digit */
        local rf = `rf' + 1
    }
    * rf is first r.eff that is not complete
    * increase clock in lowest incomplete digit
    * disp "update rf = " `rf'
    matrix M_ip[1,`rf'] = M_ip[1,`rf'] + 1
}
timer off 2
*quietly{
    *now rf too high
    *!! disp "********** level " `toplev' " complete ************"
    *a noi disp "lint" `toplev' "[" $which "] = " `lint`toplev''[$which]
    if(`toplev'>1){
        if `what'==1{
            local rf = 1
            while `rf'< $HG_tprf{
                * disp in re "setting HG_MU`rf' and HG_SD`rf'"
                * disp in re "by `sortlst': replace HG_MU`rf' = HG_E`rf'`toplev'/lint`toplev'[_N]"
                qui by `sortlst': replace ${HG_MU`rf'} = ${HG_E`rf'`toplev'}/`lint`toplev''[_N]
                qui by `sortlst': replace ${HG_SD`rf'} = ${HG_V`rf'`rf'`toplev'}/`lint`toplev''[_N] - ${HG_MU`rf'}^2
                qui replace ${HG_SD`rf'} = cond(${HG_SD`rf'}>0,sqrt(${HG_SD`rf'}),0)
                *noi disp " "
                *noi disp "HG_MU`rf'[$which] = " ${HG_MU`rf'}[$which]
                *summ ${HG_MU`rf'}
                *noi disp "HG_SD`rf'[$which] = " ${HG_SD`rf'}[$which]
                *summ ${HG_SD`rf'}
                local rf2 =  1
                while `rf2' < `rf'{ /* must use MU's that are already calculated */
                    qui by `sortlst': replace ${HG_C`rf'`rf2'} = ${HG_V`rf'`rf2'`toplev'}/`lint`toplev''[_N] - (${HG_MU`rf2'}*${HG_MU`rf'})
                    *noi disp "HG_C`rf'`rf2'[$which] = " ${HG_C`rf'`rf2'}[$which]
                    *summ  ${HG_C`rf'`rf2'}
                    local rf2 = `rf2' + 1
                }
                local rf = `rf' + 1
            }

            if $HG_adapt==1{ /* have recalculated means and covariances and need diff. ones for adaptive */
                prepadpt
            }
        }
        else if `what'>=5{
            if $HG_post{
                qui by `sortlst': replace ${HG_stat`toplev'} = ${HG_stat`toplev'}/`lint`toplev''[_N]
            }
        }
        else if `what'==2{
            local i = 1
            while `i'<= M_nip[1,2] {
                qui by `sortlst': replace ${HG_p`i'} = ${HG_p`i'}/`lint`toplev''[_N]
                local i = `i' + 1
            }
        }
        *a disp "taking log of lint" `toplev' " = " `lint`toplev''[$which]
        *a disp "subtracting " `lfac`toplev''[$which]

        ** begin junk (conditional log-likelihood for capture-recapture)
/*
        qui replace `lint`toplev'' = exp(ln(`lint`toplev'')-`lfac`toplev'')
        qui summ `lint`toplev'', meanonly
        qui replace `lint`toplev'' = ${HG_wt`toplev'}*(ln(`lint`toplev'') - ln(r(sum)))
*/
        ** end junk

        quietly replace `lint`toplev'' = (ln(`lint`toplev'')-`lfac`toplev'')* ${HG_wt`toplev'}
    }
    *noi display "lint" `toplev' "[" $which "] = " `lint`toplev''[$which]
    if `what'==3{
        if `toplev'==1{ /* by sorlst won't work for composite link */
            qui replace `res' = `lint`toplev''
        }
        else{
            qui by `sortlst': replace `res' = `lint`toplev''[_N]
        }
    }
    else if `what'==4{
        if `toplev'==1{ /* by sorlst won't work for composite link */
            qui replace `res' = `lint`toplev''/ ${HG_wt`toplev'}
        }
        else{
            qui by `sortlst': replace `res' = `lint`toplev''[_N]/ ${HG_wt`toplev'}
        }
    }
    else if `what'>=5{
        qui replace `res' = ${HG_stat`toplev'}
    }
    qui by `sortlst': replace  `extra' = cond(_n==_N,1,0)
    *mlsum `lnf' = `lint`toplev'' if `extra' == 1 /* can only use this when program called by ML */
    *noi summ `lint`toplev'' if `extra' == 1
    *list $HG_clus `lint`toplev'' if `extra' ==1 & `lint`toplev''< -300

    qui count if `extra' == 1
    local n = r(N)
    summarize `lint`toplev'', meanonly
    if `n' > r(N) {
        * noi disp "there are " r(N) " values of likelihood, should be " `n'
        * noi list $HG_clus if `extra' == 1& `lint`toplev''==.
        * noi matrix list `bo'
        * noi disp "lnf equal to missing in last step"
        scalar `lnf' = .
        exit
    }
    scalar `lnf' = r(sum)
    * display "total lnf = " `lnf'
    * capture drop lint`toplev'
    * gen double lint`toplev' = `lint`toplev''
*} /* qui */
end

program define prepadpt
*qui replace $HG_C21 = 0
*qui replace $HG_MU1 = 0
*qui replace $HG_MU2 = 0
*qui replace $HG_SD1 = 1
*qui replace $HG_SD2 = 1
*qui replace $HG_C31 = 0
*qui replace $HG_C32 = 0
*disp in re " "

    local tplv = $HG_tplv

    local i = $HG_tprf - 1
    while `i' >= 1 {
        * disp in re " "
        qui replace ${HG_MU`i'} = 0 if ${HG_MU`i'} ==.
        *qui replace ${HG_SD`i'} = 1e-05 if ${HG_SD`i'} < 1e-25 | ${HG_SD`i'} ==.
        qui replace ${HG_SD`i'} = 1e-05 if ${HG_SD`i'} < 1e-05 | ${HG_SD`i'} ==.
* variance
        * disp in re "HG_V`i'`i'`tplv' = HG_SD`i'^2"
        qui replace ${HG_V`i'`i'`tplv'} = ${HG_SD`i'}^2
        local j =  $HG_tprf - 1
        while `j' > `i'{
            * disp in re "HG_V`i'`i'`tplv' = HG_V`i'`i'`tplv' -  HG_C`j'`i'*HG_V`j'`i'`tplv' "
            ** new
            qui replace ${HG_C`j'`i'} = 0 if ${HG_V`i'`i'`tplv'} - ${HG_C`j'`i'}*${HG_V`j'`i'`tplv'} < 0
            qui replace ${HG_V`i'`i'`tplv'} = ${HG_V`i'`i'`tplv'} - ${HG_C`j'`i'}*${HG_V`j'`i'`tplv'}

            qui replace ${HG_MU`i'} = ${HG_MU`i'} - ${HG_C`j'`i'}*${HG_MU`j'}
            * disp in re "HG_MU`i' = HG_MU`i' - HG_C`j'`i'*HG_MU`j' = " in ye ${HG_MU`i'}[$which]
            * noi summ ${HG_MU`i'}
            local j = `j' - 1
        }

        * disp in re "negative variances"
        * noi summ ${HG_V`i'`i'`tplv'} if ${HG_V`i'`i'`tplv'} <0
        * noi list $HG_clus ${HG_V`i'`i'`tplv'} if ${HG_V`i'`i'`tplv'} <0
        * qui replace ${HG_V`i'`i'`tplv'} = 1e-10 if ${HG_V`i'`i'`tplv'} <0
        qui replace ${HG_SD`i'} = sqrt(${HG_V`i'`i'`tplv'})
        *disp in re "HG_MU`i' = " in ye in ye ${HG_MU`i'}[$which]
        *disp in re "HG_SD`i' = " in ye ${HG_SD`i'}[$which]
        * noi summ ${HG_SD`i'}

* covariances and betas
        local k = `i' - 1
        while `k' >= 1 {
            qui replace ${HG_V`i'`k'`tplv'} = ${HG_C`i'`k'}
            * disp in re "HG_V`i'`k'`tplv' = HG_C`k'`i'"
            local j =  $HG_tprf - 1
            while `j' > `i'{
                *disp in re "HG_V`i'`k'`tplv' = HG_V`i'`k'`tplv' - HG_C`j'`i'*HG_V`j'`k'`tplv'"

                qui replace ${HG_V`i'`k'`tplv'} = ${HG_V`i'`k'`tplv'} - ${HG_C`j'`i'}*${HG_V`j'`k'`tplv'}
                local j = `j' - 1
            }
            qui replace ${HG_C`i'`k'} = ${HG_V`i'`k'`tplv'}/${HG_V`i'`i'`tplv'}
            * noi summ ${HG_C`i'`k'}
            *disp in re "HG_C`i'`k' = " in ye ${HG_C`i'`k'}[$which]
            local k = `k' - 1
        }

        local i = `i' - 1
    }
*qui replace $HG_C21 = 0
/*
qui replace $HG_MU2 = .57022599
qui replace $HG_SD2 = .83303331
qui replace $HG_MU1 = -.55565732
qui replace $HG_SD1 = .65728191
*/
end

program define lnupdate
version 6.0
args levno lintlv lintprv lprodlv lfaclv lfacprv extra lnf prvsort sortlst first what
tempvar lpkpl
quietly{
    * disp in re "!!! update level " `levno'
    * noi matrix list M_znow

    /* set previous lint to ln(lint) */
    local lprev = `levno' - 1
    if(`levno' > 2){
        *!! disp " replace lint" `lprev' " by ln(lint" `lprev' ")"
        quietly count if `lintprv' < 1e-308
        if r(N) > 0{
            /* overflow problem */
            * noi disp "overflow at level " `lprev'
            scalar `lnf' = .
            exit
        }
        if `what'==1{
            local rf = 1
            while `rf'< M_nrfc[2,`lprev']{
                * disp "by `prvsort': replace HG_E`rf'`lprev' = HG_E`rf'`lprev'/lintprv[_N]"
                    qui by `prvsort': replace ${HG_E`rf'`lprev'} = ${HG_E`rf'`lprev'}/`lintprv'[_N]
                * disp in re "by `prvsort': replace HG_V`rf'`rf'`lprev' = HG_V`rf'`rf'`lprev'/lintprv[_N]"
                qui by `prvsort': replace ${HG_V`rf'`rf'`lprev'} = ${HG_V`rf'`rf'`lprev'}/`lintprv'[_N]
                local rf2 = `rf' + 1
                while `rf2' < $HG_tprf{
                    * disp in re "by `prvsort': replace HG_V`rf2'`rf'`lprev' = HG_V`rf2'`rf'`lprev'/lintprv[_N]"
                    qui by `prvsort': replace ${HG_V`rf2'`rf'`lprev'} = ${HG_V`rf2'`rf'`lprev'}/`lintprv'[_N]
                    local rf2 = `rf2' + 1
                }
                local rf = `rf' + 1
            }

        }
        if `what'>=5{
            if $HG_post{
                qui by `prvsort': replace ${HG_stat`lprev'} = ${HG_stat`lprev'}/`lintprv'[_N]
            }
        }
        quietly replace `lintprv' = ln(`lintprv')
        quietly replace `lintprv' = (`lintprv'-`lfacprv')*${HG_wt`lprev'}
    }

    /* sum previous lprod within cluster at current level */
    *!! disp " "
    *!! disp "by `sortlst': replace lprod" `levno' "=cond(_n==N, sum(lint" `lprev' "))"
    * noi summ `lintprv'
    quietly by `sortlst': replace `lprodlv' = cond(_n==_N,sum(`lintprv'),.)
    * noi disp " "
    * noi disp "lprod" `levno' " = " `lprodlv'[$which]

    /* accumulate terms for integral */

    /* get lpkpl: log of product of r.effs at level */
    qui $HG_lzpr `levno' `lpkpl'
    * noi disp "exp(lprod`leno'+lpkpl) = " exp(`lpkpl'[$which]+`lprodlv'[$which])

    if `first' { /* first term for lint`levno' */
        quietly replace `extra' = 0
        quietly replace `lfaclv' = -`lprodlv' - `lpkpl'
        *a noi disp " "
        * noi disp "lfac`levno' = " `lfaclv'[$which]
        * noi disp "lintlv = 1"
                qui replace `lintlv' = 1
    }
    else{
        local max = 500
        quietly replace `extra' = cond(`lprodlv'+ `lpkpl'+`lfaclv'>`max', /*
                    */ -(`lprodlv'+`lpkpl'+`lfaclv')+`max',0)
        * disp "extra = " `extra'[$which]
        quietly replace `lfaclv'=`lfaclv'+`extra'
        * noi disp "lfac`levno' = " `lfaclv'[$which]

                /* increment lint at current level using lprod at previous level */
            *a disp "increase lint" `levno' " by exp(lprodlv + lpkpl +lfaclv)"
            quietly replace `lintlv' = exp(`extra')*`lintlv' + exp(`lprodlv'+ `lpkpl'+`lfaclv')
                * noi disp "increase by " exp(`lprodlv'[$which]+`lpkpl'[$which]+`lfaclv'[$which]) " to "  `lintlv'[$which]
    }


/* posterior means and variances*/
    if `what'==1{
        local rf = 1
        while `rf'< M_nrfc[2,`levno'] {
            * noi disp "update `rf' `levno'"
            quietly by `sortlst': replace ${HG_E`rf'`levno'}=/*
            */ exp(`extra'[_N])*${HG_E`rf'`levno'}+ ${HG_E`rf'`lprev'}*exp(`lprodlv'[_N]+`lpkpl'+`lfaclv'[_N])
            * noi disp "HG_E" `rf' `levno' "[" $which "] = "  ${HG_E`rf'`levno'}[$which]
            quietly by `sortlst': replace ${HG_V`rf'`rf'`levno'}=/*
            */ exp(`extra'[_N])*${HG_V`rf'`rf'`levno'}+ ${HG_V`rf'`rf'`lprev'}*exp(`lprodlv'[_N]+`lpkpl'+`lfaclv'[_N])
            *noi disp "HG_V" `rf' `rf' `levno' "[" $which "] = "  ${HG_V`rf'`rf'`levno'}[$which]
            local rf2 = `rf' + 1
            while `rf2' < $HG_tprf {
                * noi disp "HG_V`rf2'`rf'`levno' is ${HG_V`rf2'`rf'`levno'}"
                quietly by `sortlst': replace ${HG_V`rf2'`rf'`levno'}=/*
                */ exp(`extra'[_N])*${HG_V`rf2'`rf'`levno'}+ ${HG_V`rf2'`rf'`lprev'}*exp(`lprodlv'[_N]+`lpkpl'+`lfaclv'[_N])
                * noi disp "HG_V" `rf2' `rf' `levno' "[" $which "] = "  ${HG_V`rf2'`rf'`levno'}[$which]
                local rf2 = `rf2' + 1
            }
            local rf = `rf' + 1
        }

    }
    if `what'>=5{
        if $HG_post{
            * disp in re "sortlist: `sortlst'"
            *sort `sortlst' /* removed oct 2004 */
            quietly by `sortlst': replace ${HG_stat`levno'}=/*
            */ exp(`extra'[_N])*${HG_stat`levno'}+ /*
            */ ${HG_stat`lprev'}*exp(`lprodlv'[_N]+`lpkpl'+`lfaclv'[_N])
        }
        else{
            *disp in re "HG_stat`levno'[1] = " ${HG_stat`levno'}[1] " + " ${HG_stat`lprev'}[1] " * " exp(`lpkpl')
            qui replace ${HG_stat`levno'} = ${HG_stat`levno'} + ${HG_stat`lprev'}*exp(`lpkpl')
        }
    }

    else if `what'==2{
        * noi matrix list M_ip
        local i = M_ip[1,2]
        local j = 1
        while `j'<`i'{
            quietly by `sortlst': replace ${HG_p`j'} = exp(`extra'[_N])*${HG_p`j'}
            local j = `j' + 1
        }
        quietly by `sortlst': replace ${HG_p`i'} = exp(`lprodlv'[_N]+`lpkpl'+`lfaclv'[_N])
    }
    /* reset previous lint to zero */
    if `levno'>2{
        if `what'==1{
            local rf = 1
            while `rf'<M_nrfc[2,`lprev']{
                * disp in re "replace HG_V`rf'`rf'`lprev' = 0"
                    qui replace ${HG_E`rf'`lprev'} = 0
                qui replace ${HG_V`rf'`rf'`lprev'} = 0
                local rf2 = `rf' + 1
                while `rf2' < $HG_tprf{
                    * disp in re "replace HG_V`rf2'`rf'`lprev' = 0"
                    qui replace ${HG_V`rf2'`rf'`lprev'} = 0
                    local rf2 = `rf2' + 1
                }
                local rf = `rf' + 1
            }

        }
        *!! disp "setting lint" `lprev' " to zero"
        quietly replace `lintprv' = 0
        quietly replace `lfacprv' = 0
        *!!!!new
        if `what'>=5{
            *disp in re "setting HG_stat" `lprev' " to zero"
            quietly replace ${HG_stat`lprev'} = 0
        }
    }
 } /* qui */
end

program define zipf
    version 6.0
    * updates znow
    * matrix list M_ip
    args levno

    * disp "in zip, levno is " `levno'
    local i = M_nrfc[2,`levno'-1] + 1

    *!! disp "update"
    *  same class for all random effects
    local k = M_nrfc[1,`levno']
    local k = M_ip[1,`k']
    local last = M_nrfc[2,`levno']
    while `i' <= `last'{
        local npt = M_nip[2,`i']
        local im = `i' - 1
        * disp "     "`im' "th z to " `which' "th location"
        * disp " using M_zlc`npt' "
        matrix M_znow[1,`im'] = M_zlc`npt'[1,`k']
        local i = `i' + 1
    }
end

program define zipf1
    version 6.0
    * updates znow
    * matrix list M_ip
    args levno
    tempname mzlc

    * disp "in zip, levno is " `levno'
    local i = M_nrfc[2,`levno'-1] + 1

    *!! disp "update"
    *  same class for all random effects
    local k = M_nrfc[1,`levno']
    local k = M_ip[1,`k']
    local last = M_nrfc[2,`levno']
    while `i' <= `last'{
        local npt = M_nip[2,`i']
        local im = `i' - 1
        * disp "     "`im' "th z to " `which' "th location"
        * disp " using M_zlc`npt' "
        matrix M_znow[1,`im'] = M_zlc`npt'[1,`k']
        local i = `i' + 1
    }
    local llev = `levno' - 1
    local im = M_nrfc[2,`llev']
    while `im' < `last'{
* change HG_E and HG_V
        scalar `mzlc' = M_znow[1,`im']
        qui replace ${HG_E`im'`llev'} = ${HG_MU`im'} + ${HG_SD`im'}*`mzlc'
        local im2 = $HG_tprf - 1
        while `im2' > `im'{
            scalar `mzlc' = M_znow[1,`im2']
            * disp in re "replace HG_E`im'`llev' = HG_E`im'`llev' + HG_C`im2'`im'*(HG_MU`im2' + HG_SD`im2'*`mzlc')"
            qui replace ${HG_E`im'`llev'} = ${HG_E`im'`llev'} + ${HG_C`im2'`im'}*(${HG_MU`im2'} + ${HG_SD`im2'}*`mzlc')
            local im2 = `im2' - 1
        }
        * noi disp "HG_E`im'`llev'[$which] = " ${HG_E`im'`llev'}[$which]
        qui replace ${HG_V`im'`im'`llev'} = ${HG_E`im'`llev'}^2

        local im = `im' + 1
    }
    local im = M_nrfc[2,`llev']
    while `im' < `last'{
* covariances with same level and higher level effects
        local llev2 = `llev'
        local im2 = `im' + 1
        while `llev2'<$HG_tplv{
            while `im2' < M_nrfc[2,`llev2' + 1]{
                qui replace ${HG_V`im2'`im'`llev'} = ${HG_E`im'`llev'}*${HG_E`im2'`llev2'}
                * disp in re "HG_V`im2'`im'`llev' = HG_E`im'`llev'*HG_E`im2'`llev2' = " in ye ${HG_V`im2'`im'`llev'}[$which]
                local im2 = `im2' + 1
            }
            local llev2 = `llev2' + 1
        }
        local im = `im' + 1
    }
end

program define zipg
    version 6.0
    * updates znow
    * matrix list M_ip
    args levno

    * disp "in zip, levno is " `levno'
    local i = M_nrfc[2,`levno'-1] + 1

    *!! disp "update"
    if $HG_mult{
        *local npt = M_nip[2,`i']
        local npt = M_nip[2,`levno']
        local k = M_nrfc[1,`levno']
        local k = M_ip[1,`k']
        local f = M_nrfc[2,`levno'-1]
    }

    local last = M_nrfc[2,`levno']
    while `i' <= `last'{
        local im = `i' - 1

        if $HG_mult{
            matrix M_znow[1,`im'] = M_zlc`npt'[`i'-`f',`k']
        }
        else{
            local npt = M_nip[2,`i']
            local which = M_ip[1,`i']
            * disp "     "`im' "th z to " `which' "th location"
            * disp " using M_zlc`npt' "
            matrix M_znow[1,`im'] = M_zlc`npt'[1,`which']
        }

        *!! disp M_znow[1,`im']
        local i = `i' + 1
    }
end


program define zipga
    version 6.0
    * updates znow
    * matrix list M_ip
    args levno
    tempname mzlc

    * disp "in zip, levno is " `levno'
    local llev = `levno' - 1
    local i = M_nrfc[2,`llev'] + 1

    *!! disp "update"
    if $HG_mult{
        *local npt = M_nip[2,`i']
        local npt = M_nip[2,`levno']
        local k = M_nrfc[1,`levno']
        local k = M_ip[1,`k']
        local f = M_nrfc[2,`levno'-1]
    }
    local last = M_nrfc[2,`levno']

    while `i' <= `last'{
        *local npt = M_nip[2,`i']
        local im = `i' - 1

        if $HG_mult{
            matrix M_znow[1,`im'] = M_zlc`npt'[`i'-`f',`k']
        }
        else{
            local npt = M_nip[2,`i']
            local which = M_ip[1,`i']
            * disp "     "`im' "th z to " `which' "th location"
            * disp " using M_zlc`npt' "
            matrix M_znow[1,`im'] = M_zlc`npt'[1,`which']
        }
        local i = `i' + 1
    }
    local im = M_nrfc[2,`llev']
    while `im' < `last'{

* change HG_E
        scalar `mzlc' = M_znow[1,`im']
        qui replace ${HG_E`im'`llev'} = ${HG_MU`im'} + ${HG_SD`im'}*`mzlc'
        *local im2 = M_nrfc[2,`llev']

        local im2 = $HG_tprf - 1
        while `im2' > `im'{
            scalar `mzlc' = M_znow[1,`im2']
            * disp in re "replace HG_E`im'`llev' = HG_E`im'`llev' + HG_C`im2'`im'*(HG_MU`im2' + HG_SD`im2'*`mzlc')"
            qui replace ${HG_E`im'`llev'} = ${HG_E`im'`llev'} + ${HG_C`im2'`im'}*(${HG_MU`im2'} + ${HG_SD`im2'}*`mzlc')
            local im2 = `im2' - 1
        }
        * noi disp "HG_E`im'`llev'[$which] = " ${HG_E`im'`llev'}[$which]

        *!! disp M_znow[1,`im']
        local im = `im' + 1
    }
end


program define zipg1
    version 6.0
    * updates znow
    * matrix list M_ip
    args levno
    tempname mzlc

    * noi disp "in zip, levno is " `levno'
    local llev = `levno' - 1
    local i = M_nrfc[2,`llev'] + 1

    *!! disp "update"
    if $HG_mult{
        *local npt = M_nip[2,`i']
        local npt = M_nip[2,`levno']
        local k = M_nrfc[1,`levno']
        local k = M_ip[1,`k']
        local f = M_nrfc[2,`levno'-1]
    }
    local last = M_nrfc[2,`levno']

    while `i' <= `last'{
        local im = `i' - 1
        *local npt = M_nip[2,`i']
        if $HG_mult{
            * disp in re "M_znow[1,`im'] = M_zlc`npt'[`i'-`f',`k']"
            matrix M_znow[1,`im'] = M_zlc`npt'[`i'-`f',`k']
        }
        else{
            local npt = M_nip[2,`i']
            local which = M_ip[1,`i']
            * disp "     "`im' "th z to " `which' "th location"
            * disp " using M_zlc`npt' "
            matrix M_znow[1,`im'] = M_zlc`npt'[1,`which']
        }

        local i = `i' + 1
    }
    local im = M_nrfc[2,`llev']
    while `im' < `last'{
* change HG_E and HG_V
        scalar `mzlc' = M_znow[1,`im']
        qui replace ${HG_E`im'`llev'} = ${HG_MU`im'} + ${HG_SD`im'}*`mzlc'
        local im2 = $HG_tprf - 1
        while `im2' > `im'{
            scalar `mzlc' = M_znow[1,`im2']
            * disp in re "replace HG_E`im'`llev' = HG_E`im'`llev' + HG_C`im2'`im'*(HG_MU`im2' + HG_SD`im2'*`mzlc')"
            qui replace ${HG_E`im'`llev'} = ${HG_E`im'`llev'} + ${HG_C`im2'`im'}*(${HG_MU`im2'} + ${HG_SD`im2'}*`mzlc')
            local im2 = `im2' - 1
        }
        * noi disp "HG_E`im'`llev'[$which] = " ${HG_E`im'`llev'}[$which]
        qui replace ${HG_V`im'`im'`llev'} = ${HG_E`im'`llev'}^2

        local im = `im' + 1
    }
    local im = M_nrfc[2,`llev']
    while `im' < `last'{
* covariances with same level and higher level effects
        local llev2 = `llev'
        local im2 = `im' + 1
        while `llev2'<$HG_tplv{
            while `im2' < M_nrfc[2,`llev2' + 1]{
                qui replace ${HG_V`im2'`im'`llev'} = ${HG_E`im'`llev'}*${HG_E`im2'`llev2'}
                * disp in re "HG_V`im2'`im'`llev' = HG_E`im'`llev'*HG_E`im2'`llev2' = " in ye ${HG_V`im2'`im'`llev'}[$which]
                local im2 = `im2' + 1
            }
            local llev2 = `llev2' + 1
        }
        local im = `im' + 1
    }
end


program define lzprobf
    version 6.0
    * for free masses
    * returns product of pk needed for integration at level lev for current ip
    args levno lpkpl
    tempname mzps mznow
    * disp in re "in zprob, levno is " `levno'

    local i=M_nrfc[1,`levno'-1] + 1

    *!! disp "-----------lpkpl: sum of log of"

    local npt = M_nip[2,`i']
    * disp "     prob for " `i' "th r.eff: " `which' "th weight"
    * disp " using M_zps`npt' "

    local which = M_ip[1,`i']
    if M_np[1,`levno']>0{
        qui gen double `lpkpl' = ${HG_p`npt'`which'}
    }
    else{
        qui gen double `lpkpl' =  M_zps`npt'[1,`which']
    }

    * disp in re "lpkpl[$which] = " `lpkpl'[$which]
end

program define lzprobg
    version 6.0
    * product gaussian quadrature
    * returns product of pk needed for integration at level lev for current ip
    args levno lpkpl
    tempname mzps mznow
    * disp in re "in zprob, levno is " `levno'
    qui gen double `lpkpl' = 0
    local lv = `levno' - 1
    local i=M_nrfc[1,`lv'] + 1

    *!! disp "-----------lpkpl: sum of log of"
    local last = M_nrfc[1,`levno']
    while `i' <= `last'{
        local npt = M_nip[2,`i']
        * disp "     prob for " `i' "th r.eff: " `which' "th weight"
        * disp " using M_zps`npt' "

        local which = M_ip[1,`i']
        local im = `i' - 1
        scalar `mzps' = M_zps`npt'[1,`which']
        qui replace `lpkpl' = `lpkpl'+ `mzps'
        if $HG_adapt{
            scalar `mznow' = M_znow[1,`im']
            *qui replace `lpkpl' = `lpkpl' + ln(${HG_SD`im'}) + `mznow'^2/2 - (${HG_MU`im'} + ${HG_SD`im'}*`mznow')^2/2
            qui replace `lpkpl' = `lpkpl' + ln(${HG_SD`im'}) + `mznow'^2/2 - (${HG_E`im'`lv'})^2/2
        }
        local i=`i'+1
    }
    * disp in re "lpkpl[$which] = " `lpkpl'[$which]
end

program define lzprobm
    version 6.0
    * mult version
    * returns product of pk needed for integration at level lev for current ip
    args levno lpkpl
    tempname mzps mznow
    * disp in re "in zprob, levno is " `levno'

    local lv = `levno' - 1
    local i=M_nrfc[1,`lv'] + 1

    *!! disp "-----------lpkpl: sum of log of"

    local npt = M_nip[2,`i']
    * disp "     prob for " `i' "th r.eff: " `which' "th weight"
    * disp " using M_zps`npt' "

    local which = M_ip[1,`i']
    scalar `mzps' = M_zps`npt'[1,`which']
    qui gen double `lpkpl' = `mzps'

    if $HG_adapt{
        local i=M_nrfc[2,`lv'] + 1
        local last = M_nrfc[2,`levno']
        while `i' <= `last'{
            local im = `i' - 1
            scalar `mznow' = M_znow[1,`im']
            qui replace `lpkpl' = `lpkpl' + ln(${HG_SD`im'}) + `mznow'^2/2 - (${HG_E`im'`lv'})^2/2
            *disp in re "HG_E`im'`lv'[$which] = " ${HG_E`im'`lv'}[$which]
            local i = `i' + 1
        }
    }
    *disp in re "lpkpl[$which] = " `lpkpl'[$which]
end


program define lpyz
    version 6.0
* returns log of prob of obs. given znow
    args lpyz

    * disp "-----------------called lpyz"

    tempvar zu xb mu /* linear predictor and zu: r.eff*design matrix for r.eff */

/* ----------------------------------------------------------------------------- */
*quietly{


    if $HG_tprf>1{

        matrix score double `zu' = M_znow
        if $HG_adapt{
            qui replace `zu' = $HG_zuoff + `zu'
            *qui replace `zu' = $HG_E11*$HG_s2 + $HG_E21*$HG_s3
        }


    }
    else{
        qui gen double `zu' = 0
    }

    * matrix list M_znow
    * disp "ML_y1: $ML_y1 " $ML_y1[$which]
    * matrix list M_ip
    disp " xb1 = " $HG_xb1[$which]
    disp " zu = " `zu'[$which]


    if $HG_mlog>0{
        nominal `lpyz' `zu'
    }

    if $HG_oth{
        local myand
        if "$HG_lv"~=""&($HG_nolog>0|$HG_mlog>0){
            local myand $HG_lvolo~=1
        }
        quietly gen double `mu' = 0
        *if $HG_noC|$HG_comp>0 {
    if $HG_noC {
            link "$HG_link" `mu' $HG_xb1 `zu' $HG_s1
            if $HG_comp>0 {
                compos `mu' "`myand'"
            }
            *disp " mu = " `mu'[$which]
            if $HG_comp>0{
                if "`myand'"~=""{
                    local myand `myand' & $HG_ind>0
                }
                else{
                    local myand $HG_ind>0
                }
            }
            family "$HG_famil" `lpyz' `mu' "`myand'"
        }

        else {
            if $HG_lev1 != 0 {
                local s1opt "st($HG_s1)"
            }
            if "$HG_denom" != "" {
                local denopt "denom($HG_denom)"
            }
            if "$HG_fv" != "" {
                local fvopt "fv($HG_fv)"
            }
            if "$HG_lv" != "" {
                local lvopt "lv($HG_lv)"
                local othopt "oth(M_oth)"
            }
            if "`myand'" != "" {
                local ifopt "if `myand'"
            }
            if $HG_comp > 0 { /* only got here if new Stata 8 */
                local comp "comp($HG_comp)"
                local cclus "cluster($HG_clus)"
            }
            noi _gllamm_fl `lpyz' `mu' `ifopt', `s1opt' /*
            */ link($HG_link) family($HG_famil) `denopt' `fvopt' /*
            */ `lvopt' xb($HG_xb1) zu(`zu') /*
            */ y($ML_y1) `othopt' `comp' `cclus'
        }
        disp  " mu = " `mu'[$which]
    }

    if $HG_nolog>0{
        if $HG_noC {
            ordinal `lpyz' `zu'
        }
        else {
            if $HG_lev1 != 0 {
                local stopt st($HG_s1)
            }
            if "$HG_lv"!="" {
                local lvopt lv($HG_lv)
            }
            local j 1
            while `j'<=$HG_tpff {
                local xbeta `xbeta' ${HG_xb`j'}
                local j = `j' + 1
            }
            _gllamm_ord `lpyz', y($ML_y1) xb(`xbeta') /*
            */ zu(`zu') link($HG_linko) nlog($HG_nolog) /*
            */ olog(M_olog) nresp(M_nresp) resp(M_resp) /*
            */ `stopt' `lvopt'
        }
    }

*} /* qui */
end

program define compos
    version 6.0
    args mu und

    tempvar junk mu2
    local ifs
    if "`und'"~=""{
        local ifs if `und'
    }
    gen double `junk'=0
    gen double `mu2'=.
    local i = 1
    *disp in re "in compos: HG_clus is: $HG_clus"
    while `i'<= $HG_comp{
        *disp in re "in compos: variable HG_co`i' is: ${HG_co`i'}"
        replace `junk' = `mu'*${HG_co`i'}
        qui by $HG_clus: replace `junk' = sum(`junk')
        qui by $HG_clus: replace `mu2' = `junk'[_N] if $HG_ind==`i'
        local i = `i' + 1
    }
    qui replace `mu' = `mu2' `ifs'
end

program define nominal
    version 6.0
    args lpyz zu
    tempvar mu

    if $HG_smlog{
        local s $HG_s1
    }
    else{
        local s = 1
    }
    local and
    if "$HG_lv"~=""{
        local and & $HG_lv == $HG_mlog
        local mlif if $HG_lv == $HG_mlog
    }
    disp "mlogit link `mlif'"
    if $HG_exp==1&$HG_expf==0{
        qui gen double `mu' = exp(`zu'/`s') if $ML_y1==M_respm[1,1] `and'
        local n=rowsof(M_respm)
        local i=2
        while `i'<=`n'{
            local prev = `i' - 1
            * disp "xb`prev':" ${HG_xb`prev'}[$which]
            qui replace `mu' = exp((${HG_xb`prev'} + `zu')/`s') if $ML_y1==M_respm[`i',1] `and'
            local i = `i' + 1
        }

        sort $HG_clus $HG_ind
        qui by $HG_clus: replace `lpyz'=cond(_n==_N,sum(`mu'),.) `mlif'
        qui replace `lpyz' = ln(`mu'/`lpyz') `mlif'
/*  skip sort
    qui by $HG_clus: replace `lpyz' = sum(`mu') `mlif'
    qui by $HG_clus: replace `lpyz' = cond($HG_ind==1,ln(`mu'/`lpyz'[_N]),.) `mlif'

*/
    }
    else if $HG_exp==1&$HG_expf==1{
        qui gen double `mu' = exp(($HG_xb1 + `zu')/`s') `mlif'
        sort $HG_clus $HG_ind
        * disp "sort $HG_clus $HG_ind"
        qui by $HG_clus: replace `lpyz'=cond(_n==_N,sum(`mu'),.) `mlif'
        * disp "denom = " `lpyz'[$which]
        qui replace `lpyz' = ln(`mu'/`lpyz') `mlif'
    }
    else{
        tempvar den tmp
        local n=rowsof(M_respm)
        local i = 2
        qui gen double `mu' = 1 if $ML_y1==M_respm[1,1] `mlif'
        qui gen double `den' = 1
        qui gen double `tmp' = 0
        while `i'<= `n'{
            local prev = `i' - 1
            qui replace `tmp' = exp((${HG_xb`prev'} + `zu')/`s') `mlif'
            qui replace `mu' =  `tmp' if $ML_y1==M_respm[`i',1] `mlif'
            replace `den' = `den' + `tmp' `mlif'
            local i = `i' + 1
        }
        replace `lpyz' = ln(`mu'/`den') `mlif'
    }
end

program define ordinal
    version 6.0
    args lpyz zu
    local no = 1
    local xbind = 2
    tempvar mu p1 p2
    qui gen double `p1' = 0
    qui gen double `p2' = 0
    qui gen double `mu' = 0

    while `no' <= $HG_nolog{
        local olog = M_olog[1,`no']
        local lnk: word `no' of $HG_linko

        if "`lnk'"=="ologit"{
            local func logitl
        }
        else if "`lnk'"=="oprobit"{
            local func probitl
        }
        else if "`lnk'"=="ocll"{
            local func cll
        }
        else if "`lnk'"=="soprobit"{
            local func sprobitl
        }
        local and
        if "$HG_lv"~=""&$HG_nolog>0{
            local and & $HG_lv == `olog'
        }
        * disp "ordinal link is `lnk', and = `and'"
        local n=M_nresp[1,`no']

        * disp  "HG_xb1: " $HG_xb1
        * disp  "xbind = " `xbind'
        * disp  ${HG_xb`xbind'}[$which]

        qui replace `mu' = $HG_xb1+`zu'-${HG_xb`xbind'}
        `func' `mu' `p1'
        qui replace `lpyz' = ln(1-`p1') /*
            */ if $ML_y1==M_resp[1,`no'] `and'
        qui replace `p2' = `p1'
        local i = 2
        while `i' < `n'{
            local nxt = `xbind' + `i' - 1

            * disp "nxt = " `nxt'
            * disp ${HG_xb`nxt'}[$which]

            qui replace `mu' = $HG_xb1+`zu'-${HG_xb`nxt'}
            `func' `mu' `p2'

            * disp "p1 and p2: "  `p1'[$which] " " `p2'[$which]

            qui replace `lpyz' = ln(`p1' -`p2') /*
                */ if $ML_y1==M_resp[`i',`no'] `and'
            qui replace `p1' = `p2'
            local i = `i' + 1
        }
        local xbind = `xbind' + `n' -1
        qui replace `lpyz' = ln(`p2') /*
            */ if $ML_y1==M_resp[`n',`no'] `and'
        local no = `no' + 1
    } /* next ordinal response */
    *tab $ML_y1 if `lpyz'==. `and'
    qui replace `lpyz' = -100 if `lpyz'==. `and'
end

program define logitl
    version 6.0
    args mu p
    qui replace `p' = 1/(1+exp(-`mu'))
end

program define cll
    version 6.0
    args mu p
    qui replace `p' = 1-exp(-exp(`mu'))
end

program define probitl
    version 6.0
    args mu p
    qui replace `p' = normprob(`mu')
end

program define sprobitl
    version 6.0
    args mu p
    qui replace `p' = normprob(`mu'/$HG_s1)
end


program define link
    version 6.0
* returns mu for requested link
    args which mu xb zu s1
    * disp " in link, which is `which' "

    tokenize "`which'"
    local i=1
    local ifs
    while "`1'"~=""{
        if "$HG_lv" ~= ""{
            local oth =  M_oth[1,`i']
            local ifs if $HG_lv==`oth'
        }
        * disp "`1' link `ifs'"

        if ("`1'" == "logit"){
            quietly replace `mu' = 1/(1+exp(-`xb'-`zu')) `ifs'
        }
        else if ("`1'" == "probit"){
            * disp "doing probit "
            quietly replace `mu' = normprob((`xb'+`zu')) `ifs'
        }
        else if ("`1'" == "sprobit"){
            quietly replace `mu' = normprob((`xb'+`zu')/`s1') `ifs'
        }
        else if ("`1'" == "log"){
            * disp "doing log "
            quietly replace `mu' = exp(`xb'+`zu') `ifs'
        }
        else if ("`1'" == "recip"){
            * disp "doing recip "
            quietly replace `mu' = 1/(`xb'+`zu') `ifs'
        }
        else if ("`1'" == "cll"){
            * disp "doing cll "
            quietly replace `mu' = 1 - exp(-exp(`xb'+`zu')) `ifs'
        }
        else if ("`1'" == "ll"){
            quietly replace `mu' = exp(-exp(`xb'+`zu')) `ifs'
        }
        else if ("`1'" == "ident"){
            quietly replace `mu' = `xb'+`zu' `ifs'
        }
        local i = `i' + 1
        mac shift
    }

end

program define family
    version 6.0
    args which lpyz mu und

    tokenize "`which'"
    local i=1
    * disp "in family, und = `und'"
    if "$HG_fv" == ""{
        local ifs
        if "`und'"~=""{local und if `und'}
    }
    else{
        if "`und'"~=""{local und & `und'}
    }
    while "`1'"~=""{
        if "$HG_fv" ~=""{
            local ifs if $HG_fv == `i'
        }
        if ("`1'" == "binom"){
            famb `lpyz' `mu' "`ifs'" "`und'"
        }
        else if ("`1'" == "poiss"){
            famp `lpyz' `mu' "`ifs'" "`und'"
        }
        else if ("`1'" == "gauss") {
            *disp in re "famg lpyz mu $HG_s1 `ifs' `und'"
            famg `lpyz' `mu' $HG_s1 "`ifs'" "`und'"  /* get log of conditional prob. */
        }
        else if ("`1'" == "gamma"){
            famga `lpyz' `mu' $HG_s1 "`ifs'" "`und'"
        }
        else{
            disp in re "unknown family in gllam_ll"
            exit 198
        }
        local i = `i' + 1
        mac shift
    }
end

program define famg
    version 6.0
* returns log of normal density conditional on r.effs
    args lpyz mu s1 if and
    * disp "running famg `if' `and'"
    * disp "s1 = " `s1'[$which] ", mu = " `mu'[$which] " and Y = " $ML_y1[$which]
        quietly replace `lpyz' = /*
        */ -(ln(2*_pi*`s1'^2) + (($ML_y1-`mu')/`s1')^2)/2 `if' `and'
end

program define famb
    version 6.0
* returns log of binomial density conditional on r.effs
* $HG_denom is denominator
    args lpyz mu if and
    * disp "running famb `if' `and'"
    * disp "mu = " `mu'[$which] " and Y = " $ML_y1[$which]
    qui replace `lpyz' = cond($ML_y1>0,$ML_y1*ln(`mu'),0)                 /*
        */ + cond($HG_denom-$ML_y1>0,($HG_denom-$ML_y1)*ln(1-`mu'),0) /*
        */ + cond($HG_denom>1,lngamma($HG_denom+1)-lngamma($ML_y1+1)  /*
        */ -  lngamma($HG_denom-$ML_y1+1),0) `if' `and'
    *tab $ML_y1 `if' `and' & `lpyz'==.
    qui replace `lpyz' = cond(`lpyz'==.,-100,`lpyz') `if' `and'
    * disp "done famb"
end

program define famp
    version 6.0
* returns log of poisson density conditional on r.effs
    args lpyz mu if and
    *!! disp "running famp `if'"
    * disp in re "if and: `if' `and'"
    quietly replace `lpyz' = /*
        */ $ML_y1*(ln(`mu'))-`mu'-lngamma($ML_y1+1) `if' `and'
    * qui replace `lpyz' = cond(`lpyz'==.,-100,`lpyz') `if' `and'
    * disp "done famp"
end

program define famga
    version 6.0
* returns log of gamma density conditional on r.effs
    args lpyz mu s1 if and
    *!! disp "running famg `if'"
    *!! disp "mu = " `mu'[$which]
    *!! disp "s1 = " `s1'[$which]
    qui replace `mu' = 0.0001 if `mu' <= 0
    tempvar nu
    qui gen double `nu' = `s1'^(-2)
        quietly replace `lpyz' = /*
        */ `nu'*(ln(`nu')-ln(`mu')) - lngamma(`nu')/*
        */ + (`nu'-1)*ln($ML_y1) - `nu'*$ML_y1/`mu' `if' `and'
end

program define timer
version 6.0
end