Nonparametric estimation of the measure of functional dependence

1.   Introduction and motivation

The study of association or dependence plays an important role in statistics. One of the important aspects of this is how to measure the strength of various associations among random variables. Among the measures of associations between random variables, Pearson's correlation coefficient, Spearman's $\rho$, and Kendall's $\tau$ are the most prominent ones. But they only measure linear or monotonic relationships, not suitable for a general nonlinear relationship. For example, when the relationship between two random variables is parabolic-shaped, none of the measures above will be applicable. Thus, a measure addressing this issue is desirable.

Two random variables are said to be mutually complete dependent (MCD) when they have mutual functional relationship. This concept was first introduced by Lancaster ^[1]. This is also known as the strongest dependence. In this relation, one variable is completely predictable of the other. Siburg and Stoimenov ^[2] constructed a measure of MCD for continuous random variables. Tasena and Dhompongsa ^[3] extended this measure to the multivariate case. In their papers, they proposed to measure the distance between two copulas by a modified Sobolev norm. For details about inner product and Sobolev norm on copula space, we refer to Darsow and Olsen ^[4].

Dette et al. ^[5] found that a simple modification of the measure of MCD can be used to measure the strength of functional dependence (MFD). Note that this functional dependence includes a wider range of dependence since it could be nonlinear or even nonmonotonic. The discrete form of MFD was given in Shan et al. ^[6]. Given the theoretical definition of MFD, the next question will be how to estimate it. Similar to the measure of MCD, MFD is also constructed based on copula. So a straightforward way of estimating MFD contains two steps. First, estimating the copula or its density. Second, using the estimated copula or its density to estimate the MFD. Generally in step one, there are two approaches to estimate the copula or its density, parametric way or nonparametric way. In the parametric way, one assumes a parametric model for copula, and estimates its parameters by the method of maximum likelihood (MLE). However, unlike distribution functions, the dependence structure, i.e., the copula, is usually hidden behind the data set, which makes the claim of having prior knowledge of copula family quite questionable. So in this paper, we will mainly consider the nonparametric way.

In this article, we propose estimating the MFD using the kernel method with the beta kernel. We introduce the new estimators and study their asymptotic properties in Section 3. Using Monte Carlo simulations, we investigate the finite sample performance of the proposed estimators relative to traditional measures of dependence. Our simulation results, reported in Section 4, show that the new estimators are accurate and stable with the choice of different parameters in a given model. In Section 5, a real data example is used to illustrate the new estimators. The paper is concluded by some remarks in Section 6.

2.   Preliminaries

This section introduces some notations and concepts that will be used in the remainder of the article. We will focus on bivariate continuous distributions. Consider an independent and identically distributed sample $(X_1, Y_1), \cdots, (X_n, Y_n)$ of a bivariate random vector $(\bf{X}, \bf{Y})$ with joint distribution function $H$, and marginal distribution functions $F$ and $G$, respectively. Then, by Sklar's Theorem ^[7], there exists a unique copula $C: I^2 = [0, 1]^2\to I = [0, 1]$ such that

Therefore, the copula density is given by

The measure of MCD is based on the norm of functions in the copula space $\mathbb{C}$. The Sobolev norm for copula takes the following form

The above norm has the following properties.

Proposition 1. The Sobolev norm for copulas satisfies $\| C \|^2 \in [2 / 3, 1]$ for all $C \in \mathbb{C}$. Moreover, the following properties hold:

i. $\| C \|^2 = 2 / 3$ if and only if $X$ and $Y$ are independent.

ii. $\| C \|^2 = 1$ if and only if $X$ and $Y$ are MCD.

Inspired by this proposition, Siburg and Stoimenov ^[2] proposed the measure of MCD, which can be defined as follows.

Definition 1. Given two continuous random variables $X, Y$ with copula $C$, we define

$\rho (X, Y)$ can be interpreted as a normalized Sobolev distance of $C$ from the independent copula denoted by $\Pi$:

where $C_m$ is a MCD copula. A close look at the MCD shows that it can be decomposed into two opposite functional dependencies, i.e., $Y$ is a function of $X$ and $X$ is a function of $Y$. Therefore, the measure of functional dependence (MFD) can be derived by modifying the measure for MCD, which is the main idea of Dette et al. ^[5]. Its discrete form was discussed in Shan et al. ^[6]. The construction of MFD is based on the following propositions.

