Introduction

This document presents mathematical details about the Dendrochronology Program Library in R (dplR) (Bunn 2008, 2010) which is an add-on package for R (R Core Team 2015). The first half deals with the spline smoothing function caps; the second covers the computation of Gini coefficients in gini.coef.

The original implementations of the functions covered here were not written by the author of this document. Therefore the functions were analyzed with a reverse engineering approach.

This document was first written when spline smoothing in dplR was performed by ffcsaps, a pure R function. As of dplR version 1.7.3 that role belongs to caps, a wrapper around a Fortran subroutine from Ed Cook’s ARSTAN, and ffcsaps is deprecated. The two functions parameterize the spline differently but, as the section on equivalence demonstrates, they compute the same spline. The analysis below has been rewritten around caps, with the ffcsaps parameterization retained at the end because it explains a factor of two that a reader comparing the two implementations will otherwise trip over.

Spline smoothing parameters in caps

The caps function fits a cubic smoothing spline to a given data vector. In the manual (Rd file) of the function (Bunn et al. 2025), it is stated that the frequency response of the spline is f at a wavelength (period) of nyrs years¹, where these two are parameters of the function. We aim to clarify how they relate to the single smoothing parameter of the spline and what that parameter stands for.

The manual of the caps function cites Cook and Kairiukstis (1990). On page 111, they give the following frequency (amplitude) response function for the spline:

$$u(f)=1-\frac{1}{1 + \frac{p(\cos (2\pi f) +2)}{6(\cos (2\pi f) -1)^2}} \qquad (1)$$

where $f$ is frequency and $p$ is stated to be the Lagrange multiplier of the spline, the single parameter that determines the frequency response. However, the exact definition of the optimization problem is absent. Neither is it given in Cook and Peters (1981), the reference used by Cook and Kairiukstis (1990). I did not find a copy of Peters and Cook (1981) when trying to follow the chain of references further.

Note that the relationship between frequency and period using mixed notation of caps and equation (1) is $f = 1/\mathtt{nyrs}$. Setting parameters f and nyrs in caps is equivalent to the following directive: set the smoothing parameter to a value that fulfills $u(1/\mathtt{nyrs}) = \mathtt{f}$. By making the variable substitutions and rearranging equation (1) we get the following equation for $p$:

$$p = \frac{6 \mathtt{f} (\cos (2\pi / \mathtt{nyrs}) -1)^2}{(1 - \mathtt{f})(\cos (2\pi / \mathtt{nyrs}) +2)} \qquad (2)$$

p=((1.d0/(1.d0-pct)-1.d0)*6.d0*(cos(pi*2.d0/v)-1.d0)**2)/(cos(pi*2.d0/v)+2.d0)

Empirical frequency response

Whether the fitted spline actually has the advertised frequency response is a separate question from whether the code implements equation (2) correctly, and it is worth checking. We smooth 500 independent series of 1536 i.i.d. standard normal samples, take the ratio of the modulus of the discrete Fourier transform of the smoothed series to that of the input, and average over the repeats.

Figure 1. Theoretical frequency response of the spline filter (equations 1 and 2, green line) versus the response measured with i.i.d. normal series of 1536 samples, mean of 500 repeats, using caps (blue circles). The legend on the bottom panel applies to all panels.

Figure 1 shows the result. Theory meets practice well, particularly for low frequencies. The measured response at the nominal cutoff frequency $1/\mathtt{nyrs}$ is 0.504, 0.511 and 0.530 for $\mathtt{nyrs} = 4$, 16 and 64 respectively, against a nominal $\mathtt{f} = 0.5$. It must be noted that the theoretical result does not take into account the effect of having a series of finite length, which is why the agreement degrades as nyrs grows toward the length of the series.

Equivalence of caps and the legacy ffcsaps

Users with results produced by dplR 1.7.2 or earlier will want to know whether caps changed any numbers. It did not, with one documented exception. Below, the pure R implementation of ffcsaps as it stood in dplR 1.6.9 is reproduced as ffcsaps.legacy and compared against caps.

