Confusion about Markov random fields in the mgcv package in R

Question

In order to implement a spatial analysis, I tried a simple Markov random field smoother in an example in the mgcv package in R, where the manual is here:

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

This is the example I tried:

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

However, when I took a look at estimated coefficients in b$coefficients, there are 48 estimates from the 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, the district shapes list shows 49 areas (from 0~48). When I tried my data, the same situation happened because data with 28 areas only produced 27 estimates from the Markov random field smoother.

My understanding is, the Markov random field used as a spatial function can be regarded as a structured random effect; however, the Markov random field smoother in the mgcv package in R seems to automatically set up the first area as a reference level. I am wondering whether it is just like a fixed effect but under the consideration of spatial autocorrelation?

If so, an extended problem is how to explain such output? I feel very weird in that the spatial estimate can be explained like the difference between each area and the reference area, but this interpretation is not too meaningful.

I am thinking whether we can fit a Markov random field smoother like a random effect in R. Hope anyone who is familiar with this package can provide some suggestions. Thanks!

There is an update of this question from my misunderstanding to the coefficients object in the mgcv package. Now I may recognize that the output for Markov random fields in b$coefficients is only its knots. I am still looking for where I can get the spatial estimate of each area used to draw the map shown by plot(b). If anyone knows how to extract the spatial estimates, that must be pretty thankful! — cchien, Mar 24 '15 at 16:38
The most updated status: This question has been solved by using predict(object, type="terms", se=T). — cchien, Mar 25 '15 at 17:05

cmbarbu · Answer 1 · 2015-03-23T11:15:47.563

The coefficients in a multivariate Gaussian smoothing are not and should not be interpreted as coefficients individually applied to each covariate s is a function of. They describe the correlation relationship between the covariates; correlation described by a number of coefficients to be set by the k parameter in the s R function.

By default, s sets k to its maximum, n-1. k cannot be bigger than n-1 with n the number of covariates in s as an intercept is necessary to set the average level the smoothing function will vary around and the total number of fitted parameters must not be bigger than the size of the "data".

For further details, check the paragraph on choose.k in the mgcv help file.

If you are interested by something directly applicable to each of your districts, you should check the values predicted by the smoothing function. Following gamObject help it is given by the fitted.values item.

Here I 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 to visualize the fit, it looks good and the correspondence between observed and estimated seems reasonable: plot(columb$crime,b$fitted.values)

I know values in b$fitted.values are predictions, which are definitely 49 values with no doubt in this example; however, what I mean is the values in b$coefficients which contains all estimates including linear and nonlinear estimated coefficients. In this example, a total of 49 values in b$coefficients, but the first value is the estimated intercept, and the remaining 48 estimated coefficients are from s(district). So this is why I asked that there is one missing estimate in s(district). — cchien, Mar 22 '15 at 06:56

Confusion about Markov random fields in the mgcv package in R

1 Answers1