`ghcm`

is an R package used to perform conditional independence tests for functional data.

This vignette gives a brief overview of the usage of the `ghcm`

package. We give a brief presentation of the idea behind the GHCM and the conditions under which the test is valid. Subsequently, we provide several examples of the usage of the `ghcm`

package by analysing a simulated data set.

In this section we briefly describe the idea behind the GHCM. For the full technical details and theoretical results, see [1].

Let \(X\), \(Y\) and \(Z\) be random variables of which we are given \(n\) i.i.d. observations \((X_1, Y_1, Z_1), \dots, (X_n, Y_n, Z_n)\); here \(X\), \(Y\) and \(Z\) can be either scalar or functional. Existing methods, such as the GCM [2] implemented in the `GeneralisedCovarianceMeasure`

package [3], can deal with most cases where both \(X\) and \(Y\) are scalar hence our primary interest is in the cases where at least one of \(X\) and \(Y\) are functional.

The GHCM estimates the squared Hilbertâ€“Schmidt norm of the expected conditional covariance of \(X\) and \(Y\) given \(Z\), \(\| \mathscr{K} \|_{\mathrm{HS}}^2\), and rejects the hypothesis \(X \mbox{${}\perp\mkern-11mu\perp{}$}Y \,|\,Z\) if this quantity is too large. We denote by \(\langle f, g \rangle\) the inner product between \(f\) and \(g\). We compute the GHCM as follows.

- Regress \(X\) on \(Z\) and \(Y\) on \(Z\) yielding residuals \(\hat{\varepsilon}\) and \(\hat{\xi}\), respectively.
- Compute the test statistic \[ T = \frac{1}{n} \sum_{i=1}^n \sum_{j=1}^n \langle \hat{\varepsilon}_i, \hat{\varepsilon}_j \rangle \langle \hat{\xi}_i, \hat{\xi}_j \rangle. \]
- Compute the \(n \times n\) matrix \(\Gamma\) by setting \[ \Gamma_{ij} := \langle \hat{\varepsilon}_i, \hat{\varepsilon}_j \rangle \langle \hat{\xi}_i, \hat{\xi}_j \rangle. \]
- Compute the \(d \leq n-1\) non-zero eigenvalues \((\lambda_k)_{k=1}^d\) of \[ C := \frac{1}{n-1} (\Gamma - J \Gamma - \Gamma J + J \Gamma J), \] where \(J\) is the \(n \times n\) matrix with all entries equal to \(1/n\).
- Produce a \(p\)-value by setting \[ p := \mathbb{P}\left(\sum_{k=1}^d \lambda_k \zeta_k^2 > T\right), \] where \(\zeta_1, \dots, \zeta_d\) are i.i.d. standard Gaussian random variables.

Assuming that the regression methods perform sufficiently well, the GHCM has approximately uniformly distributed \(p\)-values when the null is true. It should be noted that there are situations where \(X \mbox{${}\not\!\perp\mkern-11mu\perp{}$}Y \,|\,Z\) but the GHCM is unable to detect this dependence for any sample size, since \(\mathscr{K}\) can be zero in this case.

The GHCM as implemented in the function uses different methods to compute the inner products as required above depending on the format of the given residuals. For residuals from scalar or multivariate variables or functional variables observed on a constant, fixed grid, the inner products are computed as the usual Euclidean inner products and no further preprocessing is done. For residuals coming from irregularly observed functional variables or functions on a fixed grid with missing values, the inner products are computed as L2 inner products computed from spline interpolants of the curves. See the following sections for examples of both types of use.

To give concrete examples of the usage of the package, we perform conditional independence tests on a simulated data set consisting of both functional and scalar variables. The functional variables are observed on a common equidistant grid of \(101\) points on \([0, 1]\).

```
library(ghcm)
set.seed(111)
data(ghcm_sim_data)
<- seq(0, 1, length.out=101)
grid colnames(ghcm_sim_data)
#> [1] "Y_1" "Y_2" "X" "Z" "W"
```

`ghcm_sim_data`

consists of 500 observations of the scalar variables \(Y_1\) and \(Y_2\) and the functional variables \(X\), \(Z\) and \(W\). The curves and the mean curve for functional data can be seen in Figures 1, 2 and 3.

In all of the upcoming examples we will use functions from the `refund`

R-package [4] to perform regressions. The simulated data is obtained such that the regressions we perform are well-specified hence we do not have to worry about the performance of our regression methods. However, in actual applications of the GHCM, it is critical that the regression methods employed estimate the conditional expectations \(\mathbb{E}(X \,|\,Z)\) and \(\mathbb{E}(Y \,|\,Z)\) sufficiently well for the \(p\)-values to be valid. Any use of the GHCM should be prefaced by an analysis of the performance of the regression methods in use.

We first test whether \(Y_1\) and \(Y_2\) are conditionally independent given the functional variables. This is relevant if, say, weâ€™re trying to predict \(Y_1\) and we want to know whether including \(Y_2\) as a predictor would be helpful. A naive correlation-based approach would suggest that \(Y_2\) could be relevant since:

```
cor(ghcm_sim_data$Y_1, ghcm_sim_data$Y_2)
#> [1] 0.5762141
```

To perform the conditional independence test, we need a scalar-on-function regression method and we will use the `pfr`

function from the `refund`

package [4] with `lf`

-terms. We run the test in the following code:

```
library(refund)
<- pfr(Y_1 ~ lf(X) + lf(Z) + lf(W) , data=ghcm_sim_data)
m_1 <- pfr(Y_2 ~ lf(X) + lf(Z) + lf(W), data=ghcm_sim_data)
m_2 <- ghcm_test(resid(m_1), resid(m_2))
test print(test)
#> H0: X _||_ Y | Z, p: 0.9583592
#> Not rejected at 5 % level
```

We obtain a \(p\)-value of 0.958. It should be noted that since the asymptotic distribution of the test statistic depends on the underlying distribution, there is no way to know the \(p\)-value from just the test statistic alone, hence it is not reported.

We now test whether \(Y_1 \mbox{${}\perp\mkern-11mu\perp{}$}X \,|\,Z\). This is relevant if weâ€™re interested in modelling \(Y_1\) and want to determine whether \(X\) should be included in a model that already includes \(Z\). We can plot \(X\) and \(Z\) and color the curves based on the value of \(Y_1\) as can be seen in Figures 4 and 5 below.