[R] MRF smoothers in MGCV - specifying neighbours

[Roger Bivand]

Mark Payne <markpayneatwork <at> gmail.com> writes:

> 
> Hi,
> 
> Does anyone have an example of a Markov Random Field smoother (MRF) in MGCV
> where they have specified the neighbourhood directly, rather than supplying
> polygons? Does anyone understand how the rules should be? Based on the
> columb example, I have setup my data set and neighbourhood like so:
> 
> > head(nb.l)
> $`10/10`
> [1] 135 155 153
> 
...
> 
> > head(obs)
>             x          y      truth x.idx y.idx xy.idx
> 24  1.4835147  0.8026673  2.3605204    13    10  13/10
...
> 
> but get the following error:
> 
> > mdl <- gam(truth ~
> s(xy.idx,bs="mrf",xt=list(nb=nb.l)),data=obs,method="REML")
> Error in smooth.construct.mrf.smooth.spec(object, dk$data, dk$knots) :
>   mismatch between nb/polys supplied area names and data area names

I could reconstruct your problem with:

library(mgcv)
library(spdep)
example(columbus)
nb <- poly2nb(columbus)
names(nb) <- attr(nb, "region.id")
columbus$ID <- names(nb)
bnb <- gam(CRIME ~ s(ID, bs="mrf", xt=list(nb=nb)), data=columbus,  
method="REML")

#Error in ExtractData(object, data, knots) : 
#  'names' attribute [1] must be the same length as the vector [0]

It appears that the ID variable is happier as integer:

columbus$ID <- as.integer(names(nb))
bnb <-  gam(CRIME ~ s(ID, bs="mrf", xt=list(nb=nb)), data=columbus,
method="REML")
summary(bnb)

Using another representation:

nb <- poly2nb(columbus, row.names=columbus$NEIGNO)
names(nb) <- attr(nb, "region.id")
bnb <-  gam(CRIME ~ s(NEIGNO, bs="mrf", xt=list(nb=nb)), data=columbus,
method="REML")

also works, but NEIGNO is numeric. 

It appears that the ID variable should be numeric (or integer), but that the
names of the list of neighbour vector should be character, but otherwise
match. The call tree is not simple; the error seems to occur in
mgcv:::ExtractData.

Hope this helps,

Roger

> 
> However, there is a perfect match between the nb list names and the data
> area names:
> > all(levels(obs$xy.idx) %in% names(nb.l))
> [1] TRUE
> >
> 
> Any suggestions where to start?
> 
> Mark
> 
> 	[[alternative HTML version deleted]]
> 
>