Dear Statalisters,
I would like to point out a problem that occurs when running various rather complex Bayesian multilevel models using the bayes prefix. The problem occurs in both Stata 17/MP and Stata 18/MP for Unix (Linux 64-bit x86-64). A test data set for such a case is attached. Here I try to run a crossed effects model with two ('fixed') covariates, two crossed effects and random coefficients ('slopes') in both levels (using Isabel’s suggestions presented here Including covariates in crossed-effects models).
Here’s the code I’m using
The code crashes Stata after estimating the model, i.e. Stata exits without saving the estimates, here is the log-file.
My guess is that the problem must have to do with the size of the burn-in period: I have tried
burnin(20000), burnin(100000) and burnin(150000).
For all three sizes, everything works fine.
Has anyone else encountered similar problems and successfully fixed them?
Just in case: at the moment, I’m using Nikolay’s workaround to retrieve estimation results (described here Retrieving estimation results from "bayes [...] saving(xyz)" command):
Best,
Alex
I would like to point out a problem that occurs when running various rather complex Bayesian multilevel models using the bayes prefix. The problem occurs in both Stata 17/MP and Stata 18/MP for Unix (Linux 64-bit x86-64). A test data set for such a case is attached. Here I try to run a crossed effects model with two ('fixed') covariates, two crossed effects and random coefficients ('slopes') in both levels (using Isabel’s suggestions presented here Including covariates in crossed-effects models).
Here’s the code I’m using
Code:
*--------------------------------------------------*
cap log close
log using investigate_bayes_mixed, replace
use test_data, clear
*Including covariates in crossed-effects models as described by Isabel Canette
*https://blog.stata.com/2010/12/22/including-covariates-in-crossed-effects-models/
qui tab Group2, gen(id_Group2)
unab idvar: id_Group2*
foreach v of local idvar {
gen inter`v' =covariate2*`v'
}
*--------------------------------------------------*
*Standard frequentist model
mixed outcome covariate1 covariate2 || _all: R.Group2 || _all: inter*, cov(identity) nocons || Group1: covariate2
*--------------------------------------------------*
bayes, burnin(200000) dots(1000) rseed(42) saving(investigate_bayes_mixed_simulation, replace): ///
mixed outcome covariate1 covariate2 || _all: R.Group2 || _all: inter*, cov(identity) nocons || Group1: covariate2
estimates save investigate_bayes_mixed, replace
*--------------------------------------------------*
log close
Code:
. use test_data, clear
.
. *Including covariates in crossed-effects models as described by Isabel Canette
. *https://blog.stata.com/2010/12/22/including-covariates-in-crossed-effects-models/
.
. qui tab Group2, gen(id_Group2)
. unab idvar: id_Group2*
. foreach v of local idvar {
2. gen inter`v' =covariate2*`v'
3. }
.
. *--------------------------------------------------*
. *Standard frequentist model
. mixed outcome covariate1 covariate2 || _all: R.Group2 || _all: inter*, cov(identity) nocons || Group1: covariate2
Performing EM optimization ...
Performing gradient-based optimization:
Iteration 0: Log likelihood = 37.598276
Iteration 1: Log likelihood = 38.546812
Iteration 2: Log likelihood = 38.905076
Iteration 3: Log likelihood = 38.952799
Iteration 4: Log likelihood = 38.954077
Iteration 5: Log likelihood = 38.954091
Iteration 6: Log likelihood = 38.954091
Computing standard errors ...
Mixed-effects ML regression Number of obs = 97
Grouping information
-------------------------------------------------------------
| No. of Observations per group
Group variable | groups Minimum Average Maximum
----------------+--------------------------------------------
_all | 1 97 97.0 97
Group1 | 24 1 4.0 20
-------------------------------------------------------------
Wald chi2(2) = 14.75
Log likelihood = 38.954091 Prob > chi2 = 0.0006
------------------------------------------------------------------------------
outcome | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
covariate1 | -.0130288 .004734 -2.75 0.006 -.0223072 -.0037503
covariate2 | -.2171314 .066152 -3.28 0.001 -.3467869 -.0874759
_cons | .7379007 .0734132 10.05 0.000 .5940134 .881788
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects parameters | Estimate Std. err. [95% conf. interval]
-----------------------------+------------------------------------------------
_all: Identity |
var(R.Group2) | .0233182 .0101985 .0098949 .0549513
-----------------------------+------------------------------------------------
_all: Identity |
var(inter~21..inte~219)(1) | 2.66e-11 . . .
-----------------------------+------------------------------------------------
Group1: Independent |
var(covariate2) | .0042044 .0217997 1.62e-07 108.9365
var(_cons) | 1.68e-16 1.27e-12 0 .
-----------------------------+------------------------------------------------
var(Residual) | .0187106 .0032308 .0133385 .0262462
------------------------------------------------------------------------------
LR test vs. linear model: chi2(4) = 49.66 Prob > chi2 = 0.0000
(1) interid_Group21 interid_Group22 interid_Group23 interid_Group24 interid_Group25 interid_Group26 interid_Group27 interid_Group28 interid_Group29 interid_Group210 interid_Group211 interid_Group212 interid_Group213
interid_Group214 interid_Group215 interid_Group216 interid_Group217 interid_Group218 interid_Group219
Note: LR test is conservative and provided only for reference.
.
. *--------------------------------------------------*
. *Bayesian model, with burnin(200000)
. bayes, burnin(200000) dots(1000) rseed(42) saving(investigate_bayes_mixed_simulation, replace): ///
> mixed outcome covariate1 covariate2 || _all: R.Group2 || _all: inter*, cov(identity) nocons || Group1: covariate2
note: Gibbs sampling is used for regression coefficients and some variance components.
Burn-in 200000 aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa done
Simulation 10000 .......... done
file investigate_bayes_mixed_simulation.dta not found; file saved.
Multilevel structure
------------------------------------------------------------------------------
Group2
{U0}: random intercepts
_all
{V0}: random intercepts
{V1}: random coefficients for interid_Group21
{V2}: random coefficients for interid_Group22
{V3}: random coefficients for interid_Group23
{V4}: random coefficients for interid_Group24
{V5}: random coefficients for interid_Group25
{V6}: random coefficients for interid_Group26
{V7}: random coefficients for interid_Group27
{V8}: random coefficients for interid_Group28
{V9}: random coefficients for interid_Group29
{V10}: random coefficients for interid_Group210
{V11}: random coefficients for interid_Group211
{V12}: random coefficients for interid_Group212
{V13}: random coefficients for interid_Group213
{V14}: random coefficients for interid_Group214
{V15}: random coefficients for interid_Group215
{V16}: random coefficients for interid_Group216
{V17}: random coefficients for interid_Group217
{V18}: random coefficients for interid_Group218
{V19}: random coefficients for interid_Group219
Group1
{W0}: random intercepts
{W1}: random coefficients for covariate2
------------------------------------------------------------------------------
Model summary
------------------------------------------------------------------------------
Likelihood:
outcome ~ normal(xb_outcome,{e.outcome:sigma2})
Priors:
{outcome:covariate1 covariate2 _cons} ~ normal(0,10000) (1)
{U0} ~ normal(0,{U:sigma2}) (1)
{V0} ~ normal(0,{V:sigma2}) (1)
{V1} ~ normal(0,{V:sigma2}) (1)
{V2} ~ normal(0,{V:sigma2}) (1)
{V3} ~ normal(0,{V:sigma2}) (1)
{V4} ~ normal(0,{V:sigma2}) (1)
{V5} ~ normal(0,{V:sigma2}) (1)
{V6} ~ normal(0,{V:sigma2}) (1)
{V7} ~ normal(0,{V:sigma2}) (1)
{V8} ~ normal(0,{V:sigma2}) (1)
{V9} ~ normal(0,{V:sigma2}) (1)
{V10} ~ normal(0,{V:sigma2}) (1)
{V11} ~ normal(0,{V:sigma2}) (1)
{V12} ~ normal(0,{V:sigma2}) (1)
{V13} ~ normal(0,{V:sigma2}) (1)
{V14} ~ normal(0,{V:sigma2}) (1)
{V15} ~ normal(0,{V:sigma2}) (1)
{V16} ~ normal(0,{V:sigma2}) (1)
{V17} ~ normal(0,{V:sigma2}) (1)
{V18} ~ normal(0,{V:sigma2}) (1)
{V19} ~ normal(0,{V19:sigma2}) (1)
{W0} ~ normal(0,{W0:sigma2}) (1)
{W1} ~ normal(0,{W1:sigma2}) (1)
{e.outcome:sigma2} ~ igamma(.01,.01)
Hyperpriors:
{U:sigma2} ~ igamma(.01,.01)
{V:sigma2} ~ igamma(.01,.01)
{V19:sigma2} ~ igamma(.01,.01)
{W0:sigma2} ~ igamma(.01,.01)
{W1:sigma2} ~ igamma(.01,.01)
------------------------------------------------------------------------------
(1) Parameters are elements of the linear form xb_outcome.
Bayesian multilevel regression MCMC iterations = 210,000
Metropolis–Hastings and Gibbs sampling Burn-in = 200,000
MCMC sample size = 10,000
-------------------------------------------------------------
| No. of Observations per group
Group variable | groups Minimum Average Maximum
----------------+--------------------------------------------
_all | 1 97 97.0 97
Group1 | 1 97 97.0 97
-------------------------------------------------------------
Number of obs = 97
Acceptance rate = .5045
Efficiency: min = .002255
avg = .06496
Log marginal-likelihood max = .2841
------------------------------------------------------------------------------
| Equal-tailed
| Mean Std. dev. MCSE Median [95% cred. interval]
-------------+----------------------------------------------------------------
outcome |
covariate1 | -.0110874 .0059492 .000344 -.0109953 -.0228833 .0004529
covariate2 | -.1876556 .1218701 .006668 -.1864473 -.4312361 .05925
_cons | .6759434 .1800871 .037924 .6768317 .3181374 1.030745
-------------+----------------------------------------------------------------
Group2 |
U:sigma2 | .0267817 .0142745 .001079 .0236498 .0091008 .0632229
-------------+----------------------------------------------------------------
_all |
V:sigma2 | .0328675 .0291046 .003738 .0227667 .0049685 .1147533
V19:sigma2 | 648.4501 21147.34 396.721 .0976229 .0059268 45.86536
-------------+----------------------------------------------------------------
Group1 |
W0:sigma2 | .0076228 .005156 .000306 .0063366 .0020867 .0204744
W1:sigma2 | .0304866 .0299592 .002302 .020928 .0038986 .1133497
-------------+----------------------------------------------------------------
e.outcome |
sigma2 | .019182 .0033781 .000083 .0188363 .0135053 .0268277
------------------------------------------------------------------------------
Note: Default priors are used for model parameters.
Note: There is a high autocorrelation after 500 lags.
burnin(20000), burnin(100000) and burnin(150000).
For all three sizes, everything works fine.
Has anyone else encountered similar problems and successfully fixed them?
Just in case: at the moment, I’m using Nikolay’s workaround to retrieve estimation results (described here Retrieving estimation results from "bayes [...] saving(xyz)" command):
Code:
.
. *Make temporary directory to store simulation dataset
. capture mkdir "temp_sim"
.
. copy "investigate_bayes_mixed_simulation.dta" "temp_sim/investigate_bayes_mixed_simulation.dta"
.
. *Re-run Bayesian model with burnin(1)
. use test_data, clear
.
.
. qui tab Group2, gen(id_Group2)
. unab idvar: id_Group2*
. foreach v of local idvar {
2. gen inter`v' =covariate2*`v'
3. }
.
.
. *--------------------------------------------------*
. bayes, burnin(1) rseed(42) saving(investigate_bayes_mixed_simulation, replace): ///
> mixed outcome covariate1 covariate2 || _all: R.Group2 || _all: inter*, cov(identity) nocons || Group1: covariate2
note: Gibbs sampling is used for regression coefficients and some variance components.
Burn-in 1 done
Simulation 10000 aaaaaaaaa1000aaaaaaaaa2000aaaaa....3000.........4000.........5000.........6000.........7000.........8000.........9000.........10000 done
file investigate_bayes_mixed_simulation.dta saved.
Multilevel structure
------------------------------------------------------------------------------
Group2
{U0}: random intercepts
_all
{V0}: random intercepts
{V1}: random coefficients for interid_Group21
{V2}: random coefficients for interid_Group22
{V3}: random coefficients for interid_Group23
{V4}: random coefficients for interid_Group24
{V5}: random coefficients for interid_Group25
{V6}: random coefficients for interid_Group26
{V7}: random coefficients for interid_Group27
{V8}: random coefficients for interid_Group28
{V9}: random coefficients for interid_Group29
{V10}: random coefficients for interid_Group210
{V11}: random coefficients for interid_Group211
{V12}: random coefficients for interid_Group212
{V13}: random coefficients for interid_Group213
{V14}: random coefficients for interid_Group214
{V15}: random coefficients for interid_Group215
{V16}: random coefficients for interid_Group216
{V17}: random coefficients for interid_Group217
{V18}: random coefficients for interid_Group218
{V19}: random coefficients for interid_Group219
Group1
{W0}: random intercepts
{W1}: random coefficients for covariate2
------------------------------------------------------------------------------
Model summary
------------------------------------------------------------------------------
Likelihood:
outcome ~ normal(xb_outcome,{e.outcome:sigma2})
Priors:
{outcome:covariate1 covariate2 _cons} ~ normal(0,10000) (1)
{U0} ~ normal(0,{U:sigma2}) (1)
{V0} ~ normal(0,{V:sigma2}) (1)
{V1} ~ normal(0,{V:sigma2}) (1)
{V2} ~ normal(0,{V:sigma2}) (1)
{V3} ~ normal(0,{V:sigma2}) (1)
{V4} ~ normal(0,{V:sigma2}) (1)
{V5} ~ normal(0,{V:sigma2}) (1)
{V6} ~ normal(0,{V:sigma2}) (1)
{V7} ~ normal(0,{V:sigma2}) (1)
{V8} ~ normal(0,{V:sigma2}) (1)
{V9} ~ normal(0,{V:sigma2}) (1)
{V10} ~ normal(0,{V:sigma2}) (1)
{V11} ~ normal(0,{V:sigma2}) (1)
{V12} ~ normal(0,{V:sigma2}) (1)
{V13} ~ normal(0,{V:sigma2}) (1)
{V14} ~ normal(0,{V:sigma2}) (1)
{V15} ~ normal(0,{V:sigma2}) (1)
{V16} ~ normal(0,{V:sigma2}) (1)
{V17} ~ normal(0,{V:sigma2}) (1)
{V18} ~ normal(0,{V:sigma2}) (1)
{V19} ~ normal(0,{V19:sigma2}) (1)
{W0} ~ normal(0,{W0:sigma2}) (1)
{W1} ~ normal(0,{W1:sigma2}) (1)
{e.outcome:sigma2} ~ igamma(.01,.01)
Hyperpriors:
{U:sigma2} ~ igamma(.01,.01)
{V:sigma2} ~ igamma(.01,.01)
{V19:sigma2} ~ igamma(.01,.01)
{W0:sigma2} ~ igamma(.01,.01)
{W1:sigma2} ~ igamma(.01,.01)
------------------------------------------------------------------------------
(1) Parameters are elements of the linear form xb_outcome.
Bayesian multilevel regression MCMC iterations = 10,001
Metropolis–Hastings and Gibbs sampling Burn-in = 1
MCMC sample size = 10,000
-------------------------------------------------------------
| No. of Observations per group
Group variable | groups Minimum Average Maximum
----------------+--------------------------------------------
_all | 1 97 97.0 97
Group1 | 1 97 97.0 97
-------------------------------------------------------------
Number of obs = 97
Acceptance rate = .5247
Efficiency: min = .001728
avg = .1486
Log marginal-likelihood max = 1
------------------------------------------------------------------------------
| Equal-tailed
| Mean Std. dev. MCSE Median [95% cred. interval]
-------------+----------------------------------------------------------------
outcome |
covariate1 | -.010756 .0057898 .000262 -.0107464 -.022353 .0008853
covariate2 | -.1865147 .1113828 .005497 -.1863413 -.3999516 .0360584
_cons | .589891 .2360603 .05678 .6385582 -.1050936 .9248694
-------------+----------------------------------------------------------------
Group2 |
U:sigma2 | .0252251 .0130876 .000589 .0224263 .0083923 .0581508
-------------+----------------------------------------------------------------
_all |
V:sigma2 | .0413304 .0492257 .00763 .0251147 .0051311 .1700203
V19:sigma2 | 363.2888 14164.7 141.647 .1094709 .0061633 61.7817
-------------+----------------------------------------------------------------
Group1 |
W0:sigma2 | .0078951 .0050638 .000273 .0065598 .0022338 .0212309
W1:sigma2 | .0275205 .027656 .002792 .0185549 .0036315 .1026532
-------------+----------------------------------------------------------------
e.outcome |
sigma2 | .0192395 .0035545 .000092 .0187985 .0135397 .0274202
------------------------------------------------------------------------------
Note: Default priors are used for model parameters.
Note: Adaptation continues during simulation.
Note: There is a high autocorrelation after 500 lags.
Note: Adaptation tolerance is not met in at least one of the blocks.
. estimates save investigate_bayes_mixed, replace
file investigate_bayes_mixed.ster saved
.
. *--------------------------------------------------*
.
. *Copy simulation dataset back to current directory
. copy "temp_sim/investigate_bayes_mixed_simulation.dta" "investigate_bayes_mixed_simulation.dta", replace
.
. *Re-run model based on initial run
.
. bayes
Multilevel structure
------------------------------------------------------------------------------
Group2
{U0}: random intercepts
_all
{V0}: random intercepts
{V1}: random coefficients for interid_Group21
{V2}: random coefficients for interid_Group22
{V3}: random coefficients for interid_Group23
{V4}: random coefficients for interid_Group24
{V5}: random coefficients for interid_Group25
{V6}: random coefficients for interid_Group26
{V7}: random coefficients for interid_Group27
{V8}: random coefficients for interid_Group28
{V9}: random coefficients for interid_Group29
{V10}: random coefficients for interid_Group210
{V11}: random coefficients for interid_Group211
{V12}: random coefficients for interid_Group212
{V13}: random coefficients for interid_Group213
{V14}: random coefficients for interid_Group214
{V15}: random coefficients for interid_Group215
{V16}: random coefficients for interid_Group216
{V17}: random coefficients for interid_Group217
{V18}: random coefficients for interid_Group218
{V19}: random coefficients for interid_Group219
Group1
{W0}: random intercepts
{W1}: random coefficients for covariate2
------------------------------------------------------------------------------
Model summary
------------------------------------------------------------------------------
Likelihood:
outcome ~ normal(xb_outcome,{e.outcome:sigma2})
Priors:
{outcome:covariate1 covariate2 _cons} ~ normal(0,10000) (1)
{U0} ~ normal(0,{U:sigma2}) (1)
{V0} ~ normal(0,{V:sigma2}) (1)
{V1} ~ normal(0,{V:sigma2}) (1)
{V2} ~ normal(0,{V:sigma2}) (1)
{V3} ~ normal(0,{V:sigma2}) (1)
{V4} ~ normal(0,{V:sigma2}) (1)
{V5} ~ normal(0,{V:sigma2}) (1)
{V6} ~ normal(0,{V:sigma2}) (1)
{V7} ~ normal(0,{V:sigma2}) (1)
{V8} ~ normal(0,{V:sigma2}) (1)
{V9} ~ normal(0,{V:sigma2}) (1)
{V10} ~ normal(0,{V:sigma2}) (1)
{V11} ~ normal(0,{V:sigma2}) (1)
{V12} ~ normal(0,{V:sigma2}) (1)
{V13} ~ normal(0,{V:sigma2}) (1)
{V14} ~ normal(0,{V:sigma2}) (1)
{V15} ~ normal(0,{V:sigma2}) (1)
{V16} ~ normal(0,{V:sigma2}) (1)
{V17} ~ normal(0,{V:sigma2}) (1)
{V18} ~ normal(0,{V:sigma2}) (1)
{V19} ~ normal(0,{V19:sigma2}) (1)
{W0} ~ normal(0,{W0:sigma2}) (1)
{W1} ~ normal(0,{W1:sigma2}) (1)
{e.outcome:sigma2} ~ igamma(.01,.01)
Hyperpriors:
{U:sigma2} ~ igamma(.01,.01)
{V:sigma2} ~ igamma(.01,.01)
{V19:sigma2} ~ igamma(.01,.01)
{W0:sigma2} ~ igamma(.01,.01)
{W1:sigma2} ~ igamma(.01,.01)
------------------------------------------------------------------------------
(1) Parameters are elements of the linear form xb_outcome.
Bayesian multilevel regression MCMC iterations = 10,001
Metropolis–Hastings and Gibbs sampling Burn-in = 1
MCMC sample size = 10,000
-------------------------------------------------------------
| No. of Observations per group
Group variable | groups Minimum Average Maximum
----------------+--------------------------------------------
_all | 1 97 97.0 97
Group1 | 1 97 97.0 97
-------------------------------------------------------------
Number of obs = 97
Acceptance rate = .5247
Efficiency: min = .002255
avg = .06496
Log marginal-likelihood max = .2841
------------------------------------------------------------------------------
| Equal-tailed
| Mean Std. dev. MCSE Median [95% cred. interval]
-------------+----------------------------------------------------------------
outcome |
covariate1 | -.0110874 .0059492 .000344 -.0109953 -.0228833 .0004529
covariate2 | -.1876556 .1218701 .006668 -.1864473 -.4312361 .05925
_cons | .6759434 .1800871 .037924 .6768317 .3181374 1.030745
-------------+----------------------------------------------------------------
Group2 |
U:sigma2 | .0267817 .0142745 .001079 .0236498 .0091008 .0632229
-------------+----------------------------------------------------------------
_all |
V:sigma2 | .0328675 .0291046 .003738 .0227667 .0049685 .1147533
V19:sigma2 | 648.4501 21147.34 396.721 .0976229 .0059268 45.86536
-------------+----------------------------------------------------------------
Group1 |
W0:sigma2 | .0076228 .005156 .000306 .0063366 .0020867 .0204744
W1:sigma2 | .0304866 .0299592 .002302 .020928 .0038986 .1133497
-------------+----------------------------------------------------------------
e.outcome |
sigma2 | .019182 .0033781 .000083 .0188363 .0135053 .0268277
------------------------------------------------------------------------------
Note: Default priors are used for model parameters.
Note: There is a high autocorrelation after 500 lags.
.
. log close
Alex
Comment