daily_response() can now detrend the aggregated climate series of each
moving window with the smoothers used by dplR::detrend(). Besides the
existing "SLD", dc_method accepts "Spline", "Mean", "Ar" and
"Friedman", controlled by the new arguments dc_nyrs, dc_f,
dc_difference, dc_span, dc_bass and dc_scale.
For strictly positive series the results are identical to dplR::detrend().
The ARSTAN safeguards that dplR applies when a fitted curve is not strictly
positive are not reproduced, because they would silently replace the selected
method for climate variables that cross zero; such windows return NA
instead.
dc_difference = FALSE (the dplR default) performs ratio detrending and is
rejected for variables with negative values, such as temperature in degrees
Celsius. Use dc_difference = TRUE for these.
The returned object gains a $detrending element recording the detrending
method and settings that were used.
Results for dc_method = "SLD" are unchanged.
dc_method = "SLD" no longer recycles fitted values when the climate series
contains missing values.
$optimized_return_all keeps its year row names when detrending is applied.
Warnings about missing values are no longer triggered by scattered gaps, such as DOY 366 in non-leap years, which the window aggregation removes. They now report only the years containing a gap at least as long as the shortest window width, and name those years.
critical_r() derived its threshold from the Pearson sampling distribution
regardless of which coefficient was actually computed, and the criteria that
rely on it — "alpha", "fdr_BH", "fdr_BY" and "sidak" — inherited that
assumption. Kendall's tau has a different null distribution, so those criteria
were far too strict on tau surfaces and discarded genuine climate responses.
critical_r() gains a cor_method argument, accepting "pearson" (the
default, unchanged), "spearman" (which uses the same transformation applied
to ranks) and "kendall". For Kendall's tau the critical value now follows the
normal approximation with variance 2(2n + 5) / 9n(n − 1), the same null
distribution stats::cor.test() uses when exact = FALSE. The four response
functions pass their cor_method or pcor_method through, and the false
discovery rate criteria build their p values from the matching distribution.
Results change for anyone using cor_method = "kendall" or
pcor_method = "kendall" together with remove_insignificant = TRUE. More
windows are now retained, correctly so. The magnitude depends on sample size:
the old threshold was 41 per cent too strict at n = 24 and 47 per cent too
strict at n = 73. On one test surface the "sidak" threshold fell from 0.443,
which retained no windows at all, to 0.276, against 0.297 from the permutation
test that had always been correct.
"maxstat" and "cluster" were never affected, because they build their null
by permuting the statistic actually being computed and so make no
distributional assumption. They remain the most reliable choice for rank
correlations, and for partial correlations with Kendall's tau in particular,
where the effective-sample-size adjustment used by the analytic criteria is
itself an approximation.
The logical argument remove_insignificant_boot has been removed and replaced
by remove_insignificant_method, which selects the criterion applied when
remove_insignificant = TRUE. This keeps one switch for whether to filter and
one for how, rather than adding an argument per method.
remove_insignificant_method accepts:
"alpha" — per-window significance threshold. This is the default and
reproduces the behaviour of earlier versions."boot" — the bootstrap confidence interval excludes zero. Requires
boot = TRUE."fdr_BY" — false discovery rate, Benjamini–Yekutieli. Valid under
arbitrary dependence between windows."fdr_BH" — false discovery rate, Benjamini–Hochberg. Valid under positive
dependence and less conservative than BY."maxstat" — max-statistic permutation test."sidak" — Šidák correction on the effective number of tests."cluster" — cluster-extent inference.The "alpha" and "boot" criteria are applied one window at a time and
therefore do not account for how many windows were screened. The remaining
criteria do.
"maxstat" builds a null distribution of the largest absolute correlation
obtained anywhere in the search by repeatedly shuffling the response, and
retains only windows exceeding the resulting threshold. The number of
permutations is set by the new argument perm_n (default 999). It makes no
distributional assumptions, but requires perm_n re-evaluations of the whole
search. It is available for method = "cor" and for method = "lm" with a
single predictor, but not for method = "brnn", where the permutation would
require refitting the network at every window.
"sidak" targets the same family-wise error rate as "maxstat" but estimates
it analytically. Because neighbouring windows overlap, the number of windows
screened greatly overstates the number of independent tests. The effective
number is obtained from the eigenvalues of the correlation matrix among the
screened windows, following Li and Ji (2005), and a Šidák correction is applied
to that number rather than to the nominal count. The correlation matrix has
rank at most the number of years, so its non-zero eigenvalues are computed from
a small years-by-years matrix; no permutations are needed and the criterion is
effectively instantaneous. It is available for every method, including
method = "brnn", where "maxstat" cannot run.
"sidak" is an approximation to "maxstat". Across 54 combinations of site and
climate variable the two thresholds agreed to within 0.016 on average, with a
maximum discrepancy of 0.038, for ordinary correlations; for partial
correlations the corresponding figures were 0.020 and 0.046. The approximation
is close to exact for strongly autocorrelated variables such as temperature and
tends to be slightly permissive for weakly autocorrelated variables such as
daily precipitation, because the eigenvalue approach does not represent
temporal autocorrelation. "maxstat" is therefore preferable for final results,
while "sidak" suits exploratory work and method = "brnn".
"cluster" is cluster-extent inference, adapted from neuroimaging. The surface
is thresholded at the new argument cluster_threshold, connected regions of
supra-threshold windows are identified, and each region's mass, the sum of
absolute correlations inside it, is compared against the largest mass arising
in permuted surfaces. It controls the family-wise error rate over clusters
rather than over windows, so it establishes that a region contains signal
without identifying which windows within it carry that signal. It is available
for method = "cor", because the permuted statistic has to be laid out on the
same surface as the observed one.
"cluster" has little discriminating power on tree-ring response surfaces, and
is provided mainly for completeness and comparison. A real climate signal forms
one large connected region whose mass greatly exceeds anything the permuted
surfaces produce, so on daily surfaces every cluster that forms is retained and
the permutation step never rejects; the number of windows kept is then set by
cluster_threshold rather than by the data. On a daily search of about 3,400
windows, raising the threshold from the nominal critical correlation to 0.45
reduced the retained area from 95 to 40 per cent while the max-statistic
threshold on the same data stayed near 0.36. On the much smaller monthly
surfaces the permutation occasionally does reject a cluster, but not stably
between runs. "maxstat" or "sidak" are preferable when the question is
which windows carry signal.
For the seascorr functions the degrees of freedom account for the control variable throughout, and the effective number of tests is estimated from the primary windows.
Applies to daily_response(), monthly_response(),
daily_response_seascorr() and monthly_response_seascorr().
The returned list element $boot_ci_filter has been renamed
$significance_filter. It records which method was applied, the criterion
used, and how many matrix cells were retained or removed. "alpha", "sidak",
"maxstat" and "cluster" also report the threshold applied, as
critical_r. For "maxstat" it additionally reports the observed maximum
(observed_max_r), the global p value (p_global), the number of
permutations, and the number of windows searched. For "sidak" it reports the
effective number of tests (effective_tests), the per-test significance level
(alpha_per_test), and the number of windows searched. For "cluster" it
reports the forming threshold used (forming_threshold), the critical cluster
mass (cluster_mass_critical), how many clusters formed and how many were
retained (n_clusters_formed, n_clusters_retained), the largest observed
mass, and the max-statistic threshold from the same permutations
(maxstat_critical_r), so the two criteria can be compared directly.
When "maxstat" removes every window, the message now reports the global p
value and states that no correlation exceeded what the search alone would
produce, rather than suggesting that the threshold be relaxed.
Missing values are no longer handled silently. All four functions now default
to cor_na_use = "everything" or pcor_na_use = "everything", so a window
overlapping a missing value returns NA and the user decides explicitly how to
proceed. Previously the seascorr functions defaulted to
pcor_na_use = "pairwise.complete.obs" and quietly computed over the available
years, which was inconsistent with daily_response() and monthly_response().
All four functions now warn when the climate data contain missing values and
the chosen rule does not handle them, naming the argument to set. Previously
only daily_response() did so.
For method = "cor" and for the seascorr partial correlations, a window whose
bootstrap produces any missing replicate is now NA in all four functions,
rather than being computed from the remaining replicates. For method = "lm"
and method = "brnn" the opposite rule applies, again in all four functions: a
missing replicate there indicates a model fit that did not converge rather than
missing data, so it is dropped and the window is computed from the replicates
that succeeded.
Windows are still aggregated with na.rm = TRUE, so a window with a small
number of missing days is aggregated from the days that are present.
The bootstrap for method = "cor", and for the partial correlations in the
seascorr functions, is now computed by matrix algebra rather than by
resampling loops. Each replicate correlation is obtained from weighted sums
over the original observations, so a whole window's replicates come from a
single matrix product.
Approximate speedups on a search of about 6500 windows with 1000 replicates:
| correlation | before | after | |---|---|---| | Pearson | ~3 min | ~4 s | | Spearman | ~9 min | ~45 s | | Kendall | ~19 min | ~55 s |
Confidence intervals of type "norm", "perc" and "basic" are now computed
directly from the replicates and reproduce boot::boot.ci() exactly. The
"bca" type, and the "lm" and "brnn" methods, are unchanged and continue
to use boot::boot().
Results will differ slightly from earlier versions even with the same seed. Resample indices are now drawn once and shared by every window, instead of being redrawn for each window. Differences are of the order of Monte Carlo error (typically below 0.02 in correlation units). The shared resamples mean all windows are evaluated on the same resampled years, so bootstrap surfaces are internally comparable rather than each window carrying independent noise.
cor_na_use and pcor_na_use are honoured by the new code path, including
"everything", under which any replicate that draws a missing value returns
NA, matching stats::cor().
boot_ci_type = "stud" is now rejected, with a message explaining that
studentized intervals require a variance estimate for each bootstrap
replicate, which is not available for correlation statistics. Previously
this combination was accepted, ran the full search, and then failed part
way through with the uninformative error number of items to replace is not a multiple of replacement length. boot_ci_type is validated once at the
start of the call, so an invalid value now fails immediately rather than
after minutes of computation.
When remove_insignificant = TRUE, the confidence-interval matrices
$boot_lower and $boot_upper are now masked with the same criterion as
$calculations. Previously they retained values for windows the filter had
removed, so plotting intervals over a filtered surface could suggest
significance for windows that had in fact been rejected.
Kendall, M. G. (1938). A new measure of rank correlation. Biometrika 30, 81–93.
Li, J. and Ji, L. (2005). Adjusting multiple testing in multilocus analyses using the eigenvalues of a correlation matrix. Heredity 95, 221–227.
Nichols, T. E. and Holmes, A. P. (2002). Nonparametric permutation tests for functional neuroimaging: a primer with examples. Human Brain Mapping 15, 1–25.
quantile_prob for daily_response() and monthly_response()quantile_prob_env_data_primary and quantile_prob_env_data_control for daily_response_seascorr() and monthly_response_seascorr()daily_response(), daily_response_seascorr(), monthly_response(), and monthly_response_seascorr().number_previous_years, which allows users to include more than one previous year in response-function analyses.previous_year now acts as the main switch for previous-year analyses. If previous_year = FALSE, previous-year settings are ignored and the analysis is performed for the current year only.day_interval or month_interval while previous_year = FALSE, the interval is ignored with a warning and the analysis is performed for the current year only.Y-3, Y-2, Y-1, and Y, with year separators where appropriate.brnn method. Occasional unstable BRNN model fits are now handled more safely, returning NA for failed windows instead of interrupting the full analysis.