[R-sig-Geo] bad degrees of freedom in bptest.sarlm

[Roger Bivand]

On Fri, 24 Feb 2012, Adam Slez wrote:

> 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)

Thanks, k and n should be set after the if(), not before, I think. I'll 
check and release if this seems correct.

Roger

>   ...
>    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]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no