Proposition 2. Let $X$ and $Y$ be two random variables with copula $C$. Then,

$i. X$ and $Y$ are independent if and only if $\partial_1 C_{X, Y} (u, v) = v$ for Lebesgue almost all $(u, v) \in I^2$.

$ii. Y$ is almost surely (a.s.) a Borel function of $X$ if and only if $\partial_1 C_{X, Y} (u, v) \in \{ 0, 1 \}$ for Lebesgue almost all $(u, v) \in I^2$.

Proposition 3. For any $C_{X, Y} \in \mathbb{C}$, we have $\| \partial_1 C_{X, Y} \|^2_2 \in [1 / 3, 1 / 2]$, Moreover,

$i. \| \partial_1 C_{X, Y} \|^2_2 = 1 / 3$ if and only if $X$ and $Y$ are independent.

$ii. \| \partial_1 C_{X, Y} \|^2_2 = 1 / 2$ if and only if $Y$ is a.s. a Borel function of $X$.

Notice that $\| \partial_1 C_{X, Y} \|$ reaches its boundaries at two extreme cases. $\| \partial_2 C_{X, Y} \|$ has similar propositions which we will not reproduce here. Using these definitions and propositions, we can define

This is a standardized form of $\| \partial_1 C_{X, Y} \|$. Hence $\rho_1$ inherits similar properties:

i. $\rho_1 = 0$ if and only if $X$ and $Y$ are independent.

ii. $\rho_1 = 1$ if and only if $Y$ is a.s. a Borel function of $X$.

Similarly, $\rho_2 = 1$ indicates $X$ is a.s. a Borel function of $Y$. Those properties suggest that $\rho_i (i = 1, 2)$ can be used to measure functional dependence. We can assess the strength of dependence from the magnitude of $\rho_i$.

3.   Nonparametric estimation

The measures $\rho_1$, $\rho_2$ and $\rho$ are all constructed based on copulas. They both can be stated in terms of copulas or copula densities. In other words, (2.5) can be written as

and (2.6) can be written as

Accordingly, there are two approaches to estimate the measures, through copula or its density. We will focus on the latter in this paper.

The estimation of the copula density has been discussed in many papers. For example, Kauermann et al. ^[8] estimated copula density with B-spline. Genest et al. ^[9] estimated copula density through wavelets. It is recognized that the estimation of copula density involves more technical difficulties than usual density estimation. One of the big issues in this respect is boundary bias. Several methods have been proposed to address this issue. Omelka et al. ^[10] suggested an improved version of mirror-reflection estimator. Charpentier et al. ^[11] suggested to use transformation estimator. Chen ^[12] suggested to use beta kernels whose support matches the support of copulas. Geenens et al. ^[13] used probit transformation to reduce boundary bias effect in kernel estimation of copula density. Majdara and Nooshabadi ^[14] provides a novel method in estimating copula density in high dimension space.

Dette et al. ^[5] discussed the asymptotic behavior of the estimation $\rho_1$ based on symmetric kernels and Eq (2.5). In this paper, we will use the beta kernel in estimating $\rho_1$ and $\rho_2$. Beta kernel smoothing was considered by Harrell and Davis ^[15]. Chen ^[12], Chen ^[16] applied beta kernel smoothing in density estimation, and found that the beta estimator can reduce boundary bias and variance compared with local linear estimators for densities with finite support. Following the same idea, Charpentier et al. ^[11] proposed the beta kernel based estimator for copula density.

3.1.   Nonparametric estimation of the copula density

Let $(X_1, Y_1), (X_2, Y_2), \cdots (X_n, Y_n)$ be a sample from $H_{X, Y} (x, y)$ with unknown marginals. Denote the copula corresponding to $H_{X, Y} (x, y)$ by $C (u, v)$. We assume that both $H$ and $C$ are completely unknown. We start by the most convenient situation in which we assume that the copula $C$ is twice differentiable. Let $c$ denote the density of the copula. Usually copulas are estimated via pseudo-observations $(\hat{F} (X_i), \hat{G}(Y_i))$, where $\hat{F}$ and $\hat{G}$ are the empirical distribution functions, i.e.,

