Page 1 of 1

95%Credible intervals are smaller than OR

Posted: Sat Mar 20, 2021 6:10 pm
by ManuelDewez
Hello,

I am running the following multilevel model with a binary outcome with mcmc:
runmlwin unique_dipstick_avai2cats cons private secondaryhosp paed_mother_child centered_unique_turn19_time centered_exp_capita gov VHIOOP, level2(country2: cons) level1(id:) discrete(distribution(binomial) link(logit) denom(cons) pql2) nopause
runmlwin unique_dipstick_avai2cats cons private secondaryhosp paed_mother_child centered_unique_turn19_time centered_exp_capita gov VHIOOP, level2(country2: cons, residuals(u)) level1(id:) discrete(distribution(binomial) link(logit) denom(cons) pql2) mcmc(burnin(1000) chain(550000)) initsprevious nopause nogroup
I have 504 observations (level 1) and 29 second level units.

My outcome is whether urine dipstick are available (0=no; 1=yes). I am using mcmc because the response proportion is extreme (0=42/504; 1=462/504).

I checked the diagnostics and following the Raftery-Lewis and Brooks-Drappers diagnostics I re run the model with 550,00 iterations.

Some of the odds ratio and Standard deviations are huge for 2 of the variables (OR for gov:1156.948 , and VHIOOP:13730.1), and the 95%Credible interval is smaller than the parameter for VHIOOP: 95%CI = 0.050643 - 164.6105

Burnin = 1000
Chain = 550000
Thinning = 1
Run time (seconds) = 156
Deviance (dbar) = 214.32
Deviance (thetabar) = 195.41
Effective no. of pars (pd) = 18.91
Bayesian DIC = 233.23

Disptick_available OR Std. Dev ESS P [95% Cred.Int]

gov 1156.948 24956.07 8874 0.003 2.891191 2926.618
VHIOOP 13730.1 1120034 14575 0.307 0.050643 164.6105

Can anyone explain why is this happening, and what should I do to get correct OR and 95%CI?

Thank you!

Manuel

Re: 95%Credible intervals are smaller than OR

Posted: Mon Mar 22, 2021 2:48 pm
by ChrisCharlton
That does sound a bit strange. Have you tried plotting the chains? What do the results look like it you don't convert them to odds-ratios?

Re: 95%Credible intervals are smaller than OR

Posted: Tue Mar 23, 2021 9:32 am
by ManuelDewez
Hello Chris,

The ACF plot shows quite some autocorrelation and the MCSE plot show that I should run the model with more iterations ( which I did)
The non-ORs results look fine ( the parameters are within the 95% Credible Intervals).

Gov and VHIOOP are 2 categories of the same variable. I think the issue is that there are too few observations for this variable (only 8 for gov and 5 for VHIOOP:


unique_dip |
stick_avai | source_financing2
2cats | compshi gov vhi_oop | Total
-----------+---------------------------------+----------
0 | 29 8 5 | 42
1 | 221 211 30 | 462
-----------+---------------------------------+----------
Total | 250 219 35 | 504

And this is creating the troubles ?
I think I will just ignore this variable for this particular model (I don't have this issue with other models where there is a more balanced distribution of the observations across the outcome's binary categories). What do you think?

Thank you.

Manuel

Re: 95%Credible intervals are smaller than OR

Posted: Tue Mar 23, 2021 1:34 pm
by ChrisCharlton
The way that the credible intervals for this would be calculated is as follows:
  1. Fetch the MCMC chain for the parameter
  2. Exponentiate the values to convert to odds-ratios
  3. Sort the chain
  4. Extract the 2.5% and 97.5% quantiles
I suppose that if the chain contained a few extreme values then this could result in the mean being outside the range of the credible intervals. You might also want to look at the mode and median values for comparison. For a more visual examination of the chain distribution you can use the fiveway option in the mcmcsum command to look at the kernel density.

Re: 95%Credible intervals are smaller than OR

Posted: Tue Mar 23, 2021 2:15 pm
by ManuelDewez
Thank you Chris,

this is the syntax I am using :
runmlwin unique_dipstick_avai2cats cons private secondaryhosp paed_mother_child centered_unique_turn19_time centered_exp_capita gov VHIOOP if unique_hospital==1, level2(country2: cons) level1(id:) discrete(distribution(binomial) link(logit) denom(cons) pql2) nopause
followed by :
runmlwin unique_dipstick_avai2cats cons private secondaryhosp paed_mother_child centered_unique_turn19_time centered_exp_capita gov VHIOOP if unique_hospital==1, level2(country2: cons, residuals(u)) level1(id:) discrete(distribution(binomial) link(logit) denom(cons) pql2) mcmc(burnin(1000) chain(550000)) initsprevious nopause nogroup
to get OR I use:
runmlwin, or

this provides ORs and 95 %Cred. Int.


How do I sort and extract the 2.5 and 97.5 % quantiles?
Thank you

Re: 95%Credible intervals are smaller than OR

Posted: Tue Mar 23, 2021 2:58 pm
by ChrisCharlton
You can obtain the MCMC chains as a dataset with the following command:

Code: Select all

mcmcsum, getchains
Once you have this you can use the standard Stata data manipulation commands to process it, i.e. gen, sort and centile .

You ought to get the same results as before as this should be equivalent to what runmlwin is doing itself.

Re: 95%Credible intervals are smaller than OR

Posted: Tue Mar 23, 2021 3:56 pm
by ManuelDewez
That worked well, thank you Chris!