spatial - Confusion about Markov random fields in the mgcv package in R -

in order implement spatial analysis, tried simple markov random field smoother in example in mgcv package in r, manual here:

https://stat.ethz.ch/r-manual/r-devel/library/mgcv/html/smooth.construct.mrf.smooth.spec.html

this example tried:

library(mgcv) data(columb)         ## data frame data(columb.polys)   ## district shapes list xt <- list(polys=columb.polys) ## neighbourhood structure info mrf b <- gam(crime ~ s(district,bs="mrf",xt=xt),data=columb,method="reml") 

however, when took @ estimated coefficients in b$coefficients, there 48 estimates markov random field smoother:

> b$coefficients  (intercept)  s(district).1  s(district).2  s(district).3  s(district).4  35.12882390   -10.96490165    20.99250496    16.04968951    10.49535483   s(district).5  s(district).6  s(district).7  s(district).8  s(district).9 16.56626217    14.55352540    17.90043996    -0.60239588    13.41215603  s(district).10 s(district).11 s(district).12 s(district).13 s(district).14     18.61920671   -11.13853418    -2.95677338     7.89719220     3.04717540  s(district).15 s(district).16 s(district).17 s(district).18 s(district).19    -11.18235328    12.57473374    19.83013619    10.56130003    12.36240748  s(district).20 s(district).21 s(district).22 s(district).23 s(district).24     15.65160761    20.40965885    24.79853590     0.05312873   -14.65881026  s(district).25 s(district).26 s(district).27 s(district).28 s(district).29    -13.01294201     7.16191556    -9.36311304     3.65410713   -16.37092777  s(district).30 s(district).31 s(district).32 s(district).33 s(district).34     11.23500771    13.92036006   -14.67653893   -12.39341674    11.02216471  s(district).35 s(district).36 s(district).37 s(district).38 s(district).39    -12.93210046   -15.48924425     3.42745125    -2.54916472    -1.90604972  s(district).40 s(district).41 s(district).42 s(district).43 s(district).44    -16.25160966    -7.46491914    -4.48126353    -7.61064264    -2.91807488  s(district).45 s(district).46 s(district).47 s(district).48    -12.12765102     6.68446503     2.55883220    -0.20920888  

however, district shapes list shows 49 areas (from 0~48). when tried data, same situation happened because data 28 areas produced 27 estimates markov random field smoother.

my understanding is, markov random field used spatial function can regarded structured random effect; however, markov random field smoother in mgcv package in r seems automatically set first area reference level. wondering whether fixed effect under consideration of spatial autocorrelation?

if so, extended problem how explain such output? feel weird in spatial estimate can explained difference between each area , reference area, interpretation not meaningful.

i thinking whether can fit markov random field smoother random effect in r. hope familiar package can provide suggestions. thanks!

the coefficients in multivariate gaussian smoothing not , should not interpreted coefficients individually applied each covariate s function of. describe correlation relationship between covariates; correlation described number of coefficients set k parameter in s r function.

by default, s sets k maximum, n-1. k cannot bigger n-1 n number of covariates in s intercept necessary set average level smoothing function vary around , total number of fitted parameters must not bigger size of "data".

for further details, check paragraph on choose.k in mgcv file.

if interested directly applicable each of districts, should check values predicted smoothing function. following gamobject given fitted.values item.

here get:

   > b$fitted.values    [1] 18.81758 22.12502 30.13315 33.14305 44.11208 30.17184 20.96227 39.77438    [9] 35.64875 32.88071 54.08242 49.13961 43.58527 49.65618 47.64344 50.99036    [17] 32.48752 46.50207 51.70913 21.95138 40.98711 36.13709 21.90757 45.66465    [25] 52.92006 43.65122 45.45233 48.74153 53.49958 57.88845 18.43111 20.07698    [33] 40.25183 23.72681 36.74403 16.71899 44.32493 47.01028 18.41338 20.69650    [41] 20.15782 17.60067 36.51737 30.54075 31.18387 16.83831 25.62500 28.60866    [49] 25.47928 

the result of plot(b) allows visualize fit, looks , correspondence between observed , estimated seems reasonable: plot(columb$crime,b$fitted.values)