Hello, I'm running a hierarchical Bayesian regression. If I set off-the-shelf prior distributions like this, I have no problem:
bayesmh outcome treatment covariate1 covariate2 covariate3 covariate4 covariate5 U[clusterID],
likelihood(normal({var_resid}))
block({var_resid} {var_U}, split)
burnin(20000) mcmcsize(5000) thin(2) nchains(2) rseed(712) dots
prior({outcome:treatment}, normal(8.0, 6.7))
prior({outcome:_cons}, normal(0, 0.4))
prior({outcome:covariate1}, normal(0, 0.4))
prior({outcome:covariate2}, normal(0, 0.4))
prior({outcome:covariate3}, normal(0, 0.4))
prior({outcome:covariate4}, normal(0, 0.4))
prior({outcome:covariate5}, normal(0, 0.4))
prior({U}, normal(0,{var_U}))
prior({var_U}, density(2*({var_U}>0)*normalden(sqrt({var_U}), 0, .5)))
prior({var_resid}, igamma(0.01,0.01))
initall( {outcome:_cons} rnormal(0,0.1)
{outcome:treatment} rnormal(0,0.2)
{outcome:covariate1} rnormal(0,0.3)
{outcome:covariate2} rnormal(0,0.3)
{outcome:covariate3} rnormal(0,0.3)
{outcome:covariate4} rnormal(0,0.3)
{outcome:covariate5} rnormal(0,0.3)
{var_U} runiform(0.01, 0.02)
{var_resid} runiform(0.01, 0.02) );
Note that {var_U} gets a bespoke prior density (half-normal on the SD scale), no problem.
But if I add a density suboption to {outcome:treatment}, even if it is just the same as above, I get an error:
bayesmh outcome treatment covariate1 covariate2 covariate3 covariate4 covariate5 U[clusterID],
likelihood(normal({var_resid}))
block({var_resid} {var_U}, split)
burnin(20000) mcmcsize(5000) thin(2) nchains(2) rseed(712) dots
prior({outcome:treatment}, density(normalden({outcome:treatment}, 8.0, 6.7))) <--- this is the line that causes the trouble
prior({outcome:_cons}, normal(0, 0.4))
prior({outcome:covariate1}, normal(0, 0.4))
prior({outcome:covariate2}, normal(0, 0.4))
prior({outcome:covariate3}, normal(0, 0.4))
prior({outcome:covariate4}, normal(0, 0.4))
prior({outcome:covariate5}, normal(0, 0.4))
prior({U}, normal(0,{var_U}))
prior({var_U}, density(2*(sqrt({var_U})>0)*normalden(sqrt({var_U} ), 0, .5)))
prior({var_resid}, igamma(0.01,0.01))
initall( {outcome:_cons} rnormal(0,0.1)
{outcome:treatment} rnormal(0,0.2)
{outcome:covariate1} rnormal(0,0.3)
{outcome:covariate2} rnormal(0,0.3)
{outcome:covariate3} rnormal(0,0.3)
{outcome:covariate4} rnormal(0,0.3)
{outcome:covariate5} rnormal(0,0.3)
{var_U} runiform(0.01, 0.02)
{var_resid} runiform(0.01, 0.02) );
The error reads:
attempting to create a group named outcome but a linear equation exists with that name
r(110);
Now, I think this stems from me naming the parameter {outcome:treatment} inside the prior(...density()) suboption. I read in [BAYES] v18.0 p.153 about using the xb (sub)suboption to distinguish between {outcome:treatment} and ({outcome:_cons}+{outcome:treatment}*treatment), unless the linear combination has been previously defined. But surely it has been previously defined, because it is in the main body of the command:
bayesmh outcome treatment covariate1 covariate2 covariate3 covariate4 covariate5 U[clusterID]
Unless Stata encounters the prior definition first, defines:
{outcome:treatment} := ({outcome:_cons}+{outcome:treatment}*treatment)
and then does not know what to do with covariate1-covariate5. I can't think of any other reason why it would try to make a {outcome:} group twice. But it seems implausible that the prior options would be processed before the linear predictor expression.
Using -set trace on- does not help as there are so many commands and subcommands firing off that it is not clear to a mere mortal what generated the error. stset (of all things!!!) is the last command to appear in the trace.
Can anyone advise on how I can set a bespoke prior density for a coefficient? (I don't actually want to set a normal prior, I have a formula which is not a Stata preset but it would be too complicated and distracting to show it here). I wonder whether there is some trick about define(), which used to be redefine() and could solve some issues with random effects naming.
Many thanks
Robert
bayesmh outcome treatment covariate1 covariate2 covariate3 covariate4 covariate5 U[clusterID],
likelihood(normal({var_resid}))
block({var_resid} {var_U}, split)
burnin(20000) mcmcsize(5000) thin(2) nchains(2) rseed(712) dots
prior({outcome:treatment}, normal(8.0, 6.7))
prior({outcome:_cons}, normal(0, 0.4))
prior({outcome:covariate1}, normal(0, 0.4))
prior({outcome:covariate2}, normal(0, 0.4))
prior({outcome:covariate3}, normal(0, 0.4))
prior({outcome:covariate4}, normal(0, 0.4))
prior({outcome:covariate5}, normal(0, 0.4))
prior({U}, normal(0,{var_U}))
prior({var_U}, density(2*({var_U}>0)*normalden(sqrt({var_U}), 0, .5)))
prior({var_resid}, igamma(0.01,0.01))
initall( {outcome:_cons} rnormal(0,0.1)
{outcome:treatment} rnormal(0,0.2)
{outcome:covariate1} rnormal(0,0.3)
{outcome:covariate2} rnormal(0,0.3)
{outcome:covariate3} rnormal(0,0.3)
{outcome:covariate4} rnormal(0,0.3)
{outcome:covariate5} rnormal(0,0.3)
{var_U} runiform(0.01, 0.02)
{var_resid} runiform(0.01, 0.02) );
Note that {var_U} gets a bespoke prior density (half-normal on the SD scale), no problem.
But if I add a density suboption to {outcome:treatment}, even if it is just the same as above, I get an error:
bayesmh outcome treatment covariate1 covariate2 covariate3 covariate4 covariate5 U[clusterID],
likelihood(normal({var_resid}))
block({var_resid} {var_U}, split)
burnin(20000) mcmcsize(5000) thin(2) nchains(2) rseed(712) dots
prior({outcome:treatment}, density(normalden({outcome:treatment}, 8.0, 6.7))) <--- this is the line that causes the trouble
prior({outcome:_cons}, normal(0, 0.4))
prior({outcome:covariate1}, normal(0, 0.4))
prior({outcome:covariate2}, normal(0, 0.4))
prior({outcome:covariate3}, normal(0, 0.4))
prior({outcome:covariate4}, normal(0, 0.4))
prior({outcome:covariate5}, normal(0, 0.4))
prior({U}, normal(0,{var_U}))
prior({var_U}, density(2*(sqrt({var_U})>0)*normalden(sqrt({var_U} ), 0, .5)))
prior({var_resid}, igamma(0.01,0.01))
initall( {outcome:_cons} rnormal(0,0.1)
{outcome:treatment} rnormal(0,0.2)
{outcome:covariate1} rnormal(0,0.3)
{outcome:covariate2} rnormal(0,0.3)
{outcome:covariate3} rnormal(0,0.3)
{outcome:covariate4} rnormal(0,0.3)
{outcome:covariate5} rnormal(0,0.3)
{var_U} runiform(0.01, 0.02)
{var_resid} runiform(0.01, 0.02) );
The error reads:
attempting to create a group named outcome but a linear equation exists with that name
r(110);
Now, I think this stems from me naming the parameter {outcome:treatment} inside the prior(...density()) suboption. I read in [BAYES] v18.0 p.153 about using the xb (sub)suboption to distinguish between {outcome:treatment} and ({outcome:_cons}+{outcome:treatment}*treatment), unless the linear combination has been previously defined. But surely it has been previously defined, because it is in the main body of the command:
bayesmh outcome treatment covariate1 covariate2 covariate3 covariate4 covariate5 U[clusterID]
Unless Stata encounters the prior definition first, defines:
{outcome:treatment} := ({outcome:_cons}+{outcome:treatment}*treatment)
and then does not know what to do with covariate1-covariate5. I can't think of any other reason why it would try to make a {outcome:} group twice. But it seems implausible that the prior options would be processed before the linear predictor expression.
Using -set trace on- does not help as there are so many commands and subcommands firing off that it is not clear to a mere mortal what generated the error. stset (of all things!!!) is the last command to appear in the trace.
Can anyone advise on how I can set a bespoke prior density for a coefficient? (I don't actually want to set a normal prior, I have a formula which is not a Stata preset but it would be too complicated and distracting to show it here). I wonder whether there is some trick about define(), which used to be redefine() and could solve some issues with random effects naming.
Many thanks
Robert
Comment