Estimation methods
Given the form of bridge function \(F\), obtaining a moment-based estimation
\(\widehat \sigma_{jk}\) requires
inversion of \(F\).
latentcor implements two methods for calculation of the
inversion:
method = "original"Subsection describing original method and relevant parametertolmethod = "approx"Subsection describing approximation method and relevant parameterratio
Both methods calculate inverse bridge function applied to each
element of sample Kendall’s \(\tau\)
matrix. Because the calculation is performed point-wise (separately for
each pair of variables), the resulting point-wise estimator of
correlation matrix may not be positive semi-definite.
latentcor performs projection of the pointwise-estimator to
the space of positive semi-definite matrices, and allows for shrinkage
towards identity matrix using the parameter nu (see Subsection describing adjustment of point-wise
estimator and relevant parameter nu).
Original method (method = "original")
Original estimation approach relies on numerical inversion of \(F\) based on solving uni-root optimization
problem. Given the calculated \(\widehat
\tau_{jk}\) (sample Kendall’s \(\tau\) between variables \(j\) and \(k\)), the estimate of latent correlation
\(\widehat \sigma_{jk}\) is obtained by
calling optimize function to solve the following
optimization problem: \[
\widehat r_{jk} = \arg\min_{r} \{F(r) - \widehat \tau_{jk}\}^2.
\] The parameter tol controls the desired accuracy
of the minimizer and is passed to optimize, with the
default precision of 1e-8:
estimate = latentcor(X, types = c("con", "bin", "ter", "tru"), method = "original", tol = 1e-8)Algorithm for Original method
Input: \(F(r)=F(r, \mathbf{\Delta})\) - bridge function based on the type of variables \(j\), \(k\)
- Step 1. Calculate \(\hat{\tau}_{jk}\) using (1).
estimate$K
#> [,1] [,2] [,3] [,4]
#> [1,] 1.0000000 0.2557576 0.2456566 0.3331313
#> [2,] 0.2557576 1.0000000 0.1555556 0.2339394
#> [3,] 0.2456566 0.1555556 1.0000000 0.2183838
#> [4,] 0.3331313 0.2339394 0.2183838 1.0000000- Step 2. For binary/truncated variable \(j\), set \(\hat{\mathbf{\Delta}}_{j}=\hat{\Delta}_{j}=\Phi^{-1}(\pi_{0j})\) with \(\pi_{0j}=\sum_{i=1}^{n}\frac{I(x_{ij}=0)}{n}\). For ternary variable \(j\), set \(\hat{\mathbf{\Delta}}_{j}=(\hat{\Delta}_{j}^{1}, \hat{\Delta}_{j}^{2})\) where \(\hat{\Delta}_{j}^{1}=\Phi^{-1}(\pi_{0j})\) and \(\hat{\Delta}_{j}^{2}=\Phi^{-1}(\pi_{0j}+\pi_{1j})\) with \(\pi_{0j}=\sum_{i=1}^{n}\frac{I(x_{ij}=0)}{n}\) and \(\pi_{1j}=\sum_{i=1}^{n}\frac{I(x_{ij}=1)}{n}\).
estimate$zratios
#> [[1]]
#> [1] NA
#>
#> [[2]]
#> [1] 0.5
#>
#> [[3]]
#> [1] 0.3 0.8
#>
#> [[4]]
#> [1] 0.5- Compute \(F^{-1}(\hat{\tau}_{jk})\)
as \(\hat{r}_{jk}=argmin\{F(r)-\hat{\tau}_{jk}\}^{2}\)
solved via
optimizefunction in R with accuracytol.
estimate$Rpointwise
#> [,1] [,2] [,3] [,4]
#> [1,] 1.0000000 0.5529903 0.4480984 0.5826171
#> [2,] 0.5529903 1.0000000 0.4050223 0.5821513
#> [3,] 0.4480984 0.4050223 1.0000000 0.4653875
#> [4,] 0.5826171 0.5821513 0.4653875 1.0000000Approximation method (method = "approx")
A faster approximation method is based on multi-linear interpolation of pre-computed inverse bridge function on a fixed grid of points (Yoon, Müller, and Gaynanova 2021). This is possible as the inverse bridge function is an analytic function of at most 5 parameters:
- Kendall’s \(\tau\)
- Proportion of zeros in the 1st variable
- (Possibly) proportion of zeros and ones in the 1st variable
- (Possibly) proportion of zeros in the 2nd variable
- (Possibly) proportion of zeros and ones in the 2nd variable
In short, d-dimensional multi-linear interpolation uses a weighted
average of \(2^{d}\) neighbors to
approximate the function values at the points within the d-dimensional
cube of the neighbors, and to perform interpolation,
latentcor takes advantage of the R package
chebpol (Gaure 2019). This
approximation method has been first described in (Yoon, Müller, and Gaynanova 2021) for
continuous/binary/truncated cases. In latentcor, we
additionally implement ternary case, and optimize the choice of grid as
well as interpolation boundary for faster computations with smaller
memory footprint.
#estimate = latentcor(X, types = c("con", "bin", "ter", "tru"), method = "approx")Algorithm for Approximation method
Input: Let \(\check{g}=h(g)\), pre-computed values \(F^{-1}(h^{-1}(\check{g}))\) on a fixed grid \(\check{g}\in\check{\cal{G}}\) based on the type of variables \(j\) and \(k\). For binary/continuous case, \(\check{g}=(\check{\tau}_{jk}, \check{\Delta}_{j})\); for binary/binary case, \(\check{g}=(\check{\tau}_{jk}, \check{\Delta}_{j}, \check{\Delta}_{k})\); for truncated/continuous case, \(\check{g}=(\check{\tau}_{jk}, \check{\Delta}_{j})\); for truncated/truncated case, \(\check{g}=(\check{\tau}_{jk}, \check{\Delta}_{j}, \check{\Delta}_{k})\); for ternary/continuous case, \(\check{g}=(\check{\tau}_{jk}, \check{\Delta}_{j}^{1}, \check{\Delta}_{j}^{2})\); for ternary/binary case, \(\check{g}=(\check{\tau}_{jk}, \check{\Delta}_{j}^{1}, \check{\Delta}_{j}^{2}, \check{\Delta}_{k})\); for ternary/truncated case, \(\check{g}=(\check{\tau}_{jk}, \check{\Delta}_{j}^{1}, \check{\Delta}_{j}^{2}, \check{\Delta}_{k})\); for ternay/ternary case, \(\check{g}=(\check{\tau}_{jk}, \check{\Delta}_{j}^{1}, \check{\Delta}_{j}^{2}, \check{\Delta}_{k}^{1}, \check{\Delta}_{k}^{2})\).
Step 1 and Step 2 same as Original method.
Step 3. If \(|\hat{\tau}_{jk}|\le \mbox{ratio}\times \bar{\tau}_{jk}(\cdot)\), apply interpolation; otherwise apply Original method.
To avoid interpolation in areas with high approximation errors close
to the boundary, we use hybrid scheme in Step 3. The parameter
ratio controls the size of the region where the
interpolation is performed (ratio = 0 means no
interpolation, ratio = 1 means interpolation is always
performed). For the derivation of approximate bound for BC, BB, TC, TB,
TT cases see Yoon, Müller, and Gaynanova
(2021). The derivation of approximate bound for NC, NB, NN, NT
case is in the Appendix.
\[ \bar{\tau}_{jk}(\cdot)= \begin{cases} 2\pi_{0j}(1-\pi_{0j}) & for \; BC \; case\\ 2\min(\pi_{0j},\pi_{0k})\{1-\max(\pi_{0j}, \pi_{0k})\} & for \; BB \; case\\ 2\{\pi_{0j}(1-\pi_{0j})+\pi_{1j}(1-\pi_{0j}-\pi_{1j})\} & for \; NC \; case\\ 2\min(\pi_{0j}(1-\pi_{0j})+\pi_{1j}(1-\pi_{0j}-\pi_{1j}),\pi_{0k}(1-\pi_{0k})) & for \; NB \; case\\ 2\min(\pi_{0j}(1-\pi_{0j})+\pi_{1j}(1-\pi_{0j}-\pi_{1j}), \\ \;\;\;\;\;\;\;\;\;\;\pi_{0k}(1-\pi_{0k})+\pi_{1k}(1-\pi_{0k}-\pi_{1k})) & for \; NN \; case\\ 1-(\pi_{0j})^{2} & for \; TC \; case\\ 2\max(\pi_{0k},1-\pi_{0k})\{1-\max(\pi_{0k},1-\pi_{0k},\pi_{0j})\} & for \; TB \; case\\ 1-\{\max(\pi_{0j},\pi_{0k},\pi_{1k},1-\pi_{0k}-\pi_{1k})\}^{2} & for \; TN \; case\\ 1-\{\max(\pi_{0j},\pi_{0k})\}^{2} & for \; TT \; case\\ \end{cases} \]
By default, latentcor uses ratio = 0.9 as
this value was recommended in Yoon, Müller, and
Gaynanova (2021) having a good balance of accuracy and
computational speed. This value, however, can be modified by the
user.
#latentcor(X, types = c("con", "bin", "ter", "tru"), method = "approx", ratio = 0.99)$R
#latentcor(X, types = c("con", "bin", "ter", "tru"), method = "approx", ratio = 0.4)$R
latentcor(X, types = c("con", "bin", "ter", "tru"), method = "original")$R
#> [,1] [,2] [,3] [,4]
#> [1,] 1.0000000 0.5524373 0.4476503 0.5820345
#> [2,] 0.5524373 1.0000000 0.4046173 0.5815691
#> [3,] 0.4476503 0.4046173 1.0000000 0.4649222
#> [4,] 0.5820345 0.5815691 0.4649222 1.0000000The lower is the ratio, the closer is the approximation
method to original method (with ratio = 0 being equivalent
to method = "original"), but also the higher is the cost of
computations.
library(microbenchmark)
#microbenchmark(latentcor(X, types = c("con", "bin", "ter", "tru"), method = "approx", ratio = 0.99)$R)
#microbenchmark(latentcor(X, types = c("con", "bin", "ter", "tru"), method = "approx", ratio = 0.4)$R)
microbenchmark(latentcor(X, types = c("con", "bin", "ter", "tru"), method = "original")$R)
#> Unit: milliseconds
#> expr
#> latentcor(X, types = c("con", "bin", "ter", "tru"), method = "original")$R
#> min lq mean median uq max neval
#> 39.4301 40.65735 41.72452 41.0003 42.33695 50.8438 100Rescaled Grid for Interpolation
Since \(|\hat{\tau}|\le \bar{\tau}\), the grid does not need to cover the whole domain \(\tau\in[-1, 1]\). To optimize memory associated with storing the grid, we rescale \(\tau\) as follows: \(\check{\tau}_{jk}=\tau_{jk}/\bar{\tau}_{jk}\in[-1, 1]\), where \(\bar{\tau}_{jk}\) is as defined above.
In addition, for ternary variable \(j\), it always holds that \(\Delta_{j}^{2}>\Delta_{j}^{1}\) since \(\Delta_{j}^{1}=\Phi^{-1}(\pi_{0j})\) and \(\Delta_{j}^{2}=\Phi^{-1}(\pi_{0j}+\pi_{1j})\). Thus, the grid should not cover the the area corresponding to \(\Delta_{j}^{2}\le\Delta_{j}^{1}\). We thus rescale as follows: \(\check{\Delta}_{j}^{1}=\Delta_{j}^{1}/\Delta_{j}^{2}\in[0, 1]\); \(\check{\Delta}_{j}^{2}=\Delta_{j}^{2}\in[0, 1]\).
Speed Comparison
To illustrate the speed improvement by
method = "approx", we plot the run time scaling behavior of
method = "approx" and method = "original"
(setting types for gen_data by replicating
c("con", "bin", "ter", "tru") multiple times) with
increasing dimensions \(p = [20, 40, 100, 200,
400]\) at sample size \(n =
100\) using simulation data. Figure below summarizes the observed
scaling in a log-log plot. For both methods we observe the expected
\(O(p^2)\) scaling behavior with
dimension p, i.e., a linear scaling in the log-log plot. However,
method = "approx" is at least one order of magnitude faster
than method = "original" independent of the dimension of
the problem.
Adjustment of pointwise-estimator for positive-definiteness
Since the estimation is performed point-wise, the resulting matrix of estimated latent correlations is not guaranteed to be positive semi-definite. For example, this could be expected when the sample size is small (and so the estimation error for each pairwise correlation is larger)
set.seed("234820")
X = gen_data(n = 6, types = c("con", "bin", "ter", "tru"))$X
X
#> [,1] [,2] [,3] [,4]
#> [1,] -0.5182800 0 1 0.1021738
#> [2,] -1.3017092 0 0 0.0000000
#> [3,] 0.3145191 1 2 0.4213514
#> [4,] -0.6093291 0 1 1.2771610
#> [5,] -1.3175490 1 0 0.0000000
#> [6,] -0.7807245 1 1 0.0000000
out = latentcor(X, types = c("con", "bin", "ter", "tru"))
out$Rpointwise
#> [,1] [,2] [,3] [,4]
#> [1,] 1.0000000 -0.1477240 0.9990000 0.8548518
#> [2,] -0.1477240 1.0000000 0.3523666 -0.5030324
#> [3,] 0.9990000 0.3523666 1.0000000 0.9114307
#> [4,] 0.8548518 -0.5030324 0.9114307 1.0000000
eigen(out$Rpointwise)$values
#> [1] 2.85954424 1.29130852 0.09944544 -0.25029820latentcor automatically corrects the pointwise estimator
to be positive definite by making two adjustments. First, if
Rpointwise has smallest eigenvalue less than zero, the
latentcor projects this matrix to the nearest positive
semi-definite matrix. The user is notified of this adjustment through
the message (supressed in previous code chunk), e.g.
out = latentcor(X, types = c("con", "bin", "ter", "tru"))Second, latentcor shrinks the adjusted matrix of
correlations towards identity matrix using the parameter \(\nu\) with default value of 0.001
(nu = 0.001), so that the resulting R is
strictly positive definite with the minimal eigenvalue being greater or
equal to \(\nu\). That is \[
R = (1 - \nu) \widetilde R + \nu I,
\] where \(\widetilde R\) is the
nearest positive semi-definite matrix to Rpointwise.
out = latentcor(X, types = c("con", "bin", "ter", "tru"), nu = 0.001)
out$Rpointwise
#> [,1] [,2] [,3] [,4]
#> [1,] 1.0000000 -0.1477240 0.9990000 0.8548518
#> [2,] -0.1477240 1.0000000 0.3523666 -0.5030324
#> [3,] 0.9990000 0.3523666 1.0000000 0.9114307
#> [4,] 0.8548518 -0.5030324 0.9114307 1.0000000
out$R
#> [,1] [,2] [,3] [,4]
#> [1,] 1.0000000 -0.1053533 0.9232992 0.9048072
#> [2,] -0.1053533 1.0000000 0.2372115 -0.4244433
#> [3,] 0.9232992 0.2372115 1.0000000 0.7723678
#> [4,] 0.9048072 -0.4244433 0.7723678 1.0000000As a result, R and Rpointwise could be
quite different when sample size \(n\)
is small. When \(n\) is large and \(p\) is moderate, the difference is
typically driven by parameter nu.
set.seed("234820")
X = gen_data(n = 100, types = c("con", "bin", "ter", "tru"))$X
out = latentcor(X, types = c("con", "bin", "ter", "tru"), nu = 0.001)
out$Rpointwise
#> [,1] [,2] [,3] [,4]
#> [1,] 1.0000000 0.5528213 0.4476819 0.5823120
#> [2,] 0.5528213 1.0000000 0.4058967 0.5808889
#> [3,] 0.4476819 0.4058967 1.0000000 0.4567771
#> [4,] 0.5823120 0.5808889 0.4567771 1.0000000
out$R
#> [,1] [,2] [,3] [,4]
#> [1,] 1.0000000 0.5522684 0.4472342 0.5817297
#> [2,] 0.5522684 1.0000000 0.4054908 0.5803080
#> [3,] 0.4472342 0.4054908 1.0000000 0.4563203
#> [4,] 0.5817297 0.5803080 0.4563203 1.0000000