MCMC: calculating estimates, CrIs for combinations of terms/interactions

[amaald21]

Hello,

I'm wondering if there is a way to calculate estimates and credible intervals for different combinations of terms (including interaction terms) if the model was run using runmlwin and MCMC?

I'm looking to create a table like the one in this article (Table 1): https://doi.org/10.1093/ije/dyr218
which would present, for example, odds ratios and credible intervals for the combination of Term A + Term B + Interaction(A+B)

I would typically do this in Stata using the 'lincom' command, but it looks like lincom doesn't work because the estimates/interaction terms are variables, not e-class results, in the saved chains file.

Is there a way of doing this, or a lincom-type command that would work with the saved chains?

Thanks very for your help, much appreciated!
Amanda

[ChrisCharlton]

I think that you should be able to do this with the normal variable manipulation commands such as generate and egen. Once you have the new variable representing the result of the calculation you should be able to get the credible intervals using the summarize command, or the mcmcsum provided with runmlwin.

[amaald21]

Hi Chris,

Thanks so much - that makes perfect sense! I've tried this out using generate and mcmcsum and have got the combinations of terms and their credible intervals.

Is there a way to get the p-value as an r-class result too?

Thank you!
Amanda

[ChrisCharlton]

It should be fairly easy to modify mcmcsum to do this. You can see where the r-class is filled in on lines 535-552 or 410-427 of mcmcsum. To add Bayesian p-values you just need to decide whether you want the mean, median or mode version and add an extra line to return this, i.e.

Code: Select all

return scalar pvalue         = `pvalmean'

If you prefer more conventional p-values you can calculate these from the estimate and variance as usual.

[amaald21]

Thanks Chris. I've tried this and it seems to work, but the p-values are only to one significant figure (they display as '0'). Could you please advise how I might alter the ado so that I can get the p value to a few more significant figures e.g. 0.001?

Thank you!
Amanda

[ChrisCharlton]

The p-values reported here are one-tailed bayesian p-values, and essentially give the proportion of the chain that crosses zero (or one if the variable is exponentiated and it is aware of this). You can find the code for this in runmlwin_mcmcdiag.ado, and is as follows:

Code: Select all

		local compval 0
		if "`eform'" ~= "" {
			local compval 1
		}
		
		if `mean' > `compval' {
			quietly count if `parameter' < `compval'
		}
		else {
			quietly count if `parameter' > `compval'
		}		
		local pvalmean = r(N) / `N'

I suggest testing with a variable that contains both positive and negative numbers to check whether this is why you are only getting values of zero.