predict.STmodel: Computes Conditional Expectation at Unobserved Locations

Description

Compute the conditional expectations (i.e. predictions) at the unobserved space-time locations. Predictions are computed for the space-time locations in object and/or STdata, conditional on the observations in object and parameters given in x.

Usage

## S3 method for class 'STmodel':
predict(object, x, STdata = NULL,
    Nmax = 1000, only.pars = FALSE, nugget.unobs = 0,
    only.obs = FALSE, pred.var = TRUE, pred.covar = FALSE,
    beta.covar = FALSE, combine.data = FALSE, type = "p",
    ...)

Arguments

object

STmodel object for which to compute predictions.

Model parameters for which to compute the conditional expectation. Either as a vector/matrix or an estimateSTmodel from estimate.STmodel.

STdata

STdata/STmodel object with locations/times for which to predict. If not given predictions are computed for locations/times in object

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.

nugget.unobs

Value of nugget at unonserved locations, either a scalar or a vector with one element per unobserved site. NOTE: All sites in STdata are considered unobserved!

only.obs

Compute predictions at only locations specified by observations in STdata. Used to limit computations when doing cross-validation. only.obs=TRUE implies pred.covar=FALSE and combine.data=FALSE<

pred.var,pred.covar

Compute point-wise prediction variances; or 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 ti

beta.covar

Compute the full covariance matrix for the latent beta-fields, otherwise only the diagonal elements of V(beta|obs) are computed.

combine.data

Combine object and STdata and predict for the joint set of points, see c.STmodel.

type

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

...

Ignored additional arguments.

Value

The function returns a list containing (objects not computed will be missing):
optsCopy of options used in the function call.
parsA list with regression parameters and related variances. pars contain gamma.E and alpha.E with regression coefficients for the spatio-temporal model and land-use covaraiates; variances are found in gamma.V and alpha.V; cross-covariance between gamma and alpha in gamma.alpha.C.
betaA list with estimates of the beta-fields, including the regression mean mu, conditional expectations EX, possibly variances VX, and the full covariance matrix VX.full.
EX.mupredictions based on the regression parameters, geographic covariates, and temporal trends. I.e. only the deterministic part of the spatio-temporal model.
EX.mu.betaPredictions based on the latent-beta fields, but excluding the residual nu field.
EXFull predictions at the space-time locations in object and/or STdata.
VXPointwise variances for all locations in EX.
VX.predPointwise prediction variances for all locations in EX, i.e. including contribution from nugget.unobs.
VX.fullA list with (number of locations) elements, each element is a (number of timepoints) - by - (number of timepoints) temporal covariance matrix for the timeseries at each location.
IA vector with the locations of the observations in object or STdata. To extract predictions at the observations locations use EX[I].

Details

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.

Prediction are computed as the conditional expectation of a latent field given observations. This implies that E(X_i| Y_i) != Y_i, with the difference being due to smoothing over the nugget. Further two possible variance are given, V(X_i|Y_i) and V(X_i|Y_i)+nugget_i. Here the nugget for unobserved locations needs to be specified as an additional argument nugget.nobs. The two variances correspond, losely, to confidence and prediction intervals.

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; this implies that the covariances between temporal predictions at the same location are calculated but not, due to memory restrictions, any covariances between locations. beta.covar=TRUE gives the full covariance matrices for the latent beta-fields.

Examples

Run this code

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

##find regression parameters using GLS
x.reg <- predict(mesa.model, est.mesa.model, only.pars = TRUE)
str(x.reg$pars)

##Split data into FIXED and AQS
I.aqs <- mesa.model$locations$ID[mesa.model$locations$type=="AQS"]
I.aqs <- mesa.model$obs$ID %in% I.aqs
mesa.model.aqs <- dropObservations(mesa.model, !I.aqs)
mesa.model.fixed <- dropObservations(mesa.model, I.aqs)

##compute predictions at all locations, including beta-fields
  pred.mesa.model <- predict(mesa.model, est.mesa.model,
                             pred.var=TRUE)
  ##predict at FIXED using AQS-sites
  pred.mesa.model.obs <- predict(mesa.model.aqs, est.mesa.model,
                                 STdata=mesa.model.fixed, only.obs=TRUE)
##Let's load precomputed results instead.
data(pred.mesa.model)

##Compare the predictions at all locations and only obs
print(pred.mesa.model)
print(pred.mesa.model.obs)

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

##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=pred.mesa.model$beta$EX[,1],
     main="Temporal Intercept",
     xlab="Empirical estimate", ylab="Spatio-Temporal Model")
abline(0,1,col="grey")
##or just the regression part of the beta fields
points(x=beta[,1], y=pred.mesa.model$beta$mu[,1], col=2)

##plot predictions and observations for 4 locations
par(mfrow=c(4,1),mar=c(2.5,2.5,2,.5))
plot(pred.mesa.model, ID=1, STmodel=mesa.model, pred.var=TRUE)
plot(pred.mesa.model, ID=10, STmodel=mesa.model, pred.var=TRUE)
plot(pred.mesa.model, ID=17, STmodel=mesa.model, pred.var=TRUE)
plot(pred.mesa.model, ID=22, STmodel=mesa.model, pred.var=TRUE)

Run the code above in your browser using DataLab