with $\mathbb{1}(A)$ being the usual indicator function. The beta kernel based estimator of the copula density at point $(u, v) \in [0, 1]^2$ is

where $K (\cdot, \alpha, \beta)$ is the density of the beta distribution with parameters $\alpha$ and $\beta$, i.e.,

with $B (\alpha, \beta) = \Gamma (\alpha + \beta)/\Gamma (\alpha) \Gamma (\beta)$. For convenience, the same bandwidth is used in both kernels in $\hat{c}_h (u, v)$. Charpentier et al. ^[11] claimed the asymptotic normality of $\hat{c}_h (u, v)$by showing that, for all $(u, v) \in [0, 1]^2$,

as $nh \rightarrow \infty$ and $h \rightarrow 0$, where "$\overset{L}{\rightarrow }$" means convergence in distribution. Nagler ^[17] provided detailed proof and gave the bias and variance of $\hat{c}_h (u, v)$ in Proposition 4. He also discussed bandwidth selection for $\hat{c}_h (u, v)$.

Proposition 4. Let $c (u, v)$ be twice continuously differentiable on $(0, 1)^2$, and $h_n \rightarrow 0$ and $n h_n \rightarrow \infty$ as $n \rightarrow \infty$. Then, for all $(u, v) \in (0, 1)^2$,

Note that there is a little difference between Proposition 4, which is on the interior of $(0, 1)^2$, and Charpentier's claim ^[11], which is on the whole $[0, 1]^2$. Since the definition of functional dependence is based on the integration of copula densities, which will not be affected by the boundaries. We will consider the whole $[0, 1]^2$ in the remainder.

3.2.   Nonparametric estimation of MFD via the copula density

Since all of the three measures, $\rho_1$, $\rho_2$ and $\rho$, are constructed in similar manner, we only take $\rho_1$ as an example and show how to estimate it through estimating the copula density using the beta kernel. Let $\ell^{\infty} ([0, 1]^2)$ be the space of the collection of all uniformly bounded real-valued functions defined on $[0, 1]^2$, equipped with the uniform metric $m$ defined as

Define $\phi_i :$ $\ell^{\infty} ([0, 1]^2) \rightarrow \mathbb{R}$, $i = 1, 2$, by

Then, the three measures, $\rho^2_1$, $\rho^2_1$ and $\rho^2$, are functionals of $c(\cdot, \cdot)$. So, it suffices to show that $\phi_1$ and $\phi_2$ are Hadamard differentiable.

Theorem 1. Let $c (u, v)$ be twice continuously differentiable on $[0, 1]^2$, and $h \rightarrow 0$ and $n h \rightarrow \infty$ as $n \rightarrow \infty$. Then,

where $\phi'_1 (l (u, y)) = \int_0^1 \int_0^1 \left(2 \int_0^v c (u, y) dy \int_0^v l (u, y) dy \right) dudv$.

Proof. $\rho_1$, $\rho_2$ can be represented as a map $\phi_1$, $\phi_2$: $\ell^{\infty} ([0, 1]^2) \rightarrow \mathbb{R} ([0, 1])$ via $\rho_1 = \phi_1 (c)$ and $\rho_2 = \phi_2 (c)$, respectively. The function space $\ell^{\infty} ([0, 1]^2)$ is equipped with the uniform metric $m$. For all converging sequences $t_n \rightarrow 0$ and $l_n \rightarrow l$ such that $c + t_n l_n \in \ell^{\infty} ([0, 1]^2)$ for every $n$, we have

So, the Hadamard derivative of $\phi_1$ at $c$ is

Therefore, according to the Delta method ^[18],

This completes the proof.

The asymptotic distributions of $\hat{\rho}^2_1$ and $\hat{\rho}^2$ can be derived in exactly the same manner, so we omit their details.

3.3.   Nonparametric estimation of MFD via copula functions

In the following, we show that estimators of MFD through copula have the same asymptotic distributions as those established through copula density. As an example, let's check the asymptotic distribution of $\hat{\rho}^2$.

Mapping a copula density to an MFD can be decomposed into two steps as follows

The first map is a double-integration that is linear and continuous, and thus, it is Hadamard-differentiable. We only need to check the second map.

Let $\mathbb{D}^1_2 ([0, 1]^2)$ be the Sobolev space, $D_1, D_2 \in \mathbb{D}^1_2 ([0, 1]^2)$, and define the inner product

