They are different ways of doing the same thing. 

I would guess that Professor Bivand is concerned that it might not be possible to compute (I - pW)^-1 since W is usually sparse but (I-pW)^-1  is often dense and impossible to compute for large N (unless you have unlimited RAM). Also this way professor Bivand is making use of the native R packages, namely lm which already has an lm.fit and lm.residuals functions ! (very clever)

Professor LeSage does it assuming smaller problems or if you are reading the book "Introduction to Spatial Econometrics" he does it so that it is simpler.

But if you think of the way OLS is estimated in matrix form you will see that it is the same thing.

Finally i like Professor's Bivand solution more. If you think about it for large N (I-pW)^-1 can pose numerical difficulties in terms of precision ...

If you want a more