I am attempting to do a 2 stage residual inclusion (2SRI) model where the first stage is estimating a multinomial regression. The problem I am having is figuring out how to calculate the Pearson residual to enter into my second stage. I have seen one other unanswered post related to this question.

Let's say my data looks like this:

ID, provider, copay, opioid
1, MD1, 50, yes
2, MD2, 10, no
3, MD1, 25, no
4, MD3, 10, no
5, MD2, 30, yes
6, MD3, 14, yes

Using MD1 as the base comparison, the first stage estimation would ideally tell me how a dollar increase in copay influences the likelihood that you see MD1 versus MD2, MD1 versus MD3, etc. In Stata speak:

mlogit provider copay
predict res, rstandard

Ideally I would be able to predict the residual as above but it seems that I need to manually calculate this residual rather than rely on a postestimation option to predict the residual. And I'm pretty sure I need the Pearson residual though I am open to the fact I may be wrong on this. If someone could help me figure out the most efficient approach to the calculation, I would greatly appreciate it.

And just to complete the thought, in the second stage, I want to know how the provider choice influence opioid outcomes, and let's assume that provider choice is influenced by copay - hence the need for the first stage.

logit opioid provider res

Thank you!
Bianca