Table 1 gives the largest absolute difference between the two implementations across four test series, two values of nyrs and two values of f. The worst case is 2.5e-12, which is floating-point noise. For integer nyrs, caps and ffcsaps compute the same spline.

There is one genuine difference. caps passes nyrs to Fortran as an integer, so a fractional nyrs is truncated, whereas ffcsaps used it as given. Fitting series CAM011 of the ca533 data set with $\mathtt{nyrs} = 302.67$, two thirds of the series length, the two differ by 1.23e-04, which is 0.01% of the range of the series. Passing the truncated value $302$ to ffcsaps instead brings the difference back down to 1.1e-11, confirming that truncation is the whole of the discrepancy.

This is reachable in ordinary use, by two routes. A nyrs between 0 and 1 selects the proportion-of-series-length shorthand, and caps multiplies it by the series length, which will rarely land on a whole number. Internally, plot.crn, wavelet.plot and ssf pass a fractional nyrs of their own, computed as a fixed proportion of the series length; detrend.series and rcs apply floor first and so are unaffected. The effect on the fitted curve is small, but it is not zero.

Figure 2. Top: series CAM011 of the ca533 data set (grey) with a 32-year spline fitted by caps (blue) and by the legacy ffcsaps (vermillion, dashed); the two curves are indistinguishable. Bottom: the difference between the two fitted curves, in ring-width units.

Figure 2 shows the two fits on a real ring-width series together with their difference, which is at the level of the floating-point representation.

A note on the ffcsaps parameterization

The deprecated ffcsaps contained code lines corresponding to the equation

$$\mathtt{p.inv} = \frac{1}{\mathtt{p}} = \frac{(1 - \mathtt{f})(\cos (2\pi / \mathtt{nyrs}) +2)}{12 \mathtt{f} (\cos (2\pi / \mathtt{nyrs}) -1)^2} + 1 \qquad (4)$$

where $\mathtt{p}$ and its inverse $\mathtt{p.inv}$ are variables used in the code. Writing equation (2) for the inverse,

$$\frac{1}{p} = \frac{(1 - \mathtt{f})(\cos (2\pi / \mathtt{nyrs}) +2)}{6 \mathtt{f} (\cos (2\pi / \mathtt{nyrs}) -1)^2} \qquad (5)$$

we find that equations (5) and (4) are connected by

$$\frac{1}{p} = 2 \left(\frac{1}{\mathtt{p}} - 1\right) \qquad (6)$$

or equivalently

$$\frac{\mathtt{p}}{1 - \mathtt{p}} = 2 p \qquad (7)$$

So the variable named p in ffcsaps and the Lagrange multiplier $p$ of Cook and Kairiukstis (1990) were not the same quantity, despite sharing a name. They are two parameterizations of the same penalty, related by equation (7). This is the factor of two that a reader comparing src/capsf.f95 against ffcsaps would otherwise have to discover the hard way.

The ffcsaps form is the convex-combination parameterization, in which the spline minimizes

$$\mathtt{p} \times \text{Error} + (1 - \mathtt{p}) \times \text{Roughness} \qquad (8)$$

with $\mathtt{p} \in [0, 1]$. Following from equations (7) and (8), the splines described in Cook and Kairiukstis (1990), and hence those computed by caps, are the result of minimizing

$$2 p \times \text{Error} + \text{Roughness} \qquad (9)$$

with the same definitions of Error and Roughness, details of which are omitted here. The section above confirms empirically that the two forms give the same fitted curve.

Formulation of the Gini coefficient in gini.coef

The gini.coef function computes the Gini coefficient (Gini index) of a given data vector. The manual (Rd file) of the function has a reference to Biondi and Qeadan (2008) which uses the following formula for the Gini coefficient ($G$):

$$G = \frac{1}{2 n \sum_{i=1}^{n} x_i} \sum_{i=1}^{n} \sum_{j=1}^{n} \left| x_i - x_j \right| \qquad (10)$$

