[Dplr-commits] r813 - in pkg/dplR: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Apr 12 19:36:14 CEST 2014
Author: mvkorpel
Date: 2014-04-12 19:36:14 +0200 (Sat, 12 Apr 2014)
New Revision: 813
Modified:
pkg/dplR/DESCRIPTION
pkg/dplR/NAMESPACE
pkg/dplR/R/detrend.R
pkg/dplR/R/detrend.series.R
pkg/dplR/R/helpers.R
pkg/dplR/man/detrend.series.Rd
Log:
* Progress with the 'verbose' option of detrend.series(). Output
somewhat reformatted.
* detrend.series() has new option 'return.info' for including similar
information as printed with 'verbose = TRUE' in the return value.
This option is not (yet) implemented in detrend(). Andy: Currently
there is more information (about uncertainty of model coefficients)
returned with 'return.info = TRUE' than printed with 'verbose =
TRUE'.
Modified: pkg/dplR/DESCRIPTION
===================================================================
--- pkg/dplR/DESCRIPTION 2014-04-12 02:41:53 UTC (rev 812)
+++ pkg/dplR/DESCRIPTION 2014-04-12 17:36:14 UTC (rev 813)
@@ -3,7 +3,7 @@
Type: Package
Title: Dendrochronology Program Library in R
Version: 1.6.0
-Date: 2014-04-11
+Date: 2014-04-12
Authors at R: c(person("Andy", "Bunn", role = c("aut", "cph",
"cre", "trl"), email = "andy.bunn at wwu.edu"), person("Mikko",
"Korpela", role = c("aut", "trl")), person("Franco", "Biondi",
Modified: pkg/dplR/NAMESPACE
===================================================================
--- pkg/dplR/NAMESPACE 2014-04-12 02:41:53 UTC (rev 812)
+++ pkg/dplR/NAMESPACE 2014-04-12 17:36:14 UTC (rev 813)
@@ -42,4 +42,4 @@
S3method(print, redfit)
S3method(plot, rwl)
S3method(plot, crn)
-S3method(summary, rwl)
\ No newline at end of file
+S3method(summary, rwl)
Modified: pkg/dplR/R/detrend.R
===================================================================
--- pkg/dplR/R/detrend.R 2014-04-12 02:41:53 UTC (rev 812)
+++ pkg/dplR/R/detrend.R 2014-04-12 17:36:14 UTC (rev 813)
@@ -16,7 +16,7 @@
stop("'rwl' must be a data.frame")
rn <- row.names(rwl)
- if(!make.plot &&
+ if(!make.plot && !verbose &&
("Spline" %in% method2 || "ModNegExp" %in% method2) &&
!inherits(try(suppressWarnings(req.it <-
requireNamespace("iterators",
Modified: pkg/dplR/R/detrend.series.R
===================================================================
--- pkg/dplR/R/detrend.series.R 2014-04-12 02:41:53 UTC (rev 812)
+++ pkg/dplR/R/detrend.series.R 2014-04-12 17:36:14 UTC (rev 813)
@@ -3,34 +3,45 @@
method = c("Spline", "ModNegExp", "Mean", "Ar"),
nyrs = NULL, f = 0.5, pos.slope = FALSE,
constrain.modnegexp = c("never", "when.fail", "always"),
- verbose=FALSE)
+ verbose = FALSE, return.info = FALSE)
{
stopifnot(identical(make.plot, TRUE) || identical(make.plot, FALSE),
- identical(pos.slope, FALSE) || identical(pos.slope, TRUE))
+ identical(pos.slope, FALSE) || identical(pos.slope, TRUE),
+ identical(verbose, TRUE) || identical(verbose, FALSE),
+ identical(return.info, TRUE) || identical(return.info, FALSE))
known.methods <- c("Spline", "ModNegExp", "Mean", "Ar")
constrain2 <- match.arg(constrain.modnegexp)
method2 <- match.arg(arg = method,
choices = known.methods,
several.ok = TRUE)
-
-
-
- if(verbose){
- nyrs.tmp <- ifelse(test=is.null(nyrs),yes="NULL",nyrs)
- cat("\nVerbose output: ", y.name,
- "\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~",
- "\n Options",
- "\n make.plot =", make.plot,
- "\n method(s) =", paste(deparse(method2), sep = "\n",
- collapse = "\n"),
- "\n nyrs =", nyrs.tmp,
- "\n f =", f,
- "\n pos.slope =", pos.slope,
- "\n constrain.modnegexp =", constrain2,
- "\n verbose =", verbose,
- "\n")
- }
-
+ if (verbose) {
+ widthOpt <- getOption("width")
+ indentSize <- 1
+ indent <- function(x) {
+ paste0(paste0(rep.int(" ", indentSize), collapse = ""), x)
+ }
+ sepLine <-
+ indent(paste0(rep.int("~", max(1, widthOpt - 2 * indentSize)),
+ collapse = ""))
+ cat(gettext("Verbose output: ", domain="R-dplR"), y.name, "\n",
+ sep = "")
+ opts <- c("make.plot" = make.plot,
+ "method(s)" = deparse(method2),
+ "nyrs" = if (is.null(nyrs)) "NULL" else nyrs,
+ "f" = f,
+ "pos.slope" = pos.slope,
+ "constrain.modnegexp" = constrain2,
+ "verbose" = verbose,
+ "return.info" = return.info)
+ optNames <- names(opts)
+ optChar <- c(gettext("Options", domain="R-dplR"),
+ paste(str_pad(optNames,
+ width = max(nchar(optNames)),
+ side = "right"),
+ opts, sep = " "))
+ cat(sepLine, indent(optChar), sep = "\n")
+ }
+
## Remove NA from the data (they will be reinserted later)
good.y <- which(!is.na(y))
if(length(good.y) == 0) {
@@ -40,18 +51,40 @@
}
y2 <- y[good.y]
nY2 <- length(y2)
+
## Recode any zero values to 0.001
- if(verbose) {
- cat("\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
- if(any(y2==0)){
- cat("\n Zeros in series:")
- cat("\n ", names(y2)[y2==0])
- }
- else cat("\n Zeros in series: 0 \n")
+ if (verbose || return.info) {
+ years <- names(y2)
+ if (is.null(years)) {
+ years <- good.y
+ }
+ zeroFun <- function(x) list(zero.years = years[is.finite(x) & x == 0])
+ nFun <- function(x) list(n.zeros = length(x[[1]]))
+ zero.years.data <- zeroFun(y2)
+ n.zeros.data <- nFun(zero.years.data)
+ dataStats <- c(n.zeros.data, zero.years.data)
+ if (verbose) {
+ cat("", sepLine, sep = "\n")
+ if (n.zeros.data[[1]] > 0){
+ if (is.character(years)) {
+ cat(indent(gettext("Zero years in input series:\n",
+ domain="R-dplR")))
+ } else {
+ cat(indent(gettext("Zero indices in input series:\n",
+ domain="R-dplR")))
+ }
+ cat(indent(paste(zero.years.data[[1]], collapse = " ")),
+ "\n", sep = "")
+ } else {
+ cat(indent(gettext("No zeros in input series.\n",
+ domain="R-dplR")))
+ }
+ }
}
y2[y2 == 0] <- 0.001
resids <- list()
+ stats <- list()
if("ModNegExp" %in% method2){
## Nec or lm
@@ -60,9 +93,10 @@
a <- mean(Y[seq_len(max(1, floor(nY * 0.1)))])
b <- -0.01
k <- mean(Y[floor(nY * 0.9):nY])
- nlsForm <- Y ~ a * exp(b * seq_len(nY)) + k
+ nlsForm <- Y ~ I(a * exp(b * seq_along(Y)) + k)
nlsStart <- list(a=a, b=b, k=k)
checked <- FALSE
+ constrained <- FALSE
## Note: nls() may signal an error
if (constrain == "never") {
nec <- nls(formula = nlsForm, start = nlsStart)
@@ -71,6 +105,7 @@
lower = c(a=0, b=-Inf, k=0),
upper = c(a=Inf, b=0, k=Inf),
algorithm = "port")
+ constrained <- TRUE
} else {
nec <- nls(formula = nlsForm, start = nlsStart)
coefs <- coef(nec)
@@ -85,6 +120,7 @@
lower = c(a=0, b=-Inf, k=0),
upper = c(a=Inf, b=0, k=Inf),
algorithm = "port")
+ constrained <- TRUE
}
}
if (!checked) {
@@ -100,39 +136,49 @@
return(NULL)
}
}
- fits
+ tmpFormula <- nlsForm
+ formEnv <- new.env(parent = environment(detrend.series))
+ formEnv[["Y"]] <- Y
+ formEnv[["a"]] <- coefs["a"]
+ formEnv[["b"]] <- coefs["b"]
+ formEnv[["k"]] <- coefs["k"]
+ environment(tmpFormula) <- formEnv
+ structure(fits, constrained = constrained,
+ formula = tmpFormula, summary = summary(nec))
}
ModNegExp <- try(nec.func(y2, constrain2), silent=TRUE)
mneNotPositive <- is.null(ModNegExp)
- if(verbose) {
- cat("\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~",
- "\n Detrend by ModNegExp.",
- "\n Trying to fit nls model...")
- if(class(ModNegExp) != "try-error"){
- cat("\n nls fit but nec.func doesn't return coefs. just fits!")
- cat("\n nls coefs",
- "\n a: goes here",
- "\n b: goes here",
- "\n k: goes here",
- "\n")
- }
+
+ if (verbose) {
+ cat("", sepLine, sep = "\n")
+ cat(indent(gettext("Detrend by ModNegExp.\n", domain = "R-dplR")))
+ cat(indent(gettext("Trying to fit nls model...\n",
+ domain = "R-dplR")))
}
-
if (mneNotPositive || class(ModNegExp) == "try-error") {
- if(verbose) { cat("\n nls failed...fitting linear model...") }
- ## Straight line via linear regression
- if (mneNotPositive) {
- warning("Fits from ModNegExp are not all positive, see constrain.modnegexp argument in detrend.series")
- }
- tm <- cbind(1, seq_len(nY2))
- lm1 <- lm.fit(tm, y2)
- coefs <- lm1[["coefficients"]]
- if(verbose) {
- cat("\n linear model fit",
- "\n Intercept: ", round(coefs[1],4),
- "\n Slope: ", round(coefs[2],4))
- }
- if (all(is.finite(coefs)) && (coefs[2] <= 0 || pos.slope)) {
+ if (verbose) {
+ cat(indent(gettext("nls failed... fitting linear model...",
+ domain = "R-dplR")))
+ }
+ ## Straight line via linear regression
+ if (mneNotPositive) {
+ warning("Fits from ModNegExp are not all positive, see constrain.modnegexp argument in detrend.series")
+ }
+ x <- seq_len(nY2)
+ lm1 <- lm(y2 ~ x)
+ coefs <- coef(lm1)
+ xIdx <- names(coefs) == "x"
+ coefs <- c(coefs[!xIdx], coefs[xIdx])
+ if (verbose) {
+ cat(indent(c(gettext("Linear model fit", domain = "R-dplR"),
+ gettextf("Intercept: %s", format(coefs[1]),
+ domain = "R-dplR"),
+ gettextf("Slope: %s", format(coefs[2]),
+ domain = "R-dplR"))),
+ sep = "\n")
+ }
+ if (all(is.finite(coefs)) && (coefs[2] <= 0 || pos.slope)) {
+ tm <- cbind(1, x)
ModNegExp <- drop(tm %*% coefs)
useMean <- !isTRUE(ModNegExp[1] > 0 &&
ModNegExp[nY2] > 0)
@@ -143,14 +189,40 @@
useMean <- TRUE
}
if (useMean) {
- if(verbose) { cat("\n lm has a positive slope",
- "\n pos.slope = FALSE",
- "\n Detrend by mean.",
- "\n Mean =",mean(y2)) }
- ModNegExp <- rep.int(mean(y2), nY2)
+ theMean <- mean(y2)
+ if (verbose) {
+ cat(indent(c(gettext("lm has a positive slope",
+ "pos.slope = FALSE",
+ "Detrend by mean.",
+ domain = "R-dplR"),
+ gettextf("Mean = %s", format(theMean),
+ domain = "R-dplR"))),
+ sep = "\n")
+ }
+ ModNegExp <- rep.int(theMean, nY2)
+ mneStats <- list(method = "Mean", mean = theMean)
+ } else {
+ mneStats <- list(method = "Line", coefs = coef(summary(lm1)))
}
+ } else if (verbose || return.info) {
+ mneSummary <- attr(ModNegExp, "summary")
+ mneCoefs <- mneSummary[["coefficients"]]
+ mneCoefsE <- mneCoefs[, 1]
+ if (verbose) {
+ cat(indent(c(gettext("nls coefs", domain = "R-dplR"),
+ paste0(names(mneCoefsE), ": ",
+ format(mneCoefsE)))),
+ sep = "\n")
+ }
+ mneStats <- list(method = "ModNegExp",
+ is.constrained = attr(ModNegExp, "constrained"),
+ formula = attr(ModNegExp, "formula"),
+ coefs = mneCoefs)
+ } else {
+ mneStats <- NULL
}
resids$ModNegExp <- y2 / ModNegExp
+ stats$ModNegExp <- mneStats
do.mne <- TRUE
} else {
do.mne <- FALSE
@@ -165,18 +237,24 @@
nyrs2 <- floor(nY2 * 0.67)
else
nyrs2 <- nyrs
- if(verbose) {
- cat("\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~",
- "\n Detrend by spline.",
- "\n Spline parameters",
- "\n nyrs =", nyrs2," f =", f)
- }
+ if (verbose) {
+ cat("", sepLine, sep = "\n")
+ cat(indent(c(gettext(c("Detrend by spline.",
+ "Spline parameters"), domain = "R-dplR"),
+ paste0("nyrs = ", nyrs2, ", f = ", f))),
+ sep = "\n")
+ }
Spline <- ffcsaps(y=y2, x=seq_len(nY2), nyrs=nyrs2, f=f)
if (any(Spline <= 0)) {
warning("Spline fit is not all positive")
- Spline <- rep.int(mean(y2), nY2)
+ theMean <- mean(y2)
+ Spline <- rep.int(theMean, nY2)
+ splineStats <- list(method = "Mean", mean = theMean)
+ } else {
+ splineStats <- list(method = "Spline", nyrs = nyrs2)
}
resids$Spline <- y2 / Spline
+ stats$Spline <- splineStats
do.spline <- TRUE
} else {
do.spline <- FALSE
@@ -184,27 +262,33 @@
if("Mean" %in% method2){
## Fit a horiz line
- Mean <- rep.int(mean(y2), nY2)
- if(verbose) {
- cat("\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~",
- "\n Detrend by mean.",
- "\n Mean =",mean(y2))
+ theMean <- mean(y2)
+ Mean <- rep.int(theMean, nY2)
+ if (verbose) {
+ cat("", sepLine, sep = "\n")
+ cat(indent(c(gettext("Detrend by mean.", domain = "R-dplR"),
+ paste("Mean = ", format(theMean)))),
+ sep = "\n")
}
+ meanStats <- list(method = "Mean", mean = theMean)
resids$Mean <- y2 / Mean
+ stats$Mean <- meanStats
do.mean <- TRUE
} else {
do.mean <- FALSE
}
if("Ar" %in% method2){
## Fit an ar model - aka prewhiten
- if(verbose) {
- ar.tmp <- ar(y2)
- cat("\n ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~",
- "\n Detrend by prewhitening.")
- print(ar.tmp)
+ Ar <- ar.func(y2, model = TRUE)
+ arModel <- attr(Ar, "model")
+ if (verbose) {
+ cat("", sepLine, sep = "\n")
+ cat(indent(gettext("Detrend by prewhitening.", domain = "R-dplR")))
+ print(arModel)
}
- Ar <- ar.func(y2)
- # This will propogate NA to rwi as a result of detrending.
+ arStats <- list(method = "Ar", order = arModel[["order"]],
+ ar = arModel[["ar"]])
+ # This will propogate NA to rwi as a result of detrending.
# Other methods don't. Problem when interacting with other
# methods?
# Also, this can (and does!) produce negative RWI values.
@@ -214,32 +298,48 @@
Ar[Ar<0] <- 0
}
resids$Ar <- Ar / mean(Ar,na.rm=TRUE)
+ stats$Ar <- arStats
do.ar <- TRUE
} else {
do.ar <- FALSE
}
-
+
resids <- data.frame(resids)
+ if (verbose || return.info) {
+ zero.years <- lapply(resids, zeroFun)
+ n.zeros <- lapply(zero.years, nFun)
+ stats <- mapply(c, stats, n.zeros, zero.years)
+ if (verbose) {
+ n.zeros2 <- unlist(n.zeros, use.names = FALSE)
+ zeroFlag <- n.zeros2 > 0
+ methodNames <- names(stats)
+ if (any(zeroFlag)) {
+ cat("", sepLine, sep = "\n")
+ for (i in which(zeroFlag)) {
+ if (is.character(years)) {
+ cat(indent(gettextf("Zero years in %s series:\n",
+ methodNames[i], domain="R-dplR")))
+ } else {
+ cat(indent(gettextf("Zero indices in %s series:\n",
+ methodNames[i], domain="R-dplR")))
+ }
+ cat(indent(paste(zero.years[[i]][[1]], collapse = " ")),
+ "\n", sep = "")
+ }
+ }
+ }
+ }
if(make.plot){
op <- par(no.readonly=TRUE)
on.exit(par(op))
par(mar=c(2.1, 2.1, 2.1, 2.1), mgp=c(1.1, 0.1, 0),
tcl=0.5, xaxs='i')
- n.plots <- 1 + ncol(resids)
- if(n.plots == 5){
- mat <- matrix(c(1,1,2,3,4,5), nrow=3, ncol=2,byrow=TRUE)
- }
- if(n.plots == 4){
- mat <- matrix(c(1,2,3,4), nrow=2, ncol=2,byrow=TRUE)
- }
- if(n.plots == 3){
- mat <- matrix(c(1,1,2,3), nrow=2, ncol=2,byrow=TRUE)
- }
- if(n.plots == 2){
- mat <- matrix(c(1,2), nrow=2, ncol=1,byrow=TRUE)
- }
-
+ mat <- switch(ncol(resids),
+ matrix(c(1,2), nrow=2, ncol=1, byrow=TRUE),
+ matrix(c(1,1,2,3), nrow=2, ncol=2, byrow=TRUE),
+ matrix(c(1,2,3,4), nrow=2, ncol=2, byrow=TRUE),
+ matrix(c(1,1,2,3,4,5), nrow=3, ncol=2, byrow=TRUE))
layout(mat,
widths=rep.int(0.5, ncol(mat)),
heights=rep.int(1, nrow(mat)))
@@ -250,7 +350,7 @@
if(do.spline) lines(Spline, col="green", lwd=2)
if(do.mne) lines(ModNegExp, col="red", lwd=2)
if(do.mean) lines(Mean, col="blue", lwd=2)
-
+
if(do.spline){
plot(resids$Spline, type="l", col="green",
main=gettext("Spline", domain="R-dplR"),
@@ -296,5 +396,10 @@
resids2 <- resids2[, method2]
## Make sure names (years) are included if there is only one method
if(!is.data.frame(resids2)) names(resids2) <- names(y)
- resids2
+ if (return.info) {
+ list(series = resids2,
+ model.info = stats[method2], data.info = dataStats)
+ } else {
+ resids2
+ }
}
Modified: pkg/dplR/R/helpers.R
===================================================================
--- pkg/dplR/R/helpers.R 2014-04-12 02:41:53 UTC (rev 812)
+++ pkg/dplR/R/helpers.R 2014-04-12 17:36:14 UTC (rev 813)
@@ -82,12 +82,16 @@
}
### AR function for chron, normalize1, normalize.xdate, ...
-ar.func <- function(x) {
+ar.func <- function(x, model = FALSE) {
y <- x
idx.goody <- !is.na(y)
ar1 <- ar(y[idx.goody])
y[idx.goody] <- ar1$resid+ar1$x.mean
- y
+ if (isTRUE(model)) {
+ structure(y, model = ar1)
+ } else {
+ y
+ }
}
### Range of years. Used in cms, rcs, rwl.stats, seg.plot, spag.plot, ...
Modified: pkg/dplR/man/detrend.series.Rd
===================================================================
--- pkg/dplR/man/detrend.series.Rd 2014-04-12 02:41:53 UTC (rev 812)
+++ pkg/dplR/man/detrend.series.Rd 2014-04-12 17:36:14 UTC (rev 813)
@@ -10,7 +10,7 @@
method = c("Spline", "ModNegExp", "Mean", "Ar"),
nyrs = NULL, f = 0.5, pos.slope = FALSE,
constrain.modnegexp = c("never", "when.fail", "always"),
- verbose=FALSE)
+ verbose=FALSE, return.info=FALSE)
}
\arguments{
@@ -48,7 +48,12 @@
constrained solution, even if the unconstrained one would have been
valid. See \sQuote{Details}. }
- \item{verbose}{ logical. Write out details? }
+ \item{verbose}{ a \code{logical} flag. Write out details? }
+
+ \item{return.info}{ a \code{logical} flag. If \code{TRUE}, details
+ about models and data will be added to the return value. See
+ \sQuote{Value}. }
+
}
\details{
This detrends and standardizes a tree-ring series. The detrending is
@@ -110,9 +115,69 @@
See the references below for further details on detrending.
}
\value{
+
If several methods are used, returns a \code{data.frame} containing
- the detrended series (\code{\var{y}}) according to the methods
- used. If only one method is selected, returns a vector.
+ the detrended series (\code{\var{y}}) according to the methods used.
+ The columns are named and ordered to match \code{\var{method}}. If
+ only one method is selected, returns a vector.
+
+ If \code{\var{return.info}} is \code{TRUE}, the return value is a
+ \code{list} with three parts:
+
+ \item{series}{ the main result described above (\code{data.frame} or
+ vector) }
+
+ \item{model.info}{ Information about the models corresponding to each
+ output series. Whereas \code{\var{series}} may return a vector,
+ \code{\var{model.info}} is always a list where each top level
+ element corresponds to one selected method. Also these elements are
+ named and arranged according to the methods selected. Each element
+ is a list with some of the following sub-elements, depending on
+ which detrending method was actually used:
+
+ \describe{
+
+ \item{n.zeros}{ See \code{"data.info"} below. Always present. }
+
+ \item{zero.years}{ See \code{"data.info"}. Always present. }
+
+ \item{method}{ The method actually used for detrending. One of
+ \code{"Mean"}, \code{"Line"}, \code{"ModNegExp"},
+ \code{"Spline"} or \code{"Ar"}. Always present. }
+
+ \item{mean}{ Mean of the input series, missing values removed.
+ Only for method \code{"Mean"}. }
+
+ \item{coefs}{ Coefficients of the model. Methods \code{"Line"}
+ and \code{"ModNegExp"}.}
+
+ \item{formula}{ The \code{"ModNegExp"} \code{\link{formula}}. }
+
+ \item{is.constrained}{ A \code{logical} flag indicating whether
+ the parameters of the \code{"ModNegExp"} model were
+ constrained. Only interesting when argument
+ \code{\var{constrain.modnegexp}} is set to \code{"when.fail"}. }
+
+ \item{nyrs}{ The value of \code{\var{nyrs}} used for
+ \code{\link{ffcsaps}}. Only for method \code{"Spline"}. }
+
+ \item{order}{ The order of the autoregressive model, selected by
+ AIC (Akaike information criterion). Only for method
+ \code{"Ar"}. }
+
+ \item{ar}{ The autoregressive coefficients used by method
+ \code{"Ar"}. A \code{numeric} vector ordered by increasing
+ lag. }
+
+ }
+
+ }
+
+ \item{data.info}{ Information about the input series: number
+ (\code{"n.zeros"}) and location (\code{"zero.years"}) of zero
+ values. If the locations are in a \code{character} vector, they are
+ years. Otherwise they are indices to the input series. }
+
}
\references{
Cook, E. R. and Kairiukstis, L. A. (1990) \emph{Methods of
More information about the Dplr-commits
mailing list