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

[Roger Bivand]

On Tue, 20 Mar 2012, Adam Slez wrote:

> Roger,
>
> Out of curiosity, did you ever reach a final word on this?

Please check in 0.5-45, released at the weekend and on CRAN as a source 
and Windows binary package. The changelog lists the commit as February 24.

Roger

>
> Adam
>
> On Fri, Feb 24, 2012 at 7:18 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>
>> 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
>>
>>
>

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