https://github.com/cran/RandomFields
Tip revision: e994a4415e67fa60cbfd3f208aaab20872521c0b authored by Martin Schlather on 14 February 2019, 21:02:19 UTC
version 3.3
version 3.3
Tip revision: e994a44
RFformula.Rd
\name{RFformula}
\alias{RFformula}
\alias{RMformula}
\title{RFformula - syntax to design random field models with trend
or linear mixed models}
\description{It is described how to create a formula, which, for example, can be used as an argument of \command{\link{RFsimulate}} and
\command{\link{RFfit}} to simulate and to fit data according 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 can be adressed via
the expression \command{\link{RMfixed}} and the function
\command{\link{RMtrend}}. In simple cases, the trend can also
be given in a very simple, see the examples below.
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.
See \link{RFformulaAdvanced} for rather complicated model definitions.
}
\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 linear
model, i.e. the corresponding parameter has to be estimated
(e.g. \code{~ 1 + ...}) whereas in models not given as formula the
parameters to be estimated must be given explicitly.
}
\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.
}
}
\me
\seealso{
\command{\link{RMmodel}},
\command{\link{RFsimulate}},
\command{\link{RFfit}},
\code{\link[=RandomFields-package]{RandomFields}}.}
\keyword{spatial}
\examples{\dontshow{StartExample()}
RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
## RFoptions(seed=NA) to make them all random again
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:
model <- ~ Mean + RMexp(var=Var, scale=Sc)
simulated2 <- RFsimulate(model = model,x=x, y=y, Var=V, Sc=S, Mean=M)
plot(simulated2)
# a third way of specifying the model using the argument 'param'
# the initials of the variables do not be captical letters
model <- ~ M + RMexp(var=var, scale=sc)
simulated3 <- RFsimulate(model = model,x=x, y=y,
param=list(var=V, sc=S, M=M))
plot(simulated3)
# Estimate parameters of underlying covariance function via
# maximum likelihood
model.na <- ~ NA + RMexp(var=NA, scale=NA)
fitted <- RFfit(model=model.na, data=simulated)
# compare sample mean of data with ML estimate, which is very similar:
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()}}