Parameter significance of bayesian logistic

[ece bacaksiz]

Hi all
How can i compute parameter significances after bayes logit? In search button, i see ''bic'' recommended, but i couldnt make it work.
Thanks for help.

[Matt Warkentin]

Hi Ece,

Formal hypothesis testing is not all that common in Bayesian analysis. It is more of a frequentist approach, so if you're looking for an analogous approach you could evaluate the posterior density for any specific parameter and assess the probability that a parameter is in a certain range. If you look at the equi-tailed or high-density posterior intervals provided by Stata you can see whether it contains a null value, which can help provide some insight into the "significance" of the parameter. Or you can look at the probability that the posterior effect is greater than/less than a value.

Sometimes Bayes factors are used to compare models with and without a parameter included to test which model fits best.

Try the following code to investigate whether drug==2 has an effect on mortality:

sysuse cancer

bayesmh died i.drug , likelihood(logit) prior( {died:}, normal(0,10))
bayesstats summary, hpd

bayesgraph diag {died: 2.drug}

bayestest interval {died: 2.drug} , upper(0)

From the above code we see that the probability the effect for drug==2 is less than zero (i.e. reduces mortality) is extremely large, so we may be inclined to think this is a meaningful association.

Next we can compute the Bayes factors for the model including and excluding the entire drug variable (just as an example).

bayesmh died i.drug , likelihood(logit) prior( {died:}, normal(0,10)) saving(full)
estimates store full

bayesmh died , likelihood(logit) prior( {died:}, normal(0,10)) saving(intercept)
estimates store intercept

bayestest model full intercept

bayestest model full intercept, prior(0.8, 0.2)

Using the results from the first --bayestest-- code we can compute the Bayes factor as simply the ratio of the P(M | y) and we see that the full model including the drug variable is approx. 80 times more probable. This was when we gave equal probability to the two models. Let's say we think that the model with the drug variable is more likely based on our expert knowledge, we can change the prior to be 80% for the full model and 20% for the intercept only model. The new Bayes factor is 1286, so the full model is MUCH more likely (BF = 0.9969/0.0031 * 0.8/0.2).

Hope this was informative.

[ece bacaksiz]

Hi Matt
Thanks for interest. Explanation was great, I used the command you advised; fitted it for multivariable model, for 3 independents, then 2 independents, and 1, last. I get the log(ML) and P(M|y) calculations. The smaller ML is better fit, right ?(for the unrestricted model).Then we decide the include the variable. But this is decision only stands for including -or not- a variable. I understand there is no probability calculation?

[Matt Warkentin]

Hi Ece,

You'll have to clarify some of your question. I am unsure what you are asking in the last two lines.

Yes, a smaller log-likelihood is considered better but will improve with more covariates. You'll see in the second column in the table you get the prior for the models P(M) which will be equally distributed among the models if you don't specify your own priors. The third column is the probability of the model given the data and can be used to compute Bayes Factors. More information on BF can be found https://en.wikipedia.org/wiki/Bayes_factor and here https://www.stata.com/manuals/bayesb...esbayesstatsic.

Calculating Bayes Factors can be done automatically by Stata. For example, similar to what you've described...

sysuse cancer

bayesmh died i.drug age, likelihood(logit) prior( {died:}, normal(0,10)) ///
saving(both, replace)
estimates store both

bayesmh died age, likelihood(logit) prior( {died:}, normal(0,10)) ///
saving(age_only, replace)
estimates store age_only

bayesmh died i.drug, likelihood(logit) prior( {died:}, normal(0,10)) ///
saving(drug_only, replace)
estimates store drug_only

bayesmh died , likelihood(logit) prior( {died:}, normal(0,10)) ///
saving(intercept, replace)
estimates store intercept

bayesstats ic intercept drug_only age_only both, basemodel(intercept)
bayestest model intercept drug_only age_only both

We can see that the model containing drug only was the best as it had the largest probability. In comparing each model against the intercept only model we can see that the Bayes Factors for the drug only, age only, and both models are 77, 0.3, and 5, respectively. There are rules of thumb available to interpret the strength of these Bayes factors. These rules of thumb can be found at both links I provided.

[ece bacaksiz]

Dear Matt
I get this result.. Dependent, A and B are binary. C is continuous. Can you help how i can manage with this..

