A New Non-parametric Test for Multivariate Paired Data and Pair Matching

Let $N=2n$ be the total number of observations and let $Z_{i}=I(i\leq n)X_{i}+I(i>n)Y_{i-n}$, $i=1,\ldots,N$, where $I(\cdot)$ is the indicator function. Let $g_{i}$ be the indicator function that equals 1 when $Z_{i}$ is assigned to sample 1 under the paired-comparison permutation null distribution, and 0 if $Z_{i}$ is assigned to sample 2. It is easy to see that $\text{pr}(g_{i}=1)=0.5$. We use $a\wedge b$ to denote the minimum of $a$ and $b$, and use $a\vee b$ to denote the maximum of $a$ and $b$. Let $i$ and $i^{*}$ be the two indices belonging to pair $\mathbf{i}=(i\wedge i^{*},i\vee i^{*})$, i.e., $i^{*}=i+n$ if $i\leq n$, and $i^{*}=i-n$ if $i>n$. Then, we always have $g_{i}+g_{i^{*}}=1$. Let $j$ and $j^{*}$ be the two indices belonging to pair $\mathbf{j}=(j\wedge j^{*},j\vee j^{*})$. Since assigning $Z_{i}$ to sample 1 is independent of assigning $Z_{j}$ $(\mathbf{j}\neq\mathbf{i})$ to sample 1, $g_{i}$ and $g_{j}$ are independent.

For an edge in $G$, we denote it by the indices of the nodes connected by the edge. Then by definition, we have

Here, we do not distinguish edge $(i,j)$ from edge $(j,i)$. For graph $G$, let $|G|$ be the number of edges in the graph. Let $G_{1}$ be the subgraph in $G$ that connects observations from different pairs, i.e., $G_{1}$ contains edges $\{(i,j)\in G:j\neq i^{*}\}$. Let $G_{1,i}$ be the subgraph in $G_{1}$ that connects to node $i$. Then $|G_{1,i}|$ is the degree of node $i$ in $G_{1}$. Let $C_{1}$ be the number of pairs of edges $(i,j),~{}(i^{*},j^{*})\in G_{1}$, and $C_{2}$ be the number of pairs of edges $(i,j),~{}(i,j^{*})\in G_{1}$. The analytic expressions are provided in the following theorem.

To ensure that the proposed test statistic $D$ is well defined, $\Sigma_{R}$ needs to be invertible.

For the critical value $C(\alpha)$, it can be determined by performing the paired-comparison permutation directly, which is however time consuming. To make the test more application friendly, we study the asymptotic distribution of the statistic $D$.

Before stating the theorem, we define two additional terms on the similarity graph $G_{1}$: for edge $e\in G_{1}$, let $e_{-}$ and $e_{+}$ be the indices of the nodes connected by the edge $e$.

We use $a=O(b)$ to denote that $a$ and $b$ are of the same order, and $a=o(b)$ to denote that $a$ is of a smaller order than $b$.

To derive the asymptotic behavior of our statistic, we work under the following conditions for some $\gamma>0$:

Based on Theorem 3, it is easy to obtain the asymptotic distribution of $D$.

We reject the null hypothesis at $\alpha$ level when $D>C(\alpha)$. Based on Corollary 1, $C(\alpha)$ can be approximated by $\chi_{2}^{2}(1-\alpha)$, the $(1-\alpha)$ quantile of the $\chi_{2}^{2}$ distribution.

We next check how well the rejection region $D>\chi_{2}^{2}(1-\alpha)$ can control Type I error through numerical studies. We consider the following three settings:

• •

  Setting 1 (S1): $(X_{1},Y_{1})^{T},(X_{2},Y_{2})^{T},\ldots,(X_{n},Y_{n})^{T}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}_{2d}(\nu,\Gamma)$;

• •

  Setting 2 (S2): $(X_{1},Y_{1})^{T},(X_{2},Y_{2})^{T},\ldots,(X_{n},Y_{n})^{T}\stackrel{{\scriptstyle iid}}{{\sim}}\text{multivariate }t_{3}(\nu,\Gamma)$;

• •

  Setting 3 (S3): $\log(X_{1},Y_{1})^{T},\log(X_{2},Y_{2})^{T},\ldots,\log(X_{n},Y_{n})^{T}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}_{2d}(\nu,\Gamma)$.

