[R-sig-Geo] bad degrees of freedom in bptest.sarlm
Roger Bivand
Roger.Bivand at nhh.no
Tue Mar 20 20:39:53 CET 2012
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
More information about the R-sig-Geo
mailing list