[Matt Warkentin]

Since each model was given the same prior probability, we can just focus on column 3. The largest probability is for model AB with probability of 0.62. If you think this might be you model of choice, you can compare it to the other models and see if the evidence for this model being best is supported.

Model AB vs. A --> 0.62 / 0.22 = 2.8

Model AB vs. AC --> 0.62 / 0.06 = 10.3

Model AB vs. Full --> 0.62 / 0.11 = 5.6

These are your Bayes factors. Using the rules of thumb in the Wikipedia link I sent you, we can say there is 'strong' evidence for AB over AC, 'barely worth a mention' for AB over A, and 'substantial' evidence for AB vs. Full model. You could repeat this for other pairwise model comparisons or even composite comparisons.

[ece bacaksiz]

Matt
Is there a problem in not comparing with the model only-intercept? And do i have to declare this,,,i mean do i have to say ''only-intercept model gave 0.000 value and can't take it base to compare''. I would be glad if i can learn why this happens.
Secondly, how would it be possible to compare probabilities of parameter significance with models predicted by different methods?
Thanks a lot.

[Matt Warkentin]

Hi Ece,

It may not be that the intercept-only model had a probability of exactly 0.000, but it is likely subject to rounding error. If you extract the p-values matrix from the command you might be able to compare to that model -- but based on the table you shared, I was unable to compute that Bayes Factors (or it was equal to infinity). As far as reporting it, well that I can't advise on really except to say that more transparency and disclosure is better when reporting results. Though that being said, comparing to the intercept-only model is the lowest hanging fruit as far as showing whether a variable or set of variables are improving model fit.

Comparing parameter significance across different methods is generally unwise. Do you have a more specific context?

[ece bacaksiz]

Hi Matt
Thanks for useful informations.
I want to compare three different methods, each based on different estimation methods. There is almost no common points in these methods. I think that parameter significance can be comparable in terms of protecting type 1 error,,at least an available benchmarking tool.. Does Stata have this option?

[Joseph Coveney]

I want to compare three different methods, each based on different estimation methods. There is almost no common points in these methods. I think that parameter significance can be comparable in terms of protecting type 1 error,,at least an available benchmarking tool.. Does Stata have this option?

Why, yes, it does!

And perhaps you've had it shown to you before.

Toward that larger question, you might want to take a look at the attached log file, which shows the results of a comparison of test size and power between exlogistic, firthlogit and a Bayesian-like penalization method using log-F(1, 1) priors (as implemented in the user-written penlogit command that is available from the Stata Journal website) for a model of three indicator variable predictors for your sample size of 53 and a baseline response of 50%.

[ece bacaksiz]

Dear Sir

In post 5, in couldnt run the command starts with ''foreach var of new....'' (both in attached files) so i couldnt see the results

i wrote before.
Thanks for help.

[ece bacaksiz]

Error is in

simulate e1 = r(e1) e2 = r(e2) f1 = r(f1) f2 = r(f2) p1 = r(p1) p2 = r(p2),
> ///
> reps(3000) nodots: simem

here. It separates simulate and reps. So i couldn be able to reach your calculation.

[ece bacaksiz]

. simulate e1 = r(e1) e2 = r(e2) f1 = r(f1) f2 = r(f2) p1 = r(p
> 1) p2 = r(p2), reps(3000) nodots: simem
too few ')' or ']'
an error occurred when simulate executed simem

Could you please help. I want to make this work. And manage my study.

[Joseph Coveney]

You don't need to be able to reach my calculation. I attached two files to the post. One is a DO-file, which you may just run as is--you don't need to mess with anything in it or do whatever it is that you're doing that breaks it. The other is the log file that shows the results for that DO-file, so that you don't even need to run the DO-file at all.

If you examine the log file, you will see that all three methods--exact logistic regression, Firth's logistic regression and the Bayesian logistic regression model using the regularizing prior--have very similar frequentist operating characteristics against a dataset that mimics yours. You can choose any of the methods with confidence that none of them will be any more liable to lead you astray when used properly than the other two.

So you don't need to spend endless additional hours in your method-comparison pursuit. Instead, you can focus your efforts on getting up to speed on your Stata coding skills, and on making sure that you get the big picture for your study right.