Fix confidence intervals and improve DOF calculation in FilterSpec

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: PaleoSpec
 Title: Spectral tools for the ECUS group
-Version: 0.32
+Version: 0.33
 Authors@R: c(
   person("Thomas", "Laepple", email = "tlaepple@awi.de", role = c("aut", "cre")),
   person("Thomas", "Muench", email = "tmuench@awi.de", role = c("aut")),

NEWS.md

@@ -1,3 +1,14 @@
+# PaleoSpec 0.33
+
+* AddConfInterval calculation changed calculated confidence intervals give the 
+  interval within which the true spectrum should lie with frequency 1-p. 
+  Previous versions returned the interval which, when applied to the 'true' 
+  spectrum, would contain new spectral estimates with frequency 1-p.
+
+* Improvements to calculation of DOF in FilterSpec and FilterSpecLog. DOF 
+  calculation now accounts for non independence of ordinates from existing 
+  filter or tapering.
+
 # PaleoSpec 0.32
 
 * SpecACF can now use slepian tapers

R/AddConfInterval.R

@@ -14,33 +14,62 @@
 #' @return the input object including the new list elements \code{lim.1} and
 #'   \code{lim.2} giving the upper and lower bound of the confidence interval,
 #'   respectively.
-#' @author Thomas Laepple
+#' @description
+#' Calculated confidence intervals give the interval within which the true
+#' spectrum should lie with frequency 1-p. This behaviour changed with version
+#' 0.33.
+#'
+#' Previous versions returned the interval which, when applied to the 'true' spectrum, would contain new spectral estimates
+#' with frequency 1-p.
+#'
+#' @author Thomas Laepple, Andrew Dolman
 #' @examples
 #'
-#' N.R <- 1000
-#' N.T <- 100
-#' save.spec <- matrix(NA, N.T / 2, N.R)
-#' for (i.R in 1 : N.R) {
-#'   save.spec[, i.R] <- SpecMTM(ts(SimPowerlaw(1, N.T)))$spec
-#' }
+#' alpha <- 0.1
+#' beta <- 1
+#'
+#' spec_sim <- SpecMTM(ts(SimPLS(N = 1e03, beta = beta, alpha = alpha)))
+#'
+#' # Using a large nominal p value of 0.25 to reduce variation between random
+#' # timeseries
+#'
+#' # The true spec should be p/2 times above the lower CI, and p/2 below the upper CI
+#'
+#' nominal_p <- 0.25
+#' spec_sim <- AddConfInterval(spec_sim, pval = nominal_p)
+#' true_spec <- alpha * spec_sim$freq^-beta
+#'
+#' LPlot(spec_sim)
+#' abline(a = log10(alpha), b = -beta, lty = 2, col = "red")
 #'
-#' q.empirical <- apply(save.spec, 1, quantile, c(0.025, 0.975))
-#' testspec <- SpecMTM(ts(SimPowerlaw(1, N.T)))
+#' SpecCoverage <- function(spec, true_spec){
 #'
-#' LPlot(AddConfInterval(testspec), ylim = c(0.05, 10))
-#' lines(testspec$freq, q.empirical[1, ], col = "red")
-#' lines(testspec$freq, q.empirical[2, ],col = "red")
-#' legend("bottomleft", lwd = 2, col = c("black", "red"),
-#'        legend = c("one realization with chisq conf intervals",
-#'                   "MC confidence intervals"))
+#'   stopifnot("lim.1" %in% names(spec))
+#'
+#'   n <- length(spec$freq)
+#'   below <- sum(spec$lim.1 < true_spec)
+#'   above <- sum(spec$lim.2 > true_spec)
+#'
+#'   total <- below + above
+#'
+#'   count <- list(n = n, above = above, below = below, total = total,
+#'                 p_above = above / n, p_below = below / n,
+#'                 p_total = total / n,
+#'                 nominal_pval = spec$pval)
+#'
+#'   as.data.frame(count)
+#' }
+#'
+#' SpecCoverage(spec_sim, true_spec)
 #'
 #' @export
 AddConfInterval <- function(spec, pval = 0.05, MINVALUE = 1e-10) {
 
   is.spectrum(spec)
 
-  spec$lim.1 <- spec$spec * qchisq(c(1 - pval / 2), spec$dof) / (spec$dof)
-  spec$lim.2 <- spec$spec * qchisq(c(pval / 2), spec$dof) / (spec$dof)
+  spec$lim.1 <- spec$spec * 1 / (qchisq(c(pval / 2), spec$dof) / (spec$dof))
+  spec$lim.2 <- spec$spec * 1 / (qchisq(c(1 - pval / 2), spec$dof) / (spec$dof))
+
   spec$lim.1[spec$lim.1 < MINVALUE] <- MINVALUE
   spec$lim.2[spec$lim.2 < MINVALUE] <- MINVALUE
 

R/FilterSpec.R

@@ -1,11 +1,14 @@
 #' Filter a Power Spectrum Object
 #'
 #' @param spec A spec object
+#' @param method single integer for choosing an endpoint constraint method.
+#' Available choices are integers 0-4, see details of \code{\link{ApplyFilter}}
 #' @param keep_low_f Keep filtered (smoothed) low frequencies or replace with unfiltered
 #' @inheritParams stats::spec.pgram
 #' @inheritParams ApplyFilter
 #' @return A spec object (list)
 #' @family functions to filter / smooth spectra
+#' @author Andrew Dolman
 #' @export
 #'
 #' @examples
@@ -57,34 +60,47 @@ FilterSpec <- function(spec, spans, method = 3, keep_low_f = TRUE) {
   kernel <- stats::kernel("modified.daniell", spans %/% 2)
   filter <- kernel[-kernel$m:kernel$m]
 
-  spec_filt <- ApplyFilter(spec$spec, filter = filter, method = method)
+  spec_filt <- ApplyFilter(spec$spec, filter = filter,
+                           method = method)
 
   if (keep_low_f == FALSE) {
     # replace filtered spec with original in area where freqs have been reflected
     i <- 1:ceiling(length(filter) / 2)
     spec_filt[i] <- spec$spec[i]
-
     iend <- length(spec$freq) - (i-1)
-
     spec_filt[iend] <- spec$spec[iend]
-
   }
 
   spec$spec <- as.numeric(spec_filt)
 
-  # degrees of freedom of the kernel
-  df.kern <- stats::df.kernel(kernel)
 
-  spec$dof <- df.kern * spec$dof / 2
+  ## New degrees of freedom of the filtered spec
+  ## This can be estimated using the convolution of the new filter and an
+  ## approximate equivalent filter to the existing DOF
 
-  if (keep_low_f == FALSE) {
+  # Approximate an equivalent filter to the existing DOF
+  fl <- round(mean(spec$dof)/2)
+  if ((fl%%2) == 0) fl <- fl-1 #ensure that the filter length is odd
+
+  prior_filter <- rep(1/fl, fl)
+
+  combined_filter <- convolve(filter, prior_filter, type = "open")
+
+  m3 <- floor(length(combined_filter)/2)
+  kernel3 <- combined_filter[1:ceiling(length(combined_filter)/2)]
+  kernel3 <- kernel(coef = sort(kernel3, decreasing = TRUE), m = m3)
+
+  new.dof <- stats::df.kernel(kernel3)
+
+  spec$dof <- rep(new.dof, length(spec$freq))
 
+
+  if (keep_low_f == FALSE) {
     i <- 1:ceiling(length(filter) / 2)
     spec$dof[i] <- dof0[i]
 
     iend <- length(spec$freq) - (i-1)
     spec$dof[iend] <- dof0[iend]
-
   }
 
   # Adjust DOF in reflected filter region
@@ -116,7 +132,6 @@ FilterSpec <- function(spec, spans, method = 3, keep_low_f = TRUE) {
 
   spec <- AddConfInterval(spec)
 
-
   return(spec)
 }
 
@@ -126,11 +141,12 @@ FilterSpec <- function(spec, spans, method = 3, keep_low_f = TRUE) {
 #'
 #' @param spec A spec object
 #' @inheritParams LogSmooth
-#' @inheritParams ApplyFilter
-#'
+#' @param method single integer for choosing an endpoint constraint method.
+#' Available choices are integers 0-4, see details of \code{\link{ApplyFilter}}
 #' @return A spec object (list)
 #' @family functions to filter / smooth spectra
 #' @export
+#' @author Andrew Dolman
 #' @examples
 #' library(PaleoSpec)
 #'
@@ -219,7 +235,7 @@ FilterSpecLog <- function(spec,
   dof0 <- spec$dof
 
   # Gets the difference in delta_f for the log and standard freq axis
-  NpF <- function(freq, fw, df){
+  NpF2 <- function(freq, fw, df){
 
     posdiff <- (exp(log(freq) + df) - freq)
     negdiff <- (freq - exp(log(freq) - df))
@@ -228,9 +244,12 @@ FilterSpecLog <- function(spec,
 
     2 * fw * (fdiff/df) * 1/(2*max(freq))
   }
-  df.logkern <- NpF(spec$freq, length(filter), df = diff(freq_logspace[1:2]))
 
-  spec$dof <- spec$dof + df.mod * df.logkern * spec$dof/2
+  deltaf.logkern <- NpF2(spec$freq, length(filter), df = diff(freq_logspace[1:2]))
+
+  # this would be more accurate as the dof of the convolution of the filter not sum of the old and new dofs
+  spec$dof <- df.mod * deltaf.logkern + spec$dof
+
   spec$shape <- spec$dof/2
   spec$spans <- paste(spans, collapse = ",")
 
@@ -244,7 +263,5 @@ FilterSpecLog <- function(spec,
 
   spec <- AddConfInterval(spec)
 
-
-
   return(spec)
 }

tests/testthat/test-AddConfInterval.R

@@ -10,8 +10,10 @@ test_that("adding the confidence interval is correct.", {
 
   pval <- 0.05
 
-  lim1 <- spec * qchisq(1 - pval / 2, dof) / dof
-  lim2 <- spec * qchisq(pval / 2, dof) / dof
+  lim1 <- spec * 1 / (qchisq(c(pval / 2), dof) / (dof))
+  lim2 <- spec * 1 / (qchisq(c(1 - pval / 2), dof) / (dof))
+
+
   lim1[lim1 < MINVALUE] <- MINVALUE
   lim2[lim2 < MINVALUE] <- MINVALUE
 
@@ -23,8 +25,9 @@ test_that("adding the confidence interval is correct.", {
   MINVALUE <- 1.
   pval <- 0.1
 
-  lim1 <- spec * qchisq(1 - pval / 2, dof) / dof
-  lim2 <- spec * qchisq(pval / 2, dof) / dof
+  lim1 <- spec * 1 / (qchisq(c(pval / 2), dof) / (dof))
+  lim2 <- spec * 1 / (qchisq(c(1 - pval / 2), dof) / (dof))
+
   lim1[lim1 < MINVALUE] <- MINVALUE
   lim2[lim2 < MINVALUE] <- MINVALUE
 

tests/testthat/test-gg_spec.R

@@ -11,7 +11,7 @@ sp1 <- SpecMTM(ts1)
 sp2 <- SpecMTM(ts2)
 
 # 1 as spec_df
-sp1_df <- as.data.frame(sp1)
+sp1_df <- Spec2DF(sp1)
 
 # list of spectra, 1 is a spec_df
 sp_lst <- list(sp2 = sp2, sp1_df = sp1_df)
@@ -64,3 +64,4 @@ test_that("removeFirst and removeLast operate on a per spec_id basis", {
   expect_that(nrow(g0$data) - nrow(gl3$data), equals(6))
 
 })
+