Here, $\nu^{T}=(\nu_{1},\nu_{2})^{T}$ and $\Gamma=\begin{pmatrix}\Gamma_{1}&\Gamma_{12}\\ \Gamma_{12}&\Gamma_{2}\end{pmatrix}$ with $\nu_{1}=\nu_{2}=\mathbf{0}_{d}$, $\Gamma_{1}=\Gamma_{2}=\mathbf{I}_{d}$ and $\Gamma_{12}=0.6\mathbf{I}_{d}$. Then, in each setting, the distributions of $(X_{i},Y_{i})^{T}$ and $(Y_{i},X_{i})^{T}$ are the same.

Table 2 shows the empirical sizes of the proposed test under 0.05 and 0.1 nominal levels. The empirical size is computed as the fraction of trials (out of 1000) that the test statistic is greater than $\chi_{2}^{2}(1-\alpha)$. We see that the empirical sizes are well controlled for the proposed test even for quite small sample sizes under both low and high dimensions in all three settings.

We first examine the performance of the proposed test for data from the same family of distributions. We consider the three settings (S1), (S2) and (S3) in Section 2.2 but with different choices of $\nu$ and $\Gamma$. Let the number of pairs $n$ be fixed at the moderate size with $n=60$. We assess how the proposed statistic behaves when the dimension is comparable to or larger than the number of pairs by considering three dimensions: $d=50$, $d=100$ and $d=1000$. For each setting, we consider two alternatives for each $d$ listed below.

1. (i)

   Only $\nu_{1}$ differs from $\nu_{2}$ with $\nu_{1}=\mathbf{0}_{d}$, $\nu_{2}=0.5d^{-1/4}\mathbf{1}_{d}$, $\Gamma_{1}=\Gamma_{2}=\mathbf{I}_{d}$ and $\Gamma_{12}=0.6\mathbf{I}_{d}$.

2. (ii)

   Parameter $\nu_{1}$ differs from $\nu_{2}$ and $\Gamma_{1}$ differs from $\Gamma_{2}$ with $\nu_{1}=\mathbf{0}_{d}$, $\nu_{2}=0.5d^{-1/4}\mathbf{1}_{d}$, $\Gamma_{1}=\mathbf{I}_{d}$, $\Gamma_{2}=c_{d}\mathbf{I}_{d}$ and $\Gamma_{12}=0.6c_{d}^{1/2}\mathbf{I}_{d}$, where $c_{d}=1.15$ for $d=50$, $c_{d}=1.1$ for $d=100$ and $c_{d}=1.05$ for $d=1000$.

The results under various scenarios are summarized in Table 3. We first look at the results under moderate-dimensional settings ($d=50$). When only $\nu_{1}$ differs from $\nu_{2}$, the paired Hotelling’s $T^{2}$ test has high power for the multivariate normal and multivariate $t$ distributions, but is not good for the multivariate log-normal distribution. When $\nu_{1}\neq\nu_{2}$ and $\Gamma_{1}\neq\Gamma_{2}$, the performance of paired Hotelling’s $T^{2}$ is much worse than the proposed test $D$, especially for the multivariate log-normal distribution. For results under dimensions $d=100$ and $d=1000$, $D$ exhibits pretty good performance under all the settings, while the paired Hotelling’s $T^{2}$ cannot be applied as $d>n$.

Additionally, under setting (ii), $\nu_{1}$ and $\nu_{2}$ are set to be the same as those under setting (i), while $\Gamma_{2}$ becomes different from $\Gamma_{1}$. Thus, the difference between the two samples should be more significant. Comparing the results under setting (ii) with those under setting (i), we observe an increasing power of $D$ for the three distributions, while the power of pHT decreases for the multivariate normal and multivariate $t$ distributions.

Now we examine the performance of $D$ when $X_{i}$ and $Y_{i}$ are from different families of distributions. Suppose $X_{i},Y_{i}\in\mathbb{R}^{d}$, $X_{1},\ldots,X_{n}\stackrel{{\scriptstyle iid}}{{\sim}}F_{X}$, $Y_{1},\ldots,Y_{n}\stackrel{{\scriptstyle iid}}{{\sim}}F_{Y}$. We generate data by

