Difference (and corresponding delta-based 95% CI) between stored predictions from 2 nbreg models

[Omar Okasha]

Dear Statalist experts,

I think this is sort of basic question (looked it up here before posting but could not find the answer, perhaps I used wrong search terms, if so, my apologies).

Summary of my case: I have stored estimates from two negative binomial regression models for which I would like to estimate the difference and 95% CI using delta method.

Details:
I'm running two negative binomial regression models on a time-series data. The models are adjusted for seasonality (using sin/cos function) and sex. I am trying to predict the outcome (incidence rate for a certain disease) under two different scenarios (out-of-sample prediction?). I managed to store these predictions from the two different models, for certain time points, using prvalue command in a loop. Below are the steps used.

First, running the first model on a specific segment in time series:

nbreg n t sex cos* sin* if tin( ,2009m12), vce(robust) offset(logcohort) irr

n is the count outcome of interest, t is the time since beginning of study, and tin specifies the time segment (from beginning to December 2009) that I am using to get the predicted rates in another time segment that starts from t=91 and ends at t=144 (from January 2011 to the end of the study), as shown below:

 foreach i of num 91/144 {
 qui prvalue, x(t=`i') rest(mean) level(95) delta
 replace est=r(mu) if t==`i'
 replace est_hi=r(mu_hi) if t==`i'
 replace est_lo=r(mu_lo) if t==`i'
}

so now I've stored the predicted rates and corresponding delta-based 95% CI as est, est_lo and est_hi.

Now to the second model (only change is the model running on different time segment, therefore yielding the counterfactual scenario):

est clear
nbreg n t sex cos* sin* if tin(2011m1, ), vce(robust) offset(logcohort) irr

And the predictions from this model:

 foreach i of num 91/144 {
 qui prvalue if tin(2011m1, ), x(t=`i') rest(mean) level(95) delta
 replace exp=r(mu) if t==`i'
 replace exp_hi=r(mu_hi) if t==`i'
 replace exp_lo=r(mu_lo) if t==`i'
}
*

The predicted rates and corresponding delta-based 95% CI are now stored as exp, exp_lo and exp_hi.

To my question: I would like to estimate the difference between est and exp values, and the delta-based 95% CI for this difference. I could not figure out the way to do that.
Your advice is much appreciated.

Omar

[Jun Xu]

It's been quite a while, but I don't think prvalue, as how it was written, is able to accommodate robust standard errors. It may run through though. I would suggest that you simply use bootstrap if it is not very time consuming (less than a day?).

Jun Xu

[Omar Okasha]

Dear Jun,

I appreciate your follow up to this question. I am actually right in the middle of trying to figure our how to do so! this is a bit difficult for me cause I am not a statistician! And I am trying to find relevant references on how exactly this can be done in Stata. Bootstrapping the 95% CI (AKA Interval Estimates ?) was done in a previous study form the US CDC (Moore et al, Lancet 2015), and what they did was running a model on a time period to derive predictions in another time period, so that they can compare the actual data to those predictions (the counterfactual). Those predictions and corresponding 95% CI were based on bootstrapping technique.
I appreciate your help.
Omar

[Omar Okasha]

This is the article I am talking about..