getting estimates when using bayes prefix for melogit

[Madu Abuchi]

Hi Stata forum members,

I need some advice on how to get estimates after fitting melogit using bayesian framework. I have tried using the -parmest- command but I get an error that says "Estimates matrix e(b) must have exactly 1 row"

Below is my example code:

sysuse auto, clear
bayes: melogit foreign trunk || rep78:,
parmest,format(estimate min95 max95 %8.2f p %8.1e) list(,)

can someone help?

Thanks in anticipation.

Madu

Just to add that I use Stata/SE 15.1 and the error number is r(498);

[Joseph Coveney]

I don't know what parmest is, but you could

matrix list e(mean)
matrix list e(cri)
// etc.

to get summary statistics of the posterior distribution.

More generally,

help bayesian_postestimation

[Madu Abuchi]

Thanks Joseph.

Also, I was wondering if you/anyone can advice how I can get the predicted estimate (prevalence) of my outcome, specifically 2.5% percentile and 97.5% percentile of the posterior distribution. I don't seem to have grabbed this from the bayesian postestimation help tab.

Thanks

[Joseph Coveney]

how I can get the predicted estimate (prevalence) of my outcome, specifically 2.5% percentile and 97.5% percentile of the posterior distribution.

I don't quite follow you here.

matrix list e(cri) will get you the 2.5 and 97.5 percentiles of the posterior distribution (well, of the the sample from the posterior that you drew).

If you want a posterior prediction, say, of the outcome, then you'd typically simulate from the posterior.

[Madu Abuchi]

I am working on a multilevel data and my interest is on the posterior prediction after fitting my model.

Please, are you able to advise codewise how I can get around simulating the posterior prediction in Stata. I use Stata 15 and approaching this using the bayes prefix: Below is my sample date/model specification

sysuse auto, clear
bayes: melogit foreign trunk mpg || rep78:,

Thanks

[Joseph Coveney]

I don't think that there is any getting around it, but you can find a worked example here it it helps.

[Joseph Coveney]

According to this you don't go so far as simulating for a binomial / Bernoulli outcome variable. So, for example, if you wanted the predicted probability at the means of the predictors (including random effect) for your particular model, you could do something like the following.

version 15.1

clear *

set seed `=strreverse("1476747")'

quietly sysuse auto

// Fill in missing values of repair record
summarize rep78, meanonly
local first `r(min)'
local last `r(max)'
quietly replace rep78 = runiformint(`first', `last') if mi(rep78)

// Means of predictors

* Indicator variables for random effect
quietly tabulate rep78, generate(rid)
forvalues rep78 = `first'/`last' {
    tempname rid`rep78'
    summarize rid`rep78', meanonly
    scalar define `rid`rep78'' = r(mean)
}

* Other predictors
foreach var of varlist trunk mpg {
    tempname `var'
    summarize `var', meanonly
    scalar define ``var'' = r(mean)
}

bayes, normalprior(`=sqrt(10)') ///
    burnin(20000) adaptation(every(500)) ///
    nodots nomodelsummary nomesummary: ///
        melogit foreign c.(trunk mpg) || rep78:

quietly use "`e(filename)'", clear
quietly expand _frequency
keep eq*

generate double xb = `trunk' * eq3_p1 + `mpg' * eq3_p2 + eq3_p3
forvalues rep78 = `first'/`last' {
    quietly replace xb = xb + `rid`rep78'' * eq1_p`rep78'
}
generate double mu = invlogit(xb)

// Pr(Foreign) at means of predictors
summarize mu, meanonly
display in smcl as text r(mean)
centile mu, centile(2.5 97.5)

exit

[Madu Abuchi]

Hi Joseph,

Thank you for sharing your knowledge and the sample codes.

I think something is missing from my understanding and I am wondering if it is possible to link the predicted posterior distribution of the outcome to the dataset. For example, in the frequentist approach, the post-estimation predictions include the predicted values in the main dataset, but I noticed with the bayesian approach it creates a separate dataset. I have looked at help file to check what is included in this new dataset: It reads

"The saved dataset has the following structure. Variance index records iteration numbers. The bayes prefix saves only states (sets of parameter values) that are different from one iteration to another and the frequency of each state in variable frequency. (Some states may be repeated for discrete parameters.) As such, index may not necessarily contain consecutive integers. Remember to use frequency as a frequency weight if you need to obtain any summaries of this dataset......" This doesn't seem to me as though I can merge this with the original dataset.

So my question is if there's a way to merge this to the main dataset?? Is there any variable worth using for merging? I don't think its the _index variable. The end game of my analysis is to map the predicted posterior in a country using the longitude and latitude coordinates which already exist in the main dataset.

Any thoughts will be appreciated.

Regards

[Joseph Coveney]

This doesn't seem to me as though I can merge this with the original dataset.

So my question is if there's a way to merge this to the main dataset?? Is there any variable worth using for merging? I don't think its the _index variable. The end game of my analysis is to map the predicted posterior in a country using the longitude and latitude coordinates which already exist in the main dataset

You'd cross, not merge the posterior distribution sample with the original dataset and create summary statistics of the posterior predictions.

version 15.1

clear *

set seed `=strreverse("1477241")'

quietly sysuse auto
keep make foreign rep78 trunk mpg

// Fill in missing values of repair record
summarize rep78, meanonly
local first `r(min)'
local last `r(max)'
quietly replace rep78 = runiformint(`first', `last') if mi(rep78)

// Indicator variables for random effect
quietly tabulate rep78, generate(rid)

// Centering predictors
foreach var of varlist trunk mpg {
    summarize `var', meanonly
    generate double c_`var' = `var' - r(mean)
}

bayes, normalprior(`=sqrt(10)') block({foreign:}, split) ///
    prior({U0:sigma2}, jeffreys) ///
    burnin(25000) adaptation(every(500)) ///
    nodots nomodelsummary nomesummary: ///
        melogit foreign c.c_* || rep78:

preserve
quietly use "`e(filename)'", clear
quietly expand _frequency
keep eq*

tempfile tmpfil0
quietly save `tmpfil0'

restore
cross using `tmpfil0'

generate double xb = c_trunk * eq3_p1 + c_mpg * eq3_p2 + eq3_p3
forvalues rep78 = `first'/`last' {
    quietly replace xb = xb + rid`rep78' * eq1_p`rep78'
}
generate double mu = invlogit(xb)

sort make mu
by make: egen double mu_ba = mean(mu)
by make: generate double pct2_5 = mu[250]
by make: generate double pct97_5 = mu[9750]
quietly by make: keep if _n == 1

melogit foreign c.c_* || rep78:, nolrtest nolog
quietly predict double mu_ml, mu

keep make foreign rep78 trunk mpg mu_* pct*

format mu_* pct* %04.2f

sort rep78 foreign make
list rep78 make trunk mpg foreign mu_ml mu_ba pct*, noobs sepby(rep78) abbreviate(20)

exit

[Madu Abuchi]

Thank you heaps. Exactly what I wanted.

[Madu Abuchi]

Hi Joseph,

Assuming one is working with a cluster (say rep78) without a missing value would this step in #9 where the fitted value is replaced be necessary?:

 forvalues rep78 = `first'/`last' {
    quietly replace xb = xb + rid`rep78' * eq1_p`rep78'
}

I noticed the local macro

`first'

and

`last'

were created earlier in the workflow to fill missing values of rep78.

Thank you

[Gleasia Reid]

Hi, so I am new to Bayesian estimation in Stata and I am working on a Bayesian logit model, trying to get predicted probabilities and odds ratios using "bayesstats summary" command but I haven't been successful.