corentinchoisy/simul_confusion
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Initial commit

Browse Source
This commit is contained in:
Corentin Choisy
2024-06-04 11:19:21 +02:00
commit a808b2a204
710 changed files with 277632 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

405
Modules/ado/plus/r/remcor.ado Normal file
View File

@ -0,0 +1,405 @@
*! version 2.0.3 SRH 1 June 2003
program define remcor
version 6.0
* takes b and transform set of r.effs at same level using choleski
* returns bi where the r.effs are independent and correlation equations
* have been removed
* also takes exponential of sd parameters
* and evaluates contributions to linear predictor that don't change
args b stop
tempname b2 s1 cov t r d dd u denom mean mzps
tempvar junk
gen double `junk'=0
global HG_error=0
disp "*********in remcor:"
/* fixed effects $HG_xb1, $HG_xb2 etc. (tempnames stored in global macros-list) */
qui disp "fixed parameters: "
qui matrix list `b'
if $HG_const==1{
matrix `b2' = `b'*M_T' + M_a
matrix coleq `b2' = $HG_cole
matrix colnames `b2' = $HG_coln
noi matrix list `b2'
}
else{
matrix `b2' = `b'
}
local nffold=0
local ff = 1
local nxt = 1
while(`ff' <= $HG_tpff){
matrix `b2' = `b2'[1, `nxt'...]
local np = M_nffc[1, `ff'] - `nffold'
qui disp "np = " `np'
if `np'>0{
local nxt = `np' + 1
local nffold = M_nffc[1,`ff']
matrix `s1' = `b2'[1,1..`np']
qui matrix list `s1'
qui disp "tempname: ${HG_xb`ff'}"
matrix score double ${HG_xb`ff'} = `s1' /* nontemp */
}
else{
qui gen double ${HG_xb`ff'} = 0
}
qui disp ${HG_xb`ff'}[$which]
local ff=`ff'+1
}
if "$HG_off"~=""{qui replace $HG_xb1=$HG_xb1+$HG_off}
qui disp "HG_xb1 = " $HG_xb1[$which]
if $HG_ethr{
local j = 1
local ii = 1
while `ii'<=$HG_nolog{
local j = `j' + 1
local jm = `j' + M_nresp[1,`ii']-3
while `j' <= `jm'{
local jp = `j' + 1
* disp in re "replace HG_xb`jp' = HG_xb`j' + exp(HG_xb`jp')"
qui replace ${HG_xb`jp'} = ${HG_xb`j'} + exp(${HG_xb`jp'})
local j = `j' + 1
}
local ii = `ii' + 1
}
}
/* random effects */
/* level 1 */
local np = M_nbrf[1,1]
if `np'>0{
matrix `b2' = `b2'[1, `nxt'...]
local nxt = 1
matrix `s1' = `b2'[1,1..`np']
qui matrix list `s1'
matrix score double $HG_s1 = `s1'
if $HG_nats{
qui replace $HG_s1 = abs($HG_s1)
}
else{
qui replace $HG_s1=exp($HG_s1)
}
qui disp "s1 = $HG_s1 = " $HG_s1[$which]
local nxt = `nxt' + `np'
}
local lev = 2
local rf = 2
local nrfold = M_nrfc[2,1]
/* MASS POINTS */
if($HG_free){
tempname pdenom
gen double `pdenom' = 1.0
while(`lev'<=$HG_tplv&`rf'<=$HG_tprf){
local j1 = M_nrfc[2, `lev']
local nrf = `j1' - `nrfold'
local nip = M_nip[1, `lev']
scalar `denom' = 1 /* =exp(0) */
qui replace `pdenom' = 1.0
matrix M_zps`lev' = J(1,`nip',0)
local k = 1
qui disp "`nip' integration points at level `lev'"
local nloc = `nip'
local npar = M_np[1,`lev']
if $HG_cip|`nip'>1{ local nloc = `nloc'-1}
while `k' <= `nloc' {
local j = `nrfold'+1
while `j'<=`j1'{
qui disp "level `lev', class `k' and random effect `j'"
qui disp " nxt = " `nxt'
matrix `b2' = `b2'[1, `nxt'...]
qui matrix list `b2'
local nxt = 1
if `k'==1{
/* linear predictors come before first masspoint */
local np = M_nbrf[1,`j']-1
if `np'>0 {
qui disp "extracting coefficients for r.eff"
matrix `s1' = `b2'[1,1..`np']
matrix score double ${HG_s`j'} = `s1'
qui disp "HG_s`j' = ${HG_s`j'} = " ${HG_s`j'}[$which]
local nxt = `nxt' + `np'
}
/* first coeff fixed at one */
if M_frld[1,`rf']~=1{
matrix `s1' = (1)
local lab: colnames `b2'
local lab: word `nxt' of `lab'
matrix colnames `s1'=`lab'
qui matrix list `s1'
capture drop `junk'
matrix score double `junk' = `s1'
if `np'>0{
qui replace ${HG_s`j'}=${HG_s`j'}+`junk'
}
else{
qui gen double ${HG_s`j'}=`junk'
}
qui matrix list `s1'
}
qui disp "HG_s`j' = ${HG_s`j'} = " ${HG_s`j'}[$which]
qui disp "making M_zlc`j'"
matrix M_zlc`j' = J(1,`nip',0)
}
matrix M_zlc`j'[1,`k'] = `b2'[1,`nxt']
local nxt = `nxt' + 1
local j = `j' + 1
}
if `k'<`nip'{
if `npar'>0{
qui disp "extract probability parameters for HG_p`cl' (npar=`npar')"
matrix `s1' = `b2'[1,`nxt'..(`nxt'+`npar'-1)]
qui matrix list `s1'
capture drop `junk'
matrix score double `junk' = `s1'
qui gen double ${HG_p`lev'`k'} = `junk'
local nxt = `nxt' + `npar' - 1
qui replace `pdenom' = `pdenom' + exp(${HG_p`lev'`k'})
}
scalar `mzps' = exp(`b2'[1,`nxt'])
if `mzps' == . {
global HG_error=1
exit
}
matrix M_zps`lev'[1,`k'] = `b2'[1,`nxt']
local nxt = `nxt' + 1
scalar `denom' = `denom' + `mzps'
}
local k = `k' + 1
}
if `npar'>0{
qui gen double ${HG_p`lev'`nip'} = 0.0
local k = 1
while `k' <= `nip'{
qui disp "divide HG_p`lev'`k'"
qui replace ${HG_p`lev'`k'} = ${HG_p`lev'`k'} - ln(`pdenom')
local k = `k' + 1
}
}
local k = 1
while `k' <= `nip'{
matrix M_zps`lev'[1,`k'] = M_zps`lev'[1,`k'] - ln(`denom')
local k = `k' + 1
}
local j = `nrfold' + 1
while `j' <= `j1'{ /* define last location */
if $HG_cip == 1{
local k = 1
scalar `mean' = 0
while `k'<`nip'{
scalar `mean' = `mean' + M_zlc`j'[1,`k']*exp(M_zps`lev'[1,`k'])
local k = `k' + 1
}
scalar `mzps' = exp(M_zps`lev'[1,`nip'])
matrix M_zlc`j'[1,`nip'] = -`mean'/`mzps'
}
else if `nip'>1{
matrix M_zlc`j'[1,`nip'] = `b2'[1,`nxt']
local nxt = `nxt' + 1
}
qui disp "M_zlc`j'"
qui matrix list M_zlc`j'
local j = `j' + 1
}
qui disp "M_zps`lev'"
qui matrix list M_zps`lev'
local nrfold = `j1'
local lev = `lev' + 1
}
}/*endif HG_free */
else{
/* ST. DEVS */
qui disp "random parameters: "
if $HG_tprf>1{matrix CHmat = J($HG_tprf-1,$HG_tprf-1,0)}
while(`lev'<=$HG_tplv&`rf'<=$HG_tprf){
local np = M_nbrf[1,`rf']
qui disp "np = " `np'
local nrf = M_nrfc[2, `lev'] - `nrfold'
matrix `t' = J(`nrf',`nrf',0)
local i = 1
while (`i' <= `nrf'){
qui disp " nxt = " `nxt'
matrix `b2' = `b2'[1, `nxt'...]
local nxt = 1
qui matrix list `b2'
local np = M_nbrf[1, `rf'] - 1
qui disp `np' " loadings at random effect " `rf' ", level " `lev'
if `np'>0{
matrix `s1' = `b2'[1,1..`np']
/*
* fudge: exponentiate s1
local ij = 1
while `ij'<=`np'{
matrix `s1'[1,`ij'] = exp(`s1'[1,`ij'])
local ij = `ij' + 1
}
* end fudge
*/
qui matrix list `s1'
matrix score double ${HG_s`rf'} = `s1'
local nxt = `nxt' + `np'
}
/* first (single non-) loading fixed at one, label in st. dev */
if M_frld[1,`rf']~=1{
matrix `s1' = (1)
local lab: colnames `b2'
local lab: word `nxt' of `lab'
matrix colnames `s1' = `lab'
capture drop `junk'
tempname junk
matrix score double `junk' = `s1'
if `np'>0{
qui replace ${HG_s`rf'} = ${HG_s`rf'} + `junk'
}
else{
matrix score double ${HG_s`rf'} = `s1'
*qui replace ${HG_s`rf'} = `junk'
}
}
qui disp "HG_s`rf' = ${HG_s`rf'} = " ${HG_s`rf'}[$which]
* extract standard deviation
* fudge: take exponential
* matrix `t'[`i',`i'] = exp(`b2'[1, `nxt'])
matrix `t'[`i',`i'] = `b2'[1, `nxt']
matrix CHmat[`rf'-1,`rf'-1]=`t'[`i',`i']
local nxt = `nxt' + 1
local i = `i' + 1
local rf = `rf' + 1
}
if (`nrf'>1&$HG_cor==1){ /* deal with correlations */
/* extract correlation parameters */
local i = 2
while (`i' <= `nrf'){
local k = `i' + `nrfold' - 1
local j = 1
while (`j' < `i'){
local l = `j' + `nrfold' - 1
qui disp "i = " `i' " j = " `j' " nxt = " `nxt'
matrix `t'[`i',`j'] = `b2'[1,`nxt']
matrix CHmat[`k',`l'] = `t'[`i',`j']
local j = `j' + 1
local nxt = `nxt' + 1
}
local i = `i' + 1
}
}
qui matrix list `t'
matrix M_chol = `t'
/* unpacked parameters */
local nrfold = M_nrfc[2,`lev']
local lev = `lev' + 1
} /* loop through levels */
}/*endelse HG_free */
if "`stop'"~=""{
exit
}
local nrfold = M_nrfc[2,1]
/* use B-matrix */
if $HG_tprf>1&$HG_bmat==1{
qui disp "dealing with B-matrix"
local i = 1
matrix Bmat = J($HG_tprf-1,$HG_tprf-1,0)
while `i'<$HG_tprf{
local j = 1
while `j' < $HG_tprf{
if M_b[`i',`j']>0{
matrix Bmat[`i',`j']=`b2'[1,`nxt']
local nxt = `nxt' + 1
}
local j = `j' + 1
}
local i = `i' + 1
}
qui matrix list Bmat
/* only works if B-matrix is upper diagonal */
local i=2
while `i'<$HG_tprf{
local k = `i' + `nrfold'
local j = 1
qui disp "making s`k'"
while `j'<`i'{
local l = `j' + `nrfold'
qui replace ${HG_s`k'} = ${HG_s`k'} + Bmat[`j',`i']*${HG_s`l'}
qui disp " adding Bmat[`j',`i']s`l'"
local j = `j' + 1
}
local i = `i' + 1
}
}
/* deal with geqs */
if $HG_ngeqs>0{
qui disp "dealing with geqs"
local i = 1
while `i'<=$HG_ngeqs{
local k = M_ngeqs[1,`i']
local n = M_ngeqs[2,`i']
qui disp "random effect `k' has `n' covariates"
local nxt2 = `nxt'+`n'-1
matrix `s1' = `b2'[1,`nxt'..`nxt2']
qui matrix list `s1'
local nxt = `nxt2' + 1
capture drop `junk'
matrix score double `junk' = `s1'
qui disp "multiply " `junk'[$which] " by HG_s`k' and add to HG_xb1"
qui replace $HG_xb1 = $HG_xb1 + `junk'*${HG_s`k'}
qui disp "HG_xb1:" $HG_xb1[$which]
local i = `i' + 1
}
}
/* use inter */
if $HG_inter~=0{
local k = $HG_l + 1
local j = $HG_r + 1
qui disp "HG_s`k' = HG_s`k'*HG_s`j'
qui replace ${HG_s`k'} = ${HG_s`k'}*${HG_s`j'}
}
/* use CHmat */
if $HG_free==0&$HG_tprf>1{
qui disp "dealing with Cholesky matrix"
qui matrix list CHmat
local i = 1
while (`i'<$HG_tprf){
local k = `i' + `nrfold'
qui replace `junk'=0
local j = `i'
qui disp "making s`k'"
while `j'<$HG_tprf{
local l = `j' + `nrfold'
qui replace `junk' = `junk' + CHmat[`j',`i']*${HG_s`l'}
qui disp " adding CHmat[`j',`i']s`l'"
local j = `j' + 1
}
qui replace ${HG_s`k'}=`junk'
qui disp "s`k' = ${HG_s`k'} = " ${HG_s`k'}[$which]
local i = `i' + 1
}
}
* label M_znow
local i=2
local lab
while `i'<=$HG_tprf{
local lab "`lab' ${HG_s`i'}"
local i = `i' + 1
}
matrix colnames M_znow=`lab'
qui disp "M_znow:"
qui matrix list M_znow
end

163
Modules/ado/plus/r/reop_ll.ado Normal file
View File

@ -0,0 +1,163 @@
*! Version 1.0.1 November 22 2000, by Guillaume R. Frechette (STB-59: sg158)
program define reop_ll
version 6.0
args todo b lnf g
tempvar theta1 F F1 F2 p db db1 db2 dr C
tempname rho s2su u x w gr
local nm1 = $S_n-1
local np1 = $S_n+1
local i = 1
while `i' < $S_n {
tempname _cut`i'
tempname g_cut`i'
tempname g_cut1`i'
tempname g_cut2`i'
local i = `i'+1
}
mleval `theta1' = `b', eq(1)
local h 0
local i 1
local j 2
while `i' < $S_n {
mleval `_cut`i'' = `b', eq(`j') scalar
if `h'>0 {
if `_cut`i''<=`_cut`h'' {
scalar `_cut`h''=`_cut`i''-0.01
/* the above correction should not be needed
in general, but it might preclude `bad'
estimates if, for instance, one performs a
random search */
}
}
local h = `h'+1
local i = `i'+1
local j = `j'+1
}
mleval `rho' = `b', eq(`np1') scalar
if `rho' >= 1 {
scalar `rho' = 0.99
di "rho >= 1, set to rho = 0.99"
}
scalar `s2su' = sqrt(2*`rho'/(1 - `rho'))
quietly {
gen double `F' = . in 1
gen double `F1' = . in 1
gen double `F2' = . in 1
by $S_i: gen double `p' = cond(_n==_N,0,.)
gen double `db' = 0
gen double `db1' = 0
gen double `db2' = 0
gen double `dr' = 0
gen double `C' = .
local m 1
while `m' <= $S_quad {
scalar `x' = $S_x[1,`m']
scalar `w' = $S_w[1,`m']
scalar `u' = `s2su'*`x'
local condf "cond($S_lhs==0, normprob(`_cut1'-`theta1'-`u')"
local i 1
local j 2
while `j' < $S_n {
local condf "`condf', cond($S_lhs==`i', normprob(`_cut`j''-`theta1'-`u')-normprob(`_cut`i''-`theta1'-`u')"
local i = `i'+1
local j = `j'+1
}
local condf "`condf', 1-normprob(`_cut`nm1''-`theta1'-`u'))"
local i 2
while `i' < $S_n {
local i = `i'+1
local condf "`condf')"
}
replace `C' = `condf'
replace `C' = 0.00000001 if `C' == 0
by $S_i: replace `F' = /*
*/ cond(_n==1,`C',`C'*`F'[_n-1])
replace `p' = `p' + `w'*`F'
local condg1 "cond($S_lhs==0, 0"
local condg2 "cond($S_lhs==0, -exp(-0.5*(`_cut1'-`theta1'-`u')^2)"
local h 1
local i 2
while `i' < $S_n {
local condg1 "`condg1', cond($S_lhs==`h', exp(-0.5*(`_cut`h''-`theta1'-`u')^2)"
local condg2 "`condg2', cond($S_lhs==`h', -exp(-0.5*(`_cut`i''-`theta1'-`u')^2)"
local h = `h'+1
local i = `i'+1
}
local condg1 "`condg1', exp(-0.5*(`_cut`nm1''-`theta1'-`u')^2))"
local condg2 "`condg2', 0)"
local i 2
while `i' < $S_n {
local i = `i'+1
local condg1 "`condg1')"
local condg2 "`condg2')"
}
replace `F1' = `F'
replace `F2' = `F'
by $S_i: replace `F1' = `condg1'*`F1'[_N]/`C'
by $S_i: replace `F2' = `condg2'*`F2'[_N]/`C'
by $S_i: replace `F' = (`condg1'+`condg2')*`F'[_N]/`C'
replace `db' = `db' + `w'*`F'
replace `db1' = `db1' + `w'*`F1'
replace `db2' = `db2' + `w'*`F2'
replace `dr' = `dr' + `w'*`u'*`F'
local m = `m' + 1
}
tempname lp grps
gen double `lp' = ln(`p'/sqrt(_pi))
by $S_i: gen byte `grps'=_n==_N
sum `grps' if $ML_samp
local N = r(sum)
sum `lp' if $ML_samp, meanonly
if r(N) !=`N' {
scalar `lnf' = .
exit
}
scalar `lnf' = r(sum)
if `todo'==0|`lnf'==. { exit }
by $S_i: replace `p' = `p'[_N]
replace `db' = `db'/(sqrt(2*_pi)*`p')
replace `db1' = -`db1'/(sqrt(2*_pi)*`p')
replace `db2' = -`db2'/(sqrt(2*_pi)*`p')
matrix vecaccum `g' = `db' $S_rhs, nocons
local i 0
local j 1
while `j' < $S_n {
capture {
matrix vecaccum `g_cut1`j'' = `db1' if $S_lhs == `j'
matrix vecaccum `g_cut2`j'' = `db2' if $S_lhs == `i'
}
if _rc {
scalar `lnf' = .
exit
}
matrix `g_cut`j'' = `g_cut1`j''+`g_cut2`j''
matrix `g' = `g', `g_cut`j''
local i = `i'+1
local j = `j'+1
}
replace `dr' = sum(`dr'/`p')
scalar `gr' = `dr'[_N]/(2*sqrt(2*_pi)*`rho'*(1-`rho'))
matrix `g' = `g', `gr'
}
end

160
Modules/ado/plus/r/reopc_ll.ado Normal file
View File

@ -0,0 +1,160 @@
*! Version 1.0.1 November 22 2000, by Guillaume R. Frechette (STB-59: sg158)
program define reopc_ll
version 6.0
args todo b lnf g
tempvar F F1 F2 p db1 db2 dr C
tempname rho s2su u x w gr
local nm1 = $S_n-1
local np1 = $S_n+1
local i = 1
while `i' < $S_n {
tempname _cut`i'
tempname g_cut`i'
tempname g_cut1`i'
tempname g_cut2`i'
local i = `i'+1
}
local h 0
local i 1
local j 2
while `i' < $S_n {
mleval `_cut`i'' = `b', eq(`i') scalar
if `h'>0 {
if `_cut`i''<=`_cut`h'' {
scalar `_cut`h''=`_cut`i''-0.01
/* the above correction should not be needed
in general, but it might preclude `bad'
estimates if, for instance, one performs a
random search */
}
}
local h = `h'+1
local i = `i'+1
local j = `j'+1
}
mleval `rho' = `b', eq($S_n) scalar
if `rho' >= 1 {
scalar `rho' = 0.99
di "rho >= 1, set to rho = 0.99"
}
scalar `s2su' = sqrt(2*`rho'/(1 - `rho'))
quietly {
gen double `F' = . in 1
gen double `F1' = . in 1
gen double `F2' = . in 1
by $S_i: gen double `p' = cond(_n==_N,0,.)
gen double `db1' = 0
gen double `db2' = 0
gen double `dr' = 0
gen double `C' = .
local m 1
while `m' <= $S_quad {
scalar `x' = $S_x[1,`m']
scalar `w' = $S_w[1,`m']
scalar `u' = `s2su'*`x'
local condf "cond($S_lhs==0, normprob(`_cut1'-`u')"
local i 1
local j 2
while `j' < $S_n {
local condf "`condf', cond($S_lhs==`i', normprob(`_cut`j''-`u')-normprob(`_cut`i''-`u')"
local i = `i'+1
local j = `j'+1
}
local condf "`condf', 1-normprob(`_cut`nm1''-`u'))"
local i 2
while `i' < $S_n {
local i = `i'+1
local condf "`condf')"
}
replace `C' = `condf'
replace `C' = 0.00000001 if `C' == 0
by $S_i: replace `F' = /*
*/ cond(_n==1,`C',`C'*`F'[_n-1])
replace `p' = `p' + `w'*`F'
local condg1 "cond($S_lhs==0, 0"
local condg2 "cond($S_lhs==0, -exp(-0.5*(`_cut1'-`u')^2)"
local h 1
local i 2
while `i' < $S_n {
local condg1 "`condg1', cond($S_lhs==`h', exp(-0.5*(`_cut`h''-`u')^2)"
local condg2 "`condg2', cond($S_lhs==`h', -exp(-0.5*(`_cut`i''-`u')^2)"
local h = `h'+1
local i = `i'+1
}
local condg1 "`condg1', exp(-0.5*(`_cut`nm1''-`u')^2))"
local condg2 "`condg2', 0)"
local i 2
while `i' < $S_n {
local i = `i'+1
local condg1 "`condg1')"
local condg2 "`condg2')"
}
replace `F1' = `F'
replace `F2' = `F'
by $S_i: replace `F1' = `condg1'*`F1'[_N]/`C'
by $S_i: replace `F2' = `condg2'*`F2'[_N]/`C'
by $S_i: replace `F' = (`condg1'+`condg2')*`F'[_N]/`C'
replace `db1' = `db1' + `w'*`F1'
replace `db2' = `db2' + `w'*`F2'
replace `dr' = `dr' + `w'*`u'*`F'
local m = `m' + 1
}
tempname lp grps
gen double `lp' = ln(`p'/sqrt(_pi))
by $S_i: gen byte `grps'=_n==_N
sum `grps' if $ML_samp
local N = r(sum)
sum `lp' if $ML_samp, meanonly
if r(N) !=`N' {
scalar `lnf' = .
exit
}
scalar `lnf' = r(sum)
if `todo'==0|`lnf'==. { exit }
by $S_i: replace `p' = `p'[_N]
replace `db1' = -`db1'/(sqrt(2*_pi)*`p')
replace `db2' = -`db2'/(sqrt(2*_pi)*`p')
local i 0
local j 1
while `j' < $S_n {
capture {
matrix vecaccum `g_cut1`j'' = `db1' if $S_lhs == `j'
matrix vecaccum `g_cut2`j'' = `db2' if $S_lhs == `i'
}
if _rc {
scalar `lnf' = .
exit
}
matrix `g_cut`j'' = `g_cut1`j''+`g_cut2`j''
if `j' == 1 {
matrix `g' = `g_cut`j''
}
else matrix `g' = `g', `g_cut`j''
local i = `i'+1
local j = `j'+1
}
replace `dr' = sum(`dr'/`p')
scalar `gr' = `dr'[_N]/(2*sqrt(2*_pi)*`rho'*(1-`rho'))
matrix `g' = `g', `gr'
}
end

98
Modules/ado/plus/r/reoprob.ado Normal file
View File

@ -0,0 +1,98 @@
*! Version 1.0.0 (11/22/00), G. R. Frechette (STB-59: sg158; STB-61: sg158.1)
program define reoprob, eclass
version 6.0
if replay() {
if "`e(cmd)'" ~= "reoprob" {
error 301
}
Replay `0'
}
else Estimate `0'
end
program define Estimate, eclass
syntax varlist [if] [in] [, I(varname) Quadrat(int 12) /*
*/ Level(passthru) *]
tempvar touse x w
marksample touse
markout `touse' `i'
tokenize `varlist'
local lhs "`1'"
macro shift 1
local rhs "`*'"
/* Get points and weights for Gauss-Hermite quadrature. */
ghquadm `quadrat' `x' `w'
/* Set up macros for ml function. */
global S_i "`i'"
global S_x "`x'"
global S_w "`w'"
global S_quad "`quadrat'"
global S_rhs "`rhs'"
global S_lhs "`lhs'"
/* get starting values */
tempname b0 s0
quietly oprobit `lhs' `rhs' if `touse'
mat `b0' = e(b)
mat `b0' = [`b0', 0.5]
quietly oprobit `lhs' if `touse'
mat `s0' = e(b)
mat `s0' = [`s0', 0.5]
/* number of categories */
quietly tab1 $S_lhs if `touse'
global S_n = _result(2)
/* create our version of `lhs' that runs from 0, ..., n-1 where
n is the number of categories */
tempvar dv
rename $S_lhs `dv'
quietly egen $S_lhs = group( `dv' ) if `touse'
quietly replace $_lhs = $_lhs-1
/* estimation equations */
local meqe "($S_lhs=`rhs', nocons)"
local start "(_cut1: $S_lhs=)"
local i = 1
while ( `i' < $S_n ) {
local meqe "`meqe' /_cut`i'"
local i = `i' + 1
}
local i = 2
while ( `i' < $S_n ) {
local start "`start' /_cut`i'"
local i = `i' + 1
}
local meqe "`meqe' /rho"
local start "`start' /rho"
/* Sort data. */
sort $S_i
/* optimization */
di in green _n "Fitting constant-only model:"
ml model d1 reopc_ll `start' if `touse', /*
*/ init(`s0', copy) maximize /*
*/ search(off) /*
*/ `options'
di in green _n "Fitting full model:"
ml model d1 reop_ll `meqe' if `touse', /*
*/ continue init(`b0', copy) maximize /*
*/ search(off) /*
*/ `options' /*
*/ title("Random Effects Ordered Probit")
estimate local cmd "reoprob"
Replay, `level'
end
program define Replay
syntax [, Level(int $S_level)]
ml display, level(`level')
end

95
Modules/ado/plus/r/reoprob.hlp Normal file
View File

@ -0,0 +1,95 @@
.-
help for ^reoprob^ (STB-59: sg158; STB-61: sg158.1)
.-
Random-effects ordered probit
-----------------------------
^reoprob^ depvar varlist [^if^ exp] [^in^ range] ^,^
[ ^i(^varname^)^ ^q^uadrat^(^#^)^ ^l^evel^(^#^)^ maximize_options ]
This command shares the features of all estimation commands; see help @est@.
To reset problem-size limits, see help @matsize@.
Description
-----------
^reoprob^ estimates a random-effects ordered probit model for panel datasets
using maximum likelihood estimation. The likelihood for each unit is
approximated by Gauss-Hermite quadrature.
Options
-------
^i(^varname^)^ specifies the variable corresponding to an independent unit
(e.g., a subject id). ^i(^varname^)^ is not optional.
^quadrat(^#^)^ specifies the number of points to use for Gaussian-Hermite
quadrature. It is optional, and the default is 12. Increasing this value
improves accuracy, but also increases computation time. Computation time
is roughly proportional to its value.
^level(^#^)^ specifies the confidence level, in percent, for confidence
intervals. The default is ^level(95)^ or as set by ^set level^.
maximize_options controls the maximization process and the display of
information; see [R] maximize. ^nolog^ suppresses the display of the
likelihood iterations. Use the ^trace^ option to view parameter
convergence. The ^ltol(^#^)^ and ^tol(^#^)^ option can be used to loosen
the convergence criterion (respectively 1e-7 and 1e-6 by default) during
specification searches. ^iter(^#^)^ specifies the maximum number of
iterations.
Examples
--------
. ^reoprob y x, i(id)^
. ^reoprob y x^
. ^reoprob y x, i(id) quad(24) nolog^
. ^reoprob y x, i(id) trace^
. ^reoprob^
Method
------
^reoprob^ uses the d1 method (analytic first derviatives) of Stata's ^ml^
commands. See Butler and Moffitt (1982) for details about using Gauss-Hermite
quadrature to approximate such integrals. Also see Green (2000) for
information on how to estimate a basic ordered probit model.
Author
------
Guillaume R. Frechette
Ohio State University
Department of Economics
410 Arps Hall
1945 North High Street
Columbus, OH 43210-1172
Tel: (614) 688-4140
Fax: (614) 292-4192
e-mail: frechette.6@@osu.edu
http://www.econ.ohio-state.edu/frechette/
Reference
---------
Butler, J.S. and R. Moffitt. 1982. A computationally efficient
quadrature procedure for the one-factor multinomial probit model.
Econometrica 50: 761-764.
Green, W. H. 2000. Econometric Analysis. Prentice Hall, New Jersey.
pp. 875-878.
Also see
--------
Manual: ^[R] xt, [R] xtprobit, [R] maximize, [R] oprobit^
On-line: help for @xt@, @xtreg@

48
Modules/ado/plus/r/replist.ado Normal file
View File

@ -0,0 +1,48 @@
program def replist, rclass
*! NJC 1.0.0 19 June 2000
version 6.0
gettoken list 0 : 0, parse(",")
if "`list'" == "" | "`list'" == "," {
di in r "nothing in list"
exit 198
}
* note that 0 copies => empty list
syntax , Copies(numlist int >=0) [ Block Global(str) Noisily ]
if length("`global'") > 8 {
di in r "global name must be <=8 characters"
exit 198
}
if `copies' == 1 {
local newlist "`list'"
}
else if `copies' > 1 {
if "`block'" != "" {
local c = 1
while `c' <= `copies' {
local newlist "`newlist'`list' "
local c = `c' + 1
}
}
else {
tokenize `list'
local n : word count `list'
local i = 1
while `i' <= `n' {
local c = 1
while `c' <= `copies' {
local newlist "`newlist'``i'' "
local c = `c' + 1
}
local i = `i' + 1
}
}
}
if "`noisily'" != "" { di "`newlist'" }
if "`global'" != "" { global `global' "`newlist'" }
return local list `newlist'
end

2
Modules/ado/plus/r/replist.hlp Normal file
View File

@ -0,0 +1,2 @@
.h listutil

31
Modules/ado/plus/r/revlist.ado Normal file
View File

@ -0,0 +1,31 @@
program def revlist, rclass
*! NJC 1.2.0 7 June 2000
* NJC 1.1.0 22 Dec 1999
* NJC 1.0.0 24 Oct 1999
version 6.0
gettoken list 0 : 0, parse(",")
if "`list'" == "" | "`list'" == "," {
di in r "nothing in list"
exit 198
}
syntax [, Noisily Global(str) ]
if length("`global'") > 8 {
di in r "global name must be <=8 characters"
exit 198
}
tokenize `list'
local nwords : word count `list'
local i = `nwords'
while `i' >= 1 {
local newlist "`newlist'``i'' "
local i = `i' - 1
}
if "`noisily'" != "" { di "`newlist'" }
if "`global'" != "" { global `global' "`newlist'" }
return local list `newlist'
end

2
Modules/ado/plus/r/revlist.hlp Normal file
View File

@ -0,0 +1,2 @@
.h listutil

129
Modules/ado/plus/r/rfpr_ll1.ado Normal file
View File

@ -0,0 +1,129 @@
*! version 1.0.0 7 June 1995 sg41: STB-26
program define rfpr_ll1 /* b log_likelihood grad */
version 4.0
local b "`1'"
local f "`2'"
local g "`3'"
macro shift 3
local options "FIRSTIT LASTIT FAST(string)"
parse "`*'"
local y = trim("$S_mldepn")
local doit "$S_sample"
local i "$S_ivar"
tempname rho s2su beta u
tempvar xb F p
local rhocol = colnumb(`b',"rho:_cons")
scalar `rho' = abs(`b'[1,`rhocol'])
if `rho' >= 1 {
scalar `rho' = 0.99
di "rho >= 1, set to rho = 0.99"
}
matrix `b'[1,`rhocol'] = `rho'
scalar `s2su' = sqrt(2*`rho'/(1 - `rho'))
matrix `beta' = `b'[1,"`y':"]
quietly {
matrix score double `xb' = `beta' if `doit'
gen double `F' = . in 1
by `doit' `i': gen double `p' = cond(_n==_N,0,.) if `doit'
/* Do computation this way if only log likelihood required. */
if "`fast'" == "0" {
local m 1
while `m' <= $S_quad {
scalar `u' = `s2su'*$S_x[`m']
#delimit ;
by `doit' `i': replace `F' =
cond(_n==1,
cond(`y', normprob(`xb' + `u'),
1 - normprob(`xb' + `u')),
cond(`y', normprob(`xb' + `u'),
1 - normprob(`xb' + `u'))
*`F'[_n-1])
if `doit' ;
#delimit cr
replace `p' = `p' + $S_w[`m']*`F' if `doit'
local m = `m' + 1
}
replace `F' = sum(log(`p'/sqrt(_pi))) if `doit'
scalar `f' = `F'[_N]
exit
}
/* Do computation this way if first derivatives required. */
tempname gr
tempvar db dr
gen double `db' = 0
gen double `dr' = 0
local m 1
while `m' <= $S_quad {
scalar `u' = `s2su'*$S_x[`m']
#delimit ;
by `doit' `i': replace `F' =
cond(_n==1,
cond(`y', normprob(`xb' + `u'),
1 - normprob(`xb' + `u')),
cond(`y', normprob(`xb' + `u'),
1 - normprob(`xb' + `u'))*`F'[_n-1])
if `doit' ;
replace `p' = `p' + $S_w[`m']*`F' if `doit' ;
by `doit' `i': replace `F' =
cond(`y',1,-1)*exp(-0.5*(`xb'+`u')^2)*`F'[_N]
/cond(`y',normprob(`xb'+`u'),1-normprob(`xb'+`u'))
if `doit' ;
#delimit cr
replace `db' = `db' + $S_w[`m']*`F' if `doit'
replace `dr' = `dr' + $S_w[`m']*`u'*`F' if `doit'
local m = `m' + 1
}
/* Compute log likelihood. */
replace `F' = sum(log(`p'/sqrt(_pi))) if `doit'
scalar `f' = `F'[_N]
/* Compute first derivatives. */
by `doit' `i': replace `p' = `p'[_N] if `doit'
replace `db' = `db'/(sqrt(2*_pi)*`p') if `doit'
local nvar = colsof(`beta') - 1
if `nvar' > 0 {
matrix `beta' = `beta'[1, 1..`nvar']
local vars : colnames(`beta')
}
matrix vecaccum `g' = `db' `vars' if `doit'
replace `dr' = sum(`dr'/`p') if `doit'
matrix `gr' = (0)
matrix `gr'[1,1] = `dr'[_N]/(2*sqrt(2*_pi)*`rho'*(1-`rho'))
matrix `g' = `g' , `gr'
}
end

233
Modules/ado/plus/r/rfprobit.ado Normal file
View File

@ -0,0 +1,233 @@
*! version 1.1.0 20jun1995 sg41: STB-26
program define rfprobit
version 4.0
local options "Level(integer $S_level)"
if substr("`1'",1,1)=="," | "`*'"=="" {
if "$S_E_cmd"~="rfprobit" {
error 301
}
parse "`*'"
}
else {
local varlist "req ex"
local if "opt"
local in "opt"
local options /*
*/ "`options' I(string) Quadrat(integer 6) noCHIsq noLOg *"
parse "`*'"
parse "`varlist'", parse(" ")
local y "`1'"
macro shift
xt_iis `i'
local i "$S_1"
tempvar doit x w
mark `doit' `if' `in'
markout `doit' `y' `*' `i'
/* Check to see if outcome varies. */
quietly count if `doit'
local n = _result(1)
quietly count if `y'==0 & `doit'
local n0 = _result(1)
if `n0'==0 | `n0'==`n' {
di _n in blu "outcome does not vary"
exit
}
/* Sort data. */
sort `doit' `i'
/* Get points and weights for Gaussian-Hermite quadrature. */
ghquad double(`x' `w'), n(`quadrat')
/* Set up macros for ml function. */
global S_sample "`doit'"
global S_ivar "`i'"
global S_x "`x'"
global S_w "`w'"
global S_quad "`quadrat'"
/* Fit constant-only model. */
if "`chisq'"=="" & "`*'"~="" {
rfpr_ml `y', title("Constant-only model") /*
*/ `log' `options'
local lf0 "lf0($S_1)"
local rho "rho($S_2)"
}
/* Fit full model. */
rfpr_ml `y' `*', title("Full model") post `lf0' `rho' /*
*/ `log' `options'
}
/* Display results. */
ml mlout rfprobit, level(`level')
/* Compute LR test for rho = 0. */
local chi = 2*($S_E_ll - $S_E_rho0)
#delimit ;
di in gr "LR test of rho = 0: chi2(" in ye "1" in gr ") = "
in ye %7.2f `chi' _n
in gr " Prob > chi2 = " in ye %7.4f
chiprob(1,`chi') ;
#delimit cr
end
program define rfpr_ml /* y x, title(string) POST noLOg ml_options */
version 4.0
local varlist "req ex"
local options "TITLE(string) POST LF0(string) RHO(string) noLOg *"
parse "`*'"
parse "`varlist'", parse(" ")
local y "`1'"
macro shift
local doit "$S_sample"
if "`lf0'"~="" { local lf0 "lf0(`lf0')" }
tempname llrho0 b0 b1 lllast b ll V
tempvar mldoit
/* Get initial values. */
quietly probit `y' `*' if `doit'
scalar `llrho0' = _result(2)
if "`log'"=="" {
di _n in gr "`title'"
di in gr "rho =" in ye %4.1f 0 /*
*/ in gr " Log Likelihood = " in ye `llrho0'
}
matrix `b0' = get(_b)
matrix coleq `b0' = `y'
matrix `b1' = (0)
matrix colnames `b1' = rho:_cons
matrix `b0' = `b0' , `b1'
local rcol = colnumb(`b0',"rho:_cons")
/* Search for good starting value for rho if not supplied. */
if "`rho'"=="" {
global S_mldepn "`y'"
scalar `lllast' = `llrho0'
local rho 0.05
local rhotry 0.1
while `rhotry' < 0.91 {
matrix `b0'[1,`rcol'] = `rhotry'
rfpr_ll1 `b0' `ll' `b' , fast(0)
if "`log'"=="" {
di in gr "rho =" in ye %4.1f `rhotry' /*
*/ in gr " Log Likelihood ~ " in ye `ll'
}
if `ll' < `lllast' { /* exit loop */
local rhotry 1
}
else {
scalar `lllast' = `ll'
local rho `rhotry'
local rhotry = `rhotry' + 0.1
}
}
}
matrix `b0'[1,`rcol'] = `rho'
/* Set up ml commands. */
ml begin
ml function rfpr_ll1
ml method deriv1
eq `y': `y' `*'
eq rho:
ml model `b' = `y' rho, depv(10) from(`b0')
ml sample `mldoit' if `doit', noauto
if "`log'"=="" { noisily ml max `ll' `V', `options' }
else quietly ml max `ll' `V', `options'
global S_1 = `ll'
global S_2 = `b'[1,`rcol']
if "`post'"~="" {
ml post rfprobit, `lf0' pr2 title("Random-Effects Probit")
}
global S_E_rho0 = `llrho0'
end
/*
Routines that compute weights and points for Gaussian-Hermite
quadrature follow:
*/
* version 1.0.1 29jun1995
program define ghquad
version 4.0
local varlist "req new min(2) max(2)"
local options "N(integer 10)"
parse "`*'"
parse "`varlist'", parse(" ")
local x "`1'"
local w "`2'"
if `n' + 2 > _N {
di in red /*
*/ "`n' + 2 observations needed to compute quadrature points"
exit 2001
}
tempname xx ww
local i 1
local m = int((`n' + 1)/2)
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*`x'[1] }
else if `i' == 4 { scalar `xx' = 1.91*`xx'-0.91*`x'[2] }
else { scalar `xx' = 2*`xx'-`x'[`i'-2] }
hermite `n' `xx' `ww'
qui replace `x' = `xx' in `i'
qui replace `w' = `ww' in `i'
local i = `i' + 1
}
if mod(`n', 2) == 1 { qui replace `x' = 0 in `m' }
qui replace `x' = -`x'[`n'+1-_n] in `i'/`n'
qui replace `w' = `w'[`n'+1-_n] in `i'/`n'
end
program define hermite /* integer n, scalar x, scalar w */
version 4.0
local n "`1'"
local x "`2'"
local w "`3'"
local last = `n' + 2
tempvar p
tempname i
qui gen double `p' = .
scalar `i' = 1
while `i' <= 10 {
qui replace `p' = 0 in 1
qui replace `p' = _pi^(-0.25) in 2
qui replace `p' = `x'*sqrt(2/(_n-2))*`p'[_n-1] /*
*/ - sqrt((_n-3)/(_n-2))*`p'[_n-2] in 3/`last'
scalar `w' = sqrt(2*`n')*`p'[`last'-1]
scalar `x' = `x' - `p'[`last']/`w'
if abs(`p'[`last']/`w') < 3e-14 {
scalar `w' = 2/(`w'*`w')
exit
}
scalar `i' = `i' + 1
}
di in red "hermite did not converge"
exit 499
end

102
Modules/ado/plus/r/rfprobit.hlp Normal file
View File

@ -0,0 +1,102 @@
.-
help for ^rfprobit^ (STB-26: sg41)
.-
Random-effects probit
---------------------
^rfprobit^ depvar [indepvars] [^if^ exp] [^in^ range] ^,^ [ ^i(^varname^)^
^q^uadrat^(^#^)^ ^nochi^sq ^nolo^g ^l^evel^(^#^)^ maximize_options ]
This command shares the features of all estimation commands; see help @est@.
To reset problem-size limits, see help @matsize@.
Description
-----------
^rfprobit^ estimates a random-effects probit model for cross-sectional time-
series data sets using maximum likelihood estimation. The likelihood (for
an independent unit i) is expressed as an integral which is computed using
Gaussian-Hermite quadrature. This computational procedure is only accurate
when there are a small-to-moderate number of time periods T_i per unit i.
It is recommended that ^rfprobit^ only be used when max(T_i) <= 50.
Options
-------
^i(^varname^)^ specifies the variable corresponding to an independent unit
(e.g., a subject id). This variable represents the i in x_it. Either
this option must be specified or i must be set using the ^iis^ command;
see help for @xt@.
^quadrat(^#^)^ specifies the number of points to use for Gaussian-Hermite
quadrature. Default is 6. Increasing this value slightly improves
accuracy, but also increases computation time. Computation time is
roughly proportional to its value.
^nochisq^ omits the estimation of the constant-only model. This will reduce
computation time at the cost of not being able to calculate the model
chi-squared or pseudo R^^2.
^nolog^ suppress the display of the likelihood iterations.
^level(^#^)^ specifies the significance level, in percent, for confidence
intervals of the coefficients; see help @level@.
maximize_options control the maximization process; see [7] maximize.
Use the ^trace^ option to view parameter convergence.
The ^ltol(^#^)^ option can be used to loosen the convergence
criterion (default is 1e-6) during specification searches.
Examples
--------
. ^rfprobit y x, i(id)^
. ^iis id^
. ^rfprobit y x^
. ^rfprobit y x, i(id) nochisq^
. ^rfprobit y x, i(id) quad(8) nolog^
. ^rfprobit y x, i(id) trace^
. ^rfprobit^
Method
------
^rfprobit^ uses the deriv1 method (analytic first derviatives) of Stata's ^ml^
commands. See Butler and Moffitt (1982) for details.
Author
------
Bill Sribney
Stata Corporation
702 University Drive East
College Station, TX 77840
Phone: 409-696-4600
800-782-8272
Fax: 409-696-4601
email: tech_support@@stata.com
Reference
---------
Butler, J.S. and R. Moffitt. 1982. A computationally efficient quadrature
procedure for the one-factor multinomial probit model. Econometrica 50:
761-764.
Also see
--------
STB: STB-26 sg41
Manual: [5s] xt, [5s] xtreg, [7] maximize
On-line: help for @xt@, @xtreg@

44
Modules/ado/plus/r/rotlist.ado Normal file
View File

@ -0,0 +1,44 @@
program def rotlist, rclass
*! NJC 1.1.0 6 June 2000
* NJC 1.0.0 26 Apr 2000
version 6.0
gettoken list 0 : 0, parse(",")
if "`list'" == "" | "`list'" == "," {
di in r "nothing in list"
exit 198
}
syntax , Rot(int) [ Global(str) Noisily ]
if length("`global'") > 8 {
di in r "global name must be <=8 characters"
exit 198
}
local n : word count `list'
local rot = mod(`rot', `n')
if `rot' == 0 {
local newlist "`list'"
}
else {
tokenize `list'
local i = `rot' + 1
while `i' <= `n' {
local newlist "`newlist'``i'' "
local i = `i' + 1
}
local i = 1
while `i' <= `rot' {
local newlist "`newlist'``i'' "
local i = `i' + 1
}
}
if "`noisily'" != "" { di "`newlist'" }
if "`global'" != "" { global `global' "`newlist'" }
return local list `newlist'
end

2
Modules/ado/plus/r/rotlist.hlp Normal file
View File

@ -0,0 +1,2 @@
.h listutil