plotting predicted probabilities after poisson regression with cubic splines

[Daniela van Santen]

Hi,

I need to plot the predicted incidence after running a Poisson regression model when using cubic splines. I am using Stata 13. So basically I want my graph to show the predicted incidence (y axis) over calendar year (x axis).

I have aggregated data of individuals distributed over 19 years, with the dataset containing person years of follow up, calendar year and number of events.

example dataset:
year person_timeyears failures
2009 45.0 1
2010 225.5 6
2011 361.8 7
2012 829.8 21
2013 1059.5 11
..
..
..

My model looks as follows:

______

** First making splines of calendar year with 3 knots:
mkspline calendarsp =year, cubic knots(2009,2012,2017) displayknots
mat knots = r(knots)

** running the Poisson model
poisson failures calendarsp*, e(person_timeyears)

______

I know that when using predict I need to calculate the confidence intervals in a different way than when using a linear regression. At least that what I have been reading in Stata forums. One of the links I used for an example code is: https://www.stata.com/support/faqs/s...ence-intervals.

I tried to use the following code to make the graph, but as you can see in the figure I've attached something goes wrong with the 95%CI. I think this is because I am using splines. I was also wondering whether the red line is at least the predicted incidence. Not sure if I am obtaining the results I want.

______

** code to make a graph:

predict yhat
predict xb,xb
predict error, stdp

generate lb = xb - invnormal(0.975)*error
generate ub = xb + invnormal(0.975)*error
generate plb = invlogit(lb)
generate pub = invlogit(ub)

twoway rarea plb pub year, sort || ///
line yhat year, sort

_____

I hope someone can help me. Thank you!!

Kind regards , Daniela.

[Joseph Coveney]

Why are you using an inverse-logit transformation for the confidence bounds after fitting a generalized linear model where the link function is log?

[Daniela van Santen]

Hi Joseph,

Thank you for your question. Basically I am not sure how to calculate the 95%CIs as you can see. I've been reading forums and it just confused me on what to use exactly - I mistakenly read the logistic regression forum post and applied that to my example here. Should I just be using the following code for the predicted incidence with 95% CIs?:

predict yhat
predict se, stdp
gen lb = yhat - 1.96*se
gen ub = yhat + 1.96*se

However, is I apply this, the 95%CIs are very narrow.

Thank you again!

[Joseph Coveney]

Maybe you could try something like that below (begin at the "Begin here" comment after fitting the Poisson regression model). I don't have Stata Release 13 anymore, but the postestimation commands shown ought to work with that version. One or both of the two graphs are what you seem to be after.

version 15.1

clear *

set seed `=strreverse("1450624")'

quietly set obs 19

generate int predictor = 2008 - 9 + _n
label variable predictor Year
generate double exposure = 1000 * runiform()
generate int response = rpoisson(10)

mkspline s = predictor, cubic nknots(3)

poisson response c.s?, exposure(exposure) nolog

*
* Begin here
*
predict double xb, xb nooffset
predict double std, stdp
generate double lb = xb - invnormal(0.975) * std
generate double ub = xb + invnormal(0.975) * std
format xb ?b %3.1f

predict double ir, ir
generate double lir = exp(lb)
generate double uir = exp(ub)
format ir *ir %04.2f

graph twoway ///
    line ir predictor, lcolor(black) lpattern(solid) || ///
    line lir predictor, lcolor(black) lpattern(dash) || ///
    line uir predictor, lcolor(black) lpattern(dash) ///
        xlabel(2000(3)2018) ///
        ytitle(Incidence Rate) ylabel( , angle(horizontal) nogrid) ///
        legend(off)
graph export ir.png, replace

graph twoway ///
    line xb predictor, lcolor(black) lpattern(solid) || ///
    line lb predictor, lcolor(black) lpattern(dash) || ///
    line ub predictor, lcolor(black) lpattern(dash) ///
        xlabel(2000(3)2018) ///
        ytitle(Linear Prediction (Exposure-corrected)) ylabel( , angle(horizontal) nogrid) ///
        legend(off)
graph export xb.png, replace

exit

[Daniela van Santen]

Hi Joseph,

The code seems to be working. Thank you so much!

I just had an additional question regarding your code and I hope you don't mind me asking.

What is the second graph doing? it says exposure corrected, but I am not sure what that means. The y-axis give me negative numbers, so I don't know how to interpret this. Ive attached the graph in this post.

Kind regards,
Daniela.

[Joseph Coveney]

What is the second graph doing? it says exposure corrected, but I am not sure what that means. The y-axis give me negative numbers, so I don't know how to interpret this.

It looks like you didn't need splines at all.

Anyway, the answers to your questions are in the help file for the Poisson regression command's postestimation commands.

help poisson postestimation

and then click on the hyperlink for predict on the help-file screen that pops up in the Viewer window.