Nothing
R/temporal_SVD_incomplete_data.R
In SpatioTemporal: Spatio-Temporal Model Estimation
#######################################################################
## FUNCTIONS THAT COMPUTE SVDs FOR TEMPORAL TRENDS WITH MISSING DATA ##
#######################################################################
##Functions in this file:
## SVDmiss - EX:ok
## SVDsmooth - EX:ok
## SVDsmoothCV - EX:with SVDsmooth
## calcSmoothTrends - EX:ok
## plot.SVDcv - EX:with SVDsmooth
## boxplot.SVDcv - EX:with SVDsmooth
## print.SVDcv - EX:with SVDsmooth
## updateTrend.STdata - EX:ok
## updateTrend.STmodel - EX:with updateTrend.STdata
## updateTrend - EX:with updateTrend.STdata
## updateSTdataTrend - EX:with updateTrend.STdata
## internalUpdateTrendFnc - internal function for updateTrend
## internalUpdateTrendSmooth - internal function for updateTrend
## internalCreateTrendFnc - internal function that creates a trend function
##' Function that completes a data matrix using iterative svd as described in
##' Fuentes et. al. (2006). The function iterates between computing the svd for
##' the matrix and replacing the missing values by linear regression of the
##' columns onto the first \code{ncomp} svd components. As initial replacement
##' for the missing values regression on the column averages are used. The
##' function \emph{will fail} if entire rows and/or columns are missing from the
##' data matrix.
##'
##' @title Missing Data SVD
##' @param X Data matrix, with missing values marked by \code{NA}.
##' @param niter Maximum number of iterations to run before exiting,
##' \code{Inf} will run until the \code{conv.reldiff} criteria is met.
##' @param ncomp Number of SVD components to use in the reconstruction (>0).
##' @param conv.reldiff Assume the iterative procedure has converged when
##' the relative difference between two consecutive iterations is less than
##' \code{conv.reldiff}.
##' @return A list with the following components:
##' \item{Xfill}{The completed data matrix with missing values replaced by
##' fitting the data to the \code{ncomp} most important svd
##' components}
##' \item{svd}{The result of svd on the completed data matrix, i.e.
##' \code{svd(Xfill)}}
##' \item{status}{A vector of status variables: \code{diff}, the
##' absolute difference between the two last iterations;
##' \code{rel.diff}, the relative difference;
##' \code{n.iter}, the number of iterations; and
##' \code{max.iter}, the requested maximum number of
##' iterations.}
##'
##' @references
##' M. Fuentes, P. Guttorp, and P. D. Sampson. (2006) Using Transforms to
##' Analyze Space-Time Processes in Statistical methods for spatio-temporal
##' systems (B. Finkenstadt, L. Held, V. Isham eds.) 77-150
##'
##' @author Paul D. Sampson and Johan Lindstrom
##'
##' @example Rd_examples/Ex_SVDmiss.R
##' @family SVD for missing data
##' @family data matrix
##' @export
SVDmiss <- function(X, niter=25, ncomp=min(4,dim(X)[2]), conv.reldiff=0.001)
{
##ensure at least one iteration
niter <- max(niter,1)
##and ensure sane number of components
if( ncomp<1 ){
stop("ncomp should be >0, is: ", ncomp)
}
##First find missing values
Ina <- is.na(X)
if( all(!Ina) ){
##if X has no missing data this is simple
svd0 <- svd(X)
XF <- X
i <- diff <- reldiff <- 0
}else{
##X has missing data, use iterative method
##Iterative svd calculation with missing data.
##Initial first element of U matrix is average curve.
u1 <- rowMeans(X, na.rm = TRUE)
XM <- matrix(1, nrow(X), ncol(X))
XM[Ina] <- 0
XZ <- X
# v1 is proportional to X'u1/(u1'u1), but with calculations
# filling in zeros for missing values in the sums.
XZ[Ina] <- 0.
# Fill in missing values for initial complete SVD calculation.
# Then iterate using only complete data.
v1 <- diag(t(XZ) %*% (XM * u1))/diag(t(XM * u1) %*% (XM * u1))
XF <- X
XF[Ina] <- (matrix(u1, ncol = 1) %*% matrix(v1, nrow = 1))[Ina]
if( any(is.na(XF)) )
stop("Unable to complete matrix, too much missing data")
reldiff <- conv.reldiff+1
i <- 0
while(i<niter && reldiff>conv.reldiff){
svd0 <- svd(XF)
Xnew <- X
Xnew[Ina] <- (svd0$u[, 1:ncomp] %*%
diag(svd0$d[1:ncomp],nrow=length(svd0$d[1:ncomp])) %*%
t(svd0$v[,1:ncomp]))[Ina]
diff <- max(abs(Xnew - XF))
reldiff <- diff/max(abs(XF[Ina]))
XF <- Xnew
i <- i+1
}
}#if( all(!is.na(X)) ) ... else ...
final.diff <- c(diff,reldiff,i,niter)
names(final.diff) <- c("diff","rel.diff","n.iter","max.iter")
return(list(svd=svd0, Xfill=XF, status=final.diff))
}##function SVDmiss
##' Function that computes smooth functions for a data matrix with missing
##' values, as described in Fuentes et. al. (2006), or does cross validation to
##' determine a suitable number of basis functions. The function uses
##' \code{\link{SVDmiss}} to complete the matrix and then computes smooth basis
##' functions by applying \code{\link[stats:smooth.spline]{smooth.spline}} to the
##' SVD of the completed data matrix.
##'
##' \code{SVDsmoothCV} uses leave-one-column-out cross-validation; holding one column
##' out from \code{X}, calling \code{SVDsmooth}, and then regressing
##' the held out column on the resulting smooth functions. Cross-validation
##' statistics computed for each of these regressions include MSE, R-squared,
##' AIC and BIC. The weighted average (weighted by number of observations in the
##' colum) is then reported as CV-statistics.
##'
##' @title Smooth Basis Functions for Data Matrix with Missing Values
##'
##' @param X Data matrix, with missing values marked by \code{NA} (use
##' \code{\link{createDataMatrix}}). Rows and/or columns
##' that are completely missing will be dropped (with a message), for the rows the
##' smooths will be interpolated using
##' \code{\link[stats:predict.smooth.spline]{predict.smooth.spline}}.
##' @param n.basis Number of smooth basis functions to compute, will be passed
##' as \code{ncomp} to \code{\link{SVDmiss}}; for \code{SVDsmoothCV} a
##' vector with the different number of basis functions to evaluate (including 0).
##' @param date.ind Vector giving the observation time of each row in
##' \code{X}, used as \code{x} in \cr
##' \code{\link[stats:smooth.spline]{smooth.spline}} when computing the smooth
##' basis functions. If missing \code{\link{convertCharToDate}} is used to
##' coerce the \code{rownames(X)}.
##' @param scale If \code{TRUE}, will use \code{\link[base:scale]{scale}} to scale
##' \code{X} before calling \code{\link{SVDmiss}}.
##' @param niter,conv.reldiff Controls convergence, passed to \code{\link{SVDmiss}}.
##' @param df,spar The desired degrees of freedom/smoothing parameter for the
##' spline, \cr see \code{\link[stats:smooth.spline]{smooth.spline}}
##' @param fnc If \code{TRUE} return a function instead of the trend-matrix, see
##' Value below.
##' @param ... Additional parameters passed to \code{SVDsmooth}; i.e. \code{date.ind},
##' \code{scale}, \code{niter}, \code{conv.reldiff}, \code{df}, \code{spar},
##' and/or \code{fnc}.
##'
##' @return Depends on the function:
##' \item{SVDsmooth}{A matrix (if \code{fnc==FALSE}) where each column is a
##' smooth basis function based on the SVD of the completed
##' data matrix. The left most column contains the smooth of
##' the most important SVD. If \code{fnc==TRUE} a function
##' that will create the data matrix if called as
##' \code{fnc(date.ind)}, \code{fnc(1:dim(X)[1])}, or
##' \code{fnc(convertCharToDate( rownames(X) ))}.}
##' \item{SVDsmoothCV}{A list of class \code{SVDcv} with components:
##' \describe{
##' \item{CV.stat,CV.sd}{\code{data.frame}s with mean and standard
##' deviation of the CV statistics for each of the number of basis
##' functions evaluated.}
##' \item{MSE.all,R2.all,AIC.all,BIC.all}{\code{data.frame}s with the
##' individual MSE, R2, AIC, and BIC values for each column in the
##' data matrix and for each number of basis functions evaluated.}
##' \item{smoothSVD}{A list with \code{length(n.basis)} components. If
##' \code{fnc==FALSE} each component contains an array
##' where \code{smoothSVD[[j]][,,i]} is the result of
##' \code{SVDsmooth} applied to \code{X[,-i]} with
##' \code{n.basis[j]} smooth functions; if
##' \code{fnc==FALSE} each component contains a list
##' of functions as \code{smoothSVD[[j]][[i]]}.}
##' }
##' }
##'
##' @references
##' M. Fuentes, P. Guttorp, and P. D. Sampson. (2006) Using Transforms to
##' Analyze Space-Time Processes in Statistical methods for spatio-temporal
##' systems (B. Finkenstadt, L. Held, V. Isham eds.) 77-150
##'
##' @author Paul D. Sampson and Johan Lindstrom
##'
##' @example Rd_examples/Ex_SVDsmooth.R
##' @family SVD for missing data
##' @family data matrix
##' @family SVDcv methods
##' @importFrom stats sd
##' @importFrom stats smooth.spline
##' @export
SVDsmooth <- function(X, n.basis=min(2,dim(X)[2]), date.ind=NULL, scale=TRUE,
niter=100, conv.reldiff=0.001, df=NULL, spar=NULL,
fnc=FALSE){
##date.ind contains NA/missing, use rownames of X
if( missing(date.ind) || is.null(date.ind) || any(is.na(date.ind)) ){
date.ind <- convertCharToDate( rownames(X) )
}
if( is.null(date.ind) ){
date.ind <- 1:dim(X)[1]
}
##check length od date.ind
if( length(date.ind)!=dim(X)[1] ){
stop("length(date.ind)!=dim(X)[1]")
}
##check number of columns
if( max(n.basis)>dim(X)[2] ){
stop("Number of basis functions cannot exceed dim(X)[2]")
}
##scale data to zero mean and unit variance
if( scale ){
X <- scale(X)
}
##drop missing columns and rows
Icol <- colSums(is.na(X))!=dim(X)[1]
Irow <- rowSums(is.na(X))!=dim(X)[2]
X.reduced <- X[Irow,Icol]
if( any(!Icol) ){
message( paste("Dropping column(s):", paste(which(!Icol), collapse=", ")) )
}
if( any(!Irow) ){
message( paste("Smooth interpolated at", sum(!Irow), "time(s).") )
}
##missing data SVD
X.svd <- SVDmiss(X.reduced, niter=niter, ncomp=n.basis,
conv.reldiff=conv.reldiff)
##create function that compute trend matrix
##first compute splines for each basis-set
spline <- lapply(1:n.basis,
function(j){
if( is.null(df) && is.null(spar) ){
return(smooth.spline(date.ind[Irow], X.svd$svd$u[,j]))
}else{
return(smooth.spline(date.ind[Irow], X.svd$svd$u[,j],
df=df, spar=spar))
}
})
##second compute mean and scaling for each spline
scale.spline <- lapply(spline, function(x){ c(mean(x$y), sd(x$y)) })
##define the function
trend.fnc <- function(x=date.ind){
X.comps <- matrix(NA, length(x), length(spline))
for(i in 1:length(spline)){
##scale components to unit variance and zero mean
X.comps[,i] <- scale(predict(spline[[i]], as.double(x))$y,
center=scale.spline[[i]][1],
scale=scale.spline[[i]][2])
##Ensure that components have alternating sign
X.comps[,i] <- (-1)^i*X.comps[,i]*sign(X.comps[1,i])
}
##add names
rownames(X.comps) <- as.character(x)
colnames(X.comps) <- paste("V", 1:length(spline), sep="")
return( X.comps )
}
##clean up the environment of the function (reduces overhead)
tmp <- ls(environment(trend.fnc))
tmp <- tmp[ !(tmp %in% c("date.ind", "spline", "scale.spline")) ]
rm(tmp, pos=environment(trend.fnc))
##should we return the function or the matrix?
if( fnc ){
return( trend.fnc )
}
return( trend.fnc() )
}##function SVDsmooth
##' @rdname SVDsmooth
##' @importFrom stats extractAIC
##' @importFrom stats lm
##' @importFrom stats weighted.mean
##' @export
SVDsmoothCV <- function(X, n.basis, ...){
##check number of columns
if( max(n.basis)>(dim(X)[2]-1) )
stop("Number of basis functions cannot exceed dim(X)[2]-1.")
##matrices that holds the individual CV-statistics
CV.all <- array(NA, c(length(n.basis), dim(X)[2], 4))
dimnames(CV.all) <- list(paste("n.basis", as.character(n.basis), sep="."),
colnames(X), c("MSE", "R2", "AIC", "BIC"))
N <- matrix(colSums(!is.na(X)), 1, dim(X)[2])
colnames(N) <- colnames(X)
##loop over the number of basis functions requested
trend <- vector("list", length(n.basis))
trend.fnc <- vector("list", length(n.basis))
names(trend) <- names(trend.fnc) <- dimnames(CV.all)[[1]]
for(i in 1:length(n.basis)){
##for each basis functions do leave one out CV
trend[[i]] <- array(NA, c(dim(X)[1], n.basis[i], dim(X)[2]))
trend.fnc[[i]] <- vector("list", dim(X)[2])
for(j in 1:dim(X)[2]){
if( n.basis[i]!=0 ){
if(i==1 && j==1){
tmp.trend <- SVDsmooth(X[,-j,drop=FALSE], n.basis[i], ...)
}else{
##suppress message given above to avoid noisy CV output
suppressMessages( tmp.trend <- SVDsmooth(X[,-j,drop=FALSE],
n.basis[i], ...))
}
##did we return function or matrix?
if( is.function(tmp.trend) ){
trend.fnc[[i]][[j]] <- tmp.trend
trend[[i]][,,j] <- tmp.trend()
}else{
trend[[i]][,,j] <- tmp.trend
}
}
if( N[j] > (n.basis[i]+1) ){
##compute CV-error
if( n.basis[i]!=0 ){
tmp.lm <- lm(X[,j]~trend[[i]][,,j])
}else{
tmp.lm <- lm(X[,j]~1)
}
##compute CV-statistics
CV.all[i,j,"MSE"] <- mean(tmp.lm$residuals^2)
CV.all[i,j,"R2"] <- summary(tmp.lm)$r.squared
CV.all[i,j,"AIC"] <- extractAIC(tmp.lm)[2]
CV.all[i,j,"BIC"] <- extractAIC(tmp.lm, k=log(N[j]))[2]
}
}
#add names to the trend functions
if( n.basis[i]!=0 ){
dimnames(trend[[i]]) <- list(rownames(X), paste("V", 1:n.basis[i], sep=""),
colnames(X))
}else{
dimnames(trend[[i]]) <- list(rownames(X), NULL, colnames(X))
}
names(trend.fnc[[i]]) <- colnames(X)
}##for(i in 1:length(n.basis))
weights <- N/sum(N)
CV.stat <- apply(CV.all, c(1,3), weighted.mean, w=weights, na.rm=TRUE)
CV.sd <- CV.stat
for(i in 1:dim(CV.sd)[1]){
for(j in 1:dim(CV.sd)[2]){
CV.sd[i,j] <- weighted.mean((CV.all[i,,j] - CV.stat[i,j])^2,
w=weights, na.rm=TRUE)
}
}
CV.sd <- sqrt(CV.sd)
##return a class object
out <- list(CV.stat=as.data.frame(CV.stat),
CV.sd=as.data.frame(CV.sd),
MSE.all=as.data.frame(t(CV.all[,,"MSE"])),
R2.all=as.data.frame(t(CV.all[,,"R2"])),
AIC.all=as.data.frame(t(CV.all[,,"AIC"])),
BIC.all=as.data.frame(t(CV.all[,,"BIC"])),
smoothSVD=switch(is.function(tmp.trend)+1, trend, trend.fnc))
class(out) <- "SVDcv"
return( out )
}##function SVDsmoothCV
##' A front end function for calling \code{\link{SVDsmooth}} (and
##' \code{\link{SVDsmoothCV}}), with either a \code{STdata} object
##' or vectors containing observations, dates and locations.
##'
##' The function uses \code{\link{createDataMatrix}} to create a
##' data matrix which is passed to \code{\link{SVDsmooth}} (and
##' \code{\link{SVDsmoothCV}}). The output can be used as \cr
##' \code{STdata$trend = calcSmoothTrends(...)$trend}, or \cr
##' \code{STdata$trend = calcSmoothTrends(...)$trend.cv[[i]]}.
##' However, it is recommended to use \code{\link{updateTrend.STdata}}.
##'
##' @title Smooth Basis Functions for a STdata Object
##'
##' @param STdata A \code{STdata}/\code{STmodel} data structure containing
##' observations, see \code{\link{mesa.data.raw}}. Use either this or the \code{obs},
##' \code{date}, and \code{ID} inputs.
##' @param obs A vector of observations.
##' @param date A vector of observation times.
##' @param ID A vector of observation locations.
##' @param subset A subset of locations to extract the data matrix for. A warning
##' is given for each name not found in \code{ID}.
##' @param extra.dates Additional dates for which smooth trends should be
##' computed (any duplicates will be removed).
##' @param n.basis Number of basis functions to compute, see
##' \code{\link{SVDsmooth}}.
##' @param cv Also compute smooth functions using leave one out
##' cross-validation, \cr see \code{\link{SVDsmoothCV}}.
##' @param ... Additional parameters passed to \code{\link{SVDsmooth}}
##' and \code{\link{SVDsmoothCV}}; except \code{fnc}, which is always
##' \code{TRUE}.
##'
##' @return Returns a list with
##' \item{trend}{A data.frame containing the smooth trends and the dates.
##' This can be used as the \code{trend} in \code{STdata$trend}.}
##' \item{trend.cv}{If \code{cv==TRUE} a list of data.frames; each one containing
##' the smooth trend obtained when leaving one site out.
##' Similar to \cr \code{SVDsmoothCV(data)$smoothSVD[[1]]}).}
##' \item{trend.fnc,trend.fnc.cv}{Functions that produce the content of the above
##' data.frames, see \code{\link{SVDsmooth}}.}
##'
##' @author Johan Lindstrom and Paul D. Sampson
##'
##' @example Rd_examples/Ex_calcSmoothTrends.R
##' @family SVD for missing data
##' @family STdata
##' @export
calcSmoothTrends <- function(STdata=NULL, obs=STdata$obs$obs,
date=STdata$obs$date, ID=STdata$obs$ID,
subset=NULL, extra.dates=NULL, n.basis=2,
cv=FALSE, ...){
##add extra dates
date <- c(date, extra.dates)
##and expand obs and ID
obs <- c(obs, rep(NA,length(extra.dates)) )
ID <- c(ID, rep(ID[1],length(extra.dates)) )
##create data matrix (this removes duplicated extra.dates)
data <- createDataMatrix(obs=obs, date=date, ID=ID, subset=subset)
##internal function
extractTrend <- function(x){
x <- as.data.frame(x)
x$date <- convertCharToDate( rownames(x) )
rownames(x) <- NULL
return(x)
}
##now let's do SVD
data.comps.fnc <- SVDsmooth(data, n.basis, fnc=TRUE, ...)
data.comps <- extractTrend(data.comps.fnc())
##and cross-validation
if(cv){
##Message of interpolation already displayed above...
suppressMessages( svd.cv <- SVDsmoothCV(data, n.basis, fnc=TRUE, ...) )
svd.fnc <- svd.cv$smoothSVD[[1]]
svd.tmp <- vector("list", length(svd.fnc))
for(i in 1:length(svd.fnc)){
svd.tmp[[i]] <- extractTrend( svd.fnc[[i]]() )
}
names(svd.tmp) <- names(svd.fnc)
}else{
svd.tmp <- NULL
svd.fnc <- NULL
}
return( list(trend=data.comps, trend.cv=svd.tmp, trend.fnc=data.comps.fnc,
trend.fnc.cv=svd.fnc) )
}##function calcSmoothTrends
##########################
## S3-METHODS FOR SVDcv ##
##########################
##' \code{\link[graphics:plot]{plot}} and
##' \code{\link[graphics:boxplot]{boxplot}} methods for class \code{SVDcv}.
##' Plots summary statistics for the cross-validation. Plots include
##' RMSE, R2, BIC, and scatter plots of BIC for each column.
##'
##' @title Plot and Boxplot cross-validation statistics for \code{SVDcv} object
##' @param x \code{SVDcv} object to plot.
##' @param y Which CV-statistic to plot. For pairs \code{"all"} implies
##' \code{"BIC"}.
##' @param pairs \code{TRUE}/\code{FALSE} plot cross-validation statistics,
##' or scatter plot of individual BIC:s.
##' @param sd \code{TRUE}/\code{FALSE} add uncertainty to each CV-statistic.
##' @param ... Additional parameters passed to \code{\link[graphics:plot]{plot}} or
##' \code{\link[graphics:plot]{pairs}}.
##' @return Nothing
##'
##' @examples
##' ##See SVDsmooth example
##'
##' @author Johan Lindstrom
##'
##' @family SVDcv methods
##' @family SVD for missing data
##' @importFrom graphics abline
##' @importFrom graphics lines
##' @importFrom graphics par
##' @importFrom graphics points
##' @method plot SVDcv
##' @export
plot.SVDcv <- function(x, y=c("all","MSE","R2","AIC","BIC"),
pairs=FALSE, sd=FALSE, ...){
stCheckClass(x, "SVDcv", "'x'")
##we have to use y, cast to resonable name
y <- match.arg(y)
##basis
n.basis <- sapply(x$smoothSVD, dim)[2,]
if( !pairs ){
##plot summary of cross-validation statistics
if( y=="all" ){
par(mfrow=c(2,2),mar=c(4,4,.5,.5))
y <- c("MSE", "R2", "AIC", "BIC")
}
for(i in y){
if( sd ){
ylim <- range(x$CV.stat[[i]] + x$CV.sd[[i]],
x$CV.stat[[i]] - x$CV.sd[[i]])
}else{
ylim <- range(x$CV.stat[[i]])
}
args <- internalPlotFixArgs(list(...),
default=list(xlab="n.basis", ylab=i),
add=list(x=n.basis, y=x$CV.stat[[i]],
type="l", ylim=ylim))
do.call(plot, args)
if( sd ){
lines(n.basis, x$CV.stat[[i]] + x$CV.sd[[i]], lty=3)
lines(n.basis, x$CV.stat[[i]] - x$CV.sd[[i]], lty=3)
}
}
}else{
##plot the CV.stats for each column
if( y=="all" ){ y <- "BIC" }
y <- paste(y, "all", sep=".")
pairs(x[[y]], panel=function(x,y){points(x,y); abline(0,1)}, ...)
}
return(invisible())
}##function plot.SVDcv
##' @rdname plot.SVDcv
##' @importFrom graphics boxplot
##' @importFrom graphics par
##' @method boxplot SVDcv
##' @export
boxplot.SVDcv <- function(x, y=c("all","MSE","R2","AIC","BIC"), ...){
stCheckClass(x, "SVDcv", "'x'")
##we have to use y, cast to resonable name
y <- match.arg(y)
##basis
n.basis <- sapply(x$smoothSVD, dim)[2,]
##plot summary of cross-validation statistics
if( y=="all" ){
par(mfrow=c(2,2),mar=c(4,4,.5,.5))
y <- c("MSE", "R2", "AIC", "BIC")
}
for(i in y){
args <- internalPlotFixArgs(list(...),
default=list(ylab=i),
add=list(at=n.basis,
x=x[[paste(i, "all", sep=".")]]))
do.call(boxplot, args)
}
return(invisible())
}##function boxplot.SVDcv
##' \code{\link[base:print]{print}} and \code{\link[base:summary]{summary}}
##' methods for class \code{SVDcv}, prints cross-validation statistics.
##'
##' @title Print details for \code{SVDcv} object
##' @param x \code{SVDcv} object to print information for.
##' @param ... ignored additional arguments.
##' @return Nothing
##'
##' @examples
##' ##See SVDsmooth example
##'
##' @author Johan Lindstrom
##'
##' @family SVDcv methods
##' @family SVD for missing data
##' @method print SVDcv
##' @export
print.SVDcv <- function(x, ...){
stCheckClass(x, "SVDcv", "'x'")
cat("Result of SVDsmoothCV, average of CV-statistics:\n")
print(x$CV.stat)
return(invisible())
}##function print.SVDcv
##' @rdname print.SVDcv
##' @param object \code{SVDcv} object to compute summary for.
##' @method summary SVDcv
##' @export
summary.SVDcv <- function(object, ...){
stCheckClass(object, "SVDcv", "'object'")
cat("Individual MSE:s by column:\n")
print( apply(object$MSE.all, 2, summary) )
cat("\nIndividual R2:s by column:\n")
print( apply(object$R2.all, 2, summary) )
cat("\nIndividual AIC:s by column:\n")
print( apply(object$AIC.all, 2, summary) )
cat("\nIndividual BIC:s by column:\n")
print( apply(object$BIC.all, 2, summary) )
return(invisible())
}##function summary.SVDcv
############################################
## S3-METHODS that update temporal trends ##
############################################
##' Updates/sets the temporal trend for \code{STdata} or \code{STmodel}
##' objects. It also checks that the spatio-temporal covariate exists for all
##' dates in the trend, mainly an issue if \code{extra.dates!=NULL} adds
##' additional times at which to do predictions.
##'
##' If \code{n.basis} is given this will use \code{\link{calcSmoothTrends}} to
##' compute smoothed SVDs of data for use as temporal trends. If \code{fnc} is
##' given, \code{n.basis} is ignored and \code{fnc} should be a function that,
##' given a vector of dates, returns an object that can be coerced to a
##' data.frame with \emph{numeric} temporal trends; recall that an intercept is
##' \strong{always} added.
##'
##' For a \code{STmodel} object the new trend \emph{must have no more} components
##' than the existing trend; if a function is given \code{colnames} of the new
##' trend \emph{must match} those of the existing trend.
##' In both cases the returned \code{STdata} or \code{STmodel} object will have
##' both a \code{$trend} and \code{$trend.fnc} field.
##'
##' Function \code{updateSTdataTrend} is deprecated and will be removed in
##' future versions of the package.
##'
##' @title Update Trend in \code{STdata} or \code{STmodel} Object
##'
##' @param object A \code{STdata} or \code{STmodel} object, see
##' \code{\link{mesa.data.raw}}.
##' @param n.basis number of basis functions for the temporal trend
##' @param extra.dates Additional dates for which smooth trends should be
##' computed (otherwise only those in \code{object$obs$date} are used);
##' \emph{only} for \code{STdata}.
##' @param fnc Function that defines the trend, see Details and Example.
##' @param ... Additional parameters passed to \code{\link{calcSmoothTrends}}.
##'
##' @return Returns a modfied version of the input, with an added/altered
##' trend.
##'
##' @example Rd_examples/Ex_updateTrend.R
##'
##' @author Johan Lindstrom
##' @family STdata functions
##' @family SVD for missing data
##' @method updateTrend STdata
##' @export
updateTrend.STdata <- function(object, n.basis=0, fnc=NULL,
extra.dates=NULL, ...){
##check class belonging
stCheckClass(object, "STdata", name="object")
##is there an existing trend
if( !is.null(object$trend) ){
message("Replacing existing trend.")
}
##dates
if( length(object$obs$date)==0 ){
##special case since c(..., Date) casts to numeric
##(an issue when object$obs$date is empty)
dates <- sort(unique(extra.dates))
}else{
dates <- sort(unique(c(object$obs$date,extra.dates)))
}
if( !is.null(fnc) ){
if( !missing(n.basis) ){
warning("fnc defined, ignoring n.basis")
}
##call function and update the trend
object <- internalUpdateTrendFnc(object, fnc, dates)
}else{
object <- internalUpdateTrendSmooth(object, n.basis, dates, ...)
}
##ensure that we have all dates (mainly checks that added dates are also in ST-covariate.
trend.dates <- object$trend$date
ST.dates <- convertCharToDate( rownames(object$SpatioTemporal) )
if( any(!(trend.dates %in% ST.dates)) && !is.null(ST.dates) ){
stop( paste("Trend dates not found in rownames(object$SpatioTemporal):",
paste(trend.dates[!(trend.dates %in% ST.dates)],
collapse=", ")) )
}
##return modified object
return(object)
}##function updateTrend.STdata
##' @rdname updateTrend.STdata
##' @family STmodel functions
##' @method updateTrend STmodel
##' @export
updateTrend.STmodel <- function(object, n.basis=0, fnc=NULL, ...){
##check class belonging
stCheckClass(object, "STmodel", name="object")
##can only use exactly the existing dates
dates <- object$trend$date
##keep the old trend for comparison
trend.old <- object$trend
F.old <- object$F
if( !is.null(fnc) ){
if( !missing(n.basis) ){
warning("fnc defined, ignoring n.basis")
}
object <- internalUpdateTrendFnc(object, fnc, dates)
}else{
object <- internalUpdateTrendSmooth(object, n.basis, dates, ...)
}
##ensure that we have exatcly the same dates as before
if( !isTRUE(all.equal(dates, object$trend$date)) ){
stop("Dates for the trend have changed (this should NOT happen)")
}
##check that we dont expand the size of the trend
if( dim(trend.old)[2] < dim(object$trend)[2] ){
stop("New trend must have no more components than existing trend.")
}
##and update the F
object$F <- internalSTmodelCreateF(object)
Ind <- match( colnames(object$F), colnames(F.old) )
if( any(is.na(Ind)) ){
stop("colnames of new trend not found in the old: ",
paste(colnames(object$F)[is.na(Ind)], collapse=", "))
}
##Lets pick out components relevant to the new trend (LUR and beta-fields)
object$LUR <- object$LUR[Ind]
object$LUR.all <- object$LUR.all[Ind]
object$LUR.list <- object$LUR.list[Ind]
object$cov.beta$covf <- object$cov.beta$covf[Ind]
object$cov.beta$nugget <- object$cov.beta$nugget[Ind]
##return modified object
return(object)
}##function updateTrend.STdata
##' @rdname updateTrend.STdata
##' @export
updateTrend <- function(object, n.basis=0, fnc=NULL, ...){
UseMethod("updateTrend")
}
##' @rdname updateTrend.STdata
##' @export
updateSTdataTrend <- function(object, n.basis=0, extra.dates=NULL, fnc=NULL, ...){
message("Deprecated, use updateTrend instead. NOTE change in parameter order!")
if( missing(n.basis) ){
return( updateTrend.STdata(object=object, fnc=fnc,
extra.dates=extra.dates, ...) )
}else{
return( updateTrend.STdata(object=object, n.basis=n.basis, fnc=fnc,
extra.dates=extra.dates, ...) )
}
}##function updateSTdataTrend
########################################
## Internal functions for updateTrend ##
########################################
internalUpdateTrendFnc <- function(object, fnc, dates){
##call function
object$trend <- as.data.frame(fnc(dates))
object$trend.fnc <- fnc
##drop rownames
rownames(object$trend) <- NULL
##sanity checks
if( !is.data.frame(object$trend) || dim(object$trend)[1]!=length(dates) ){
stop("fnc(x) should return a data.frame with dim(.)[1]=length(x).")
}
if( any(!sapply(object$trend, is.numeric)) ){
stop("fnc(x) should return a data.frame with 'numeric' columns")
}
if( "date" %in% names(object$trend) ){
stop("Non of the columns of fnc(x) may be named 'date'.")
}
##add a column of dates
object$trend$date <- dates
return( object )
}##function internalUpdateTrendFnc
internalUpdateTrendSmooth <- function(object, n.basis, dates, ...){
if(n.basis==0){
##no temporal trend, just add dates for prediction
object$trend <- data.frame( date=dates )
##constant, just hide the trend.fnc
object$trend.fnc <- NULL
}else{
tmp <- calcSmoothTrends(STdata=object, n.basis=n.basis,
extra.dates=dates, ...)
object$trend <- tmp$trend
object$trend.fnc <- tmp$trend.fnc
}
return( object )
}##function internalUpdateTrendFnc
##############################################
## Internal functions for creation of trend ##
##############################################
##' @importFrom stats spline
internalCreateTrendFnc <- function(trend){
X.trend <- trend
Y.trend <- X.trend[,colnames(X.trend)!="date",drop=FALSE]
X.trend <- X.trend$date
trend.fnc <- function(x){
out <- matrix(NA, length(x), dim(Y.trend)[2])
colnames(out) <- colnames(Y.trend)
for(i in 1:dim(out)[2]){
out[,i] <- spline(X.trend, Y.trend[,i], xout=x)$y
}
return( out )
}
return( trend.fnc )
}##function internalCreateTrendFnc
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.
#######################################################################
## FUNCTIONS THAT COMPUTE SVDs FOR TEMPORAL TRENDS WITH MISSING DATA ##
#######################################################################
##Functions in this file:
## SVDmiss - EX:ok
## SVDsmooth - EX:ok
## SVDsmoothCV - EX:with SVDsmooth
## calcSmoothTrends - EX:ok
## plot.SVDcv - EX:with SVDsmooth
## boxplot.SVDcv - EX:with SVDsmooth
## print.SVDcv - EX:with SVDsmooth
## updateTrend.STdata - EX:ok
## updateTrend.STmodel - EX:with updateTrend.STdata
## updateTrend - EX:with updateTrend.STdata
## updateSTdataTrend - EX:with updateTrend.STdata
## internalUpdateTrendFnc - internal function for updateTrend
## internalUpdateTrendSmooth - internal function for updateTrend
## internalCreateTrendFnc - internal function that creates a trend function
##' Function that completes a data matrix using iterative svd as described in
##' Fuentes et. al. (2006). The function iterates between computing the svd for
##' the matrix and replacing the missing values by linear regression of the
##' columns onto the first \code{ncomp} svd components. As initial replacement
##' for the missing values regression on the column averages are used. The
##' function \emph{will fail} if entire rows and/or columns are missing from the
##' data matrix.
##'
##' @title Missing Data SVD
##' @param X Data matrix, with missing values marked by \code{NA}.
##' @param niter Maximum number of iterations to run before exiting,
##' \code{Inf} will run until the \code{conv.reldiff} criteria is met.
##' @param ncomp Number of SVD components to use in the reconstruction (>0).
##' @param conv.reldiff Assume the iterative procedure has converged when
##' the relative difference between two consecutive iterations is less than
##' \code{conv.reldiff}.
##' @return A list with the following components:
##' \item{Xfill}{The completed data matrix with missing values replaced by
##' fitting the data to the \code{ncomp} most important svd
##' components}
##' \item{svd}{The result of svd on the completed data matrix, i.e.
##' \code{svd(Xfill)}}
##' \item{status}{A vector of status variables: \code{diff}, the
##' absolute difference between the two last iterations;
##' \code{rel.diff}, the relative difference;
##' \code{n.iter}, the number of iterations; and
##' \code{max.iter}, the requested maximum number of
##' iterations.}
##'
##' @references
##' M. Fuentes, P. Guttorp, and P. D. Sampson. (2006) Using Transforms to
##' Analyze Space-Time Processes in Statistical methods for spatio-temporal
##' systems (B. Finkenstadt, L. Held, V. Isham eds.) 77-150
##'
##' @author Paul D. Sampson and Johan Lindstrom
##'
##' @example Rd_examples/Ex_SVDmiss.R
##' @family SVD for missing data
##' @family data matrix
##' @export
SVDmiss <- function(X, niter=25, ncomp=min(4,dim(X)[2]), conv.reldiff=0.001)
{
##ensure at least one iteration
niter <- max(niter,1)
##and ensure sane number of components
if( ncomp<1 ){
stop("ncomp should be >0, is: ", ncomp)
}
##First find missing values
Ina <- is.na(X)
if( all(!Ina) ){
##if X has no missing data this is simple
svd0 <- svd(X)
XF <- X
i <- diff <- reldiff <- 0
}else{
##X has missing data, use iterative method
##Iterative svd calculation with missing data.
##Initial first element of U matrix is average curve.
u1 <- rowMeans(X, na.rm = TRUE)
XM <- matrix(1, nrow(X), ncol(X))
XM[Ina] <- 0
XZ <- X
# v1 is proportional to X'u1/(u1'u1), but with calculations
# filling in zeros for missing values in the sums.
XZ[Ina] <- 0.
# Fill in missing values for initial complete SVD calculation.
# Then iterate using only complete data.
v1 <- diag(t(XZ) %*% (XM * u1))/diag(t(XM * u1) %*% (XM * u1))
XF <- X
XF[Ina] <- (matrix(u1, ncol = 1) %*% matrix(v1, nrow = 1))[Ina]
if( any(is.na(XF)) )
stop("Unable to complete matrix, too much missing data")
reldiff <- conv.reldiff+1
i <- 0
while(i<niter && reldiff>conv.reldiff){
svd0 <- svd(XF)
Xnew <- X
Xnew[Ina] <- (svd0$u[, 1:ncomp] %*%
diag(svd0$d[1:ncomp],nrow=length(svd0$d[1:ncomp])) %*%
t(svd0$v[,1:ncomp]))[Ina]
diff <- max(abs(Xnew - XF))
reldiff <- diff/max(abs(XF[Ina]))
XF <- Xnew
i <- i+1
}
}#if( all(!is.na(X)) ) ... else ...
final.diff <- c(diff,reldiff,i,niter)
names(final.diff) <- c("diff","rel.diff","n.iter","max.iter")
return(list(svd=svd0, Xfill=XF, status=final.diff))
}##function SVDmiss
##' Function that computes smooth functions for a data matrix with missing
##' values, as described in Fuentes et. al. (2006), or does cross validation to
##' determine a suitable number of basis functions. The function uses
##' \code{\link{SVDmiss}} to complete the matrix and then computes smooth basis
##' functions by applying \code{\link[stats:smooth.spline]{smooth.spline}} to the
##' SVD of the completed data matrix.
##'
##' \code{SVDsmoothCV} uses leave-one-column-out cross-validation; holding one column
##' out from \code{X}, calling \code{SVDsmooth}, and then regressing
##' the held out column on the resulting smooth functions. Cross-validation
##' statistics computed for each of these regressions include MSE, R-squared,
##' AIC and BIC. The weighted average (weighted by number of observations in the
##' colum) is then reported as CV-statistics.
##'
##' @title Smooth Basis Functions for Data Matrix with Missing Values
##'
##' @param X Data matrix, with missing values marked by \code{NA} (use
##' \code{\link{createDataMatrix}}). Rows and/or columns
##' that are completely missing will be dropped (with a message), for the rows the
##' smooths will be interpolated using
##' \code{\link[stats:predict.smooth.spline]{predict.smooth.spline}}.
##' @param n.basis Number of smooth basis functions to compute, will be passed
##' as \code{ncomp} to \code{\link{SVDmiss}}; for \code{SVDsmoothCV} a
##' vector with the different number of basis functions to evaluate (including 0).
##' @param date.ind Vector giving the observation time of each row in
##' \code{X}, used as \code{x} in \cr
##' \code{\link[stats:smooth.spline]{smooth.spline}} when computing the smooth
##' basis functions. If missing \code{\link{convertCharToDate}} is used to
##' coerce the \code{rownames(X)}.
##' @param scale If \code{TRUE}, will use \code{\link[base:scale]{scale}} to scale
##' \code{X} before calling \code{\link{SVDmiss}}.
##' @param niter,conv.reldiff Controls convergence, passed to \code{\link{SVDmiss}}.
##' @param df,spar The desired degrees of freedom/smoothing parameter for the
##' spline, \cr see \code{\link[stats:smooth.spline]{smooth.spline}}
##' @param fnc If \code{TRUE} return a function instead of the trend-matrix, see
##' Value below.
##' @param ... Additional parameters passed to \code{SVDsmooth}; i.e. \code{date.ind},
##' \code{scale}, \code{niter}, \code{conv.reldiff}, \code{df}, \code{spar},
##' and/or \code{fnc}.
##'
##' @return Depends on the function:
##' \item{SVDsmooth}{A matrix (if \code{fnc==FALSE}) where each column is a
##' smooth basis function based on the SVD of the completed
##' data matrix. The left most column contains the smooth of
##' the most important SVD. If \code{fnc==TRUE} a function
##' that will create the data matrix if called as
##' \code{fnc(date.ind)}, \code{fnc(1:dim(X)[1])}, or
##' \code{fnc(convertCharToDate( rownames(X) ))}.}
##' \item{SVDsmoothCV}{A list of class \code{SVDcv} with components:
##' \describe{
##' \item{CV.stat,CV.sd}{\code{data.frame}s with mean and standard
##' deviation of the CV statistics for each of the number of basis
##' functions evaluated.}
##' \item{MSE.all,R2.all,AIC.all,BIC.all}{\code{data.frame}s with the
##' individual MSE, R2, AIC, and BIC values for each column in the
##' data matrix and for each number of basis functions evaluated.}
##' \item{smoothSVD}{A list with \code{length(n.basis)} components. If
##' \code{fnc==FALSE} each component contains an array
##' where \code{smoothSVD[[j]][,,i]} is the result of
##' \code{SVDsmooth} applied to \code{X[,-i]} with
##' \code{n.basis[j]} smooth functions; if
##' \code{fnc==FALSE} each component contains a list
##' of functions as \code{smoothSVD[[j]][[i]]}.}
##' }
##' }
##'
##' @references
##' M. Fuentes, P. Guttorp, and P. D. Sampson. (2006) Using Transforms to
##' Analyze Space-Time Processes in Statistical methods for spatio-temporal
##' systems (B. Finkenstadt, L. Held, V. Isham eds.) 77-150
##'
##' @author Paul D. Sampson and Johan Lindstrom
##'
##' @example Rd_examples/Ex_SVDsmooth.R
##' @family SVD for missing data
##' @family data matrix
##' @family SVDcv methods
##' @importFrom stats sd
##' @importFrom stats smooth.spline
##' @export
SVDsmooth <- function(X, n.basis=min(2,dim(X)[2]), date.ind=NULL, scale=TRUE,
niter=100, conv.reldiff=0.001, df=NULL, spar=NULL,
fnc=FALSE){
##date.ind contains NA/missing, use rownames of X
if( missing(date.ind) || is.null(date.ind) || any(is.na(date.ind)) ){
date.ind <- convertCharToDate( rownames(X) )
}
if( is.null(date.ind) ){
date.ind <- 1:dim(X)[1]
}
##check length od date.ind
if( length(date.ind)!=dim(X)[1] ){
stop("length(date.ind)!=dim(X)[1]")
}
##check number of columns
if( max(n.basis)>dim(X)[2] ){
stop("Number of basis functions cannot exceed dim(X)[2]")
}
##scale data to zero mean and unit variance
if( scale ){
X <- scale(X)
}
##drop missing columns and rows
Icol <- colSums(is.na(X))!=dim(X)[1]
Irow <- rowSums(is.na(X))!=dim(X)[2]
X.reduced <- X[Irow,Icol]
if( any(!Icol) ){
message( paste("Dropping column(s):", paste(which(!Icol), collapse=", ")) )
}
if( any(!Irow) ){
message( paste("Smooth interpolated at", sum(!Irow), "time(s).") )
}
##missing data SVD
X.svd <- SVDmiss(X.reduced, niter=niter, ncomp=n.basis,
conv.reldiff=conv.reldiff)
##create function that compute trend matrix
##first compute splines for each basis-set
spline <- lapply(1:n.basis,
function(j){
if( is.null(df) && is.null(spar) ){
return(smooth.spline(date.ind[Irow], X.svd$svd$u[,j]))
}else{
return(smooth.spline(date.ind[Irow], X.svd$svd$u[,j],
df=df, spar=spar))
}
})
##second compute mean and scaling for each spline
scale.spline <- lapply(spline, function(x){ c(mean(x$y), sd(x$y)) })
##define the function
trend.fnc <- function(x=date.ind){
X.comps <- matrix(NA, length(x), length(spline))
for(i in 1:length(spline)){
##scale components to unit variance and zero mean
X.comps[,i] <- scale(predict(spline[[i]], as.double(x))$y,
center=scale.spline[[i]][1],
scale=scale.spline[[i]][2])
##Ensure that components have alternating sign
X.comps[,i] <- (-1)^i*X.comps[,i]*sign(X.comps[1,i])
}
##add names
rownames(X.comps) <- as.character(x)
colnames(X.comps) <- paste("V", 1:length(spline), sep="")
return( X.comps )
}
##clean up the environment of the function (reduces overhead)
tmp <- ls(environment(trend.fnc))
tmp <- tmp[ !(tmp %in% c("date.ind", "spline", "scale.spline")) ]
rm(tmp, pos=environment(trend.fnc))
##should we return the function or the matrix?
if( fnc ){
return( trend.fnc )
}
return( trend.fnc() )
}##function SVDsmooth
##' @rdname SVDsmooth
##' @importFrom stats extractAIC
##' @importFrom stats lm
##' @importFrom stats weighted.mean
##' @export
SVDsmoothCV <- function(X, n.basis, ...){
##check number of columns
if( max(n.basis)>(dim(X)[2]-1) )
stop("Number of basis functions cannot exceed dim(X)[2]-1.")
##matrices that holds the individual CV-statistics
CV.all <- array(NA, c(length(n.basis), dim(X)[2], 4))
dimnames(CV.all) <- list(paste("n.basis", as.character(n.basis), sep="."),
colnames(X), c("MSE", "R2", "AIC", "BIC"))
N <- matrix(colSums(!is.na(X)), 1, dim(X)[2])
colnames(N) <- colnames(X)
##loop over the number of basis functions requested
trend <- vector("list", length(n.basis))
trend.fnc <- vector("list", length(n.basis))
names(trend) <- names(trend.fnc) <- dimnames(CV.all)[[1]]
for(i in 1:length(n.basis)){
##for each basis functions do leave one out CV
trend[[i]] <- array(NA, c(dim(X)[1], n.basis[i], dim(X)[2]))
trend.fnc[[i]] <- vector("list", dim(X)[2])
for(j in 1:dim(X)[2]){
if( n.basis[i]!=0 ){
if(i==1 && j==1){
tmp.trend <- SVDsmooth(X[,-j,drop=FALSE], n.basis[i], ...)
}else{
##suppress message given above to avoid noisy CV output
suppressMessages( tmp.trend <- SVDsmooth(X[,-j,drop=FALSE],
n.basis[i], ...))
}
##did we return function or matrix?
if( is.function(tmp.trend) ){
trend.fnc[[i]][[j]] <- tmp.trend
trend[[i]][,,j] <- tmp.trend()
}else{
trend[[i]][,,j] <- tmp.trend
}
}
if( N[j] > (n.basis[i]+1) ){
##compute CV-error
if( n.basis[i]!=0 ){
tmp.lm <- lm(X[,j]~trend[[i]][,,j])
}else{
tmp.lm <- lm(X[,j]~1)
}
##compute CV-statistics
CV.all[i,j,"MSE"] <- mean(tmp.lm$residuals^2)
CV.all[i,j,"R2"] <- summary(tmp.lm)$r.squared
CV.all[i,j,"AIC"] <- extractAIC(tmp.lm)[2]
CV.all[i,j,"BIC"] <- extractAIC(tmp.lm, k=log(N[j]))[2]
}
}
#add names to the trend functions
if( n.basis[i]!=0 ){
dimnames(trend[[i]]) <- list(rownames(X), paste("V", 1:n.basis[i], sep=""),
colnames(X))
}else{
dimnames(trend[[i]]) <- list(rownames(X), NULL, colnames(X))
}
names(trend.fnc[[i]]) <- colnames(X)
}##for(i in 1:length(n.basis))
weights <- N/sum(N)
CV.stat <- apply(CV.all, c(1,3), weighted.mean, w=weights, na.rm=TRUE)
CV.sd <- CV.stat
for(i in 1:dim(CV.sd)[1]){
for(j in 1:dim(CV.sd)[2]){
CV.sd[i,j] <- weighted.mean((CV.all[i,,j] - CV.stat[i,j])^2,
w=weights, na.rm=TRUE)
}
}
CV.sd <- sqrt(CV.sd)
##return a class object
out <- list(CV.stat=as.data.frame(CV.stat),
CV.sd=as.data.frame(CV.sd),
MSE.all=as.data.frame(t(CV.all[,,"MSE"])),
R2.all=as.data.frame(t(CV.all[,,"R2"])),
AIC.all=as.data.frame(t(CV.all[,,"AIC"])),
BIC.all=as.data.frame(t(CV.all[,,"BIC"])),
smoothSVD=switch(is.function(tmp.trend)+1, trend, trend.fnc))
class(out) <- "SVDcv"
return( out )
}##function SVDsmoothCV
##' A front end function for calling \code{\link{SVDsmooth}} (and
##' \code{\link{SVDsmoothCV}}), with either a \code{STdata} object
##' or vectors containing observations, dates and locations.
##'
##' The function uses \code{\link{createDataMatrix}} to create a
##' data matrix which is passed to \code{\link{SVDsmooth}} (and
##' \code{\link{SVDsmoothCV}}). The output can be used as \cr
##' \code{STdata$trend = calcSmoothTrends(...)$trend}, or \cr
##' \code{STdata$trend = calcSmoothTrends(...)$trend.cv[[i]]}.
##' However, it is recommended to use \code{\link{updateTrend.STdata}}.
##'
##' @title Smooth Basis Functions for a STdata Object
##'
##' @param STdata A \code{STdata}/\code{STmodel} data structure containing
##' observations, see \code{\link{mesa.data.raw}}. Use either this or the \code{obs},
##' \code{date}, and \code{ID} inputs.
##' @param obs A vector of observations.
##' @param date A vector of observation times.
##' @param ID A vector of observation locations.
##' @param subset A subset of locations to extract the data matrix for. A warning
##' is given for each name not found in \code{ID}.
##' @param extra.dates Additional dates for which smooth trends should be
##' computed (any duplicates will be removed).
##' @param n.basis Number of basis functions to compute, see
##' \code{\link{SVDsmooth}}.
##' @param cv Also compute smooth functions using leave one out
##' cross-validation, \cr see \code{\link{SVDsmoothCV}}.
##' @param ... Additional parameters passed to \code{\link{SVDsmooth}}
##' and \code{\link{SVDsmoothCV}}; except \code{fnc}, which is always
##' \code{TRUE}.
##'
##' @return Returns a list with
##' \item{trend}{A data.frame containing the smooth trends and the dates.
##' This can be used as the \code{trend} in \code{STdata$trend}.}
##' \item{trend.cv}{If \code{cv==TRUE} a list of data.frames; each one containing
##' the smooth trend obtained when leaving one site out.
##' Similar to \cr \code{SVDsmoothCV(data)$smoothSVD[[1]]}).}
##' \item{trend.fnc,trend.fnc.cv}{Functions that produce the content of the above
##' data.frames, see \code{\link{SVDsmooth}}.}
##'
##' @author Johan Lindstrom and Paul D. Sampson
##'
##' @example Rd_examples/Ex_calcSmoothTrends.R
##' @family SVD for missing data
##' @family STdata
##' @export
calcSmoothTrends <- function(STdata=NULL, obs=STdata$obs$obs,
date=STdata$obs$date, ID=STdata$obs$ID,
subset=NULL, extra.dates=NULL, n.basis=2,
cv=FALSE, ...){
##add extra dates
date <- c(date, extra.dates)
##and expand obs and ID
obs <- c(obs, rep(NA,length(extra.dates)) )
ID <- c(ID, rep(ID[1],length(extra.dates)) )
##create data matrix (this removes duplicated extra.dates)
data <- createDataMatrix(obs=obs, date=date, ID=ID, subset=subset)
##internal function
extractTrend <- function(x){
x <- as.data.frame(x)
x$date <- convertCharToDate( rownames(x) )
rownames(x) <- NULL
return(x)
}
##now let's do SVD
data.comps.fnc <- SVDsmooth(data, n.basis, fnc=TRUE, ...)
data.comps <- extractTrend(data.comps.fnc())
##and cross-validation
if(cv){
##Message of interpolation already displayed above...
suppressMessages( svd.cv <- SVDsmoothCV(data, n.basis, fnc=TRUE, ...) )
svd.fnc <- svd.cv$smoothSVD[[1]]
svd.tmp <- vector("list", length(svd.fnc))
for(i in 1:length(svd.fnc)){
svd.tmp[[i]] <- extractTrend( svd.fnc[[i]]() )
}
names(svd.tmp) <- names(svd.fnc)
}else{
svd.tmp <- NULL
svd.fnc <- NULL
}
return( list(trend=data.comps, trend.cv=svd.tmp, trend.fnc=data.comps.fnc,
trend.fnc.cv=svd.fnc) )
}##function calcSmoothTrends
##########################
## S3-METHODS FOR SVDcv ##
##########################
##' \code{\link[graphics:plot]{plot}} and
##' \code{\link[graphics:boxplot]{boxplot}} methods for class \code{SVDcv}.
##' Plots summary statistics for the cross-validation. Plots include
##' RMSE, R2, BIC, and scatter plots of BIC for each column.
##'
##' @title Plot and Boxplot cross-validation statistics for \code{SVDcv} object
##' @param x \code{SVDcv} object to plot.
##' @param y Which CV-statistic to plot. For pairs \code{"all"} implies
##' \code{"BIC"}.
##' @param pairs \code{TRUE}/\code{FALSE} plot cross-validation statistics,
##' or scatter plot of individual BIC:s.
##' @param sd \code{TRUE}/\code{FALSE} add uncertainty to each CV-statistic.
##' @param ... Additional parameters passed to \code{\link[graphics:plot]{plot}} or
##' \code{\link[graphics:plot]{pairs}}.
##' @return Nothing
##'
##' @examples
##' ##See SVDsmooth example
##'
##' @author Johan Lindstrom
##'
##' @family SVDcv methods
##' @family SVD for missing data
##' @importFrom graphics abline
##' @importFrom graphics lines
##' @importFrom graphics par
##' @importFrom graphics points
##' @method plot SVDcv
##' @export
plot.SVDcv <- function(x, y=c("all","MSE","R2","AIC","BIC"),
pairs=FALSE, sd=FALSE, ...){
stCheckClass(x, "SVDcv", "'x'")
##we have to use y, cast to resonable name
y <- match.arg(y)
##basis
n.basis <- sapply(x$smoothSVD, dim)[2,]
if( !pairs ){
##plot summary of cross-validation statistics
if( y=="all" ){
par(mfrow=c(2,2),mar=c(4,4,.5,.5))
y <- c("MSE", "R2", "AIC", "BIC")
}
for(i in y){
if( sd ){
ylim <- range(x$CV.stat[[i]] + x$CV.sd[[i]],
x$CV.stat[[i]] - x$CV.sd[[i]])
}else{
ylim <- range(x$CV.stat[[i]])
}
args <- internalPlotFixArgs(list(...),
default=list(xlab="n.basis", ylab=i),
add=list(x=n.basis, y=x$CV.stat[[i]],
type="l", ylim=ylim))
do.call(plot, args)
if( sd ){
lines(n.basis, x$CV.stat[[i]] + x$CV.sd[[i]], lty=3)
lines(n.basis, x$CV.stat[[i]] - x$CV.sd[[i]], lty=3)
}
}
}else{
##plot the CV.stats for each column
if( y=="all" ){ y <- "BIC" }
y <- paste(y, "all", sep=".")
pairs(x[[y]], panel=function(x,y){points(x,y); abline(0,1)}, ...)
}
return(invisible())
}##function plot.SVDcv
##' @rdname plot.SVDcv
##' @importFrom graphics boxplot
##' @importFrom graphics par
##' @method boxplot SVDcv
##' @export
boxplot.SVDcv <- function(x, y=c("all","MSE","R2","AIC","BIC"), ...){
stCheckClass(x, "SVDcv", "'x'")
##we have to use y, cast to resonable name
y <- match.arg(y)
##basis
n.basis <- sapply(x$smoothSVD, dim)[2,]
##plot summary of cross-validation statistics
if( y=="all" ){
par(mfrow=c(2,2),mar=c(4,4,.5,.5))
y <- c("MSE", "R2", "AIC", "BIC")
}
for(i in y){
args <- internalPlotFixArgs(list(...),
default=list(ylab=i),
add=list(at=n.basis,
x=x[[paste(i, "all", sep=".")]]))
do.call(boxplot, args)
}
return(invisible())
}##function boxplot.SVDcv
##' \code{\link[base:print]{print}} and \code{\link[base:summary]{summary}}
##' methods for class \code{SVDcv}, prints cross-validation statistics.
##'
##' @title Print details for \code{SVDcv} object
##' @param x \code{SVDcv} object to print information for.
##' @param ... ignored additional arguments.
##' @return Nothing
##'
##' @examples
##' ##See SVDsmooth example
##'
##' @author Johan Lindstrom
##'
##' @family SVDcv methods
##' @family SVD for missing data
##' @method print SVDcv
##' @export
print.SVDcv <- function(x, ...){
stCheckClass(x, "SVDcv", "'x'")
cat("Result of SVDsmoothCV, average of CV-statistics:\n")
print(x$CV.stat)
return(invisible())
}##function print.SVDcv
##' @rdname print.SVDcv
##' @param object \code{SVDcv} object to compute summary for.
##' @method summary SVDcv
##' @export
summary.SVDcv <- function(object, ...){
stCheckClass(object, "SVDcv", "'object'")
cat("Individual MSE:s by column:\n")
print( apply(object$MSE.all, 2, summary) )
cat("\nIndividual R2:s by column:\n")
print( apply(object$R2.all, 2, summary) )
cat("\nIndividual AIC:s by column:\n")
print( apply(object$AIC.all, 2, summary) )
cat("\nIndividual BIC:s by column:\n")
print( apply(object$BIC.all, 2, summary) )
return(invisible())
}##function summary.SVDcv
############################################
## S3-METHODS that update temporal trends ##
############################################
##' Updates/sets the temporal trend for \code{STdata} or \code{STmodel}
##' objects. It also checks that the spatio-temporal covariate exists for all
##' dates in the trend, mainly an issue if \code{extra.dates!=NULL} adds
##' additional times at which to do predictions.
##'
##' If \code{n.basis} is given this will use \code{\link{calcSmoothTrends}} to
##' compute smoothed SVDs of data for use as temporal trends. If \code{fnc} is
##' given, \code{n.basis} is ignored and \code{fnc} should be a function that,
##' given a vector of dates, returns an object that can be coerced to a
##' data.frame with \emph{numeric} temporal trends; recall that an intercept is
##' \strong{always} added.
##'
##' For a \code{STmodel} object the new trend \emph{must have no more} components
##' than the existing trend; if a function is given \code{colnames} of the new
##' trend \emph{must match} those of the existing trend.
##' In both cases the returned \code{STdata} or \code{STmodel} object will have
##' both a \code{$trend} and \code{$trend.fnc} field.
##'
##' Function \code{updateSTdataTrend} is deprecated and will be removed in
##' future versions of the package.
##'
##' @title Update Trend in \code{STdata} or \code{STmodel} Object
##'
##' @param object A \code{STdata} or \code{STmodel} object, see
##' \code{\link{mesa.data.raw}}.
##' @param n.basis number of basis functions for the temporal trend
##' @param extra.dates Additional dates for which smooth trends should be
##' computed (otherwise only those in \code{object$obs$date} are used);
##' \emph{only} for \code{STdata}.
##' @param fnc Function that defines the trend, see Details and Example.
##' @param ... Additional parameters passed to \code{\link{calcSmoothTrends}}.
##'
##' @return Returns a modfied version of the input, with an added/altered
##' trend.
##'
##' @example Rd_examples/Ex_updateTrend.R
##'
##' @author Johan Lindstrom
##' @family STdata functions
##' @family SVD for missing data
##' @method updateTrend STdata
##' @export
updateTrend.STdata <- function(object, n.basis=0, fnc=NULL,
extra.dates=NULL, ...){
##check class belonging
stCheckClass(object, "STdata", name="object")
##is there an existing trend
if( !is.null(object$trend) ){
message("Replacing existing trend.")
}
##dates
if( length(object$obs$date)==0 ){
##special case since c(..., Date) casts to numeric
##(an issue when object$obs$date is empty)
dates <- sort(unique(extra.dates))
}else{
dates <- sort(unique(c(object$obs$date,extra.dates)))
}
if( !is.null(fnc) ){
if( !missing(n.basis) ){
warning("fnc defined, ignoring n.basis")
}
##call function and update the trend
object <- internalUpdateTrendFnc(object, fnc, dates)
}else{
object <- internalUpdateTrendSmooth(object, n.basis, dates, ...)
}
##ensure that we have all dates (mainly checks that added dates are also in ST-covariate.
trend.dates <- object$trend$date
ST.dates <- convertCharToDate( rownames(object$SpatioTemporal) )
if( any(!(trend.dates %in% ST.dates)) && !is.null(ST.dates) ){
stop( paste("Trend dates not found in rownames(object$SpatioTemporal):",
paste(trend.dates[!(trend.dates %in% ST.dates)],
collapse=", ")) )
}
##return modified object
return(object)
}##function updateTrend.STdata
##' @rdname updateTrend.STdata
##' @family STmodel functions
##' @method updateTrend STmodel
##' @export
updateTrend.STmodel <- function(object, n.basis=0, fnc=NULL, ...){
##check class belonging
stCheckClass(object, "STmodel", name="object")
##can only use exactly the existing dates
dates <- object$trend$date
##keep the old trend for comparison
trend.old <- object$trend
F.old <- object$F
if( !is.null(fnc) ){
if( !missing(n.basis) ){
warning("fnc defined, ignoring n.basis")
}
object <- internalUpdateTrendFnc(object, fnc, dates)
}else{
object <- internalUpdateTrendSmooth(object, n.basis, dates, ...)
}
##ensure that we have exatcly the same dates as before
if( !isTRUE(all.equal(dates, object$trend$date)) ){
stop("Dates for the trend have changed (this should NOT happen)")
}
##check that we dont expand the size of the trend
if( dim(trend.old)[2] < dim(object$trend)[2] ){
stop("New trend must have no more components than existing trend.")
}
##and update the F
object$F <- internalSTmodelCreateF(object)
Ind <- match( colnames(object$F), colnames(F.old) )
if( any(is.na(Ind)) ){
stop("colnames of new trend not found in the old: ",
paste(colnames(object$F)[is.na(Ind)], collapse=", "))
}
##Lets pick out components relevant to the new trend (LUR and beta-fields)
object$LUR <- object$LUR[Ind]
object$LUR.all <- object$LUR.all[Ind]
object$LUR.list <- object$LUR.list[Ind]
object$cov.beta$covf <- object$cov.beta$covf[Ind]
object$cov.beta$nugget <- object$cov.beta$nugget[Ind]
##return modified object
return(object)
}##function updateTrend.STdata
##' @rdname updateTrend.STdata
##' @export
updateTrend <- function(object, n.basis=0, fnc=NULL, ...){
UseMethod("updateTrend")
}
##' @rdname updateTrend.STdata
##' @export
updateSTdataTrend <- function(object, n.basis=0, extra.dates=NULL, fnc=NULL, ...){
message("Deprecated, use updateTrend instead. NOTE change in parameter order!")
if( missing(n.basis) ){
return( updateTrend.STdata(object=object, fnc=fnc,
extra.dates=extra.dates, ...) )
}else{
return( updateTrend.STdata(object=object, n.basis=n.basis, fnc=fnc,
extra.dates=extra.dates, ...) )
}
}##function updateSTdataTrend
########################################
## Internal functions for updateTrend ##
########################################
internalUpdateTrendFnc <- function(object, fnc, dates){
##call function
object$trend <- as.data.frame(fnc(dates))
object$trend.fnc <- fnc
##drop rownames
rownames(object$trend) <- NULL
##sanity checks
if( !is.data.frame(object$trend) || dim(object$trend)[1]!=length(dates) ){
stop("fnc(x) should return a data.frame with dim(.)[1]=length(x).")
}
if( any(!sapply(object$trend, is.numeric)) ){
stop("fnc(x) should return a data.frame with 'numeric' columns")
}
if( "date" %in% names(object$trend) ){
stop("Non of the columns of fnc(x) may be named 'date'.")
}
##add a column of dates
object$trend$date <- dates
return( object )
}##function internalUpdateTrendFnc
internalUpdateTrendSmooth <- function(object, n.basis, dates, ...){
if(n.basis==0){
##no temporal trend, just add dates for prediction
object$trend <- data.frame( date=dates )
##constant, just hide the trend.fnc
object$trend.fnc <- NULL
}else{
tmp <- calcSmoothTrends(STdata=object, n.basis=n.basis,
extra.dates=dates, ...)
object$trend <- tmp$trend
object$trend.fnc <- tmp$trend.fnc
}
return( object )
}##function internalUpdateTrendFnc
##############################################
## Internal functions for creation of trend ##
##############################################
##' @importFrom stats spline
internalCreateTrendFnc <- function(trend){
X.trend <- trend
Y.trend <- X.trend[,colnames(X.trend)!="date",drop=FALSE]
X.trend <- X.trend$date
trend.fnc <- function(x){
out <- matrix(NA, length(x), dim(Y.trend)[2])
colnames(out) <- colnames(Y.trend)
for(i in 1:dim(out)[2]){
out[,i] <- spline(X.trend, Y.trend[,i], xout=x)$y
}
return( out )
}
return( trend.fnc )
}##function internalCreateTrendFnc
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.