In equation (10), the Gini coefficient is defined in terms of pairwise differences between all pairs of observations ($x_i,\ i \in 1, \dots, n$). More specifically, the Gini coefficient is one half of the relative mean difference, which is defined as the mean of the absolute pairwise distances divided by the mean of the observations.

The C source code of the gini.coef function uses the following formula for the Gini index:

$$G = \left(X_n (n - 1) - 2 \sum_{i=1}^{n-1}X_i\right) / (X_n n) \qquad (11)$$

where $n$ is the number of observations and $X_i$ is the $i$th cumulative sum

$$X_i = \sum_{j=1}^{i} x_j \qquad (12)$$

of sorted observations $x_j$:

$$\forall i: i < j \Rightarrow x_i \leq x_j \qquad (13)$$

Equation (11) can be reformulated as

$$G = 1 - \frac{1}{n} - \frac{2}{X_n n} \sum_{i=1}^{n-1}X_i \qquad (14)$$

or as

$$G = \left(\frac{1}{2} - \left(\frac{1}{2n} + \frac{1}{X_n n} \sum_{i=1}^{n-1}X_i\right)\right) / \frac{1}{2} \qquad (15)$$

When we assign

$$A + B = \frac{1}{2} \qquad (16)$$

and

$$B = B_1 + B_2 = \frac{1}{2n} + \frac{1}{X_n n} \sum_{i=1}^{n-1}X_i \qquad (17)$$

equation (15) becomes

$$G = A / (A + B) \qquad (18)$$

or equivalently

$$G = 1 - 2 B \qquad (19)$$

Figure 3 is a graphical representation of the Gini coefficient using an example data set of the following six observed values: $\{0.2, 0.4, 0.75, 0.95, 1.2, 2.5\}$. It shows the definition of the Gini coefficient as the ratio of the area above the Lorenz curve (Lorenz 1905) to the total area of the triangle (Xu 2003). The Lorenz curve is defined by the cumulative distribution function of the empirical probability distribution of the observations. The sides of the triangle corresponding to the axes are normalized to length 1.

Figure 3. Graphical representation of the Gini coefficient based on areas defined by the Lorenz curve (n = 6). See equations (16), (17), (18) and (19).

Comparing Figure 3 to equation (17), $B_2 = \sum_{i=1}^{n-1}X_i / (X_n n)$ is the sum of the areas of the cyan bars. Summing the areas of the teal triangles, we get

$$\sum_{i=1}^{n}\left( \frac{1}{2} \frac{1}{n} \frac{x_i}{X_n} \right) = \frac{1}{2 n X_n}\sum_{i=1}^{n} x_i = \frac{1}{2 n} = B_1 \qquad (20)$$

Note that $B_1$ only depends on the number of observations, not on their values. From equations (17) and (19) we find that the value of the Gini coefficient at maximum inequality (winner takes all) is $G_{\text{max}}(n)=1 - 1 / n$. When all observed values are equal, the Lorenz curve matches the line of equality, and the Gini coefficient is $G_{\text{min}}=0$. We have assumed that all values $x_i$ are non-negative.

The equivalence of different definitions of the Gini coefficient is reviewed in Xu (2003). One of the results shown in the paper is that the geometric definition (18) used by the gini.coef function is equivalent to the definition based on the relative mean difference (10). This can be experimentally verified by comparing the results of the following R function to those of gini.coef.

## Gini index is one half of relative mean difference.
## x should not have NA values.
gini.rmd <- function(x) {
    mean(abs(outer(x, x, "-"))) / mean(x) * 0.5
}

giniMax <- max(abs(vapply(ca533, function(x) {
    x <- x[!is.na(x)]
    gini.rmd(x) - gini.coef(x)
}, numeric(1))))
giniMax

## [1] 2.542411e-14

Over all 34 series of the ca533 data set the two agree to 2.5e-14.