[Analogue-commits] r122 - in pkg: R man

[noreply at r-forge.r-project.org]

Author: gsimpson
Date: 2009-04-26 16:06:22 +0200 (Sun, 26 Apr 2009)
New Revision: 122

Added:
   pkg/R/stdError.R
   pkg/man/stdError.Rd
Log:
Adds stdError() to compute weighted standard deviations of environment values for the k-closest analogues in MAT models.

Added: pkg/R/stdError.R
===================================================================
--- pkg/R/stdError.R	                        (rev 0)
+++ pkg/R/stdError.R	2009-04-26 14:06:22 UTC (rev 122)
@@ -0,0 +1,85 @@
+## compute the standard error of MAT reconstructed values
+## following the ideas of ter Braak 1995 to use the weighted
+## variance of the k-closest analogues
+
+`stdError` <- function(object, ...) {
+    UseMethod("stdError")
+}
+
+`stdError.mat` <- function(object, k, ...) {
+    getOrd <- function(dis) {
+        nas <- is.na(dis)
+        order(dis[!nas])
+    }
+    getWts <- function(i, dis, ords, k.seq) {
+        nas <- is.na(dis[,i])
+        dis[!nas, i][ords[,i]][k.seq]
+    }
+    getEnv <- function(i, dis, ords, k.seq, y){
+        nas <- is.na(dis[,i])
+        y[!nas][ords[,i]][k.seq]
+    }
+    if(missing(k)) {
+        k <- getK(object, ...)
+        auto <- object$auto
+    } else {
+        auto <- FALSE
+    }
+    ## create k sequence
+    k.seq <- seq_len(k)
+    ## ordering of objects in terms of dissim
+    ords <- apply(object$Dij, 2, getOrd)
+    SEQ <- seq_len(ncol(ords))
+    ## weights = 1/Dij
+    wts <- 1 / sapply(SEQ, getWts, object$Dij, ords, k.seq)
+    ## produce matrix of Env data for each site
+    env <- sapply(SEQ, getEnv, object$Dij, ords, k.seq,
+                  object$orig.y)
+    ## mean of env of k closest analogues
+    ybar <- colMeans(env)
+    wtdSD <- sqrt(colSums(wts * sweep(env, 2, ybar, "-")^2) /
+                  colSums(wts))
+    names(wtdSD) <- names(object$orig.y)
+    class(wtdSD) <- "stdError"
+    return(wtdSD)
+}
+
+`stdError.predict.mat` <- function(object, k, ...) {
+    getOrd <- function(dis) {
+        nas <- is.na(dis)
+        order(dis[!nas])
+    }
+    getWts <- function(i, dis, ords, k.seq) {
+        nas <- is.na(dis[,i])
+        dis[!nas, i][ords[,i]][k.seq]
+    }
+    getEnv <- function(i, dis, ords, k.seq, y){
+        nas <- is.na(dis[,i])
+        y[!nas][ords[,i]][k.seq]
+    }
+    if(missing(k)) {
+        k <- getK(object, ...)
+        auto <- object$auto
+    } else {
+        auto <- FALSE
+    }
+    ## create k sequence
+    k.seq <- seq_len(k)
+    ## ordering of objects in terms of dissim
+    ords <- apply(object$Dij, 2, getOrd)
+    SEQ <- seq_len(ncol(ords))
+    ## weights = 1/Dij
+    wts <- 1 / sapply(SEQ, getWts, object$Dij, ords, k.seq)
+    ## produce matrix of Env data for each site
+    env <- sapply(SEQ, getEnv, object$Dij, ords, k.seq,
+                  object$observed)
+    ## mean of env of k closest analogues
+    ybar <- colMeans(env)
+    wtdSD <- sqrt(colSums(wts * sweep(env, 2, ybar, "-")^2) /
+                  colSums(wts))
+    names(wtdSD) <- colnames(object$predictions$model$predicted)
+    class(wtdSD) <- "stdError"
+    attr(wtdSD, "k") <- k
+    attr(wtdSD, "auto") <- object$auto
+    return(wtdSD)
+}

Added: pkg/man/stdError.Rd
===================================================================
--- pkg/man/stdError.Rd	                        (rev 0)
+++ pkg/man/stdError.Rd	2009-04-26 14:06:22 UTC (rev 122)
@@ -0,0 +1,71 @@
+\name{stdError}
+\alias{stdError}
+\alias{stdError.mat}
+\alias{stdError.predict.mat}
+\title{Standard error of MAT fitted and predicted values}
+\description{
+  Computes the weighted standard deviation of the environment for the
+  \emph{k}-closest analogues for each sample. This measure is proposed
+  as a measure of reconstruction uncertainty for MAT models.
+}
+\usage{
+stdError(object, ...)
+
+\method{stdError}{mat}(object, k, ...)
+
+\method{stdError}{mat}(object, k, ...)
+}
+\arguments{
+  \item{object}{Object for which the uncertainty measure is to be
+    computed. Currently methods for \code{\link{mat}} and
+    \code{\link{predict.mat}}.}
+  \item{k}{numeric; how many analogues to take? If missing, the default,
+    \code{k} is chosen using \code{\link{getK}}.}
+  \item{\dots}{Additional arguments passed to other methods. Currently
+    not used.}
+}
+%\details{
+%  TODO
+%}
+\value{
+  A named numeric vector of weighted standard deviations of the
+  environment for the \emph{k} closest analogues used to compute the MAT
+  predicted values.
+
+  The returned vector has attributes \code{"k"} and \code{"auto"},
+  indicating the number of analogues used and whether this was
+  determined from \code{object} or supplied by the user.
+}
+%\references{ ~put references to the literature/web site here ~ }
+\author{Gavin L. Simpson}
+%\note{ ~~further notes~~ 
+%
+%}
+\seealso{\code{\link{minDC}}, \code{\link{mat}},
+  \code{\link{predict.mat}}.}
+\examples{
+## Imbrie and Kipp Sea Surface Temperature
+data(ImbrieKipp)
+data(SumSST)
+data(V12.122)
+
+## merge training set and core samples
+dat <- join(ImbrieKipp, V12.122, verbose = TRUE)
+
+## extract the merged data sets and convert to proportions
+ImbrieKipp <- dat[[1]] / 100
+ImbrieKippCore <- dat[[2]] / 100
+
+## fit the MAT model using the squared chord distance measure
+ik.mat <- mat(ImbrieKipp, SumSST, method = "SQchord")
+
+## standard errors
+stdError(ik.mat)
+
+## reconstruct for the V12-122 core data
+coreV12.mat <- predict(ik.mat, V12.122, k = 3)
+## standard errors
+stdError(coreV12.mat)
+}
+\keyword{methods}
+\keyword{univar}
\ No newline at end of file