where $\bigtriangledown$ is gradient. The Sobolev norm induced by the inner product in $\mathbb{D}^1_2 ([0, 1]^2)$ is

for $D \in \mathbb{D}^1_2 ([0, 1]^2)$.

Let $\mathfrak{C} \subset \mathbb{D}^1_2 ([0, 1]^2)$ be the copula space and $C \in \mathfrak{C}$ is a copula. Define $\psi : \mathfrak{C} \mapsto \mathbb{R}$ by $\psi (C) = |C|^2$. Then, the derivative of $\psi$ at $C$ along $D$ is:

The last step follows from Theorem 2.3 in ^[19]. This result shows the convergence of $\hat{\rho}^2$, and straightforward calculations will show that it is consistent with the asymptotic distribution in Theorem 1.

4.   Implementation and simulations

4.1.   Choosing the evaluation grid

In copula density estimation, Nagler ^[20] suggested using a grid that is equally spaced after a transformation by the inverse Gaussian cdf, which is shown in Figure 1. Our simulation results below show that evaluating copula density at a set of grid points in a similar pattern will improve the accuracy of estimators of MCD. To compare the impact of the choice of grid, we considered two copula families, the Gaussian copula with parameters ${0, 0.1, 0.2, 0.5, 0.8, 0.9}$ and the Gumbel copula with parameters ${1, 10/9, 10/7, 10/3, 5, 10}$. Two samples of sizes 200 and 1000 were taken from each copula, respectively. First, copula densities are estimated from each sample based on the KDEcopula package. Then, the estimated copula density was evaluated on two sets of grid points: the usual grid with equally spaced points and a normalized grid. From the discretized copula density, we calculate the estimate of MFD. Figures 2 and 3 show the mean absolute error (MAE) of estimators with sample size 200 and 1000 under 500 replication. In each case, we find the MAE of estimators based on an equally spaced grid, labeled "equal", are significantly higher than the MAE of the same estimators based on the transformed grid, labeled "norm".

4.2.   Simulation

In this section, we explore the finite sample performance of the proposed estimators using the mean squared error (MSE). To put all estimators on the same scale, we standardize MFD. In other words, we use $\rho_1$, $\rho_2$ and $\rho$. The corresponding estimators will be denoted by $\hat{\rho}_1 (Y \mid X)$, $\hat{\rho}_2 (X \mid Y)$ and $\hat{\rho} (X, Y)$. For two-dimensional density estimation, using cross-validation to choose the bandwidth is computationally expensive. Therefore, in all simulations, a rule-of-thumb bandwidth is used. More precisely, the bandwidth is selected based on the asymptotic mean integrated squared error (AMISE)-optimality with respect to the Frank copula. For further details on bandwidth selection in this context, we refer to Nagler ^[20]. In all simulations reported here, the integration is calculated over a grid of $30 \times 30$ points. For the choice of the grid, we adopt the method in Nagler ^[20]. That is, we apply the Gaussian cumulative distribution function to equally spaced 30 knots on a line segment $[-3, 3]$. The final two-dimensional grid is shown in Figure 1. This choice takes into account the fact that copula densities usually have high fluctuation on the boundary and corners. Putting more evaluating points on those regions will reduce approximation errors.

In Tables 1 and 2, we present the simulated MSE of the estimators of $\rho_1$ and $\rho_2$ for samples size 50, 100 and 200 for Gaussian copula. These results are based on 1000 replications. Both copulas are generated by R package "copula". We find that our estimators have reasonable precision in all cases. As the sample size increases, MSE is getting smaller. And there is no significant difference in MSE for different $\theta$ values, which indicates our estimator is stable for the choice of $\theta$'s.

4.3.   Comparison of measures

In the second part of the simulation, we will compare the performance of MFDs with other measures of dependence, e.g., linear correlation coefficient $r$, Spearman's $\rho$ and Kendall's $\tau$ under several different types of relationships. We choose three different dependence structures: elliptical distributions, monotonic dependence and regressional dependence, represented by normal copula, cubic function and quadratic function, respectively.

The first example is a quadratic function. 500 data are generated from the following model,

