*! runmlwin.ado, George Leckie and Chris Charlton, 17Jul2025
program runmlwin, eclass
if c(stata_version) >= 15 local user user
if _caller() >= 12 version 12.0, `user'
if _caller() <= 9 version 9.0, `user'
display " " // display an empty line between the model command and the model output
if replay() { // replay the results
if ("`e(cmd)'" ~= "runmlwin") error 301
syntax [, Level(cilevel) CFORMAT(string) PFORMAT(string) SFORMAT(string) noHEADer noGRoup noCONTrast noFETable noRETable SD CORrelations OR IRr RRr MOde MEdian Zratio *]
runmlwin_display, level(`level') cformat(`cformat') pformat(`pformat') sformat(`sformat') `header' `group' `contrast' `fetable' `retable' `sd' `correlations' `or' `irr' `rrr' `mode' `median' `zratio' // display the results
makecns, displaycns
}
else { // fit the model
syntax anything [if] [in], [Level(cilevel) CFORMAT(string) PFORMAT(string) SFORMAT(string) noHEADer noGRoup noCONTrast noFETable noRETable SD CORrelations OR IRr RRr MOde MEdian Zratio *]
timer clear 99
timer on 99
Estimates `0'
if "`e(mcmcnofit)'" == "1" { // If the user specifies the nofit option in MCMC then runmlwin needs to exit, but without error
exit
}
timer off 99
quietly timer list
ereturn scalar time = r(t99)
timer clear 99
runmlwin_display, level(`level') cformat(`cformat') pformat(`pformat') sformat(`sformat') `header' `group' `contrast' `fetable' `retable' `sd' `correlations' `or' `irr' `rrr' `mode' `median' `zratio' // display the results
}
end
program Estimates, eclass sortpreserve
if _caller() >= 12 version 12.0
if _caller() <= 9 version 9.0
* Locate and record highest level
local runmlwin_cmdline `0'
* Separate off options
local maxlevels 1
gettoken comma tmpstr : 0, parse(",") bind
while strpos(`"`tmpstr'"', "level") != 0 {
local tmpstr = substr(`"`tmpstr'"', strpos(`"`tmpstr'"', "level") + 5, .)
local lev = substr(`"`tmpstr'"', 1, strpos(`"`tmpstr'"', "(") - 1)
* Check this really is a level option
if `=real("`lev'")' != . {
if `lev' > `maxlevels' {
local maxlevels `lev'
}
}
}
if ("`verbose'"~="") display as text "Highest level: `maxlevels'"
local levargs
forvalues l = 2/`maxlevels' {
local levargs LEVEL`l'(string) `levargs'
}
#delimit ;
syntax anything(name=eqlist id="equations" equalok) [if] [in],
[
`levargs' // min=2 as first var is id, second var is, e.g., cons
]
LEVEL1(string) // min = 1 as in discrete response models no variables are specified at level 1. Also in multivariate response models there are no covariates with random effects at level 1.
[
Weights(string)
Constraints(numlist >0 integer)
IGLS RIGLS TOLerance(numlist >0 integer min=1 max=1) MAXIterations(numlist >0 integer min=1 max=1)
FPSandwich RPSandwich
INITSPrevious INITSB(namelist min=1 max=1) INITSV(namelist min=1 max=1) INITSModel(namelist min=1 max=1) // Initial values matrix
DISCRETE(string)
MCMC(string)
SEED(numlist integer min=1 max=1) SIMulate(namelist min=1 max=1)
MLWINPATH(string) MLWINSCRIPTPATH(string)
VIEWMacro SAVEMacro(string) SAVEWorksheet(string) SAVEStata(string)
USEWorksheet(string)
Level(cilevel) CFORMAT(string) PFORMAT(string) SFORMAT(string) OR IRr RRr SD CORrelations MOde MEdian Zratio noHEADer noGRoup noCONTrast noFETable noRETable
noDrop FORCESort FORCERecast noMLWIN noPause noVERSIONCHECK BATCh noSORT
PLUGIN Verbose VIEWFULLMacro SAVEFULLMacro(string) SAVEEQuation(string)
MLWINSETTINGS(string)
] ;
#delimit cr
local doublevar 1
if c(mode) == "batch" local batch = c(mode)
if "`mlwinpath'" == "" & "`mlwinscriptpath'" ~= "" & "`batch'" ~= "" local mlwinpath `mlwinscriptpath'
if "`mlwinpath'" == "" & "$MLwiNScript_path" ~= "" & "`batch'" ~= "" local mlwinpath $MLwiNScript_path
if "`mlwinpath'" == "" & "$MLwiN_path" ~= "" local mlwinpath $MLwiN_path
// will need to allow for fact that users store MLwiN in different locations. User can put global in profile.do to allow this
// slightly unfortunate that the installation folder changes with each update from 2.20 to 2.21 etc at user has to keep on changing this to run this command
if "`mlwin'"~="nomlwin" & "`mlwinpath'" ~= "" {
capture confirm file "`mlwinpath'"
if _rc == 601 {
display as error "`mlwinpath' does not exist." _n
exit 198
}
if "`versioncheck'" ~= "noversioncheck" {
quietly capture runmlwin_verinfo `mlwinpath'
if _rc == 198 {
display as error "`mlwinpath' is not a valid version of MLwiN"
exit 198
}
local majorver `r(ver1)'
local minorver : display %02.0f `r(ver2)'
local versionok = 1
local versionold = 0
if (`majorver' < 2) | (`majorver' == 2 & `minorver' < 36) local versionok = 0
if (`majorver' < 3) | (`majorver' == 3 & `minorver' < 14) local versionold = 1
if `versionok' == 0 {
display as error "runmlwin assumes MLwiN version 2.36 or higher. You can download the latest version of MLwiN at:" _n "https://www.bristol.ac.uk/cmm/software/mlwin/download/upgrades.html." _n "If you want to ignore this warning and attempt to continue anyway you can use the noversioncheck option"
exit 198
}
if `versionold' == 1 display as error "WARNING: Your version of MLwiN is out of date. You can download the latest version of MLwiN at:" _n "https://www.bristol.ac.uk/cmm/software/mlwin/download/upgrades.html"
}
local mlwinversion `majorver'.`minorver'
capture confirm number `mlwinversion' // Check version number is valid
if _rc { // if not replace with missing
local mlwinversion .
}
if `mlwinversion' < 3 { // Checking for double->single rounding could be removed if we stopped supporting MLwiN prior to v3.00
local doublevar 0
}
}
* Mark the sample (Defines a 0/1 to-use variable that records which observations are to be used in subsequent code)
marksample touse, novarlist // This gives the full sample which is sent to MLwiN. Not a listwise delete sample.
* Discrete syntax
if ("`discrete'"~="") {
local 0 , `discrete'
syntax , Distribution(string) [Link(namelist min=1 max=1) DEnominator(varlist numeric) Extra Offset(varname numeric) Proportion(varname) Basecategory(numlist integer min=1 max=1) MQL1 MQL2 PQL1 PQL2]
local validdistributions normal binomial poisson nbinomial multinomial
local checkdistribution :list distribution & validdistributions
if "`checkdistribution'"=="" {
display as error "Invalid distribution(). Valid distributions are: normal, binomial, poisson, nbinomial, multinomial"
exit 198
}
* Method of linearization
if ("`mql2'"=="" & "`pql1'"=="" & "`pql2'"=="") local linearization MQL1
if ("`mql2'"~="") local linearization MQL2
if ("`pql1'"~="") local linearization PQL1
if ("`pql2'"~="") local linearization PQL2
* Check that the link function is correctly specified.
if "`link'"=="" {
if ("`distribution'"=="binomial") {
di as error "You must specify the link() function. Valid link functions for the binomial distribution are: logit, probit and cloglog."
exit 198
}
if ("`distribution'"=="multinomial") {
di as error "You must specify the link() function. Valid link functions for the multinomial distribution are: mlogit, ologit, oprobit, ocloglog."
exit 198
}
if ("`distribution'"=="poisson" | "`distribution'"=="nbinomial") local link log
}
if "`link'" ~= "" {
local 0 , `link'
syntax , [Identity Logit Probit Cloglog Mlogit OLogit OProbit OCloglog LOG]
if ("`distribution'"=="binomial") {
if ~inlist("`link'","logit","probit","cloglog") {
display as error "Invalid link() function. Valid link functions for the binomial distribution are: logit, probit and cloglog." _n
exit 198
}
}
if ("`distribution'"=="multinomial") {
if ~inlist("`link'","mlogit","ologit","oprobit","ocloglog") {
display as error "Invalid link() function. Valid link functions for the multinomial distribution are: mlogit, ologit, oprobit, ocloglog." _n
exit 198
}
}
if ("`distribution'"=="poisson") {
if ~inlist("`link'","log") {
display as error "Invalid link() function. Valid link functions for the poisson distribution are: log." _n
exit 198
}
}
if "`identity'" ~= "" local link identity
if "`logit'" ~= "" local link logit
if "`probit'" ~= "" local link probit
if "`cloglog'" ~= "" local link cloglog
if "`mlogit'" ~= "" local link mlogit
if "`ologit'" ~= "" local link ologit
if "`oprobit'" ~= "" local link oprobit
if "`ocloglog'" ~= "" local link ocloglog
if "`log'" ~= "" local link log
}
}
else {
local link identity
}
******************************************************************************
* (1A) PARSE RESPONSE AND FIXED PART PREDICTOR LIST
******************************************************************************
tempname o
.`o' = ._eqlist.new, eqopts(NOConstant) eqargopts(EQ) numdepvars(1) noneedvarlist // For properties of this class see C:\Program Files (x86)\Stata11\ado\base\_eqspec.class
.`o'.parse `eqlist'
local numstataeqns = `.`o'.eq count'
if ("`verbose'"~="") display as text "number of eqns: `numstataeqns'"
* Derive model type
if (`numstataeqns'==1 & "`distribution'" ~= "multinomial") local mtype univariate
if ("`distribution'" == "multinomial") local mtype multinomial
if (`numstataeqns'>1 & "`distribution'" ~= "multinomial") local mtype multivariate
assert"`mtype'"~=""
local response
local fp
if "`mtype'"=="multivariate" {
// Count the number of responses
local nummlwineqns = 0
local valideqns
forvalues e = 1/`numstataeqns' {
local 0 , `.`o'.eq options `e''
syntax [, EQuation(numlist integer)]
local eq `equation'
if `:list sizeof eq' == 1 {
local valideqns `valideqns' `eq'
local name`eq' `.`o'.eq name `e''
local response`eq' `.`o'.eq depvars `e''
// Add to complete list
local response `response' `response`eq''
local ++nummlwineqns
}
}
* Reserve the first nummlwineqns for separate coefficients
local numfpbrackets = `nummlwineqns'
forvalues e = 1/`numstataeqns' {
local 0 , `.`o'.eq options `e''
syntax [, EQuation(numlist integer)]
local eq `equation'
if "`eq'" == "" {
display as error "equation option not set for `e'"
exit 198
}
if `:list sizeof eq' == 1 {
if ("`verbose'"~="") display as text "Separate coefficients"
local fp`eq' `.`o'.eq indepvars `e''
local numfpvars`eq' :list sizeof fp`e'
local dups :list dups fp`eq'
if "`dups'" ~= "" {
display as error "`dups' are duplicated in the fixed part equation `eq'"
exit 198
}
foreach var of local fp`eq' {
local fpname`eq' `fpname`eq'' `var'.`eq'
}
local pat
forvalues a = 1/`nummlwineqns' {
if `:list a in eq' local pat `pat' 1
else local pat `pat' 0
}
local rpat`eq' `pat'
}
else {
if ("`verbose'"~="") display as text "Common coefficients"
// NOTE: This may give unnecessary errors if the common equation comes before a corresponding separate one
if "`:list eq - valideqns'" ~= "" {
display as error "invalid equation(s) specified"
exit 198
}
local ++numfpbrackets
local fp`numfpbrackets' `.`o'.eq depvars `e'' `.`o'.eq indepvars `e''
local numfpvars`numfpbrackets' :list sizeof fp`e'
local dups :list dups fp`numfpbrackets'
if "`dups'" ~= "" {
display as error "`dups' are duplicated in the fixed part equation `eq'"
exit 198
}
foreach var of local fp`e' {
local fpname`numfpbrackets' `fpname`numfpbrackets'' `var'.`=subinstr("`eq'", " ", "", .)'
}
local pat
forvalues a = 1/`nummlwineqns' {
if `:list a in eq' local pat `pat' 1
else local pat `pat' 0
}
local rpat`numfpbrackets' `pat'
}
}
forvalues e = 1/`numfpbrackets' {
// Add to complete list
local fp `fp' `fp`e''
}
local numfpvars :list sizeof fp
}
* Multinomial response
if "`mtype'"=="multinomial" {
local response1 `.`o'.eq depvars 1'
capture tab `response1'
local nummlwineqns `=`r(r)' - 1'
local numfpbrackets = `nummlwineqns'
local valideqns
forvalues e = 1/`nummlwineqns' {
local valideqns `valideqns' `e'
}
if ("`verbose'"~="") display as text "Valid Equations: `valideqns'"
forvalues e = 2/`numstataeqns' {
local 0 , `.`o'.eq options `e''
syntax [, CONtrast(numlist integer)]
local cat `contrast'
if ("`verbose'"~="") display as text "`cat'"
if "`cat'" == "" {
display as error "cat option not set for `e'"
exit 198
}
if "`:list cat - valideqns'" ~= "" {
display as error "Invalid contrast(s) specified. The contrast(s) specified in contrasts() must tally with the contrast(s) in the model, in this case: 1,2,...,`nummlwineqns'" _n
exit 198
}
if `:list sizeof cat' == 1 {
if ("`verbose'"~="") display as text "Separate coefficients"
local fp`cat' `.`o'.eq depvars `e'' `.`o'.eq indepvars `e''
local dups :list dups fp`cat'
if "`dups'" ~= "" {
display as error "`dups' are duplicated in the fixed part contrast `contrast'"
exit 198
}
}
else {
if ("`verbose'"~="") display as text "Common coefficients"
local ++numfpbrackets
local fp`numfpbrackets' `.`o'.eq depvars `e'' `.`o'.eq indepvars `e''
local numfpvars`numfpbrackets' :list sizeof fp`numfpbrackets'
local fpname`numfpbrackets'
local dups :list dups fp`numfpbrackets'
if "`dups'" ~= "" {
display as error "`dups' are duplicated in the fixed part contrast `contrast'"
exit 198
}
capture levelsof `response1', local(responsecats)
local catvals
local respnum = 1
foreach catn of local responsecats {
if "`catn'" ~= "`basecategory'" {
if `:list respnum in cat' local catvals `catvals' `catn'
local ++respnum
}
}
foreach var of local fp`numfpbrackets' {
local fpname`numfpbrackets' `fpname`numfpbrackets'' `var'.`=subinstr("`catvals'", " ", "", .)'
}
local pat
forvalues a = 1/`nummlwineqns' {
if `:list a in cat' local pat `pat' 1
else local pat `pat' 0
}
local rpat`numfpbrackets' `pat'
}
}
// Check there are at least nummlwineqns equations
capture levelsof `response1', local(responsecats)
local e = 1
foreach cat of local responsecats {
if "`cat'" ~= "`basecategory'" {
// Add to complete list
local response`e' `response1'
local response `response' `response`e''
local fp`e' `.`o'.eq indepvars 1' `fp`e''
local numfpvars`e' :list sizeof fp`e'
local dups :list dups fp`e'
if "`dups'" ~= "" {
display as error "`dups' are duplicated in the fixed part"
exit 198
}
local fpname`e'
foreach var of local fp`e' {
local fpname`e' `fpname`e'' `var'.`cat'
}
local pat
forvalues a = 1/`nummlwineqns' {
if `:list a in e' local pat `pat' 1
else local pat `pat' 0
}
local rpat`e' `pat'
local ++e
}
}
forvalues e = 1/`numfpbrackets' {
// Add to complete list
local fp `fp' `fp`e''
}
local numfpvars :list sizeof fp
}
if "`mtype'" == "univariate" {
local numfpbrackets = `numstataeqns'
local nummlwineqns = `numstataeqns'
* Parse each equation
forvalues e = 1/`numfpbrackets' {
local name`e' `.`o'.eq name `e''
local response`e' `.`o'.eq depvars `e''
if _caller() >= 11 {
local 0 `.`o'.eq indepvars `e''
syntax varlist(fv)
fvexpand `varlist'
local fpname`e' `r(varlist)'
fvrevar `r(varlist)'
local fp`e' `r(varlist)'
local newvars :list fp`e' - varlist
capture _rmcoll `fp`e'', noconstant
local fp`e' `r(varlist)'
if "`newvars'" ~= "" {
quietly compress `newvars'
}
local i = 1
foreach x in `fp`e'' {
local xname :word `i' of `fpname`e''
if substr("`x'", 1, 2) == "o." {
local fp`e' : list fp`e' - x
local fpname`e' : list fpname`e' - xname
display as text "note: `xname' omitted because of collinearity"
}
else {
if "`:type `x''" == "double" & `:list x in newvars' & `doublevar' == 0{
quietly recast float `x', force
display as error "`xname' has more precision that MLwiN can handle, forcing to float"
}
local ++i
}
}
}
else {
local fpname`e' `.`o'.eq indepvars `e''
local fp`e' `.`o'.eq indepvars `e''
}
local dups :list dups fp`e'
if "`dups'" ~= "" {
display as error "`dups' are duplicated in the fixed part"
exit 198
}
local numfpvars`e' :list sizeof fp`e'
// Add to complete list
local response `response' `response`e''
local fp `fp' `fp`e''
}
local numfpvars :list sizeof fp
}
local fpname
forvalues e = 1/`numfpbrackets' {
local fpname `fpname' `fpname`e''
}
forvalues e = 1/`numfpbrackets' {
if ("`=word("`distribution'", `e')'"=="binomial") {
if ("`extra'" ~= "") {
capture assert (`denominator' == 1) & (`response' == 1 | `response' == 0)
if ~_rc {
display as error "Extra option applies to binomial responses, not binary"
exit 198
}
}
}
}
if ("`verbose'"~="") {
* Check that the syntax has been parsed properly
display _n as result "MODEL " as text "Num. of responses: `numfpbrackets' " as text "Model type: `mtype'" _n
forvalues e = 1/`numfpbrackets' {
display as result "FP EQUATION `e' " as text "Name: `name`e'', Response: `response`e'', Predictors: `fp`e'', Num. of predictors: `numfpvars`e'', Options: `opts`e''" _n
}
}
******************************************************************************
* (1B) PARSE LEVELS
******************************************************************************
* Level parsing (Level X identifier, random part covariates and options)
local residualsall
local factorson = 0
* Total number of random part variables (not total number of random part parameters)
local numrpvars = 0
local emptylevels
forvalues l = `maxlevels'(-1)1 {
if "`level`l''"~="" {
gettoken lev`l'id 0 : level`l', parse(":") // get everything before colon and put in lev`l'id, put rest of string in 0
if `:list sizeof lev`l'id' > 1 { // if you have specified more than one variable before the colon then issue an error message
display as error "Only one level `l' ID is allowed; a colon should appear after the level `l' ID and before the variable list" _n
exit 198
}
gettoken colon rp`l': 0, parse(":") // Put : in colon and the variable list and any options in rp`l'
if strpos("`rp`l''", ",") > 0 { // If there is a comma is the variable list and option list (i.e. options have been specified after a comma as expected) then ...
gettoken rp`l' options : rp`l', parse(",") bind // put the true var list in rp`l' and put the comma and the options in options
if "`rp`l''" == "," local rp`l' "" // No RP variables specified
else gettoken comma options: options, parse(",") // Put , in comma and the true option list in options
}
local 0 `rp`l'', lev`l'id(`lev`l'id') `options' // Define new local called 0 which has contents in "". We can prob remove `varlist' but need to check.
if "`l'" ~= "1" { // If level is level 2,3,4,5,...
// syntax interpret whatever is in `0'
#delimit ;
syntax [anything(name=rpeqlist id="equations" equalok)],
[LEV`l'id(varname)]
[Diagonal]
[ELements(namelist min=1 max=1)]
[DESIGN(string)]
[RESET(namelist min=1 max=1)]
[Residuals(string)]
[Weightvar(varname)]
[MMIds(varlist)]
[MMWeights(varlist)]
[CARIds(varlist)]
[CARWeights(varlist)]
[PAREXpansion]
[FLInits(string)]
[FLConstraints(string)]
[FVInits(string)]
[FVConstraints(string)]
[FScores(string)];
#delimit cr
}
else { // If level is level 1
#delimit ;
syntax [anything(name=rpeqlist id="equations" equalok)],
LEV`l'id(varname)
[Diagonal]
[ELements(namelist min=1 max=1)]
[DESIGN(string)]
[RESET(namelist min=1 max=1)]
[Residuals(string)]
[Weightvar(varname)]
[MMIds(varlist)]
[MMWeights(varlist)]
[CARIds(varlist)]
[CARWeights(varlist)]
[FLInits(string)]
[FLConstraints(string)]
[FVInits(string)]
[FVConstraints(string)]
[FScores(string)];
#delimit cr
}
* Parse the random part equations and variable specification
tempname o
.`o' = ._eqlist.new, eqopts(NOConstant) eqargopts(EQ) numdepvars(1) noneedvarlist // For properties of this class see C:\Program Files (x86)\Stata11\ado\base\_eqspec.class
.`o'.parse `rpeqlist' // capture because in some models we may choose to include no predictors in the random part at a particular level (e.g. a bivariate normal response model fitted as a univariate response model has no random part at level 1)
local numstataeqns`l' = `.`o'.eq count'
local numrpbrackets`l' = `nummlwineqns'
if ("`verbose'"~="") display as text "Valid Equations: `valideqns'"
forvalues e = 1/`numstataeqns`l'' { // loop over the equations
local 0 , `.`o'.eq options `e''
syntax [, EQuation(numlist integer) CONtrast(numlist integer) *]
if ("`contrast'" ~= "") local eq `contrast'
else local eq `equation'
if ("`verbose'"~="") display as text "`eq'"
if "`mtype'" == "univariate" {
local rp`l'_`e' `.`o'.eq depvars `e'' `.`o'.eq indepvars `e''
local rp`l'name`e' `rp`l'_`e''
local dups :list dups rp`l'_`e'
if "`dups'" ~= "" {
display as error "`dups' are duplicated in the random part level `l' equation `e'"
exit 198
}
}
else {
if "`eq'" == "" {
if `e' == 1 {
forvalues a = 1/`nummlwineqns' {
local name`a' `.`o'.eq name `e''
local rp`l'_`a' `.`o'.eq depvars `e'' `.`o'.eq indepvars `e''
local rpeq`l'_`a' `a'
local dups :list dups rp`l'_`a'
if "`dups'" ~= "" {
display as error "`dups' are duplicated in the random part level `l' contrast `a'"
exit 198
}
}
}
else {
display as error "equation option not set for `e'"
exit 198
}
}
else {
if "`:list eq - valideqns'" ~= "" {
display as error "invalid contrast(s) or equation(s) specified"
exit 198
}
if `:list sizeof eq' == 1 {
// Separate
local name`eq' `name`eq'' `.`o'.eq name `e''
local rp`l'_`eq' `rp`l'_`eq'' `.`o'.eq depvars `e'' `.`o'.eq indepvars `e''
local rpeq`l'_`eq' `eq'
local dups :list dups rp`l'_`eq'
if "`dups'" ~= "" {
display as error "`dups' are duplicated in the random part level `l' equation `eq'"
exit 198
}
}
else {
// Common
local ++numrpbrackets`l'
local name`numrpbrackets`l'' `.`o'.eq name `e''
local rp`l'_`numrpbrackets`l'' `.`o'.eq depvars `e'' `.`o'.eq indepvars `e''
local numrpvars`numrpbrackets`l'' :list sizeof rp`l'_`numrpbrackets`l''
local rpeq`l'_`numrpbrackets`l'' `eq'
local rp`l'name`numrpbrackets`l''
local dups :list dups rp`l'_`numrpbrackets`l''
if "`dups'" ~= "" {
display as error "`dups' are duplicated in the random part level `l' equation `eq'"
exit 198
}
if "`mtype'" == "multivariate" {
foreach var of local rp`l'_`numrpbrackets`l'' {
local rp`l'name`numrpbrackets`l'' `rp`l'name`numrpbrackets`l'' ' `var'.`=subinstr("`eq'", " ", "", .)'
}
}
if "`mtype'" == "multinomial" {
capture levelsof `response1', local(responsecats)
local eqvals
local respnum = 1
foreach cat of local responsecats {
if "`cat'" ~= "`basecategory'" {
if `:list respnum in eq' local eqvals `eqvals' `cat'
local ++respnum
}
}
foreach var of local rp`l'_`numrpbrackets`l'' {
local rp`l'name`numrpbrackets`l'' `rp`l'name`numrpbrackets`l'' ' `var'.`=subinstr("`eqvals'", " ", "", .)'
}
}
local pat
forvalues a = 1/`nummlwineqns' {
if `:list a in eq' local pat `pat' 1
else local pat `pat' 0
}
local rppat`l'_`numrpbrackets`l'' `pat'
}
}
}
}
if "`mtype'" == "multivariate" {
forvalues eq = 1/`nummlwineqns' {
if `l' == 1 & "`discrete'" ~= "" & "`=word("`distribution'", `eq')'"=="binomial" local rp`l'_`eq' bcons
local numrpvars`eq' :list sizeof rp`l'_`eq'
local rp`l'name`eq'
foreach var of local rp`l'_`eq' {
local rp`l'name`eq' `rp`l'name`eq'' `var'.`=subinstr("`eq'", " ", "", .)'
}
}
}
if "`mtype'" == "multinomial" {
capture levelsof `response1', local(responsecats)
local eq = 1
foreach cat of local responsecats {
if "`cat'" ~= "`basecategory'" {
local numrpvars`eq' :list sizeof rp`l'_`eq'
local rp`l'name`eq'
foreach var of local rp`l'_`eq' {
local rp`l'name`eq' `rp`l'name`eq'' `var'.`=subinstr("`cat'", " ", "", .)'
}
local ++eq
}
}
}
if "`mtype'" == "multivariate" || "`mtype'" == "multinomial" {
forvalues eq = 1/`nummlwineqns' {
local pat
forvalues a = 1/`nummlwineqns' {
if `:list a in eq' local pat `pat' 1
else local pat `pat' 0
}
local rppat`l'_`eq' `pat'
}
}
local rp`l'
forvalues e = 1/`numrpbrackets`l'' { // loop over the equations
local rp`l' `rp`l'' `rp`l'_`e''
if "`mtype'" == "univariate" local rp`l'name`e' `rp`l'_`e''
}
forvalues e = 1/`numrpbrackets`l'' { // loop over the equations
* Ensure random part names are in the correct order
local origrp`l'name `origrp`l'name ' `rp`l'name`e''
local tmpname `rp`l'name`e''
local rp`l'name `rp`l'name' `:list fpname & tmpname'
local tmpname `:list tmpname - fpname'
forvalues lev = `maxlevels'(-1)`=`e'+1' {
local rp`l'name `rp`l'name' `:list rp`lev'name & tmpname'
local tmpname `:list tmpname - rp`lev'name'
}
local rp`l'name `rp`l'name' `tmpname'
}
// NOTE: The following lines probably aren't needed
if "`numfpbrackets'" == "1"{
local rp`l'_1 `.`o'.eq depvars 1' `.`o'.eq indepvars 1'
local rp`l' `rp`l'_1'
}
if ("`verbose'"~="") {
display as result "LEVEL `l'"
forvalues e = 1/`numfpbrackets' {
display as result "RP EQUATION `e', " as text "Name: `name`e'', " as text "Predictors: `rp`l'_`e'', " as text "Num. of predictors: `numrpvars`e'', " as text "Options: `opts`e''" _n
}
}
local numrp`l'vars :list sizeof rp`l'
local numrpvars = `numrpvars' + `numrp`l'vars'
if `numrp`l'vars' == 0 {
if "`diagonal'"~="" | "`elements'"~="" | ("`weightvar'"~="" & `l' > 1) | "`mmids'"~="" | "`mmweights'"~="" | "`carids'"~="" | "`carweights'"~="" || "`parexpansion'"~="" {
display as error "Cannot specify random part options at level `l' unless random part predictor variables are specified"
exit 198
}
if "`reset'"~="" | "`residuals'"~="" {
if ~(`l' == 1 & inlist("`distribution`e''","binomial","poisson","nbinomial")) {
display as error "Cannot specify reset or residual options at level `l', as no random part predictor variables are specified"
exit 198
}
}
if `l' == 1 & "`mtype'" == "univariate" & "`discrete'" == "" {
display as error "WARNING: In univariate continuous models, one or more variables should generally be included in the random part at level 1"
}
}
* Residuals option
local numresi`l' 0
if ("`residuals'"~="") {
local 0 `residuals'
syntax name(name = residname) [, VARiances STANDardised LEVerage INFluence DELetion SAMPling REFlated noRECODE RPAdjust noAdjust SAVEChains(string)]
if "`residname'" ~= "" { // If a variable stub has been given for the residual options, then...
if "`variances'" ~= "" & ("`leverage'" ~= "" | "`influence'" ~= "" | "`deletion'" ~= ""){
display as error "variances option cannot be used with leverage, deletion or influence"
exit 198
}
if "`:list residname & residualsall'"~="" {
display as error "The new residual variables stub `residname' is defined at more than one level. Specify unique stubs at each level." _n
exit 198
}
local resivars`l' `rp`l''
if `l' == 1 & "`discrete'" ~= "" { // Note, this does not currently handle multinomial
if "`mtype'"=="multinomial" {
if "`link'"=="mlogit" local resivars`l' bcons.1 bcons.2
else local resivars`l' bcons.1
}
else {
forvalues a = 1/`nummlwineqns' {
local resivars`l' `resivars`l'' bcons.`a'
if "`=word("`distribution'", `a')'" == "nbinomial" local resivars`l' `resivars`l'' bcons2.`a'
}
}
}
local numresi`l' :list sizeof resivars`l'
forvalues q = 0/`=`:list sizeof resivars`l''-1' {
capture confirm variable `residname'`q', exact
if !_rc {
display as error "The new residual variables stub `residname' is attempting to generate a new variable `residname'`q', but this variable already exists." _n
exit 198
}
capture confirm variable `residname'`q'se, exact
if !_rc {
display as error "The new residual variables stub `residname' is attempting to generate a new variable `residname'`q'se, but this variable already exists." _n
exit 198
}
}
if "`rpadjust'" ~= "" & "`adjust'" == "noadjust" {
display as error "Cannot specify both RPAdjust and NOAjust at level `l'"
exit 198
}
local residuals`l' `residname'
local residstd`l' `standardised'
local residlever`l' `leverage'
local residinfl`l' `influence'
local residdel`l' `deletion'
local residsamp`l' `sampling'
local residref`l' `reflated'
local residvar`l' `variances'
local residrpadj`l' `rpadjust'
local residnoadj`l' `adjust'
local residrecode`l' `recode'
local residualsall :list residualsall | residname
}
local saveresiduals`l' `savechains'
}
if "`flinits'" ~= "" | "`flconstraints'" ~= "" | "" ~= "`fvinits'" | "`fvconstraints'" ~= "" | "`fscores'" ~= "" local factorson = 1
if `factorson' == 1 {
if "`mtype'" ~= "multivariate" | "`mcmc'" == "" | "`diagonal'" == "" {
display as error "Factors can only be specified for multivariate MCMC models with a diagonal covariance matrix"
exit 198
}
if "`flinits'" == "" | "`flconstraints'" == "" | "" == "`fvinits'" | "`fvconstraints'" == "" {
display as error "flinits, flconstraints, fvinits and fvconstraints must all be specified for factor models"
exit 198
}
if rowsof(`flinits') ~= `:list sizeof rp`l'' | rowsof(`flconstraints') ~= `:list sizeof rp`l'' | colsof(`flinits') ~= colsof(`flconstraints') | colsof(`flinits') ~= rowsof(`fvinits') | colsof(`flinits') ~= colsof(`fvinits') | colsof(`flinits') ~= rowsof(`fvconstraints') | colsof(`flinits') ~= colsof(`fvconstraints') {
display as error "Factor matrices are wrong dimensions"
exit 198
}
tempname fvtempmat
capture mat `fvtempmat' = cholesky(`fvinits')
if c(rc) == 506 {
display as error "Factor variance matrix is not positive definite"
exit 198
}
local flinit`l' `flinits'
local flconstraint`l' `flconstraints'
local fvinit`l' `fvinits'
local fvconstraint`l' `fvconstraints'
local fscores`l' `fscores'
}
if ("`weightvar'"~="") {
if "`mtype'" ~= "univariate" | "`mcmc'" ~= "" {
display as error "Weights are only valid for univariate models estimated using (R)IGLS"
exit 198
}
local weight`l' `weightvar'
}
if ("`mmids'" ~= "" & "`carids'" ~= "") | ("`mmweights'" ~= "" & "`carweights'" ~= "") | ("`mmids'" ~= "" & "`carweights'" ~= "") | ("`mmweights'" ~= "" & "`carwids'" ~= ""){
display as error "MM and CAR cannot both be specified at the same level"
exit 198
}
local car
if "`carids'" ~= "" | "`carweights'" ~= "" {
local mmids `carids'
local mmweights `carweights'
local car car
}
if ("`mmids'" ~= "") {
local mmids`l' `mmids'
}
if ("`mmweights'" ~= "") {
if "`mcmc'" == "" | "`mmids`l''" == ""{
display as error "Multiple membership weights are only allowed in MCMC models where multiple membership has been defined"
exit 198
}
else {
local mmweights`l' `mmweights'
}
}
if "`car'" ~= "" {
if "`mmids'" == "" | "`mmweights'" == "" {
display as error "CAR is only valid for multiple membership models"
exit 198
}
else local car`l' `car'
}
if "`parexpansion'" ~= "" {
if "`mcmc'" == "" {
display as error "Parameter expansion is only valid for MCMC estimation"
exit 198
}
local parexpansion`l' `parexpansion'
local parexpansion
}
* Diagonal option
if ("`diagonal'"~="" & "`elements'"~="") { // If an on/off vector has been given for the elements options, then...
display as error "Specify either diagonal or elements(), but not both options." _n
exit 198
}
if ("`diagonal'"~="") local diagonal`l' `diagonal' // If an on/off vector has been given for the elements options, then...
* Elements option
if ("`elements'"~="") local elements`l' `elements' // If an on/off vector has been given for the elements options, then...
* Reset option
if ("`reset'"~="") { // If reset has been specified, then ...
local 0 , `reset'
syntax , [ALL VARiances NONE]
if "`all'" == "" & "`variances'" == "" & "`none'" == "" {
display as error "The reset() option must be one of: all, variances, none." _n
exit 198
}
if "`all'" ~= "" local reset`l' all
if "`variances'" ~= "" local reset`l' variances
if "`none'" ~= "" local reset`l' none
}
* Design option
if ("`design'"~="") {
local 0 `design'
syntax using [, FILE(string) Keep(string)]
preserve
use "`using'", clear
if "`keep'"~="" {
keep `keep'
order `keep' // note that we sort the variables by the order that they are specified in the keep() option.
}
unab vars: *
if "`file'"=="" tempfile file
saveold "`file'", nolabel replace
restore
local filedesign`l' `file'
local design`l' `vars'
local numdesign`l' :list sizeof vars
if ("`verbose'"~="") display as text "Number of DESIGN vectors at level `l': `numdesign`l'' (`design`l'')"
}
else local numdesign`l' = 0
if `l' > 1 & `numrp`l'vars' == 0 local emptylevels `emptylevels' `l'
}
else {
local numrp`l'vars 0
local emptylevels `emptylevels' `l'
}
}
if ("`linearization'" == "MQL2" | "`linearization'" == "PQL1" | "`linearization'" == "PQL2" ) & "`emptylevels'" ~= "" {
display as error "Only MQL1 can be specified if there are empty higher levels"
exit 198
}
* Generate an indicator for whether any design matrices have been specified at any levels
forvalues l = 1/`maxlevels' {
if ("`design`l''"~="") local anydesign yes
}
* Total number of variables (not total number of parameters) (not total number of columns used in MLwiN either as we have denom columns for example)
local numvars = `numfpvars' + `numrpvars'
* Level order
local levelorder
forvalues l = 1/`maxlevels' {
if "`lev`l'id'"~="" local levelorder `lev`l'id' `levelorder'
}
* Number of levels
local numlevels :list sizeof levelorder
******************************************************************************
* (1C) PARSE GENERAL OPTIONS
******************************************************************************
local hasweights = 0
forvalues l = 1/`maxlevels' {
if "`weight`l''" ~= "" local hasweights = 1
}
if "`weights'" ~= "" {
if `hasweights' == 0 {
display as error "weights options specified but no weights defined"
exit 198
}
local 0, `weights'
syntax , [noStandardisation noFPSandwich noRPSandwich]
if "`standardisation'" == "" local weighttype standardised
else local weighttype raw
}
// Default weight type is standardised
if "`weighttype'" == "" & `hasweights' == 1 local weighttype standardised
// Is this check correct?
if "`weighttype'" == "standardised" {
forvalues l = 2/`maxlevels' {
if "`weight`l''" ~= "" & "`weight`=`l'-1''" == ""{
display as error "Standardised weights have been requested for level `l', this requires non-equal standardised weights for all levels below `l'. Please define weights for these levels"
exit 198
}
}
}
* Assert that tolerance is integer
local iglstolerance `tolerance'
* If special option plugin is specified then do macro quietly (i.e. without pauses)
if ("`plugin'"~="") local pause no
* If savemacro option is specified then check the file extension
if ("`savemacro'"~="") {
local 0 `savemacro'
syntax anything(name=savemacro), [REPLACE]
if "`replace'"=="" confirm new file "`savemacro'"
}
* If savefullmacro option is specified then check the file extension
if ("`savefullmacro'"~="") {
local 0 `savefullmacro'
syntax anything(name=savefullmacro), [REPLACE]
if "`replace'"=="" confirm new file "`savefullmacro'"
}
if ("`saveequation'"~="") {
if ("`plugin'" ~= "") {
display as error "Cannot capture equation image in plugin mode"
exit 198
}
local 0 `saveequation'
syntax anything(name=saveequation), [REPLACE]
if "`replace'"=="" confirm new file "`saveequation'"
}
* If saveworksheet option is specified then check the file extension
if ("`saveworksheet'"~="") {
local 0 `saveworksheet'
syntax anything(name=saveworksheet), [REPLACE]
if length("`saveworksheet'") < 5 local saveworksheet `saveworksheet'.wsz
else if lower(substr("`saveworksheet'", -4, .)) ~= ".wsz" local saveworksheet `saveworksheet'.wsz
if "`replace'"=="" confirm new file "`saveworksheet'"
}
if "`initsmodel'" ~= "" {
quietly estimates dir
local savedmodels `r(names)'
if `:list initsmodel in savedmodels' {
tempname currentmodel
// NOTE: capture in case there is no current model (is there a better check?)
quietly capture estimates store `currentmodel', nocopy
quietly estimates restore `initsmodel'
tempname matPREVIOUS_b
matrix `matPREVIOUS_b' = e(b)
local inits `matPREVIOUS_b'
tempname matPREVIOUS_V
matrix `matPREVIOUS_V' = e(V)
local initsV `matPREVIOUS_V'
quietly capture estimates restore `currentmodel'
quietly capture estimates drop `currentmodel'
}
else {
display as error "saved model `initsmodel' not found"
exit 198
}
}
if "`initsv'"~="" & "`initsb'" == "" {
display as error "initsb must be specified if initsv is used"
exit 198
}
* Parse initial values
if ("`initsprevious'"~="" & "`initsb'"~="") {
display as error "Specify either initsprevious or initsb() and/or initsv(), but not both options." _n
exit 198
}
if ("`initsprevious'"~="") {
tempname matPREVIOUS_b
matrix `matPREVIOUS_b' = e(b)
local inits `matPREVIOUS_b'
tempname matPREVIOUS_V
matrix `matPREVIOUS_V' = e(V)
local initsV `matPREVIOUS_V'
}
if ("`initsb'"~="") {
local inits `initsb'
if rowsof(`inits')~=1 {
display as error "The initsb() option expects a row vector as the first argument." _n
exit 198
}
}
if ("`initsv'"~="") local initsV `initsv'
if ("`seed'"~="") local standardseed = `seed'
if "`mlwinsettings'" ~= "" {
local 0, `mlwinsettings'
syntax , [Size(numlist >0 integer min=1 max=1) Levels(numlist >=3 integer min=1 max=1) Columns(numlist >=1500 integer min=1 max=1) Variables(numlist >=1 integer min=1 max=1) TEMPMAT OPTIMAT RNGVersion(integer 10)]
if "`verbose'" ~= "" display "Maximum number of levels: `levels', Number of columns: `columns'"
if "`levels'" ~= "" local toplevel `levels'
if "`columns'" ~= "" local numwscolumns `columns'
}
******************************************************************************
* (1D) PARSE DISCRETE SUBOPTIONS
******************************************************************************
* Response type
if ("`discrete'"=="") {
local distribution
forvalues e = 1/`nummlwineqns' {
local distribution `distribution' normal
}
}
else {
if ("`nummlwineqns'" == "1" & "`distribution'" == "normal") {
display as error "The discrete option is unnecessary for univariate normal models"
exit 198
}
}
* Discrete moved higher from here
if ("`verbose'"~="") {
if ("`mtype'" == "multinomial") {
// NOTE: This assumes that the values coming out of levelsof is in the same order as the value labels, is this always true?
capture levelsof `response1', local(responsecats)
display as text "Response names:"
local count = 1
foreach cat of local responsecats {
if "`count'" ~= "`basecategory'" display as text "`:label (`response1') `cat''"
local ++count
}
}
if ("`mtype'" == "multivariate") {
display as text "Response names:"
forvalues a = 1/`nummlwineqns' {
display as text "`response`a''"
}
}
}
if ("`discrete'"~="") {
local tmpbin = "binomial"
* Check that the denom function is correctly specified.
if "`denominator'"=="" {
if `:list tmpbin in distribution' {
di as error "You must specify the denominator() function for the binomial distribution."
exit 198
}
if inlist("`distribution'","multinomial") {
di as error "You must specify the denominator() function for the multinomial distribution."
exit 198
}
}
else {
if ~`:list tmpbin in distribution' & "`distribution'" ~= "multinomial" {
display as error "Remove denominator(). denominator() is only valid for the binomial and multinomial distributions." _n
exit 198
}
}
* Check that the offset function is correctly specified.
if "`offset'" ~= "" {
local tmppois "poisson nbinomial"
if "`:list tmppois & distribution'" == "" {
display as error "Remove offset(). offset() is only valid for the Poisson or negative binomial distribution." _n
exit 198
}
}
}
if ("`mtype'" ~= "multinomial") {
if (`:list sizeof distribution' != `nummlwineqns') {
display as error "Incorrect number of distributions specified." _n
exit 198
}
local tmpbin = "binomial"
if ((`:list tmpbin in distribution') & (`:list sizeof denominator' != `nummlwineqns')) {
display as error "Incorrect number of denominators specified." _n
exit 198
}
local e = 0
foreach var of local distribution {
local ++e
local distribution`e' `var'
}
local e = 0
foreach var of local denominator {
local ++e
local denominator`e' `var'
}
* Assert that the response is well behaved
forvalues e = 1/`nummlwineqns' {
if ("`distribution`e''"=="binomial") {
capture assert (`response`e''>=0 & `response`e''<=1) | `response`e''>=.
if _rc~=0 {
di as error "Binary response variables must take the value 0 or 1. If you have a proportion response variable, then it must lie between 0 and 1."
exit 198
}
}
if ("`distribution`e''"=="poisson") {
capture assert `response`e''>=0
if _rc~=0 {
di as error "The count response variable must be 0 or a positive integer."
exit 198
}
}
if ("`distribution`e''"=="nbinomial") {
capture assert `response`e''>=0
if _rc~=0 {
di as error "The count response variable must be 0 or a positive integer."
exit 198
}
}
if ~inlist("`distribution`e''","normal","binomial","poisson","nbinomial") {
display as error "The `distribution`e'' response distribution is not recognised by MLwiN." _n
exit 198
}
}
* Check that no level 1 variables have been specified
if inlist("`distribution`e''","binomial","poisson","nbinomial") & "`rp1'"~="" {
display as error "You cannot make covariates random at level 1 when using the `distribution`e'' response distribution." _n
exit 198
}
* Check that the denominator has been specified
if inlist("`distribution`e''","binomial") & "`denominator`e''"=="" {
display as error "You must specify a denominator variable when using the `distribution`e'' response distribution. See the denom() option." _n
exit 198
}
}
else {
qui inspect `response1'
if r(N_unique)<=1 {
display as error "In multinomial response models, the response must have two or more categories" _n
exit 198
}
* Check that no level 1 variables have been specified
if "`rp1'" ~= "" {
display as error "You cannot make covariates random at level 1 when using the multinomial response distribution." _n
exit 198
}
* Check that the denominator has been specified
if "`denominator'"=="" {
display as error "You must specify a denominator variable when using the multinomial response distribution. See the denom() option." _n
exit 198
}
* Check that the base category has been specified
if "`basecategory'"=="" {
display as error "You must specify a (numeric) base category when using the multinomial response distribution. See the basecategory() option." _n
exit 198
}
* Check that a valid base category has been specified for models with unordered response variables
qui levelsof `response1', local(responsevalues)
local responsevalues = subinstr("`responsevalues'"," ",",",.)
if ~inlist(`basecategory',`responsevalues') {
display as error "The base category `basecategory' does not exist for the response variable `response1'. See the basecategory() option." _n
exit 198
}
}
******************************************************************************
* (1E) PARSE MCMC SUBOPTIONS
******************************************************************************
* Estimation model
if ("`mcmc'"=="" & "`rigls'"=="") local emode IGLS
if ("`mcmc'"=="" & "`rigls'"~="") local emode RIGLS
if ("`mcmc'"~="") local emode MCMC
if ("`mcmc'"~="") {
if "`initsprevious'"=="" & "`initsb'"=="" & "`initsmodel'" == "" {
display as error "When fitting models using MCMC, initial values must be specified; see initsprevious, initsmodel and initsb()/initsv() options."
exit 198
}
* The syntax for the burn-in is as follows:
*
* MCMC 0
* burn in for <value 1> iterations,
* use adaptation method <value 2>,
* scale factor <value 3>,
* %acceptance <value 4>,
* tolerance <value 5>
* {residual starting values in C1,
* s.e. residual starting values) in C2}
* {priors in C3}
* fixed effect method <value 6>
* residual method <value 7>
* level 1 variance method <value 8>
* other levels variance method <value 9>
* default prior <value 10>
* model type <value 11>.
local 0 , `mcmc'
#delimit ;
syntax , [
ON
LOGformulation
CC
MSUBScripts
CORResiduals(namelist min=1 max=1)
ME(string)
PRIORMatrix(namelist min=1 max=1)
RPPriors(string)
SAVEWinbugs(string)
Burnin(integer 500)
Chain(integer 5000)
Thinning(integer 1)
Refresh(integer 50)
SCALE(real 5.8)
noADAPTation
ACCEPTance(real 0.5)
TOLerance(real 0.1)
CYCLES(integer 1)
FEMethod(namelist min=1 max=1)
REMethod(namelist min=1 max=1)
LEVELONEVARMETHOD(namelist min=1 max=1)
HIGHERLEVELSVARMETHOD(namelist min=1 max=1)
SMCMC
SMVN
ORTHogonal
HCentring(numlist min=1 max=1 >0 integer)
CARCentre
SEED(integer 1)
noDIAGnostics
SAVEChains(string)
SIMRESiduals /* Undocumented option. The default is to calculate starting values for higher level residuals using the residuals window. This is a potential alternative which simulates residual starting values from the MVN residual differences (i.e. a drawnorm). */
IMPUTEIterations(numlist >0 integer ascending)
IMPUTESummaries
];
#delimit cr
if `thinning' > 0.5*`chain' {
display as error "Thinning must be less than half the chain length."
exit 198
}
if ("`fpsandwich'"=="fpsandwich" | ("`weighttype'"~="" & "`fpsandwich'"=="")) | ("`rpsandwich'"=="rpsandwich" | ("`weighttype'"~="" & "`rpsandwich'"=="")){
display as error "Sandwich estimates are not available for MCMC"
exit 198
}
if "`hcentring'" ~= "" {
if `:list hcentring in emptylevels' {
display as error "Hierarchical centring cannot be specified for empty levels"
exit 198
}
if `hcentring' < 2 | `hcentring' > `maxlevels' {
display as error "Hierarchical centring must be at one of levels 2 to `maxlevels'"
exit 198
}
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") local ++hcentring
}
if "`imputeiterations'" ~= "" {
foreach val of local imputeiterations {
if mod(`val', `refresh') ~= 0 {
display as error "`val' is not a valid iteration to impute at"
exit 198
}
}
}
if ("`me'"~="") {
if "`mtype'" == "multivariate" {
display as error "Measurement error not implemented for multivariate in MLwiN." _n
exit 198
}
local 0 `me'
syntax varlist(min=1), variances(numlist min=1 >=0)
if `:list sizeof varlist' != `:list sizeof variances' {
display as error "incorrect number of measurement error variances"
exit 198
}
local mevars `varlist'
}
if ("`discrete'"=="") local temp gibbs
if ("`discrete'"~="") local temp univariatemh // univariate MH // NOTE SOMETIMES THIS WILL NEET TO BE MULTIVARITE MH
if ("`femethod'"=="") local femethod `temp'
if ("`remethod'"=="") local remethod `temp'
if ("`levelonevarmethod'"=="") local levelonevarmethod `remethod' // Looks like the GUI does this as well
if ("`higherlevelsvarmethod'"=="") local higherlevelsvarmethod gibbs
macro drop temp
if ("`savechains'"~="") {
local 0 `savechains'
syntax anything(name=savechains), [REPLACE]
// Add .dta extension if this has been omitted
if lower(substr("`savechains'", -4, .)) ~= ".dta" {
local savechains `savechains'.dta
}
if "`replace'" == "" confirm new file "`savechains'"
}
else {
tempfile savechains
}
forvalues l = 1/`maxlevels' {
if "`saveresiduals`l''" ~= "" {
local 0 `saveresiduals`l''
syntax anything(name=saveresiduals`l'), [REPLACE]
// Add .dta extension if this has been omitted
if lower(substr("`saveresiduals`l''", -4, .)) ~= ".dta" {
local saveresiduals`l' `saveresiduals`l''.dta
}
if "`replace'" == "" confirm new file "`saveresiduals`l''"
}
}
local defaultprior 1
local gama 0.001
local gamb 0.001
if "`rppriors'" ~= "" {
local 0 , `rppriors'
syntax , [Uniform Gamma(numlist min=2 max=2 >=0)]
if "`uniform'" ~= "" & "`gamma'" ~= "" {
display as error "Cannot specify both uniform and gamma prior"
exit 198
}
if "`uniform'" == "" & "`gamma'" == "" {
display as error "You must specify uniform or gamma(a b) for priors"
exit 198
}
if "`uniform'" ~= "" local defaultprior 0
if "`gamma'" ~= "" {
local defaultprior 1
local gama :word 1 of `gamma'
local gamb :word 2 of `gamma'
}
}
if "`corresiduals'" ~= "" {
local 0 , `corresiduals'
syntax , [UNstructured EXchangeable AREQvars EQCORRSINDEPvars ARINDEPvars]
if "`unstructured'"=="" & "`exchangeable'"=="" & "`areqvars'"=="" & "`eqcorrsindepvars'"=="" & "`arindepvars'"=="" {
display as error "Invalid level 1 residual correlation structure specified."
exit 198
}
}
// really we want method(femethod(string) remethod(string) scale(real 5.8) adaptive(acceptance(integer 50) tolerance(integer 10)) cycles(integer 1))
local mcmc_burnin = `burnin' // MLwiN MCMC command option 1
if "`adaptation'" == "noadaptation" local mcmc_adaptation = 0 // MLwiN MCMC command option 2
else local mcmc_adaptation = 1
local mcmc_scale = `scale' // MLwiN MCMC command option 3
if `acceptance' < 0 | `acceptance' > 1 {
display as error "Acceptance rate must be between zero and one inclusive"
exit 198
}
local mcmc_acceptance = ceil(`acceptance'*100) // MLwiN MCMC command option 4
if `tolerance' < 0 | `tolerance' > 1 {
display as error "Tolerance rate must be between zero and one inclusive"
exit 198
}
local mcmc_tolerance = ceil(`tolerance'* 100) // MLwiN MCMC command option 5
local 0 , `femethod'
syntax , [GIBBS UNIvariatemh MULTIvariatemh]
if "`gibbs'" ~= "" local mcmc_femethod = 1 // 1 = gibbs
if "`univariatemh'" ~= "" local mcmc_femethod = 2 // 2 = univariate MH
if "`multivariatemh'" ~= "" local mcmc_femethod = 3 // 3 = multivariate MH
if "`gibbs'" == "" & "`univariatemh'" == "" & "`multivariatemh'" == "" {
display as error "`femethod' is not valid"
exit 198
}
local 0 , `remethod'
syntax , [GIBBS UNIvariatemh MULTIvariatemh]
if "`gibbs'" ~= "" local mcmc_remethod = 1 // 1 = gibbs
if "`univariatemh'" ~= "" local mcmc_remethod = 2 // 2 = univariate MH
if "`multivariatemh'" ~= "" local mcmc_remethod = 3 // 3 = multivariate MH
if "`gibbs'" == "" & "`univariatemh'" == "" & "`multivariatemh'" == "" {
display as error "`remethod' is not valid"
exit 198
}
* Looks like you can't set this through the GUI
local 0 , `levelonevarmethod'
syntax , [GIBBS UNIvariatemh MULTIvariatemh]
if "`gibbs'" ~= "" local mcmc_levelonevarmethod = 1 // 1 = gibbs
if "`univariatemh'" ~= "" local mcmc_levelonevarmethod = 2 // 2 = univariate MH
if "`multivariatemh'" ~= "" local mcmc_levelonevarmethod = 3 // 3 = multivariate MH
if "`gibbs'" == "" & "`univariatemh'" == "" & "`multivariatemh'" == "" {
display as error "`levelonevarmethod' is not valid"
exit 198
}
local 0 , `higherlevelsvarmethod'
syntax , [GIBBS UNIvariatemh MULTIvariatemh]
if "`gibbs'" ~= "" local mcmc_higherlevelsvarmethod = 1 // 1 = gibbs
if "`univariatemh'" ~= "" local mcmc_higherlevelsvarmethod = 2 // 2 = univariate MH
if "`multivariatemh'" ~= "" local mcmc_higherlevelsvarmethod = 3 // 3 = multivariate MH
if "`gibbs'" == "" & "`univariatemh'" == "" & "`multivariatemh'" == "" {
display as error "`higherlevelsvarmethod' is not valid"
exit 198
}
local mcmc_defaultprior = `defaultprior' // MLwiN MCMC command option 10
if "`mtype'" == "univariate" & "`distribution'"=="normal" local mcmc_modeltype = 1 // Normal models
if "`mtype'" == "univariate" & "`distribution'"=="binomial" local mcmc_modeltype = 2 // Binomial models
if "`mtype'" == "univariate" & "`distribution'"=="poisson" local mcmc_modeltype = 3 // Poisson models
if "`mtype'" == "univariate" & "`distribution'"=="nbinomial" local mcmc_modeltype = 8 // Negative-binomial models
if "`mtype'" == "multivariate" & "`discrete'" == "" local mcmc_modeltype = 4 // Multivariate normal models
if "`mtype'" == "multivariate" & "`discrete'" ~= "" local mcmc_modeltype = 5 // Multivariate mixed models (must be normal/binomial with probit link)
if "`mtype'" == "multinomial" & "`link'"=="mlogit" local mcmc_modeltype = 6 // Unordered multinomial models
if "`mtype'" == "multinomial" & ("`link'"=="ologit"|"`link'"=="oprobit"|"`link'"=="ocloglog") local mcmc_modeltype = 7 // Ordered multinomial models
capture assert inlist(`mcmc_modeltype',1,2,3,4,5,6,7,8)
if _rc {
display as error "The selected distribution is not available for MCMC"
exit 198
}
if "`verbose'"~="" {
local i = 1
foreach option in mcmc_burnin mcmc_adaptation mcmc_scale mcmc_acceptance mcmc_tolerance mcmc_femethod mcmc_remethod mcmc_levelonevarmethod mcmc_higherlevelsvarmethod mcmc_defaultprior mcmc_modeltype {
di "MLwiN MCMC command option `i': ``option''" _col(50) "(`option')"
local ++i
}
}
local resptype
local mixed 0
if "`distribution'" == "binomial" local resptype 0
if "`distribution'" == "poisson" local resptype 1
if "`distribution'" == "nbinomial" local resptype 2
if "`distribution'" == "normal" local resptype 3
if "`mtype'" == "multinomial" & "`link'"=="mlogit" local resptype 4
if "`mtype'" == "multinomial" & ("`link'"=="ologit"|"`link'"=="oprobit"|"`link'"=="ocloglog") local resptype 5
if "`mtype'" == "multivariate" {
local resptype 3
local firstdist :word 1 of `distribution'
foreach dist in `distribution' { // If one of the response is nonlinear set response type to this
if "`dist'" ~= "`firstdist'" local mixed 1
if "`dist'" == "binomial" {
local resptype 0
if "`link'" ~= "probit" {
display as error "Only the probit link is allowed for binomial responses in multivariate models for MCMC estimation"
exit 198
}
}
if "`dist'" == "poisson" {
local resptype 1
display as error "Poisson is currently unavailable in multivariate models for MCMC estimation"
exit 198
}
if "`dist'" == "nbinomial" {
local resptype 2
display as error "Poisson is currently unavailable in multivariate models for MCMC estimation"
exit 198
}
}
}
local allnormal 1
local resptypes `distribution'
forvalues n = 1/`:list sizeof resptypes' {
if "`=word("`resptypes'", `n')'" != "normal" local allnormal 0
}
local mcmc_chain = `chain'
local mcmc_thinning = `thinning'
local mcmc_refresh = `refresh'
local mcmc_seed = `seed'
forvalues l = 2/`maxlevels' {
if "`mmids`l''"~="" |"`mmweights`l''"~="" {
if `:list sizeof mmweights`l'' ~= `:list sizeof mmids`l'' {
display as error "Number of variables in mmids`l'() must equal number of variables in mmweights`l'()."
exit 198
}
mata: checkmm("`mmids`l''", "`mmweights`l''")
local numweights`l' = `:list sizeof mmwewights`l''
local cc on
local mm on
}
}
/*
Valid MCMC options:
smcmc smvn orth hcen parexpansion
1 0 0 0 0
0 1 0 0 0
0 0 1 1 1
*/
local anyparexpansion = 0
forvalues l = 2/`maxlevels' {
if "`parexpansion`l''" ~= "" local anyparexpansion = 1
}
if "`smcmc'" ~= "" & ("`smvn'" ~= "" | "`orthogonal'" ~= "" | "`hcentring'" ~= "" | `anyparexpansion' > 0) {
display as error "SMCMC is exclusive to SMVN, ORTHogonal, HCEN and PAREXpansion"
exit 198
}
if "`smvn'" ~= "" & ("`smcmc'" ~= "" | "`orthogonal'" ~= "" | "`hcentring'" ~= "" | `anyparexpansion' > 0) {
display as error "SMVN is exclusive to SMCMC, ORTHogonal, HCEN and PAREXpansion"
exit 198
}
if "`smcmc'" ~= "" & ("`mtype'" ~= "univariate" | "`discrete'" ~= "" | "`numrp1vars'" ~= "1" | "`maxlevels'" ~= "2") {
display as error "SMCMC is only valid for two-level normal response models with no complex level 1."
exit 198
}
}
local carlev = 0
forvalues l = 1/`maxlevels' {
if ("`car`l''"~="" & "`carcentre'"=="") {
if `carlev' ~= 0 {
display as error "CAR can only be specified at one level"
exit 198
}
else local carlev = `l'
}
}
******************************************************************************
* (2) GENERATE THE EMPTY VECTOR OF PARAMETER ESTIMATES
******************************************************************************
* Fixed part parameter vector(s)
forvalues r = 1/`numfpbrackets' {
tempname FP_b`r'
mat `FP_b`r'' = J(1,`numfpvars`r'',.)
mat coleq `FP_b`r'' = FP`r'
local fname
foreach name of local fpname`r' {
if _caller() >= 11 {
local fname `fname' `=strtoname("`name'")'
}
else {
mata: st_local("name", validname("`name'"))
local fname `fname' `=abbrev("`name'", 32)'
}
}
mat colnames `FP_b`r'' = `fname'
}
tempname matFP_b
forvalues r = 1/`numfpbrackets' {
mat `matFP_b' = (nullmat(`matFP_b'),`FP_b`r'')
}
* Random part parameter vector (some terms may still need to be deleted)
forvalues l = `maxlevels'(-1)1 {
tempname RP`l'_b // define matrix names as locals
if `:list sizeof rp`l'name'>=1 {
runmlwin_rplpars , rpname(`rp`l'name')
local rp`l'pars `r(rplong)'
runmlwin_rplpars , rpname(`origrp`l'name')
local origrp`l'pars `r(rplong)'
local numrp`l'pars: word count `rp`l'pars'
if `numrp`l'pars' > 0 {
mat `RP`l'_b' = J(1,`numrp`l'pars',.)
mat coleq `RP`l'_b' = RP`l'
mat colnames `RP`l'_b' = `rp`l'pars'
}
* Remove parameters from the variance-covariance matrices
if "`elements`l''"~="" {
* Assert that length of `elements`l'' is same as length of `rp`l'pars'
if colsof(`elements`l'')~=`numrp`l'pars' {
display as error "The matrix `elements`l'' should contain `numrp`l'pars' elements (i.e. one for each parameter in the level `l' var-cov matrix)" _n
exit 198
}
* Assert that the elements of `elements`l'' are 0 or 1
local names : colfullnames `elements`l''
local c = 1
foreach p of local names {
if ~inlist(`elements`l''[1,`c'],0,1) {
display as error "All elements in matrix `elements`l'' should be 0 (= remove parameter from model) or 1 (= retain parameter in model)" _n
exit 198
}
local ++c
}
* Remove elements from RP_b
tempname matTEMP
local names : colnames `RP`l'_b'
local keepnames
forvalues i = 1/`=colsof(`elements`l'')' {
if (`elements`l''[1,`i']==1) local keepnames `keepnames' `:word `i' of `origrp`l'pars''
}
foreach p of local names {
if `:list p in keepnames' matrix `matTEMP' = (nullmat(`matTEMP'), `RP`l'_b'[1,"RP`l':`p'"])
}
mat `RP`l'_b' = `matTEMP'
}
if "`diagonal`l''"~="" {
* Work out what the diagonal matrix will look like, then convert the lower diagonal to a vector
tempname matTEMP1 matTEMP2
mat `matTEMP1' = I(`numrp`l'vars')
local C = `numrp`l'vars'
local R = `numrp`l'vars'
forvalues r=1/`R' {
forvalues c=1/`C' {
if `c'<=`r' mat `matTEMP2' = (nullmat(`matTEMP2'), `matTEMP1'[`r',`c'])
}
}
* Remove elements from RP_b
// This code is same as that in the elements`l' statement above. Code could be made more efficient
tempname matTEMP3
local c = 1
local names : colfullnames `RP`l'_b'
foreach p of local names {
if (`matTEMP2'[1,`c']==1) matrix `matTEMP3' = (nullmat(`matTEMP3'), `RP`l'_b'[1,"`p'"])
local ++c
}
mat `RP`l'_b' = `matTEMP3'
}
}
}
tempname matRP_b
forvalues l = `maxlevels'(-1)1 {
if `:list sizeof rp`l'name'>=1 {
if `numrp`l'pars' > 0 {
mat `matRP_b' = (nullmat(`matRP_b'), `RP`l'_b')
}
}
}
if "`distribution'" == "binomial" | "`distribution'" == "poisson" | "`distribution'" == "nbinomial" | "`distribution'" == "multinomial"{
tempname mat_OD
if "`distribution'" == "nbinomial" local odvars bcons_1 bcons2_1
else local odvars bcons_1
if "`distribution'" == "multinomial" & "`link'"=="mlogit" & "`mcmc'" == "" local odvars `odvars' bcons_2
mat `mat_OD' = J(1, `:list sizeof odvars', .)
mat colnames `mat_OD' = `odvars'
mat coleq `mat_OD' = OD
mat `matRP_b' = (nullmat(`matRP_b'), `mat_OD')
}
* Random part design
forvalues l = `maxlevels'(-1)1 {
if "`design`l''"~="" {
mat `RP`l'_b'design = J(1,`numdesign`l'',.)
mat coleq `RP`l'_b'design = RP`l'D
mat colnames `RP`l'_b'design = `design`l''
mat `matRP_b' = (nullmat(`matRP_b'), `RP`l'_b'design)
}
}
* Factors
tempname mat_RP_b_fact
forvalues l = `maxlevels'(-1)1 {
if "`flinit`l''"~="" {
tempname matTEMP
forvalues fact = 1/`=colsof(`flinit`l'')' {
local factnames
forvalues resp = 1/`=rowsof(`flinit`l'')' {
local factnames `factnames' f`fact'_`resp'
}
mat `matTEMP' = J(1, rowsof(`flinit`l''), .)
mat coleq `matTEMP' = RP`l'FL
mat colnames `matTEMP' = `factnames'
mat `mat_RP_b_fact' = (nullmat(`mat_RP_b_fact'), `matTEMP')
}
}
}
tempname mat_RP_b_factvar
forvalues l = 1/`maxlevels' {
if "`flinit`l''" ~= "" {
tempname matTEMP
local factvarnames
forvalues fact1 = 1/`=colsof(`flinit`l'')' {
forvalues fact2 = 1/`fact1' {
if `fact1' == `fact2' local factvarnames `factvarnames' var(f`fact1')
else local factvarnames `factvarnames' cov(f`fact1'\f`fact2')
}
}
mat `matTEMP' = J(1, `:list sizeof factvarnames', .)
mat coleq `matTEMP' = RP`l'FV
mat colnames `matTEMP' = `factvarnames'
mat `mat_RP_b_factvar' = (nullmat(`mat_RP_b_factvar'), `matTEMP')
}
}
local C = colsof(`matFP_b')
local namesV : colfullnames `matFP_b'
tempname matFP_V
mat `matFP_V' = J(`C',`C',0)
matrix rownames `matFP_V' = `namesV'
matrix colnames `matFP_V' = `namesV'
* Put MLwiN values back into RP_V
local C = colsof(`matRP_b')
local namesV : colfullnames `matRP_b'
tempname matRP_V
mat `matRP_V' = J(`C',`C',0)
matrix rownames `matRP_V' = `namesV'
matrix colnames `matRP_V' = `namesV'
* Full parameter matrix (this is what Stata expects MLwiN to return)
* Create matb with the correct dimensions
tempname matb
mat `matb' = (`matFP_b',`matRP_b', nullmat(`mat_RP_b_fact'), nullmat(`mat_RP_b_factvar'))
mat rowname `matb' = `response1'
local names : colfullnames `matb'
tempname matV
mat `matV' = J(colsof(`matb'), colsof(`matb'), 0)
local namesV : colfullnames `matb'
mat rownames `matV' = `namesV'
mat colnames `matV' = `namesV'
if ("`constraints'"~="") {
* Parse contraint(s)
local C = colsof(`matb')
forvalues c=1/`C' {
mat `matb'[1,`c'] = 0
}
ereturn post `matb' // Need to create e(b) for the current model in order to set constraints
tempname matb
mat `matb' = e(b)
if ("`verbose'"~="") mat list e(b)
forvalues c=1/`C' {
mat `matb'[1,`c'] = .
}
tempname consT
tempname consa
tempname consC
makecns `constraints', displaycns
matcproc `consT' `consa' `consC'
if "`verbose'" ~= "" display "Constraints:"
if "`verbose'" ~= "" mat list `consC'
* Extract the fixed part constraint submatrix (note that some rows will refer to random part constraints, we delete these below)
tempname FPCtemp
mat `FPCtemp' = (`consC'[.,1..`numfpvars'],`consC'[.,`=colsof(`consC')'..`=colsof(`consC')'])
* Extract the random part constraint submatrix (note that some rows will refer to fixed part constraints, we delete these below)
local temp = `numfpvars' + 1
tempname RPCtemp
mat `RPCtemp' = (`consC'[.,`temp'...])
if "`verbose'" ~= "" display "FPCtemp"
if "`verbose'" ~= "" mat list `FPCtemp'
if "`verbose'" ~= "" display "RPCtemp"
if "`verbose'" ~= "" mat list `RPCtemp'
local numconstraints = rowsof(`consC')
* Work out which constraints and how many constraints are FP constraints
tempname isFPconstraint
mat `isFPconstraint' = J(`numconstraints',1,.)
forvalues r=1/`=rowsof(`FPCtemp')'{
local flag = 0
forvalues c=1/`=colsof(`FPCtemp')-1' {
if `FPCtemp'[`r',`c'] ~= 0 local flag = 1
}
mat `isFPconstraint'[`r',1]==`flag'
}
tempname mat_sum
mat `mat_sum' = J(1,`numconstraints',1)*`isFPconstraint'
local numFPconstraints = `mat_sum'[1,1]
* Work out which constraints and how many constraints are RP constraints
tempname isRPconstraint
mat `isRPconstraint' = J(`numconstraints',1,.)
forvalues r=1/`=rowsof(`RPCtemp')' {
local flag = 0
forvalues c=1/`=colsof(`RPCtemp')-1' {
if `RPCtemp'[`r',`c'] ~= 0 local flag = 1
}
mat `isRPconstraint'[`r',1]==`flag'
}
mat `mat_sum' = J(1,`numconstraints',1)*`isRPconstraint'
local numRPconstraints = `mat_sum'[1,1]
if "`verbose'" ~= "" display "num fp vars: `numfpvars', num rp vars: `numrpvars'"
if "`verbose'" ~= "" display "isFPconstraint:"
if "`verbose'" ~= "" mat list `isFPconstraint'
if "`verbose'" ~= "" display "isRPconstraint:"
if "`verbose'" ~= "" mat list `isRPconstraint'
if "`verbose'" ~= "" display "num constraints: `numconstraints', num FP constraints: `numFPconstraints', num RP constraints: `numRPconstraints'"
* Check that each constraint relates to only the fixed part or only the random part
tempname matTEMP
mat `matTEMP' = `isFPconstraint' + `isRPconstraint'
forvalues r=1/`numconstraints' {
if (`matTEMP'[`r',1]>1) {
mat list `consC'
display as error "Row `r' of the constraint matrix is invalid as it involves both fixed part and random part parameters." _n
exit 198
}
}
* Strip out redundant rows of the fixed part constraint matrix (these relate to the random part)
tempname FPC
if `numFPconstraints'>0 {
forvalues r=1/`numconstraints' {
if (`isFPconstraint'[`r',1]==1) mat `FPC' = (nullmat(`FPC') \ `FPCtemp'[`r',1...])
}
local fpconstraints `FPC'
}
* Strip out redundant rows of the random part constraint matrix (these relate to the fixed part)
tempname RPC
if `numRPconstraints'>0 {
forvalues r=1/`numconstraints' {
if (`isRPconstraint'[`r',1]==1) mat `RPC' = (nullmat(`RPC') \ `RPCtemp'[`r',1...])
}
local rpconstraints `RPC'
}
ereturn clear // As we don't want an e(b) kicking around
}
* Store the fixed part constraint(s) as strings in locals ready to insert into MLwiN macro
if ("`fpconstraints'"~="") {
local numfpconstraints = rowsof(`fpconstraints')
forvalues r = 1/`numfpconstraints' {
local C = colsof(`fpconstraints')
local fpconstraint`r' = `fpconstraints'[`r',1] // put first element of constraint into the string
forvalues c=2/`C' {
local temp = `fpconstraints'[`r',`c']
local fpconstraint`r' `fpconstraint`r'' `temp' // put second, third,... elements of constraint into the string
}
}
if "`verbose'" ~= "" display("Fixed part constraints:")
if "`verbose'" ~= "" mat list `fpconstraints'
}
* Store the random part constraint(s) as strings in locals ready to insert into MLwiN macro
if "`rpconstraints'"~="" {
local numrpconstraints = rowsof(`rpconstraints')
forvalues r = 1/`numrpconstraints' {
local C = colsof(`rpconstraints')
local rpconstraint`r' = `rpconstraints'[`r',1] // put first element of constraint into the string
forvalues c=2/`C' {
local temp = `rpconstraints'[`r',`c']
local rpconstraint`r' `rpconstraint`r'' `temp' // put second, third,... elements of constraint into the string
}
}
if "`verbose'" ~= "" display("Random part constraints:")
if "`verbose'" ~= "" mat list `rpconstraints'
}
* Substitute initial values into the parameter vector
if "`inits'"~="" | "`initsmodel'" ~= "" {
if "`initsprevious'"~="" | "`initsmodel'" ~= "" {
tempname srcrowidx
tempname srccolidx
tempname destrowidx
tempname destcolidx
matrix `srcrowidx' = J(1,1,1)
matrix `destrowidx' = J(1,1,1)
local match 0
foreach p of local names {
local idx = colnumb(`inits',"`p'")
if "`idx'" ~= "." {
local match 1
matrix `srccolidx' = (nullmat(`srccolidx'), `idx')
matrix `destcolidx' = (nullmat(`destcolidx'), `=colnumb(`matb',"`p'")')
}
}
if `match' == 1 {
mata: copymatrix("`inits'", "`srcrowidx'", "`srccolidx'", "`matb'", "`destrowidx'", "`destcolidx'");
}
else {
if "`initsprevious'"~="" display as error "The initsprevious option is invalid: No parameters in the previous model match those specified in the current model"
if "`initsmodel'" ~= "" display as error "The initsmodel option is invalid: No parameters in the stored model match those specified in the current model"
exit 198
}
matrix drop `srcrowidx'
matrix drop `destrowidx'
matrix drop `srccolidx'
matrix drop `destcolidx'
local match 0
foreach p1 of local namesV {
local cidx = colnumb(`initsV',"`p1'")
local ridx = rownumb(`initsV',"`p1'")
if "`ridx'" ~= "." && "`cinx'" ~= "." {
local match 1
matrix `srcrowidx' = (nullmat(`srcrowidx'), `ridx')
matrix `destrowidx' = (nullmat(`destrowidx'), `=rownumb(`matV', "`p1'")')
matrix `srccolidx' = (nullmat(`srccolidx'), `cidx')
matrix `destcolidx' = (nullmat(`destcolidx'), `=colnumb(`matV', "`p1'")')
}
}
if `match' == 1 {
mata: copymatrix("`initsV'", "`srcrowidx'", "`srccolidx'", "`matV'", "`destrowidx'", "`destcolidx'");
}
else {
if "`initsprevious'"~="" display as error "The initsprevious option is invalid: No parameters in the previous model match those specified in the current model"
if "`initsmodel'" ~= "" display as error "The initsmodel option is invalid: No parameters in the stored model match those specified in the current model"
exit 198
}
matrix drop `srcrowidx'
matrix drop `destrowidx'
matrix drop `srccolidx'
matrix drop `destcolidx'
}
if "`initsb'"~="" {
if colsof(`inits') ~= `=colsof(`matFP_b')+colsof(`matRP_b')' {
display as error "The row vector `inits' of initial values in initsb() is of length `=colsof(`inits')', but -runmlwin- expects the row vector to be of length `=colsof(`matFP_b')+colsof(`matRP_b')'."
exit 198
}
matrix `matb'[1,1] = `inits'
* Note this is not currently used as we always check there are two inits matrices
if "`initsV'" == "" {
forvalues i = 1/`=colsof(`matb')' {
matrix `matV'[`i', `i'] = (`inits'[1,`i']/2)*(`inits'[1,`i']/2)
}
}
}
if "`initsv'"~="" {
if colsof(`initsV') ~= `=colsof(`matFP_b')+colsof(`matRP_b')' | rowsof(`initsV') ~= `=colsof(`matFP_b')+colsof(`matRP_b')' {
display as error "The matrix of initial variance values in initsv() is of a different length to that expected by -runmlwin-."
exit 198
}
matrix `matV'[1,1] = `initsV'
}
}
if "`mlwin'" == "nomlwin" {
ereturn clear
ereturn matrix matb = `matb'
ereturn matrix matV = `matV'
}
******************************************************************************
* (3) PREPARE DATA FOR MLWIN
******************************************************************************
* Preserve the data
preserve
* Remove data label
label data
* Generate unique identifier variable
label var `_sortindex' "Observation number in the current Stata data set"
* Keep the estimation sample
qui keep if `touse'
local _nobs = _N
* List of values of the response variable
if ("`mtype'"=="multinomial") {
qui levelsof `response1', local(responsecats)
if "`proportion'" ~= "" {
capture assert `proportion' >= 0 & `proportion' <= 1
if _rc {
display as error "Values in the multinomial proportion variable should lie between zero and one"
exit 198
}
tempvar propsum
if ("`cc'"~="") local hier `lev1id'
else {
local hier
local levelused = 0
forvalues l = `maxlevels'(-1)2 {
if `numrp`l'vars' > 0 local levelused = 1
if `levelused' == 1 local hier `hier' `lev`l'id'
}
local hier `hier' `lev1id'
}
bysort `hier': egen `propsum' = total(`proportion')
capture assert inrange(`propsum', 0.995, 1.005)
if _rc {
display as error "The sum of multinomial proportions for each unit should equal one"
exit 198
}
drop `propsum'
sort `_sortindex'
drop `_sortindex'
quietly reshape wide `proportion', i(`lev1id') j(`response1')
gen `_sortindex' = _n
gen `response1' = 1
local propvars
local propvarnames
local i = 1
foreach r of local responsecats {
if `r' ~= `basecategory' local propvars `propvars' `proportion'`r'
if `r' ~= `basecategory' local propvarnames `propvarnames' '`proportion'`r''
quietly replace `response1' = `r' in `i'
local ++i
}
}
tempname newvaluelabel
foreach r of local responsecats {
label define `newvaluelabel' `r' "`: label (`response1') `r''", add
}
label values `response1' `newvaluelabel'
if "`link'"=="ologit"|"`link'"=="oprobit"|"`link'"=="ocloglog" {
local firstvalue = real(word("`responsecats'",1))
local lastvalue = real(word("`responsecats'",-1))
if ~inlist(`basecategory', `firstvalue',`lastvalue') {
di as error "The basecategory must be either the first value (`response1' = `firstvalue') or last value (`response1' = `lastvalue') of the ordered response variable." _n
exit 198
}
}
}
* List of ID variables
local idvars
forvalues l = `maxlevels'(-1)1 {
local idvars `idvars' `lev`l'id'
// MM IDs are added consecutively in case MM IDs share variables with level IDs (what if ID is used in more than one level?)
if "`mmids`l''" ~= "" local idvars `idvars' `mmids`l''
}
* List of variables which appear in the random part of the model
local rp
forvalues l = `maxlevels'(-1)1 {
local rp `rp' `rp`l''
}
* List of variable which appear as denominators in discrete response models
local denomvars
if "`mtype'" == "multivariate" {
forvalues e = 1/`nummlwineqns' {
local denomvars `denomvars' `denominator`e''
}
}
else {
local denomvars `denominator'
}
* List of multiple membership weight variables
local mmweightvars
forvalues l = 2/`maxlevels' {
if "`mmweights`l''" ~= "" local mmweightvars `mmweightvars' `mmweights`l''
}
* List of (R)IGLS sampling weight variables
local weightvars
forvalues l = 1/`maxlevels' {
if "`weight`l''" ~= "" local weightvars `weightvars' `weight`l''
}
* Keep variables that appear in model
local allvarsformlwin ///
`idvars' ///
`response' ///
`fp' ///
`rp' ///
`denomvars' `offset' `fpconst' `rpconst' `mmweightvars' `mmweightids' `weightvars' `_sortindex' `propvars'
local data_has_missing_values = 0
foreach var of local allvarsformlwin {
if "`var'" ~= "bcons" {
quietly count if missing(`var')
if r(N) > 0 {
if "`mtype'" == "multinomial" {
display as error "For multinomial response models (i.e. ordered and unordered categorical responses), the data must be manually listwise deleted (i.e. cases with missing values for one or more variables included in the model must be dropped) prior to calling runmlwin."
exit 198
}
local data_has_missing_values = 1
}
}
}
* If a multivariate response model, create and add bcons variable to the list of variables sent to MLwiN
if "`mtype'" == "multivariate" gen byte bcons = 1
* Keep variables
if "`drop'"=="" keep `allvarsformlwin'
order `allvarsformlwin'
* Strip off the value labels of all fixed or random part covariates
local temp `fp' `rp'
foreach var of local temp {
label values `var'
}
* Calculate number of columns that will be used in the worksheet
qui describe
local numcolumns = r(k)
local N = r(N)
if "`toplevel'" == "" {
local toplevel `maxlevels'
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") local ++toplevel
local ++toplevel // Add an extra level, as NLEV command gives the wrong answer if all available levels are in use and this can cause the discrete/multicat macros to give the wrong answers
if `toplevel' < 3 local toplevel = 3
}
if "`numwscolumns'" == "" local numwscolumns 1500
if ("`verbose'"~="") display as text "Expanded variables: Count: `numexpand', Length: `longN' Estimated sheet size: `sheetsize'"
qui compress
* Assert that the data is sorted according to the model hierarchy
if "`cc'"=="" && "`sort'" == "" {
if "`mtype'" == "univariate" local bottomlev 1
else local bottomlev 0
if `maxlevels'>`bottomlev' {
local idvars2andhigher
local levelused = 0
forvalues l = `maxlevels'(-1)`=`bottomlev' + 1' {
if `numrp`l'vars' > 0 local levelused = 1
if `levelused' == 1 local idvars2andhigher `idvars2andhigher' `lev`l'id'
}
* Check whether correct sort order is set, otherwise examine the data
if `=strpos("`:sortedby'", "`idvars2andhigher'")' ~= 1 {
tempvar sortorder1
tempvar sortorder2
quietly gen `sortorder1' = _n
sort `idvars2andhigher' `sortorder1'
quietly gen `sortorder2' = _n
capture assert `sortorder2'==`sortorder1'
if _rc {
if "`forcesort'" ~= "" {
display as error "Warning: The data have been resorted in the order of the model hierarchy: `idvars2andhigher' `lev1id'." _n
}
else {
display as error "The data must be sorted according to the order of the model hierarchy: `idvars2andhigher' `lev1id'." _n
exit 198
}
}
drop `sortorder1' `sortorder2'
}
}
}
foreach var of varlist * {
if "`:type `var''" == "double" & `doublevar' == 0 {
if "`forcerecast'" ~= "" {
recast float `var', force
display as error "Warning: `var' has been recast to float, this has reduced precision"
}
else {
display as error in smcl "`var' is held to more precision than MLwiN can handle, to reduce the precision use {stata recast float `var', force}"
exit 198
}
}
if "`:type `var''" == "long" {
if `doublevar' == 0 {
capture recast float `var'
}
else {
capture recast double `var'
}
if `r(N)' > 0 {
if "`forcerecast'" ~= "" {
if `doublevar' == 0 {
recast float `var', force
display as error "Warning: `var' has been recast to float, if this is an ID variable its meaning may have changed"
}
else {
recast double `var', force
display as error "Warning: `var' has been recast to double, if this is an ID variable its meaning may have changed"
}
}
else {
display as error "`var' is held to more precision than MLwiN can handle, it must be recoded so that values lie in the range +/- 16,777,215"
exit 198
}
}
}
}
tempfile filedata
qui saveold "`filedata'", replace
******************************************************************************
* (4) WRITE MLWIN SUB MACRO TO FIT MODEL
******************************************************************************
tempname macro1
if "`savemacro'" ~= "" local filemacro1 `savemacro'
else tempfile filemacro1
if "`batch'" ~= "" local pause nopause
qui file open `macro1' using "`filemacro1'", write replace
if ("`pause'"=="nopause" | "`plugin'" ~= "") file write `macro1' "ECHO 0" _n
* Give details about runmlwin and date and time stamp the macro
file write `macro1' "NOTE ***********************************************************************" _n
file write `macro1' "NOTE MLwiN macro created by runmlwin Stata command: `c(current_date)', `c(current_time)'" _n
file write `macro1' "NOTE See: www.bristol.ac.uk/cmm/runmlwin for help" _n
file write `macro1' "NOTE ***********************************************************************" _n(2)
if "`verbose'" ~= "" display as text "Number of levels: `toplevel'"
if "`verbose'" ~= "" file write `macro1' "ECHO 1" _n
file write `macro1' "NOTE Initialise MLwiN storage" _n "INIT `toplevel' 10000 `numwscolumns' 150 30" _n(2)
file write `macro1' "MARK 0" _n
if "`batch'" ~= "" {
file write `macro1' "NOTE divert error messages to a file" _n "ERRO 0" _n
tempfile errlog
file write `macro1' "LOGO '`errlog'' 1" _n(2)
}
* Open the equations window
if ("`pause'"=="") {
file write `macro1' "NOTE Open the equations window" _n
file write `macro1' "WSET 15 1" _n // open the equations window
file write `macro1' "EXPA 2" _n // expand from response equation to include first model equation then random effects equations
file write `macro1' "ESTM 1" _n(2) // switch from black symbols to coloured symbols (but not point estimates)
}
if "`useworksheet'" ~= "" {
file write `macro1' "NOTE Load previously saved model state into MLwiN" _n "ZRET '`useworksheet''" _n(2)
}
else {
* Load the data set
file write `macro1' "NOTE Import the Stata data set into MLwiN" _n "RSTA '`filedata''" _n(2)
if "`rngversion'" ~= "" file write `macro1' "NOTE Set the random number generator version" _n "RNGV `rngversion'" _n
* Load design vectors if specified
forvalues l = `maxlevels'(-1)1 {
if "`design`l''"~="" {
file write `macro1' "NOTE Load the design matrix for the random part at level `l'" _n
file write `macro1' "LINK `numdesign`l'' G21" _n "GSTA '`filedesign`l''' G21" _n "LINK 0 G21" _n(2)
}
}
* Convergence tolerance
if "`iglstolerance'" ~= "" file write `macro1' "NOTE Specify the (R)IGLS convergence tolerance" _n "TOLE `iglstolerance'" _n(2)
* Specify the univariate response variable
file write `macro1' "NOTE Specify the response variable(s)" _n
if ("`mtype'"~="multivariate" & "`mtype'"~="multinomial") file write `macro1' "RESP '`response1''" _n
* Specify the multivariate response variable
local responsenames
forvalues e = `nummlwineqns'(-1)1 {
if "response`e'"~="" local responsenames '`response`e''' `responsenames'
}
if ("`mtype'"=="multivariate") file write `macro1' "MVAR 1 `responsenames'" _n
* Specify the multinomial response variable
// assert that response is categorical
if ("`distribution'"=="multinomial") {
local cresp = `numcolumns' + 1
local cresp_indicator = `numcolumns' + 2
* Expand the data
file write `macro1' "NAME c`cresp' 'resp' c`cresp_indicator' 'resp_indicator'" _n
if "`link'" == "mlogit" file write `macro1' "MNOM 0 '`response1'' 'resp' 'resp_indicator' `basecategory'" _n
else file write `macro1' "MNOM 1 '`response1'' 'resp' 'resp_indicator' `basecategory'" _n
file write `macro1' "RESP 'resp'" _n
if ("`link'"=="mlogit") file write `macro1' "RDIS 1 4" _n
if ("`link'"=="ologit"|"`link'"=="oprobit"|"`link'"=="ocloglog") file write `macro1' "RDIS 1 5" _n
if "`link'"=="ologit" | "`link'"=="mlogit" file write `macro1' "LFUN 0" _n
if "`link'"=="oprobit" file write `macro1' "LFUN 1" _n
if "`link'"=="ocloglog" file write `macro1' "LFUN 2" _n
file write `macro1' "DOFF 1 '`denominator''" _n
if ("`extra'" ~= "") file write `macro1' "EXTR 1" _n
}
if "`proportion'" ~= "" file write `macro1' "VECT `=`:list sizeof responsecats'-1' `propvarnames' 'resp'" _n
if ("`pause'"=="") file write `macro1' "PAUS 1" _n
file write `macro1' _n
* Loop over responses (equations)
forvalues e = 1/`numfpbrackets' {
* Binary response
// assert that response lies in the interval [0,1]
if ("`distribution`e''"=="binomial") {
file write `macro1' "NOTE Response distribution, link function and denominator" _n
file write `macro1' "RDIS `e' 0" _n
if ("`link'"=="probit") file write `macro1' "LFUN 1" _n
if ("`link'"=="cloglog") file write `macro1' "LFUN 2" _n
if ("`link'"=="logit" | "`link'"=="") file write `macro1' "LFUN 0" _n
file write `macro1' "DOFF `e' '`denominator`e'''" _n
if ("`extra'" ~= "") file write `macro1' "EXTR 1" _n
if ("`pause'"=="") file write `macro1' "PAUS 1" _n
file write `macro1' _n
}
* Poison response
// assert that response is positive integer?
if ("`distribution'"=="poisson") {
file write `macro1' "NOTE Response distribution, link function and denominator" _n
file write `macro1' "RDIS `e' 1" _n "LFUN 3" _n
if ("`extra'" ~= "") file write `macro1' "EXTR 1" _n
if ("`offset'"~="") file write `macro1' "DOFF `e' '`offset''" _n
if ("`pause'"=="") file write `macro1' "PAUS 1" _n
file write `macro1' _n
}
* Negative binomial response
// assert that response is positive integer?
if ("`distribution'"=="nbinomial") {
file write `macro1' "NOTE Response distribution, link function and denominator" _n
file write `macro1' "RDIS `e' 2" _n "LFUN 3" _n
if ("`extra'" ~= "") file write `macro1' "EXTR 1" _n
if ("`offset'"~="") file write `macro1' "DOFF `e' '`offset''" _n
if ("`pause'"=="") file write `macro1' "PAUS 1" _n
file write `macro1' _n
}
}
* Levels/Classification identifiers
file write `macro1' "NOTE Specify the level identifier(s)" _n
forvalues l = `maxlevels'(-1)1 {
local ll = `l'
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") local ++ll
if ("`lev`l'id'"~="") file write `macro1' "IDEN `ll' '`lev`l'id''" _n
}
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") file write `macro1' "IDEN 1 'resp_indicator'" _n
if ("`pause'"=="") file write `macro1' "PAUS 1" _n
file write `macro1' _n
* Specify covariate(s) used anywhere in the model
file write `macro1' "NOTE Specify covariate(s) used anywhere in the model" _n
local allvars :list uniq fp
forvalues l = `maxlevels'(-1)1 {
local allvars :list allvars|rp`l'
}
* Remove duplicates (for example in Multivariate)
local allvars :list uniq allvars
if ("`mtype'"=="multinomial" | "`mtype'"=="multivariate") {
forvalues e = 1/`numfpbrackets' {
if "`mtype'"=="multivariate" & "`=word("`distribution'", `e')'" == "binomial" {
// Take out and put back in bcons.* so they appear in the correct position in the terms
file write `macro1' "NOTE Ensure bcons is in the correct place for output" _n "NEXP 0 'bcons.`e''" _n "ERAS 'bcons.`e''" _n
}
if "`fp`e''" ~= "" {
file write `macro1' "RPAT `rpat`e''" _n
foreach var of varlist `fp`e''{
file write `macro1' "ADDT '`var''" _n
}
file write `macro1' "RPAT" _n
}
if `e' <= `nummlwineqns' {
forvalues l = `maxlevels'(-1)1 {
if "`rp`l'_`e''" ~= "" {
file write `macro1' "RPAT `rppat`l'_`e'' " _n
forvalues n = 1/`:list sizeof rp`l'name`e'' {
if `:list posof "`=word("`rp`l'name`e''", `n')'" in fpname' == 0 {
file write `macro1' "ADDT '`=word("`rp`l'_`e''", `n')''" _n "FPAR 0 '`=word("`rp`l'name`e''", `n')''" _n
}
}
file write `macro1' "RPAT" _n
}
}
}
}
forvalues l = `maxlevels'(-1)1 {
forvalues e = `=`nummlwineqns'+1'/`numrpbrackets`l'' {
if "`rp`l'_`e''" ~= "" {
file write `macro1' "RPAT `rppat`l'_`e'' " _n
forvalues n = 1/`:list sizeof rp`l'name`e'' {
if `:list posof "`=word("`rp`l'name`e''", `n')'" in fpname' == 0 {
file write `macro1' "ADDT '`=word("`rp`l'_`e''", `n')''" _n "FPAR 0 '`=word("`rp`l'name`e''", `n')''" _n
}
}
file write `macro1' "RPAT" _n
}
}
}
}
if ("`mtype'"=="univariate") {
forvalues e = 1/`numfpbrackets' {
foreach var of varlist `fp`e''{
file write `macro1' "ADDT '`var''" _n
}
}
* Add to the model any variables in the random part which are not in the fixed part
* Then remove from the fixed part of the model (but not entirely from the model itself!) any variables which are only in the random part
local nonfp :list allvars - fp
if "`nonfp'"~="" {
foreach var of varlist `nonfp' {
file write `macro1' "ADDT '`var''" _n
}
foreach var of varlist `nonfp' {
file write `macro1' "FPAR 0 '`var''" _n
}
}
}
if ("`pause'"=="") file write `macro1' "PAUS 1" _n
file write `macro1' _n
forvalues l = `maxlevels'(-1)1 {
local ll = `l'
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") local ++ll
* Level `l' random part covariates
if "`rp`l''"~="" {
file write `macro1' "NOTE Specify level `l' random part covariate(s)" _n
forvalues e = 1/`numrpbrackets`l'' {
if "`rp`l'_`e''" != "" {
foreach var of local rp`l'name`e' {
if ("`diagonal`l''"=="" & "`elements`l''"=="") {
file write `macro1' "SETV `ll' '`var''" _n
}
if ("`diagonal`l''"~="") {
file write `macro1' "SETE `ll' '`var'' '`var''" _n
}
}
}
}
* Add/Remove terms from variance covariance matrices
if "`elements`l''"~="" {
local elenames :colnames `RP`l'_b'
foreach p of local rp`l'pars { // for each var specified as random at level `l'
if regexm("`p'", "(cov)[\(]([a-zA-Z0-9_~]+)[\\]([a-zA-Z0-9_~]+)[\)]") == 1 {
local partype `=regexs(1)'
local rowvar `=regexs(2)'
local colvar `=regexs(3)'
capture unab rowvar: `rowvar'
capture unab colvar: `colvar'
}
if regexm("`p'", "(var)[\(]([a-zA-Z0-9_]+)[\)]") == 1 {
local partype `=regexs(1)'
local rowvar `=regexs(2)'
local colvar `=regexs(2)'
capture unab rowvar: `rowvar'
capture unab colvar: `colvar'
}
assert inlist("`partype'","var","cov")
if `:list p in elenames' {
if "`mtype'" == "multivariate" {
local rowvar = reverse(subinstr(reverse("`rowvar'"), "_", ".", 1))
local colvar = reverse(subinstr(reverse("`colvar'"), "_", ".", 1))
}
file write `macro1' "SETE `ll' '`rowvar'' '`colvar''" _n
}
}
}
if ("`pause'"=="") file write `macro1' "PAUS 1" _n
file write `macro1' _n
}
}
* Apply random part design matrix if set
forvalues l = `maxlevels'(-1)1 {
if "`design`l''"~="" {
file write `macro1' "NOTE Apply the design matrix for the random part at level `l'" _n
foreach alphaname of local design`l' {
file write `macro1' "SETD `l' '`alphaname''" _n
}
file write `macro1' _n
}
}
* Need to fit the ordered model for at least two iterations before we can impose starting values
if ("`mcmc'"~="" & "`cc'" == "") | ("`link'"=="ologit" | "`link'"=="oprobit" | "`link'"=="ocloglog" | "`link'"=="mlogit") { // CMJC added check for mlogit to get model in the correct state for lowest level residuals (e.g. adds P/-P parameter which will then be removed later)
file write `macro1' "NOTE Fit model for two iterations to correctly setup the model (this is only required for the multilevel ordered multinomial model)" _n
if "`linearization'"=="MQL1" file write `macro1' "LINE 0 1" _n
if "`linearization'"=="MQL2" file write `macro1' "LINE 0 2" _n
if "`linearization'"=="PQL1" file write `macro1' "LINE 1 1" _n
if "`linearization'"=="PQL2" file write `macro1' "LINE 1 2" _n
file write `macro1' "MAXI 2" _n "BATCH 1" _n "STAR" _n(2)
// set the maximum iterations back to the default if not set
if "`maxiterations'" == "" local maxiterations 20
local modelstarted 1
}
* Fixed part constraint(s)
if "`fpconstraints'"~="" {
file write `macro1' "NOTE Specify fixed part constraint(s)" _n
file write `macro1' "FCON b1001" _n
forvalues r = 1/`numfpconstraints' {
file write `macro1' "JOIN cb1001 `fpconstraint`r'' cb1001" _n
}
}
* Random part constraint(s)
if "`rpconstraints'"~="" {
file write `macro1' "NOTE Specify random part constraint(s)" _n
file write `macro1' "RCON b1002" _n
forvalues r = 1/`numrpconstraints' {
file write `macro1' "JOIN cb1002 `rpconstraint`r'' cb1002" _n
}
}
}
* Edit vector of initial values
if "`inits'"~="" {
runmlwin_setinit, fpb(`matFP_b') fpv(`matFP_V') rpb(`matRP_b') rpv(`matRP_V') initb(`matb') initv(`matV') mtype(`mtype') link(`link') mcmc(`mcmc') macro1(`macro1')
}
if "`inits'"~="" | "`useworksheet'" ~= "" {
file write `macro1' "ESTM 2" _n // switch from coloured symbols to point estimates
if ("`pause'"=="") file write `macro1' "PAUS 1" _n // update equations window
file write `macro1' _n
}
*************************
* Fit model using (R)IGLS
*************************
if "`mcmc'"=="" { // i.e. (R)IGLS
if "`useworksheet'" == "" {
* Specify estimation method (R)IGLS
if "`rigls'"~="" {
file write `macro1' "NOTE Set estimation method to be RIGLS" _n "METH 0" _n(2)
}
else {
file write `macro1' "NOTE Set estimation method to be IGLS" _n "METH 1" _n(2)
}
* Specify maximum number of iterations
if "`maxiterations'"~="" {
file write `macro1' "NOTE Set maximum number of (R)IGLS iterations" _n "MAXI `maxiterations'" _n(2)
}
// Why is MQL2 the default in MLwiN for poisson models?
// What about negbin is there the same problem?
if "`distribution'"=="poisson" & "`linearization'"=="MQL1" file write `macro1' "NOTE Set estimation method to be MQL1" _n "LINE 0 1" _n(2)
if "`linearization'"=="MQL2" file write `macro1' "NOTE Set estimation method to be MQL2" _n "LINE 0 2" _n(2)
if "`linearization'"=="PQL1" file write `macro1' "NOTE Set estimation method to be PQL1" _n "LINE 1 1" _n(2)
if "`linearization'"=="PQL2" file write `macro1' "NOTE Set estimation method to be PQL2" _n "LINE 1 2" _n(2)
* Allow negative variances
local reseton = 0
forvalues l = `maxlevels'(-1)1 {
if "`reset`l''"~="" local reseton = 1
}
if `reseton' == 1 {
file write `macro1' "NOTE Specify element reset options for variance-covariance matrices" _n
forvalues l = `maxlevels'(-1)1 {
local ll = `l'
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") local ++ll
if "`reset`l''"=="all" file write `macro1' "RESE `ll' 0" _n
if "`reset`l''"=="variances" file write `macro1' "RESE `ll' 1" _n
if "`reset`l''"=="none" file write `macro1' "RESE `ll' 2" _n
}
file write `macro1' "" _n
}
* Weights
forvalues l = `maxlevels'(-1)1 {
if "`weight`l''" ~= "" | "`weighttype'" ~= "" {
if "`weight`l''" ~= "" {
file write `macro1' "NOTE Specify sampling weights at level `l'" _n "WEIG `l' 1 '`weight`l'''" _n(2)
}
else {
file write `macro1' "NOTE Specify equal weights at level `l'" _n "WEIG `l' 1" _n(2)
}
if "`weighttype'" == "standardised" {
file write `macro1' "NOTE Standardised weighting" _n
file write `macro1' "LINK 1 G21" _n
file write `macro1' "FILL G21" _n
file write `macro1' "WEIG `l' 2 G21[1]" _n
file write `macro1' "LINK 0 G21" _n
file write `macro1' "WEIG 2" _n(2)
}
if "`weighttype'" == "raw" {
file write `macro1' "NOTE Raw weighting" _n
// MLwiN appears to require this to be set, even if raw weights are used (I think this is a bug)
file write `macro1' "LINK 1 G21" _n
file write `macro1' "FILL G21" _n
file write `macro1' "WEIG `l' 2 G21[1]" _n
file write `macro1' "LINK 0 G21" _n
file write `macro1' "WEIG 1" _n(2)
}
}
}
if "`weighttype'"~="" { // if weights have been specified
if "`weighttype'"=="standardised" {
file write `macro1' "NOTE Create the standardised weights" _n "WEIG" _n(2)
}
}
if "`fpsandwich'"=="fpsandwich" | ("`weighttype'"~="" & "`fpsandwich'"==""){
file write `macro1' "NOTE Turn on sandwich estimators for the fixed part parameter standard errors" _n "FSDE 2" _n(2)
}
if "`rpsandwich'"=="rpsandwich" | ("`weighttype'"~="" & "`rpsandwich'"==""){
file write `macro1' "NOTE Turn on sandwich estimators for the random part parameter standard errors" _n "RSDE 2" _n(2)
}
}
* Pause
if ("`pause'"=="") {
file write `macro1' "NOTE Pause the macro to allow the user to examine the model specification" _n "PAUS" _n(2)
}
* Fit the model
file write `macro1' "NOTE Fit the model" _n
if "`inits'"=="" & "`useworksheet'" == "" & "`modelstarted'" == "" file write `macro1' "STAR" _n // do one iteration as batch 1 not yet set (this is to get sensible numbers for the equations window)
else {
if ("`initsprevious'"~="") display as text "Model fitted using initial values specified as parameter estimates from previous model" _n
if ("`initsmodel'"~="") display as text "Model fitted using initial values specified as parameter estimates from saved estimates in `initsmodel'" _n
if ("`initsb'"~=""|"`initsv'"~="") display as text "Model fitted using initial values specified in matrices: " as result "`initsb' `initsv'" _n
if ("`useworksheet'"~="") display as text "Model fitted using state from previously saved worksheet" _n
}
if ("`pause'"=="") {
file write `macro1' "ESTM 2" _n // switch from coloured symbols to point estimates
file write `macro1' "PAUS 1" _n // update equations window
}
file write `macro1' "BATC 1" _n
if ("`inits'"~="" & ("`link'"=="ologit"|"`link'"=="oprobit"|"`link'"=="ocloglog")) file write `macro1' "SUMM" _n // This is a fix as the ordered model forgets its structure when initial values have been specified
file write `macro1' "NEXT" _n
file write `macro1' "ITNU 0 b21" _n
file write `macro1' "CONV b22" _n
if ("`pause'"=="") file write `macro1' "PAUS 1" _n
file write `macro1' _n
}
**********************
* Fit model using MCMC
**********************
if ("`mcmc'"~="") { // i.e. MCMC
if "`useworksheet'" == "" {
* Model specification: Cross-classifications
file write `macro1' "NOTE Set estimation method to MCMC" _n
file write `macro1' "EMODe 3" _n
if ("`pause'"=="") file write `macro1' "EXPA 3" _n // Expand to show priors
file write `macro1' _n
* MM Weights
if "`carcentre'" ~= "" file write `macro1' "CARC 1" _n
forvalues l = 2/`maxlevels' {
local ll = `l'
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") local ++ll
if "`mmweights`l''" ~= "" {
file write `macro1' "NOTE Declare weights at level `l'" _n
if "`car`l''"~="" {
file write `macro1' "MULM `ll' `:list sizeof mmweights`l'' '`=word("`mmweights`l''", 1)'' '`=word("`mmids`l''", 1)''" _n "CARP `ll' 1" _n
}
else {
if "`=word("`mmids`l''", 1)'" ~= "`lev`l'id'" {
display as error "The first var in mmids`l'() must match the level `l' ID specified in level`l'()."
exit 198
}
if "`mtype'"~="multivariate" {
file write `macro1' "MULM `ll' `:list sizeof mmweights`l'' '`=word("`mmweights`l''", 1)''" _n
}
else {
file write `macro1' "MULM `ll' `:list sizeof mmweights`l'' '`=word("`mmweights`l''", 1)'' '`=word("`mmids`l''", 1)''" _n
}
}
}
}
* Cross-classified models
if ("`cc'"~="") {
file write `macro1' "NOTE Declare model to be a cross-classifed model and switch to classification notation" _n "XCLA 1" _n /* switch on cross-classifications */ "INDE 1" _n(2) /* switch on classification notation */
}
* Define measurement error if defined
if "`me'"~="" {
* measurement error
file write `macro1' "NOTE Specify measurement error" _n
local mecols
local i = 1
foreach p of local mevars {
local mecols `mecols' '`p'' `:word `i' of `variances''
local ++i
}
file write `macro1' "MERR `:list sizeof mevars' `mecols'" _n
file write `macro1' _n
}
* Correlated residuals
if "`corresiduals'" ~= "" {
file write `macro1' "NOTE Set up correlated residuals" _n
if "`corresiduals'" == "exchangeable" file write `macro1' "NOTE all correlations equal and all variances equal" _n "MCCO 1" _n
if "`corresiduals'" == "areqvars" file write `macro1' "NOTE an AR1 structure with all variances equal" _n "MCCO 2" _n
if "`corresiduals'" == "eqcorrsindepvars" file write `macro1' "NOTE all correlations equal but independent variances" _n "MCCO 3" _n
if "`corresiduals'" == "arindepvars" file write `macro1' "NOTE an AR1 structure with independent variances" _n "MCCO 4" _n
}
* Use multiple subscripts notation
if ("`msubscripts'"~="") {
file write `macro1' "NOTE Use multiple subscripts notation" _n
file write `macro1' "INDE 1" _n(2)
}
* Set MCMC seed
file write `macro1' "NOTE Set MCMC seed" _n "MCRS `mcmc_seed'" _n(2)
if "`smcmc'" ~= "" {
file write `macro1' "NOTE Turn on structured MCMC" _n
file write `macro1' "SMCM 1" _n(2)
}
if "`smvn'" ~= "" {
file write `macro1' "NOTE Turn on structured MVN" _n
file write `macro1' "SMVN 1" _n(2)
}
if "`orthogonal'" ~= "" {
file write `macro1' "NOTE Turn on orthogonal parameterisation" _n
file write `macro1' "ORTH 1" _n(2)
}
if "`hcentring'" ~= "" {
file write `macro1' "NOTE Turn on hierarchical centring" _n
file write `macro1' "HCEN 1 `hcentring'" _n(2)
}
if "`logformulation'" ~= "" {
file write `macro1' "NOTE Use log formulation" _n
file write `macro1' "LCLO 1" _n(2)
}
if `defaultprior' == 1 {
file write `macro1' "NOTE Set prior distribution parameters" _n
file write `macro1' "GAMP `gama' `gamb'" _n(2)
}
forvalues l = 2/`maxlevels' {
local ll = `l'
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") local ++ll
if "`parexpansion`l''" ~= "" {
file write `macro1' "NOTE Turn on parameter expansion at level `l'" _n
file write `macro1' "PAEX `ll' 1" _n(2)
}
}
* Informative priors
if "`priormatrix'" ~= "" {
local priorcol c1092
runmlwin_writepriors, maxlevels(`maxlevels') priormat(`priormatrix') fpb(`matFP_b') rpb(`matRP_b') macro1(`macro1') priorcol(`priorcol')
}
* Add factors if specified
local numfact = 0
local numcorr = 0
forvalues l = 1/`maxlevels' {
if "`flinit`l''" ~= "" {
local numfact = `numfact' + colsof(`flinit`l'')
forvalues i = 1/`=colsof(`flinit`l'')' {
forvalues j = 1/`=`i'-1' {
if `fvinit`l''[`i',`j'] ~= . {
local ++numcorr
}
}
}
}
}
local factstr ""
forvalues l = 1/`maxlevels' {
if "`flinit`l''" ~= "" {
forvalues fact = 1/`=colsof(`flinit`l'')' {
local factstr `factstr' `=`l'+1'
forvalues resp = 1/`=rowsof(`flinit`l'')' {
local factstr `factstr' `=`flinit`l''[`resp',`fact']' `=`flconstraint`l''[`resp', `fact']'
}
local factstr `factstr' `=`fvinit`l''[`fact', `fact']' `=`fvconstraint`l''[`fact', `fact']'
}
}
}
forvalues l = 1/`maxlevels' {
if "`flinit`l''" ~= "" {
forvalues fact1 = 1/`=colsof(`flinit`l'')' {
forvalues fact2 = 1/`=`fact1'-1' {
if `fvinit`l''[`fact1', `fact2'] ~= . {
local factstr `factstr' `fact1' `fact2' `=`fvinit`l''[`fact1', `fact2']' `=`fvconstraint`l''[`fact1', `fact2']'
}
}
}
}
}
if "`factstr'" ~= "" {
file write `macro1' "NOTE set up factors" _n
file write `macro1' "FACT `nummlwineqns' `numfact' `numcorr' `factstr'" _n
}
local useresid "1" // Use IGLS residuals as starting values
local residcol
local residsecol
// Simulate residuals from the covariance matrix
if "`simresiduals'" ~= "" {
file write `macro1' "LINK 2 G27" _n // Starting residuals columns
file write `macro1' "FILL G27" _n
local hasres = 0
local residcol G27[1]
local residsecol G27[2]
local useresid = ""
forvalues l = 1/`maxlevels' {
local ll = `l'
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") local ++ll
file write `macro1' "LINK 1 G22" _n
file write `macro1' "LINK `numresi`l'' G23" _n
if "`numrp`l'pars'" ~= "" {
file write `macro1' "OMEGa `ll' G22[1]" _n
file write `macro1' "CALC G22[1] = HSYM(G22[1])" _n
file write `macro1' "NOBS `ll' b31 b32" _n
file write `macro1' "MRAN b31 G22[1] G23" _n
forvalues i = 1/`numresi`l'' {
if `hasres' == 0 { // remove zero that fill command has put into this column
file write `macro1' "ERAS G27" _n
local hasres = 1
}
file write `macro1' "JOIN G27[1] G23 G27[1]" _n
}
}
file write `macro1' "ERASe G22" _n
file write `macro1' "ERASe G23" _n
file write `macro1' "LINK 0 G22" _n "LINK 0 G23" _n
}
file write `macro1' "CALC G27[2] = G27[1] * 0" _n
}
if "`useresid'" ~= "" & "`cc'" == "" & (`maxlevels' > 1 | "`mtype'" == "multivariate" | ("`mtype'" == "multinomial" & "`link'"=="mlogit")){
* Calculate MCMC starting values for level `l' residuals
if `allnormal' ~= 1 {
if "`mtype'" == "multinomial" {
file write `macro1' "SETV 1 'bcons.1'" _n
}
if "`linearization'" == "PQL2" & "`maxlevels'" == "1" {
forvalues lev = `maxlevels'/5 {
file write `macro1' "OFFS `lev'" _n
}
display as error "can not apply PQL or 2nd order to single level models, switching to MQL/1"
local linearization MQL1
}
// NOTE at this point MLwiN sets EXTR 0, LINE 0 1 if these are not set
if "`linearization'" == "" | "`linearization'" == "MQL1" file write `macro1' "LINE 0 1" _n // cmjc added for poisson
if "`linearization'"=="MQL2" file write `macro1' "LINE 0 2" _n
if "`linearization'"=="PQL1" file write `macro1' "LINE 1 1" _n
if "`linearization'"=="PQL2" file write `macro1' "LINE 1 2" _n
// Extra is not valid for MCMC estimation
file write `macro1' "EXTR 0" _n // cmjc added for poisson
local llev = 1
if "`mtype'"=="multivariate" | "`mtype'"=="multinomial" local ++llev
forvalues lev = `llev'/`maxlevels' {
file write `macro1' "OFFS `lev'" _n
}
if "`mtype'" ~= "multivariate" { // i.e. univariate response models
local b
foreach var of local rp1name {
if "`var'" ~= "P" & "`var'" ~= "-P" local b `b' '`var''
}
if "`distribution'" == "nbinomial" local b 'bcons.1' 'bcons2.1'
else local b 'bcons.1'
if `:list sizeof b' == 0 & "`mtype'" ~= "multinomial" display as error "Level 1 random structure must be specified for discrete response models"
if `:list sizeof b' ~= 0 file write `macro1' "LINK `b' g9" _n
if "`mtype'"=="multinomial" & ~inlist("`link'","mlogit","ologit","oprobit","ocloglog") {
display as error "link function must be mlogit, ologit, oprobit or ocloglog"
exit 198
}
if "`mtype'"=="multinomial" {
if "`link'"=="mlogit" file write `macro1' "SET b10 0" _n
if "`link'"=="ologit"|"`link'"=="oprobit"|"`link'"=="ocloglog" file write `macro1' "SET b10 1" _n
}
else file write `macro1' "SET b10 `resptype'" _n // i.e. a univariate response which is not multinomial
file write `macro1' "SET b15 1" _n
file write `macro1' "SET b16 0" _n
}
else { // i.e. multivariate response models
local b
local c
forvalues eq = 1/`nummlwineqns' {
foreach var of local rp1name`eq' {
if "`=word("`distribution'", `eq')'" ~= "normal" local b `b' '`var''
else local c `c' '`var''
}
foreach var of local fpname`eq' {
if "`=word("`distribution'", `eq')'" == "normal" local c `c' '`var''
}
}
if "`b'" == "" display as error Level 2 random structure must be specified for discrete responses"
file write `macro1' "LINK `b' g9" _n
file write `macro1' "SET b10 `resptype'" _n
file write `macro1' "SET b15 2" _n
if `mixed' == 0 file write `macro1' "SET b16 0" _n
else {
file write `macro1' "SET b16 1" _n
file write `macro1' "LINK `c' g11" _n
}
}
// LINK (backup/restore of model columns before residuals is no longer needed in MLwiN 3.13 and later)
file write `macro1' "LINK -5 g19 1" _n
file write `macro1' "LINK 0 g18" _n
file write `macro1' "CALC g18=g19" _n
file write `macro1' "ITNU 1 2" _n // Set R(IGLS) iteration to 2 so that residuals are calculated correctly
}
file write `macro1' "MISR 0" _n // Don't set small residuals to system missing
file write `macro1' "LINK 2 G27" _n // Starting residuals columns
file write `macro1' "FILL G27" _n
local hasres = 0
local residcol G27[1]
local residsecol G27[2]
forvalues l = 1/`maxlevels' {
local ll = `l'
if "`mtype'"=="multivariate" | "`mtype'"=="multinomial" local ++ll
local numresi = `numrp`l'vars'
if "`mtype'"=="multinomial" & "`link'"=="mlogit" & `l' == 1 local numresi = `numresi' + 2
if (`numresi'>=1 & `ll' != 1) {
file write `macro1' "LINK `numresi' G22" _n
file write `macro1' "FILL G22" _n
file write `macro1' "LINK 1 G23" _n
file write `macro1' "FILL G23" _n
file write `macro1' "NOTE Calculate MCMC starting values for level `l' residuals" _n
file write `macro1' "RLEV `ll'" _n
file write `macro1' "RFUN" _n
file write `macro1' "RCOV 2" _n
file write `macro1' "ROUT G22 G23" _n
file write `macro1' "RESI" _n
if `hasres' == 0 { // remove zero that fill command has put into this column
file write `macro1' "ERAS G27" _n
local hasres 1
}
file write `macro1' "JOIN G27[1] G22 G27[1]" _n
file write `macro1' "JOIN G27[2] G23 G27[2]" _n
file write `macro1' "ERAS G22 G23" _n
file write `macro1' "LINK 0 G22" _n "LINK 0 G23" _n(2)
}
}
file write `macro1' "MISR 1" _n
if `allnormal' != 1 {
// UNLINK (backup/restore of model columns before residuals is no longer needed in MLwiN 3.13 and later)
file write `macro1' "CALC g19=g18" _n
file write `macro1' "ERASe g18" _n
file write `macro1' "LINK 0 g18" _n
}
}
}
* Pause
// AT THIS POINT THE MULTIVARIATE WISHART DATA INFORMED PRIOR FOR THE LEVEL 2 MATRIX OF A TWO-LEVEL RANDOM SLOPES MODEL HAS NOT BEEN SET
if ("`pause'"=="") {
file write `macro1' "NOTE Pause the macro to allow the user to examine the model specification" _n
file write `macro1' "PAUS" _n(2)
}
* MCMC Estimation algorithm
if ("`distribution'"=="normal") file write `macro1' "NOTE Fit the model in MCMC using Gibbs algorithm (response is normal)" _n // But probit is fit with GIBBS
if ("`distribution'"~="normal") file write `macro1' "NOTE Fit the model in MCMC using Metropolis Hastings algorithm (response is discrete)" _n
if ("`savewinbugs'"~="") {
local 0 , `savewinbugs'
syntax , Model(string) Inits(string) Data(string) [noFit]
local 0 `inits'
syntax anything(name=inits), [REPLACE]
if "`replace'" == "" confirm new file "`inits'"
local 0 `data'
syntax anything(name=data), [REPLACE]
if "`replace'" == "" confirm new file "`data'"
local 0 `model'
syntax anything(name=model), [REPLACE]
if "`replace'" == "" confirm new file "`model'"
local linkcode 0 // needs to be set (0 for logit, 1 for probit, 2 for complementary log-log and 3 for log link)
if "`link'" == "logit" local linkcode 0
if "`link'" == "probit" local linkcode 1
if "`link'" == "cloglog" local linkcode 2
if "`link'" == "log" local linkcode 3
if "`mtype'" == "multivariate" local linkcode 0
if "`mcmc_modeltype'" == "5"{
display as error "WinBUGS in unable to estimate multivariate models with responses other than normal"
exit 198
}
file write `macro1' "BUGO 6 `mcmc_modeltype' `linkcode' `residcol' `priorcol' '`model'' '`inits'' '`data''" _n // MCMC chains
}
local chaincols 0
forvalues l = 1/`maxlevels' {
if "`saveresiduals`l''" ~= "" local ++chaincols
}
if "`factstr'" ~= "" local chaincols = `chaincols' + 2
if `chaincols' > 0 file write `macro1' "LINK `chaincols' G22" _n
local i = 1
forvalues l = 1/`maxlevels' {
local ll = `l'
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") local ++ll
if "`saveresiduals`l''" ~= "" {
file write `macro1' "NOTE store residual chain for level `l'" _n
file write `macro1' "NAME G22[`i'] 'RESID_CHAIN_LEVEL`l''" _n
file write `macro1' "DESC G22[`i'] '\Value'" _n
file write `macro1' "SMRE `ll' 'RESID_CHAIN_LEVEL`l''" _n(2)
local ++i
}
}
if "`factstr'" ~= "" {
file write `macro1' "NAME G22[`i'] '_FACT_LOAD_CHAIN'" _n
local ++i
file write `macro1' "NAME G22[`i'] '_FACT_VAR_CHAIN'" _n
local ++i
file write `macro1' "SMFA 1 '_FACT_LOAD_CHAIN'" _n
file write `macro1' "SMFA 2 '_FACT_VAR_CHAIN'" _n
file write `macro1' "LINK 0 G22" _n(2)
}
if "`fit'" ~= "nofit" {
if "`useworksheet'" == "" {
if "`mtype'" == "multinomial" file write `macro1' "CLRV 2" _n // I'm not sure why this is needed
file write `macro1' "MTOT `mcmc_chain'" _n
file write `macro1' "MCMC 0 `mcmc_burnin' `mcmc_adaptation' `mcmc_scale' `mcmc_acceptance' `mcmc_tolerance' `residcol' `residsecol' `priorcol' `mcmc_femethod' `mcmc_remethod' `mcmc_levelonevarmethod' `mcmc_higherlevelsvarmethod' `mcmc_defaultprior' `mcmc_modeltype'" _n
if "`residcol'" ~= "" file write `macro1' "ERAS G27" _n
if "`residcol'" ~= "" file write `macro1' "LINK 0 G27" _n
}
local total_stored = ceil(`mcmc_chain' / `mcmc_thinning')
if ("`batch'"=="") { // split MCMC chain into mcmc_refresh (stored) iterations at a time then update
file write `macro1' _n // (MREF command was introduced in v3.09, so conditions could be merged if this was the minimum version supported)
local curriter 0 // Unthinned iterations
if "`imputeiterations'" ~= "" {
local nextimpute `:word 1 of `imputeiterations''
local imputeiterations `:list imputeiterations - nextimpute'
}
while `curriter' + (`mcmc_refresh' * `mcmc_thinning') <= `mcmc_chain' {
local curriter = `curriter' + (`mcmc_refresh' * `mcmc_thinning')
if "`nextimpute'" ~= "" {
local impiter = `curriter' - (`mcmc_refresh' * `mcmc_thinning')
while (`nextimpute' > `curriter' - (`mcmc_refresh' * `mcmc_thinning')) & (`nextimpute' <= `curriter') {
local nextstop = (`nextimpute' / `mcmc_thinning') - ceil(`impiter' / `mcmc_thinning')
local impiter = `impiter' + (`nextstop' * `mcmc_thinning')
file write `macro1' "MCMC 1 `nextstop' `mcmc_thinning' c1090 c1091 c1003 c1004 1 `mcmc_modeltype'" _n
file write `macro1' "LINK 1 G21" _n
file write `macro1' "DAMI 0 G21[1]" _n // Impute
file write `macro1' "LINK `nummlwineqns' G22" _n
file write `macro1' "SPLI G21[1] 'resp_indicator' G22" _n
forvalues e = 1/`nummlwineqns' {
file write `macro1' "NAME G22[`e'] '_mi_`response`e'''" _n
}
tempfile impfile
file write `macro1' "PSTA '`impfile'' '`_sortindex'' G22" _n
local imputefiles `" `imputefiles' "`impfile'" "'
file write `macro1' "ERAS G21 G22" _n
file write `macro1' "LINK 0 G21" _n "LINK 0 G22" _n
file write `macro1' "PUPN c1003 c1004" _n
file write `macro1' "AVER c1091 b99 b100" _n
file write `macro1' "PAUS 1" _n(2)
local nextimpute `:word 1 of `imputeiterations''
local imputeiterations `:list imputeiterations - nextimpute'
if "`nextimpute'" == "" continue, break
}
local nextstop = ceil((`curriter' - `impiter') / `mcmc_thinning')
if `nextstop' ~= 0 {
file write `macro1' "MCMC 1 `nextstop' `mcmc_thinning' c1090 c1091 c1003 c1004 1 `mcmc_modeltype'" _n
file write `macro1' "PUPN c1003 c1004" _n "AVER c1091 b99 b100" _n "PAUS 1" _n(2)
}
}
else {
file write `macro1' "MCMC 1 `mcmc_refresh' `mcmc_thinning' c1090 c1091 c1003 c1004 1 `mcmc_modeltype'" _n
file write `macro1' "PUPN c1003 c1004" _n "AVER c1091 b99 b100" _n "PAUS 1" _n(2)
}
}
if `mcmc_chain' - `curriter' ~= 0 {
local remainder = ceil((`mcmc_chain' - `curriter') / `mcmc_thinning')
file write `macro1' "MCMC 1 `remainder' `mcmc_thinning' c1090 c1091 c1003 c1004 1 `mcmc_modeltype'" _n
file write `macro1' "PUPN c1003 c1004" _n "AVER c1091 b99 b100" _n "PAUS 1" _n(2)
}
}
else {
if "`imputeiterations'" ~= "" {
local curriter 0 // Thinned iterations
while "`imputeiterations'" ~= "" {
local nextimpute `:word 1 of `imputeiterations''
local imputeiterations `:list imputeiterations - nextimpute'
local nextstop = (`nextimpute' / `mcmc_thinning') - `curriter'
local curriter = `curriter' + `nextstop'
file write `macro1' "MCMC 1 `nextstop' `mcmc_thinning' c1090 c1091 c1003 c1004 1 `mcmc_modeltype'" _n(2)
file write `macro1' "LINK 1 G21" _n
file write `macro1' "DAMI 0 G21[1]" _n // Impute
file write `macro1' "LINK `nummlwineqns' G22" _n
file write `macro1' "SPLI G21[1] 'resp_indicator' G22" _n
forvalues e = 1/`nummlwineqns' {
file write `macro1' "NAME G22[`e'] 'imp`response`e'''" _n
}
tempfile impfile
file write `macro1' "PSTA '`impfile'' '`_sortindex'' G22" _n
local imputefiles `" `imputefiles' "`impfile'" "'
file write `macro1' "ERAS G21 G22" _n
file write `macro1' "LINK 0 G21" _n "LINK 0 G22" _n
file write `macro1' "PUPN c1003 c1004" _n
file write `macro1' "AVER c1091 b99 b100" _n(2)
}
local nextstop = `total_stored' - `curriter'
if `nextstop' ~= 0 {
file write `macro1' "MCMC 1 `nextstop' `mcmc_thinning' c1090 c1091 c1003 c1004 1 `mcmc_modeltype'" _n(2)
file write `macro1' "PUPN c1003 c1004" _n
file write `macro1' "AVER c1091 b99 b100" _n(2)
}
}
else {
file write `macro1' "MCMC 1 `total_stored' `mcmc_thinning' c1090 c1091 c1003 c1004 1 `mcmc_modeltype'" _n(2)
file write `macro1' "PUPN c1003 c1004" _n
file write `macro1' "AVER c1091 b99 b100" _n(2)
}
}
if ("`pause'"=="") file write `macro1' "PAUS 1" _n(2)
}
}
if "`fit'" ~= "nofit" {
if "`imputesummaries'" ~= "" {
file write `macro1' "LINK 2 G21" _n
file write `macro1' "DAMI 2 G21[1] G21[2]" _n // Impute
file write `macro1' "LINK `nummlwineqns' G22" _n
file write `macro1' "SPLI G21[1] 'resp_indicator' G22" _n
forvalues e = 1/`nummlwineqns' {
file write `macro1' "NAME G22[`e'] '_simean_`response`e'''" _n
}
file write `macro1' "LINK `nummlwineqns' G23" _n
file write `macro1' "SPLI G21[2] 'resp_indicator' G23" _n
forvalues e = 1/`nummlwineqns' {
file write `macro1' "NAME G23[`e'] '_sisd_`response`e'''" _n
}
tempfile imputesummaries
file write `macro1' "PSTA '`imputesummaries'' '`_sortindex'' G22 G23" _n
file write `macro1' "ERAS G21 G22 G23" _n
file write `macro1' "LINK 0 G21" _n "LINK 0 G22" _n "LINK 0 G23" _n
}
if "`batch'" == "" {
* Open the equations window
if ("`pause'"=="nopause") {
file write `macro1' "NOTE Open the equations window" _n
file write `macro1' "WSET 15 1" _n // open the equations window to indicate that the model has finished iterating but don't display anything
file write `macro1' "EXPA 3" _n
file write `macro1' "ESTM 2" _n
}
* If any random part design matrices have been specified then open the output window and display the random part parameters
if "`anydesign'"=="yes"{
file write `macro1' "WSET 20 1" _n
file write `macro1' "ECHO 1" _n
file write `macro1' "FIXE" _n
file write `macro1' "RAND" _n
file write `macro1' "ECHO 0" _n
}
}
*******************
* Calculate Residuals
*******************
forvalues l = `maxlevels'(-1)1 {
if ("`residuals`l''"~=""){
tempfile residfile`l'
local lids
forvalues lev = `maxlevels'(-1)1 {
local lids `lids' lev`lev'id(`lev`lev'id')
}
runmlwin_calcresiduals , level(`l') levxid(`lev`l'id') numrpxvars(`numresi`l'') rpxvars(`resivars`l'') residualsx(`residuals`l'') mtype(`mtype') `residstd`l'' `residlever`l'' `residinfl`l'' `residdel`l'' `residsamp`l'' `residref`l'' `residvar`l'' `residrecode`l'' `residrpadj`l'' `residnoadj`l'' macro(`macro1') mcmc(`mcmc') residfile(`residfile`l'') cc(`cc') maxlevels(`maxlevels') mmid(`mmids`l'') `lids'
}
}
*******************
* Calculate Factor values
*******************
if `factorson' == 1 {
file write `macro1' "LINK 3 G21" _n
file write `macro1' "JOIN 1 G21[3] G21[3]" _n
file write `macro1' "NOBS 2 b33 b34" _n
file write `macro1' "DAFA G21[1] G21[2]" _n
file write `macro1' "SET b35 1" _n
forvalues l = 1/`maxlevels' {
if ("`fscores`l''"~="") {
tempfile fscoresfile`l'
local count = `=colsof(`flinit`l'')' * 2
file write `macro1' "LINK `count' G22" _n
}
local i = 1
forvalues fact = 1/`=colsof(`flinit`l'')' { {
if ("`fscores`l''"~=""){
file write `macro1' "NOBS `=`l'+1' b31 b32" _n
file write `macro1' "NAME G22[`i'] '`fscores`l''`fact''" _n
file write `macro1' "CALC b36 = b35 + b31 - 1" _n
file write `macro1' "PICK b35 b36 G21[1] G22[`i']" _n
local ++i
file write `macro1' "NAME G22[`i'] '`fscores`l''`fact'sd'" _n
file write `macro1' "PICK b35 b36 G21[2] G22[`i']" _n
local ++i
}
file write `macro1' "CALC b35=b35+b33" _n
}
if ("`fscores`l''"~=""){
file write `macro1' "GENErate 1 b31 1 G21[3]" _n
file write `macro1' "NAME G21[3] '_fscoreid'" _n
file write `macro1' "PSTA '`fscoresfile`l''' '_fscoreid' G22" _n
file write `macro1' "ERASE '_fscoreid' G22" _n
file write `macro1' "LINK 0 G22" _n
}
}
file write `macro1' "ERASE G21" _n
file write `macro1' "LINK 0 G21" _n
}
**********************
* Set Seed
**********************
if ("`standardseed'"~="") {
// Note if user does not specify seed option then seed is set to zero and runmlwin does not write the seed command to the macro
// If the user specifies the seed option then runmlwin will write the seed command to the macro (unless that is if the user specifies the seed to be 0!)
file write `macro1' "NOTE Set seed" _n
file write `macro1' "SEED `standardseed'" _n(2)
}
**********************
* Simulate new response variable
**********************
if ("`simulate'"~="") {
tempfile simfile
runmlwin_simresp, simfile(`simfile') simulate(`simulate') simfile(`simfile') macro1(`macro1')
}
****************
* Save the model
****************
if "`saveworksheet'" ~= "" {
file write `macro1' "NOTE Save the MLwiN worksheet" _n
file write `macro1' "ZSAV '`saveworksheet''" _n(2)
}
* Pause
if ("`pause'"=="") {
file write `macro1' "NOTE Pause the macro to allow the user to examine the model results" _n
file write `macro1' "PAUS" _n(2)
}
if "`saveequation'" ~= "" {
file write `macro1' "NOTE Save an image of the equations window" _n
file write `macro1' "WCAP 15 0 '`saveequation''" _n
}
}
file write `macro1' "NOTE ***********************************************************************" _n(2)
file close `macro1'
if ("`viewmacro'"~="") view "`filemacro1'"
******************************************************************************
******************************************************************************
* MLwiN sub macro to export model results
******************************************************************************
******************************************************************************
// GENERATE COLUMNS TO HOLD MCMC CHAINS AND STATISTICS
// MCMC chains for parameters can be stored here
// MCMC chains for residuals can be stored here
tempname macro2
tempfile filemacro2
qui file open `macro2' using "`filemacro2'", write replace
* Export the model results and exit
file write `macro2' "NOTE ***********************************************************************" _n
file write `macro2' "NOTE Export the model results to Stata" _n
file write `macro2' "NOTE ***********************************************************************" _n
if "`fit'" ~= "nofit" {
file write `macro2' "LINK 1 G30" _n
* Put deviance in c1300
file write `macro2' "NAME G30[1] '_Stats'" _n
if ("`discrete'"=="" & "`mcmc'"=="") {
file write `macro2' "LIKE b100" _n
file write `macro2' "EDIT 3 G30[1] b100" _n // Deviance
}
if ("`mcmc'"~="" & "`me'"=="" & `factorson' == 0 & ("`mtype'" ~= "multivariate" | "`discrete'" == "")) { // BDIC doesn't work for models with measurement error or mixed responses
file write `macro2' "BDIC b1 b2 b3 b4" _n
file write `macro2' "EDIT 3 G30[1] b1" _n // DIC
file write `macro2' "EDIT 4 G30[1] b2" _n // pD - effective number of parameters
file write `macro2' "EDIT 5 G30[1] b3" _n // Mean deviance Dbar
file write `macro2' "EDIT 6 G30[1] b4" _n // Deviance at mean parameter values D(thetabar)
}
file write `macro2' "EDIT 7 G30[1] b21" _n // Number of iterations
file write `macro2' "EDIT 8 G30[1] b22" _n // Converged
file write `macro2' "NAME c1098 '_FP_b'" _n
file write `macro2' "NAME c1099 '_FP_v'" _n
file write `macro2' "NAME c1096 '_RP_b'" _n
file write `macro2' "NAME c1097 '_RP_v'" _n
file write `macro2' "NAME c1094 '_esample'" _n
file write `macro2' "SUM '_esample' b1" _n
file write `macro2' "EDIT 9 G30[1] b1" _n // Number of missing rows
if `factorson' == 1 {
file write `macro2' "LINK 4 G29" _n
file write `macro2' "NAME G29[1] '_FACT_b'" _n
file write `macro2' "NAME G29[2] '_FACT_v'" _n
file write `macro2' "DAFL G29[1] G29[2]" _n
file write `macro2' "NAME G29[3] '_FACTVAR_b'" _n
file write `macro2' "NAME G29[4] '_FACTVAR_v'" _n
file write `macro2' "DAFV G29[3] G29[4]" _n
}
if "`savestata'" == "" {
tempfile fileresults
}
else {
local 0 `savestata'
syntax anything(name=savestata), [REPLACE]
if "`replace'" == "" confirm new file "`savestata'"
local fileresults `savestata'
}
if `factorson' == 1 {
file write `macro2' "PSTA '`fileresults'' '_FP_b' '_FP_v' '_RP_b' '_RP_v' '_Stats' '_FACT_b' '_FACT_v' '_FACTVAR_b' '_FACTVAR_v'" _n
file write `macro2' "ERAS '_Stats' '_FACT_b' '_FACT_v' '_FACTVAR_b' '_FACTVAR_v'" _n
file write `macro2' "LINK 0 G29" _n
}
else {
file write `macro2' "PSTA '`fileresults'' '_FP_b' '_FP_v' '_RP_b' '_RP_v' '_Stats'" _n
file write `macro2' "ERAS '_Stats'" _n
}
file write `macro2' "LINK 0 G30" _n
tempfile fileesample
runmlwin_saveesample, fileesample(`fileesample') nummlwineqns(`nummlwineqns') sortindex(`_sortindex') macro2(`macro2')
if "`mcmc'" ~= "" {
file write `macro2' "LINK 0 G21" _n /* Residuals */ "LINK 0 G22" _n /* Factors */ "LINK 0 G23" _n /* Factor variances */
local paramlength = colsof(`matFP_b') + colsof(`matRP_b') // numvars is 4 here
local chainnames `:colfullnames `matFP_b'' `:colfullnames `matRP_b''
if "`mtype'"=="multinomial" & ("`link'"=="ologit"|"`link'"=="oprobit"|"`link'"=="ocloglog") {
local paramlength = `paramlength' - 1
local tmp OD:bcons_1
local chainnames `:list chainnames - tmp'
}
// Export parameter chains
file write `macro2' "NOTE export parameter chain" _n
file write `macro2' "NAME c1091 'deviance'" _n
file write `macro2' "NAME c1090 'mcmcchains'" _n
file write `macro2' "COUNt 'deviance' b40" _n // Thinned iterations
file write `macro2' "CALC b41 = b40 * `thinning'" _n // Unthinned iterations
local i = 1
file write `macro2' "LINK `:list sizeof chainnames' G21" _n
file write `macro2' "FILL G21" _n
foreach var of local chainnames {
if _caller() >= 11 {
file write `macro2' "NAME G21[`i'] '`=strtoname("`var'")''" _n
}
else {
mata: st_local("tmpname", validname("`var'"))
file write `macro2' "NAME G21[`i'] '`=abbrev("`tmpname'", 32)''" _n
}
file write `macro2' "DESC G21[`i'] '`var''" _n
local ++i
}
file write `macro2' "LINK 2 G25" _n
file write `macro2' "NAME G25[1] 'parnum'" _n
file write `macro2' "NAME G25[2] 'iteration'" _n
file write `macro2' "DESC G25[2] '\Iteration'" _n
file write `macro2' "LINK 0 G25" _n
file write `macro2' "CODE `paramlength' 1 b40 'parnum'" _n
file write `macro2' "GENE 1 b41 `thinning' 'iteration'" _n
file write `macro2' "SPLIt 'mcmcchains' 'parnum' G21" _n
if "`mtype'"=="multinomial" & ("`link'"=="ologit"|"`link'"=="oprobit"|"`link'"=="ocloglog") {
file write `macro2' "LINK 1 G25" _n
file write `macro2' "PUT b40 1 G25[1]'" _n
if _caller() >= 11 {
file write `macro2' "NAME G25[1] '`=strtoname("OD:bcons_1")''" _n
}
else {
mata: st_local("tmpname", validname("OD:bcons_1"))
file write `macro2' "NAME G25[1] '`tmpname''" _n
}
file write `macro2' "DESC G25[1] 'OD:bcons_1'" _n
file write `macro2' "GSET 2 G21 G25 G21" _n
file write `macro2' "LINK 0 G25" _n
}
if `factorson' == 1 {
file write `macro2' "LINK 2 G25" _n
file write `macro2' "NAME G25[1] 'factnum'" _n
file write `macro2' "NAME G25[2] 'factvarnum'" _n
file write `macro2' "CODE `=`numfact'*`nummlwineqns'' 1 b40 'factnum'" _n
file write `macro2' "CODE `=(`numfact'*(`numfact'+1))/2' 1 b40 'factvarnum'" _n
file write `macro2' "LINK 0 G25" _n
local count = 0
forvalues l = 1/`maxlevels' {
if "`flinit`l''" ~= "" {
local count = `count' + (`=colsof(`flinit`l'')' * `=rowsof(`flinit`l'')')
}
}
file write `macro2' "LINK `count' G25" _n
local i = 1
forvalues l = 1/`maxlevels' {
if "`flinit`l''" ~= "" {
forvalues fact = 1/`=colsof(`flinit`l'')' {
forvalues resp = 1/`=rowsof(`flinit`l'')' {
file write `macro2' "NAME G25[`i'] 'RP`l'FL_f`fact'_`resp''" _n
file write `macro2' "DESC G25[`i'] 'RP`l'FL:f`fact'_`resp''" _n
local ++i
}
}
}
}
file write `macro2' "SPLIt '_FACT_LOAD_CHAIN' 'factnum' G25" _n
file write `macro2' "GSET 2 G21 G25 G21" _n
file write `macro2' "LINK 0 G25" _n
local count = 0
local basefact = 0
forvalues l = 1/`maxlevels' {
if "`flinit`l''" ~= "" {
forvalues fact1 = 1/`=colsof(`flinit`l'')' {
local count = `count' + `basefact' + `fact1'
}
local basefact = `basefact' + colsof(`flinit`l'')
}
}
file write `macro2' "LINK `count' G25" _n
local i = 1
local basefact = 0
forvalues l = 1/`maxlevels' {
if "`flinit`l''" ~= "" {
forvalues fact1 = 1/`=colsof(`flinit`l'')' {
local i = `i' + `basefact'
forvalues fact2 = 1/`fact1' {
if `fact1' == `fact2' {
file write `macro2' "NAME G25[`i'] 'RP`l'FV_var_f`fact1'_'" _n
file write `macro2' "DESC G25[`i'] 'RP`l'FV:var(f`fact1')'" _n
}
else {
file write `macro2' "NAME G25[`i'] 'RP`l'FV_cov_f`fact1'_f`fact2'_'" _n
file write `macro2' "DESC G25[`i'] 'RP`l'FV:cov(f`fact1',f`fact2')'" _n
}
local ++i
}
}
local basefact = `basefact' + colsof(`flinit`l'')
}
}
file write `macro2' "SPLIt '_FACT_VAR_CHAIN' 'factvarnum' G25" _n
file write `macro2' "GSET 2 G21 G25 G21" _n
file write `macro2' "LINK 0 G25" _n
file write `macro2' "ERAS 'factnum' 'factvarnum'" _n
}
if ("`savechains'"=="") {
tempfile savechains
}
file write `macro2' "PSTA '`savechains'' 'iteration' 'deviance' G21" _n
file write `macro2' "ERAS 'parnum' 'iteration' G21" _n
file write `macro2' "LINK 0 G21" _n
// Export residual chains
forvalues l = 1/`maxlevels' {
if "`saveresiduals`l''" ~= "" {
file write `macro2' "NOTE export residual chain for level `l'" _n
local ll = `l'
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") local ++ll
file write `macro2' "NOBS `ll' b31 b32" _n
file write `macro2' "CALC b33 = b31 * `numresi`l''" _n
file write `macro2' "CALC b34 = b31 * b40" _n
file write `macro2' "LINK 4 G25" _n
// Generate ID column the same length as the residuals
if ("`cc'"=="") {
local hier
forvalues lev = `l'/`maxlevels' {
local hier '`lev`lev'id'' `hier'
}
file write `macro2' "COMB `hier' G25[4]" _n
file write `macro2' "TAKE G25[4] '`lev`l'id'' G25[3] G25[2]" _n
file write `macro2' "ERAS G25[3] G25[4]" _n
}
else if "`mmids`l''" ~= "" {
file write `macro2' "ERAS G25[1]" _n
foreach varname of varlist `mmids`l'' {
file write `macro2' "APPE G25[1] '`varname'' G25[1]" _n
}
file write `macro2' "UNIQ G25[1] G25[2]" _n
file write `macro2' "OMIT 0 G25[2] G25[2]" _n
}
else file write `macro2' "UNIQ '`lev`l'id'' G25[2]" _n
file write `macro2' "ERAS G25[1]" _n
file write `macro2' "REPE `numresi`l'' G25[2] G25[3]" _n
file write `macro2' "ERAS G25[2]" _n
file write `macro2' "CODE G25[3] 1 b40 G25[1]" _n
file write `macro2' "ERAS G25[3]" _n
// Temporary rename level ID column so there isn't a duplicate
file write `macro2' "NAME '`lev`l'id'' '_`lev`l'id''" _n
file write `macro2' "NAME G25[1] '`lev`l'id''" _n
file write `macro2' "DESC G25[1] '\\`:variable label `lev`l'id'''" _n
file write `macro2' "CODE b40 b33 1 G25[2]" _n // iteration number
file write `macro2' "CALC G25[2] = ((G25[2] - 1) * `thinning') + 1" _n
file write `macro2' "NAME G25[2] 'iteration'" _n
file write `macro2' "DESC G25[2] '\Iteration'" _n
file write `macro2' "CODE `numresi`l'' 1 b34 G25[3]" _n // residual number
file write `macro2' "CALC G25[3] = G25[3] - 1" _n // Renumber to start at zero
file write `macro2' "NAME G25[3] 'residual'" _n
file write `macro2' "DESC G25[3] '\Residual'" _n
file write `macro2' "CODE b31 `numresi`l'' b40 G25[4]" _n // id number
file write `macro2' "NAME G25[4] 'idnum'" _n
file write `macro2' "DESC G25[4] 'Unit number'" _n
file write `macro2' "NAME 'RESID_CHAIN_LEVEL`l'' 'value'" _n
file write `macro2' "PSTA '`saveresiduals`l''' 'iteration' 'residual' 'idnum' '`lev`l'id'' 'value'" _n
file write `macro2' "ERAS 'iteration' 'residual' 'idnum' '`lev`l'id'' 'value'" _n
file write `macro2' "LINK 0 G25" _n
// Restore level ID column name
file write `macro2' "NAME '_`lev`l'id'' '`lev`l'id''" _n
}
}
}
}
if ("`plugin'"=="") file write `macro2' "EXIT" _n
file close `macro2'
* MLwiN macro (appends macro1 and macro2 together)
tempname macro3
if "`savefullmacro'" ~= "" local filemacro3 `savefullmacro'
else tempfile filemacro3
qui file open `macro3' using "`filemacro3'", write replace
file open `macro1' using "`filemacro1'", read
file open `macro2' using "`filemacro2'", read
file read `macro1' temp
file write `macro3' "`temp'" _n
while r(eof)==0 {
file read `macro1' temp
file write `macro3' "`temp'" _n
}
file read `macro2' temp
file write `macro3' "`temp'" _n
while r(eof)==0 {
file read `macro2' temp
file write `macro3' "`temp'" _n
}
file close `macro1'
file close `macro2'
file close `macro3'
if ("`viewfullmacro'"~="") view "`filemacro3'"
if "`mlwin'"=="nomlwin" {
di as error "The nomlwin option caused runmlwin to end prematurely"
exit 198
}
* Call either MLwiN.exe or MLwiN.plugin
if "`plugin'"~="" quietly mlncommand OBEY '`filemacro3''
else {
if "`mlwinpath'"=="" {
di as error "You must specify the file address for MLwiN.exe using either:" _n(2)
di as error "(1) the mlwinpath() option; for example mlwinpath(C:\Program Files\MLwiN v3.07\mlwin.exe)" _n
di as error "(2) a global called MLwiN_path; for example: . global MLwiN_path C:\Program Files\MLwiN v3.07\mlwin.exe" _n(2)
di as error "We recommend (2) and that the user places this global macro command in profile.do, see help profile." _n(2)
di as error "IMPORTANT: Make sure that you are using the latest version of MLwiN. This is available at: https://www.bristol.ac.uk/cmm/MLwiN/index.shtml" _n
exit 198
}
if "`batch'" == "" quietly runmlwin_qshell "`mlwinpath'" /run "`filemacro3'"
else {
quietly runmlwin_qshell "`mlwinpath'" /nogui /run "`filemacro3'"
display as text " --- Begin MLwiN error log --- "
type `errlog'
display as text " --- End MLwiN error log --- " _n
}
}
if "`fit'" == "nofit" {
di as txt "The savewinbugs nofit option caused runmlwin to end prematurely."
ereturn scalar mcmcnofit = 1
exit
}
******************************************************************************
* IMPORT MLWIN MODEL RESULTS BACK INTO STATA
******************************************************************************
capture use "`fileresults'", clear
if _rc {
if ("`pause'"=="") display as error "The model did not run properly in MLwiN. You most likely clicked the 'Abort Macro' button in MLwiN, rather than the 'Resume Macro' button." _n
if ("`pause'"=="nopause") display as error "The model did not run properly in MLwiN. Re-run the model without the nopause option to debug the model in MLwiN." _n
exit 198
}
local missrows = _Stats[9]
qui gen _FP_name = ""
qui gen _RP_name = ""
order _FP_name _FP_b _FP_v _RP_name _RP_b _RP_v
* Check the MLwiN columns are the expected length
qui sum _FP_b
capture assert r(N)==colsof(`matFP_b')
if _rc~=0 {
di as error "runmlwin has encountered an error importing the model results from MLwiN. Check that the model has run properly in MLwiN." _n
exit 198
}
qui sum _RP_b
capture assert r(N)==colsof(`matRP_b')
if _rc~=0 {
di as error "runmlwin has encountered an error importing the model results from MLwiN. Check that the model has run properly in MLwiN." _n
exit 198
}
* Put MLwiN values back into FP_b
forvalues c = 1/`=colsof(`matFP_b')' {
mat `matb'[1,`c'] = _FP_b[`c']
}
local startrow = colsof(`matFP_b')
* Put MLwiN values back into RP_b
forvalues c = 1/`=colsof(`matRP_b')' {
mat `matb'[1,`=`startrow'+`c''] = _RP_b[`c']
}
if `factorson' == 1 {
local startrow = `startrow' + colsof(`matRP_b')
forvalues c = 1/`=`numfact'*`nummlwineqns'' {
mat `matb'[1,`=`startrow'+`c''] = _FACT_b[`c']
}
local i = 1
local startrow = `startrow' + `=`numfact'*`nummlwineqns''
local basefact = 0
forvalues l = 1/`maxlevels' {
if "`flinit`l''" ~= "" {
forvalues fact1 = 1/`=colsof(`flinit`l'')' {
forvalues fact2 = 1/`fact1' {
mat `matb'[1,`=`startrow'+`i''] = _FACTVAR_b[`=((`basefact'+`fact1' - 1) * `basefact'+`fact1' / 2) + `basefact'+`fact2'']
local ++i
}
}
local basefact = `basefact' + colsof(`flinit`l'')
}
}
}
* Stack b
// MLwiN has precision problems with constraints so round covariance matrix to 6s.f
local sf = 6
* Put MLwiN values back into FP_V
local i = 1
forvalues r = 1/`=colsof(`matFP_b')' {
forvalues c = 1/`r' {
if ("`constraints'"=="") local val = _FP_v[`i']
else local val = round(_FP_v[`i'], 1/10^(`sf'-trunc(log10(abs(_FP_v[`i'])))))
if `val' == . local val = 0
mat `matV'[`r',`c'] = `val'
mat `matV'[`c',`r'] = `val'
local ++i
}
}
local startrow = colsof(`matFP_b')
* Put MLwiN values back into RP_V
local i = 1
forvalues r = 1/`=colsof(`matRP_b')' {
forvalues c = 1/`r' {
if ("`constraints'"=="") local val = _RP_v[`i']
else local val = round(_RP_v[`i'], 1/10^(`sf'-trunc(log10(abs(_RP_v[`i'])))))
if `val' == . local val = 0
mat `matV'[`=`startrow'+`r'',`=`startrow'+`c''] = `val'
mat `matV'[`=`startrow'+`c'',`=`startrow'+`r''] = `val'
local ++i
}
}
if `factorson' == 1 {
local startrow = `startrow' + colsof(`matRP_b')
forvalues c = 1/`=`numfact'*`nummlwineqns'' {
mat `matV'[`=`startrow'+`c'',`=`startrow'+`c''] = (_FACT_v[`c'])^2
}
local i = 1
local startrow = `startrow' + `=`numfact'*`nummlwineqns''
local basefact = 0
forvalues l = 1/`maxlevels' {
if "`flinit`l''" ~= "" {
forvalues fact1 = 1/`=colsof(`flinit`l'')' {
forvalues fact2 = 1/`fact1' {
mat `matV'[`=`startrow'+`i'',`=`startrow'+`i''] = (_FACTVAR_v[`=((`basefact'+`fact1' - 1) * `basefact'+`fact1' / 2) + `basefact'+`fact2''])^2
local ++i
}
}
local basefact = `basefact' + colsof(`flinit`l'')
}
}
}
* Stack V
* List MLwiN model results
ereturn clear
if ("`verbose'"~="") display "b and V matrices imported from MLwiN:"
if ("`verbose'"~="") mat list `matb'
if ("`verbose'"~="") mat list `matV'
if "`consC'" == "" capture noisily ereturn post `matb' `matV'
else capture noisily ereturn post `matb' `matV' `consC'
if ("`verbose'"~="") display "ereturn error code was: " _rc
if _rc {
display as error "Warning - runmlwin has experienced difficulties importing the standard errors, all standard errors have been set to zero"
tempname zeromat
mat `zeromat' = J(rowsof(`matV'), rowsof(`matV'), 0)
mat rowname `zeromat' = `:rowfullnames `matV''
mat colname `zeromat' = `:colfullnames `matV''
if "`consC'" == "" ereturn post `matb' `zeromat'
else ereturn post `matb' `zeromat' `consC'
}
tempname rptmp
mat `rptmp' = e(b)
forvalues l = `maxlevels'(-1)1 {
if `:list sizeof rp`l'name' > 0 {
mat RP`l' = J(`:list sizeof rp`l'name', `:list sizeof rp`l'name', 0)
local matnames = subinstr("`rp`l'name'", ".", "_", .)
mat rownames RP`l' = `matnames'
mat colnames RP`l' = `matnames'
local r = 1
local c = 1
foreach row of local rp`l'name {
if _caller() < 11 mata: st_local("row", validname("`row'"))
foreach col of local rp`l'name {
if _caller() < 11 {
mata: st_local("col", validname("`col'"))
}
if `c' < `r' {
if _caller() >= 11 {
local pname RP`l':cov(`=abbrev(strtoname("`col'"),13)'\\`=abbrev(strtoname("`row'"),13)')
}
else {
local pname RP`l':cov(`=abbrev("`col'",13)'\\`=abbrev("`row'",13)')
}
if `:list pname in names' {
mat RP`l'[`r', `c'] = `rptmp'[1, "`pname'"]
mat RP`l'[`c', `r'] = `rptmp'[1, "`pname'"]
}
}
if `c'==`r' {
if _caller() >= 11 {
local pname RP`l':var(`=abbrev(strtoname("`row'"), 27)')
}
else {
local pname RP`l':var(`=abbrev("`row'", 27)')
}
if `:list pname in names' mat RP`l'[`r', `c'] = `rptmp'[1, "`pname'"]
}
local ++c
}
local ++r
local c = 1
}
ereturn matrix RP`l' = RP`l'
}
}
ereturn scalar numlevels = `maxlevels'
ereturn scalar k_f = colsof(`matFP_b')
ereturn scalar k_r = colsof(`matRP_b')
ereturn scalar k = colsof(`matFP_b') + colsof(`matRP_b')
ereturn local cmd runmlwin
ereturn local cmdline `e(cmd)' `runmlwin_cmdline'
ereturn local version `mlwinversion'
ereturn scalar maxlevels = `toplevel'
ereturn scalar columns = `numwscolumns'
if ("`cc'"~="") local classtype "(cross-classified)"
else local classtype "(hierarchical)"
if ("`distribution'"~="multinomial") ereturn local depvars `response'
if ("`distribution'"=="multinomial") ereturn local depvars :list uniq response
ereturn local distribution =rtrim(ltrim("`distribution'"))
ereturn local link `link'
ereturn local denominators `denominator'
ereturn local offsets `offset'
ereturn local properties b V
ereturn local method `emode'
ereturn local linearization `linearization'
if "`extra'" ~= "" ereturn scalar extrabinomial = 1
else ereturn scalar extrabinomial = 0
if "`mtype'" == "multinomial" {
ereturn local basecategory `basecategory'
ereturn local respcategories `responsecats'
}
if ("`e(distribution)'"=="normal") ereturn local title Normal response model `classtype'
if ("`e(distribution)'"=="binomial") ereturn local title Binomial `e(link)' response model `classtype'
if ("`e(distribution)'"=="multinomial" & ("`e(link)'"=="mlogit")) ereturn local title Unordered multinomial logit response model `classtype'
if ("`e(distribution)'"=="multinomial" & ("`e(link)'"=="ologit")) ereturn local title Ordered multinomial logit response model `classtype'
if ("`e(distribution)'"=="multinomial" & ("`e(link)'"=="oprobit")) ereturn local title Ordered multinomial probit response model `classtype'
if ("`e(distribution)'"=="multinomial" & ("`e(link)'"=="ocloglog")) ereturn local title Ordered multinomial cloglog response model `classtype'
if ("`e(distribution)'"=="poisson") ereturn local title Poisson response model `classtype'
if ("`e(distribution)'"=="nbinomial") ereturn local title Negative binomial response model `classtype'
if (wordcount("`e(distribution)'")>1) ereturn local title Multivariate response model `classtype'
* Level IDs
ereturn local level1var `lev1id'
local ivars
forvalues l = `maxlevels'(-1)2 {
local ivars `ivars' `lev`l'id'
}
ereturn local ivars `ivars'
* (R)IGLS Weighting options
local weightvar
local weighttype
forvalues l = `maxlevels'(-1)1 {
if "`weight`l''" ~= "" local weightvar `weightvar' `weight`l''
else local weightvar `weightvar' .
if "`weighttype`l''" ~= "" local weighttype `weighttype' `weighttype`l''
else local weighttype `weighttype' standardised
}
ereturn local weightvar `weightvar'
ereturn local weighttype `weighttype'
* Model fit
if ("`discrete'"=="" & "`mcmc'"=="") ereturn scalar ll = -_Stats[3]/2
if ("`discrete'"=="" & "`mcmc'"=="") ereturn scalar deviance = _Stats[3]
ereturn scalar iterations = _Stats[7]
if ("`mcmc'"=="") {
ereturn scalar converged = (_Stats[8]==1)
}
* Save MCMC chains and calculate 95% credible interval and ESS and store as matrices
if ("`mcmc'" ~= "") {
ereturn scalar burnin = `mcmc_burnin'
ereturn scalar chain = `mcmc_chain'
ereturn scalar thinning = `mcmc_thinning'
ereturn scalar dic = _Stats[3] // DIC
ereturn scalar pd = _Stats[4] // pD - effective number of parameters
ereturn scalar dbar = _Stats[5] // Dbar - Mean deviance
ereturn scalar dthetabar = _Stats[6] // D(thetabar) - Deviance at mean parameter values
ereturn scalar converged = 1
ereturn scalar mcmcdiagnostics = ("`diagnostics'"~="nodiagnostics")
if ("`diagnostics'"~="nodiagnostics") {
local chainnames `:colfullnames `matFP_b'' `:colfullnames `matRP_b''
if `factorson' == 1 {
local chainnames `chainnames' `:colfullnames `mat_RP_b_fact'' `:colfullnames `mat_RP_b_factvar''
}
local temp = `:list sizeof chainnames'
runmlwin_savechains, parnames(`chainnames') chainresults(`savechains')
local savechainnames
foreach parname of local chainnames {
if _caller() >= 11 {
local savechainnames `savechainnames' `=strtoname("`parname'")'
}
else {
mata: st_local("tmpname", validname("`parname'"))
local savechainnames `savechainnames' `=abbrev("`tmpname'", 32)'
}
}
use `savechainnames' using "`savechains'", clear
foreach stat in meanmcse sd mode ess lb ub rl1 rl2 bd pvalmean pvalmedian pvalmode {
matrix `stat' = J(1, `temp', .)
matrix colnames `stat' = `chainnames'
}
local acfpoints = 100
local pacfpoints = 10
local mcsepoints = 1000
local kdpoints = 1000
matrix ACF = J(`acfpoints', `=1 + `temp'', .)
matrix colnames ACF = point `chainnames'
tempname tmpacf
matrix PACF = J(`pacfpoints', `=1 + `temp'', .)
matrix colnames PACF = point `chainnames'
tempname tmppacf
matrix quantiles = J(9, `=1 + `temp'', .)
matrix colnames quantiles = quant `chainnames'
tempname tmpquantiles
tempname mataMCSE
mata: `mataMCSE' = J(`mcsepoints',`=1 + `temp'',.)
tempname mataKD1
mata: `mataKD1' = J(`kdpoints',`temp',.)
tempname mataKD2
mata: `mataKD2' = J(`kdpoints',`temp',.)
local i = 1
foreach par of local savechainnames {
// set posit(1) for the variance parameters
local posit = 0
if regexm("`par'", "RP[0-9]+_var_([a-zA-Z0-9_]+)_") == 1 local posit = 1
quietly runmlwin_mcmcdiag `par', level(`level') thinning(`thinning') posit(`posit') acfpoints(`acfpoints') pacfpoints(`pacfpoints') mcsepoints(`mcsepoints') kernalpoints(`kdpoints') `eform' // Need to add prior option
matrix meanmcse[1,`i'] = r(meanmcse)
matrix sd[1,`i'] = r(sd)
matrix mode[1,`i'] = r(mode)
matrix ess[1,`i'] = r(ESS)
matrix lb[1,`i'] = r(lb)
matrix ub[1,`i'] = r(ub)
matrix rl1[1,`i'] = r(RL1)
matrix rl2[1,`i'] = r(RL2)
matrix bd[1,`i'] = r(BD)
matrix pvalmean[1,`i'] = r(pvalmean)
matrix pvalmedian[1,`i'] = r(pvalmedian)
matrix pvalmode[1,`i'] = r(pvalmode)
matrix `tmpacf' = r(ACF)
if `i' == 1 matrix ACF[1, 1] = `tmpacf'[1..., 2]
matrix ACF[1, `=1 + `i''] = `tmpacf'[1..., 1]
matrix `tmppacf' = r(PACF)
if `i' == 1 matrix PACF[1, 1] = `tmppacf'[1..., 2]
matrix PACF[1, `=1 + `i''] = `tmppacf'[1..., 1]
matrix `tmpquantiles' = r(quantiles)
if `i' == 1 matrix quantiles[1, 1] = `tmpquantiles'[1..., 1]
matrix quantiles[1, `=1 + `i''] = `tmpquantiles'[1..., 2]
if `i' == 1 mata: `mataMCSE'[1..., 1] = st_matrix("r(MCSE)")[1..., 2]
mata: `mataMCSE'[1..., `=1 + `i''] = st_matrix("r(MCSE)")[1..., 1]
mata: `mataKD1'[1..., `i'] = st_matrix("r(KD)")[1..., 1]
mata: `mataKD2'[1..., `i'] = st_matrix("r(KD)")[1..., 2]
local ++i
}
ereturn scalar level = `level'
mata: st_matrix("MCSE", `mataMCSE')
mata: mata drop `mataMCSE'
matrix colnames MCSE = point `chainnames'
mata: st_matrix("KD1", `mataKD1')
mata: mata drop `mataKD1'
matrix colnames KD1 = `chainnames'
mata: st_matrix("KD2", `mataKD2')
mata: mata drop `mataKD2'
matrix colnames KD2 = `chainnames'
ereturn local chains "runmlwin_chains"
foreach stat in meanmcse sd mode ess lb ub rl1 rl2 bd pvalmean pvalmedian pvalmode {
ereturn matrix `stat' = `stat'
}
ereturn matrix ACF = ACF
ereturn matrix PACF = PACF
ereturn matrix quantiles = quantiles
ereturn matrix MCSE = MCSE
ereturn matrix KD1 = KD1
ereturn matrix KD2 = KD2
}
}
* Multiple Imputation - remove _mi_ prefix from response
if `:list sizeof imputefiles' ~= 0 {
foreach impfile of local imputefiles {
use "`impfile'", clear
if _caller() >= 12 rename _mi_* *
else renpfix "_mi_"
quietly saveold "`impfile'", replace
}
}
* Restore the original data set
restore
* Merge in esample
if `missrows' > 0 {
if _caller() >= 11 quietly merge m:1 `_sortindex' using "`fileesample'", assert(master match) nogenerate
else {
quietly merge `_sortindex' using "`fileesample'", uniqusing sort
assert inlist(_merge,1,3)
drop _merge
}
qui recode _esample (.=0)
}
else qui gen _esample = `touse'
ereturn repost, esample(_esample)
* Calculate the sample size
qui count if e(sample)
ereturn scalar N = r(N)
* Calculate the min, mean and max number of units at each level
if `maxlevels' > 1 & "`sort'" == "" {
matrix N_g = J(1, `=`maxlevels'-1', .)
matrix g_min = J(1, `=`maxlevels'-1', .)
matrix g_avg = J(1, `=`maxlevels'-1', .)
matrix g_max = J(1, `=`maxlevels'-1', .)
local hier
forvalues l = `maxlevels'(-1)2 {
if ("`cc'"=="") local hier `hier' `lev`l'id'
if ("`cc'"~="") local hier `lev`l'id'
tempvar esample temp1 temp2
quietly gen `esample' = (e(sample)==1)
qui bysort `hier' `esample': gen byte `temp1' = 1 if _n==1 & `esample'==1
qui bysort `hier' `esample': gen `temp2' = _N if `esample'==1
qui sum `temp2' if `temp1'==1 & `esample'==1
drop `esample' `temp1' `temp2'
matrix N_g[1, `=`maxlevels' -`l'+1'] = r(N)
matrix g_min[1, `=`maxlevels' -`l'+1'] = r(min)
matrix g_avg[1, `=`maxlevels' -`l'+1'] = r(mean)
matrix g_max[1, `=`maxlevels' -`l'+1'] = r(max)
}
sort `_sortindex'
ereturn matrix N_g = N_g
ereturn matrix g_min = g_min
ereturn matrix g_avg = g_avg
ereturn matrix g_max = g_max
}
* Merge residuals at each level into the original data set
forvalues l = `maxlevels'(-1)1 {
local hier
forvalues t = `maxlevels'(-1)`l' {
local hier `hier' `lev`t'id'
}
local ll = `l'
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") local ++ll
if ("`residuals`l''"~="") runmlwin_mergeresiduals , numlevels(`ll') touse(`touse') hier(`hier') cc(`cc') levelid(`lev`l'id') residfile(`residfile`l'') sort(`sort') mmid(`mmids`l'')
if ("`fscores`l''"~="") runmlwin_mergefscores , numlevels(`ll') touse(`touse') hier(`hier') cc(`cc') levelid(`lev`l'id') fscoresfile(`fscoresfile`l'') sort(`sort')
}
* Merge simulated response into the original data set
if ("`simulate'"~="") runmlwin_mergesimresp, simfile(`simfile')
if "`imputesummaries'" ~= "" {
if _caller() >= 11 qui merge m:1 `_sortindex' using "`imputesummaries'", assert(match) nogenerate
else {
quietly merge `_sortindex' using "`imputesummaries'", uniqusing sort
assert inlist(_merge,3)
drop _merge
}
}
* Append multiple imputed data sets
if `:list sizeof imputefiles' ~= 0 {
if _caller() >= 11 {
mi import flongsep tmpimpfile, using(`imputefiles') id(`_sortindex') imputed(`response') clear
mi convert mlong
quietly mi erase tmpimpfile
}
else {
tempvar tag
gen `tag' = 1 if missing(`=subinstr("`response'"," ",",",.)')
foreach resp of local response {
rename `resp' `resp'_mi_m0
}
local i = 1
foreach imputedfile of local imputefiles {
merge `_sortindex' using `imputedfile', unique sort
drop _merge
foreach resp of local response {
rename `resp' `resp'_mi_m`i'
}
local ++i
}
quietly reshape long `=subinstr("`response' "," ","_mi_m ",.)', i(`_sortindex') j(_mi_m)
foreach resp of local response {
rename `resp'_mi_m `resp'
}
quietly drop if `tag'~=1 & _mi_m>0
drop `tag'
gen _mi_id = `_sortindex'
quietly gen _mi_miss = 1 if missing(`=subinstr("`response'"," ",",",.)')
unab vars : *
local mivar _mi_m
local vars `:list vars - mivar'
local vars `vars' _mi_m
order `vars'
}
}
* Sort the data by original sort order
if `: list sizeof imputefiles' ~= 0 { // Because data set has extra rows!
sort `_sortindex' _mi_m
quietly replace `_sortindex' = _n // :sortorder should still be okay, as the data is sorted the same as before, but expanded upwards by the imputed data sets
}
end
******************************************************************************
* CREATES PARAMETER NAMES FOR RP_B (I.E. SIGMA2_U0, SIGMA_U01, SIGMA2_U1, SIGMA2_E)
******************************************************************************
// Note we can't use "," in for example c(cons,standlrt). Doing so causes problems with estimates table when you try to display more than one model
// The "(" and ")" will prob also be problematic
// What we need is to display c(cons,standlrt) in the output window but to store the parameter with different notation without the ",", "(" and ")"
// A separate problem is that instead of c(cons,standlrt) we have c(standlrt,cons)
capture program drop runmlwin_rplpars
program define runmlwin_rplpars, rclass
if _caller() >= 12 version 12.0
if _caller() <= 9 version 9.0
syntax , RPName(string)
local rplong
local r = 1
local c = 1
foreach row of local rpname {
if _caller() < 11 {
mata: st_local("row", validname("`row'"))
}
foreach col of local rpname {
if _caller() < 11 {
mata: st_local("col", validname("`col'"))
}
if `c' < `r' {
if _caller() >= 11 {
local rplong `rplong' cov(`=abbrev(strtoname("`col'"),13)'\\`=abbrev(strtoname("`row'"),13)') // Note that `col' and `row' are deliberately the wrong way around in order to follow MLwiN convention
}
else {
local rplong `rplong' cov(`=abbrev("`col'",13)'\\`=abbrev("`row'",13)') // Note that `col' and `row' are deliberately the wrong way around in order to follow MLwiN convention
}
local ++c
}
if `c'==`r' {
if _caller() >= 11 {
local rplong `rplong' var(`=abbrev(strtoname("`row'"), 27)')
}
else {
local rplong `rplong' var(`=abbrev("`row'", 27)')
}
local ++c
}
}
local ++r
local c = 1
}
return local rplong `rplong'
end
capture program drop runmlwin_saveesample
program define runmlwin_saveesample
if _caller() >= 12 version 12.0
if _caller() <= 9 version 9.0
syntax, FILEESAMPLE(string) NUMMLWINEQNS(string) SORTINDEX(string) MACRO2(string)
local _sortindex `sortindex'
tempname fh
file open `fh' using "`fileesample'", write
file close `fh'
file write `macro2' "NOTE generate esample for Stata if there a missing values" _n
file write `macro2' "SWIT b1" _n
file write `macro2' "CASE 0:" _n
file write `macro2' "LEAVE" _n
file write `macro2' "CASE:" _n
file write `macro2' "CALC '_esample' = abso('_esample' - 1)" _n // Because in MLwiN 1 = not in estimation sample
if `nummlwineqns'>1 {
file write `macro2' "NOTE unexpand esample for multinomial/multivariate" _n
file write `macro2' "LINK 1 G21" _n
file write `macro2' "COUNT '`_sortindex'' b1" _n
file write `macro2' "CODE b1 `nummlwineqns' 1 G21[1]" _n
file write `macro2' "MLAV G21[1] '_esample' '_esample'" _n
file write `macro2' "TAKE G21[1] '_esample' G21[1] '_esample'" _n
file write `macro2' "CALC '_esample' = '_esample' > 0" _n
file write `macro2' "ERAS G21[1]" _n
file write `macro2' "LINK 0 G21" _n
}
file write `macro2' "PSTA '`fileesample'' '`_sortindex'' '_esample'" _n
file write `macro2' "ENDS" _n
end
capture program drop runmlwin_simresp
program define runmlwin_simresp
if _caller() >= 12 version 12.0
if _caller() <= 9 version 9.0
syntax, SIMFILE(string) SIMULATE(string) SIMFILE(string) MACRO1(string)
tempname fh
file open `fh' using "`simfile'", write
file close `fh'
file write `macro1' "NOTE Simulate new response variable" _n
file write `macro1' "LINK 2 G21" _n
file write `macro1' "SIMU G21[1]" _n
file write `macro1' "NAME G21[1] '`simulate''" _n
file write `macro1' "COUN G21[1] b30" _n
file write `macro1' "GENE 1 b30 1 G21[2]" _n
file write `macro1' "NAME G21[2] '_simprespid'" _n
file write `macro1' "PSTA '`simfile'' '`simulate'' '_simprespid'" _n
file write `macro1' "ERAS '`simulate'' '_simprespid'" _n
file write `macro1' "LINK 0 G21" _n
end
capture program drop runmlwin_setinit
program define runmlwin_setinit
if _caller() >= 12 version 12.0
if _caller() <= 9 version 9.0
syntax , FPB(string) FPV(string) RPB(string) RPV(string) INITB(string) INITV(string) MTYPE(string) LINK(string) [MCMC(string)] MACRO1(string)
local matFP_b `fpb'
local matFP_V `fpv'
local matRP_b `rpb'
local matRP_V `rpv'
local matb `initb'
local matV `initv'
tempname matINIT
* Fixed part initial values
file write `macro1' "NOTE Specify fixed part initial values" _n
local names : colfullnames `matFP_b'
local i = 1
foreach p of local names {
file write `macro1' "NOTE `p'" _n
mat `matINIT' = `matb'[1,"`p'"]
file write `macro1' "EDIT `i' c1098 `=`matINIT'[1,1]'" _n
local ++i
}
file write `macro1' _n
* Fixed part initial values sampling (co)variances
file write `macro1' "NOTE Specify fixed part initial sampling (co)variances values" _n
local names : colfullnames `matFP_V'
local i = 1
local e = 1
file write `macro1' "PUT `=(`=rowsof(`matFP_V')'*(`=rowsof(`matFP_V')' + 1))/2' 0 c1099" _n
foreach p1 of local names {
local j = 1
foreach p2 of local names {
if `i' >= `j' {
file write `macro1' "NOTE `p1', `p2'" _n
mat `matINIT' = `matV'["`p1'","`p2'"]
file write `macro1' "EDIT `e' c1099 `=`matINIT'[1,1]'" _n
local ++e
}
local ++j
}
local ++i
}
file write `macro1' _n
* Random part initial values
// Need to check that the correct number of initial values have been specified
// what if no random part e.g. single level logit
file write `macro1' "NOTE Specify random part initial values" _n
local names : colfullnames `matRP_b'
local P = colsof(`matRP_b')
local i = 1
foreach p of local names {
if `i'<=`P' {
file write `macro1' "NOTE `p'" _n
mat `matINIT' = `matb'[1,"`p'"]
file write `macro1' "EDIT `i' c1096 `=`matINIT'[1,1]'" _n
local ++i
}
}
file write `macro1' _n
* Random part initial values standard errors
// Need to check that the correct number of initial values have been specified
// what if no random part e.g. single level logit
file write `macro1' "NOTE Specify random part initial sampling (co)variances values" _n
local names : colfullnames `matRP_V'
local P = colsof(`matRP_V')
if ("`mtype'"=="multinomial" & "`link'"=="mlogit") local P = `P' - 1 // Note we might have to add unordered to the if statement as I think ordered models do not have this problem
local i = 1
local e = 1
file write `macro1' "PUT `=(`=rowsof(`matRP_V')'*(`=rowsof(`matRP_V')' + 1))/2' 0 c1097" _n
foreach p1 of local names {
local j = 1
foreach p2 of local names {
if `i' >= `j' & `i'<=`P' & `j' <= `P' {
file write `macro1' "NOTE `p1', `p2'" _n
mat `matINIT' = `matV'["`p1'","`p2'"]
file write `macro1' "EDIT `e' c1097 `=`matINIT'[1,1]'" _n
local ++e
}
local ++j
}
local ++i
}
if "`mtype'"=="multinomial" & "`link'"=="mlogit" & "`mcmc'" ~= "" file write `macro1' "EDIT `e' c1097 0" _n
end
capture program drop runmlwin_writepriors
program define runmlwin_writepriors
if _caller() >= 12 version 12.0
if _caller() <= 9 version 9.0
syntax , MAXLevels(integer) PRIORMAT(string) FPB(string) RPB(string) MACRO1(string) PRIORCOL(string)
local priormatrix `priormat'
local matFP_b `fpb'
local matRP_b `rpb'
tempname matprior
tempname matpriorSD
tempname mateq
if rowsof(`priormatrix')~=2 {
display "The prior matrix must have two rows. The 1st row is for the means. The 2nd row is for the SDs."
exit 198
}
local r = 1
local names :colfullnames `matFP_b'
foreach var of local names {
mat `matprior' = `priormatrix'[1, "`var'"]
mat `matpriorSD' = `priormatrix'[2, "`var'"]
local prior = `matprior'[1,1]
local priorsd = `matpriorSD'[1,1]
file write `macro1' "NOTE fixed part prior for `var'" _n
if `prior' ~= . {
file write `macro1' "EDIT `r' `priorcol' 1" _n
local ++r
file write `macro1' "EDIT `r' `priorcol' `prior'" _n
local ++r
file write `macro1' "EDIT `r' `priorcol' `priorsd'" _n
local ++r
}
else {
file write `macro1' "EDIT `r' `priorcol' 0" _n
local ++r
}
}
forvalues l = `maxlevels'(-1)1 {
local hasprior = 0
local eqnames : coleq `matRP_b'
local levname "RP`l'"
local match : list levname in eqnames
if `l' == 1 & `match' == 0 {
local levname "OD"
local match : list levname in eqnames
}
if `match' == 1 {
matrix `mateq' = `matRP_b'[1, "`levname':"]
local names :colfullnames `mateq'
local first = 1
foreach var of local names {
mat `matprior' = `priormatrix'[1, "`var'"]
mat `matpriorSD' = `priormatrix'[2, "`var'"]
if (`matprior'[1,1] == . & `matpriorSD'[1,1] ~= .) | (`matprior'[1,1] ~= . & `matpriorSD'[1,1] == .) {
display as error "Either both mean and SD have to be specified, nor neither"
exit 198
}
if `first' ~= 1 {
if (`hasprior' == 1 & `matprior'[1,1] == .) | (`hasprior' == 0 & `matprior'[1,1] ~= .) {
display as error "Whole matrix must be specified for random part priors"
exit 198
}
if `priorsd' ~= `matpriorSD'[1,1] {
display as error "Inconsistent sample size"
exit 198
}
}
else {
local priorsd = `matpriorSD'[1,1]
}
if `matprior'[1,1] ~= . local hasprior = 1
local first = 0
}
file write `macro1' "NOTE random part prior for level `l'" _n
if `hasprior' == 1 {
file write `macro1' "EDIT `r' `priorcol' 1" _n
local ++r
foreach var of local names {
mat `matprior' = `priormatrix'[1, "`var'"]
local prior = `matprior'[1,1]
file write `macro1' "EDIT `r' `priorcol' `prior'" _n
local ++r
}
// Check this is a positive integer
if floor(`priorsd') ~= `priorsd' | `priorsd' < 0 {
display as error "Sample size must be a positive integer"
exit 198
}
file write `macro1' "EDIT `r' `priorcol' `priorsd'" _n
local ++r
}
else {
file write `macro1' "EDIT `r' `priorcol' 0" _n
local ++r
}
}
}
file write `macro1' "PRIO `priorcol'" _n
end
******************************************************************************
* CALCULATE RESIDUALS IN MLWIN
******************************************************************************
capture program drop runmlwin_calcresiduals
program define runmlwin_calcresiduals
if _caller() >= 12 version 12.0
if _caller() <= 9 version 9.0
syntax , LEVEL(integer) LEVXID(string) NUMRPXVARS(string) RPXVARS(string) RESIDUALSX(string) MTYPE(string) [STANDardised LEVerage INFluence DELetion SAMPling REFlated VARiances noRECode RPAdjust noAdjust MCMC(string) CC(string) MMID(string) *] MACRO1(string) RESIDFile(string) MAXLEVELS(integer)
local lids
forvalues lev = `maxlevels'(-1)1 {
local lids `lids' lev`lev'id(string)
}
local 0 ,`options'
syntax, `lids'
if "`numrpxvars'" == "0" local numrpxvars "1"
file write `macro1' "NOTE Store level `level' residuals" _n
if "`recode'" ~= "" {
file write `macro1' "MISR 0" _n
}
file write `macro1' "LINK 0 G21" _n // Residual Estimates
file write `macro1' "LINK 0 G22" _n // Residual Standard Errors/Variance
file write `macro1' "LINK 0 G23" _n // Standardised Residual Estimates
file write `macro1' "LINK 0 G24" _n // Leverage Residuals
file write `macro1' "LINK 0 G25" _n // Reflated Residuals
file write `macro1' "LINK 0 G26" _n // Deletion Residuals
file write `macro1' "LINK 0 G27" _n // Influence Residuals
file write `macro1' "LINK 0 G28" _n // Sampling Residuals
file write `macro1' "LINK 0 G29" _n // Temp
file write `macro1' "LINK 0 G30" _n // Temp
file write `macro1' "LINK `numrpxvars' G21" _n
file write `macro1' "FILL G21" _n
file write `macro1' "LINK `numrpxvars' G22" _n
file write `macro1' "FILL G22" _n
forvalues i = 1/`numrpxvars' {
* Name the columns of residual estimates
file write `macro1' "NAME G21[`i'] '`residualsx'`=`i' - 1''" _n
file write `macro1' "DESC G21[`i'] '`residualsx'`=`i' - 1' residual estimate'" _n
}
if "`variances'" == "" {
forvalues i = 1/`numrpxvars' {
* Name the columns of residual standard errors
file write `macro1' "NAME G22[`i'] '`residualsx'`=`i' - 1'se'" _n
file write `macro1' "DESC G22[`i'] '`residualsx'`=`i' - 1'se residual standard error'" _n
}
}
else {
local residual_var
forvalues i = 1/`numrpxvars' {
* Name the columns of residual standard errors
file write `macro1' "NAME G22[`i'] '`residualsx'`=`i' - 1'var'" _n
file write `macro1' "DESC G22[`i'] '`residualsx'`=`i' - 1'var residual variance'" _n
}
}
file write `macro1' "RFUN" _n
if "`variances'" == "" file write `macro1' "ROUT G21 G22" _n // Designate column c300 for residuals and column C301 for their associated variances
else file write `macro1' "ROUT G21 G22" _n // Designate column c300 for residuals and column C301 for their associated variances
local ll = `level'
if ("`mtype'"=="multivariate" | "`mtype'"=="multinomial") local ++ll
file write `macro1' "RLEV `ll'" _n // Specify the level at which variances are to be calculated
file write `macro1' "RCOV 1" _n // Output residuals and their variances
if "`standardised'" ~= "" | "`deletion'" ~= "" | "`leverage'" ~= ""{ // deletion requires standardised to be calculated, leverage requires standard error calculated here
file write `macro1' "LINK `numrpxvars' G23" _n
file write `macro1' "FILL G23" _n
forvalues i = 1/`numrpxvars' {
* Name the columns of standardised residual estimates
file write `macro1' "NAME G23[`i'] '`residualsx'`=`i' - 1'std'" _n
file write `macro1' "DESC G23[`i'] '`residualsx'`=`i' - 1'std standardised residual'" _n
}
file write `macro1' "RTYP 0" _n // Compute diagnostic variances
if "`mcmc'" ~= "" {
file write `macro1' "MCRE" _n // Save the residuals and the variances
}
else {
file write `macro1' "RESI" _n // Save the residuals and the variances
}
local ccount = 1
forvalues i = 1/`numrpxvars' {
file write `macro1' "CALC G23[`i'] = G21[`i']/sqrt(G22[`i'])" _n // Convert the variances to standard errors
}
}
// NOTE leverage depends on residuals with rtype 0
if "`leverage'" ~= "" | "`influence'" ~= "" { // influence requires leverage to be calculated
file write `macro1' "LINK `numrpxvars' G24" _n
file write `macro1' "FILL G24" _n
forvalues i = 1/`numrpxvars' {
* Name the columns of leverage residual estimates
file write `macro1' "NAME G24[`i'] '`residualsx'`=`i' - 1'lev'" _n
file write `macro1' "DESC G24[`i'] '`residualsx'`=`i' - 1'lev leverage residual'" _n
}
forvalues i = 1/`numrpxvars' {
file write `macro1' "LINK 1 G30" _n
file write `macro1' "OMEGa `ll' '`=word("`rpxvars'", `i')'' G30[1]" _n // retrieve variance for corresponding random parameter
file write `macro1' "PICK 1 G30[1] b50" _n
file write `macro1' "ERASe G30[1]" _n
file write `macro1' "LINK 0 G30" _n
file write `macro1' "CALC G24[`i'] = 1 - sqrt(G22[`i']) / sqrt(b50)" _n
}
if "`standardised'" == "" & "`deletion'" == "" {
file write `macro1' "ERASE G23" _n
file write `macro1' "LINK 0 G23" _n
}
}
if "`adjust'" == "noadjust" {
file write `macro1' "RTYP 4" _n // No adjustment
}
else if "`rpadjust'" ~= "" {
file write `macro1' "RTYP 2" _n // Fixed and Random part adjustment
}
else {
file write `macro1' "RTYP 1" _n // Fixed part adjustment
}
if "`mcmc'" ~= "" {
file write `macro1' "MCRE" _n // Save the residuals and the variances
}
else {
file write `macro1' "RESI" _n // Save the residuals and the variances
}
if "`reflated'" ~= "" {
file write `macro1' "LINK `numrpxvars' G25" _n
file write `macro1' "FILL G25" _n
forvalues i = 1/`numrpxvars' {
* Name the columns of reflated residual estimates
file write `macro1' "NAME G25[`i'] '`residualsx'`=`i' - 1'ref'" _n
file write `macro1' "DESC G25[`i'] '`residualsx'`=`i' - 1' reflated residual estimate'" _n
}
}
if "`reflated'" ~= "" {
file write `macro1' "REFLate G21 G25" _n
}
if "`variances'" == "" {
forvalues i = 1/`numrpxvars' {
file write `macro1' "CALC G22[`i'] = sqrt(G22[`i'])" _n // Convert the variances to standard errors
}
}
if "`deletion'" ~= "" | "`influence'" ~= "" { // influence requires deletion to be calculated
forvalues i = 1/`numrpxvars' {
file write `macro1' "LINK `numrpxvars' G26" _n
file write `macro1' "FILL G26" _n
* Name the columns of deletion residual estimates
file write `macro1' "NAME G26[`i'] '`residualsx'`=`i' - 1'del'" _n
file write `macro1' "DESC G26[`i'] '`residualsx'`=`i' - 1'del deletion residual'" _n
}
file write `macro1' "NOBS `ll' b31 b32" _n
forvalues i = 1/`numrpxvars' {
file write `macro1' "CALC G26[`i'] = G23[`i'] / sqrt((b31 - 1 - G23[`i'] ^2)/(b31 - 2))" _n
}
if "`standardised'" == "" & "`leverage'" == ""{
file write `macro1' "ERASE G23" _n
file write `macro1' "LINK 0 G23" _n
}
}
if "`influence'" ~= "" {
forvalues i = 1/`numrpxvars' {
file write `macro1' "LINK `numrpxvars' G27" _n
file write `macro1' "FILL G27" _n
* Name the columns of influence residual estimates
file write `macro1' "NAME G27[`i'] '`residualsx'`=`i' - 1'inf'" _n
file write `macro1' "DESC G27[`i'] '`residualsx'`=`i' - 1'inf influence residual'" _n
}
forvalues i = 1/`numrpxvars' {
file write `macro1' "SUM G24[`i'] b50" _n
file write `macro1' "CALC G27[`i'] = G24[`i'] / b50" _n
file write `macro1' "CALC G27[`i'] = sqrt(G27[`i'] / (1 - G27[`i'])) * abso(G26[`i'])" _n
}
if "`deletion'" == "" {
file write `macro1' "ERASE G26" _n
file write `macro1' "LINK 0 G26" _n
}
if "`leverage'" == "" {
file write `macro1' "ERASE G24" _n
file write `macro1' "LINK 0 G24" _n
}
}
if "`sampling'" ~= "" {
file write `macro1' "LINK `=(`numrpxvars'*(`numrpxvars'+1))/2' G28" _n
local numcombs = 0
forvalues i = 1/`numrpxvars' {
forvalues j = 1/`i' {
local ++numcombs
if `i' == `j' {
file write `macro1' "NAME G28[`numcombs'] '`residualsx'`=`j'-1'var'" _n
file write `macro1' "DESC G28[`numcombs'] '`residualsx'`=`j'-1'var sampling variance'" _n
}
else {
file write `macro1' "NAME G28[`numcombs'] '`residualsx'`=`j'-1'`residualsx'`=`i'-1'cov'" _n
file write `macro1' "DESC G28[`numcombs'] '`residualsx'`=`j'-1'`residualsx'`=`i'-1'cov sampling covariance'" _n
}
}
}
file write `macro1' "LINK `numrpxvars' G29" _n
file write `macro1' "FILL G29" _n
file write `macro1' "LINK 2 G30" _n
file write `macro1' "FILL G30" _n
file write `macro1' "RFUN " _n
file write `macro1' "ROUT G29 G30[2]" _n
file write `macro1' "RCOV 2" _n
file write `macro1' "RESI" _n
// NOTE: This is square rooted, as the residual covariances are sometimes negative
file write `macro1' "NOBS `ll' b31 b32" _n
file write `macro1' "CODE `numcombs' 1 b31 G30[1]" _n
file write `macro1' "SPLIt G30[2] G30[1] G28" _n
file write `macro1' "ERAS G30" _n
file write `macro1' "LINK 0 G30" _n
file write `macro1' "ERAS G29" _n
file write `macro1' "LINK 0 G29" _n
}
file write `macro1' _n
file write `macro1' "LINK 1 G30" _n
file write `macro1' "NOBS `ll' b30 b31" _n
file write `macro1' "GENE 1 b30 1 G30[1]" _n
file write `macro1' "NAME G30[1] '_residualid'" _n
file write `macro1' "LINK 4 G29" _n
// Generate ID column the same length as the residuals
if ("`cc'"=="") {
local hier
forvalues lev = `level'/`maxlevels' {
local hier '`lev`lev'id'' `hier'
}
file write `macro1' "COMB `hier' G29[4]" _n
file write `macro1' "TAKE G29[4] '`lev`level'id'' G29[3] G29[2]" _n
file write `macro1' "ERAS G29[3] G29[4]" _n
}
else if "`mmid'" ~= "" {
file write `macro1' "ERAS G29[1]" _n
foreach varname of varlist `mmid' {
file write `macro1' "APPE G29[1] '`varname'' G29[1]" _n
}
file write `macro1' "UNIQ G29[1] G29[2]" _n
file write `macro1' "OMIT 0 G29[2] G29[2]" _n
}
else file write `macro1' "UNIQ '`lev`level'id'' G29[2]" _n
file write `macro1' "OMIT missing G29[2] G29[2]" _n
// Temporary rename level ID column so there isn't a duplicate
file write `macro1' "NAME '`lev`level'id'' '_`lev`level'id''" _n
file write `macro1' "NAME G29[2] '`lev`level'id''" _n
file write `macro1' "DESC G29[2] '\\`:variable label `lev`level'id'''" _n
if ("`cc'"~="") { // Ensure the data is sorted by ID variable, as the residual ID is based on it this should already be the case
file write `macro1' "SORT 1 '`lev`level'id'' '_residualid' G21 G22 G23 G24 G25 G26 G27 G28 '`lev`level'id'' '_residualid' G21 G22 G23 G24 G25 G26 G27 G28" _n
}
file write `macro1' "PSTA '`residfile'' '_residualid' '`lev`level'id'' G21 G22 G23 G24 G25 G26 G27 G28" _n
file write `macro1' "ERAS '_residualid' G21 G22 G23 G24 G25 G26 G27 G28 G29" _n
file write `macro1' "LINK 0 G21" _n
file write `macro1' "LINK 0 G22" _n
file write `macro1' "LINK 0 G23" _n
file write `macro1' "LINK 0 G24" _n
file write `macro1' "LINK 0 G25" _n
file write `macro1' "LINK 0 G26" _n
file write `macro1' "LINK 0 G27" _n
file write `macro1' "LINK 0 G28" _n
file write `macro1' "LINK 0 G29" _n
// Restore level ID column name
file write `macro1' "NAME '_`lev`level'id'' '`lev`level'id''" _n
if "`recode'" ~= "" {
file write `macro1' "MISR 1" _n
}
end
******************************************************************************
* MERGE THE SAVED MLWIN RESIDUALS BACK TO THE ORIGINAL STATA DATA SET
******************************************************************************
capture program drop runmlwin_mergeresiduals
program define runmlwin_mergeresiduals, sortpreserve
if _caller() >= 12 version 12.0
if _caller() <= 9 version 9.0
syntax , NUMLevels(int) TOUSE(string) HIER(string) [CC(string) LEVELID(string) SORT(string) MMID(string)] RESIDFile(string)
if "`sort'" ~= "" & "`cc'" == "" & `numlevels' > 1 {
display as error "Residuals cannot be reliably merged back when nosort is specified"
exit 198
}
local levXid_and_higher_levels `hier' // use this if hierarchical model
local levXid `levelid' // Use this if XC model
if `numlevels'==1 { // In single level models the observation/row number is the unique ID
quietly gen _residualid = _n if `touse' == 1
if _caller() >= 11 {
qui merge m:1 _residualid using "`residfile'", assert(master match) nogenerate noreport
}
else {
mata: setsortorder("`residfile'", "_residualid")
sort _residualid
quietly merge _residualid using "`residfile'", uniqusing //sort
assert inlist(_merge,1,3)
drop _merge
}
drop _residualid
}
else {
if ("`cc'"~="") { // cross-classified models in MCMC can have non-nested IDs in which case user must specify unique IDs
if "`mmid'" ~= "" {
local mmnum 1
foreach var of varlist `mmid' {
rename `levelid' _`levelid'
if "`levelid'" == "`var'" {
gen `levelid' = _`levelid'
}
else {
gen `levelid' = `var'
}
unab oldvars: *
if _caller() >= 11 {
qui merge m:1 `levelid' using "`residfile'", keep(master match) nogenerate noreport
}
else {
mata: setsortorder("`residfile'", "`levelid'")
sort `levelid'
quietly merge `levelid' using "`residfile'", uniqusing nokeep //sort
drop _merge
}
drop _residualid
unab allvars: *
local newvars :list allvars - oldvars
if _caller() >= 12 {
rename (`newvars') (=_`mmnum')
}
else {
foreach newvar of varlist `newvars' {
rename `newvar' `newvar'_`mmnum'
}
}
local ++mmnum
drop `levelid'
rename _`levelid' `levelid'
}
}
else {
if _caller() >= 11 {
qui merge m:1 `levelid' using "`residfile'", assert(master match) nogenerate noreport
}
else {
sort `levelid'
preserve
use "`residfile'", clear
quietly sort `levelid' // Make sure using data set is sorted
quietly compress
quietly save, replace
restore
quietly merge `levelid' using "`residfile'", uniqusing
assert inlist(_merge,1,3)
drop _merge
}
drop _residualid
}
}
else { // hierarchical models
quietly egen _residualid = group(`levXid_and_higher_levels') if `touse' == 1, missing
if _caller() >= 11 {
qui merge m:1 _residualid using "`residfile'", assert(master match) nogenerate noreport
}
else {
mata: setsortorder("`residfile'", "_residualid")
sort _residualid
quietly merge _residualid using "`residfile'", uniqusing //sort
assert inlist(_merge,1,3)
drop _merge
}
drop _residualid
}
}
end
******************************************************************************
* MERGE THE SAVED MLWIN FACTOR VALUES BACK TO THE ORIGINAL STATA DATA SET
******************************************************************************
capture program drop runmlwin_mergefscores
program define runmlwin_mergefscores, sortpreserve
if _caller() >= 12 version 12.0
if _caller() <= 9 version 9.0
syntax , NUMLevels(int) TOUSE(string) HIER(string) [CC(string) LEVELID(string) SORT(string)] FSCORESFile(string)
if "`sort'" ~= "" & "`cc'" == "" & `numlevels' > 1 {
display as error "Factor scores cannot be reliably merged back when nosort is specified"
exit 198
}
local levXid_and_higher_levels `hier' // use this if hierarchical model
local levXid `levelid' // Use this if XC model
if `numlevels'==1 { // In single level models the observation/row number is the unique ID
quietly gen _fscoreid = _n if `touse' == 1
}
else {
if ("`cc'"~="") { // cross-classified models in MCMC can have non-nested IDs in which case user must specify unique IDs
quietly egen _fscoreid = group(`levXid') if `touse' == 1, missing
}
else { // hierarchical models
quietly egen _fscoreid = group(`levXid_and_higher_levels') if `touse' == 1, missing
}
}
if _caller() >= 11 {
qui merge m:1 _fscoreid using "`fscoresfile'", assert(master match) nogenerate
}
else {
mata: setsortorder("`fscoresfile'", "_fscoreid")
sort _fscoreid
quietly merge _fscoreid using "`fscoresfile'", uniqusing //sort
assert inlist(_merge,1,3)
drop _merge
}
drop _fscoreid
end
******************************************************************************
* MERGE THE SAVED MLWIN SIMULATED RESPONSE BACK TO THE ORIGINAL STATA DATA SET
******************************************************************************
capture program drop runmlwin_mergesimresp
program define runmlwin_mergesimresp
if _caller() >= 12 version 12.0
if _caller() <= 9 version 9.0
syntax , SIMFile(string)
quietly gen _simprespid = _n
if _caller() >= 11 {
qui merge m:1 _simprespid using "`simfile'", assert(match) nogenerate
}
else {
quietly merge _simprespid using "`simfile'", uniqusing sort
assert inlist(_merge,3)
drop _merge
}
drop _simprespid
end
******************************************************************************
* SAVE THE MCMC CHAIN RESIDUALS AT EACH LEVEL TO STATA DATA SETS
******************************************************************************
capture program drop runmlwin_savechains
program define runmlwin_savechains
if _caller() >= 12 version 12.0
if _caller() <= 9 version 9.0
syntax , PARNames(string) CHAINResults(string)
use "`chainresults'", clear
local colnames
foreach var of local parnames {
if _caller() >= 11.1 {
local colnames `colnames' `=strtoname("`var'")'
}
else {
mata: st_local("tmpname", validname("`var'"))
local colnames `colnames', "`=abbrev("`tmpname'", 32)'"
}
}
if _caller() >= 11.1 {
qui putmata runmlwin_chains = (iteration `colnames' deviance), replace
}
else {
mata: runmlwin_chains = st_data(., ("iteration" `colnames', "deviance"))
}
end
******************************************************************************
* DISPLAY MODEL RESULTS
******************************************************************************
capture program drop runmlwin_display
program runmlwin_display
if _caller() >= 12 version 12.0
if _caller() <= 9 version 9.0
syntax [, OR IRr RRr SD CORrelations level(cilevel) CFORMAT(string) PFORMAT(string) SFORMAT(string) noHEADer noGRoup noCONTrast noFETable noRETable MOde MEdian ZRatio * ]
if "`or'" ~= "" & !inlist("`e(link)'", "logit", "ologit") {
display as error "Odds ratio options is only valid for univariate logit models"
exit 198
}
if "`rrr'" ~= "" & "`e(link)'" ~= "mlogit" {
display as error "Relative-rate ratio options is only valid for univariate unordered multinomial logit models"
exit 198
}
local tmppois "poisson"
local tmpnbin "nbinomial"
local distribution `e(distribution)'
if "`irr'" ~= "" & ~`:list tmppois in distribution' & ~`:list tmpnbin in distribution' {
display as error "Incidence-rate ratio options is only valid for univariate poisson or negative binomial models"
exit 198
}
if "`mode'" ~= "" & "`median'" ~= "" {
display as error "Only specify one of mode or median"
exit 198
}
if "`mode'" ~= "" {
local esttype mode
}
else if "`median'" ~= "" {
local esttype median
}
else {
local esttype mean
}
local eform = 0
if "`or'" ~= "" | "`irr'" ~= "" | "`rrr'" ~= "" local eform = 1
* Univariate or Multinomial model
local mtype univariate
if "`e(distribution)'" == "multinomial" {
local mtype multinomial
local respcats `e(respcategories)'
}
local responses `e(depvars)'
if `:list sizeof responses' > 1 local mtype multivariate
* Check whether all the responses in a multivariate response model are normal or not
local allnormal 1
local resptypes `e(distribution)'
forvalues n = 1/`:list sizeof resptypes' {
if "`=word("`resptypes'", `n')'" != "normal" {
local allnormal 0
}
}
* Info required for both the FP and RP tables
local multiplier = invnormal(1 - ((100 - `level')/2)/100) // returns 1.96 if 95% confidence intervals are requested
local kk = length("`level'") // as in 95 is two characters but 7 is only 1 character
* Numeric formatting
if `"`cformat'"' == "" {
local cformat `c(cformat)'
}
if `"`cformat'"' == "" {
local cformat %9.0g
}
if fmtwidth(`"`cformat'"') > 9 {
local cformat %9.0g
display as text "note: invalid cformat(), using default"
}
if `"`pformat'"' == "" {
local pformat `c(pformat)'
}
if `"`pformat'"' == "" {
local pformat %5.3f
}
if fmtwidth(`"`pformat'"') > 5 {
local pformat %5.3f
display as text "note: invalid pformat(), using default"
}
if `"`sformat'"' == "" {
local sformat `c(sformat)'
}
if `"`sformat'"' == "" {
local sformat %8.2f
}
if fmtwidth(`"`sformat'"') > 8 {
local sformat %8.2f
display as text "note: invalid sformat(), using default"
}
* Header
if "`header'" == "" {
* Display estimation algorithm
quietly runmlwin_verinfo $MLwiN_path // This may be inaccurate if the user changes version between estimation and display (or uses the plugin)
display as text "MLwiN `e(version)' multilevel model" _c
display _col(49) as text "Number of obs" _col(68) "=" _col(70) as result %9.0g e(N)
di as txt "`e(title)'"
display as text "Estimation algorithm: " as result "`e(method)'" _c
if ("`e(linearization)'"~="" & "`e(method)'"~="MCMC") display as text ", " as result "`e(linearization)'" _c
if "`group'" == "" {
display ""
* Display number of units at each level of the model hierarchy
if e(numlevels)>1 {
if "`e(N_g)'" == "" | "`e(g_min)'" == "" | "`e(g_avg)'" == "" | "`e(g_max)'" == "" {
display as error "Grouping information is not available when nosort option is specified"
}
else {
local ivars `e(ivars)'
local levels : list uniq ivars
tempname Ng min avg max
mat `Ng' = e(N_g)
mat `min' = e(g_min)
mat `avg' = e(g_avg)
mat `max' = e(g_max)
di
di as txt "{hline 16}{c TT}{hline 42}
di as txt _col(17) "{c |}" _col(21) "No. of" _continue
di as txt _col(34) "Observations per Group"
di as txt _col(2) "Level Variable" _col(17) "{c |}" _continue
di as txt _col(21) "Groups" _col(31) "Minimum" _continue
di as txt _col(42) "Average" _col(53) "Maximum"
di as txt "{hline 16}{c +}{hline 42}"
local i 1
foreach k of local levels {
local lev = abbrev("`k'",12)
local p = 16 - length("`lev'")
di as res _col(`p') "`lev'" _continue
di as txt _col(17) "{c |}" _continue
di as res _col(19) %8.0g `Ng'[1,`i'] _continue
di as res _col(29) %9.0g `min'[1,`i'] _continue
di as res _col(40) %9.1f `avg'[1,`i'] _continue
di as res _col(51) %9.0g `max'[1,`i']
local ++i
}
di as txt "{hline 16}{c BT}{hline 42}"
}
}
}
* Weights
// 1 2 3 4 5 6 7
//123456789012345678901234567890123456789012345678901234567890123456789012345678
//-------------------------------------------
// | Weights
// Level Variable | Variable Type
//----------------+--------------------------
// Unit ID | ############ standardised
//-------------------------------------------
local maxlevels = e(numlevels)
local hasweights = 0
forvalues l = `maxlevels'(-1)1 {
if "`e(weighttype`l')'" ~= "" {
local hasweights = 1
}
}
if "`hasweights'" == "1" {
local ivars `e(ivars)'
local levels : list uniq ivars
di
di as txt "{hline 16}{c TT}{hline 26}
di as txt _col(17) "{c |}" _col(29) "Weights"
di as txt _col(2) "Level Variable" _col(17) "{c |}" _col(23) "Variable" _col(40) "Type"
di as txt "{hline 16}{c +}{hline 26}"
forvalues l = `maxlevels'(-1)1 {
if "`e(weighttype`l')'" ~= "" {
if "`l'" ~= "1" {
local lev = abbrev(word("`levels'", `=`l'-1') ,12)
}
else {
local lev `e(level1var)'
}
local weightvar = abbrev("`e(weightvar`l')'",12)
local p1 = 16 - length("`lev'")
local p2 = 31 - length("`weightvar'")
local p3 = 44 - length("`e(weighttype`l')'")
di as res _col(`p1') "`lev'" as txt _col(17) "{c |}" as txt _col(`p2') "`weightvar'" _col(`p3') "`e(weighttype`l')'"
}
}
di as txt "{hline 16}{c BT}{hline 26}"
}
* Key for the contrasts in multinomial models
if "`contrast'" == "" {
if "`e(distribution)'"=="multinomial" {
display ""
local respcategories `e(respcategories)'
local basecategory `e(basecategory)'
local num_respcategories = wordcount("`respcategories'")
tokenize "`respcategories'"
local p = length("`respcategories'") + 4
if (`p'<20) local p = 20
di as txt "{hline 13}{c TT}{hline `p'}
di as txt _col(5) "Contrast" _col(14) "{c |}" ///
_col(16) "Log-odds"
di as txt "{hline 13}{c +}{hline 20}"
* Unordered
if "`e(link)'"=="mlogit" {
local contrast_denom = `basecategory'
local contrast_numers : list respcategories - contrast_denom
forvalues contrast = 1/`=`num_respcategories' - 1' {
local respcategory : word `contrast' of `contrast_numers'
local p = 13 - length("`contrast'")
di as res _col(`p') "`contrast'" as txt _col(14) "{c |}" _c
di as res " `respcategory'" as txt " vs. " as res "`basecategory'"
}
}
* Ordered
if "`e(link)'"=="ologit"|"`e(link)'"=="oprobit"|"`e(link)'"=="ocloglog" {
local first ""
forvalues contrast = 1/`=`num_respcategories' - 1' {
local first `first' `1'
macro shift
local rest `*'
* Lowest response value is basecategory
if (`basecategory'==real(word("`respcategories'",1))) {
local contrast_numer `rest'
local contrast_denom `first'
}
* Highest response value is basecategory
if (`basecategory'==real(word("`respcategories'",-1))) {
local contrast_numer `first'
local contrast_denom `rest'
}
local p = 13 - length("`contrast'")
di as res _col(`p') "`contrast'" as txt _col(14) "{c |}" _c
di as res " `contrast_numer'" as txt " vs. " as res "`contrast_denom'"
}
}
di as txt "{hline 13}{c BT}{hline 20}"
}
}
* Display model fit statistics
if e(method)=="IGLS" | e(method)=="RIGLS" {
display ""
if e(method)=="IGLS" & `allnormal'==1 {
display as text "Run time (seconds)" _col(22) "= " as result %10.2f e(time)
display as text "Number of iterations" _col(22) "= " as result %10.0f e(iterations)
if `e(ll)' == . {
display as error "Caution. MLwiN was unable to calculate the log likelihood, there may be a problem with your model"
}
display as text "Log likelihood" _col(22) "= " as result %10.0g e(ll)
display as text "Deviance" _col(22) "= " as result %10.0g -2*e(ll)
}
if e(method)=="RIGLS" & `allnormal'==1 {
display as text "Run time (seconds)" _col(27) "= " as result %10.2f e(time)
display as text "Number of iterations" _col(27) "= " as result %10.0f e(iterations)
if `e(ll)' == . {
display as error "Caution. MLwiN was unable to calculate the log likelihood, there may be a problem with your model"
}
display as text "Log restricted-likelihood" _col(27) "= " as result %10.0g e(ll)
display as text "Restricted-deviance" _col(27) "= " as result %10.0g -2*e(ll)
}
if `allnormal'==0 {
display as text "Run time (seconds)" _col(22) "= " as result %10.2f e(time)
display as text "Number of iterations" _col(22) "= " as result %10.0f e(iterations)
}
}
if e(method)=="MCMC" {
display ""
display as text "Burnin" _col(28) "= " as result %10.0f e(burnin)
display as text "Chain" _col(28) "= " as result %10.0f e(chain)
display as text "Thinning" _col(28) "= " as result %10.0f e(thinning)
display as text "Run time (seconds)" _col(28) "= " as result %10.3g e(time)
display as text "Deviance (dbar)" _col(28) "= " as result %10.2f e(dbar)
display as text "Deviance (thetabar)" _col(28) "= " as result %10.2f e(dthetabar)
display as text "Effective no. of pars (pd)" _col(28) "= " as result %10.2f e(pd)
display as text "Bayesian DIC" _col(28) "= " as result %10.2f e(dic)
}
}
* FE table
// IGLS FP
// 1 2 3 4 5 6 7
//123456789012345678901234567890123456789012345678901234567890123456789012345678
//------------------------------------------------------------------------------
// normexam | Coef. Std. Err. z P>|z| [95% Conf. Interval]
//-------------+----------------------------------------------------------------
// cons | ######### ######### ######## ###### ######### #########
//------------------------------------------------------------------------------
// MCMC FP
// 1 2 3 4 5 6 7
//123456789012345678901234567890123456789012345678901234567890123456789012345678
//------------------------------------------------------------------------------
// normexam | Mean Std. Dev. z ESS [95% Cred. Interval]
//-------------+----------------------------------------------------------------
// cons | ######### ######### ######## ###### ######### #########
//------------------------------------------------------------------------------
if "`fetable'" == "" {
* Display fixed part of the model
display as txt "{hline 13}{c TT}{hline 64}
if ("`mtype'"=="univariate") {
local p = 13 - length("`=abbrev("`e(depvars)'", 12)'")
if e(method)~="MCMC" {
display as txt _col(`p') "`=abbrev("`e(depvars)'", 12)'" _continue
display as txt _col(14) "{c |}" _continue
if "`or'" ~= "" display as txt _col(16) "Odds Ratio" _continue
if "`irr'" ~= "" display as txt _col(23) "IRR" _continue
if "`rrr'" ~= "" display as txt _col(23) "RRR" _continue
if "`or'" == "" & "`irr'" == "" & "`rrr'" == "" display as txt _col(21) "Coef." _continue
display as txt _col(29) "Std. Err." _continue
display as txt _col(44) "z" _continue
display as txt _col(49) "P>|z|" _continue
display as txt _col(`=61 -`kk'') `"[`=strsubdp("`level'")'% Conf. Interval]"'
}
if e(method)=="MCMC" {
display as txt _col(`p') "`=abbrev("`e(depvars)'", 12)'" _continue
display as txt _col(14) "{c |}" _continue
if "`or'" ~= "" display as txt _col(16) "Odds Ratio" _continue
if "`irr'" ~= "" display as txt _col(23) "IRR" _continue
if "`rrr'" ~= "" display as txt _col(23) "RRR" _continue
if "`or'" == "" & "`irr'" == "" & "`rrr'" == "" {
if "`esttype'" == "mode" {
display as txt _col(21) "Mode" _continue
}
else if "`esttype'" == "median" {
display as txt _col(19) "Median" _continue
}
else {
display as txt _col(21) "Mean" _continue
}
}
display as txt _col(29) "Std. Dev." _continue
if "`zratio'" ~= "" {
display as txt _col(44) "z" _continue
display as txt _col(49) "P>|z|" _continue
}
else {
display as txt _col(43) "ESS" _continue
display as txt _col(51) "P" _continue
}
display as txt _col(`=61 -`kk'') `"[`=strsubdp("`level'")'% Cred. Interval]"'
}
}
else {
if e(method)~="MCMC" {
display as txt _col(14) "{c |}" _continue
if "`or'" ~= "" display as txt _col(16) "Odds Ratio" _continue
if "`irr'" ~= "" display as txt _col(23) "IRR" _continue
if "`rrr'" ~= "" display as txt _col(23) "RRR" _continue
if "`or'" == "" & "`irr'" == "" & "`rrr'" == "" display as txt _col(21) "Coef." _continue
display as txt _col(29) "Std. Err." _continue
display as txt _col(44) "z" _continue
display as txt _col(49) "P>|z|" _continue
display as txt _col(`=61 -`kk'') `"[`=strsubdp("`level'")'% Conf. Interval]"'
}
if e(method)=="MCMC" {
display as txt _col(14) "{c |}" _continue
if "`or'" ~= "" display as txt _col(16) "Odds Ratio" _continue
if "`irr'" ~= "" display as txt _col(23) "IRR" _continue
if "`rrr'" ~= "" display as txt _col(23) "RRR" _continue
if "`or'" == "" & "`irr'" == "" & "`rrr'" == "" {
if "`esttype'" == "mode" {
display as txt _col(21) "Mode" _continue
}
else if "`esttype'" == "median" {
display as txt _col(19) "Median" _continue
}
else {
display as txt _col(21) "Mean" _continue
}
}
display as txt _col(29) "Std. Dev." _continue
if "`zratio'" ~= "" {
display as txt _col(44) "z" _continue
display as txt _col(49) "P>|z|" _continue
}
else {
display as txt _col(43) "ESS" _continue
display as txt _col(51) "P" _continue
}
display as txt _col(`=61 -`kk'') `"[`=strsubdp("`level'")'% Cred. Interval]"'
}
}
local eqlist :coleq e(b)
local eqlist :list uniq eqlist
local fpeqlist
foreach eq of local eqlist {
if (substr("`eq'",1,2)=="FP") local fpeqlist `fpeqlist' `eq'
}
local eqnum 0
foreach eq in `fpeqlist' {
local ++eqnum
display as txt "{hline 13}{c +}{hline 64}"
if "`mtype'" == "multivariate" && `:list sizeof responses' >= `eqnum'{
display as res "`=word("`responses'", `eqnum')'" as txt _col(14) "{c |}"
}
if "`mtype'" == "multinomial" && `:list sizeof respcats' > `eqnum'{
display as res "Contrast `eqnum'" as txt _col(14) "{c |}"
}
local r = 1
tempname matb matV matb_eq matV_eq
matrix `matb' = e(b)
matrix `matV' = e(V)
matrix `matb_eq' = `matb'[1, "`eq':"]
matrix `matV_eq' = `matV'["`eq':", "`eq':"]
if e(method)=="MCMC" & e(mcmcdiagnostics)==1 {
tempname matress matrlb matrub matrmode matrmedian matp matess_eq matlb_eq matub_eq matmode_eq matmedian_eq matp_eq
matrix `matress' = e(ess)
matrix `matrlb' = e(lb)
matrix `matrub' = e(ub)
matrix `matrmode' = e(mode)
matrix `matrmedian' = e(quantiles)
if "`esttype'" == "mode" {
matrix `matp' = e(pvalmode)
}
else if "`esttype'" == "median" {
matrix `matp' = e(pvalmedian)
}
else {
matrix `matp' = e(pvalmean)
}
matrix `matess_eq' = `matress'[1, "`eq':"]
matrix `matlb_eq' = `matrlb'[1, "`eq':"]
matrix `matub_eq' = `matrub'[1, "`eq':"]
matrix `matmode_eq' = `matrmode'[1, "`eq':"]
* NOTE: This assumes that the median is always the fifth row
matrix `matmedian_eq' = `matrmedian'[5, "`eq':"]
matrix `matp_eq' = `matp'[1, "`eq':"]
}
local names :colnames `matb_eq'
tempname matest matse matess matub matlb matmode matmedian
foreach par of local names {
* Exponential form
if `eform' == 1 {
* MCMC
if e(method)=="MCMC" & e(mcmcdiagnostics)==1 {
tempname expchain
local colnames :colfullnames e(b)
mata: `expchain' = exp(runmlwin_chains[ , `=`:list posof "`eq':`par'" in colnames' + 1'])
preserve
drop _all
label drop _all
if _caller() >= 11.1 {
getmata (`expchain') = `expchain'
}
else {
mata: (void) st_addvar("double", "`expchain'")
mata: st_addobs(length(`expchain'))
mata: st_store(., "`expchain'", `expchain')
}
mata: mata drop `expchain'
local posit = 0
if regexm("`par'", "(var)[\(]([a-zA-Z0-9_]+)[\)]") == 1 {
local posit = 1
}
quietly runmlwin_mcmcdiag `expchain', level(`level') thinning(e(thinning)) posit(`posit') `eform'
mat `matest' = r(mean)
mat `matse' = r(sd)*r(sd)
mat `matess' = r(ESS)
mat `matlb' = r(lb)
mat `matub' = r(ub)
mat `matmode' = r(mode)
mat `matmedian' = r(quantiles)
* NOTE: This assumes that the median is always the fifth row
mat `matmedian' = `matmedian'[5, 2]
restore
}
* (R)IGLS
else {
capture quietly nlcom (eform: exp([`eq']`par'))
if !_rc {
mat `matlb' = exp(`matb_eq'[1, `r'] - `multiplier'*sqrt(`matV_eq'[`r',`r']))
mat `matub' = exp(`matb_eq'[1, `r'] + `multiplier'*sqrt(`matV_eq'[`r',`r']))
mat `matest' = r(b)
mat `matse' = r(V)
mat `matest' = `matest'[1,1]
mat `matse' = `matse'[1,1]
}
else {
mat `matest' = J(1, 1, .)
mat `matse' = J(1, 1, .)
mat `matlb' = J(1, 1, .)
mat `matub' = J(1, 1, .)
}
}
}
* Original estimation metric
else {
matrix `matest' = `matb_eq'[1, `r']
matrix `matse' = `matV_eq'[`r',`r']
matrix `matlb' = `matest'[1, 1] - `multiplier'*sqrt(`matse'[1, 1])
matrix `matub' = `matest'[1, 1] + `multiplier'*sqrt(`matse'[1, 1])
* MCMC
if e(method)=="MCMC" & e(mcmcdiagnostics)==1 {
if e(level) == `level' {
mat `matess' = `matess_eq'[1, `r']
mat `matlb' = `matlb_eq'[1, `r']
mat `matub' = `matub_eq'[1, `r']
mat `matmode' = `matmode_eq'[1, `r']
mat `matmedian' = `matmedian_eq'[1, `r']
}
else {
tempname estchain
local colnames :colfullnames e(b)
mata: `estchain' = runmlwin_chains[ , `=`:list posof "`eq':`par'" in colnames' + 1']
preserve
drop _all
label drop _all
if _caller() >= 11.1 {
getmata (`estchain') = `estchain'
}
else {
mata: (void) st_addvar("double", "`estchain'")
mata: st_addobs(length(`estchain'))
mata: st_store(., "`estchain'", `estchain')
}
mata: mata drop `estchain'
local posit = 0
if regexm("`par'", "(var)[\(]([a-zA-Z0-9_]+)[\)]") == 1 {
local posit = 1
}
quietly runmlwin_mcmcdiag `estchain', level(`level') thinning(e(thinning)) posit(`posit') `eform'
mat `matlb' = r(lb)
mat `matub' = r(ub)
restore
}
}
}
local varname `=abbrev("`par'", 12)'
local k = length("`varname'")
local p = 13 - `k'
if e(method)~="MCMC" {
display as txt _col(`p') "`varname'" _continue
display as txt _col(14) "{c |}" _continue
display as res _col(17) `cformat' `matest'[1, 1] _continue
display as res _col(28) `cformat' sqrt(`matse'[1, 1]) _continue
display as res _col(36) `sformat' `matb_eq'[1, `r']/sqrt(`matV_eq'[`r',`r']) _continue
display as res _col(49) `pformat' 2*normal(-abs(`matb_eq'[1, `r']/sqrt(`matV_eq'[`r',`r']))) _continue
display as res _col(58) `cformat' `matlb'[1, 1] _continue
display as res _col(70) `cformat' `matub'[1, 1]
local ++r
}
if e(method)=="MCMC" {
display as txt _col(`p') "`varname'" _continue
display as txt _col(14) "{c |}" _continue
if "`esttype'" == "mode" {
display as res _col(17) `cformat' `matmode'[1, 1] _continue
}
else if "`esttype'" == "median" {
display as res _col(17) `cformat' `matmedian'[1, 1] _continue
}
else {
display as res _col(17) `cformat' `matest'[1, 1] _continue
}
display as res _col(28) `cformat' sqrt(`matse'[1, 1]) _continue
if "`zratio'" == "" {
if e(mcmcdiagnostics)==1 {
display as res _col(36) %9.0f `matess'[1, 1] _continue
display as res _col(49) `pformat' `matp_eq'[1, `r'] _continue
}
}
else {
display as res _col(36) `sformat' `matb_eq'[1, `r']/sqrt(`matV_eq'[`r',`r']) _continue
display as res _col(49) `pformat' 2*normal(-abs(`matb_eq'[1, `r']/sqrt(`matV_eq'[`r',`r']))) _continue
}
if e(mcmcdiagnostics)==1 {
display as res _col(58) `cformat' `matlb'[1, 1] _continue
display as res _col(70) `cformat' `matub'[1, 1] _continue
}
display
local ++r
}
}
}
display as txt "{hline 13}{c BT}{hline 64}"
display ""
}
* RE Table
// IGLS RP
// 1 2 3 4 5 6 7
//123456789012345678901234567890123456789012345678901234567890123456789012345678
//------------------------------------------------------------------------------
// Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval]
//-----------------------------+------------------------------------------------
//Level 1: |
// var(cons) | ######### ######### ######### #########
//------------------------------------------------------------------------------
//
// MCMC RP
// 1 2 3 4 5 6 7
//123456789012345678901234567890123456789012345678901234567890123456789012345678
//------------------------------------------------------------------------------
// Random-effects Parameters | Mean Std. Dev. ESS [95% Cred. Int.]
//-----------------------------+------------------------------------------------
//Level 1: |
// var(cons) | ######### ######### ###### ######### #########
//------------------------------------------------------------------------------
if "`retable'" == "" {
* Display random part of the model
local eqlist :coleq e(b)
local eqlist :list uniq eqlist
* Random part design
* Change display order so parameter estimates are reported at their appropriate level
local ordlist
local i = 1
local maxlevels = e(numlevels)
forvalues l = `maxlevels'(-1)1 {
if `:list posof "RP`l'" in eqlist' > 0 {
local ordlist `ordlist' RP`l'
local ordlist`i'des Level `l': `e(level`l')'
if `l' == 1 {
local ordlist`i'des `ordlist`i'des' `e(level1var)'
}
else {
local levnames `e(ivars)'
local namepos = `:list sizeof levnames' - `=`l' - 2'
local ordlist`i'des `ordlist`i'des' `:word `namepos' of `e(ivars)''
}
local ++i
}
if `:list posof "RP`l'D" in eqlist' > 0 {
local ordlist `ordlist' RP`l'D
local ordlist`i'des Level `l': `e(level`l')' (Design)
local ++i
}
if `:list posof "RP`l'FL" in eqlist' > 0 {
local ordlist `ordlist' RP`l'FL
local ordlist`i'des Level `l' factors:
local ++i
}
if `:list posof "RP`l'FV" in eqlist' > 0 {
local ordlist `ordlist' RP`l'FV
local ordlist`i'des Level `l' factor covariances:
local ++i
}
}
local i = 1
local iLAST = wordcount("`ordlist'")
// Here we want to display the random part only if variables have been specified
if `:list sizeof ordlist' > 0 {
display as txt "{hline 29}{c TT}{hline 48}"
if e(method)~="MCMC" {
display as txt _col(4) "Random-effects Parameters" _continue
display as txt _col(30) "{c |}" _continue
display as txt _col(34) "Estimate" _continue
display as txt _col(45) "Std. Err." _continue
display as txt _col(`=61 -`kk'') `"[`=strsubdp("`level'")'% Conf. Interval]"'
}
if e(method)=="MCMC" {
display as txt _col(4) "Random-effects Parameters" _continue
display as txt _col(30) "{c |}" _continue
if "`esttype'" == "mode" {
display as txt _col(36) "Mode" _continue
}
else if "`esttype'" == "median" {
display as txt _col(34) "Median" _continue
}
else {
display as txt _col(36) "Mean" _continue
}
display as txt _col(43) "Std. Dev." _continue
display as txt _col(55) "ESS" _continue
display as txt _col(`=65 -`kk'') `"[`=strsubdp("`level'")'% Cred. Int]"'
}
foreach eq in `ordlist' {
local r = 1
display as txt "{hline 29}{c +}{hline 48}"
display as res "`ordlist`i'des'" _continue
display as txt _col(30) "{c |}"
tempname matb matv matb_eq matV_eq
matrix `matb' = e(b)
matrix `matV' = e(V)
matrix `matb_eq' = `matb'[1, "`eq':"]
matrix `matV_eq' = `matV'["`eq':", "`eq':"]
if e(method)=="MCMC" & e(mcmcdiagnostics)==1 {
tempname matress matrlb matrub matrmode matrmedian matess_eq matlb_eq matub_eq matmode_eq matmedian_eq
matrix `matress' = e(ess)
matrix `matrlb' = e(lb)
matrix `matrub' = e(ub)
matrix `matrmode' = e(mode)
matrix `matrmedian' = e(quantiles)
matrix `matess_eq' = `matress'[1, "`eq':"]
matrix `matlb_eq' = `matrlb'[1, "`eq':"]
matrix `matub_eq' = `matrub'[1, "`eq':"]
matrix `matmode_eq' = `matrmode'[1, "`eq':"]
* NOTE: This assumes that the median is always the fifth row
matrix `matmedian_eq' = `matrmedian'[5, "`eq':"]
}
tempname matest matse matess matlb matub matmode matmedian
local names :colnames `matb_eq'
foreach par of local names {
// Don't display when constrained to one
local partype
if regexm("`par'", "(var)[\(]([a-zA-Z0-9_]+)[\)]") == 1 {
local partype `=regexs(1)'
local par1 `=regexs(2)'
local par2 `=regexs(2)'
}
if regexm("`par'", "(cov)[\(]([a-zA-Z0-9_~]+)[\\]([a-zA-Z0-9_~]+)[\)]") == 1 {
local partype `=regexs(1)'
local par1 `=regexs(2)'
local par2 `=regexs(3)'
}
* SD specified and parameter is a variance
if "`sd'" ~= "" & "`partype'" == "var" {
* MCMC
if e(method)=="MCMC" & e(mcmcdiagnostics)==1 {
tempname sdchain
local colnames :colfullnames e(b)
mata: `sdchain' = sqrt(runmlwin_chains[ , `=`:list posof "`eq':`par'" in colnames' + 1'])
preserve
drop _all
label drop _all
if _caller() >= 11.1 {
getmata (`sdchain') = `sdchain'
}
else {
mata: (void) st_addvar("double", "`sdchain'")
mata: st_addobs(length(`sdchain'))
mata: st_store(., "`sdchain'", `sdchain')
}
mata: mata drop `sdchain'
quietly runmlwin_mcmcdiag `sdchain', level(`level') thinning(e(thinning)) posit(1) `eform'
mat `matest' = r(mean)
mat `matse' = r(sd)*r(sd)
mat `matess' = r(ESS)
mat `matlb' = r(lb)
mat `matub' = r(ub)
mat `matmode' = r(mode)
mat `matmedian' = r(quantiles)
* NOTE: This assumes that the median is always the fifth row
mat `matmedian' = `matmedian'[5, 2]
restore
}
* (R)IGLS
else {
capture quietly nlcom (sd: sqrt([`eq']`par'))
if !_rc {
mat `matlb' = sqrt(`matb_eq'[1, `r'] - `multiplier'*sqrt(`matV_eq'[`r',`r']))
mat `matub' = sqrt(`matb_eq'[1, `r'] + `multiplier'*sqrt(`matV_eq'[`r',`r']))
mat `matest' = r(b)
mat `matse' = r(V)
mat `matest' = `matest'[1,1]
mat `matse' = `matse'[1,1]
}
else {
mat `matest' = J(1, 1, .)
mat `matse' = J(1, 1, .)
mat `matlb' = J(1, 1, .)
mat `matub' = J(1, 1, .)
}
}
local par sd(`=abbrev("`par1'", 24)')
}
* CORR specified and parameter is a covariance
else if "`correlations'" ~= "" & "`partype'" == "cov" {
// Search for corresponding variances
local possnames :colnames `matb_eq'
local posspar1
local posspar2
foreach poss of local possnames {
if regexm("`poss'", "(var)[\(]([a-zA-Z0-9_]+)[\)]") == 1 {
local poss `=regexs(2)'
if "`=abbrev("`poss'", 12)'" == "`par1'" local posspar1 `posspar1' "`poss'"
if "`=abbrev("`poss'", 12)'" == "`par2'" local posspar2 `posspar2' "`poss'"
}
}
if `:list sizeof posspar1' != 1 | `:list sizeof posspar2' != 1 {
display as error "Unable to identify unique variance parameters for `par'. Don't use the corr option."
exit 198
}
else {
local par1 `posspar1'
local par2 `posspar2'
}
* MCMC
if e(method)=="MCMC" & e(mcmcdiagnostics)==1 {
// Recalculate correlations from chain (There is an argument that the estimates and covariance should always be estimated from chains on display too, otherwise there is potential for wasted calculations)
tempname corrchain
local colnames :colfullnames e(b)
mata: `corrchain' = runmlwin_chains[ , `=`:list posof "`eq':`par'" in colnames' + 1'] :/ sqrt(runmlwin_chains[ , `=`:list posof "`eq':var(`=abbrev("`par1'", 24)')" in colnames' + 1'] :* runmlwin_chains[ , `=`:list posof "`eq':var(`=abbrev("`par2'", 24)')" in colnames' + 1'])
preserve
drop _all
label drop _all
if _caller() >= 11.1 {
getmata (`corrchain') = `corrchain'
}
else {
mata: (void) st_addvar("double", "`corrchain'")
mata: st_addobs(length(`corrchain'))
mata: st_store(., "`corrchain'", `corrchain')
}
mata: mata drop `corrchain'
quietly runmlwin_mcmcdiag `corrchain', level(`level') thinning(e(thinning)) posit(0) `eform'
mat `matest' = r(mean)
mat `matse' = r(sd)*r(sd)
mat `matess' = r(ESS)
mat `matlb' = r(lb)
mat `matub' = r(ub)
mat `matmode' = r(mode)
mat `matmedian' = r(quantiles)
* NOTE: This assumes that the median is always the fifth row
mat `matmedian' = `matmedian'[5, 2]
restore
}
* (R)IGLS
else {
capture quietly nlcom (corr: [`eq']`par' / sqrt([`eq']var(`=abbrev("`par1'", 24)') * [`eq']var(`=abbrev("`par2'", 24)')))
if !_rc {
mat `matest' = r(b)
mat `matse' = r(V)
mat `matest' = `matest'[1,1]
mat `matse' = `matse'[1,1]
mat `matlb' = `matest'[1, 1] - `multiplier'*sqrt(`matse'[1, 1]) // This is correct! This is what xtmixed does when reporting covariance parameters (as opposed to the default of correlation parameters)
mat `matub' = `matest'[1, 1] + `multiplier'*sqrt(`matse'[1, 1]) // So we can't simply transform the lb and up of the correlation parameter on the original estimation metric because the correlation is a non-linear combination of several other parameters.
}
else {
mat `matest' = J(1, 1, .)
mat `matse' = J(1, 1, .)
mat `matlb' = J(1, 1, .)
mat `matub' = J(1, 1, .)
}
}
local par corr(`=abbrev("`par1'", 12)',`=abbrev("`par2'", 12)')
}
* All other cases
else {
if "`partype'" == "var" {
local par var(`=abbrev("`par1'", 24)')
}
if "`partype'" == "cov" {
local par cov(`=abbrev("`par1'", 12)',`=abbrev("`par2'", 12)')
}
mat `matest' = `matb_eq'[1, `r']
mat `matse' = `matV_eq'[`r',`r']
mat `matlb' = `matest'[1, 1] - `multiplier'*sqrt(`matse'[1, 1])
mat `matub' = `matest'[1, 1] + `multiplier'*sqrt(`matse'[1, 1])
if e(method)=="MCMC" & e(mcmcdiagnostics)==1 {
mat `matess' = `matess_eq'[1, `r']
mat `matmode' = `matmode_eq'[1, `r']
mat `matmedian' = `matmedian_eq'[1, `r']
if e(level) == `level' {
mat `matlb' = `matlb_eq'[1, `r']
mat `matub' = `matub_eq'[1, `r']
}
else {
tempname estchain
local colnames :colfullnames e(b)
mata: `estchain' = runmlwin_chains[ , `=`:list posof "`eq':`par'" in colnames' + 1']
preserve
drop _all
label drop _all
if _caller() >= 11.1 {
getmata (`estchain') = `estchain'
}
else {
mata: (void) st_addvar("double", "`estchain'")
mata: st_addobs(length(`estchain'))
mata: st_store(., "`estchain'", `estchain')
}
mata: mata drop `estchain'
local posit = 0
if "`partype'" == "var" {
local posit = 1
}
quietly runmlwin_mcmcdiag `estchain', level(`level') thinning(e(thinning)) posit(`posit') `eform'
mat `matlb' = r(lb)
mat `matub' = r(ub)
restore
}
}
}
local k = length("`par'")
local p = 29 - `k'
if e(method)~="MCMC" {
display as txt _col(`p') "`par'" _continue
display as txt _col(30) "{c |}" _continue
display as res _col(33) `cformat' `matest'[1, 1] _continue
display as res _col(44) `cformat' sqrt(`matse'[1,1]) _continue
display as res _col(58) `cformat' `matlb'[1, 1] _continue
display as res _col(70) `cformat' `matub'[1, 1]
local ++r
}
if e(method)=="MCMC" {
display as txt _col(`p') "`par'" _continue
display as txt _col(30) "{c |}" _continue
if "`esttype'" == "mode" {
display as res _col(32) `cformat' `matmode'[1, 1] _continue
}
else if "`esttype'" == "median" {
display as res _col(32) `cformat' `matmedian'[1, 1] _continue
}
else {
display as res _col(32) `cformat' `matest'[1, 1] _continue
}
display as res _col(42) `cformat' sqrt(`matse'[1, 1]) _continue
if e(mcmcdiagnostics)==1 {
display as res _col(52) %6.0f `matess'[1, 1] _continue
display as res _col(60) `cformat' `matlb'[1, 1] _continue
display as res _col(70) `cformat' `matub'[1, 1] _continue
}
display
local ++r
}
}
local ++i
}
* Horizontal line
display as txt "{hline 29}{c BT}{hline 48}
}
}
local eqlist :coleq e(b)
local eqlist :list uniq eqlist
if `:list posof "OD" in eqlist' > 0 & ("`e(extrabinomial)'" == "1" | "`e(distribution)'" == "nbinomial") {
tempname matodb_eq
tempname matodV_eq
matrix `matodb_eq' = `matb'[1, "OD:"]
matrix `matodV_eq' = `matV'["OD:", "OD:"]
if e(method)=="MCMC" & e(mcmcdiagnostics)==1 {
tempname matress matrlb matrub matrmode matrmedian matodess_eq matodlb_eq matodub_eq matodmode_eq matodmedian_eq
matrix `matress' = e(ess)
matrix `matrlb' = e(lb)
matrix `matrub' = e(ub)
matrix `matrmode' = e(mode)
matrix `matrmedian' = e(quantiles)
matrix `matodess_eq' = `matress'[1, "OD:"]
matrix `matodlb_eq' = `matrlb'[1, "OD:"]
matrix `matodub_eq' = `matrub'[1, "OD:"]
matrix `matodmode_eq' = `matrmode'[1, "OD:"]
* NOTE: This assumes that the median is always the fifth row
matrix `matodmedian_eq' = `matrmedian'[5, "OD:"]
}
local names :colnames `matodb_eq'
display as txt "{hline 29}{c TT}{hline 48}"
if e(method)~="MCMC" {
display as txt _col(4) "Overdispersion Parameters" _continue
display as txt _col(30) "{c |}" _continue
display as txt _col(34) "Estimate" _continue
display as txt _col(45) "Std. Err." _continue
display as txt _col(`=61 -`kk'') `"[`=strsubdp("`level'")'% Conf. Interval]"'
}
if e(method)=="MCMC" {
display as txt _col(4) "Overdispersion Parameters" _continue
display as txt _col(30) "{c |}" _continue
if "`esttype'" == "mode" {
display as txt _col(36) "Mode" _continue
}
else if "`esttype'" == "median" {
display as txt _col(34) "Median" _continue
}
else {
display as txt _col(36) "Mean" _continue
}
display as txt _col(43) "Std. Dev." _continue
display as txt _col(55) "ESS" _continue
display as txt _col(`=65 -`kk'') `"[`=strsubdp("`level'")'% Cred. Int]"'
}
* Horizontal line
display as txt "{hline 29}{c +}{hline 48}
local r = 1
foreach par of local names {
if "`r'" == "1" && "`e(distribution)'" == "nbinomial" && "`e(extrabinomial)'" == "0" {
local ++r
continue
}
mat `matest' = `matodb_eq'[1, `r']
mat `matse' = `matodV_eq'[`r',`r']
mat `matlb' = `matest'[1, 1] - `multiplier'*sqrt(`matse'[1, 1])
mat `matub' = `matest'[1, 1] + `multiplier'*sqrt(`matse'[1, 1])
if e(method)=="MCMC" & e(mcmcdiagnostics)==1 {
mat `matess' = `matodess_eq'[1, `r']
mat `matmode' = `matodmode_eq'[1, `r']
mat `matmedian' = `matodmedian_eq'[1, `r']
if e(level) == `level' {
mat `matlb' = `matodlb_eq'[1, `r']
mat `matub' = `matodub_eq'[1, `r']
}
else {
tempname estchain
local colnames :colfullnames e(b)
mata: `estchain' = runmlwin_chains[ , `=`:list posof "`eq':`par'" in colnames' + 1']
preserve
drop _all
label drop _all
if _caller() >= 11.1 {
getmata (`estchain') = `estchain'
}
else {
mata: (void) st_addvar("double", "`estchain'")
mata: st_addobs(length(`estchain'))
mata: st_store(., "`estchain'", `estchain')
}
mata: mata drop `estchain'
local posit = 0
if "`partype'" == "var" {
local posit = 1
}
quietly runmlwin_mcmcdiag `estchain', level(`level') thinning(e(thinning)) posit(`posit') `eform'
mat `matlb' = r(lb)
mat `matub' = r(ub)
restore
}
}
local k = length("`par'")
local p = 29 - `k'
if e(method)~="MCMC" {
display as txt _col(`p') "`par'" _continue
display as txt _col(30) "{c |}" _continue
display as res _col(33) `cformat' `matest'[1, 1] _continue
display as res _col(44) `cformat' sqrt(`matse'[1,1]) _continue
display as res _col(58) `cformat' `matlb'[1, 1] _continue
display as res _col(70) `cformat' `matub'[1, 1]
local ++r
}
if e(method)=="MCMC" {
display as txt _col(`p') "`par'" _continue
display as txt _col(30) "{c |}" _continue
if "`esttype'" == "mode" {
display as res _col(32) `cformat' `matmode'[1, 1] _continue
}
else if "`esttype'" == "median" {
display as res _col(32) `cformat' `matmedian'[1, 1] _continue
}
else {
display as res _col(32) `cformat' `matest'[1, 1] _continue
}
display as res _col(42) `cformat' sqrt(`matse'[1, 1]) _continue
if e(mcmcdiagnostics)==1 {
display as res _col(52) %6.0f `matess'[1, 1] _continue
display as res _col(60) `cformat' `matlb'[1, 1] _continue
display as res _col(70) `cformat' `matub'[1, 1] _continue
}
display
local ++r
}
}
* Horizontal line
display as txt "{hline 29}{c BT}{hline 48}
}
* Notes
if e(method)=="IGLS" & e(converged)==0 display _n as error "WARNING: IGLS algorithm failed to converge. Increase the number of iterations. See the maxiterations() option."
if e(method)=="RIGLS" & e(converged)==0 display _n as error "WARNING: RIGLS algorithm failed to converge. Increase the number of iterations. See the maxiterations() option."
end
mata:
void copymatrix(string scalar srcmat, string scalar srcrowidxmat, string scalar srccolidxmat, string scalar destmat, string scalar destrowidxmat, string scalar destcolidxmat) {
real colvector srcrowidx;
real colvector srccolidx;
real colvector destrowidx;
real colvector destcolidx;
real matrix temp;
src = st_matrix(srcmat);
srcrowidx = st_matrix(srcrowidxmat);
srccolidx = st_matrix(srccolidxmat);
destrowidx = st_matrix(destrowidxmat);
destcolidx = st_matrix(destcolidxmat);
temp = st_matrix(destmat);
temp[destrowidx, destcolidx] = src[srcrowidx, srccolidx];
_editmissing(temp, 0);
st_replacematrix(destmat, temp);
}
void checkmm(string scalar mmids, string scalar mmweights) {
real matrix idmat;
real matrix weightmat;
real scalar i;
real scalar j;
if (callersversion() >= 11) {
st_view(idmat, ., mmids);
st_view(weightmat, ., mmweights);
}
else {
st_view(idmat, ., tokens(mmids));
st_view(weightmat, ., tokens(mmweights));
}
for (i = 1; i <= rows(idmat); i++) {
for (j = 1; j <= cols(idmat); j++) {
if (idmat[i, j] != 0 && weightmat[i, j] == 0) {
errprintf("Error, the MM ID variable %g for observation %g is present but a zero MM weight has been specified for it\n", j, i);
exit(198);
}
if (idmat[i, j] == 0 && weightmat[i, j] != 0) {
errprintf("Error, the MM ID variable %g for observation %g is absent but a positive MM weight has been specified for it\n", j, i);
exit(198);
}
for (k = 1; k < j; k++) {
if (idmat[i, j] == idmat[i, k] & idmat[i, j] != 0 ) {
errprintf("Error, the MM ID variable %g for observation %g is duplicated in ID variable %g\n", j, i, k);
exit(198)
}
}
}
}
}
string scalar validname(string scalar name) {
real rowvector codes;
real rowvector valid;
codes = ascii(name);
if (codes[1] >=48 && codes[1] <= 57) {
codes = 95,codes;
}
valid = ((codes:<65 :| codes:>90) :& (codes:<97 :| codes:>122) :& (codes:<48 :| codes:>57))
for (i = 1; i <= length(codes); i++) {
if (valid[i] == 1) {
codes[i] = 95;
}
}
// The following code would truncate to 32 characters
/*
if (length(codes) > 32) {
codes = codes[|1\32|];
}
*/
name = char(codes);
return(name);
}
void setsortorder(string scalar filename, string scalar sortorder) {
string rowvector order;
colvector C;
real scalar fh;
real scalar ver;
real scalar byteorder;
real scalar filetype;
real scalar reserved;
real scalar nvar;
real scalar nobs;
real scalar varlablen;
real scalar colnamelen;
real scalar recsize;
real scalar fmtlen;
string scalar data_label;
string scalar time_stamp;
real matrix typlist;
string matrix varlist;
real matrix srtlist;
order = tokens(sortorder);
if (length(order) ~= length(uniqrows(order))) {
printf("variables in sort order are not unique");
exit(error(198));
}
C = bufio()
fh = fopen(filename, "rw");
ver = fbufget(C, fh, "%1bu");
byteorder = fbufget(C, fh, "%1bu");
filetype = fbufget(C, fh, "%1bu");
reserved = fbufget(C, fh, "%1bu");
bufbyteorder(C, byteorder);
if (filetype != 1)
{
display("Invalid file");
exit();
}
nvar = fbufget(C, fh, "%2bu");
nobs = fbufget(C, fh, "%4bu");
if (ver == 105) {
varlablen = 32;
colnamelen = 9;
recsize = 16;
fmtlen = 12;
}
else if (ver == 108) {
varlablen = 81;
colnamelen = 9;
recsize = 16;
fmtlen = 12;
}
else if (ver == 110 || ver == 111 || ver == 113) {
varlablen = 81;
colnamelen = 33;
recsize = 32;
fmtlen = 12;
}
else if (ver == 114 || ver == 115) {
varlablen = 81;
colnamelen = 33;
recsize = 32;
fmtlen = 49;
}
else {
printf("unrecognised version: %f", ver);
exit(error(198));
}
data_label = fbufget(C, fh, "%" + strofreal(varlablen) + "s");
time_stamp = fbufget(C, fh, "%18s");
typlist = J(1, nvar, .);
for (i = 1; i <= nvar; i++) {
typlist[1, i] = fbufget(C, fh, "%1bu");
}
varlist = J(1, nvar, "");
for (i = 1; i <= nvar; i++) {
varlist[1, i] = fbufget(C, fh, "%" + strofreal(colnamelen) + "s");
}
srtlist = J(1, nvar + 1, 0);
for (i = 1; i <= length(order); i++) {
found = 0;
for (j = 1; j <= nvar; j++) {
if (order[i] == varlist[1, j]) {
found = 1;
srtlist[1, i] = j
break;
}
}
if (found == 0) {
printf("variable %s not found", order[i]);
exit(error(198));
}
}
for (i = 1; i <= nvar + 1; i++) {
fbufput(C, fh, "%2bu", srtlist[1, i]);
}
fclose(fh);
}
end
******************************************************************************
exit