Getting predicted values from xtpoisson (fixed effects)

[Bob Reed]

Hi,

I want to get fitted values from xtpoisson, fe. I can best illustrate my problem with the following data and code. Open up the dataset "ForStata". Then run the program below, and compare yvar with the 5 predicted values: yhat1-yhat5.

yvar is the dependent variable.

yhat1 comes from a GLM/family(poisson) link(log) model. Note how close yhat1 is to the dep. variable = yvar.

yhat2 and yhat3 come from an xtpoisson model (random effects). Note that the estimates are the same as the GLM estimates, and that yhat2=yhat3=yhat1.

yhat4 and yhat5 come from an xtpoisson model (fixed effects). Note that yhat4 and yhat5 are very different from the previous predicted values, and are not close to the dep variable.

QUESTION: Is there any way of getting predicted values from the xtpoisson model (fixed effects) so I can see how well it does in predicting the original dep variable and examine the residuals?

PROGRAM:

gen yvar = exports_ttl
glm exports_ttl ln_RGDP i.bothin i.year i.countrycode RER sd_RER, family(poisson) link(log)
predict yhat1
xtpoisson exports_ttl ln_RGDP i.bothin i.year i.countrycode RER sd_RER
predict yhat2, nu0
predict yhat3, iru0
xtpoisson exports_ttl ln_RGDP i.bothin i.year i.countrycode RER sd_RER, fe
predict yhat4, nu0
predict yhat5, iru0
// Compare yvar (the dep var) with the four predicted values in the data editor (browser)

[Joao Santos Silva]

Dear Bob,

So many people ask that question that I prepared an example illustrating how to do it. The example uses the auto dataset from Stata, but you can easily modify it to fit your case:

sysuse auto, clear
g id=rep78
drop if id==.
xtset id
xtpoisson price mpg, fe
predict double fitted, xb
gen double yhat=exp(fitted)
 egen meany=mean(price) , by(id)
 egen meanyhat=mean(yhat), by(id)
 gen double exp_alpha=meany/meanyhat if meanyhat>0
 replace yhat=yhat*exp_alpha if meanyhat>0
su price yhat

All the best,

Joao

[Bob Reed]

Thank you Joao. That is really, really helpful. :-)

[Stephen Jenkins]

I agree: very helpful. Bob: I added the following to Joao's code, looking to visualise your "problem"

. ge d = yhat - price

. scatter d price

Clearly a regression line through the scatterplot would have a negative slope. Whether this, and the magnitudes shown, are bothersome in your particular context is your call, I think

[Bob Reed]

Thanks Stephen. Plotting residuals against explanatory variables seems to be a useful diagnostic. Can you recommend a good reference for how to conduct diagnostics with nonlinear models such as xtpoisson?

[Joao Santos Silva]

Dear Bob,

There is a literature on that, including the goodness-of-link test, but I would start with the simple RESET test. We have a description of how to implement it in the "log of gravity" page.

Best regards,

Joao

[Bob Reed]

Thanks, Joao. And suppose all the usual estimators (OLS with logged dep variable, GLM with different families/links, xtpoisson) fail the RESET test (which I think would not be especially unusual). Any words of wisdom what to do then? :-)

[Joao Santos Silva]

Bob,

If you are estimating a gravity equation you really need an exponential model as estimated by the commands in the Poisson family. If the model fails the RESET (with clustered robust SEs) you should try to change the set of regressors or include interactions, non-linear transformations, etc.

Best of luck,

Joao

[Bob Reed]

JOAO: Thanks so much!

[Bob Reed]

Dear Joao, I seem to have run into a problem implementing your technique. The problem is that when I check your procedure, I get contradictory answers, as the program below show.

As before, open up the dataset "ForStata." Then run the program below, and compare yhat1 with yhat2. yhat1 is the estimate I get from Stata using the unconditional FE version of xtpoisson. yhat2 is what I get when I apply your procedure to the conditional FE version. After you run the program, go to the data browser and compare yval with yhat1 and yhat2 (note that yval is the dependent variable that I am trying to predict).

--------------------

// First I get a prediction using unconditional fixed effects
// I call these yhat1
xtpoisson exports_ttl ln_RGDP i.countrycode i.year, vce(robust)
gen bob = e(sample)
predict double yhat1 , nu

