RFformula: 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 RFsimulate and RFfit to simulate and to fit data accordingly to the model described by the formula.

In general, the created formula serves two purposes:

to describe models in theLinear Mixed Models-framework including fixed and random effects
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 RMfixed and the function RMtrend; the covariance structures of the zero-mean multivariate normally distributed random effects and random field components are adressed by objects of class RMmodel, which allow for a very flexible covariance specification.

Arguments

Details

The formula should be of the type $$response ~ fixed effects + random effects + error term$$ or $$response ~ trend + zero-mean random field + nugget effect,$$ respectively. Thereby:

response optional; name of response variable
fixed effects/trend: optional, should be a sum (using+)%%check link of components either of the formX@RMfixed(beta)orRMtrend(...)with$X$being a design matrix and$\beta$being a vector of coefficients (seeRMfixedandRMtrend). Note that a fixed effect of the form$X$is interpreted asX@RMfixed(beta=NA)by default (and$\beta$is estimated provided that the formula is used inRFfit).
random effects/zero-mean random field: optional, should be a sum (using+)%%check link of components of the formZ@modelwhere$Z$is a design matrix andmodelis an object of classRMmodel. Z@modeldescribes a vector of random effects which is normally distributed with zero mean and covariance matrix$Z \Sigma Z^T$where$Z^T$is the transpose of$Z$and$\Sigma$is the covariance matrix according tomodel. Note that a random effect/random fluctuation of the formmodelis viewed asI@modelwhere$I$is the identity matrix of corresponding dimension.
error term/nugget effect optional, should be of the formRMnugget(...).RMnuggetdescribes a vector of iid Gaussian random variables.
Please note that the character@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@is used to access slots of corresponding objects.

References

Chiles, J.-P. and P. Delfiner (1999)Geostatistics. Modeling Spatial Uncertainty.New York, Chichester: John Wiley & Sons.
McCulloch, C. E., Searle, S. R. and Neuhaus, J. M. (2008)Generalized, linear, and mixed models.Hoboken, NJ: John Wiley & Sons.
Ruppert, D. and Wand, M. P. and Carroll, R. J. (2003)Semiparametric regression.Cambridge: Cambridge University Press.

Examples

Run this code

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
model <- RMexp(var=V, scale=S)
x <- y <- seq(1, 3, by= if (interactive()) 0.1 else 1)

simulated <- RFsimulate(model = model, x=x, y=y)
# add overall mean/trend of 3 to the response/simulated data 
simulated@data <- simulated@data + 3
#returns an object of class RFspatialPointsDataFrame including
#coordinates and the simulated response vector
plot(simulated)


# an alternative code to the above code:
simulated2 <- RFsimulate(model = ~ 1@RMfixed(beta=3) +
               RMexp(var=V, scale=S),x=x, y=y, V=V, S=S)
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@ml@trend


\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
# 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)
}


\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)
}



\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)
}



RFoptions(modus_operandi="normal")

FinalizeExample()

Run the code above in your browser using DataCamp Workspace

Description

Arguments

Details

References

See Also

Examples