bad degrees of freedom in bptest.sarlm

Am I correct in thinking that the degrees of freedom associated with the
bptest.sarlm, much like the bptest code from which it came, should be
driven by the number of parameters in the auxiliary regression NOT the
model used to generate the residuals?  The way the code for bptest.sarlm is
written, the degrees of freedom are calculated from the original
regression, as opposed to the auxiliary model.  The easiest way to see this
is to do something like this:

library(spdep)
data(columbus)
wgts<-nb2listw(col.gal.nb)
model1<-lagsarlm(CRIME~INC, data = columbus, listw = wgts)
model2<-lagsarlm(CRIME~INC+HOVAL, data = columbus, listw = wgts)
bp1<-bptest.sarlm(model1, ~EW, data = columbus, studentize = F)
bp2<-bptest.sarlm(model2, ~EW, data = columbus, studentize = F)
bp3<-bptest.sarlm(model2, ~EW+CP+NSA, data = columbus, studentize = F)

All I'm doing here is running spatially adjusted BP tests on two different
spatial lag models.  In the first two cases, I'm just testing for groupwise
heteroskedasticity with regimes defined in terms of the EW variable.  If I
understand the logic of the BP test, the degrees of freedom associated with
bp1 and bp2 should be the same, while the number of degrees of freedom
associated with bp2 and bp3 should differ.  This doesn't seem to be the
case.

If we repeat this exercise for non-spatial models and tests, changes in the
degrees of freedom work as expected (recognizing of course that the actual
degrees of freedom for a non-spatial BP tests differs from the degrees of
freedom for a spatially-adjusted BP test):

library(lmtest)
model3<-lm(CRIME~INC, data = columbus)
model4<-lm(CRIME~INC+HOVAL, data = columbus)
bp4<-bptest(model3, ~EW, data = columbus, studentize = F)
bp5<-bptest(model4, ~EW, data = columbus, studentize = F)
bp6<-bptest(model4, ~EW+CP+NSA, data = columbus, studentize = F)

The bptest.sarlm routine is cribbed directly from the bptest code.  If we
take a look a the first part of the bptest.sarlm code, we can see the
problem (I'm using ... to denote skipped code):

function (object, varformula = NULL, studentize = TRUE, data = list()) {
   if (!inherits(object, "sarlm"))
       stop("not sarlm object")
   Z <- object$tarX
   k <- ncol(Z)
   n <- nrow(Z)
   if (!is.null(varformula))
       Z <- model.matrix(varformula, data = data)

  ...
   df <- k - 1
  ...
}

It is pretty clear that df is driven by k which is unaffected by Z <-
model.matrix(varformula, data = data), even when varformula is not null.  I
can't see any reason why the degrees of freedom for a spatially-adjusted BP
test would follow from the regression producing the residuals as opposed to
the auxiliary regression.  Any help would be much appreciated.

Adam

	[[alternative HTML version deleted]]