Back to the main page.


This example illustrates the usage of the proposed multivariate permutation approach based on the discrete optimal transport for a \(d\)-dimensional test statistic \(\mathbf{T}\) for \(d=3\). For simplicity, we use the same three sample problem as in the first example, but now we will consider all three pairwise comparisons.

I. Data generation and computation of \(\mathbf{T}_0\)

We use the same setup as in the first example, i.e. observations \(X_{j,i}\sim N(\mu_j,1)\), \(i=1,\dots,n_j\) for samples \(j\in\{1,2,3\}\) and for \(n_1=n_2=n_3=5\) and \(\mu_1=\mu_2=0\), \(\mu_3=2\). The data are generated in the same way with the same seed. The obtained values for the three samples are plotted below:

The test of the null hypothesis \[ \widetilde{H}_0:\mu_1=\mu_2=\mu_3 \] will be based on a three-variate test statistic \(\mathbf{T}=(T_1,T_2,T_3)^\top\), where \(T_1\), \(T_2\) and \(T_3\) are the classical t-statistics for testing \[ H_{0,1}:\mu_1=\mu_2, \ H_{0,2}:\mu_1=\mu_3 \text{ and } H_{0,3}:\mu_2=\mu_3, \] respectively.

First, we compute the test statistics from the data to get the vector \(\mathbf{T}_0\). In the code below n1=5 is the sample size of each of the samples.

T1.0=t.test(x[1:n1],x[(n1+1):(2*n1)])$stat
T2.0=t.test(x[1:n1],x[(2*n1+1):(3*n1)])$stat
T3.0=t.test(x[(n1+1):(2*n1)],x[(2*n1+1):(3*n1)])$stat
(T.0=c(T1.0,T2.0,T3.0))
##          t          t          t 
##  0.3750625 -4.5776556 -3.9617182

II. Generation of a grid

This part provides code for construction of both types of grids considered in the paper. Part A. below shows how to obtain a product grid grid, while section B. computes a non-product grid.

A. Product grid

Recall that a product grid \(G_{B+1}\) is a set of \(B+1=n_R\cdot n_S\) points is \(\mathbf{g}_{i,j}\) of the form \[ \mathbf{g}_{ij}=\frac{i}{n_R+1}\cdot \mathbf{s}_j, \quad i=1,\dots,n_R, \ j=1,\dots,n_S, \] where \(\mathbf{s}_j=(s_{j,1},s_{j,2},s_{j,3})^\top\) are vectors from the unit sphere obtained as a suitable transformation of low discrepancy points \((y_{j,1},y_{j,2})^\top\), \(j=1,\dots,n_S\), in \([0,1]^2\). In particular, it is shown in Chapter 1 of (Fang and Wang 1994) that if \(\mathbf{Y}= (Y_1,Y_2)^\top\) has a uniform distribution on \([0,1]^2\), then \(\mathbf{S}=(S_1,S_2,S_3)^\top\) defined as \[\begin{align*} S_1&=1-2 Y_1,\\ S_2&=2 \sqrt{Y_1(1-Y_1)}\cos(2\pi Y_2), \\ S_3& = 2\sqrt{Y_1(1-Y_1)}\sin(2\pi Y_2) \end{align*}\] has a uniform distribution on the unit sphere in \(\mathbb{R}^3\).

In this example we take \(B=799\), so \(B+1=800=n_R n_S\) for \(n_R=20\) and \(n_S=40\). The \(n_S\) vectors \((y_{i,1},y_{i,2})^\top\) from \([0,1]^2\) are taken as a Halton sequence, computed by the function halton from package randtoolbox.

n.R=20
n.S=40
d=3
N=n.R*n.S

library(randtoolbox)
y=halton(n.S,dim=2) # Halton sequence of n.S points

# transformation to the unit sphere  
s=matrix(NA,nrow=n.S,ncol=3)


s[,1]=(1-2*y[,1])
s[,2]=2*sqrt(y[,1]*(1-y[,1]))*cos(2*pi*y[,2])
s[,3]=2*sqrt(y[,1]*(1-y[,1]))*sin(2*pi*y[,2])
  

# final product grid
  
grid<-matrix(rep(0,d*N),ncol=d)
  radius=(1:n.R)/(n.R+1)
  
  for(i in 1:n.R) 
    for(j in 1:n.S)
      grid[(i-1)*n.S+j,]=radius[i]*s[j,]

The resulting grid is plotted in the following interactive plot:

Alternatively, one can use a Sobol sequence of points in \([0,1]^2\) (function sobol from randtoolbox) or a GLP sequence in \([0,1]^2\) (see the part B below).

B. Non-product grid

A non-product grid of \(B+1\) points is \(G_{B+1}=\{\mathbf{g}_i\}_{i=1}^{B+1}\), where \[ \mathbf{g}_i = y_{i,1} \mathbf{s}_i, \] where \((y_{i,1},y_{i,2},y_{i,3})^\top\), \(i=1,\dots,B+1\), form a low discrepancy sequence in \([0,1]^3\) and \((y_{i,2},y_{i,3})^\top\) are transformed to vectors \(\mathbf{s}_i=(s_{i,1},s_{i,2},s_{i,3})^\top\) from the unit sphere in \(\mathbb{R}^3\) using the same transformation as in part A.

