The R Project SVN R-packages

\name{smooth.construct.re.smooth.spec}
\alias{smooth.construct.re.smooth.spec}
\alias{Predict.matrix.random.effect}
%- Also NEED an `\alias' for EACH other topic documented here.
\title{Simple random effects in GAMs}

\description{\code{\link{gam}} can deal with simple independent random effects, by exploiting the link 
between smooths and random effects to treat random effects as smooths. \code{s(x,bs="re")} implements 
this. Such terms can can have any number of predictors, which can be any mixture of numeric or factor 
variables. The terms produce a parametric interaction of the predictors, and penalize the corresponding 
coefficients with a multiple of the identity matrix, corresponding to an assumption of i.i.d. normality.
See details.

}

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

\arguments{ 
\item{object}{For the \code{smooth.construct} method a smooth specification object, 
usually generated by a term \code{s(x,...,bs="re",)}. For the \code{predict.Matrix} method 
an object of class \code{"random.effect"} 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}{generically a list containing any knots supplied for basis setup --- unused at present.}

}

\value{ An object of class \code{"random.effect"} or a matrix mapping the coefficients of the random effect to the random effects themselves.
}

\details{Exactly how the random effects are implemented is best seen by example. Consider the model 
term \code{s(x,z,bs="re")}. This will result in the model matrix component corresponding to \code{~x:z-1} 
being added to the model matrix for the whole model. The coefficients associated with the model matrix 
component are assumed i.i.d. normal, with unknown variance (to be estimated). This assumption is 
equivalent to an identity penalty matrix (i.e. a ridge penalty) on the coefficients. Because such a 
penalty is full rank, random effects terms do not require centering constraints. 

If the nature of the random effect specification is not clear, consider a couple more examples: 
\code{s(x,bs="re")} results in \code{model.matrix(~x-1)} being appended to the overall model matrix, 
while  \code{s(x,v,w,bs="re")} would result in  \code{model.matrix(~x:v:w-1)} being appended to the model 
matrix. In both cases the corresponding model coefficients are assumed i.i.d. normal, and are hence 
subject to ridge penalties.

Some models require differences between the coefficients corresponding to different levels of the same random effect. See  \code{\link{linear.functional.terms}} for how to implement this.

If the random effect precision matrix is of the form \eqn{\sum_j \lambda_j S_j}{sum_j p_j S_j} for known matrices \eqn{S_j}{S_j} and unknown parameters \eqn{\lambda_j}{p_j}, then a list containing the \eqn{S_j}{S_j} can be supplied in the \code{xt} argument of 
\code{\link{s}}. In this case an array \code{rank} should also be supplied in \code{xt} giving the ranks of the \eqn{S_j}{S_j} matrices. See simple example below. 

  
Note that smooth \code{id}s are not supported for random effect terms. Unlike most smooth terms, side 
conditions are never applied to random effect terms in the event of nesting (since they are identifiable 
without side conditions).

Random effects implemented in this way do not exploit the sparse structure of many random effects, and 
may therefore be relatively inefficient for models with large numbers of random effects, when \code{gamm4}
or \code{\link{gamm}} may be better alternatives. Note also that \code{\link{gam}} will not support 
models with more coefficients than data. 

The situation in which factor variable random effects intentionally have unobserved levels 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.

The implementation is designed so that supplying random effect factor levels to \code{\link{predict.gam}} that were not levels of
the factor when fitting, will result in the corresponding random effect (or interactions involving it) being set to zero (with zero standard error) for prediction. See \code{\link{random.effects}} for an example. This is achieved by the \code{Predict.matrix} method zeroing any rows of the prediction matrix involving factors that are \code{NA}. \code{\link{predict.gam}} will set any factor observation to \code{NA} if it is a level not present in the fit data. 

}

\references{
Wood, S.N. (2008) Fast stable direct fitting and smoothness
 selection for generalized additive models. Journal of the Royal
 Statistical Society (B) 70(3):495-518
}

\author{ Simon N. Wood \email{simon.wood@r-project.org}}

\seealso{\code{\link{gam.vcomp}}, \code{\link{gamm}}}

\examples{
## see ?gam.vcomp

require(mgcv)
## simulate simple random effect example
set.seed(4)
nb <- 50; n <- 400
b <- rnorm(nb)*2 ## random effect
r <- sample(1:nb,n,replace=TRUE) ## r.e. levels
y <- 2 + b[r] + rnorm(n)
r <- factor(r)
## fit model....
b <- gam(y ~ s(r,bs="re"),method="REML")
gam.vcomp(b)

## example with supplied precision matrices...
b <- c(rnorm(nb/2)*2,rnorm(nb/2)*.5) ## random effect now with 2 variances
r <- sample(1:nb,n,replace=TRUE) ## r.e. levels
y <- 2 + b[r] + rnorm(n)
r <- factor(r)
## known precision matrix components...
S <- list(diag(rep(c(1,0),each=nb/2)),diag(rep(c(0,1),each=nb/2)))
b <- gam(y ~ s(r,bs="re",xt=list(S=S,rank=c(nb/2,nb/2))),method="REML")
gam.vcomp(b)
summary(b)
}
\keyword{models} \keyword{regression}%-- one or more ..