Nothing
###############################################
## FUNCTIONS THAT MANIPULATE STmodel OBJECTS ##
###############################################
##Functions in this file:
## updateCovf EX:ok
## createLUR EX:INTERNAL
## createST EX:INTERNAL
#####################################################################
## Helper function that updates covariance function specifications ##
#####################################################################
##' Updates/sets the covariance functions for \code{STmodel} objects. Used by
##' \code{\link{createSTmodel}}.
##'
##' The covariance function is specified using lists for \code{cov.beta} and
##' \code{cov.nu}. The lists should contain the following elements:
##' \describe{
##' \item{covf}{The type of covariance function(s), see
##' \code{\link{namesCovFuns}}.}
##' \item{nugget}{For the beta-fields: a vector of \code{TRUE}/\code{FALSE}
##' indicating if each beta-field should contain a nugget. \cr
##' For the nu-field: Either \code{TRUE}/\code{FALSE} for constant nugget/no
##' nugget; a formula; or length=1 character vector. For the latter two the
##' nugget is allowed to vary as \code{exp(B*theta)} where: \cr
##' \code{nugget = as.formula(paste("~", paste(cov.nu$nugget, collapse="+")))} \cr
##' \code{ covars = model.frame(nugget, covars, drop.unused.levels=TRUE)} \cr
##' \code{ B=model.matrix(nugget, covars)} \cr
##' \code{ B=as.matrix(B) } \cr
##' The resulting regression matrix is stored as \code{STmodel$cov.nu$nugget.matrix}
##' giving nugget for the observed locations. Unobserved locations are assumed
##' to have a zero nugget.}
##' \item{random.effect}{Only used for \code{cov.nu}, \code{TRUE}/\code{FALSE}
##' indicating if a random.effect for the mean value should be included, see
##' \code{\link{makeSigmaNu}}.}
##' }
##'
##' @title Update Covariance Functions in \code{STmodel} Objects
##' @param STmodel \code{STmodel} object with observations, covariates, trends, etc;
##' see \code{\link{mesa.model}}.
##' @param cov.beta,cov.nu Covariance specification for the beta- and nu-fields
##' should contain fields \code{covf}, \code{nugget}, and \code{random.effect}
##' (for the nu field); see details for description of these fields.
##' For \code{cov.beta} the fields should contain one element for each
##' smooth-temporal trend / beta-field if the fields have only one element,
##' these elements are repeated implying the same covariance for all beta-fields.
##' @return updated version of \code{STmodel} with new covariance specifications.
##'
##' @example Rd_examples/Ex_updateCovf.R
##'
##' @author Johan Lindstrom
##' @family STmodel functions
##' @family covariance functions
##' @importFrom stats as.formula
##' @importFrom stats model.frame
##' @importFrom stats model.matrix
##' @export
updateCovf <- function(STmodel, cov.beta=STmodel$cov.beta,
cov.nu=STmodel$cov.nu){
##check class belonging
stCheckClass(STmodel, "STmodel", name="STmodel")
##check for fields
stCheckFields(cov.beta, c("covf","nugget"), name="cov.beta")
stCheckFields(cov.nu, c("covf","nugget","random.effect"), name="cov.nu")
##find the number of temporal trends (including intercept)
nt <- dim(STmodel$trend)[2]
##expand elements in cov.beta
for(i in 1:length(cov.beta) ){
if( length(cov.beta[[i]])==1 ){
cov.beta[[i]] <- rep(cov.beta[[i]], nt)
}else if( length(cov.beta[[i]])!=nt ){
stop( sprintf("Covariance specification, %s, for beta has %d!=%d elements.",
names(cov.beta)[i], length(cov.beta[[i]]), nt) )
}
}##for(i in 1:length(cov.beta) )
##check that covariance functions are valid
tmp <- unlist(lapply( parsCovFuns( c(cov.beta$covf, cov.nu$covf) ), is.null))
if( any(tmp) ){
stop( paste("Unknown covariance specification(s):",
paste(names(tmp)[tmp], collapse=", ")) )
}
##check random effect
if( !is.logical(cov.nu$random.effect) ){
stop("'random.effect' specification must be 'logical'.")
}
##check nugget - beta
if( !is.logical(cov.beta$nugget) ){
stop("'nugget' specification for beta-fields must be 'logical'.")
}
##check iid and nugget - beta
if( any(cov.beta$covf=="iid" & !cov.beta$nugget) ){
stop("beta-fields: Covariance model iid requires a nugget.")
}
##check nugget - nu
##... but first drop unobserved locations from the covars
covars <- STmodel$locations
Ind <- STmodel$locations$ID %in% unique(STmodel$obs$ID)
covars <- covars[Ind,,drop=FALSE]
if( is.logical(cov.nu$nugget) ){
##either constant or nothing.
if( cov.nu$nugget[1]==TRUE ){
cov.nu$nugget <- as.formula("~1", env=.GlobalEnv)
}else{
cov.nu$nugget <- FALSE
cov.nu$nugget.matrix <- matrix(NULL, dim(covars)[1], 0)
}
}
if( class(cov.nu$nugget)=="formula" || is.character(cov.nu$nugget) ){
if( is.character(cov.nu$nugget) ){
cov.nu$nugget <- as.formula(paste("~", paste(cov.nu$nugget, collapse="+")),
env=.GlobalEnv)
}
covars.tmp <- model.frame(cov.nu$nugget, covars, drop.unused.levels=TRUE)
##A formula has been specified, use this
tmp <- tryCatch(model.matrix(cov.nu$nugget, covars.tmp),
error = function(e){
tmp <- model.matrix(cov.nu$nugget, covars);
##find constant columns
test <- apply(tmp, 2, function(x){length(unique(x))})
test <- test!=1 | names(test)=="(Intercept)"
warning(
paste("Some covariates are constant or factors with only one level",
"in the specification of the nugget. This will result in",
"unidentifiable parameters. To avoid this the following elemets",
"of the formula have been dropped:",
paste(names(test[!test]),collapse=", "), sep="\n") )
##drop constant columns
tmp <- tmp[,test, drop=FALSE];
return(tmp)})
##store nugget matrix.
cov.nu$nugget.matrix <- as.matrix( tmp )
}else if( !is.logical(cov.nu$nugget) ){
stop( paste("Unknown specification of nugget for nu-field:", cov.nu$nugget) )
}
##check iid and nugget - beta
if( dim(cov.nu$nugget.matrix)[2]==0 && cov.nu$covf=="iid" ){
stop("nu-field: Covariance model iid requires a nugget.")
}
##add rownames
rownames(cov.nu$nugget.matrix) <- covars$ID
##return
STmodel$cov.beta <- cov.beta
STmodel$cov.nu <- cov.nu
return( STmodel )
}##function updateCovf
#####################################################
## Helper functions that add LUR and ST covariates ##
#####################################################
##' Extracts the requested geographic and spatio-temporal covariates from
##' a \code{STmodel} object and formats them into suitable matrices.
##' For \strong{INTERNAL} use by \code{\link{createSTmodel}}
##'
##' @title Add Covariate Fields to \code{STdata} Object.
##' @param STdata \code{STdata} object with observations, covariates, trends,
##' etc; see \code{\link{mesa.data.raw}}.
##' @param LUR.list specification of covariates; e.g. output from
##' \code{\link{processLUR}}.
##'
##' @return \code{STdata} with added fields \code{LUR}, or \code{ST} and
##' \code{ST.all}.
##'
##' @author Johan Lindstrom
createLUR <- function(STdata, LUR.list){
##check class belonging
stCheckClass(STdata, "STdata", name="STdata")
##extract important fields
covars <- STdata$covars
scale.covars <- STdata$scale.covars
##should we scale the covariates
if( !is.null(scale.covars) ){
##check that we have all the covariates
stCheckFields(scale.covars$mean, names(covars), name="scale.covars$mean")
stCheckFields(scale.covars$sd, names(covars), name="scale.covars$sd")
##ensure right order
scale.covars$mean <- scale.covars$mean[ names(covars) ]
scale.covars$sd <- scale.covars$sd[ names(covars) ]
##scale all non NA:s
I <- !is.na(scale.covars$mean)
n1 <- dim(covars)[1]
n2 <- sum(I)
covars[,I] <- scale(covars[,I,drop=FALSE], center=scale.covars$mean[I],
scale=scale.covars$sd[I])
}
##create a list of LUR:s, one for each temporal trend
STdata$LUR.all <- vector("list", length(LUR.list))
names(STdata$LUR.all) <- names(LUR.list)
for( i in 1:length(LUR.list) ){
STdata$LUR.all[[i]] <- model.matrix(LUR.list[[i]], covars)
rownames(STdata$LUR.all[[i]]) <- covars$ID
}
##also create a reduced LUR for observations only
if( length(STdata$obs$obs)!=0 ){
STdata$LUR <- STdata$LUR.all
I.idx <- 1:max(match(STdata$obs$ID,STdata$covars$ID))
for( i in 1:length(STdata$LUR) ){
STdata$LUR[[i]] <- STdata$LUR[[i]][I.idx,,drop=FALSE]
}
}
return(STdata)
}##function createLUR
##' @param ST.list specification of spatio-temporal covariates; e.g. output from
##' \code{\link{processST}}.
##' @rdname createLUR
createST <- function(STdata, ST.list){
##check class belonging
stCheckClass(STdata, "STdata", name="STdata")
##extract important fields
SpatioTemporal <- STdata$SpatioTemporal
obs <- STdata$obs
if( length(ST.list)!=0 ){
##ST-covariate for all locations
STdata$ST.all <- SpatioTemporal[,,ST.list,drop=FALSE]
##and now for the observations
STdata$ST <- matrix(0, dim(obs)[1], length(ST.list))
##dates and ID:s for the spatio-temporal covaraites
ST.dates <- convertCharToDate( rownames(SpatioTemporal) )
ST.ID <- colnames(SpatioTemporal)
##correspinding index
I <- (match(obs$ID,ST.ID)-1)*length(ST.dates) + match(obs$date,ST.dates)
##extract these points from the Spatio temporal covariate
for(i in 1:length(ST.list)){
tmp <- SpatioTemporal[,,ST.list[i],drop=FALSE]
STdata$ST[,i] <- tmp[I]
}
colnames(STdata$ST) <- ST.list
STdata$ST <- data.matrix(STdata$ST)
}
return( STdata )
}##function createST
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.