Margins after menbreg with exposure

[Caroline Thirukumaran]

Hello,

I am estimating the following multilevel negative binomial model with an exposure term.

 menbreg surgery i.year##c.prop_race $covars, exposure(population) || countynum:, irr

In this model, surgery represents the number of surgeries, year is the year in which they were performed (2013-2016), prop_race is the proportion of a particular race in the county, population is the total population of the county, countynum is the county ID, and $covars is a set of county covariates.

1. After estimating this model, I would like to obtain the number of surgeries or their rates for counties at, say, the 25th percentile of prop_race for each year. How could I program this such that the code takes into account the exposure variable. Will I be able to interpret the output in counts of surgeries or as a rates? While I am not certain, I am thinking that the following code may help me get to some of it

 margins i.year, at(prop_race=(0.06))

if 0.06 is the value of prop_race at the 25th percentile.

2. Is there a way I could use the -predict- command to obtain numbers similar to what I will obtain from margins.

I am using Stata/MP 15.1. I have also looked through some previous discussions which examine margins after negative binomial models with exposure (not multilevel though), or after multilevel negative binomial models (but without the exposure term).

Thank you,
Caroline

[Clyde Schechter]

It looks to me as if you can get exactly what you are asking for from your -margins- command. You say it may help you get some of it: what does it not provide that you are looking for. By the way, you don't have to pre-calculate the 25th percentile yourself. -margins- can do it for you:

margins i.year, at((p25) prop_race)

The outputs you get will be expected numbers of surgeries, all adjusted to the observed sample distribution of the covariates. And the exposure will be automatically included in the calculation--you don't have to do anything to make that happen.

As for your second question, it depends on what you mean. In one sense, the answer is clearly no: -predict- gives you observation by observation modeled results, whereas the output of -margins- is a summary statistic for the entire sample and it estimates a population-level parameter. In a different sense, -margins- is actually a wrapper program for -predict-. It uses -predict- to calculate individual predictions and then generates its summary statistics from those.

[Caroline Thirukumaran]

Thank you for your response, Dr. Schecter. These points are very helpful.

- Following-up on #1: If I needed to present the adjusted rates (i.e. surgeries *1000 /population) instead of numbers of surgeries - are there options that I could add to the margins command that would help me get the rates?

- Following-up on #2: If I were to use the -predict- command and obtain adjusted numbers per-county per-year, which commands would help me get the marginal estimates that the -margins- command gives me?

Best,
Caroline

[Clyde Schechter]

#1.

margins i.year, at((p25) prop_race) expression(1000*predict(mu)/population)

#2. The gist of the code would be this:

levelsof year, local(years)
centile prop_race, centile(25)
local set_prop_race `r(c_1)'
preserve
replace prop_race = `set_prop_race'
foreach y of local years {
    replace year = `y'
    predict result, mu
    replace result = 1000*result/population
    summ result, meanonly
    display "Adjusted rate for `y': `r(mean)'"
    drop result
}
restore

I have not tested this code, so it may contain errors. But this is the gist of it. You loop over all the values of year, and in each iteration you replace year by the current index year, and prop_race is always replaced by the 25th percentile. You run -predict- for mu, then multiply by 1000 and divide b population. Summarize that and output the mean.

Note that this does not give you the standard errors: just the adjusted rates themselves. I do not know how to code for the standard errors: that is more complicated.

[Caroline Thirukumaran]

Dr. Schechter - many thanks for your response! These are so helpful! I will run your code for # 2 and will work through any errors if they arise. Thanks again!

Best,
Caroline