Revision 513ad2c4d40f9e25102134188eb154b18140f6a2 authored by Martin Schlather on 16 April 2017, 09:57:35 UTC, committed by cran-robot on 16 April 2017, 09:57:35 UTC
1 parent 1b4f8cd
RFformula.Rd
\name{RFformula}
\alias{RFformula}
\title{RFformula - syntax to design random field models with trend
or linear mixed models}
\description{It is described how to create a formula, which can
e.g. be used as an argument of \command{\link{RFsimulate}} and
\command{\link{RFfit}} to simulate and to fit data accordingly to the
model described by the formula.
In general, the created formula serves two purposes:
\itemize{
\item to describe models in the \dQuote{Linear Mixed
Models}-framework including fixed and random effects
\item to define models for random fields including trend surfaces from
a geostatistical point of view.
}
Thereby, fixed effects and trend surfaces are adressed via
the expression \command{\link{RMfixed}} and the function
\command{\link{RMtrend}}; the covariance structures of the zero-mean
multivariate normally distributed random effects and random field
components are adressed by objects of class \code{\link[=RMmodel-class]{RMmodel}}, which
allow for a very flexible covariance specification.
}
\details{
The formula should be of the type
\deqn{response ~ fixed effects + random effects + error term}
or
\deqn{response ~ trend + zero-mean random field + nugget effect,}
respectively.
Thereby:
\itemize{
\item response\cr
optional; name of response variable
\item fixed effects/trend:\cr
optional, should be a sum (using \command{\link[=RMplus]{+}})%%check link
of components either of the form \code{X@RMfixed(beta)} or
\code{\link{RMtrend}(...)} with \eqn{X} being a design matrix
and \eqn{\beta} being a vector of coefficients (see
\command{\link{RMfixed}} and \command{\link{RMtrend}}).\cr
Note that a fixed effect of the form \eqn{X} is interpreted as
\code{X@RMfixed(beta=NA)} by default (and \eqn{\beta} is estimated
provided that the formula is used in \command{\link{RFfit}}).
\item random effects/zero-mean random field:\cr
optional, should be a sum (using \command{\link[=RMplus]{+}})%%check link
of components of the form \code{Z@model}
where \eqn{Z} is a design matrix and \code{model} is an object of
class \code{\link[=RMmodel-class]{RMmodel}}.\cr
\code{Z@model} describes a vector of random effects which is
normally distributed with zero mean and covariance matrix \eqn{Z
\Sigma Z^T} where \eqn{Z^T} is the transpose of \eqn{Z} and
\eqn{\Sigma} is the covariance matrix according to \code{model}.\cr
Note that a random effect/random fluctuation of the form
\code{model} is viewed as \code{I@model} where \eqn{I} is the
identity matrix of corresponding dimension.
\item error term/nugget effect\cr
optional, should be of the form \code{\link{RMnugget}(...)}.
\command{\link{RMnugget}} describes a vector
of iid Gaussian random variables.
Please note that the character \dQuote{@} in the RFformula-context can only
be used to multiply design-matrices with corresponding vectors of
fixed or random effects, whereas in the context of S4-classes \dQuote{@} is
used to access slots of corresponding objects.
}
}
\section{IMPORTANT}{
Note that in formula constants are interpreted as part of a mixed
model, i.e. the corresponding parameter has to be estimated
(e.g. \code{~ 1 + ...}) whereas models not given as formula the
parameters to be estimated must be given explicitely.
}
\note{ %to do: make the following point cleares
(additional) argument names should always start with a capital
letter. Small initial letters are reserved for \command{\link{RFoptions}}.
}
\references{
\itemize{
\item Chiles, J.-P. and P. Delfiner (1999) \emph{Geostatistics. Modeling
Spatial Uncertainty.} New York, Chichester: John Wiley & Sons.
\item McCulloch, C. E., Searle, S. R. and Neuhaus, J. M. (2008)
\emph{Generalized, linear, and mixed models.} Hoboken, NJ: John
Wiley & Sons.
\item Ruppert, D. and Wand, M. P. and Carroll, R. J. (2003)
\emph{Semiparametric regression.} Cambridge: Cambridge University
Press.
}
}
\author{Martin Schlather, \email{schlather@math.uni-mannheim.de}
}
\seealso{
\command{\link{RMmodel}},
\command{\link{RFsimulate}},
\command{\link{RFfit}},
\code{\link[=RandomFields-package]{RandomFields}}.}
\keyword{spatial}
\examples{
RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
## RFoptions(seed=NA) to make them all random again
\dontshow{StartExample()}
RFoptions(modus_operandi="sloppy")
##############################################################
#
# Example : Simulation and fitting of a two-dimensional
# Gaussian random field with exponential covariance function
#
###############################################################
V <- 10
S <- 0.3
M <- 3
model <- RMexp(var=V, scale=S) + M
x <- y <- seq(1, 3, 0.1)
simulated <- RFsimulate(model = model, x=x, y=y)
plot(simulated)
# an alternative code to the above code:
simulated2 <- RFsimulate(model = ~ 1@RMfixed(beta=M) +
RMexp(var=V, scale=S),x=x, y=y, V=V, S=S, M=M)
plot(simulated2)
# Estimate parameters of underlying covariance function via
# maximum likelihood
model.na <- ~ 1@RMfixed(beta=NA) + RMexp(var=NA, scale=NA)
fitted <- RFfit(model=model.na, data=simulated)
# compare sample mean of data with ML estimate:
mean(simulated@data[,1])
fitted
\dontshow{\dontrun{
##############################################################
#
# Example 2: Fitting a standard linear mixed model using maximum
# likelihood to estimate the variance components
#
###############################################################
# Y = W*beta + Z*u + e
# where u ~ N(0, (sigma_u)^2 A) and e ~ N(0, (sigma_e)^2)
W <- rep(1,times=10)
Z <- matrix(rnorm(150) ,nrow=10, ncol=15)
A <- RFcovmatrix(0:14, model=RMexp())
response <- W*5 + Z%*%matrix(1:225, nrow=15)%*%rnorm(15, sd=10) + rnorm(10, sd=3)
# Estimate beta, (sigma_u)^2 and (sigma_e)^2:
model <- response ~ W@RMfixed(beta=NA) +
Z@RMfixcov(A, var=NA) +
RMnugget(var=NA)
fitted <- RFfit(model=model, data=response, W=W, Z=Z, A=A)
}}
\dontshow{\dontrun{
#### THIS EXAMPLE IS NOT PROGRAMMED YET
###########################################################
#
# Example 3: Simulate and fit a geostatistical model
#
###########################################################
# Simulate a Gaussian random field with trend m((x,y)) = 2 + 1.5 x - 3 y
# with vector of coordinates (x,y)
# and covariance function C(s,t) = 3*exp(-||(2*(s_1-t_1),s_2-t_2)||) +
# 1.5*exp(-||(2*(s_1-t_1),s_2-t_2)||^2)
# for s=(s_1,s_2) and t=(t_1,t_2)
model <- ~ RMtrend(mean=2) +
RMtrend(arbitraryfct=function(x) x, fctcoeff=1.5) +
RMtrend(arbitraryfct=function(y) y, fctcoeff=-3) +
RMplus(RMexp(var=3), RMgauss(var=1.5),
Aniso=matrix(c(2,0,0,1),nrow=2))
simulated <- RFsimulate(model=model,x=seq(0,2,0.1),y=seq(-1,3,0.2))
# equivalent model
model <- RMtrend(polydeg=1, polycoeff=c(2, 1.5, .3)) +
RMplus(RMexp(var=3), RMgauss(var=1.5),
Aniso=matrix(c(2,0,0,1), nrow=2))
simulated <- RFsimulate(model=model, x=seq(0,2,0.1), y=seq(-1,3,0.2))
# Estimate trend (polynomial of degree 1) and variance components such
# that var_exp = 2*var_gauss as in the true model used for simulation
model.na <- ~ RMtrend(polydeg=1) +
RMplus(RMexp(var=2),RMgauss(var=2),var=NA,
Aniso=matrix(c(NA,0,0,NA),nrow=2))
fit <- RFfit(model=model.na, data=simulated)
}}
\dontshow{\dontrun{
#### THIS EXAMPLE IS NOT PROGRAMMED YET
####################################################################
#
# Example 4: Simulate and fit a multivariate geostatistical model
#
####################################################################
# Simulate a bivariate Gaussian random field with trend
# m_1((x,y)) = x + 2*y and m_2((x,y)) = 3*x + 4*y
# where m_1 is a hyperplane describing the trend for the first response
# variable and m_2 is the trend for the second one;
# the covariance function is a multivariate nugget effect
x <- y <- 0:10
model <- ~ RMtrend(plane=matrix(c(1,2,3,4), ncol=2)) +
RMnugget(var=0.5, vdim=2)
multi.simulated <- RFsimulate(model=model, x=x, y=y)
# Fit the Gaussian random field with unknown trend for the second
# response variable and unknown variances
model.na <- ~ RMtrend(plane=matrix(c(NA,NA,NA,NA), ncol=2)) +
RMnugget(var=NA, vdim=2)
fit <- RFfit(model=model.na, data=multi.simulated)
}}
\dontshow{RFoptions(modus_operandi="normal")}
\dontshow{FinalizeExample()}
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...