[R] MRF smoothers in MGCV - specifying neighbours

[Simon Wood]

Hi Mark,

I'm not sure what is happening here - there is no chance that nb.l 
contains a neighbourhood not in the levels of obs$xy.idx, I suppose? 
i.e. is

all(names(nb.l)%in%levels(obs$xy.idx))

also TRUE? Here is some code illustrating what nb should look like (and 
in response to Roger Bivand's suggestion I also tried this replacing all 
the labels with things like "x/y", but it still works).


## example mrf fit using polygons....
library(mgcv)
## Load Columbus Ohio crime data (see ?columbus for details and credits)
data(columb)       ## data frame
data(columb.polys) ## district shapes list
xt <- list(polys=columb.polys) ## neighbourhood structure info for MRF
par(mfrow=c(2,2))
## First a full rank MRF...
b0 <- gam(crime ~ s(district,bs="mrf",xt=xt),data=columb,method="REML")

## same fit based on direct neighbour spec...
nb <- mgcv:::pol2nb(columb.polys)$nb
xt <- list(nb=nb)
b <- gam(crime ~ s(district,bs="mrf",xt=xt),data=columb,method="REML")

best,
Simon



On 08/05/14 01:58, Mark Payne wrote:
> 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
>
> $`10/2`
> [1] 27  8  6
>
> $`10/3`
> [1] 48  7 28 26
>
> $`10/4`
> [1] 69 27 49 47
>
> $`10/5`
> [1] 48 70 68
>
> $`10/7`
> [1] 115  95  93
>
>> head(obs)
>              x          y      truth x.idx y.idx xy.idx
> 24  1.4835147  0.8026673  2.3605204    13    10  13/10
> 26  1.0452111  0.4673685  1.8316741    11     8   11/8
> 43  2.1514977 -0.2640058 -2.8812026    17     5   17/5
> 46  2.8473951  0.5445714  3.6347799    20     9   20/9
> 53  1.7983253 -0.6905912 -2.5473984    15     3   15/3
> 86 -0.1839814 -0.7824026 -0.5776616     5     2    5/2
>>
>
> 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
>
> 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]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>


-- 
Simon Wood, Mathematical Science, University of Bath BA2 7AY UK
+44 (0)1225 386603               http://people.bath.ac.uk/sw283