where $\Delta$ is a constant vector, and $\alpha_{i},\epsilon_{i},\tau_{i}$ are independent with each other. Let $\alpha_{i}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}_{d}(\mathbf{0}_{d},\Omega_{1})$ and $\epsilon_{i}\stackrel{{\scriptstyle iid}}{{\sim}}t_{3}(\mathbf{0}_{d},\Omega_{2}/3)$. We consider three types of distributions for $\tau_{i}$: (i) $\tau_{i}\stackrel{{\scriptstyle iid}}{{\sim}}\mathcal{N}_{d}(\mathbf{0}_{d},\Omega_{2})$, (ii) $\tau_{i}\stackrel{{\scriptstyle iid}}{{\sim}}\text{multivariate skew normal distribution with mean }\mathbf{0}_{d}$, variance $\Omega_{2}$, skewness 1 and (iii) $\tau_{i}\stackrel{{\scriptstyle iid}}{{\sim}}\text{multivariate Laplace distribution with mean }\mathbf{0}_{d}$, variance $\Omega_{2}$. Therefore, $E(X_{i})=\mathbf{0}_{d}$, $E(Y_{i})=\Delta$, $\text{var}(X_{i})=\text{var}(Y_{i})=\Omega_{1}+\Omega_{2}$ and $\text{cov}(X_{i},Y_{i})=\Omega_{1}$. Let $n=60$ and $d=50,100,1000$. Let $\Omega_{1}=(\Omega_{ij})$ with the $(i,j)$ element of $\Omega$ being $\Omega_{ij}=0.5^{|i-j|}$ and $\Omega_{2}=\mathbf{I}_{d}$. For each $d$, we consider two settings for $\Delta$: $\Delta=\mathbf{0}_{d}$ and $\Delta=0.2\mathbf{1}_{d}$.

Table 3 gives the results. We first take a look at the results when there is no mean difference between the two samples ($\Delta=\mathbf{0}_{d}$). When $\Delta=\mathbf{0}_{d}$, the means and variances of $X_{i}$ and $Y_{i}$ are the same, while the two distributions $F_{X}$ and $F_{Y}$ are different, and we would expect a powerful test to reject the null hypothesis. However, the paired Hotelling’s $T^{2}$ test has almost no power under $d=50$ and cannot be applied to high-dimensional scenarios ($d=100,1000$), so it cannot detect the shape difference between two distributions. The proposed test $D$ performs well under all these scenarios. When there is a mean difference between the two samples ($\Delta=0.2\mathbf{1}_{d}$), the paired Hotelling’s $T^{2}$ exhibits some power under $d=50$, though the proposed test $D$ has a higher power. For the high dimensional scenarios ($d=100,1000$), the paired Hotelling’s $T^{2}$ is inapplicable, while the power of $D$ is very high.

Following a similar simulation setting as in Franklin et al. (2014), we assume that 20-dimensional covariates $\mathbf{X}_{(1)}=(X_{1},\cdots,X_{20})^{T}$ are observed, where $X_{j}~{}(j=1,\cdots,20)$ independent identically follows Laplace distribution with mean 0 and variance 0.65. The exposure/treatment, $T$, depends on these 20 covariates as well as three unobserved transformations of them ($\mathbf{X}_{(2)}=(X_{21},X_{22},X_{23})^{T}=(\sin(X_{1})/4,\cos(X_{2})/4,X_{3}^{4})^{T}$). We simulate $T$ as a binary variable via the logistic model $\text{logit}\{P(T=1)\}=\alpha_{0}+\bm{\alpha}_{1}^{T}\mathbf{X}_{(1)}+\bm{\alpha}_{2}^{T}\mathbf{X}_{(2)}$, where $\bm{\alpha}_{1}$ and $\bm{\alpha}_{2}$ define the log-odds ratios between exposures and controls in the pre-matched data set. We generated 1000 subjects and determined whether each of them is exposed or not. Then the control subjects $(T=0)$ are matched to exposed subjects $(T=1)$ through their propensity scores on $\mathbf{X}_{(1)}$.

We consider the following four scenarios.

1. (i)

   Zero coefficients for both $\mathbf{X}_{(1)}$ and $\mathbf{X}_{(2)}$: $\bm{\alpha}_{1}=\mathbf{0}_{20}$ and $\bm{\alpha}_{2}=\mathbf{0}_{3}$.

2. (ii)

   Nonzero coefficients for observed covariates $\mathbf{X}_{(1)}$ only: $\bm{\alpha}_{1}=0.2\mathbf{1}_{20}$ and $\bm{\alpha}_{2}=\mathbf{0}_{3}$.

3. (iii)

   Nonzero coefficients for unobserved covariates $\mathbf{X}_{(2)}$ only: $\bm{\alpha}_{1}=\mathbf{0}_{20}$ and $\bm{\alpha}_{2}=0.33\mathbf{1}_{3}$.

4. (iv)

   Nonzero coefficients for both $\mathbf{X}_{(1)}$ and $\mathbf{X}_{(2)}$: $\bm{\alpha}_{1}=0.2\mathbf{1}_{20}$ and $\bm{\alpha}_{2}=0.33\mathbf{1}_{3}$.

