Six outcomes, six separate models or one?

[Fouziah Almouqati]

Dear all,

I want to calculate the adjusted prediction of the count of CT scans for each year from 2015 to 2022, per the anatomical area of the scan. The outcomes I am analysing are total_headCT, total_neckCT, total_chestCT, total_abdopelCT, total_spineCT, and total_otherCT. My initial plan is to conduct six separate negative binomial (or Poisson) models for each outcome, estimate the margins, and then combine these into a single marginsplot.

My question is: Is this approach acceptable? I reviewed similar literature where they used negative binomial models for count data, but there were no details on whether the analyses were conducted in one combined model or in separate models. Additionally, if I run six separate models, how can I create a margins plot ( a line ) that includes all the trends in one graph?

Description of the data:
Observations are per individual patient. For each patient, I have their sociodemographic characteristics, diagnosis, arrival mode, and other factors that I wish to adjust for in order to calculate the count for each year.

One patient may experience more than one of the outcomes. For example, they might have a head CT, a neck CT and a spine CT.

---

Below are examples of two codes I used to model head and neck CT scans separately:

nbreg total_headCT ib(first).pre_year_cat ib(first).age_grp ib(first).p_sex ib(first).ind_stat3 ib(first).claim_grp ib(first).triage_code ib(first).PD_grp4 ib(first).treat_clin_exp2 ib(first).rf_source4 ib(first).arri_means_grp4 ib(first).pre_shift ib(first).pre_day_type , dispersion(mean) irr allbaselevels
margins pre_year_cat

nbreg total_neckCT ib(first).pre_year_cat ib(first).age_grp ib(first).p_sex ib(first).ind_stat3 ib(first).claim_grp ib(first).triage_code ib(first).PD_grp4 ib(first).treat_clin_exp2 ib(first).rf_source4 ib(first).arri_means_grp4 ib(first).pre_shift ib(first).pre_day_type , dispersion(mean) irr allbaselevels
margins pre_year_cat

---

I hope my explanation is clear, and I would greatly appreciate any advice, help, or replies.

[Jeff Wooldridge]

There are many others who can provide better guidance about the graphing question. As far as the estimation, it's fine to estimate a separate model for each outcome. Moreover, because it looks like the explanatory variables are the same in each equation, there would be little or no efficiency gain to estimating the equations jointly. (In the linear case, there would be no efficiency gain using OLS on each equation or GLS.)

I assume there is an entry for each patient, with lots of zeros for some of the outcomes. This doesn't change anything. The Poisson regression estimator is, technically, more robust for estimating the conditional mean.

You want to plot the predicted values against one of the variables in the model?

[Clyde Schechter]

if I run six separate models, how can I create a margins plot ( a line ) that includes all the trends in one graph?

I don't think you will be able to do this using -marginsplot-. You can run -margins- following each regression and save the results with the (undocumented) -generate()- option, creating a new set of variables (one for each value of pre_year_cat) for each of your outcomes. I would then move those new variables into a new frame, reshape the data long, and draw the graphs. You don't give example data, but here's an illustration of the overall approach using the auto.dta. You can adapt it to your situation.

clear*
sysuse auto, clear

foreach v of varlist mpg headroom trunk {
    regress `v' i.rep78
    margins rep78,  generate(`v'_margin)
}

frame put *margin*, into(for_graphs)
frame change for_graphs
gen `c(obs_t)' obs_no = _n
rename *_margin# *#
reshape long mpg headroom trunk, i(obs_no) j(rep78)

graph twoway line mpg headroom trunk rep78, sort

Added: Sorry, the above approach is too simple and will not work in your model. It would work correctly if you had no covariates in the mode. Because of the covariates, you need to one additional command: - collapse (mean) mpg headroom trunk), by(rep78)- between the -reshape- and the -graph-. (As a practical matter, after the -collapse- you may also want to relabel the variables so that the graph legend isn't cluttered up with "(mean)".