The R Project SVN R-packages

\name{smooth.construct.mrf.smooth.spec}
\alias{smooth.construct.mrf.smooth.spec}
\alias{Predict.matrix.mrf.smooth}
\alias{mrf}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Markov Random Field Smooths}

\description{For data observed over discrete spatial units, a simple Markov random field 
smoother is sometimes appropriate. These functions provide such a smoother class for \code{mgcv}. 
See details for how to deal with regions with missing data.
}

\usage{
\method{smooth.construct}{mrf.smooth.spec}(object, data, knots)
\method{Predict.matrix}{mrf.smooth}(object, data)
}

\arguments{ 
\item{object}{For the \code{smooth.construct} method a smooth specification object, 
usually generated by a term \code{s(x,...,bs="mrf",xt=list(polys=foo))}. \code{x} is a factor variable giving labels 
for geographic districts, and the \code{xt} argument 
is obligatory: see details. For the \code{Predict.Matrix} method 
an object of class \code{"mrf.smooth"} produced by the \code{smooth.construct} method.}

\item{data}{a list containing just the data (including any \code{by} variable) required by this term, 
            with names corresponding to \code{object$term} (and \code{object$by}). The \code{by} variable 
            is the last element.} 

\item{knots}{If there are more geographic areas than data were observed for, then this argument is used to 
             provide the labels for all the areas (observed and unobserved). }

}

\value{ An object of class \code{"mrf.smooth"} or a matrix mapping the coefficients of the MRF smooth 
to the predictions for the areas listed in \code{data}.
}

\details{A Markov random field smooth over a set of discrete areas is defined using a set of area labels, and 
a neighbourhood structure for the areas. The covariate of the smooth is the vector of area labels 
corresponding to each obervation. This covariate should be a factor, or capable of being coerced to a factor.

The neighbourhood structure is supplied in the \code{xt} argument to \code{s}. This must contain at least one of 
the elements \code{polys}, \code{nb} or \code{penalty}. 
\describe{
\item{polys}{contains the polygons defining the geographic areas. 
             It is a list with as many elements as there are geographic areas. 
             \code{names(polys)} must correspond to
             the levels of the argument of the smooth, in any order (i.e. it gives the area labels). 
             \code{polys[[i]]} is a 2 column matrix the rows of which specify the vertices of the polygon(s) 
             defining the boundary of the ith area. A boundary may be made up of several closed loops: these must
             be separated by \code{NA} rows. A polygon within another is treated as a hole. The first polygon in 
             any \code{polys[[i]]}  should not be a hole.  An example 
             of the structure is provided by \code{\link{columb.polys}} (which contains an artificial hole
             in its second element, for illustration). Any list elements with duplicate names are combined into a 
             single NA separated matrix.
             
             Plotting of the smooth is not possible without a \code{polys} object.

             If \code{polys} is the only element of \code{xt} provided, then the neighbourhood structure is 
             computed from it automatically. To count as neigbours, polygons must exactly share one of more 
             vertices.
      }
\item{nb}{is a named list defining the neighbourhood structure. \code{names(nb)} must correspond to the 
          levels of the covariate of the smooth (i.e. the area labels), but can be in any order. \code{nb[[i]]} 
          is a numeric vector indexing the neighbours of the ith area (and should not include \code{i}).
      All indices are relative to \code{nb} itself, but can be translated using \code{names(nb)}. See example code.
      As an alternative each \code{nb[[i]]} can be an array of the names of the neighbours, but these will be
      converted to the arrays of numeric indices internally. 

          If no \code{penalty} is provided then it is computed automatically from this list. The ith row of 
          the penalty matrix will be zero everwhere, except in the ith column, which will contain the number 
          of neighbours of the ith geographic area, and the columns corresponding to those geographic 
          neighbours, which will each contain -1.  
         }
\item{penalty}{ if this is supplied, then it is used as the penalty matrix. It should be positive semi-definite.
Its row and column names should correspond to the levels of the covariate.}
}

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 (2017) Section 5.4.2. See Wood (2017, section 5.8.1). 

Note that smooths of this class have a built in plot method, and that the utility function \code{\link{in.out}} 
can be useful for working with discrete area data. The plot method has two schemes, \code{scheme==0} is colour, 
\code{scheme==1} is grey scale.

The situation in which there are areas with no data requires special handling. You should set \code{drop.unused.levels=FALSE} in 
the model fitting function, \code{\link{gam}}, \code{\link{bam}} or \code{\link{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. 

}

\references{
Wood S.N. (2017) Generalized additive models: an introduction with R (2nd edition). CRC.
}

\author{ Simon N. Wood \email{simon.wood@r-project.org} and Thomas Kneib
(Fabian Scheipl prototyped the low rank MRF idea) } 

\seealso{\code{\link{in.out}}, \code{\link{polys.plot}}}

\examples{
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)

## Examine an example neighbourhood list - this one auto-generated from
## 'polys' above.

nb <- b$smooth[[1]]$xt$nb 
head(nb) 
names(nb) ## these have to match the factor levels of the smooth
## look at the indices of the neighbours of the first entry,
## named '0'...
nb[['0']] ## by name
nb[[1]]   ## same by index 
## ... and get the names of these neighbours from their indices...
names(nb)[nb[['0']]]
b1 <- gam(crime ~ s(district,bs="mrf",k=20,xt=list(nb=nb))+s(income),
         data=columb,method="REML")
b1 ## fit unchanged
plot(b1) ## but now there is no information with which to plot the mrf
}
\keyword{models} \keyword{regression}%-- one or more ..