
mgcv
.
See details for how to deal with regions with missing data.## S3 method for class 'mrf.smooth.spec':
smooth.construct(object, data, knots)
## S3 method for class 'mrf.smooth':
Predict.matrix(object, data)
smooth.construct
method a smooth specification object,
usually generated by a term s(x,...,bs="mrf",xt=list(polys=foo))
. x
is a factor variable giving labels
for geographic districts, and the xt
by
variable) required by this term,
with names corresponding to object$term
(and object$by
). The by
variable
is the last element."mrf.smooth"
or a matrix mapping the coefficients of the MRF smooth
to the predictions for the areas listed in data
.The neighbourhood structure is supplied in the xt
argument to s
. This must contain at least one of
the elements polys
, nb
or penalty
.
[object Object],[object Object],[object Object]
If no basis dimension is supplied then the constructor produces a full rank MRF, with a coefficient for each geographic area. Otherwise a low rank approximation is obtained based on truncation of the parameterization given in Wood (2006) Section 4.10.4.
Note that smooths of this class have a built in plot method, and that the utility function in.out
can be useful for working with discrete area data. The plot method has two schemes, scheme==0
is colour,
scheme==1
is grey scale.
The situation in which there are areas with no data requires special handling. You should set drop.unused.levels=FALSE
in
the model fitting function, gam
, bam
or gamm
, having first ensured that any fixed effect
factors do not contain unobserved levels. Also make sure that the basis dimension is set to ensure that the total number of
coefficients is less than the number of observations.
in.out
, polys.plot
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...
b <- gam(crime ~ s(district,bs="mrf",xt=xt),data=columb,method="REML")
plot(b,scheme=1)
## Compare to reduced rank version...
b <- gam(crime ~ s(district,bs="mrf",k=20,xt=xt),data=columb,method="REML")
plot(b,scheme=1)
## An important covariate added...
b <- gam(crime ~ s(district,bs="mrf",k=20,xt=xt)+s(income),
data=columb,method="REML")
plot(b,scheme=c(0,1))
## plot fitted values by district
par(mfrow=c(1,1))
fv <- fitted(b)
names(fv) <- as.character(columb$district)
polys.plot(columb.polys,fv)
Run the code above in your browser using DataLab