In this example we illustrate the use of a GLP sequence for generating the low discrepancy sequence in \([0,1]^3\). Alternatively, one can use a Halton sequence (see part A above) or a Sobol sequence.

We will first generate a set of \(B+1=1010\) GLP points in \([0,1]^3\) using a generating vector \((h_1,h_2,h_3)^\top=c(1,140,237)^\top\). The costruction is based on Chapter 1 of (Fang and Wang 1994).

N=1010;h1=1;h2=140;h3=237 # generating vector of glp set
fract=function(x){x-floor(x)}

k=1:N
glp=matrix(NA,ncol=3,nrow=N)


glp[,1]=fract((2*k*h1-1)/(2*N))
glp[,2]=fract((2*k*h2-1)/(2*N))
glp[,3]=fract((2*k*h3-1)/(2*N))

The GLP points from \([0,1]^3\) are visualized in the plot below:

As a next step, the points are transformed into the non-product grid \(G_{B+1}\):

# transformation to s
s=matrix(NA,ncol=3,nrow=N)
y=glp[,2:3] # the first column is left as the radius
s[,1]=1-2*y[,1]
s[,2]=2*sqrt(y[,1]*(1-y[,1]))*cos(2*pi*y[,2])
s[,3]=2*sqrt(y[,1]*(1-y[,1]))*sin(2*pi*y[,2])
 
# resulting grid:
grid.GLP=glp[,1]*s

The resulting grid is shown in the interactive plot below.

The red points correspond to \(\mathbf{g}_b\) with the most extreme norm, where by the most extreme we mean \[ \|\mathbf{g}_b\|\geq Q_{0.05}, \] where \(Q_{0.05}\) is such that \[ \frac{1}{B+1}\sum_{b=1}^{B+1} \mathrm{I}\Bigl[\|\mathbf{g}_b\|\geq Q_{0.05}\Bigr]\leq 0.05. \] In other words, the null hypothesis is rejected on level \(\alpha=0.05\) if and only if \(\mathbf{T}_0\) is transported to one of the red points.

III. Permutation test

In this step, the vectors \(\mathbf{T}_1,\dots,\mathbf{T}_B\) of permutation test statistics are computed. The procedure is the same as for \(d=2\).

set.seed(1234)
B=n.R*n.S-1

#perm. test statistics T_b (in rows):
T.b=matrix(NA,ncol=3,nrow=B) 


for(b in 1:B){
  x.b=sample(x,size=length(x),replace=FALSE)
  T.b[b,1]=t.test(x.b[1:n1],x.b[(n1+1):(2*n1)])$stat
  T.b[b,2]=t.test(x.b[1:n1],x.b[(2*n1+1):(3*n1)])$stat
  T.b[b,3]=t.test(x.b[(n1+1):(2*n1)],x.b[(2*n1+1):(3*n1)])$stat
}

Dimension \(d=3\) still allows us to visualize the obtained permutation test statistics using a 3D plot. The blue points correspond to \(\mathbf{T}_1,\dots,\mathbf{T}_B\) and the red point is \(\mathbf{T}_0\). Obviously, it lies on the edge of the cloud. The next interactive plot enables rotation and so a good inspection of the whole situation.

IV. Computation of the transport and the final p-value

In this step we will compute the transportation \(F_{\pm}^{(B)}\) of \(\mathbf{T}_0,\mathbf{T}_1,\dots,\mathbf{T}_B\) to the target grid. We choose the product type grid from section II.A. for the transport. The optimization is performed using the function solve_LSAP from the package clue, similarly as for \(d=2\) in the first example.

require(clue)
## Loading required package: clue
sqdist<-function(a,b){
  return(abs(a-b)^2)
}

 d=3
 y=rbind(T.0,T.b)
 costes<-t(outer(grid[,1],y[,1],"sqdist"))
 
 for(j in 2:d) costes<-costes+t(outer(grid[,j],y[,j],"sqdist")) 
 asignacion<-solve_LSAP(costes)
 FB=grid[asignacion,]

Here, again, FB[1,] is \(F_{\pm}^{(B+1)}(\mathbf{T}_0)\) and FB[b+1,] corresponds to \(F_{\pm}^{(B+1)}(\mathbf{T}_b)\).

The value corresponding to \(F_{\pm}^{(B)}(\mathbf{T}_0)\) is stressed by a red color in the next interactive plot.

The red point lies on the most outer orbit of the grid. Since we chose \(n_R=20\) as in the example for \(d=2\), the same p-value (either \(\widehat{p}_e\) or \(\widehat{p}_a\)) is obtained.

vnorm=function(x) sqrt(sum(x^2))

(p.a=1-vnorm(FB[1,]))
## [1] 0.04761905
FB.norms=round(apply(FB,1,vnorm),dig=6)

(p.e=1/(B+1)*(1+sum(FB.norms[-1]>=FB.norms[1])))
## [1] 0.05

V. Contributions of the partial test statistics

To evaluate the partial contributions, we compute the directional vector \(\mathbf{D}\):

D=FB[1,]/vnorm(FB[1,])

# contributions in percentages
D^2*100
## [1]  3.515625 45.437154 51.047221

The contribution of the first, second and third test statistic to the rejection of \(H_0\) is 4 %, 45 % and 51 % respectively.

References

Fang, K. T., and Y. Wang. 1994. Number-Theoretic Methods in Statistics. Chapman & Hall.