[Dplr-commits] r1120 - in pkg/dplR: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 22 21:46:33 CEST 2018
Author: andybunn
Date: 2018-06-22 21:46:33 +0200 (Fri, 22 Jun 2018)
New Revision: 1120
Added:
pkg/dplR/R/sss.R
pkg/dplR/man/sss.Rd
Modified:
pkg/dplR/ChangeLog
pkg/dplR/NAMESPACE
pkg/dplR/man/rwi.stats.running.Rd
Log:
Adding sss as a new function.
Modified: pkg/dplR/ChangeLog
===================================================================
--- pkg/dplR/ChangeLog 2018-06-20 20:12:39 UTC (rev 1119)
+++ pkg/dplR/ChangeLog 2018-06-22 19:46:33 UTC (rev 1120)
@@ -1,5 +1,10 @@
* CHANGES IN dplR VERSION 1.6.9
+File: sss.R and .Rd
+----------------
+
+- Adding subsample signal strength as a stand alone function
+
File: spag.plot.R and .Rd
----------------
@@ -29,7 +34,7 @@
File: rwl.stats.R and .Rd
----------------
-- Removing sens1 and sens2 from this function. At world dendro I'm still seeing these being used by dplR users. It's a terrible stat. I'm taking it out. Help file ammended.
+- Removing sens1 and sens2 from this function. At world dendro I'm still seeing these being used by dplR users. It's a terrible stat. I'm taking it out. Help file amended.
File: detrend.series.R and .Rd
----------------
Modified: pkg/dplR/NAMESPACE
===================================================================
--- pkg/dplR/NAMESPACE 2018-06-20 20:12:39 UTC (rev 1119)
+++ pkg/dplR/NAMESPACE 2018-06-22 19:46:33 UTC (rev 1120)
@@ -55,7 +55,7 @@
write.tucson, plot.rwl, interseries.cor, summary.rwl, plot.crn,
insert.ring, delete.ring, xskel.ccf.plot, xskel.plot, latexify,
latexDate, rasterPlot, treeMean, rwl.report, print.rwl.report,
- plotRings,time.rwl,time.crn,csv2rwl,pass.filt,as.rwl)
+ plotRings,time.rwl,time.crn,csv2rwl,pass.filt,as.rwl,sss)
S3method(print, redfit)
S3method(plot, rwl)
Added: pkg/dplR/R/sss.R
===================================================================
--- pkg/dplR/R/sss.R (rev 0)
+++ pkg/dplR/R/sss.R 2018-06-22 19:46:33 UTC (rev 1120)
@@ -0,0 +1,37 @@
+sss <- function(rwi,ids=NULL){
+ # rwi.stats is robust enough to have a single call regardless of
+ # whether ids is passed in because if no ids are passed in
+ # rbar.eff == rbar.bt and n.cores=n.trees
+ # Thus, if ids is NULL assume one core / tree use and
+ # we use rbar.bt (which is the same as rbar.eff). If there are ids
+ # we use rbar.eff which is not the same as rbar.bt. Ditto for getting
+ # N. Satisfy yourself via:
+ #ca533.rwi <- detrend(ca533,method="Spline")
+ #ca533.ids <- autoread.ids(ca533)
+ #rwi.stats(ca533.rwi)
+ #rwi.stats(ca533.rwi,ca533.ids)
+
+ rwiVars <- rwi.stats(rwi, ids=ids)
+ rbar <- rwiVars$rbar.eff
+ N <- rwiVars$n.trees
+
+
+ if(is.null(ids)){
+ # n is is the number of cores in subsample
+ n <- rowSums(!is.na(rwi))
+ }
+ else {
+ # n is the number of trees in the subsample.
+ # calculating n is kind of tedious:
+ # we need n trees, not n cores in a year
+ colnames.rwi <- colnames(rwi)
+ n <- rep(NA,nrow(rwi))
+ for(i in 1:nrow(rwi)){
+ cols.with.data <- c(!is.na(rwi[i,]))
+ trees.this.year <- ids$tree[rownames(ids) %in% colnames.rwi[cols.with.data]]
+ n[i] <- length(unique(trees.this.year))
+ }
+ }
+ res <- (n*(1+(N-1)*rbar)) / (N*(1+(n-1)*rbar))
+ res
+}
Property changes on: pkg/dplR/R/sss.R
___________________________________________________________________
Added: svn:eol-style
+ native
Modified: pkg/dplR/man/rwi.stats.running.Rd
===================================================================
--- pkg/dplR/man/rwi.stats.running.Rd 2018-06-20 20:12:39 UTC (rev 1119)
+++ pkg/dplR/man/rwi.stats.running.Rd 2018-06-22 19:46:33 UTC (rev 1120)
@@ -132,6 +132,8 @@
Fritts (2001) for further details for computational details on the
output. The signal-to-noise ratio is calculated following Cook and
Pederson (2011).
+
+ Note that Buras (2017) cautions against using the expressed population signal as a statistic to determine the whether a chronology correctly represents the population signal of a data set. He reccomends the use of subsample signal strength (\code{\link{sss}}) over EPS.
If desired, the \code{\var{rwi}} can be filtered in the same manner
as the family of cross-dating functions using \code{\var{prewhiten}} and
@@ -198,6 +200,8 @@
}
\references{
+ Buras, A. (2017) A comment on the Expressed Population Signal. Dendrochronologia 44:130-132.
+
Cook, E. R. and Kairiukstis, L. A., editors (1990) \emph{Methods of
Dendrochronology: Applications in the Environmental Sciences}.
Springer. \acronym{ISBN-13}: 978-0-7923-0586-6.
@@ -236,51 +240,5 @@
rwi.stats(gp.rwi, gp.ids, period="common")
rwi.stats.legacy(gp.rwi, gp.ids) # rwi.stats prior to dplR 1.5.3
-\dontrun{
- library(graphics)
- def.par <- par(no.readonly=TRUE)
- ## Plot the chronology showing a potential cutoff year based on EPS
- eps.cut <- 0.92 # An arbitrary EPS cutoff for demonstration
- gp.crn <- chron(gp.rwi)
- ## Running stats on the rwi with an window
- foo <- rwi.stats.running(gp.rwi, gp.ids, window.length = 80)
- yrs <- time(gp.crn)
- bar <- data.frame(yrs = c(min(yrs), foo$mid.year, max(yrs)),
- eps = c(NA, foo$eps, NA))
- par(mar = c(2, 2, 2, 2), mgp = c(1.1, 0.1, 0), tcl = 0.25,
- mfcol = c(2, 1), xaxs='i')
- plot(yrs, gp.crn[, 1], type = "n", xlab = "Year",
- ylab = "RWI", axes=FALSE)
- cutoff <- max(bar$yrs[bar$eps < eps.cut], na.rm = TRUE)
- xx <- c(500, 500, cutoff, cutoff)
- yy <- c(-1, 3, 3, -1)
- polygon(xx, yy, col = "grey80")
- abline(h = 1, lwd = 1.5)
- lines(yrs, gp.crn[, 1], col = "grey50")
- lines(yrs, ffcsaps(gp.crn[, 1], nyrs = 32), col = "red", lwd = 2)
- axis(1); axis(2); axis(3);
- par(new = TRUE)
- ## Add EPS
- plot(bar$yrs, bar$eps, type = "b", xlab = "", ylab = "",
- axes = FALSE, pch = 20, col = "blue")
- axis(4, at = pretty(foo$eps))
- mtext("EPS", side = 4, line = 1.1)
- box()
- ## Second plot is the chronology after the cutoff only
- ## Chronology is rebuilt using just years after cutoff but
- ## that difference is essentially nil.
- yr.mask <- yrs > cutoff
- yrs2 <- yrs[yr.mask]
- gp.crn2 <- chron(gp.rwi[yr.mask, ])
- plot(yrs2, gp.crn2[, 1], type = "n",
- xlab = "Year", ylab = "RWI", axes=FALSE)
- abline(h = 1, lwd = 1.5)
- lines(yrs2, gp.crn2[, 1], col = "grey50")
- lines(yrs2, ffcsaps(gp.crn2[, 1], nyrs = 32),
- col = "red", lwd = 2)
- axis(1); axis(2); axis(3); axis(4)
- box()
- par(def.par)
}
-}
\keyword{ misc }
Added: pkg/dplR/man/sss.Rd
===================================================================
--- pkg/dplR/man/sss.Rd (rev 0)
+++ pkg/dplR/man/sss.Rd 2018-06-22 19:46:33 UTC (rev 1120)
@@ -0,0 +1,96 @@
+\encoding{UTF-8}
+\name{sss}
+\alias{sss}
+\title{ Subsample Signal Strength }
+\description{
+ Calculate subsample signal strength on a
+ \code{data.frame} of (usually) ring-width indices.
+}
+\usage{
+sss(rwi, ids = NULL)
+}
+\arguments{
+
+ \item{rwi}{ a \code{data.frame} with detrended and standardized ring
+ width indices as columns and years as rows such as that produced by
+ \code{\link{detrend}}. }
+
+ \item{ids}{ an optional \code{data.frame} with column one named
+ \code{"tree"} giving a \code{numeric} \acronym{ID} for each tree and
+ column two named \code{"core"} giving a \code{numeric} \acronym{ID}
+ for each core. Defaults to one core per tree as\cr
+ \code{data.frame(tree=1:ncol(\var{rwi}), core=rep(1, ncol(\var{rwi})))}. }
+}
+\details{
+ This calculates subsample signal strength (sss) following equation 3.50 in Cook and Kairiukstis (1990) but using notation from Buras (2017) because writing the prime unicode symbol seems too difficult. The function calls \code{\link{rwi.stats}} and passes it the arguments \code{ids} and \code{prewhiten}.
+
+ To make better use of variation in growth within and between series, an appropriate mask (parameter \code{\var{ids}}) should be provided that identifies each series with a tree as it is common for dendrochronologists to take more than one core per tree. The function \code{\link{read.ids}} is helpful for creating a mask based on the series \acronym{ID}.
+
+ Subsample signal strength is calculated as \eqn{\frac{n[1+(N-1)\bar{r}]}{N[1+(n-1)\bar{r}]}}{n*(1+(N-1)*rbar) / N*(1+(n-1)*rbar)} where \code{n} and \code{N} are the number of cores or trees in the subsample and sample respectively and \code{rbar} is mean interseries correlation. If there is only one core per tree \code{n} is the sample depth in a given year (\code{rowSums(!is.na(rwi))}), \code{N} is the number of cores (\code{n.cores} as given by \code{\link{rwi.stats}}), and \code{rbar} is the mean interseries correlation between all series (\code{r.bt} as given by \code{\link{rwi.stats}}). If there are multiple cores per tree \code{n} is the number of trees present in a given year, \code{N} is the number of trees (\code{n.trees} as given by \code{\link{rwi.stats}}), and \code{rbar} is the effective mean interseries correlation (\code{r.eff} as given by \code{\link{rwi.stats}}).
+
+Readers interested in the differences between subsample signal strength and the more commonly used (running) expressed population signal should look at Buras (2017) on the common misuse of the expressed population signal as well as Cook and Pederson (2011) for a more general appraoch to categorizing variability in tree-ring data.
+
+}
+
+\value{ A \code{numeric} containing the subsample signal strength that is the same as number if rows of\code{rwi}.
+}
+\references{
+
+ Buras, A. (2017) A comment on the Expressed Population Signal. Dendrochronologia 44:130-132.
+
+ Cook, E. R. and Kairiukstis, L. A., editors (1990) \emph{Methods of
+ Dendrochronology: Applications in the Environmental Sciences}.
+ Springer. \acronym{ISBN-13}: 978-0-7923-0586-6.
+
+ Cook, E. R. and Pederson, N. (2011) Uncertainty, Emergence, and
+ Statistics in Dendrochronology. In Hughes, M. K., Swetnam, T. W., and
+ Diaz, H. F., editors, \emph{Dendroclimatology: Progress and
+ Prospects}, pages 77\enc{–}{--}112. Springer. \acronym{ISBN-13}:
+ 978-1-4020-4010-8.
+
+}
+
+\author{Andy Bunn. Patched and improved by Mikko Korpela. }
+\seealso{ \code{\link{rwi.stats}}, \code{\link{read.ids}} }
+\examples{
+data(ca533)
+ca533.rwi <- detrend(ca533,method="Spline")
+# assuming 1 core / tree
+ca533.sss <- sss(ca533.rwi)
+
+ca533.ids <- autoread.ids(ca533)
+# done properly with >=1 core / tree as per the ids
+ca533.sss2 <- sss(ca533.rwi,ca533.ids)
+
+yr <- time(ca533)
+plot(yr,ca533.sss,type="l",ylim=c(0.4,1),
+ col="darkblue",lwd=2,xlab="Year",ylab="SSS")
+lines(yr,ca533.sss2,lty="dashed",
+ col="darkgreen",lwd=2)
+
+# Plot the chronology showing a potential cutoff year based on SSS
+ca533.crn <- chron(ca533.rwi)
+def.par <- par(no.readonly=TRUE)
+par(mar = c(2, 2, 2, 2), mgp = c(1.1, 0.1, 0), tcl = 0.25, xaxs='i')
+plot(yr, ca533.crn[, 1], type = "n", xlab = "Year",
+ ylab = "RWI", axes=FALSE)
+cutoff <- max(yr[ca533.sss2 < 0.85])
+xx <- c(500, 500, cutoff, cutoff)
+yy <- c(-1, 3, 3, -1)
+polygon(xx, yy, col = "grey80")
+abline(h = 1, lwd = 1.5)
+lines(yr, ca533.crn[, 1], col = "grey50")
+lines(yr, ffcsaps(ca533.crn[, 1], nyrs = 32), col = "red", lwd = 2)
+axis(1); axis(2); axis(3);
+par(new = TRUE)
+## Add EPS
+plot(yr, ca533.sss2, type = "l", xlab = "", ylab = "",
+ axes = FALSE, col = "blue")
+abline(h=0.85,col="blue",lty="dashed")
+axis(4, at = pretty(ca533.sss2))
+mtext("SSS", side = 4, line = 1.1, lwd=1.5)
+box()
+par(def.par)
+
+}
+\keyword{ misc }
Property changes on: pkg/dplR/man/sss.Rd
___________________________________________________________________
Added: svn:eol-style
+ native
More information about the Dplr-commits
mailing list