I went to look into the code in detail and actually my answer was not correct. Sorry i hasted and i should not have.
The residuals that each of them calculating are different. Here is why:
The model is as follows:
y = pWy + xb + e with e ~ n(0,1)
Now if we play arround with it we get:
y = (I - pW)^-1(xb + e)
Now what Prof. LeSage does ie:
y - (I-p_hat * W)^-1 * xb_hat = (I-pW)^-1*e
So what you are getting it the residual with the auto correlation.
On the other hand, by transforming y:
y - p_hatWy = xb + e
Estimating, xb and calculation the residuals, what Bivand is doing is giving you e instead of (I-pW)^-1*e
Which one is preferred will depend on your application!