where $\varepsilon \sim N (0, \sigma)$ and $\sigma = 1$, 5, and 10, respectively (see Figure 4). To obtain the copula data, we apply the empirical marginal distributions to the data, i.e., apply a rank transformation as shown in (3.3) to the data generated by the model (4.1). Then beta kernel estimation is applied to get the estimations of $\rho_1$ and $\rho_2$.

Obviously, neither Spearman's $\rho$ nor Kendall's $\tau$ is suitable for this situation. The simulation results in Table 3 also showed that they are almost $0$ in all cases. $\hat{\rho} (X, Y)$, on the other hand, is much higher than both Spearman's $\rho$ and Kendall's $\tau$, especially for small $\sigma$. This indicates that the type of dependence is functional, not monotonic. And the magnitude of $\hat{\rho} (X, Y)$ tells that the strength of dependence is getting weaker as $\sigma$ increases. In a comparison of $\hat{\rho}_1 (Y \mid X)$ and $\hat{\rho}_2 (X \mid Y)$, we find that functional dependence is stronger in the $Y$ to $X$ direction than the other direction since $\hat{\rho}_1 (Y \mid X)$ is higher. Again, as $\sigma$ increases, the strength of dependence in this direction is also getting weaker.

To compare the performance of MFD with Kendall's $\tau$ and Spearman's $\rho$ in monotonic dependence. 500 data are generated from the following model,

where $\varepsilon \sim N (0, \sigma)$ and $\sigma = 1$, 5 and 10. The scatter plot of model (4.2) is in Figure 5 and the simulation results are in Table 4. As shown in Table 4, the values of MFD has a similar decreasing pattern as the other two measures when $\sigma$ increases. Indeed, a cubic function is one type of functional dependence, so MFDs are capable of measuring the strength of monotonic dependence. The values of $\hat{\rho}_1 (Y \mid X)$ and $\hat{\rho}_2 (X \mid Y)$, which measure the strength of functional dependence in two directions ($Y$ to $X$ and $X$ to $Y$) separately, are close to each other, obviously this is because model (4.2) is symmetric.

Next, we take into account the Pearson's correlation coefficient. We take normal copulas as an example, which is

with $r = 0.1, 0.5~ {\operatorname{{and}}} ~0.9$. The scatter plots of Gaussian copulas are in Figure 6. Table 5 shows the simulation results. As expected, those measures show no significant difference in measuring the dependence of elliptical distributions.

The comparison of MFDs with other measures in models 4.1–4.3 shows that MFDs have good adaptability for different types of relationships.

5.   Application

The measurement of functional relationships has many possible applications (^[21,22]). The Communities and Crime Data Set ^[23] contains community crime rate of 1994 communities with 123 possibly related variables. We will use functional dependence measure as a criteria for variable selection to choose the variables which have the most impact on community crime rate. We calculate the measures $\rho$, $\rho_1$, and $\rho_2$ given in Eqs (2.4)–(2.6), respectively by using beta kernel estimation for community crime rate and each of other variables. Variables with higher values of the measures have a greater impact on community crime rate. The Table 6 shows 9 variables with highest functional dependence measures and Table 7 gives the explanation of the abbreviations. Notice that the measures can detect strong non-linear relationships. As shown in Figure 7, two of the selected variables (PctIlleg and racePctWhite), showed clear nonlinear relations with the crime rate.

6.   Discussion

This paper showed that, compared with Spearman's $\rho$ or Kendall's $\tau$, the measures of functional relationship could not only measure the strength of a relationship but also indicate the direction of a possible functional relationship. We provide a novel method to estimate the measures of functional relationships. The simulation results showed that they have fairly good accuracy.

Although MFD can quantify the strength of functional dependence, it doesn't suggest any specific form of the function. So one possible application of this measure is in variable selection. We use MFD to filter out the less correlated variables, then use parametric or nonparametric methods to construct a predicting model. In the community crime data example, we showed that MFD could detect nonlinear relationship, but as for how many variables should be retained, in other words, how to set up the threshold of MFD in variable selection is a question that needs to be discussed and may involve some subjective opinions. After the desired number of variables are chosen, people may use either parametric or non-parametric methods to set up the model.

Acknowledgments

This work was supported in part by the Education Department of Jiangxi Province under Grant GJJ190253 and Grant GJJ190259.

Conflict of interest

The authors have no conflicts of interest to declare.