cond.expectation: Computes Conditional Expectation for Unobserved Locations

Description

Compute the conditional expectations (i.e. predictions) at the unobserved space-time locations.

Usage

cond.expectation(x, mesa.data.model, mesa.data = NA, Nmax = 1000,
    only.pars = FALSE, compute.beta = FALSE, no.nugget = FALSE,
    only.obs = FALSE, pred.var = TRUE, pred.covar = FALSE,
    combine.data = FALSE, type = "p")

Arguments

Model parameters for which to compute the conditional expectation.

mesa.data.model

Data structure holding observations and information regarding where to compute predictions. Predictions are computed conditional on the observations in mesa.data.model. See further

mesa.data

Data structure holding additional locations, geographic covariates, and spatio-temporal covariates, smooth temporal trends, and optionally additional observations. See further mesa.data

Nmax

Limits the size of matrices constructed when computing expectations. Use a smaller value if memory becomes a problem.

only.pars

Compute only the regression parameters (using GLS) along with the related variance.

only.pars=TRUE is incompatible with type="f"

compute.beta

Compute the latent beta-fields at all locations in mesa.data.model and/or mesa.data.

no.nugget

Set the nugget for the residual nu field to zero when computing the conditional expectations. Makes no difference at unobserved locations but implies smoothing at locations with observations.

only.obs

Compute predictions at only the observations specified in mesa.data. Used to limit computations when doing cross-validation.

This requires that mesa.data!=NA and that

pred.var

If TRUE, compute point-wise variances for the predictions.

pred.covar

If TRUE, compute covariance matrices for the predicted time series at each location. pred.covar = TRUE implies pred.var = TRUE and sets Nmax equal to the number of timepoints.

combine.data

If TRUE will combine locations and covariates from mesa.data.model and mesa.data, but will use temporal trends from

type

A single character indicating the type of log-likelihood to use. Valid options are "f", "p", and "r", for full, profile or restricted maximum likelihood (REML).

Value

The function returns a list containing:
parsA list with regression parameters and related variances. The variances are set to zero if type="f" and the input x contains the regression parameters.
pars contain gamma.E and alpha.E which are regression coefficients for the spatio-temporal model and land-use covaraiates respectively. The associated variances are found in gamma.V and alpha.V, with the cross-covariance between gamma and alpha in gamma.alpha.C
EX.beta.muA (number of locations) - by - (number of smooth trends) matrix with mean value component of the latent beta-fields (i.e. covariates * alpha). NULL if compute.beta=FALSE.
EX.betaA (number of locations) - by - (number of smooth trends) matrix with predictions of the latent beta-fields.
NULL if compute.beta=FALSE.
VX.betaA (number of locations) - by - (number of smooth trends) matrix or a (number of locations) - by - (number of locations) - by - (number of smooth trends) 3D-array containing the prediction (co)variances for the beta-fields.
NULL if compute.beta=FALSE or pred.var=FALSE, a matrix if pred.covar=FALSE, and a 3D-array if pred.covar=TRUE.
EX.muA (number of timepoints) - by - (number of locations) matrix with predictions based on only the regression parameters, geographic covariates, and temporal trends. I.e. only the deterministic part of the spatio-temporal model. Predictions are computed at the space-time locations defined by mesa.data.model and/or mesa.data. If only.obs=TRUE this is a vector with predictions at the space-time locations of the observations in mesa.data$obs.
EX.mu.betaPredictions based on the latent-beta fields, but not on the residual nu field; i.e. the predictions in EX.mu with kriging for the geographic covariates. Same size as EX.mu above; NULL if compute.beta=FALSE.
EXFull predictions at the space-time locations defined by mesa.data.model and/or mesa.data.
Same size as EX.mu above.
VXPointwise prediction variances for all locations in EX above. NULL if pred.var=FALSE.
VX.fullA list with (number of locations) elements, each element is a (number of timepoints) - by - (number of timepoints) temporal covariance matrix for the predicted timeseries at that location. NULL if pred.covar=FALSE.
IA vector with the locations of the observations in mesa.data.model or mesa.data. To extract predictions at these observations locations use EX[I]. NULL if only.obs=TRUE.

encoding

latin1

Details

Predictions are computed for the space-time locations in mesa.data.model and/or mesa.data, conditional on the observations in mesa.data.model and parameters given in x.

Depending on the the model requested in type the parameters specified in x can be either the log-covariances or regression parameters and log-covariances parameters.

