Nothing
R/STmodel_predict.R
In SpatioTemporal: Spatio-Temporal Model Estimation
############################################
## S3-METHOD THAT PREDICTS FROM A STmodel ##
############################################
##Functions in this file:
## predict.STmodel EX:ok
## print.predictSTmode EX:ok
## plot.predictSTmodel EX:ok
##' Compute the conditional expectations (i.e. predictions) at the unobserved
##' space-time locations. Predictions are computed for the space-time locations in
##' \code{object} and/or \code{STdata}, conditional on the observations (and
##' temporal trends) in \code{object} and parameters given in \code{x}.
##'
##' 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 \code{E(X_i| Y_i) != Y_i}, with the
##' difference being due to smoothing over the nugget. Further two possible
##' variance can be computed (see below), \code{V(X_i|Y_i)} and
##' \code{V(X_i|Y_i)+nugget_i}. Here the nugget for unobserved locations needs
##' to be specified as an additional argument \code{nugget.nobs}. The two
##' variances correspond, losely, to confidence and prediction intervals.
##'
##' Variances are computed if \code{pred.var=TRUE} point-wise variances for the
##' predictions (and the latent beta-fields) are
##' computed. If instead \code{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 \emph{not}, due
##' to memory restrictions, any covariances between locations.
##' \code{beta.covar=TRUE} gives the full covariance matrices for the latent
##' beta-fields.
##'
##' If \code{transform!="none"} the field is assumed to be log-Gaussian and
##' expectations are transformed, and if \code{pred.var=TRUE} the mean squared
##' prediction errors are given.
##'
##' @title Computes Conditional Expectation at Unobserved Locations
##'
##' @param object \code{STmodel} object for which to compute predictions.
##' @param x Model parameters for which to compute the conditional
##' expectation. Either as a vector/matrix or an \code{estimateSTmodel} from
##' \code{\link{estimate.STmodel}}.
##' @param STdata \code{STdata}/\code{STmodel} object with locations/times for
##' which to predict. If not given predictions are computed for locations/times
##' in \code{object}
##' @param Nmax Limits the size of matrices constructed when computing
##' expectations. Use a smaller value if memory becomes a problem.
##' @param only.pars Compute only the regression parameters (using GLS) along
##' with the related variance.
##' @param nugget.unobs Value of nugget at unonserved locations, either a scalar
##' or a vector with one element per unobserved site. \strong{NOTE:} All sites in
##' \code{STdata} are considered unobserved!
##' @param only.obs Compute predictions at only locations specified by
##' observations in \code{STdata}. Used to limit computations when doing
##' cross-validation.
##' \code{only.obs=TRUE} \emph{implies} \code{pred.covar=FALSE} and
##' \code{combine.data=FALSE}.
##' Further \code{\link{createSTmodel}} will be called on any \code{STdata}
##' input, possibly \emph{reordering the observations.}
##' @param pred.var,pred.covar Compute point-wise prediction variances; or
##' compute covariance matrices for the predicted time series at each location.
##' \code{pred.covar=TRUE} \emph{implies} \code{pred.var=TRUE} and sets
##' \code{Nmax} equal to the number of timepoints.
##' @param beta.covar Compute the full covariance matrix for the latent
##' beta-fields, otherwise only the diagonal elements of V(beta|obs) are
##' computed.
##' @param combine.data Combine \code{object} and \code{STdata} and predict for
##' the joint set of points, see \code{\link{c.STmodel}}.
##' @param type A single character indicating the type of prediction to
##' compute. Valid options are "f", "p", and "r", for \emph{full},
##' \emph{profile} or \emph{restricted maximum likelihood} (REML). For profile
##' and full the predictions are computed assuming that \emph{both} covariance
##' parameters and regression parameters are known,
##' e.g. \code{E(X|Y,cov_par,reg_par)}; for REML predictions are compute
##' assuming \emph{only} covariance parameters known,
##' e.g. \code{E(X|Y,cov_par)}. The main difference is that REML will have
##' \emph{larger} variances due to the additional uncertainty in the
##' regression parameters.
##' @param transform Regard field as log-Gaussian and apply exponential
##' transformation to predictions. For the final expectations two options
##' exist, either a unbiased prediction or the (biased) mean-squared error
##' predictions.
##' @param LTA Compute long-term temporal averages. Either a logical value or a
##' list; if \code{TRUE} then averages at each location (and variances if
##' \code{pred.var=TRUE}) are computed; otherwise this should be a list with
##' elements named after locations and each element containing a vector (or
##' list of vectors) with dates over which to compute averages. If
##' \code{only.obs=TRUE} averages are computed over only the observations.
##' @param ... Ignored additional arguments.
##'
##' @return The function returns a list containing (objects not computed
##' will be missing):
##' \item{opts}{Copy of options used in the function call.}
##' \item{pars}{A list with regression parameters and related variances.
##' \code{pars} contain \code{gamma.E} and \code{alpha.E} with
##' regression coefficients for the spatio-temporal model and
##' land-use covaraiates; variances are found in \code{gamma.V}
##' and \code{alpha.V}; cross-covariance between gamma and alpha in
##' \code{gamma.alpha.C}.}
##' \item{beta}{A list with estimates of the beta-fields, including the
##' regression mean \code{mu}, conditional expectations \code{EX},
##' possibly variances \code{VX}, and the full covariance matrix
##' \code{VX.full}.}
##' \item{EX.mu}{predictions based on the regression parameters, geographic
##' covariates, and temporal trends. I.e. only the deterministic
##' part of the spatio-temporal model.}
##' \item{EX.mu.beta}{Predictions based on the latent-beta fields, but excluding
##' the residual nu field.}
##' \item{EX}{Full predictions at the space-time locations in
##' \code{object} and/or \code{STdata}.}
##' \item{EX.pred}{Only for \code{transform!="none"}, full predictions
##' including bias correction for prediction error.}
##' \item{VX,VX.pred}{Pointwise variances and prediction variances (i.e. incl.
##' contribution from \code{nugget.unobs}) for all locations in \code{EX}.}
##' \item{VX.full}{A 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.}
##' \item{MSPE,MSPE.pred}{Pointwise mean-square prediction errors for the
##' log-Gaussian fields.}
##' \item{log.EX,log.VX.pred,log.VX}{Pointwise predictions and variances for
##' the un-transformed fields when \code{transform!="none"}}
##' \item{LTA}{A data.frame with temporal averages for locations specified by
##' \code{LTA}. }
##' \item{I}{A vector with the locations of the observations in \code{object} or
##' \code{STdata}. To extract predictions at the observations locations use
##' \code{EX[I]}.}
##'
##' @example Rd_examples/Ex_predict_STmodel.R
##'
##' @author Johan Lindstrom
##'
##' @family STmodel methods
##' @family predictSTmodel methods
##' @importFrom stats predict
##' @method predict STmodel
##' @export
predict.STmodel <- function(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", LTA=FALSE,
transform=c("none","unbiased","mspe"), ...){
##################################
### INITIAL SETUP AND CHECKING ###
##check class belongings
stCheckClass(object, "STmodel", name="object")
if( !is.null(STdata) ){
stCheckClass(STdata, c("STdata","STmodel"), name="STdata")
}
##first ensure that type is lower case
type <- tolower(type)
##check if type is valid
stCheckType(type)
##args for transform
transform <- match.arg(transform)
if( transform=="none" ){ transform <- NULL }
if( inherits(x,"estimateSTmodel") ){
x <- coef(x,"all")$par
}##if( inherits(x,"estimateSTmodel") )
##figure out a bunch of dimensions
##dimensions$L!=0 is used to check for spatio-temporal covariates
dimensions <- loglikeSTdim(object)
##check length of input parameters
if( type=="f" && length(x)!=dimensions$nparam ){
stop( paste("type=f, requires", dimensions$nparam,
"parameters but length(x) =", length(x)) )
}
if( type!="f" && length(x)==dimensions$nparam ){
##drop the regression parameters
x <- x[(dimensions$nparam-dimensions$nparam.cov+1):dimensions$nparam]
}
if( type!="f" && length(x)!=dimensions$nparam.cov ){
stop( paste("type!=f, requires", dimensions$nparam.cov,
"parameters but length(x) =", length(x)) )
}
##only.pars is strange when type=="f"
if( only.pars && type=="f"){
warning("only.pars=TRUE and type=(f)ull only returns KNOWN parameters.",
immediate. = TRUE)
}
##only.pars, and we can ignore a bunch of things
if( only.pars ){
only.obs <- pred.covar <- beta.covar <- combine.data <- LTA <-FALSE
transform <- NULL
}else{
##check if only.obs is valid
if( only.obs && is.null(STdata) ){
stop("only.obs=TRUE requires STdata.")
}
##do we have an object to combine with
if( combine.data && is.null(STdata) ){
warning("No data to combine with; predicting for 'object'", immediate. = TRUE)
combine.data <- FALSE
}
##can't combine data if we're predicting at only obs.
if(only.obs && combine.data){
warning("only.obs=TRUE implies combine.data=FALSE.", immediate. = TRUE)
combine.data <- FALSE
}
##computing prediction covariates require !only.obs and pred.var
if(pred.covar && !pred.var){
warning("pred.covar=TRUE implies pred.var=TRUE.", immediate. = TRUE)
pred.var <- TRUE
}
if(pred.covar && only.obs){
warning("only.obs=TRUE implies pred.covar=FALSE.", immediate. = TRUE)
pred.covar <- FALSE
}
##computing beta covariances requires pred.var
if(beta.covar && !pred.var){
warning("beta.covar=TRUE implies pred.var=TRUE.", immediate. = TRUE)
pred.var <- TRUE
}
if( !is.null(transform) ){
if( type!="r" && transform!="unbiased" ){
warning("transform 'unbiased' and 'mspe' are equivalent when type!='r'",
immediate. = TRUE)
transform <- "unbiased"
}
if(pred.covar){
warning("transform implies pred.covar=FALSE.", immediate. = TRUE)
pred.covar <- FALSE
}
}##if( !is.null(transform) )
}##if( only.pars ){...}else{...}
##create the STdata used for predictions
if( is.null(STdata) ){
##pure copy, predict in the data set
STdata <- object
}else if(combine.data){
##combine the two datasets, and use for predictions (expands object with
##ONLY locations from STdata)
STdata <- c(object, STdata)
}else{
if( is.null(object$trend.fnc) && dim(object$trend)[2]!=1 ){
##trend.fnc missing and trend is more than a constant
##(backwards compability)
object$trend.fnc <- internalCreateTrendFnc(object$trend)
}
if( !inherits(STdata,"STmodel") ){
##predict only at STdata, and STdata not of class STmodel: need to cast.
##First drop ST covariate if no covariate in object
if(dimensions$L==0){
STdata$SpatioTemporal <- NULL
}
##and fix the trend (right no of trends, names and dates)
suppressMessages(STdata <- updateTrend(STdata, fnc=object$trend.fnc,
extra.dates=STdata$trend$date))
##since we're not using covariances for the prediction locations
##(specified seperately in nugget.unobs), we just pick a simple covariance
##structure, avoiding problems with missing covariates/levels.
cov.nu <- object$cov.nu
cov.nu$nugget <- TRUE
##Create an STmodel from STdata
STdata <- createSTmodel(STdata, LUR=object$LUR.list, ST=object$ST.list,
cov.beta=object$cov.beta, cov.nu=cov.nu,
locations=object$locations.list,
scale=!is.null(object$scale.covars),
scale.covars=object$scale.covars)
}else{
##STdata is an STmodel object ->
##test for consistent covariates and scaling (only equal scaling allowed).
areSTmodelsConsistent(object, STdata, "STdata")
}##if( !inherits(STdata,"STmodel") ){...}{else{...}
##STdata is now of type STmodel, lets fix the trend
suppressMessages(STdata <- updateTrend.STmodel(STdata, fnc=object$trend.fnc))
}##if( is.null(STdata) ){...}else if(...){...}else{...}
##Should we compute LTAs?
if( !is.list(LTA) && LTA ){
##yes, create list
if( only.obs ){
##averages over observtions only,
##i.e. split observation dates by their locations.
LTA.list <- split(STdata$obs$date, STdata$obs$ID)
}else{
##averages over all predictions
LTA.list <- lapply(STdata$locations$ID,
function(x){ return(STdata$trend$date) })
names(LTA.list) <- STdata$locations$ID
}
}else if( is.list(LTA) ){
##yes, list given
LTA.list <- LTA
LTA <- TRUE
}else{
##no, empty list and set LTA=FALSE
LTA.list <- NULL
LTA <- FALSE
}##if( !is.list(LTA) && LTA ){...}else if(...){...}else{...}
##check the validity of LTA.list
if( LTA ){
##check validity of sites
LTA.missing <- !(names(LTA.list) %in% STdata$locations$ID)
if( any(LTA.missing) ){
warning("Removed ", sum(LTA.missing),
" elements from LTA not found in STdata$locations$ID",
immediate. = TRUE)
LTA.list <- LTA.list[ !LTA.missing ]
}
if( length(LTA.list)!=0 ){
##convert components to lists
LTA.list <- lapply(LTA.list,
function(x){ switch(is.list(x)+1, list(x), x) })
##check validity of dates
LTA.tot <- 0
for(i in 1:length(LTA.list)){
for(j in length(LTA.list[[i]])){
LTA.missing <- !(LTA.list[[i]][[j]] %in% STdata$trend$date)
LTA.list[[i]][[j]] <- LTA.list[[i]][[j]][ !LTA.missing ]
LTA.tot <- LTA.tot + sum(LTA.missing)
}
##drop empty points
LTA.list[[i]] <- LTA.list[[i]][ sapply(LTA.list[[i]], length)!=0 ]
}##for(i in 1:length(LTA.list))
if( LTA.tot!=0 ){
warning("Removed ", LTA.tot,
" dates not found in STdata$trend$date from all LTA elements",
immediate. = TRUE)
}
##drop empty points
LTA.list <- LTA.list[ sapply(LTA.list, length)!=0 ]
}##if( length(LTA.list)!=0 )
##nothing left -> don't compute LTA's
if( length(LTA.list)==0 ){
LTA <- FALSE
LTA.list <- NULL
}
}##if( LTA )
#########################################
### COMPUTE COVARIANCE MATRICES FOR Y ###
##extract parameters from x
tmp <- loglikeSTgetPars(x, object)
if(type=="f"){
gamma <- tmp$gamma
alpha <- tmp$alpha
}
cov.pars.beta <- tmp$cov.beta
cov.pars.nu <- tmp$cov.nu
##the following only needed for type!="f", or predictions
if(type!="f" || !only.pars){
##nugget for unobserved sites
nugget.unobs <- internalFixNuggetUnobs(nugget.unobs, STdata,
cov.pars.nu$nugget)
##extract sparse matrices
Fobs <- expandF(object$F, object$obs$idx, n.loc=dimensions$n.obs)
##Create the Xtilde = [M FX] matrix
Xtilde <- as.matrix( Fobs %*% bdiag(object$LUR) )
##Add the spatio-temporal covariate (if it exists)
if( dimensions$L!=0 ){
Xtilde <- cbind(object$ST, Xtilde)
}
##create covariance matrices, beta-field
i.sigma.B <- makeSigmaB(cov.pars.beta$pars, dist = object$D.beta,
type = object$cov.beta$covf,
nugget = cov.pars.beta$nugget)
##and nu-field
i.sigma.nu <- makeSigmaNu(cov.pars.nu$pars, dist = object$D.nu,
type = object$cov.nu$covf,
nugget = cov.pars.nu$nugget,
random.effect = cov.pars.nu$random.effect,
blocks1 = object$nt, ind1 = object$obs$idx,
sparse=TRUE)
##calculate block cholesky factor of the matrices, in-place to conserve memory
i.sigma.B <- makeCholBlock(i.sigma.B, n.blocks=dimensions$m)
i.sigma.nu <- chol(i.sigma.nu)
##invert the matrices, in-place to conserve memory
i.sigma.B <- invCholBlock(i.sigma.B, n.blocks=dimensions$m)
i.sigma.nu <- chol2inv(i.sigma.nu)
##F'*inv(sigma.nu)
tF.iS <- t(Fobs) %*% i.sigma.nu
##F'*inv(sigma.nu)*F
tF.iS.F <- as.matrix(tF.iS %*% Fobs)
##inv(sigma.B|Y) = F'*inv(sigma.nu)*F + inv(sigma.B)
R.i.sigma.B.Y <- tF.iS.F + i.sigma.B
##calculate cholesky factor of inv(sigma.B|Y)
R.i.sigma.B.Y <- chol(R.i.sigma.B.Y)
##compute iSigma.nu*Xtilde and iSigma.nu*Xtilde
iS.X <- i.sigma.nu %*% Xtilde
iS.Y <- i.sigma.nu %*% object$obs$obs
##compute F'*iSigma.nu*Xtilde and F'*iSigma.nu*Y
tF.iS.X <- as.matrix(t(Fobs) %*% iS.X)
tF.iS.Y <- as.matrix(t(Fobs) %*% iS.Y)
##compute inv(sigma.B|Y)*F'*iSigma.nu*Xtilde
iSBY.tF.iS.X <- solveTriBlock(R.i.sigma.B.Y, tF.iS.X, transpose=TRUE)
iSBY.tF.iS.X <- solveTriBlock(R.i.sigma.B.Y, iSBY.tF.iS.X, transpose=FALSE)
}#if(type!="f" || !only.pars)
###############################################
### REGRESISON PARAMETERS (alpha and gamma) ###
if(type=="f"){
##regression parameters given
gamma.E <- c(tmp$gamma)
alpha.E <- unlist(tmp$alpha)
##known parameters, set variances to zero
gamma.V <- matrix(0,length(gamma.E),length(gamma.E))
alpha.V <- matrix(0,length(alpha.E),length(alpha.E))
gamma.alpha.C <- matrix(0,length(gamma.E),length(alpha.E))
##create a combined vector
gamma.alpha <- as.matrix(c(gamma.E,alpha.E))
}else{
##compute inv(Xtilde'*Sigma^-1*Xtilde) (and force symmetric)
i.XSX <- as.matrix( t(Xtilde) %*% iS.X )
i.XSX <- symmpart(i.XSX - t(tF.iS.X) %*% iSBY.tF.iS.X)
##compute t(Xtilde'*Sigma^-1*Y)
tmp2 <- as.matrix( object$obs$obs %*% iS.X )
tmp2 <- tmp2 - t(tF.iS.Y) %*% iSBY.tF.iS.X
##compute alpha and gamma
gamma.alpha <- solve(i.XSX, t(tmp2))
##extract computed parameters
if(dimensions$L!=0){
gamma.E <- c(gamma.alpha[1:dimensions$L])
alpha.E <- c(gamma.alpha[-(1:dimensions$L)])
}else{
gamma.E <- double(0)
alpha.E <- c(gamma.alpha)
}
##compute the variances
i.XSX <- solve(i.XSX)
if(dimensions$L!=0){
gamma.V <- i.XSX[1:dimensions$L,1:dimensions$L,drop=FALSE]
alpha.V <- i.XSX[-c(1:dimensions$L),-c(1:dimensions$L),drop=FALSE]
gamma.alpha.C <- i.XSX[1:dimensions$L,-c(1:dimensions$L),drop=FALSE]
}else{
gamma.V <- matrix(0,0,0)
alpha.V <- i.XSX
gamma.alpha.C <- matrix(0,length(gamma.E),length(alpha.E))
}
}##if(type=="f"){...}else{...}
##force to matrices
gamma.E <- as.matrix(gamma.E)
alpha.E <- as.matrix(alpha.E)
##add names
names.tmp <- loglikeSTnames(object, TRUE)
names.tmp <- names.tmp[1:(dimensions$nparam-dimensions$nparam.cov)]
if(dimensions$L!=0){
##first for gamma
rownames(gamma.E) <- names.tmp[1:dimensions$L]
colnames(gamma.V) <- rownames(gamma.V) <- rownames(gamma.E)
rownames(gamma.alpha.C) <- rownames(gamma.E)
rownames(alpha.E) <- names.tmp[-(1:dimensions$L)]
}else{
rownames(alpha.E) <- names.tmp
}
##and then for alpha
colnames(alpha.V) <- rownames(alpha.V) <- rownames(alpha.E)
colnames(gamma.alpha.C) <- rownames(alpha.E)
###############################
### CONSTRUCT RETURN OBJECT ###
out <- list()
##set class
class(out) <- "predictSTmodel"
##save options used for predict.
out$opts <- list(only.pars=only.pars, nugget.unobs=nugget.unobs,
only.obs=only.obs, pred.var=pred.var, pred.covar=pred.covar,
beta.covar=beta.covar, combine.data=combine.data, type=type,
transform=transform, LTA=LTA, LTA.list=LTA.list)
##collect the regression results
out$pars <- list(gamma.E=gamma.E, alpha.E=alpha.E,
gamma.V=gamma.V, alpha.V=alpha.V,
gamma.alpha.C=gamma.alpha.C)
if( out$opts$only.pars ){
##return the parameters
return( out )
}
##remove variables not needed, regression parameters...
rm(gamma.E, alpha.E, gamma.V, alpha.V, gamma.alpha.C, names.tmp)
##and options stored elsewhere
rm(only.pars, nugget.unobs, only.obs, pred.var, pred.covar,
beta.covar, combine.data, type, LTA, LTA.list, transform)
####################################
### COMPUTE inv(Sigma_oo)*(Y-mu) ###
##compute F'*iSigma.nu*(C-mu) and iSigma.nu*(C-mu)
tF.iS.C <- tF.iS.Y - tF.iS.X %*% gamma.alpha
iS.C <- iS.Y - iS.X %*% gamma.alpha
##compute inv(sigma.B|Y)*F'*iSigma.nu*(C-mu)
iSBY.tF.iS.C <- solveTriBlock(R.i.sigma.B.Y, tF.iS.C, transpose=TRUE)
iSBY.tF.iS.C <- solveTriBlock(R.i.sigma.B.Y, iSBY.tF.iS.C, transpose=FALSE)
##compute iSigma.tilde*(C-mu)
##this is the same for all the conditionals...
iSoo.C <- as.matrix( iS.C - t( tF.iS ) %*% iSBY.tF.iS.C )
##for REML we also need iSigma.tilde*Xtilde
if( out$opts$type=="r" ){
iSoo.Xtilde <- as.matrix( iS.X - t(tF.iS) %*% iSBY.tF.iS.X )
}
##remove variables not needed
rm(iS.X, iS.Y, tF.iS.X, tF.iS.Y, iSBY.tF.iS.X, tF.iS.C, iS.C, iSBY.tF.iS.C)
###################################
### DEFINE PREDICTION LOCAITONS ###
##create matrices with site index and temporal trends for all
##OBS date has to be INCREASING, i.e. same sorting as for STdata$obs...
if( out$opts$only.obs ){
##unobserved points
idx.unobs <- STdata$obs$idx
##time of observation for each site
T1 <- STdata$obs$date
##trend functions for all points matrix
Funobs <- STdata$F
if( dimensions$L!=0 ){
##extract the spatio-temporal trends
ST.unobs <- STdata$ST
}
##size of the unobserved locations
N.unobs <- length( unique(idx.unobs) )
T.unobs <- length( unique(T1) )
}else{
##size of the unobserved locations
N.unobs <- dim(STdata$locations)[1]
T.unobs <- dim(STdata$trend)[1]
##unobserved points
idx.unobs <- rep(1:N.unobs, each=T.unobs)
##time of observation for each site
T1 <- rep(STdata$trend$date, N.unobs)
##trend functions for all points matrix
Funobs <- matrix(1, (T.unobs*N.unobs), dimensions$m)
for(i in (1:dimensions$m)){
##if not constant, find relevant row in $trend
if(colnames(STdata$F)[i] != "const"){
Funobs[,i] <- rep(STdata$trend[,colnames(STdata$F)[i]], N.unobs)
}
}##for(i in (1:dimensions$m))
if( dimensions$L!=0 ){
##extract the spatio-temporal trends
ST.unobs <- matrix(STdata$ST.all, (T.unobs*N.unobs), dimensions$L)
}##if( dimensions$L!=0 )
}##if( out$opts$only.obs ) ... else ...
##compute number of observations for each time-point,
##taking into account that there might be missing dates.
date.all <- sort(unique( c(object$trend$date, STdata$trend$date) ))
nt.unobs <- nt.obs <- double(length(date.all))
for(i in c(1:length(date.all))){
nt.obs[i] <- sum( object$obs$date==date.all[i] )
}
##extract sparse matrices
Funobs <- expandF(Funobs, idx.unobs, n.loc=N.unobs)
##compute LUR times the temporal trends [M F*X] for unobserved
Xtilde.unobs <- as.matrix( Funobs %*% bdiag(STdata$LUR.all) )
if( dimensions$L!=0 ){
Xtilde.unobs <- cbind(ST.unobs, Xtilde.unobs)
}
###########################
### COMPUTE BETA FIELDS ###
##compute LUR*alpha (mean values for the beta fields)
out$beta <- list()
out$beta$mu <- matrix(bdiag(STdata$LUR.all) %*% out$pars$alpha.E,
ncol=length(STdata$LUR.all))
colnames(out$beta$mu) <- names(STdata$LUR.all)
rownames(out$beta$mu) <- rownames( STdata$LUR.all[[1]] )
##create distance matrix between unobserved and observed elements
##unobserved locations
loc.unobs.nu <- STdata$locations[, c("x.nu","y.nu"), drop=FALSE]
loc.unobs.beta <- STdata$locations[, c("x.beta","y.beta"), drop=FALSE]
##observed locations
I.obs <- match( colnames(object$D.nu), object$locations$ID)
loc.obs.nu <- object$locations[I.obs, c("x.nu","y.nu"), drop=FALSE]
loc.obs.beta <- object$locations[I.obs, c("x.beta","y.beta"), drop=FALSE]
##We also have that in the stripped data the indecies may not match so we need a
##second vector giving which idx in striped matches original idx.
Ind.2.1 <- match(object$locations$ID[I.obs], STdata$locations$ID, nomatch=0)
##precompute the cross-covariance for the beta-fields
sigma.B.C <- makeSigmaB(cov.pars.beta$pars,
dist = crossDist(loc.unobs.beta, loc.obs.beta),
type = object$cov.beta$covf,
nugget = cov.pars.beta$nugget,
ind2.to.1=Ind.2.1, sparse=TRUE)
##The right most part is F'*iSigma.tilde*(C-mu)
out$beta$EX <- out$beta$mu + matrix(sigma.B.C %*% (t(Fobs) %*% iSoo.C),
ncol=dim(out$beta$mu)[2])
##and variances
if( out$opts$pred.var ){
Sby.iSb.Sou <- as.matrix( i.sigma.B %*% t(sigma.B.C) )
Sby.iSb.Sou <- solveTriBlock(R.i.sigma.B.Y, Sby.iSb.Sou, transpose=TRUE)
Sby.iSb.Sou <- solveTriBlock(R.i.sigma.B.Y, Sby.iSb.Sou, transpose=FALSE)
##contribution from REML estimate
if( out$opts$type=="r" ){
var.beta.REML <- (cbind(rep(0,dimensions$L), bdiag(STdata$LUR.all)) -
sigma.B.C %*% (t(Fobs)%*%iSoo.Xtilde))
var.beta.REML <- as.matrix(var.beta.REML)
if( out$opts$beta.covar ){
##full matrix
V.REML <- (var.beta.REML %*% i.XSX) %*% t(var.beta.REML)
}else{
##only diagonal elements
V.REML <- rowSums( (var.beta.REML %*% i.XSX) * var.beta.REML )
}
}else{
V.REML <- 0
}##if( out$opts$type=="r" ){...}else{...}
if( out$opts$beta.covar ){
##sigma.B for unobserved locations
sigma.B.uu <- makeSigmaB(cov.pars.beta$pars,
dist = crossDist(loc.unobs.beta),
type = object$cov.beta$covf,
nugget = cov.pars.beta$nugget)
##compute variance matrix
tmp <- sigma.B.uu - sigma.B.C %*% tF.iS.F %*% Sby.iSb.Sou + V.REML
out$beta$VX.full <- list()
for(i in 1:dim(out$beta$EX)[2]){
Ind <- (1:dim(out$beta$EX)[1]) + (i-1)*dim(out$beta$EX)[1]
out$beta$VX.full[[i]] <- as.matrix(tmp[Ind, Ind, drop=FALSE])
rownames(out$beta$VX.full[[i]]) <- rownames(out$beta$EX)
colnames(out$beta$VX.full[[i]]) <- rownames(out$beta$EX)
}
names(out$beta$VX.full) <- colnames(out$beta$EX)
tmp <- diag(tmp)
}else{
##compute only diagonal elements of the STATIONARY sigma.beta matrix
sigma.B.uu <- makeSigmaB(cov.pars.beta$pars, dist = matrix(0,1,1),
type = object$cov.beta$covf,
nugget = cov.pars.beta$nugget)
##expand to full matrix
sigma.B.uu <- matrix(diag(sigma.B.uu), ncol=dim(sigma.B.uu)[1],
nrow=dim(loc.unobs.beta)[1], byrow=TRUE)
##and compute the relevant covariance
tmp <- (c(sigma.B.uu) - rowSums(sigma.B.C * t(tF.iS.F %*% Sby.iSb.Sou)) +
V.REML)
}##if( out$opts$beta.covar ){...}else{...}
out$beta$VX <- matrix(tmp, ncol=dim(out$beta$EX)[2])
dimnames(out$beta$VX) <- dimnames(out$beta$EX)
##remove variables not needed
rm(tmp, Sby.iSb.Sou, sigma.B.uu, V.REML)
}## if( out$opts$pred.var )
##remove variables not needed
rm(i.sigma.B, tF.iS.F)
######################################
### MATRICES HOLDING RETURN VALUES ###
##mean value of the field
out$EX.mu <- as.matrix( Xtilde.unobs %*% gamma.alpha )
##compute predicitons based on the beta fields
out$EX.mu.beta <- as.matrix( Funobs %*% c(out$beta$EX) )
##Add the spatio-temporal covariates, if any
if( dimensions$L!=0 ){
out$EX.mu.beta <- out$EX.mu.beta + (ST.unobs %*% out$pars$gamma.E)
}
if( !out$opts$only.obs ){
##reshape out$EX.mu and out$EX.mu.beta to match T-by-N
dim(out$EX.mu.beta) <- dim(out$EX.mu) <- c(T.unobs, N.unobs)
##and add names
colnames(out$EX.mu) <- STdata$locations$ID
rownames(out$EX.mu) <- as.character(STdata$trend$date)
dimnames(out$EX.mu.beta) <- dimnames(out$EX.mu)
}else{
dimnames(out$EX.mu.beta) <- dimnames(out$EX.mu) <- NULL
}
##lets save the contribution from the mean.
out$EX <- out$EX.mu
##additional fields for log-transformation
if( !is.null(out$opts$transform) ){
##temporary variables for transformed data
EX.trans <- out$EX.mu
EX.trans.pred <- out$EX.mu
##reserve places in the out structure
out$log.EX <- out$EX.pred <- NA
}
##Create matrices containing pointwise variances
##needed if a) variances requested OR b) log-transform
if( out$opts$pred.var || !is.null(out$opts$transform) ){
out$VX <- matrix(NA, dim(out$EX)[1], dim(out$EX)[2])
out$VX.pred <- matrix(NA, dim(out$EX)[1], dim(out$EX)[2])
##add names
if( !out$opts$only.obs ){
dimnames(out$VX) <- dimnames(out$VX.pred) <- dimnames(out$EX)
}
##and a list of covariances.
if( out$opts$pred.covar ){
out$VX.full <- list()
}
##additional fields for log-transformation
##(only if variances requested AND log-transform)
if( out$opts$pred.var && !is.null(out$opts$transform) ){
out$MSPE <- out$MSPE.pred <- out$VX
}
}##if( out$opts$pred.var )
##list holding LTA elements
if( out$opts$LTA ){
out$LTA <- vector("list", length(out$opts$LTA.list))
names(out$LTA) <- names(out$opts$LTA.list)
}
################################
### COMPUTE PREDICITONS OF X ###
##cross distance for the nu coordinates
cross.D.nu <- crossDist(loc.unobs.nu, loc.obs.nu)
## sigma.B.C * F'
sigma.B.C.tF <- sigma.B.C %*% t(Fobs)
##now we need to split the conditional expectation into parts to
##reduce the memory footprint
if( out$opts$pred.covar || out$opts$LTA ){
Ind.list <- split(1:length(out$EX), idx.unobs)
}else{
Ind.list <- split(1:length(out$EX),
rep(1:ceiling(length(out$EX)/Nmax),
each=Nmax)[1:length(out$EX)])
}
##loop over the parts
for( i in 1:length(Ind.list) ){
##index of the points we want to get conditional expectations for
Ind <- Ind.list[[i]]
##compute number of unobserved per block for this subset
T1.Ind <- T1[Ind]
for(j in c(1:length(date.all))){
nt.unobs[j] <- sum( T1.Ind==date.all[j] )
}
##order T1.Ind for computation of cross.covariance
T1.order <- order(T1.Ind)
##create full cross-covariance matrix for all the observations
sigma.B.full.C <- as.matrix(Funobs[Ind,,drop=FALSE] %*% sigma.B.C.tF)
##parameter order is sill, nugget, range (always predict with nugget=0)
sigma.nu.C <- makeSigmaNu(cov.pars.nu$pars, dist = cross.D.nu,
type = object$cov.nu$covf, nugget = 0,
random.effect = cov.pars.nu$random.effect,
ind1 = (idx.unobs[Ind])[T1.order],
ind2 = object$obs$idx,
blocks1 = nt.unobs, blocks2 = nt.obs,
ind2.to.1=Ind.2.1)
sigma.nu.C[T1.order,] <- sigma.nu.C
sigma.nu.C <- sigma.nu.C + sigma.B.full.C
##calculate conditional expectation
out$EX[Ind] <- out$EX.mu[Ind] + sigma.nu.C %*% iSoo.C
if( out$opts$LTA ){
##which site are we at and does it have a matching LTA request?
##asked here to avoid unnescesary VX.full computations.
ID.tmp <- STdata$locations$ID[unique(idx.unobs[Ind])]
##sanity check
if( length(ID.tmp)!=1 ){
stop("Something wrong with LTA-computations.")
}
LTA.tmp <- out$opts$LTA.list[[ID.tmp]]
}else{
LTA.tmp <- NULL
}
##calculate pointwise variance
if( out$opts$pred.var || !is.null(out$opts$transform) ){
##pick out the observed locations in this iteration
I.loc <- sort(unique( idx.unobs[Ind] ))
I.loc2 <- (rep(I.loc,dimensions$m) +
rep((0:(dimensions$m-1)) * N.unobs, each=length(I.loc)))
##compute relevant part of the covariance matrix
unobs.D.beta <- crossDist(loc.unobs.beta[I.loc,,drop=FALSE])
sigma.B.uu <- makeSigmaB(cov.pars.beta$pars, dist = unobs.D.beta,
type = object$cov.beta$covf,
nugget = cov.pars.beta$nugget, sparse=TRUE)
##first the unobserved covariance matrix
V.uu <- (Funobs[Ind,I.loc2,drop=FALSE] %*%
Matrix(sigma.B.uu %*% t(Funobs[Ind,I.loc2,drop=FALSE])))
##distance matrices for the unobserved nu locations
unobs.D.nu <- crossDist(loc.unobs.nu[I.loc,,drop=FALSE])
V.uu <- V.uu + makeSigmaNu(cov.pars.nu$pars, dist = unobs.D.nu,
type = object$cov.nu$covf, nugget = 0,
random.effect = cov.pars.nu$random.effect,
ind1 = idx.unobs[Ind]-min(I.loc)+1,
blocks1 = nt.unobs)
##compute iSigma.nu*Sigma.ou
iS.Sou <- i.sigma.nu %*% t(sigma.nu.C)
##compute F'*iSigma.nu*Sigma.ou
tF.iS.Sou <- as.matrix(tF.iS %*% t(sigma.nu.C))
##compute chol(inv(sigma.B.Y))' * (F'*iSigma.nu*Sigma.ou)
Sby.tF.iS.Sou <- solveTriBlock(R.i.sigma.B.Y, tF.iS.Sou, transpose=TRUE)
##contribution from REML estimate
if( out$opts$type=="r" ){
tmp <- Xtilde.unobs[Ind,,drop=FALSE] - sigma.nu.C %*% iSoo.Xtilde
if( out$opts$pred.covar || !is.null(LTA.tmp) ){
##full matrix
V.REML <- (tmp %*% i.XSX) %*% t(tmp)
if( !is.null(out$opts$transform) ){
##also compute LAMBDA correction for transform
lambda.full <- (tmp %*% i.XSX) %*% t(Xtilde.unobs[Ind,,drop=FALSE])
lambda <- diag(lambda.full)
}
}else{
##only diagonal elements
V.REML <- rowSums( (tmp %*% i.XSX) * tmp)
if( !is.null(out$opts$transform) ){
##also compute LAMBDA correction for transform
lambda <- rowSums((tmp %*% i.XSX) * (Xtilde.unobs[Ind,,drop=FALSE]))
}
}
}else{
V.REML <- 0
lambda <- 0
}
##combine parts, either to covariance matrix or VX for each estimate.
if( out$opts$pred.covar || !is.null(LTA.tmp) ){
##full matrix
tmp <- - sigma.nu.C %*% iS.Sou + t(Sby.tF.iS.Sou) %*% Sby.tF.iS.Sou
V.cond <- V.cond.0 <- as.matrix(V.uu + tmp + V.REML)
##add nugget to the diagonal
diag(V.cond) <- diag(V.cond) + out$opts$nugget.unobs[idx.unobs[Ind]]
##extract diagonal elements and store in output structure
out$VX[Ind] <- diag(V.cond.0)
out$VX.pred[Ind] <- diag(V.cond)
if(out$opts$pred.covar){
out$VX.full[[i]] <- V.cond.0
}
}else{
##only diagonal elements
tmp <- (- colSums(t(sigma.nu.C) * iS.Sou) +
colSums(Sby.tF.iS.Sou * Sby.tF.iS.Sou))
out$VX[Ind] <- diag(V.uu) + tmp + V.REML
out$VX.pred[Ind] <- out$VX[Ind] + out$opts$nugget.unobs[idx.unobs[Ind]]
}
##ensure positive variances
out$VX[Ind] <- pmax(out$VX[Ind], 0)
out$VX.pred[Ind] <- pmax(out$VX.pred[Ind], 0)
}else{
V.cond <- V.cond.0 <- NULL
}##if( out$opts$pred.var || !is.null(out$opts$transform) ){...}else{...}
##compute log-transform, store in EX.trans, EX.trans.pred
if( !is.null(out$opts$transform) ){
##use -lambda or -2*lambda, lambda=0 if type!="r"
if( out$opts$transform=="unbiased" ){ c <- 1 }else{ c <- 2 }
##expectations
EX.trans[Ind] <- exp( out$EX[Ind] + out$VX[Ind]/2 - c*lambda )
EX.trans.pred[Ind] <- exp( out$EX[Ind] + out$VX.pred[Ind]/2 - c*lambda )
##compute MSPE
if( out$opts$pred.var ){
##first compute/extract the unbiased estimate for each location
if( out$opts$transform=="mspe" ){
EX.ub <- exp( out$EX[Ind] + out$VX[Ind]/2 - lambda )
EX.ub.pred <- exp( out$EX[Ind] + out$VX.pred[Ind]/2 - lambda )
}else{
EX.ub <- EX.trans[Ind]
EX.ub.pred <- EX.trans.pred[Ind]
}
##prediction variance for unobserved locations (add nugget)
V.uu.pred <- V.uu + diag(out$opts$nugget.unobs[idx.unobs[Ind]])
if( !is.null(LTA.tmp) ){
##LTA to be computed, we'll need the full matrix
##then find the adjustment factor for the second term, only relevant for
##"unbiased" and type="r" (set to 1 o.w.)
if(out$opts$transform=="unbiased" && out$opts$type=="r"){
Z <- exp(lambda.full)
Z <- Z+t(Z) - exp(lambda.full + t(lambda.full))
}else{
Z <- 1
}
MSPE.part <- exp(V.uu) * (1 - exp(-V.cond.0)*Z)
MSPE.part.pred <- exp(V.uu.pred) * (1 - exp(-V.cond)*Z)
##Compute MSPE for EX and EX.pred
out$MSPE[Ind] <- (EX.ub*EX.ub * diag(MSPE.part))
out$MSPE.pred[Ind] <- (EX.ub.pred*EX.ub.pred * diag(MSPE.part.pred))
}else{
##then find the adjustment factor for the second term, only relevant for
##"unbiased" and type="r" (set to 1 o.w.)
if(out$opts$transform=="unbiased" && out$opts$type=="r"){
Z <- 2*exp(lambda)-exp(2*lambda)
}else{
Z <- 1
}
##Compute MSPE for EX and EX.pred
out$MSPE[Ind] <- (EX.ub*EX.ub * exp( diag(V.uu) ) *
( 1 - exp(-out$VX[Ind])*Z ))
out$MSPE.pred[Ind] <- (EX.ub.pred*EX.ub.pred * exp( diag(V.uu.pred) ) *
( 1 - exp(-out$VX.pred[Ind])*Z ))
}##if( !is.null(LTA.tmp) ){...}else{...}
##ensure positive MSPE
out$MSPE[Ind] <- pmax(out$MSPE[Ind], 0)
out$MSPE.pred[Ind] <- pmax(out$MSPE.pred[Ind], 0)
}else{
MSPE.part <- MSPE.part.pred <- NULL
}##if( out$opts$pred.var ){...}else{...}
if( !is.null(LTA.tmp) ){
##observatiosn for this location
EX.tmp <- cbind(exp(out$EX.mu[Ind]), exp(out$EX.mu.beta[Ind]),
EX.trans[Ind], EX.trans.pred[Ind])
colnames(EX.tmp) <- c("EX.mu", "EX.mu.beta", "EX", "EX.pred")
##compute mean value
out$LTA[[ID.tmp]] <- internalComputeLTA(LTA.tmp, EX.tmp, T1[Ind],
V=MSPE.part,
V.pred=MSPE.part.pred,
E.vec=EX.ub,
E.vec.pred=EX.ub.pred)
}##if( !is.null(LTA.tmp) )
}else{
##if not transform, seperate code for LTA (differences are too big...)
if( !is.null(LTA.tmp) ){
##observatiosn for this location
EX.tmp <- cbind(out$EX.mu[Ind], out$EX.mu.beta[Ind], out$EX[Ind])
colnames(EX.tmp) <- c("EX.mu", "EX.mu.beta", "EX")
##compute mean value
out$LTA[[ID.tmp]] <- internalComputeLTA(LTA.tmp, EX.tmp, T1[Ind],
V=V.cond.0, V.pred=V.cond)
}##if( !is.null(LTA.tmp) )
}##if( !is.null(out$opts$transform) ){...} else {...}
}##for( i in 1:length(Ind.list) )
##combine LTA results
if(out$opts$LTA){
##bind results
out$LTA <- do.call(cbind, out$LTA)
##and add names
tmp <- sapply(out$opts$LTA.list, length)
if( all(tmp==1) ){
colnames(out$LTA) <- names(tmp)
}else{
colnames(out$LTA) <- paste(rep(names(tmp),times=tmp),
unlist(lapply(tmp,seq)), sep=".")
}
##convert to data.frame
out$LTA <- as.data.frame(t(out$LTA))
}##if(out$opts$LTA)
if( !is.null(out$opts$transform) ){
##move EX of log-field
out$log.EX <- out$EX
##store EX of transformed field
out$EX <- EX.trans
out$EX.pred <- EX.trans.pred
##transform partial fields
out$EX.mu <- exp(out$EX.mu)
out$EX.mu.beta <- exp(out$EX.mu.beta)
##rename VX/VX.pred to log-fields
I <- grep("^VX",names(out))
names(out)[I] <- paste("log.",names(out)[I],sep="")
}
##add names to full prediction variances
if( out$opts$pred.covar ){
names(out$VX.full) <- colnames(out$EX)
for(i in 1:length(out$VX.full))
colnames(out$VX.full[[i]]) <- rownames(out$VX.full[[i]]) <- rownames(out$EX)
}##if(out$opts$pred.covar)
##index for the unobserved values into the predicted values
if( length(STdata$obs$obs)!=0 ){
if( out$opts$only.obs ){
I <- 1:length(out$EX)
}else{
I <- (match(STdata$obs$ID,colnames(out$EX))-1)*dim(out$EX)[1] +
match(STdata$obs$date, STdata$trend$date)
}
out$I <- data.frame(I=I, date=STdata$obs$date, ID=STdata$obs$ID,
stringsAsFactors=FALSE)
}##if( length(STdata$obs$obs)!=0 )
return( out )
}##function predict.STmodel
###################################
## S3 methods for predictSTmodel ##
###################################
##' \code{\link[base:print]{print}} method for class \code{predictSTmodel}.
##'
##' @title Print details for \code{predictSTmodel} object
##' @param x \code{predictSTmodel} object to print information for.
##' @param ... Ignored additional arguments.
##' @return Nothing
##'
##' @author Johan Lindstrom
##'
##' @examples
##' ##load data
##' data(pred.mesa.model)
##' print(pred.mesa.model)
##'
##'
##' @family predictSTmodel methods
##' @importFrom utils str
##' @method print predictSTmodel
##' @export
print.predictSTmodel <- function(x, ...){
##check class belonging
stCheckClass(x, "predictSTmodel", name="x")
cat("Prediction for STmodel.\n\n")
if( x$opts$only.pars ){
cat("Only computed regression parameters:\n")
}else{
cat("Regression parameters:\n")
}
if( length(x$pars$gamma.E)==0 ){
cat("\t","No spatio-temporal covariate.\n")
}else{
cat("\t", length(x$pars$gamma), "Spatio-temporal covariate(s).\n")
}
cat("\t", length(x$pars$alpha.E),
"beta-fields regression parameters in x$pars.\n\n")
if( x$opts$only.pars ){
return(invisible())
}
if( x$opts$type=="r" ){
cat("Regression parameters are assumed to be unknown and\n")
cat("\tprediction variances include uncertainties\n")
cat("\tin regression parameters.\n\n")
}else{
cat("Regression parameters are assumed to be known and\n")
cat("\tprediction variances do NOT include\n")
cat("\tuncertainties in regression parameters.\n\n")
}
cat("Prediction of beta-fields, (x$beta):\n")
str(x$beta,1)
cat("\n")
if( !is.null(x$opts$transform) ){
cat("Predictions for log-Gaussian field of type:",
x$opts$transform, "\n\n")
}
if( x$opts$only.obs ){
cat("Predictions only for", length(x$EX), "observations.\n")
}else{
cat("Predictions for", dim(x$EX)[1], "times at",
dim(x$EX)[2], "locations.\n")
}
str(x[grep("EX",names(x))],1)
cat("\n")
if( length(grep("VX",names(x)))!=0 ){
cat("Variances have been computed.\n")
str(x[grep("VX",names(x))],1)
}else{
cat("Variances have NOT been computed.\n")
}
cat("\n")
if( !is.null(x$opts$transform) ){
if( length(grep("MSPE",names(x)))!=0 ){
cat("Mean squared prediciton errors have been computed.\n")
str(x[grep("MSPE",names(x))],1)
}else{
cat("Mean squared prediciton errors NOT been computed.\n")
}
cat("\n")
}
if( isTRUE(x$opts$LTA) ){
cat(dim(x$LTA)[1], "temporal averages have been compute.\n")
str(x["LTA"],1)
}
cat("\n")
return(invisible())
}##function print.predictSTmodel
##' \code{\link[graphics:plot]{plot}} method for classes \code{predictSTmodel}
##' and \code{predCVSTmodel}. Provides several different plots of the
##' data.
##'
##' @title Plots for \code{predictSTmodel} and \code{predCVSTmodel} Objects
##'
##' @param x \code{predictSTmodel} or \code{predCVSTmodel} object to plot.
##' @param y Plot predictions as a function of either \code{"time"} or
##' \code{"obs"}ervations.
##' @param STmodel \code{STdata}/\code{STmodel} object containing observations
##' with which to compare the predictions (not used for
##' \code{plot.predCVSTmodel}).
##' @param ID The location for which we want to plot predictions. A
##' string matching names in \code{colnames(x$EX)} (or \code{x$I$ID},
##' number(s) which are used as \code{ID = colnames(x$EX)[ID]}, or
##' \code{"all"} in which case all predictions are used.
##' If several locations are given (or \code{"all"}) then
##' \code{y} must be \code{"obs"}.
##' @param col A vector of three colours: The first is the colour of the
##' predictions, second for the observations and third for the polygon
##' illustrating the confidence bands. \cr For \code{y="obs"} the colours are
##' 1) colour of the points, 2) colour of the 1-1 line, and 3) colour of the
##' polygon. If \code{ID="all"}, picking \code{col[1]="ID"} will colour code
##' the observations-prediction points by site ID.
##' @param pch,cex,lty,lwd Vectors with two elements giving the point type,
##' size, line type and line width to use when plotting the predictions and
##' observations respectively. Setting a value to \code{NA} will give no
##' points/lines for the predictions/observations. \cr
##' When plotting predictions
##' as a function of observations \code{lty[2]} is used for the addition of
##' \code{\link[graphics:abline]{abline}(0,1, lty=lty[2], col=col[2],
##' lwd=lwd[2])}; \code{pch[2]} and \code{cex[2]} are ignored.
##' @param p Width of the plotted confidence bands (as coverage percentage,
##' used to find appropriate two-sided normal quantiles).
##' @param pred.type Which type of prediction to plot, one of
##' \code{"EX"}, \code{"EX.mu"}, \code{"EX.mu.beta"}, or \code{"EX.pred"};
##' see the output from \code{\link{predict.STmodel}}
##' @param pred.var Should we plot confidence bands based on prediction (TRUE)
##' or confidence intrevalls (FALSE), see \code{\link{predict.STmodel}}.
##' Only relevant if \code{pred.type="EX"} or \code{pred.type="EX.pred"}. \cr
##' \strong{NOTE:} \emph{The default differs for \code{plot.predictSTmodel}
##' and \code{plot.predCVSTmodel}!}
##' @param add Add to existing plot?
##' @param ... Additional parameters passed to
##' \code{\link[graphics:plot]{plot}}.
##'
##' @return Nothing
##'
##' @example Rd_examples/Ex_plot_predictSTmodel.R
##'
##' @author Johan Lindstrom
##'
##' @family predictSTmodel methods
##' @method plot predictSTmodel
##' @export
plot.predictSTmodel <- function(x, y="time", STmodel=NULL, ID=x$I$ID[1],
col=c("black","red","grey"), pch=c(NA,NA),
cex=c(1,1), lty=c(1,1), lwd=c(1,1), p=0.95,
pred.type="EX", pred.var=FALSE,
add=FALSE, ...){
##check class belonging
stCheckClass(x, "predictSTmodel", name="x")
if( !is.null(STmodel) ){
stCheckClass(STmodel, c("STmodel","STdata"), name="STmodel")
}
##check for observations
if( x$opts$only.pars ){
message("Prediction structure contains only parameters.")
##return
return(invisible())
}
##we have to use y, cast to resonable name
tmp <- internalPlotPredictChecks(y, pred.type, pred.var, x$opts$transform)
plot.type <- tmp$plot.type
pred.type <- tmp$pred.type
pred.var <- tmp$pred.var
##check ID
ID <- internalFindIDplot(ID, colnames(x$EX))
ID.all <- internalCheckIDplot(ID, y)
##pick out data for plotting
if( !ID.all ){
##pick out predictions, two cases
if( x$opts$only.obs ){
I1 <- x$I$ID==ID
I2 <- 1
date <- x$I$date[I1]
}else{
I1 <- 1:dim(x[[ pred.type[1] ]])[1]
I2 <- ID
date <- rownames(x[[ pred.type[1] ]])
}
sd <- x[[pred.var]][I1,I2]
if( is.null(sd) ){
sd <- NA
}else{
sd <- sqrt(sd)
}
pred <- data.frame(x=x[[ pred.type[1] ]][I1, I2], sd=sd, date=date,
x.log=NA, stringsAsFactors=FALSE)
if( length(pred.type)==2 ){
pred$x.log <- x[[ pred.type[2] ]][I1, I2]
}
##pick out observations, if any
obs <- STmodel$obs[STmodel$obs$ID==ID,,drop=FALSE]
if( !is.null(obs) ){
tmp <- obs
obs <- matrix(NA, dim(pred)[1], 1)
rownames(obs) <- as.character(pred$date)
obs[as.character(tmp$date),] <- tmp$obs
}
}else{
##all only relevant when we have plot by observations
##pick out predictions, two cases
if(length(ID)==1 && ID=="all"){
ID <- unique(x$I$ID)
}
I.ID <- x$I$ID %in% ID
I <- x$I$I[I.ID]
sd <- x[[pred.var]][I]
if( is.null(sd) ){
sd <- NA
}else{
sd <- sqrt(sd)
}
pred <- data.frame(x=x[[ pred.type[1] ]][I], sd=sd,
date=x$I$date[I.ID], ID=x$I$ID[I.ID],
x.log=NA, stringsAsFactors=FALSE)
if( length(pred.type)==2 ){
pred$x.log <- x[[ pred.type[2] ]][I]
}
##pick out observations, if any
obs <- createDataMatrix(STmodel)
I <- (match(x$I$ID[I.ID],colnames(obs))-1)*dim(obs)[1] +
match(as.character(x$I$date[I.ID]), rownames(obs))
obs <- obs[I]
if( dim(pred)[1]!=length(obs) ){
stop("Number of observed locations do not match no. predicted.")
}
}##if( !ID.all ){...}else{...}
if( !is.null(x$opts$transform) && !is.null(obs) ){
##log-field, first fix observations
obs <- exp(obs)
}
internalPlotPredictions(plot.type, ID, pred, obs, col, pch, cex, lty, lwd, p,
add, x$opts$transform, ...)
##return
return(invisible())
}##function plot.predictSTmodel
Try the SpatioTemporal package in your browser
Any scripts or data that you put into this service are public.
SpatioTemporal documentation built on May 2, 2019, 8:49 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
############################################
## S3-METHOD THAT PREDICTS FROM A STmodel ##
############################################
##Functions in this file:
## predict.STmodel EX:ok
## print.predictSTmode EX:ok
## plot.predictSTmodel EX:ok
##' Compute the conditional expectations (i.e. predictions) at the unobserved
##' space-time locations. Predictions are computed for the space-time locations in
##' \code{object} and/or \code{STdata}, conditional on the observations (and
##' temporal trends) in \code{object} and parameters given in \code{x}.
##'
##' 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 \code{E(X_i| Y_i) != Y_i}, with the
##' difference being due to smoothing over the nugget. Further two possible
##' variance can be computed (see below), \code{V(X_i|Y_i)} and
##' \code{V(X_i|Y_i)+nugget_i}. Here the nugget for unobserved locations needs
##' to be specified as an additional argument \code{nugget.nobs}. The two
##' variances correspond, losely, to confidence and prediction intervals.
##'
##' Variances are computed if \code{pred.var=TRUE} point-wise variances for the
##' predictions (and the latent beta-fields) are
##' computed. If instead \code{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 \emph{not}, due
##' to memory restrictions, any covariances between locations.
##' \code{beta.covar=TRUE} gives the full covariance matrices for the latent
##' beta-fields.
##'
##' If \code{transform!="none"} the field is assumed to be log-Gaussian and
##' expectations are transformed, and if \code{pred.var=TRUE} the mean squared
##' prediction errors are given.
##'
##' @title Computes Conditional Expectation at Unobserved Locations
##'
##' @param object \code{STmodel} object for which to compute predictions.
##' @param x Model parameters for which to compute the conditional
##' expectation. Either as a vector/matrix or an \code{estimateSTmodel} from
##' \code{\link{estimate.STmodel}}.
##' @param STdata \code{STdata}/\code{STmodel} object with locations/times for
##' which to predict. If not given predictions are computed for locations/times
##' in \code{object}
##' @param Nmax Limits the size of matrices constructed when computing
##' expectations. Use a smaller value if memory becomes a problem.
##' @param only.pars Compute only the regression parameters (using GLS) along
##' with the related variance.
##' @param nugget.unobs Value of nugget at unonserved locations, either a scalar
##' or a vector with one element per unobserved site. \strong{NOTE:} All sites in
##' \code{STdata} are considered unobserved!
##' @param only.obs Compute predictions at only locations specified by
##' observations in \code{STdata}. Used to limit computations when doing
##' cross-validation.
##' \code{only.obs=TRUE} \emph{implies} \code{pred.covar=FALSE} and
##' \code{combine.data=FALSE}.
##' Further \code{\link{createSTmodel}} will be called on any \code{STdata}
##' input, possibly \emph{reordering the observations.}
##' @param pred.var,pred.covar Compute point-wise prediction variances; or
##' compute covariance matrices for the predicted time series at each location.
##' \code{pred.covar=TRUE} \emph{implies} \code{pred.var=TRUE} and sets
##' \code{Nmax} equal to the number of timepoints.
##' @param beta.covar Compute the full covariance matrix for the latent
##' beta-fields, otherwise only the diagonal elements of V(beta|obs) are
##' computed.
##' @param combine.data Combine \code{object} and \code{STdata} and predict for
##' the joint set of points, see \code{\link{c.STmodel}}.
##' @param type A single character indicating the type of prediction to
##' compute. Valid options are "f", "p", and "r", for \emph{full},
##' \emph{profile} or \emph{restricted maximum likelihood} (REML). For profile
##' and full the predictions are computed assuming that \emph{both} covariance
##' parameters and regression parameters are known,
##' e.g. \code{E(X|Y,cov_par,reg_par)}; for REML predictions are compute
##' assuming \emph{only} covariance parameters known,
##' e.g. \code{E(X|Y,cov_par)}. The main difference is that REML will have
##' \emph{larger} variances due to the additional uncertainty in the
##' regression parameters.
##' @param transform Regard field as log-Gaussian and apply exponential
##' transformation to predictions. For the final expectations two options
##' exist, either a unbiased prediction or the (biased) mean-squared error
##' predictions.
##' @param LTA Compute long-term temporal averages. Either a logical value or a
##' list; if \code{TRUE} then averages at each location (and variances if
##' \code{pred.var=TRUE}) are computed; otherwise this should be a list with
##' elements named after locations and each element containing a vector (or
##' list of vectors) with dates over which to compute averages. If
##' \code{only.obs=TRUE} averages are computed over only the observations.
##' @param ... Ignored additional arguments.
##'
##' @return The function returns a list containing (objects not computed
##' will be missing):
##' \item{opts}{Copy of options used in the function call.}
##' \item{pars}{A list with regression parameters and related variances.
##' \code{pars} contain \code{gamma.E} and \code{alpha.E} with
##' regression coefficients for the spatio-temporal model and
##' land-use covaraiates; variances are found in \code{gamma.V}
##' and \code{alpha.V}; cross-covariance between gamma and alpha in
##' \code{gamma.alpha.C}.}
##' \item{beta}{A list with estimates of the beta-fields, including the
##' regression mean \code{mu}, conditional expectations \code{EX},
##' possibly variances \code{VX}, and the full covariance matrix
##' \code{VX.full}.}
##' \item{EX.mu}{predictions based on the regression parameters, geographic
##' covariates, and temporal trends. I.e. only the deterministic
##' part of the spatio-temporal model.}
##' \item{EX.mu.beta}{Predictions based on the latent-beta fields, but excluding
##' the residual nu field.}
##' \item{EX}{Full predictions at the space-time locations in
##' \code{object} and/or \code{STdata}.}
##' \item{EX.pred}{Only for \code{transform!="none"}, full predictions
##' including bias correction for prediction error.}
##' \item{VX,VX.pred}{Pointwise variances and prediction variances (i.e. incl.
##' contribution from \code{nugget.unobs}) for all locations in \code{EX}.}
##' \item{VX.full}{A 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.}
##' \item{MSPE,MSPE.pred}{Pointwise mean-square prediction errors for the
##' log-Gaussian fields.}
##' \item{log.EX,log.VX.pred,log.VX}{Pointwise predictions and variances for
##' the un-transformed fields when \code{transform!="none"}}
##' \item{LTA}{A data.frame with temporal averages for locations specified by
##' \code{LTA}. }
##' \item{I}{A vector with the locations of the observations in \code{object} or
##' \code{STdata}. To extract predictions at the observations locations use
##' \code{EX[I]}.}
##'
##' @example Rd_examples/Ex_predict_STmodel.R
##'
##' @author Johan Lindstrom
##'
##' @family STmodel methods
##' @family predictSTmodel methods
##' @importFrom stats predict
##' @method predict STmodel
##' @export
predict.STmodel <- function(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", LTA=FALSE,
transform=c("none","unbiased","mspe"), ...){
##################################
### INITIAL SETUP AND CHECKING ###
##check class belongings
stCheckClass(object, "STmodel", name="object")
if( !is.null(STdata) ){
stCheckClass(STdata, c("STdata","STmodel"), name="STdata")
}
##first ensure that type is lower case
type <- tolower(type)
##check if type is valid
stCheckType(type)
##args for transform
transform <- match.arg(transform)
if( transform=="none" ){ transform <- NULL }
if( inherits(x,"estimateSTmodel") ){
x <- coef(x,"all")$par
}##if( inherits(x,"estimateSTmodel") )
##figure out a bunch of dimensions
##dimensions$L!=0 is used to check for spatio-temporal covariates
dimensions <- loglikeSTdim(object)
##check length of input parameters
if( type=="f" && length(x)!=dimensions$nparam ){
stop( paste("type=f, requires", dimensions$nparam,
"parameters but length(x) =", length(x)) )
}
if( type!="f" && length(x)==dimensions$nparam ){
##drop the regression parameters
x <- x[(dimensions$nparam-dimensions$nparam.cov+1):dimensions$nparam]
}
if( type!="f" && length(x)!=dimensions$nparam.cov ){
stop( paste("type!=f, requires", dimensions$nparam.cov,
"parameters but length(x) =", length(x)) )
}
##only.pars is strange when type=="f"
if( only.pars && type=="f"){
warning("only.pars=TRUE and type=(f)ull only returns KNOWN parameters.",
immediate. = TRUE)
}
##only.pars, and we can ignore a bunch of things
if( only.pars ){
only.obs <- pred.covar <- beta.covar <- combine.data <- LTA <-FALSE
transform <- NULL
}else{
##check if only.obs is valid
if( only.obs && is.null(STdata) ){
stop("only.obs=TRUE requires STdata.")
}
##do we have an object to combine with
if( combine.data && is.null(STdata) ){
warning("No data to combine with; predicting for 'object'", immediate. = TRUE)
combine.data <- FALSE
}
##can't combine data if we're predicting at only obs.
if(only.obs && combine.data){
warning("only.obs=TRUE implies combine.data=FALSE.", immediate. = TRUE)
combine.data <- FALSE
}
##computing prediction covariates require !only.obs and pred.var
if(pred.covar && !pred.var){
warning("pred.covar=TRUE implies pred.var=TRUE.", immediate. = TRUE)
pred.var <- TRUE
}
if(pred.covar && only.obs){
warning("only.obs=TRUE implies pred.covar=FALSE.", immediate. = TRUE)
pred.covar <- FALSE
}
##computing beta covariances requires pred.var
if(beta.covar && !pred.var){
warning("beta.covar=TRUE implies pred.var=TRUE.", immediate. = TRUE)
pred.var <- TRUE
}
if( !is.null(transform) ){
if( type!="r" && transform!="unbiased" ){
warning("transform 'unbiased' and 'mspe' are equivalent when type!='r'",
immediate. = TRUE)
transform <- "unbiased"
}
if(pred.covar){
warning("transform implies pred.covar=FALSE.", immediate. = TRUE)
pred.covar <- FALSE
}
}##if( !is.null(transform) )
}##if( only.pars ){...}else{...}
##create the STdata used for predictions
if( is.null(STdata) ){
##pure copy, predict in the data set
STdata <- object
}else if(combine.data){
##combine the two datasets, and use for predictions (expands object with
##ONLY locations from STdata)
STdata <- c(object, STdata)
}else{
if( is.null(object$trend.fnc) && dim(object$trend)[2]!=1 ){
##trend.fnc missing and trend is more than a constant
##(backwards compability)
object$trend.fnc <- internalCreateTrendFnc(object$trend)
}
if( !inherits(STdata,"STmodel") ){
##predict only at STdata, and STdata not of class STmodel: need to cast.
##First drop ST covariate if no covariate in object
if(dimensions$L==0){
STdata$SpatioTemporal <- NULL
}
##and fix the trend (right no of trends, names and dates)
suppressMessages(STdata <- updateTrend(STdata, fnc=object$trend.fnc,
extra.dates=STdata$trend$date))
##since we're not using covariances for the prediction locations
##(specified seperately in nugget.unobs), we just pick a simple covariance
##structure, avoiding problems with missing covariates/levels.
cov.nu <- object$cov.nu
cov.nu$nugget <- TRUE
##Create an STmodel from STdata
STdata <- createSTmodel(STdata, LUR=object$LUR.list, ST=object$ST.list,
cov.beta=object$cov.beta, cov.nu=cov.nu,
locations=object$locations.list,
scale=!is.null(object$scale.covars),
scale.covars=object$scale.covars)
}else{
##STdata is an STmodel object ->
##test for consistent covariates and scaling (only equal scaling allowed).
areSTmodelsConsistent(object, STdata, "STdata")
}##if( !inherits(STdata,"STmodel") ){...}{else{...}
##STdata is now of type STmodel, lets fix the trend
suppressMessages(STdata <- updateTrend.STmodel(STdata, fnc=object$trend.fnc))
}##if( is.null(STdata) ){...}else if(...){...}else{...}
##Should we compute LTAs?
if( !is.list(LTA) && LTA ){
##yes, create list
if( only.obs ){
##averages over observtions only,
##i.e. split observation dates by their locations.
LTA.list <- split(STdata$obs$date, STdata$obs$ID)
}else{
##averages over all predictions
LTA.list <- lapply(STdata$locations$ID,
function(x){ return(STdata$trend$date) })
names(LTA.list) <- STdata$locations$ID
}
}else if( is.list(LTA) ){
##yes, list given
LTA.list <- LTA
LTA <- TRUE
}else{
##no, empty list and set LTA=FALSE
LTA.list <- NULL
LTA <- FALSE
}##if( !is.list(LTA) && LTA ){...}else if(...){...}else{...}
##check the validity of LTA.list
if( LTA ){
##check validity of sites
LTA.missing <- !(names(LTA.list) %in% STdata$locations$ID)
if( any(LTA.missing) ){
warning("Removed ", sum(LTA.missing),
" elements from LTA not found in STdata$locations$ID",
immediate. = TRUE)
LTA.list <- LTA.list[ !LTA.missing ]
}
if( length(LTA.list)!=0 ){
##convert components to lists
LTA.list <- lapply(LTA.list,
function(x){ switch(is.list(x)+1, list(x), x) })
##check validity of dates
LTA.tot <- 0
for(i in 1:length(LTA.list)){
for(j in length(LTA.list[[i]])){
LTA.missing <- !(LTA.list[[i]][[j]] %in% STdata$trend$date)
LTA.list[[i]][[j]] <- LTA.list[[i]][[j]][ !LTA.missing ]
LTA.tot <- LTA.tot + sum(LTA.missing)
}
##drop empty points
LTA.list[[i]] <- LTA.list[[i]][ sapply(LTA.list[[i]], length)!=0 ]
}##for(i in 1:length(LTA.list))
if( LTA.tot!=0 ){
warning("Removed ", LTA.tot,
" dates not found in STdata$trend$date from all LTA elements",
immediate. = TRUE)
}
##drop empty points
LTA.list <- LTA.list[ sapply(LTA.list, length)!=0 ]
}##if( length(LTA.list)!=0 )
##nothing left -> don't compute LTA's
if( length(LTA.list)==0 ){
LTA <- FALSE
LTA.list <- NULL
}
}##if( LTA )
#########################################
### COMPUTE COVARIANCE MATRICES FOR Y ###
##extract parameters from x
tmp <- loglikeSTgetPars(x, object)
if(type=="f"){
gamma <- tmp$gamma
alpha <- tmp$alpha
}
cov.pars.beta <- tmp$cov.beta
cov.pars.nu <- tmp$cov.nu
##the following only needed for type!="f", or predictions
if(type!="f" || !only.pars){
##nugget for unobserved sites
nugget.unobs <- internalFixNuggetUnobs(nugget.unobs, STdata,
cov.pars.nu$nugget)
##extract sparse matrices
Fobs <- expandF(object$F, object$obs$idx, n.loc=dimensions$n.obs)
##Create the Xtilde = [M FX] matrix
Xtilde <- as.matrix( Fobs %*% bdiag(object$LUR) )
##Add the spatio-temporal covariate (if it exists)
if( dimensions$L!=0 ){
Xtilde <- cbind(object$ST, Xtilde)
}
##create covariance matrices, beta-field
i.sigma.B <- makeSigmaB(cov.pars.beta$pars, dist = object$D.beta,
type = object$cov.beta$covf,
nugget = cov.pars.beta$nugget)
##and nu-field
i.sigma.nu <- makeSigmaNu(cov.pars.nu$pars, dist = object$D.nu,
type = object$cov.nu$covf,
nugget = cov.pars.nu$nugget,
random.effect = cov.pars.nu$random.effect,
blocks1 = object$nt, ind1 = object$obs$idx,
sparse=TRUE)
##calculate block cholesky factor of the matrices, in-place to conserve memory
i.sigma.B <- makeCholBlock(i.sigma.B, n.blocks=dimensions$m)
i.sigma.nu <- chol(i.sigma.nu)
##invert the matrices, in-place to conserve memory
i.sigma.B <- invCholBlock(i.sigma.B, n.blocks=dimensions$m)
i.sigma.nu <- chol2inv(i.sigma.nu)
##F'*inv(sigma.nu)
tF.iS <- t(Fobs) %*% i.sigma.nu
##F'*inv(sigma.nu)*F
tF.iS.F <- as.matrix(tF.iS %*% Fobs)
##inv(sigma.B|Y) = F'*inv(sigma.nu)*F + inv(sigma.B)
R.i.sigma.B.Y <- tF.iS.F + i.sigma.B
##calculate cholesky factor of inv(sigma.B|Y)
R.i.sigma.B.Y <- chol(R.i.sigma.B.Y)
##compute iSigma.nu*Xtilde and iSigma.nu*Xtilde
iS.X <- i.sigma.nu %*% Xtilde
iS.Y <- i.sigma.nu %*% object$obs$obs
##compute F'*iSigma.nu*Xtilde and F'*iSigma.nu*Y
tF.iS.X <- as.matrix(t(Fobs) %*% iS.X)
tF.iS.Y <- as.matrix(t(Fobs) %*% iS.Y)
##compute inv(sigma.B|Y)*F'*iSigma.nu*Xtilde
iSBY.tF.iS.X <- solveTriBlock(R.i.sigma.B.Y, tF.iS.X, transpose=TRUE)
iSBY.tF.iS.X <- solveTriBlock(R.i.sigma.B.Y, iSBY.tF.iS.X, transpose=FALSE)
}#if(type!="f" || !only.pars)
###############################################
### REGRESISON PARAMETERS (alpha and gamma) ###
if(type=="f"){
##regression parameters given
gamma.E <- c(tmp$gamma)
alpha.E <- unlist(tmp$alpha)
##known parameters, set variances to zero
gamma.V <- matrix(0,length(gamma.E),length(gamma.E))
alpha.V <- matrix(0,length(alpha.E),length(alpha.E))
gamma.alpha.C <- matrix(0,length(gamma.E),length(alpha.E))
##create a combined vector
gamma.alpha <- as.matrix(c(gamma.E,alpha.E))
}else{
##compute inv(Xtilde'*Sigma^-1*Xtilde) (and force symmetric)
i.XSX <- as.matrix( t(Xtilde) %*% iS.X )
i.XSX <- symmpart(i.XSX - t(tF.iS.X) %*% iSBY.tF.iS.X)
##compute t(Xtilde'*Sigma^-1*Y)
tmp2 <- as.matrix( object$obs$obs %*% iS.X )
tmp2 <- tmp2 - t(tF.iS.Y) %*% iSBY.tF.iS.X
##compute alpha and gamma
gamma.alpha <- solve(i.XSX, t(tmp2))
##extract computed parameters
if(dimensions$L!=0){
gamma.E <- c(gamma.alpha[1:dimensions$L])
alpha.E <- c(gamma.alpha[-(1:dimensions$L)])
}else{
gamma.E <- double(0)
alpha.E <- c(gamma.alpha)
}
##compute the variances
i.XSX <- solve(i.XSX)
if(dimensions$L!=0){
gamma.V <- i.XSX[1:dimensions$L,1:dimensions$L,drop=FALSE]
alpha.V <- i.XSX[-c(1:dimensions$L),-c(1:dimensions$L),drop=FALSE]
gamma.alpha.C <- i.XSX[1:dimensions$L,-c(1:dimensions$L),drop=FALSE]
}else{
gamma.V <- matrix(0,0,0)
alpha.V <- i.XSX
gamma.alpha.C <- matrix(0,length(gamma.E),length(alpha.E))
}
}##if(type=="f"){...}else{...}
##force to matrices
gamma.E <- as.matrix(gamma.E)
alpha.E <- as.matrix(alpha.E)
##add names
names.tmp <- loglikeSTnames(object, TRUE)
names.tmp <- names.tmp[1:(dimensions$nparam-dimensions$nparam.cov)]
if(dimensions$L!=0){
##first for gamma
rownames(gamma.E) <- names.tmp[1:dimensions$L]
colnames(gamma.V) <- rownames(gamma.V) <- rownames(gamma.E)
rownames(gamma.alpha.C) <- rownames(gamma.E)
rownames(alpha.E) <- names.tmp[-(1:dimensions$L)]
}else{
rownames(alpha.E) <- names.tmp
}
##and then for alpha
colnames(alpha.V) <- rownames(alpha.V) <- rownames(alpha.E)
colnames(gamma.alpha.C) <- rownames(alpha.E)
###############################
### CONSTRUCT RETURN OBJECT ###
out <- list()
##set class
class(out) <- "predictSTmodel"
##save options used for predict.
out$opts <- list(only.pars=only.pars, nugget.unobs=nugget.unobs,
only.obs=only.obs, pred.var=pred.var, pred.covar=pred.covar,
beta.covar=beta.covar, combine.data=combine.data, type=type,
transform=transform, LTA=LTA, LTA.list=LTA.list)
##collect the regression results
out$pars <- list(gamma.E=gamma.E, alpha.E=alpha.E,
gamma.V=gamma.V, alpha.V=alpha.V,
gamma.alpha.C=gamma.alpha.C)
if( out$opts$only.pars ){
##return the parameters
return( out )
}
##remove variables not needed, regression parameters...
rm(gamma.E, alpha.E, gamma.V, alpha.V, gamma.alpha.C, names.tmp)
##and options stored elsewhere
rm(only.pars, nugget.unobs, only.obs, pred.var, pred.covar,
beta.covar, combine.data, type, LTA, LTA.list, transform)
####################################
### COMPUTE inv(Sigma_oo)*(Y-mu) ###
##compute F'*iSigma.nu*(C-mu) and iSigma.nu*(C-mu)
tF.iS.C <- tF.iS.Y - tF.iS.X %*% gamma.alpha
iS.C <- iS.Y - iS.X %*% gamma.alpha
##compute inv(sigma.B|Y)*F'*iSigma.nu*(C-mu)
iSBY.tF.iS.C <- solveTriBlock(R.i.sigma.B.Y, tF.iS.C, transpose=TRUE)
iSBY.tF.iS.C <- solveTriBlock(R.i.sigma.B.Y, iSBY.tF.iS.C, transpose=FALSE)
##compute iSigma.tilde*(C-mu)
##this is the same for all the conditionals...
iSoo.C <- as.matrix( iS.C - t( tF.iS ) %*% iSBY.tF.iS.C )
##for REML we also need iSigma.tilde*Xtilde
if( out$opts$type=="r" ){
iSoo.Xtilde <- as.matrix( iS.X - t(tF.iS) %*% iSBY.tF.iS.X )
}
##remove variables not needed
rm(iS.X, iS.Y, tF.iS.X, tF.iS.Y, iSBY.tF.iS.X, tF.iS.C, iS.C, iSBY.tF.iS.C)
###################################
### DEFINE PREDICTION LOCAITONS ###
##create matrices with site index and temporal trends for all
##OBS date has to be INCREASING, i.e. same sorting as for STdata$obs...
if( out$opts$only.obs ){
##unobserved points
idx.unobs <- STdata$obs$idx
##time of observation for each site
T1 <- STdata$obs$date
##trend functions for all points matrix
Funobs <- STdata$F
if( dimensions$L!=0 ){
##extract the spatio-temporal trends
ST.unobs <- STdata$ST
}
##size of the unobserved locations
N.unobs <- length( unique(idx.unobs) )
T.unobs <- length( unique(T1) )
}else{
##size of the unobserved locations
N.unobs <- dim(STdata$locations)[1]
T.unobs <- dim(STdata$trend)[1]
##unobserved points
idx.unobs <- rep(1:N.unobs, each=T.unobs)
##time of observation for each site
T1 <- rep(STdata$trend$date, N.unobs)
##trend functions for all points matrix
Funobs <- matrix(1, (T.unobs*N.unobs), dimensions$m)
for(i in (1:dimensions$m)){
##if not constant, find relevant row in $trend
if(colnames(STdata$F)[i] != "const"){
Funobs[,i] <- rep(STdata$trend[,colnames(STdata$F)[i]], N.unobs)
}
}##for(i in (1:dimensions$m))
if( dimensions$L!=0 ){
##extract the spatio-temporal trends
ST.unobs <- matrix(STdata$ST.all, (T.unobs*N.unobs), dimensions$L)
}##if( dimensions$L!=0 )
}##if( out$opts$only.obs ) ... else ...
##compute number of observations for each time-point,
##taking into account that there might be missing dates.
date.all <- sort(unique( c(object$trend$date, STdata$trend$date) ))
nt.unobs <- nt.obs <- double(length(date.all))
for(i in c(1:length(date.all))){
nt.obs[i] <- sum( object$obs$date==date.all[i] )
}
##extract sparse matrices
Funobs <- expandF(Funobs, idx.unobs, n.loc=N.unobs)
##compute LUR times the temporal trends [M F*X] for unobserved
Xtilde.unobs <- as.matrix( Funobs %*% bdiag(STdata$LUR.all) )
if( dimensions$L!=0 ){
Xtilde.unobs <- cbind(ST.unobs, Xtilde.unobs)
}
###########################
### COMPUTE BETA FIELDS ###
##compute LUR*alpha (mean values for the beta fields)
out$beta <- list()
out$beta$mu <- matrix(bdiag(STdata$LUR.all) %*% out$pars$alpha.E,
ncol=length(STdata$LUR.all))
colnames(out$beta$mu) <- names(STdata$LUR.all)
rownames(out$beta$mu) <- rownames( STdata$LUR.all[[1]] )
##create distance matrix between unobserved and observed elements
##unobserved locations
loc.unobs.nu <- STdata$locations[, c("x.nu","y.nu"), drop=FALSE]
loc.unobs.beta <- STdata$locations[, c("x.beta","y.beta"), drop=FALSE]
##observed locations
I.obs <- match( colnames(object$D.nu), object$locations$ID)
loc.obs.nu <- object$locations[I.obs, c("x.nu","y.nu"), drop=FALSE]
loc.obs.beta <- object$locations[I.obs, c("x.beta","y.beta"), drop=FALSE]
##We also have that in the stripped data the indecies may not match so we need a
##second vector giving which idx in striped matches original idx.
Ind.2.1 <- match(object$locations$ID[I.obs], STdata$locations$ID, nomatch=0)
##precompute the cross-covariance for the beta-fields
sigma.B.C <- makeSigmaB(cov.pars.beta$pars,
dist = crossDist(loc.unobs.beta, loc.obs.beta),
type = object$cov.beta$covf,
nugget = cov.pars.beta$nugget,
ind2.to.1=Ind.2.1, sparse=TRUE)
##The right most part is F'*iSigma.tilde*(C-mu)
out$beta$EX <- out$beta$mu + matrix(sigma.B.C %*% (t(Fobs) %*% iSoo.C),
ncol=dim(out$beta$mu)[2])
##and variances
if( out$opts$pred.var ){
Sby.iSb.Sou <- as.matrix( i.sigma.B %*% t(sigma.B.C) )
Sby.iSb.Sou <- solveTriBlock(R.i.sigma.B.Y, Sby.iSb.Sou, transpose=TRUE)
Sby.iSb.Sou <- solveTriBlock(R.i.sigma.B.Y, Sby.iSb.Sou, transpose=FALSE)
##contribution from REML estimate
if( out$opts$type=="r" ){
var.beta.REML <- (cbind(rep(0,dimensions$L), bdiag(STdata$LUR.all)) -
sigma.B.C %*% (t(Fobs)%*%iSoo.Xtilde))
var.beta.REML <- as.matrix(var.beta.REML)
if( out$opts$beta.covar ){
##full matrix
V.REML <- (var.beta.REML %*% i.XSX) %*% t(var.beta.REML)
}else{
##only diagonal elements
V.REML <- rowSums( (var.beta.REML %*% i.XSX) * var.beta.REML )
}
}else{
V.REML <- 0
}##if( out$opts$type=="r" ){...}else{...}
if( out$opts$beta.covar ){
##sigma.B for unobserved locations
sigma.B.uu <- makeSigmaB(cov.pars.beta$pars,
dist = crossDist(loc.unobs.beta),
type = object$cov.beta$covf,
nugget = cov.pars.beta$nugget)
##compute variance matrix
tmp <- sigma.B.uu - sigma.B.C %*% tF.iS.F %*% Sby.iSb.Sou + V.REML
out$beta$VX.full <- list()
for(i in 1:dim(out$beta$EX)[2]){
Ind <- (1:dim(out$beta$EX)[1]) + (i-1)*dim(out$beta$EX)[1]
out$beta$VX.full[[i]] <- as.matrix(tmp[Ind, Ind, drop=FALSE])
rownames(out$beta$VX.full[[i]]) <- rownames(out$beta$EX)
colnames(out$beta$VX.full[[i]]) <- rownames(out$beta$EX)
}
names(out$beta$VX.full) <- colnames(out$beta$EX)
tmp <- diag(tmp)
}else{
##compute only diagonal elements of the STATIONARY sigma.beta matrix
sigma.B.uu <- makeSigmaB(cov.pars.beta$pars, dist = matrix(0,1,1),
type = object$cov.beta$covf,
nugget = cov.pars.beta$nugget)
##expand to full matrix
sigma.B.uu <- matrix(diag(sigma.B.uu), ncol=dim(sigma.B.uu)[1],
nrow=dim(loc.unobs.beta)[1], byrow=TRUE)
##and compute the relevant covariance
tmp <- (c(sigma.B.uu) - rowSums(sigma.B.C * t(tF.iS.F %*% Sby.iSb.Sou)) +
V.REML)
}##if( out$opts$beta.covar ){...}else{...}
out$beta$VX <- matrix(tmp, ncol=dim(out$beta$EX)[2])
dimnames(out$beta$VX) <- dimnames(out$beta$EX)
##remove variables not needed
rm(tmp, Sby.iSb.Sou, sigma.B.uu, V.REML)
}## if( out$opts$pred.var )
##remove variables not needed
rm(i.sigma.B, tF.iS.F)
######################################
### MATRICES HOLDING RETURN VALUES ###
##mean value of the field
out$EX.mu <- as.matrix( Xtilde.unobs %*% gamma.alpha )
##compute predicitons based on the beta fields
out$EX.mu.beta <- as.matrix( Funobs %*% c(out$beta$EX) )
##Add the spatio-temporal covariates, if any
if( dimensions$L!=0 ){
out$EX.mu.beta <- out$EX.mu.beta + (ST.unobs %*% out$pars$gamma.E)
}
if( !out$opts$only.obs ){
##reshape out$EX.mu and out$EX.mu.beta to match T-by-N
dim(out$EX.mu.beta) <- dim(out$EX.mu) <- c(T.unobs, N.unobs)
##and add names
colnames(out$EX.mu) <- STdata$locations$ID
rownames(out$EX.mu) <- as.character(STdata$trend$date)
dimnames(out$EX.mu.beta) <- dimnames(out$EX.mu)
}else{
dimnames(out$EX.mu.beta) <- dimnames(out$EX.mu) <- NULL
}
##lets save the contribution from the mean.
out$EX <- out$EX.mu
##additional fields for log-transformation
if( !is.null(out$opts$transform) ){
##temporary variables for transformed data
EX.trans <- out$EX.mu
EX.trans.pred <- out$EX.mu
##reserve places in the out structure
out$log.EX <- out$EX.pred <- NA
}
##Create matrices containing pointwise variances
##needed if a) variances requested OR b) log-transform
if( out$opts$pred.var || !is.null(out$opts$transform) ){
out$VX <- matrix(NA, dim(out$EX)[1], dim(out$EX)[2])
out$VX.pred <- matrix(NA, dim(out$EX)[1], dim(out$EX)[2])
##add names
if( !out$opts$only.obs ){
dimnames(out$VX) <- dimnames(out$VX.pred) <- dimnames(out$EX)
}
##and a list of covariances.
if( out$opts$pred.covar ){
out$VX.full <- list()
}
##additional fields for log-transformation
##(only if variances requested AND log-transform)
if( out$opts$pred.var && !is.null(out$opts$transform) ){
out$MSPE <- out$MSPE.pred <- out$VX
}
}##if( out$opts$pred.var )
##list holding LTA elements
if( out$opts$LTA ){
out$LTA <- vector("list", length(out$opts$LTA.list))
names(out$LTA) <- names(out$opts$LTA.list)
}
################################
### COMPUTE PREDICITONS OF X ###
##cross distance for the nu coordinates
cross.D.nu <- crossDist(loc.unobs.nu, loc.obs.nu)
## sigma.B.C * F'
sigma.B.C.tF <- sigma.B.C %*% t(Fobs)
##now we need to split the conditional expectation into parts to
##reduce the memory footprint
if( out$opts$pred.covar || out$opts$LTA ){
Ind.list <- split(1:length(out$EX), idx.unobs)
}else{
Ind.list <- split(1:length(out$EX),
rep(1:ceiling(length(out$EX)/Nmax),
each=Nmax)[1:length(out$EX)])
}
##loop over the parts
for( i in 1:length(Ind.list) ){
##index of the points we want to get conditional expectations for
Ind <- Ind.list[[i]]
##compute number of unobserved per block for this subset
T1.Ind <- T1[Ind]
for(j in c(1:length(date.all))){
nt.unobs[j] <- sum( T1.Ind==date.all[j] )
}
##order T1.Ind for computation of cross.covariance
T1.order <- order(T1.Ind)
##create full cross-covariance matrix for all the observations
sigma.B.full.C <- as.matrix(Funobs[Ind,,drop=FALSE] %*% sigma.B.C.tF)
##parameter order is sill, nugget, range (always predict with nugget=0)
sigma.nu.C <- makeSigmaNu(cov.pars.nu$pars, dist = cross.D.nu,
type = object$cov.nu$covf, nugget = 0,
random.effect = cov.pars.nu$random.effect,
ind1 = (idx.unobs[Ind])[T1.order],
ind2 = object$obs$idx,
blocks1 = nt.unobs, blocks2 = nt.obs,
ind2.to.1=Ind.2.1)
sigma.nu.C[T1.order,] <- sigma.nu.C
sigma.nu.C <- sigma.nu.C + sigma.B.full.C
##calculate conditional expectation
out$EX[Ind] <- out$EX.mu[Ind] + sigma.nu.C %*% iSoo.C
if( out$opts$LTA ){
##which site are we at and does it have a matching LTA request?
##asked here to avoid unnescesary VX.full computations.
ID.tmp <- STdata$locations$ID[unique(idx.unobs[Ind])]
##sanity check
if( length(ID.tmp)!=1 ){
stop("Something wrong with LTA-computations.")
}
LTA.tmp <- out$opts$LTA.list[[ID.tmp]]
}else{
LTA.tmp <- NULL
}
##calculate pointwise variance
if( out$opts$pred.var || !is.null(out$opts$transform) ){
##pick out the observed locations in this iteration
I.loc <- sort(unique( idx.unobs[Ind] ))
I.loc2 <- (rep(I.loc,dimensions$m) +
rep((0:(dimensions$m-1)) * N.unobs, each=length(I.loc)))
##compute relevant part of the covariance matrix
unobs.D.beta <- crossDist(loc.unobs.beta[I.loc,,drop=FALSE])
sigma.B.uu <- makeSigmaB(cov.pars.beta$pars, dist = unobs.D.beta,
type = object$cov.beta$covf,
nugget = cov.pars.beta$nugget, sparse=TRUE)
##first the unobserved covariance matrix
V.uu <- (Funobs[Ind,I.loc2,drop=FALSE] %*%
Matrix(sigma.B.uu %*% t(Funobs[Ind,I.loc2,drop=FALSE])))
##distance matrices for the unobserved nu locations
unobs.D.nu <- crossDist(loc.unobs.nu[I.loc,,drop=FALSE])
V.uu <- V.uu + makeSigmaNu(cov.pars.nu$pars, dist = unobs.D.nu,
type = object$cov.nu$covf, nugget = 0,
random.effect = cov.pars.nu$random.effect,
ind1 = idx.unobs[Ind]-min(I.loc)+1,
blocks1 = nt.unobs)
##compute iSigma.nu*Sigma.ou
iS.Sou <- i.sigma.nu %*% t(sigma.nu.C)
##compute F'*iSigma.nu*Sigma.ou
tF.iS.Sou <- as.matrix(tF.iS %*% t(sigma.nu.C))
##compute chol(inv(sigma.B.Y))' * (F'*iSigma.nu*Sigma.ou)
Sby.tF.iS.Sou <- solveTriBlock(R.i.sigma.B.Y, tF.iS.Sou, transpose=TRUE)
##contribution from REML estimate
if( out$opts$type=="r" ){
tmp <- Xtilde.unobs[Ind,,drop=FALSE] - sigma.nu.C %*% iSoo.Xtilde
if( out$opts$pred.covar || !is.null(LTA.tmp) ){
##full matrix
V.REML <- (tmp %*% i.XSX) %*% t(tmp)
if( !is.null(out$opts$transform) ){
##also compute LAMBDA correction for transform
lambda.full <- (tmp %*% i.XSX) %*% t(Xtilde.unobs[Ind,,drop=FALSE])
lambda <- diag(lambda.full)
}
}else{
##only diagonal elements
V.REML <- rowSums( (tmp %*% i.XSX) * tmp)
if( !is.null(out$opts$transform) ){
##also compute LAMBDA correction for transform
lambda <- rowSums((tmp %*% i.XSX) * (Xtilde.unobs[Ind,,drop=FALSE]))
}
}
}else{
V.REML <- 0
lambda <- 0
}
##combine parts, either to covariance matrix or VX for each estimate.
if( out$opts$pred.covar || !is.null(LTA.tmp) ){
##full matrix
tmp <- - sigma.nu.C %*% iS.Sou + t(Sby.tF.iS.Sou) %*% Sby.tF.iS.Sou
V.cond <- V.cond.0 <- as.matrix(V.uu + tmp + V.REML)
##add nugget to the diagonal
diag(V.cond) <- diag(V.cond) + out$opts$nugget.unobs[idx.unobs[Ind]]
##extract diagonal elements and store in output structure
out$VX[Ind] <- diag(V.cond.0)
out$VX.pred[Ind] <- diag(V.cond)
if(out$opts$pred.covar){
out$VX.full[[i]] <- V.cond.0
}
}else{
##only diagonal elements
tmp <- (- colSums(t(sigma.nu.C) * iS.Sou) +
colSums(Sby.tF.iS.Sou * Sby.tF.iS.Sou))
out$VX[Ind] <- diag(V.uu) + tmp + V.REML
out$VX.pred[Ind] <- out$VX[Ind] + out$opts$nugget.unobs[idx.unobs[Ind]]
}
##ensure positive variances
out$VX[Ind] <- pmax(out$VX[Ind], 0)
out$VX.pred[Ind] <- pmax(out$VX.pred[Ind], 0)
}else{
V.cond <- V.cond.0 <- NULL
}##if( out$opts$pred.var || !is.null(out$opts$transform) ){...}else{...}
##compute log-transform, store in EX.trans, EX.trans.pred
if( !is.null(out$opts$transform) ){
##use -lambda or -2*lambda, lambda=0 if type!="r"
if( out$opts$transform=="unbiased" ){ c <- 1 }else{ c <- 2 }
##expectations
EX.trans[Ind] <- exp( out$EX[Ind] + out$VX[Ind]/2 - c*lambda )
EX.trans.pred[Ind] <- exp( out$EX[Ind] + out$VX.pred[Ind]/2 - c*lambda )
##compute MSPE
if( out$opts$pred.var ){
##first compute/extract the unbiased estimate for each location
if( out$opts$transform=="mspe" ){
EX.ub <- exp( out$EX[Ind] + out$VX[Ind]/2 - lambda )
EX.ub.pred <- exp( out$EX[Ind] + out$VX.pred[Ind]/2 - lambda )
}else{
EX.ub <- EX.trans[Ind]
EX.ub.pred <- EX.trans.pred[Ind]
}
##prediction variance for unobserved locations (add nugget)
V.uu.pred <- V.uu + diag(out$opts$nugget.unobs[idx.unobs[Ind]])
if( !is.null(LTA.tmp) ){
##LTA to be computed, we'll need the full matrix
##then find the adjustment factor for the second term, only relevant for
##"unbiased" and type="r" (set to 1 o.w.)
if(out$opts$transform=="unbiased" && out$opts$type=="r"){
Z <- exp(lambda.full)
Z <- Z+t(Z) - exp(lambda.full + t(lambda.full))
}else{
Z <- 1
}
MSPE.part <- exp(V.uu) * (1 - exp(-V.cond.0)*Z)
MSPE.part.pred <- exp(V.uu.pred) * (1 - exp(-V.cond)*Z)
##Compute MSPE for EX and EX.pred
out$MSPE[Ind] <- (EX.ub*EX.ub * diag(MSPE.part))
out$MSPE.pred[Ind] <- (EX.ub.pred*EX.ub.pred * diag(MSPE.part.pred))
}else{
##then find the adjustment factor for the second term, only relevant for
##"unbiased" and type="r" (set to 1 o.w.)
if(out$opts$transform=="unbiased" && out$opts$type=="r"){
Z <- 2*exp(lambda)-exp(2*lambda)
}else{
Z <- 1
}
##Compute MSPE for EX and EX.pred
out$MSPE[Ind] <- (EX.ub*EX.ub * exp( diag(V.uu) ) *
( 1 - exp(-out$VX[Ind])*Z ))
out$MSPE.pred[Ind] <- (EX.ub.pred*EX.ub.pred * exp( diag(V.uu.pred) ) *
( 1 - exp(-out$VX.pred[Ind])*Z ))
}##if( !is.null(LTA.tmp) ){...}else{...}
##ensure positive MSPE
out$MSPE[Ind] <- pmax(out$MSPE[Ind], 0)
out$MSPE.pred[Ind] <- pmax(out$MSPE.pred[Ind], 0)
}else{
MSPE.part <- MSPE.part.pred <- NULL
}##if( out$opts$pred.var ){...}else{...}
if( !is.null(LTA.tmp) ){
##observatiosn for this location
EX.tmp <- cbind(exp(out$EX.mu[Ind]), exp(out$EX.mu.beta[Ind]),
EX.trans[Ind], EX.trans.pred[Ind])
colnames(EX.tmp) <- c("EX.mu", "EX.mu.beta", "EX", "EX.pred")
##compute mean value
out$LTA[[ID.tmp]] <- internalComputeLTA(LTA.tmp, EX.tmp, T1[Ind],
V=MSPE.part,
V.pred=MSPE.part.pred,
E.vec=EX.ub,
E.vec.pred=EX.ub.pred)
}##if( !is.null(LTA.tmp) )
}else{
##if not transform, seperate code for LTA (differences are too big...)
if( !is.null(LTA.tmp) ){
##observatiosn for this location
EX.tmp <- cbind(out$EX.mu[Ind], out$EX.mu.beta[Ind], out$EX[Ind])
colnames(EX.tmp) <- c("EX.mu", "EX.mu.beta", "EX")
##compute mean value
out$LTA[[ID.tmp]] <- internalComputeLTA(LTA.tmp, EX.tmp, T1[Ind],
V=V.cond.0, V.pred=V.cond)
}##if( !is.null(LTA.tmp) )
}##if( !is.null(out$opts$transform) ){...} else {...}
}##for( i in 1:length(Ind.list) )
##combine LTA results
if(out$opts$LTA){
##bind results
out$LTA <- do.call(cbind, out$LTA)
##and add names
tmp <- sapply(out$opts$LTA.list, length)
if( all(tmp==1) ){
colnames(out$LTA) <- names(tmp)
}else{
colnames(out$LTA) <- paste(rep(names(tmp),times=tmp),
unlist(lapply(tmp,seq)), sep=".")
}
##convert to data.frame
out$LTA <- as.data.frame(t(out$LTA))
}##if(out$opts$LTA)
if( !is.null(out$opts$transform) ){
##move EX of log-field
out$log.EX <- out$EX
##store EX of transformed field
out$EX <- EX.trans
out$EX.pred <- EX.trans.pred
##transform partial fields
out$EX.mu <- exp(out$EX.mu)
out$EX.mu.beta <- exp(out$EX.mu.beta)
##rename VX/VX.pred to log-fields
I <- grep("^VX",names(out))
names(out)[I] <- paste("log.",names(out)[I],sep="")
}
##add names to full prediction variances
if( out$opts$pred.covar ){
names(out$VX.full) <- colnames(out$EX)
for(i in 1:length(out$VX.full))
colnames(out$VX.full[[i]]) <- rownames(out$VX.full[[i]]) <- rownames(out$EX)
}##if(out$opts$pred.covar)
##index for the unobserved values into the predicted values
if( length(STdata$obs$obs)!=0 ){
if( out$opts$only.obs ){
I <- 1:length(out$EX)
}else{
I <- (match(STdata$obs$ID,colnames(out$EX))-1)*dim(out$EX)[1] +
match(STdata$obs$date, STdata$trend$date)
}
out$I <- data.frame(I=I, date=STdata$obs$date, ID=STdata$obs$ID,
stringsAsFactors=FALSE)
}##if( length(STdata$obs$obs)!=0 )
return( out )
}##function predict.STmodel
###################################
## S3 methods for predictSTmodel ##
###################################
##' \code{\link[base:print]{print}} method for class \code{predictSTmodel}.
##'
##' @title Print details for \code{predictSTmodel} object
##' @param x \code{predictSTmodel} object to print information for.
##' @param ... Ignored additional arguments.
##' @return Nothing
##'
##' @author Johan Lindstrom
##'
##' @examples
##' ##load data
##' data(pred.mesa.model)
##' print(pred.mesa.model)
##'
##'
##' @family predictSTmodel methods
##' @importFrom utils str
##' @method print predictSTmodel
##' @export
print.predictSTmodel <- function(x, ...){
##check class belonging
stCheckClass(x, "predictSTmodel", name="x")
cat("Prediction for STmodel.\n\n")
if( x$opts$only.pars ){
cat("Only computed regression parameters:\n")
}else{
cat("Regression parameters:\n")
}
if( length(x$pars$gamma.E)==0 ){
cat("\t","No spatio-temporal covariate.\n")
}else{
cat("\t", length(x$pars$gamma), "Spatio-temporal covariate(s).\n")
}
cat("\t", length(x$pars$alpha.E),
"beta-fields regression parameters in x$pars.\n\n")
if( x$opts$only.pars ){
return(invisible())
}
if( x$opts$type=="r" ){
cat("Regression parameters are assumed to be unknown and\n")
cat("\tprediction variances include uncertainties\n")
cat("\tin regression parameters.\n\n")
}else{
cat("Regression parameters are assumed to be known and\n")
cat("\tprediction variances do NOT include\n")
cat("\tuncertainties in regression parameters.\n\n")
}
cat("Prediction of beta-fields, (x$beta):\n")
str(x$beta,1)
cat("\n")
if( !is.null(x$opts$transform) ){
cat("Predictions for log-Gaussian field of type:",
x$opts$transform, "\n\n")
}
if( x$opts$only.obs ){
cat("Predictions only for", length(x$EX), "observations.\n")
}else{
cat("Predictions for", dim(x$EX)[1], "times at",
dim(x$EX)[2], "locations.\n")
}
str(x[grep("EX",names(x))],1)
cat("\n")
if( length(grep("VX",names(x)))!=0 ){
cat("Variances have been computed.\n")
str(x[grep("VX",names(x))],1)
}else{
cat("Variances have NOT been computed.\n")
}
cat("\n")
if( !is.null(x$opts$transform) ){
if( length(grep("MSPE",names(x)))!=0 ){
cat("Mean squared prediciton errors have been computed.\n")
str(x[grep("MSPE",names(x))],1)
}else{
cat("Mean squared prediciton errors NOT been computed.\n")
}
cat("\n")
}
if( isTRUE(x$opts$LTA) ){
cat(dim(x$LTA)[1], "temporal averages have been compute.\n")
str(x["LTA"],1)
}
cat("\n")
return(invisible())
}##function print.predictSTmodel
##' \code{\link[graphics:plot]{plot}} method for classes \code{predictSTmodel}
##' and \code{predCVSTmodel}. Provides several different plots of the
##' data.
##'
##' @title Plots for \code{predictSTmodel} and \code{predCVSTmodel} Objects
##'
##' @param x \code{predictSTmodel} or \code{predCVSTmodel} object to plot.
##' @param y Plot predictions as a function of either \code{"time"} or
##' \code{"obs"}ervations.
##' @param STmodel \code{STdata}/\code{STmodel} object containing observations
##' with which to compare the predictions (not used for
##' \code{plot.predCVSTmodel}).
##' @param ID The location for which we want to plot predictions. A
##' string matching names in \code{colnames(x$EX)} (or \code{x$I$ID},
##' number(s) which are used as \code{ID = colnames(x$EX)[ID]}, or
##' \code{"all"} in which case all predictions are used.
##' If several locations are given (or \code{"all"}) then
##' \code{y} must be \code{"obs"}.
##' @param col A vector of three colours: The first is the colour of the
##' predictions, second for the observations and third for the polygon
##' illustrating the confidence bands. \cr For \code{y="obs"} the colours are
##' 1) colour of the points, 2) colour of the 1-1 line, and 3) colour of the
##' polygon. If \code{ID="all"}, picking \code{col[1]="ID"} will colour code
##' the observations-prediction points by site ID.
##' @param pch,cex,lty,lwd Vectors with two elements giving the point type,
##' size, line type and line width to use when plotting the predictions and
##' observations respectively. Setting a value to \code{NA} will give no
##' points/lines for the predictions/observations. \cr
##' When plotting predictions
##' as a function of observations \code{lty[2]} is used for the addition of
##' \code{\link[graphics:abline]{abline}(0,1, lty=lty[2], col=col[2],
##' lwd=lwd[2])}; \code{pch[2]} and \code{cex[2]} are ignored.
##' @param p Width of the plotted confidence bands (as coverage percentage,
##' used to find appropriate two-sided normal quantiles).
##' @param pred.type Which type of prediction to plot, one of
##' \code{"EX"}, \code{"EX.mu"}, \code{"EX.mu.beta"}, or \code{"EX.pred"};
##' see the output from \code{\link{predict.STmodel}}
##' @param pred.var Should we plot confidence bands based on prediction (TRUE)
##' or confidence intrevalls (FALSE), see \code{\link{predict.STmodel}}.
##' Only relevant if \code{pred.type="EX"} or \code{pred.type="EX.pred"}. \cr
##' \strong{NOTE:} \emph{The default differs for \code{plot.predictSTmodel}
##' and \code{plot.predCVSTmodel}!}
##' @param add Add to existing plot?
##' @param ... Additional parameters passed to
##' \code{\link[graphics:plot]{plot}}.
##'
##' @return Nothing
##'
##' @example Rd_examples/Ex_plot_predictSTmodel.R
##'
##' @author Johan Lindstrom
##'
##' @family predictSTmodel methods
##' @method plot predictSTmodel
##' @export
plot.predictSTmodel <- function(x, y="time", STmodel=NULL, ID=x$I$ID[1],
col=c("black","red","grey"), pch=c(NA,NA),
cex=c(1,1), lty=c(1,1), lwd=c(1,1), p=0.95,
pred.type="EX", pred.var=FALSE,
add=FALSE, ...){
##check class belonging
stCheckClass(x, "predictSTmodel", name="x")
if( !is.null(STmodel) ){
stCheckClass(STmodel, c("STmodel","STdata"), name="STmodel")
}
##check for observations
if( x$opts$only.pars ){
message("Prediction structure contains only parameters.")
##return
return(invisible())
}
##we have to use y, cast to resonable name
tmp <- internalPlotPredictChecks(y, pred.type, pred.var, x$opts$transform)
plot.type <- tmp$plot.type
pred.type <- tmp$pred.type
pred.var <- tmp$pred.var
##check ID
ID <- internalFindIDplot(ID, colnames(x$EX))
ID.all <- internalCheckIDplot(ID, y)
##pick out data for plotting
if( !ID.all ){
##pick out predictions, two cases
if( x$opts$only.obs ){
I1 <- x$I$ID==ID
I2 <- 1
date <- x$I$date[I1]
}else{
I1 <- 1:dim(x[[ pred.type[1] ]])[1]
I2 <- ID
date <- rownames(x[[ pred.type[1] ]])
}
sd <- x[[pred.var]][I1,I2]
if( is.null(sd) ){
sd <- NA
}else{
sd <- sqrt(sd)
}
pred <- data.frame(x=x[[ pred.type[1] ]][I1, I2], sd=sd, date=date,
x.log=NA, stringsAsFactors=FALSE)
if( length(pred.type)==2 ){
pred$x.log <- x[[ pred.type[2] ]][I1, I2]
}
##pick out observations, if any
obs <- STmodel$obs[STmodel$obs$ID==ID,,drop=FALSE]
if( !is.null(obs) ){
tmp <- obs
obs <- matrix(NA, dim(pred)[1], 1)
rownames(obs) <- as.character(pred$date)
obs[as.character(tmp$date),] <- tmp$obs
}
}else{
##all only relevant when we have plot by observations
##pick out predictions, two cases
if(length(ID)==1 && ID=="all"){
ID <- unique(x$I$ID)
}
I.ID <- x$I$ID %in% ID
I <- x$I$I[I.ID]
sd <- x[[pred.var]][I]
if( is.null(sd) ){
sd <- NA
}else{
sd <- sqrt(sd)
}
pred <- data.frame(x=x[[ pred.type[1] ]][I], sd=sd,
date=x$I$date[I.ID], ID=x$I$ID[I.ID],
x.log=NA, stringsAsFactors=FALSE)
if( length(pred.type)==2 ){
pred$x.log <- x[[ pred.type[2] ]][I]
}
##pick out observations, if any
obs <- createDataMatrix(STmodel)
I <- (match(x$I$ID[I.ID],colnames(obs))-1)*dim(obs)[1] +
match(as.character(x$I$date[I.ID]), rownames(obs))
obs <- obs[I]
if( dim(pred)[1]!=length(obs) ){
stop("Number of observed locations do not match no. predicted.")
}
}##if( !ID.all ){...}else{...}
if( !is.null(x$opts$transform) && !is.null(obs) ){
##log-field, first fix observations
obs <- exp(obs)
}
internalPlotPredictions(plot.type, ID, pred, obs, col, pch, cex, lty, lwd, p,
add, x$opts$transform, ...)
##return
return(invisible())
}##function plot.predictSTmodel
Try the SpatioTemporal package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.