// Next I get prediction using conditional fixed effects
// Check to make sure that the estimates are identical with the unconditional FE xtpoisson model above (they are)
xtpoisson exports_ttl ln_RGDP i.year, fe vce(robust)
predict double fitted , nu0
gen double yhat2=exp(fitted)
egen meany=mean(exports_ttl) if bob == 1 , by(countrycode)
egen meanyhat=mean(yhat2) if bob == 1, by(countrycode)
gen double exp_alpha=meany/meanyhat if meanyhat>0 & bob == 1
replace yhat2=yhat2*exp_alpha if meanyhat>0
gen yval = exports_ttl if bob == 1

---------------

You will see that yhat1 is different from yhat2. :-(

Now you might ask me, then why even estimate the conditional FE model. Just go with the unconditional FE model and get my estimates from there. The problem is that there are a number of cases which won't converge when I attempt to estimate the unconditional FE model, but do converge when I estimate the conditional FE model.

So I want to estimate the conditional FE model and use your procedure to get predicted values, but I am bothered by the fact the differences between yhat1 and yhat2.

Any help you could give would be greatly appreciated!

[Joao Santos Silva]

Bob,

You were not using my code! Check it like this:

// First I get a prediction using unconditional fixed effects
// I call these yhat1
poisson exports_ttl ln_RGDP i.countrycode i.year, vce(robust)
gen bob = e(sample)
predict double yhat1

// Next I get prediction using conditional fixed effects
// Check to make sure that the estimates are identical with the unconditional FE xtpoisson model above (they are)
xtpoisson exports_ttl ln_RGDP i.year, fe vce(robust)
predict double fitted, xb
gen double yhat2=exp(fitted)
egen meany=mean(exports_ttl) if bob == 1 , by(countrycode)
egen meanyhat=mean(yhat2) if bob == 1, by(countrycode)
gen double exp_alpha=meany/meanyhat if meanyhat>0 & bob == 1
replace yhat2=yhat2*exp_alpha if meanyhat>0
gen yval = exports_ttl if bob == 1

Joao

[Bob Reed]

Joao, Guilty as charged! I changed it from what you originally had to fit my problem (or so I thought), but now see that I shouldn't have. Thank you very much!

[giorgioconti]

Dear Bob,

So many people ask that question that I prepared an example illustrating how to do it. The example uses the auto dataset from Stata, but you can easily modify it to fit your case:

sysuse auto, clear
g id=rep78
drop if id==.
xtset id
xtpoisson price mpg, fe
predict double fitted, xb
gen double yhat=exp(fitted)
egen meany=mean(price) , by(id)
egen meanyhat=mean(yhat), by(id)
gen double exp_alpha=meany/meanyhat if meanyhat>0
replace yhat=yhat*exp_alpha if meanyhat>0
su price yhat

All the best,

Joao

Dear Joao,

I have the same question as Bob, but my model is random effects negative binomial (xtnbreg, re). Could you please tell me whether the code for obtaining fitted values (you provided above) will be the same for xtnbreg, re, i.e.

xtnbreg price mpg, re
predict double fitted, xb
gen double yhat=exp(fitted)
egen meany=mean(price), by(id)
egen meanyhat=mean(yhat), by(id)
gen double exp_alpha=meany/meanyhat if meanyhat>0
replace yhat=yhat*exp_alpha if meanyhat>0

The problem with my data is that when I try to predict yhat with predict yhat, nu0, Stata gives me the following error
"warning: estimated value of r is less than one; expected value undefined."
When I use your syntax, I get the predicted number of events without any problem.

I would be very grateful for your reply.

Best regards,
Giorgio

[Joao Santos Silva]

Dear Giorgio,

The code I provided is for FE only. Prediction with random effects is quite different and I am not familiar with what -xtnbreg, fe- does, so I am afraid I cannot help. I have a question, however: are you sure you need a NegBin model with RE?

Best wishes,

Joao

[giorgioconti]

Dear Joao,

Thank you for your reply! I have many panels with few time periods, so I was thinking of using random effects in order not to reduce degree of freedom.

Perhaps the options may be using cluster-robust standard errors instead of RE. This line of syntax works fine:

nbreg price mpg, vce(cluster id)
predict fitted, n