If type="f" the user has the freedom to specify all model parameters -- including the regression coefficients -- rather than just the log-covariance parameters. However, if type="f" and only log-covariance parameters are given then the regression parameters are inferred using generalised least squares (GLS). If type!="f" the regression parameters are always computed using GLS (any input is ignored). The GLS can be obtained through cond.expectation(x, mesa.data.model, only.pars = TRUE, type="p").

The space-time locations at which to compute predictions are specified through mesa.data.model and mesa.data. If mesa.data=NA predictions will be computed at all space-time locations given in mesa.data.model. If mesa.data!=NA three options exist:

Ifcombine.data=FALSEpredictions will be given atonlythe space-time locations inmesa.data. Ifmesa.datadoes not have a smooth temporal trend, the trend in mesa.data.modelwill be used (with a warning). Predictions are still computed conditional on the observations inmesa.data.modeland any observations inmesa.dataareignored.
Ifcombine.data=FALSEandonly.obs=TRUEpredictions are computed as in (1) above, but only at space-time locations given inmesa.data$obs. This is primarily used to limit the number of computations needed when doing cross-validation.
Ifcombine.data=TRUElocations and covariates from bothmesa.data.modelandmesa.dataare combined usingcombineMesaData. The smooth temporal trends are taken frommesa.data.modelwith trends inmesa.databeingignored. I.e. predictions are computed conditional ononlythe observations inmesa.data.model, at locations given in bothmesa.data.modelandmesa.dataand atonlythe times specified inmesa.data.model$trend.

In addition to computing the conditional expectation at a number of space-time locations the function also computes predictions based on only the regression part of the model as well as the latent beta-fields (the latter only if compute.beta=TRUE).

NOTE: When computing the latent beta-fields, some matrix multiplications occur in a different order. this can lead to differences on the order of 1e-13 in some conditional expectations and variances.

Predictions variances can also be computed. If pred.var=TRUE point-wise variances for the predictions (and the latent beta-fields) are computed. If instead pred.covar=TRUE the full covariance matrices for each predicted time series is computed (as well as full covariance matrices for the latent beta-fields); this implies that the covariances between temporal predictions at the same location are calculated but not, due to memory restrictions, any covariances between locations.

Examples

Run this code

##load data
data(mesa.data)
data(mesa.data.model)
data(mesa.data.res)

##extract estimated parameters
x <- mesa.data.res$par.est$res.best$par

##find regression parameters using GLS
x.reg <- cond.expectation(x, mesa.data.model, only.pars = TRUE)
str(x.reg)

##compute predictions at all locations, including beta-fields
EX <- cond.expectation(x, mesa.data.model, compute.beta=TRUE)
##compute predictions at only observations locations
EX.obs <- cond.expectation(x, mesa.data.model, mesa.data=mesa.data,
                           only.obs=TRUE, compute.beta=FALSE)
##Let's load precomputed results instead.
EX <- mesa.data.res$EX
EX.obs <- mesa.data.res$EX.obs

##Compare the predictions at all locations and only obs
dim(EX$EX)
dim(EX.obs$EX)

##estimate beta from the observations for reference
##create data matrix
D <- create.data.matrix(mesa.data)
beta <- matrix(NA,dim(D)[2], dim(mesa.data$trend)[2])
##extact the temporal trends
F <- mesa.data$trend
##drop the date column
F$date <- NULL
##estimate the beta-coeficients at each location
for(i in 1:dim(D)[2])
  beta[i,] <- summary(lm(D[,i] ~ as.matrix(F)))$coefficients[,1]

##Study the results
##Start by comparing beta fields
par(mfcol=c(1,1), mar=c(4.5,4.5,2,.5), pty="s")
plot(x=beta[,1], y=EX$EX.beta[,1], main="Temporal Intercept",
     xlab="Empirical estimate", ylab="Spatio-Temporal Model")
abline(0,1,col="grey")

##plot predictions and observations for 4 locations
par(mfrow=c(4,1),mar=c(2.5,2.5,2,.5))
plotPrediction(EX,  1, mesa.data.model)
plotPrediction(EX, 10, mesa.data.model)
plotPrediction(EX, 17, mesa.data.model)
plotPrediction(EX, 22, mesa.data.model)

##compare the only obs predictions with what we can extract from EX
par(mfcol=c(1,1), mar=c(4.5,4.5,2,.5), pty="s")
plot(EX$EX[EX$I], EX.obs$EX)

##minor (1e-14) numerical differences in the results
print(range(EX$EX[EX$I]-EX.obs$EX))

Run the code above in your browser using DataLab