negative binomial in GEE

Does anyone have any suggestions for a way round this problem for the
future (I had to resort to using Stata), or maybe more realistically, how
much work it would take to build an extendible version of the gee
"algorithm"?

I don't think it would take that much work [and I've done it a couple of 
times, once in XLISP-Stat and once in SPIDA].  An entirely interpreted 
implementation will be relatively slow, though it might be possible to 
improve that quite a bit with sparse matrix techniques.

 	-thomas