Here, scenario (i) is used to examine the empirical size since the covariates of treated subjects and controls are generated from the same distribution. For each scenario, we simulate 1000 data sets. Define the standardized mean difference as $\text{SD1}=(\bar{x}_{1}-\bar{x}_{0})/\sqrt{(s_{1}^{2}+s_{0}^{2})/2}$ where $\bar{x}_{m}$ and $s_{m}$ are the sample mean and variance for treated subjects $(m=1)$ and controls $(m=0)$. For each of the 23 covariates, we calculate the standardized mean difference between the treatment group and control group before and after matching. Figure 2 shows the boxplots of standardized mean difference of 1000 data set with the left panels denoting before matching and right panels denoting after matching. We see that the standardized mean differences of the observed covariates $\mathbf{X}_{(1)}$ are relatively close to 0 after matching under all scenarios. However, the unobserved covariates $\mathbf{X}_{(2)}$ are very unbalanced after matching when $\bm{\alpha}_{2}\neq\mathbf{0}_{3}$ (scenarios (iii) and (iv)), which indicates the distribution of ($X_{1},X_{2},X_{3}$) is not well balanced. Since we aim to test whether the joint distributions of the covariates in the matched control group and the treatment group are the same, we would expect that a good test rejects the null hypothesis under scenarios (iii) and (iv), while it accepts the null hypothesis under scenarios (i) and (ii).

Denote the method of combined differences and crossmatch test by CD and CM, respectively. To apply the CrossNN and CrossMST tests, we consider 1-MST and 5-MST as the similarity graph and denote these tests by $\text{NN}1$, $\text{NN}2$, $\text{MST}1$ and $\text{MST}2$, respectively. We present the proportion of trials that the tests reject the null hypothesis at 0.05 significance level in Table 4. We see all tests control empirical size well under scenario (i) and they deem the covariates well balanced under the reasonably balanced scenario (ii). For the results under scenarios (iii) and (iv), we observe the proposed test $D$ performs best with the highest power as expected. The crossmatch, CrossNN and CrossMST tests also show some power, while it is much less effective than the proposed test $D$. The paired Hotelling’s test and the method of combined differences do not have power under scenarios (iii) and (iv).

We first study the participants’ neuropsychologic performance over time, so we focused on the Neuropsychological Battery Summary Scores in Form C1, which contains 22 neuropsychology measurement variables. A description of these 22 variables are provided in LABEL:app:variable. We group the participants into three groups according to the CDR Dementia Staging Instrument in their first visits: no dementia (Group I), very mild dementia (Group II), and mild dementia (Group III). We test whether their neuropsychologic performances in five years differ from those in their initial visits. After removing missing data, the sample sizes ($n$) of three groups are 1747, 539, and 41, respectively.

Table 5 presents the results of the paired Hotelling’s $T^{2}$ test (pHT) and the proposed test $D$. We first check the results for Group I in testing whether their neuropsychologic measures in five years differ from those in their initial visits. We see that both the paired Hotelling’s $T^{2}$ test and our new test $D$ reject the null hypothesis with extremely small $p$-value. Thus, it is clear there is a significant difference in the neuropsychologic measures between the two visits. Similarly, we would reject the null hypothesis for Group II.

For Group III, as the sample size is relatively small, the power of the tests decreases. We see that the paired Hotelling’s $T^{2}$ test cannot reject the null hypothesis, while the proposed test $D$ rejects the null hypothesis with a small $p$-value. To see which result is more reliable, we check the data in more details. In particular, we perform the paired $t$-test to each of the 22 neuropsychologic measurement variables. It turns out that the paired $t$-test has $p$-value less than 0.05 for 11 out of these 22 variables and the smallest $p$-value is 0.001 (Table 1). Then, by the Bonferroni correction, the paired $t$-test would reject the null hypothesis at 0.05 significance level ($0.05/22=0.002>0.001$), supporting the result from $D$.

These results reveal that for all the participants from the three groups, neuropsychologic measures after several years would become significantly different from those in the initial visits. It suggests that researchers pay attention to the changes of cognitive performance even for no dementia or (very) mild dementia participants.

