nbreg only converges when starting from poisson estimates

[Molly Jeffery]

I am using nbreg (Stata MP 15.1) to analyze the rate of days on which patients use two different drugs. I don't think the format of the model matters here, but for completeness it is this:

nbreg overlap_days independent_vars if drug1_days>0, exposure(drug1_days) vce(cluster patient_id) difficult

The model does not converge. I allowed it to iterate overnight. The last several hundred iterations gave the same log pseudolikelihood and the message (not concave). Note that the Poisson model and the constant-only model do converge. It's the full model that does not.

However, a Poisson model (same as above, but substitute poisson for nbreg) of the same data does converge (quickly).

Interestingly, if I run the nbreg model starting from the poisson estimates, the nbreg model converges with no problems.

poisson overlap_days independent_vars if drug1_days>0, exposure(drug1_days) vce(cluster patient_id) difficult
matrix b0=e(b)
matrix b1=J(1,126,1)
matrix b1[1,1]=b0
nbreg overlap_days independent_vars if drug1_days>0, exposure(drug1_days) vce(cluster patient_id) difficult from(b1, copy)

Is this a legitimate approach? There are two large-ish coefficients in the final estimate that might be affecting the convergence, but both are around 25.

Note that there are a large number of variables in the model, hence the size of the b1 matrix, but I have 13 million observations.

[Clyde Schechter]

Yes, this is a legitimate approach. The likelihood function of a negative binomial model is quite complicated and may defy direct attempts at maximization from default starting estimates. A perfectly good way to resolve the non-convergence is to provide starting estimates derived from a simplified approximate model, such as Poisson. Had you not done it that way and just posted about a non-convergence problem you might well have gotten advice to try just this approach.

Coefficients around 25, if they are coefficients of 0/1 variables, or coefficients of variables that are not on a pretty small scale, are suspect. Remember that these models have a log link. So a unit change in one of these variables corresponds, in this model, to multiplying the outcome by exp(25), which is an enormous number, about 7 x 10¹⁰. If, however the scale of these variables is small enough that nothing approaching a unit change in them is possible, then there is no issue. You have to decide based on your knowledge of the content area whether a unit change in a variable leading to multiplying the expected outcome by a factor of exp(25) is reasonable--in most domains this would be quite unreasonable. If you conclude it's unreasonable, then you need to look into why you are getting such results. Most commonly it reflects errors in the data; less commonly it is a modeling problem and the variables need to be eliminated or transformed, or handled in some other way.

If you need more concrete advice, please post actual outputs so that the details of what is happening can be seen.

[Molly Jeffery]

Thank you so much for your help. Is this a problem in some circumstances? I.e., is there a reason they don't default to using the coefficients from the Poisson model to estimate the full model?

On the large coefficients: I don't interpret those coefficients. I have included the entire population in the model so that their data can contribute to the estimates of other variables. The problem coefficients are on a factor variable measuring whether the patient received both drugs from the same physician. The levels are: Yes, No, Unknown, and No overlap in month. So the "no overlap" level is perfectly correlated with the outcome variable. (the coefficients are actually around -25).

A two-part model would be the ideal model, but I'm not aware of a command that supports count models in the second equation. I need the count model because I want to model the rate of drug_1 days on which drug_2 is also available, which I'm doing with the exposure variable in the original model.

I suppose an alternative would be to directly calculate the proportion of overlap days and use it as the dependent variable in a hurdle model with the churdle command, but that does not take into account the limited nature of a proportion (bounded between 0 and 1; variance is mean-dependent).

Thank you again for working through this with me. I appreciate your help.

[Joao Santos Silva]

Dear Molly Jeffery

Adding to the useful advice Clyde provided, I wonder whether you have perfect predictors and poisson and nbreg with the Poisson initial values are giving you false convergence results. I suggest that you try estimating your model with ppml (available from SSC) and see if any variables are dropped. If no variables are dropped, all should be fine.

Best wishes,

Joao

[John Mullahy]

Molly: Adding two more thoughts to those already helpfully provided by Clyde and Joao:

1. A general consideration with nonconvergence of NB algorithms is that the data may be (conditional on x) underdispersed relative to Poisson in which case the NB probability model is problematic. This may not be relevant in your case (esp. since you did ultimately obtain convergence), but it is something to consider in general when estimating NB models.

2. Regarding your query about two-part models in #3, one could follow Cragg's (Econometrica, 1971) approach that involves using truncated distributions for part-2, and with count data use truncated-outcome commands like tpoisson and tnbreg. It's not obvious whether this is a good solution in general, or in your specific case, and interpretation and inference may be challenging. But it is one way to consider a two-part approach. (Link to Cragg's paper: https://www.jstor.org/stable/pdf/1909582.pdf )

[Carlo Lazzaro]

As far as two-parts model are concerned, I would also mention John' s pivotal publication: Specification and testing of some modified count data models. Journal of Econometrics 1986;33(3):341-65.

[Molly Jeffery]

Thank you all very much for your advice, which gives me a lot to think about. Joao and Carlo, thank you for the suggestions--I will take a look at those papers.

John--I think the churdle command (which implements Cragg 1971) may get me most of the way there, with the option for the exponential model for part 2. The goal is to place this in a medical journal, so I have been using predictive margins to ease interpretation anyway.

Thank you all again
Best--
Molly