#' Prediction based on DS or CCA model
#'
#' Apply an empirical-statistical downscaling model to new data
#'
#' \code{predict} is similar to the predict function in R
#'
#' \code{project} returns projection of climate
#'
#' @aliases predict.ds predict.ds.eof predict.ds.comb predict.mvr
#' predict.ds.pca predict.ds.station predict.cca project project.ds
#'
#' @param x A ds object
#' @param newdata An eof object containing the new data sets on which the prediction is made.
#' @param addnoise If TRUE, will add an attribute called "noise" to the ouput based on WG
#' @param n Number of runs to be generated, used only if addnoise is set to TRUE
#'
#' @return Predicted ds values.
#' @seealso \code{\link{DS}}
#'
#' @examples
#'
#' # Get predictor
#' ## Get reanalysis
#' X <- t2m.DNMI(lon=c(-40,50),lat=c(40,75))
#' ## Get Gcm output
#' Y <- t2m.NorESM.M(lon=c(-40,50),lat=c(40,75))
#' ## Combine
#' XY <- combine(X,Y)
#' # Compute common eof for January
#' ceof <- EOF(XY,it='jan')
#' # Get predictand
#' data(Oslo)
#' # Do the downscaling
#' ds <- DS(Oslo,ceof)
#' # Plot ds results
#' plot(ds)
#' # Do the prediction based on the calibration (or the fitted values)
#' ds.pre <- predict(ds)
#' # Plot predicted results based on ds object
#' plot(ds.pre)
#' # Display the attribute "aspect"
#' attr(ds.pre, "aspect")
#' ## Extract the projected results
#' plot(project(ds))
#'
#' @exportS3Method
#' @export
predict.ds <- function(x,...,newdata=NULL,addnoise=FALSE,n=100,verbose=FALSE) {
if (verbose) print(paste("predict.ds",paste(class(x),collapse='-')))
stopifnot(!missing(x),inherits(x,"ds"))
if ( inherits(x,c('eof','comb')) & (is.null(newdata) | is.logical(newdata)) ) {
if(verbose) print("no new predictor data is provided")
if(verbose) print("predictand is an EOF")
if (inherits(x,'comb')) {
if(verbose) print("predictand is a combined EOF")
y <- predict.ds.comb(x,newdata=newdata,addnoise=addnoise,
n=n,verbose=verbose)
} else {
if(verbose) print("predictand is an EOF, but use PCA call")
y <- predict.ds.pca(x,newdata=newdata,addnoise=addnoise,
n=n,verbose=verbose)
}
} else if (inherits(x,'field')) {
if (verbose) print("predictand is a field object")
if (inherits(x,'station')) {
y <- predict.ds.eof(x,newdata=newdata,addnoise=addnoise,
n=n,verbose=verbose)
} else {
y <- predict.ds.pca(x,newdata=newdata,addnoise=addnoise,
n=n,verbose=verbose)
}
} else if (inherits(x,'eof')) {
if (verbose) print("predictand is an EOF")
if(verbose) print("new predictor data is provided")
if (inherits(newdata,'comb')) {
if (verbose) print("new predictor data is a combined EOF")
y <- predict.ds.comb(x,newdata=newdata,addnoise=addnoise,
n=n,verbose=verbose)
} else if (inherits(newdata,'eof')) {
if (verbose) print("new predictor data is an EOF")
# y <- predict.ds.eof(x,newdata=newdata,addnoise=addnoise,
# n=n,verbose=verbose)
y <- predict.ds.pca(x,newdata=newdata,addnoise=addnoise,
n=n,verbose=verbose)
}
} else if (inherits(x,'pca')) {
if(verbose) print("predictand is PCA")
y <- predict.ds.pca(x,newdata=newdata,addnoise=addnoise,n=n,verbose=verbose)
} else if (inherits(x,'station')) {
y <- predict.ds.station(x,newdata=newdata,addnoise=addnoise,n=n,verbose=verbose)
}
y <- attrcp(attr(x,'original_data'),y)
attr(y,'history') <- history.stamp(x)
class(y) <- class(x)
invisible(y)
}
#' @export predict.ds.station
predict.ds.station <- function(x,...,newdata=NULL,addnoise=FALSE,n=100,verbose=FALSE) {
if (verbose) print("predict.ds.pca")
if(verbose) print(paste(class(x),collapse='-'))
if (is.null(names(newdata))) names(newdata) <- paste('X',1:dim(newdata)[2],sep='.')
t <- index(as.eof(newdata))
newdata <- data.frame(coredata(as.eof(newdata)))
## Fudge fix for when the names of the data are x.1, x.2,,, instead of X.1, X.2,..
if (nchar(names(newdata))[1]!=3) names(newdata) <- paste('X',1:dim(newdata)[2],sep='.')
names(newdata) <- toupper(names(newdata))
d <- dim(newdata)
if (is.null(d)) d <- c(length(x),1)
model <- attr(x,'model')
y <- predict(model,newdata)
y <- zoo(y, order.by=t)
# KMP 2018-12-30: Changed to include model in attrcp. Is there any reason not to copy model to y?
y <- attrcp(x, y, ignore=c("name", "n.apps", "appendix", "dimnames","aspect"))
class(y) <- class(x)[-1]
invisible(y)
}
#' @export predict.ds.eof
predict.ds.eof <- function(x,...,newdata=NULL,addnoise=FALSE,n=100,verbose=FALSE) {
if (verbose) print("predict.ds.eof")
stopifnot(!missing(x),inherits(x,"ds"))
if (verbose) print(paste(paste(class(x),collapse='-')))
X <- as.eof(x)
if (verbose) print(paste(class(X),collapse='-'))
#print(dim(X))
if (is.null(newdata)) neofs <- length(attr(X,'eigenvalues')) else {
if (inherits(newdata,'eof'))
neofs <- length(attr(newdata,'eigenvalues')) else {
if (!is.null(dim(newdata))) neofs <- dim(newdata)[2] else
neofs <- 1
}
}
# For some reason, the column names of newdata is muddled here,
# and hence Xnames is used to enforce names 'X.1', 'X.2', 'X.3', ...
Xnames <- paste("X.",1:neofs,sep="")
if (verbose) print(Xnames)
if (is.null(newdata)) {
if (verbose) print('Use calibration data')
newdata <- data.frame(X=coredata(X))
src <- attr(X,'source')
idx <- index(X)
} else {
if (verbose) print('Use new data')
idx <- index(newdata)
src <- attr(newdata,'source')
newdata <- as.data.frame(newdata)
if ((length(attr(X,'eigenvalues'))) != neofs)
warning(paste('Warning: newdata and X have different number of EOFs:',length(attr(X,'eigenvalues')),neofs))
}
#print(summary(newdata))
names(newdata) <- Xnames
model <- attr(x,'model')
if (verbose) {print('names(newdata):'); print(Xnames)}
## AM 04-04-2015 model is always a list object - Quick fix here ...
## KMP 19-11-2015 the if (!is.list(model)) solution
## does not work when there is only one model
## which is a list of coefficients, residuals, ...
if (is.model(model,verbose=verbose)) {
# if (names(model)[1]=="coefficients") { ## REB 2016-01-12: changed to the line above.
y <- predict(model,newdata=newdata) + attr(x,'mean')
} else {
# if (!is.null(newdata)) {
y <- lapply(model,predict,newdata)
# } else {
# y <- lapply(model,predict)
# }
y <- matrix(unlist(y),nrow=length(idx),ncol=length(model))
}
## predict - phase scramble of residual
## There is a bug in the following lines, works only if model is not a list object -- need fixes here ...
residual <- model$residuals
if (addnoise) {
if (verbose) print('add noise')
l <- length(index(x))
noise <- matrix(rep(NA,n*l),n,l)
for (i in 1:n)
noise[i,] <- FTscramble(noise)
noise <- zoo(t(noise),order.by=index(x))
attr(y,'noise') <- noise
}
y <- zoo(y,order.by=idx)
attr(y,'source') <- src
if (!is.null(residual)) attr(y,'residual.mean') <- mean(residual,na.rm=TRUE)
if (!is.null(residual)) attr(y,'residual.sd') <- sd(residual,na.rm=TRUE)
class(y) <- class(x)[-1] ## AM remove 'ds' from output class
# KMP 2018-12-30: Changed to include model in attrcp. Is there any reason not to copy model to y?
y <- attrcp(x, y, ignore=c("name", "n.apps", "appendix", "dimnames","aspect"))
if (verbose) print('predict.ds.eof complete')
invisible(y)
}
#' @export predict.ds.pca
predict.ds.pca <- function(x,...,newdata=NULL,addnoise=FALSE,n=100,verbose=FALSE) {
if (verbose) print("predict.ds.pca")
if(verbose) print(paste(class(x),collapse="-"))
if (is.null(newdata)) {
newdata <- data.frame(coredata(as.eof(x)))
t <- index(as.eof(x))
} else {
if (is.null(names(newdata))) names(newdata) <- paste('X',1:dim(newdata)[2],sep='.')
t <- index(as.eof(newdata))
newdata <- data.frame(coredata(as.eof(newdata)))
}
## Fudge fix for when the names of the data are x.1, x.2,,, instead of X.1, X.2,..
if (nchar(names(newdata))[1]!=3) names(newdata) <- paste('X',1:dim(newdata)[2],sep='.')
names(newdata) <- toupper(names(newdata))
d <- dim(newdata)
if (is.null(d)) d <- c(length(x),1)
model <- attr(x,'model')
y <- lapply(model,predict,newdata)
y <- matrix(unlist(y),nrow=d[1],ncol=length(model))
## Replace
y <- zoo(y, order.by=t)
y <- attrcp(x, y, ignore=c("name", "n.apps", "appendix", "dimnames","aspect"))
class(y) <- class(x)[-1]
invisible(y)
}
#' @export predict.ds.comb
predict.ds.comb <- function(x,...,newdata=NULL,addnoise=FALSE,n=100,verbose=FALSE) {
if (verbose) print("predict.ds.comb")
stopifnot(!missing(x),inherits(x,"ds"))
## If newdata is set as NULL, reassign it to FALSE
if (is.null(newdata)) newdata <- FALSE
if (is.logical(newdata)) {
if (!newdata) {
X <- attr(x,'eof')
} else {
X <- attr(x,'appendix.1')
}
} else {
X <- newdata ## newdata must be an eof object
}
if (is.null(X)) neofs <- 1 else neofs <- dim(X)[2]
model <- attr(x,'model')
## If newdata is provided as a zoo object or an EOF:
if (inherits(newdata,'zoo')) Y <- predict.ds.eof(x=x,newdata=newdata,verbose=verbose)
## In newdata is set as TRUE or FALSE -------------------------------------------------
if (is.logical(newdata)) {
# For some reason, the column names of newdata are muddled here,
# and hence Xnames is used to enforce names 'X.1', 'X.2', 'X.3', ...
#Xnames <- paste("X.",1:neofs,sep="")
if (newdata) {
## newdata <- data.frame(X=coredata(X))
n.app <- attr(x,'n.apps')
} else {
## X <- newdata
n.app <- 1
}
rownm <- rep("",n.app+1)
rownm[1] <- attr(x,'source')
if (!newdata) {
Y <- predict.ds.eof(x=x,verbose=verbose)
attr(Y,'aspect') <- 'fitted'
} else {
for (i in 1:n.app) {
X <- attr(x,paste('appendix.',i,sep=""))
if (inherits(x,c('pca','eof')))
y <- predict.ds.eof(x=x,newdata=X,verbose=verbose) else
y <- X
if (i==1)
Y <- y
else
Y <- merge(y,Y,all=TRUE)
}
attr(Y,'aspect') <- 'predicted'
}
} ## End of the section dealing with predict if newdata is TRUE/FALSE ----------------
residual <- model$residuals
if (addnoise) {
l <- length(residual)
noise <- matrix(rep(NA,n*l),n,l)
for (i in 1:n)
noise[i,] <- FTscramble(noise)
attr(Y,'noise') <- noise
names(Y) <- rownm
}
Y <- zoo(Y,order.by=index(Y))
attr(Y,'source') <- attr(X,'source')
attr(Y,'residual.mean') <- mean(residual,na.rm=TRUE)
attr(Y,'residual.sd') <- sd(residual,na.rm=TRUE)
## KMP 2018-12-30: Changed to include model in attrcp. Is there any reason not to copy model to y?
Y <- attrcp(x,Y,ignore = c("name", "n.apps", "appendix", "dimnames","aspect"))
attr(Y,'aspect') <- 'predicted'
invisible(Y)
}
#' @export
project <- function(x,...) {
UseMethod("project")
}
#' @exportS3Method
#' @export project.ds
project.ds <- function(x,...,newdata=NULL,addnoise=FALSE,n=100,verbose=FALSE) {
if (verbose) print("project.ds")
stopifnot(!missing(x),inherits(x,"ds"))
if (is.null(newdata)) {
X <- attr(x,'eof')
} else {
X <- newdata ## newdata must be an eof object
}
neofs <- length(attr(X,'eigenvalues'))
# For some reason, the column names of newdata are muddled here,
# and hence Xnames is used to enforce names 'X.1', 'X.2', 'X.3', ...
Xnames <- paste("X.",1:neofs,sep="")
if (verbose) print(Xnames)
if (is.null(newdata)) {
n.app <- attr(X,'n.apps')
} else {
## X <- newdata
n.app <- attr(newdata,'n.apps')
}
model <- attr(x,'model')
if (verbose) {print(summary(model)); print(n.app)}
rownm <- rep("",n.app+1)
rownm[1] <- attr(x,'source')
for (i in 1:n.app) {
if(is.null(newdata)) newdata <- attr(X,paste('appendix.',i,sep=""))
if(verbose) {
print("Data for obs:")
print(summary(newdata))
}
names(newdata) <- Xnames
newdata <- attrcp(X,newdata)
y <- predict.ds.eof(x=x,newdata=newdata,addnoise=FALSE,n=100,verbose=verbose)
names(newdata) <- Xnames
if (verbose) {
print("Data for GCM:")
print(summary(newdata))
}
if (verbose) {
print("DS values:")
print(summary(predict(model,newdata=newdata)))
}
if (i==1) {
Y <- y
} else {
Y <- merge(y,Y,all=TRUE)
}
}
residual <- model$residuals
if (addnoise) {
l <- length(residual)
noise <- matrix(rep(NA,n*l),n,l)
for (i in 1:n)
noise[i,] <- FTscramble(noise)
attr(Y,'noise') <- noise
}
Y <- zoo(Y,order.by=index(Y))
#print(dim(Y)); print(rownm)
#print(names(Y))
names(Y) <- rownm
attr(Y,'source') <- attr(X,'source')
attr(Y,'residual.mean') <- mean(residual,na.rm=TRUE)
attr(Y,'residual.sd') <- sd(residual,na.rm=TRUE)
# KMP 2018-12-30: Changed to include model in attrcp. Is there any reason not to copy model to y?
y <- attrcp(x, y, ignore=c("name", "n.apps", "appendix", "dimnames","aspect"))
attr(Y,'aspect') <- 'projected'
Y <- as.station(Y)
invisible(Y)
}
# To get one predictor pattern, use predict with newdata set to
# a vector where most variables are set to zero apart from one
# variable set to unity for the identification of teleconnection pattern.
#' @exportS3Method
#' @export
predict.mvr <- function(object, ..., newdata=NULL, verbose=FALSE) {
if(verbose) print("predict.mvr")
if (is.null(newdata)) newdata <- object$data
x <- newdata
psi <- object$model
Z <- object$fitted.values
if (inherits(newdata,'zoo')) {
Yhat <- zoo(coredata(x) %*% psi,order.by=index(x))
} else if (is.vector(x)) {
Yhat <- t(x) %*% psi
}
nattr <- softattr(Z)
for (i in 1:length(nattr)) {
attr(Yhat,nattr[i]) <- attr(Z,nattr[i])
}
class(Yhat) <- class(object$fitted.values)
if (inherits(Yhat,'eof')) Yhat <- eof2field(Yhat)
#attr(Yhat,'history') <- c('predict.MVR',attr(Z,'history'))
#attr(Yhat,'date-stamp') <- date()
#attr(Yhat,'call') <- match.call()
attr(Yhat,'history') <- history.stamp(x)
invisible(Yhat)
}
#' @exportS3Method
#' @export
predict.cca <- function(x, ..., newdata=NULL, verbose=FALSE) {
if(verbose) print("predict.cca")
if (!is.null(newdata)) X <- newdata else X <- x$X
Y.hat <- Psi(x) %*% X
Y.hat <- attrcp(Y.hat,X)
return(Y.hat)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.