Several papers have previously reported a sex effect on Alzheimer’s disease (Payami et al., 1996; Podcasy and Epperson, 2016; Grimm et al., 2016). In this study, we focus on the data from the initial visit and only consider the white participant who has no stroke, no transient ischemic attack, no serious heart problem (e.g., atrial fibrillation, cardiac bypass procedure, congestive heart failure), no diabetes, no brain trauma, no Parkinson’s disease, no seizures, no other neurological and psychiatric disorder, do not abuse substances (except alcohol), and meanwhile, no Frontotemporal lobar degeneration mutation and Alzheimer’s disease mutation in her/his family, no cognitive impairment in her/his first-degree family. There are 13 covariates in total under consideration (Table 6). Here, each of TOBAC30, TOBAC100, CVHATT and DEP2YRS is a 0-1 variable, which is transformed to one dummy variable. CVANGIO and ALCOHOL are categorical variables with three categories: absent, recent/active, remote/inactive, and both of them are transformed to two dummy variables, respectively. We match 105 male participants with 105 female participants, among 153 female participants, expecting the covariates in the control (female) group are balanced with the exposed/treated (male) group. We adopt pairmatch() function in the optmatch package in R (Hansen, 2007) on treated to control distances with the distance computed from the match_on() using the logit propensity scores.

To measure the balance of covariates in the control group and the treatment group, we use the standardized mean difference (SD1), and define the standardized variance difference (SD2) and standardized third central moment difference (SD3) as

where $s_{m}$ and $\nu_{m}$ are the sample variance and the third sample central moment for treated subjects $(m=1)$ and controls $(m=0)$. Table 6 lists means, SD1, SD2 and SD3 before and after matching. We see that the covariates are not balanced before matching. For example, the values of SD1, SD2 and SD3 for CVHATT are 0.422, 1.534 and 4.741. After matching, it seems that the covariates are still not balanced well. For example, SD3 of the four covariates, CVHATT, CVANGIO1, CVANGIO2 and ALCOHOL1, are larger than 4.

Now we apply the paired Hotelling’s $T^{2}$ (pHT), the method of combined differences (CD), the crossmatch (CM), the CrossNN (NN1, NN2), the CrossMST (MST1, MST2) and the proposed test ($D$) to the paired data after matching. Only our proposed test gives small $p$-value 0.005. By contrast, the $p$-values of pHT, CD, CM, NN1, NN2, MST1 and MST2 are 0.087, 0.100, 0.669, 0.541, 0.219, 0.249, 0.266, respectively, which indicates that the other tests cannot reject the null hypothesis of balanced covariates at 0.05 significance level. To explore which result is more reliable, we apply the paired $t$-test to each of the 13 covariates. It turns out that the smallest $p$-value of the paired $t$-test is 0.001. So the paired $t$-test would reject the null hypothesis at 0.05 significance level ($0.05/13=0.004>0.001$) by the Bonferroni correction, supporting the result of $D$.

Another question of interest is to compare the neuropsychologic performances between well-matched female participants and male participants. To make the matching easier, we focus on the data from the initial visit and only consider the white participant who has no stroke, no transient ischemic attack, no atrial fibrillation, no heart problem, no angioplasty/endarterectomy/stent, no active depression in the last two years, no diabetes, no brain trauma, no Parkinson’s disease, no seizures, no alcohol abuse, no other neurological and psychiatric disorder, and meanwhile, no Frontotemporal lobar degeneration mutation and Alzheimer’s disease mutation in her/his family, no cognitive impairment in her/his first-degree family. In addition, only four covariates, NACCAGE, NACCBMI, SMOKYRS and TOBAC30, are under consideration now (Table 6). Adopting the same matching method as in Section 4.2, we match 76 male participants with 76 female participants, among 120 female participants.

The means, SD1, SD2 and SD3 before and after matching are provided in Table 6. We see the covariates of the data after matching are much more balanced than those in Section 4.2. We apply the proposed test $D$ and other existing tests to the matched covariates. Since the $p$-values of pHT, CD, CM, NN1, NN2, MST1, MST2 and $D$ are 0.372, 0.366, 0.804, 0.280, 0.688, 0.082, 0.163 and 0.056, respectively, all these tests cannot reject the null hypothesis at 0.05 significance level. So we deem the covariates reasonably balanced jointly.

After obtaining well-matched participants, we apply the paired Hotelling’s $T^{2}$ test (pHT) and the proposed test ($D$) to test whether the neuropsychologic performances in the female group differ from those in the male group. Here, 22 neuropsychology measurement variables are considered as in Section 4.1. The $p$-values of pHT and $D$ are 0.068 and 0.008, respectively, so the proposed test rejects the null hypothesis of equal neuropsychologic performances between the female group and the male group, while pHT accepts the null hypothesis at 0.05 significance level. We further explore the data by applying the paired $t$-test to each of the 22 covariates. We see that the smallest $p$-value of the paired $t$-test is 0.001. The paired $t$-test would reject the null hypothesis at 0.05 significance level ($0.05/22=0.002>0.001$) by the Bonferroni correction. Therefore, the proposed test $D$ is more reliable and the sex may have an effect on Alzheimer’s disease.