# Asymptotic Theory in Model Diagnostic for General Multivariate Spatial Regression

Authors: Wayan Somayasa; Gusti N. Adhi Wibawa; La Hamimu; La Ode Ngkoimani
Journal: International Journal of Mathematics and Mathematical Sciences
 (2016)

---

## Abstract

We establish an asymptotic approach for checking the appropriateness of an assumed multivariate spatial regression model by considering the set-indexed partial sums process of the least squares residuals of the vector of observations. In this work, we assume that the components of the observation, whose mean is generated by a certain basis, are correlated. By this reason we need more effort in deriving the results. To get the limit process we apply the multivariate analog of the well-known Prohorov’s theorem. To test the hypothesis we define tests which are given by Kolmogorov-Smirnov (KS) and Cramér-von Mises (CvM) functionals of the partial sums processes. The calibration of the probability distribution of the tests is conducted by proposing bootstrap resampling technique based on the residuals. We studied the finite sample size performance of the KS and CvM tests by simulation. The application of the proposed test procedure to real data is also discussed.


---

## Body

## 1. Introduction

As mentioned in the literatures of model checks for multivariate regression, the appropriateness of an assumed model is mostly verified by analyzing the least squares residual of the observations; see, for example, Zellner [1], Christensen [2], pp. 1–22, Anderson [3], pp. 187–191, and Johnson and Wichern [4], pp. 395–398. A common feature of these works is the comparison between the length of the matrix of the residuals under the null hypothesis and that of the residuals under a proposed alternative.Instead of considering the residuals directly MacNeill [5] and MacNeill [6] proposed a method in model check for univariate polynomial regression based on the partial sums process of the residuals. These popular approaches are generalized to the spatial case by MacNeill and Jandhyala [7] for ordinary partial sums and Xie and MacNeill [8] for set-indexed partial sums process of the residuals. Bischoff and Somayasa [9] and Somayasa et al. [10] derived the limit process in the spatial case by a geometric method generalizing a univariate approach due to Bischoff [11] and Bischoff [12]. These results can be used to establish asymptotic test of Cramér-von Mises and Kolmogorov-Smirnov type for model checks and change-point problems. Model checks for univariate regression with random design using the empirical process of the explanatory variable marked by the residuals was established in Stute [13] and Stute et al. [14]. In the papers mentioned above the limit processes were explicitly expressed as complicated functions of the univariate Brownian motion (sheet).The purpose of the present article is to study the application of set-indexed partial sums technique to simultaneously check the goodness-of-fit of a multivariate spatial linear regression defined on high-dimensional compact rectangle. In contrast to the normal multivariate model studied in the standard literatures such as in Christensen [2], Anderson [3], and Johnson and Wichern [4] or in the references of model selection such as in Bedrick and Tsai [15] and Fujikoshi and Satoh [16], in this paper we will consider a multivariate regression model in which the components of the mean of the response vector are assumed to lie in different spaces and the underlying distribution model of the vector of random errors is unknown.To see the problem in more detail letn
1
≥
1
,
…
,
n
d
≥
1 be fixed. Let a p-dimensional random vector Z
≔
(
Z
i
)
i
=
1
p be observed independently over an experimental design given by a regular lattice: (1)
Ξ
n
1
⋯
n
d
≔
j
k
n
k
k
=
1
d
∈
I
d
:
1
≤
j
k
≤
n
k
,
k
=
1
,
…
,
d
,where I
d
≔
∏
k
=
1
d
0,1 is the d-dimensional unit cube. Let g
≔
(
g
i
)
i
=
1
p be the true but unknown R
p-valued regression function on I
d which represents the mean function of the observations. Let Z
j
1
⋯
j
d
≔
(
Z
i
,
j
1
⋯
j
d
)
i
=
1
p and g
≔
(
g
i
,
j
1
⋯
j
d
)
i
=
1
p be the observation and the corresponding mean in the experimental condition (
j
k
/
n
k
)
k
=
1
d. Under the null hypothesis H
0 we assume that Z
j
1
⋯
j
d follows a multivariate linear model. That is, we assume a model (2)
H
0
:
Z
j
1
⋯
j
d
=
∑
w
=
1
d
i
β
i
w
f
i
w
j
k
n
k
k
=
1
d
i
=
1
p
+
E
j
1
⋯
j
d
,
1
≤
j
1
≤
n
1
,
…
,
1
≤
j
d
≤
n
d
,where, for i
=
1
,
…
,
p, β
i
≔
(
β
i
w
)
w
=
1
d
i
∈
R
d
i is a d
i-dimensional vector of unknown parameters; f
i
≔
(
f
i
w
)
w
=
1
d
i is a d
i-dimensional vector of known regression functions whose components are assumed to be square integrable with respect to the Lebesgue measure λ
I
d on I
d, that is, f
i
w
∈
L
2
(
λ
I
d
), for all i and w. E
j
1
⋯
j
d
≔
(
ε
i
,
j
1
⋯
j
d
)
i
=
1
p is the mutually independent p-dimensional vector of random errors defined on a common probability space (
Ω
,
F
,
P
). We assume that, for all 1
≤
j
1
≤
n
1
,
…
,
1
≤
j
d
≤
n
d, E
(
E
j
1
⋯
j
d
)
=
0
∈
R
p, and C
o
v
(
E
j
1
⋯
j
d
)
=
Σ
=
(
σ
u
v
)
u
,
v
=
1
p
,
p. Let Z
n
1
⋯
n
d
≔
(
Z
j
1
⋯
j
d
)
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d be the p-dimensional pyramidal array of random observations and let E
n
1
⋯
n
d
≔
(
E
j
1
⋯
j
d
)
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d be the p-dimensional pyramidal array of random errors taking values in the Euclidean space ∏
i
=
1
p
R
n
1
×
⋯
×
n
d. Then under H
0 the observations can be represented by (3)
Z
n
1
⋯
n
d
=
∑
w
=
1
d
i
β
i
w
f
i
w
Ξ
n
1
⋯
n
d
i
=
1
p
+
E
n
1
⋯
n
d
,where f
i
w
(
Ξ
n
1
⋯
n
d
)
≔
f
i
w
(
j
k
/
n
k
)
k
=
1
d
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d. Under the alternative H
1 a multivariate nonparametric regression model (4)
Z
n
1
⋯
n
d
=
g
Ξ
n
1
⋯
n
d
+
E
n
1
⋯
n
dis assumed, where g
(
Ξ
n
1
⋯
n
d
)
≔
g
(
j
k
/
n
k
)
k
=
1
d
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d
∈
∏
i
=
1
p
R
n
1
×
⋯
×
n
d. By applying the similar argument as in Christensen [2] and Johnson and Wichern [4], the p-dimensional array of the least squares residuals of the observations is given by the following component-wise projection: (5)
R
n
1
⋯
n
d
≔
r
j
1
⋯
j
d
j
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d
≔
Z
n
1
⋯
n
d
-
p
r
∏
i
=
1
p
W
i
,
n
1
⋯
n
d
Z
n
1
⋯
n
d
,with r
j
1
⋯
j
d
≔
(
r
i
,
j
1
⋯
j
d
)
i
=
1
p, for 1
≤
j
k
≤
n
k, and k
=
1
,
…
,
d. Thereby, for i
=
1
,
…
,
p, we define W
i
,
n
1
⋯
n
d
≔
[
f
i
1
(
Ξ
n
1
⋯
n
d
)
,
…
,
f
i
d
i
(
Ξ
n
1
⋯
n
d
)
] as the subspace of R
n
1
×
⋯
×
n
d spanned by the arrays {
f
i
1
(
Ξ
n
1
⋯
n
d
)
,
…
,
f
i
d
i
(
Ξ
n
1
⋯
n
d
)
}. It is worth mentioning that the Euclidean space R
n
1
×
⋯
×
n
d is furnished with the inner product denoted by 〈
·
,
·
〉
R
n
1
×
⋯
×
n
d and defined by(6)
A
n
1
⋯
n
d
,
B
n
1
⋯
n
d
R
n
1
×
⋯
×
n
d
≔
∑
j
1
=
1
n
1
⋯
∑
j
d
=
1
n
d
a
j
1
⋯
j
d
b
j
1
⋯
j
d
,for every A
n
1
⋯
n
d
≔
(
a
j
1
⋯
j
d
)
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d, and B
n
1
⋯
n
d
≔
(
b
j
1
⋯
j
d
)
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d
∈
R
n
1
×
⋯
×
n
d.Next we define the set-indexed partial sums operator. LetA be the family of convex subset of I
d, and let η
λ
I
d be the Lebesgue pseudometric on A defined by η
λ
I
d
(
A
1
,
A
2
)
≔
λ
I
d
(
A
1
Δ
A
2
), for A
1
,
A
2
∈
A. Let C
(
A
) be the set of continuous functions on A under η
λ
I
d. We embed the array of the residual R
n
1
⋯
n
d into a p-dimensional stochastic process indexed by A by using the component-wise set-indexed partial sums operator (7)
V
n
1
⋯
n
d
:
∏
i
=
1
p
R
n
1
×
⋯
×
n
d
⟼
C
p
A
≔
∏
i
=
1
p
C
A
,such that, for any B
∈
A,(8)
V
n
1
⋯
n
d
R
n
1
⋯
n
d
B
≔
∑
j
1
=
1
n
1
⋯
∑
j
d
=
1
n
d
n
1
⋯
n
d
λ
I
d
B
∩
C
j
1
⋯
j
d
r
j
1
⋯
j
d
,where, for 1
≤
j
1
≤
n
1
,
…
,
1
≤
j
d
≤
n
d, C
j
1
⋯
j
d
≔
∏
k
=
1
d
(
(
j
k
-
1
)
/
n
k
,
j
k
/
n
k
]. Let us call this process the p-dimensional set-indexed least squares residual partial sums process. The space C
p
(
A
) is furnished with the uniform topology induced by the metric φ defined by (9)
φ
u
,
w
≔
∑
i
=
1
p
u
i
-
w
i
A
=
∑
i
=
1
p
sup
A
∈
A
⁡
u
i
A
-
w
i
A
, for u
≔
(
u
i
)
i
=
1
p and w
≔
(
w
i
)
i
=
1
p
∈
C
p
(
A
).We notice that, in the works of Bischoff and Somayasa [9], Bischoff and Gegg [17], and Somayasa and Adhi Wibawa [18], the limit process of the partial sums process of the least squares residuals has been investigated by applying the existing geometric method of Bischoff [11, 12]. However, the method becomes not applicable anymore in deriving the limit process of V
n
1
⋯
n
d
(
R
n
1
⋯
n
d
) as the dimension of W
i
,
n
1
⋯
n
d varies. Therefore, in this work, we attempt to adopt the vectorial analog of Prohorov’s theorem; see, for example, Theorem 6.1 in Billingsley [19] for obtaining the limit process. For our result we need to extend the ordinary partial sums formula to p-dimensional case defined on I
d as follows. Let K
≔
{
1,2
,
…
,
d
} and C
k
K be the set of all k-combinations of the set K, with k
=
1
,
…
,
d. For a chosen value of k, we denote the ith element of C
k
K by a k-tuple (
i
(
k
1
)
,
i
(
k
2
)
,
…
,
i
(
k
k
)
), for 1
≤
i
≤
C
k
K, where C
k
K is the number of elements of C
k
K which is clearly given by (10)
C
k
K
=
d
d
-
1
d
-
2
⋯
d
-
k
+
1
k
k
-
1
k
-
2
⋯
1
.For example, let K
=
{
1,2
,
3
}. Then, for k
=
1, we denote the elements of C
1
K as 1
(
k
1
)
≔
1, 2
(
k
1
)
≔
2, and 3
(
k
1
)
≔
3. In a similar way, we denote the elements of C
2
K which consists of {
(
1,2
)
,
(
1,3
)
,
(
2,3
)
}, respectively, by (
1
(
k
1
)
,
1
(
k
2
)
)
≔
(
1,2
), (
2
(
k
1
)
,
2
(
k
2
)
)
≔
(
1,3
), and (
3
(
k
1
)
,
3
(
k
2
)
)
≔
(
2,3
). Finally the element of C
3
K
=
{
(
1,2
,
3
)
} is sufficiently written by (
1
(
k
1
)
,
1
(
k
2
)
,
1
(
k
3
)
). Hence the p-dimensional ordinary partial sums operator transforms any p-dimensional array A
n
1
⋯
n
d
=
(
a
j
1
⋯
j
d
)
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d
∈
∏
i
=
1
p
R
n
1
×
⋯
×
n
d to a continuous function on I
d defined by(11)
T
n
1
⋯
n
d
A
n
1
⋯
n
d
t
≔
1
n
1
⋯
n
d
∑
j
1
=
1
⋯
j
d
=
1
n
1
t
1
⋯
n
d
t
d
a
j
1
⋯
j
d
+
∑
k
=
1
d
∑
i
=
1
C
2
K
∏
u
=
1
k
t
i
k
u
-
n
i
k
u
t
i
k
u
n
i
k
u
n
1
⋯
n
d
n
1
⋯
n
d
/
n
i
k
1
⋯
n
i
k
d
∑
j
1
=
1
,
…
,
j
d
=
1
/
j
i
k
1
,
…
,
j
i
k
d
n
1
t
1
,
…
,
n
d
t
d
a
j
1
⋯
j
d
∣
j
i
k
1
=
n
i
k
1
t
i
k
1
+
1
,
…
,
j
i
k
d
=
n
i
k
d
t
i
k
d
+
1
,for every t
≔
(
t
1
,
…
,
t
d
)
⊤
∈
I
d, where for positive integers b
u
,
b
u
+
1
,
…
,
b
u
+
m
∈
Z
+ we define a notation (12)
a
j
1
⋯
j
d
∣
j
u
=
b
u
,
…
,
j
u
+
m
=
b
u
+
m
≔
a
j
1
,
…
,
j
u
-
1
,
b
u
,
…
,
b
u
+
m
,
j
u
+
m
+
1
,
…
,
j
d
.It is clear that the partial sums process of the residuals obtained using (11) is a special case of (8) since for every t
≔
(
t
1
,
…
,
t
d
)
⊤
∈
I
d it holds (13)
T
n
1
⋯
n
d
A
n
1
⋯
n
d
t
=
V
n
1
⋯
n
d
A
n
1
⋯
n
d
∏
k
=
1
d
0
,
t
k
.It is worth noting that the extension of the study from univariate to multivariate model and also the expansion of the dimension of the lattice points are strongly motivated by the prediction problem in mining industry and geosciences. As for an example recently Tahir [20] presented data provided by PT Antam Tbk (a mining industry in Southeast Sulawesi). The data consist of a joint measurement of the percentage of several chemical elements and substances such as Ni, Co, Fe, MgO, S
i
O
2, and CaO which are recorded in every point of a three-dimensional lattice defined over the exploration region of the company. Hence, by the inherent existence of the correlation among the variables, the statistical analysis for the involved variables must be conducted simultaneously.There have been many methods proposed in the literatures for testingH
0. Most of them have been derived for the case of W
1
=
⋯
=
W
p under normally distributed random error. Generalized likelihood ratio test which has been leading to Wilk’s lambda statistic or variant of it can be found in Zellner [1], Christensen [2], pp. 1–22, Anderson [3], pp. 187–191, and Johnson and Wichern [4], pp. 395–398. Mardia and Goodall [21] derived the maximum likelihood estimation procedure for the parameters of the general normal spatial multivariate model with stationary observations. This approach can be straightforwardly extended for obtaining the associated likelihood ratio test in model check for the model. Unfortunately, in the practice especially when dealing with mining data, normal distribution is sometimes found to be not suitable for describing the distribution model of the observations, so that the test procedures mentioned above are consequently no longer applicable. Asymptotic method established in Arnold [22] for multivariate regression with W
1
=
⋯
=
W
p can be generalized in such a way that it is valid for the general model. As a topic in statistics it must be well known. However, we cannot find literatures where the topic has been studied.The rest of the paper is organized as follows. In Section2 we show that when H
0 is true Σ
-
1
/
2
V
n
1
⋯
n
d
(
R
n
1
⋯
n
d
) converges weakly to a projection of the p-dimensional set-indexed Brownian sheet. The limit process is shown to be useful for testing H
0 asymptotically based on the Kolmogorov-Smirnov (KS-test) and Cramér-von Mises (CvM-test) functionals of the set-indexed p-dimensional least squares residual partial sum processes, defined, respectively, by (14)
KS
n
1
⋯
n
d
,
A
≔
sup
A
∈
A
⁡
Σ
-
1
/
2
V
n
1
⋯
n
d
R
n
1
⋯
n
d
A
R
p
CvM
n
1
⋯
n
d
,
A
≔
1
n
1
⋯
n
d
∑
A
∈
A
Σ
-
1
/
2
V
n
1
⋯
n
d
R
n
1
⋯
n
d
A
R
p
2
.For both tests the rejection of H
0 is for large value of the KS and CvM statistics, respectively. Under localized alternative the above sequence of random processes converges weakly to the above limit process with an additional deterministic trend (see Section 3). In Section 4, we define a consistent estimator for Σ. In Section 5 we investigate a residual based bootstrap method for the calibration of the tests. Monte Carlo simulation for the purpose of studying the finite sample behavior of the KS and CvM tests is reported in Section 6. Application of the test procedure in real data is presented in Section 7. The paper is closed in Section 8 with conclusion and some remarks for future research. Auxiliary results needed for deriving the limit process are presented in the appendix. We note that all convergence results derived throughout this paper which hold for n
1
,
…
,
n
d simultaneously go to infinity, that is, for n
k
→
∞, for all k
=
1
,
…
,
d; otherwise they will be stated in some way. The notion of convergence in distribution and convergence in probability will be conventionally denoted by →
D and →
P, respectively.


## 2. The Limit of Vn
1
⋯
n
d
(Rn
1
⋯
n
d
) under H
0

LetB
≔
{
B
(
A
)
:
A
∈
A
} be the one-dimensional set-indexed Brownian sheet having sample path in C
(
A
). We refer the reader to Pyke [23], Bass and Pyke [24], and Alexander and Pyke [25] for the definition and the existence of B. Let us consider a subspace of C
(
A
) which is closely related to B, defined by(15)
H
B
≔
h
v
:
A
⟶
R
,
∃
v
∈
L
2
λ
I
d
,
s.t.
h
v
A
≔
∫
A
v
t
λ
d
t
.Under the inner product and the norm defined by (16)
h
v
,
h
w
H
B
≔
∫
I
d
v
t
w
t
λ
I
d
t
,
h
v
2
≔
∫
I
d
v
2
t
λ
I
d
t
,it is clear that H
B and L
2
(
λ
I
d
) are isometric. For our result we need to define subspaces W
i and W
i
H
B associated with the regression functions f
i
1
,
…
,
f
i
d
i, where W
i
≔
[
f
i
1
,
…
,
f
i
d
i
]
⊂
L
2
(
λ
I
d
) and W
i
H
B
≔
[
h
f
i
1
,
…
,
h
f
i
d
i
]
⊂
H
B, for i
=
1
,
…
,
p.Now we are ready to state the limit process of the sequence ofp-dimensional set-indexed residual partial sums processes for the model specified under H
0.Theorem 1.
Fori
=
1
,
…
,
p, let {
f
i
1
,
…
,
f
i
d
i
} be an orthonormal basis (ONB) of W
i. We assume that W
i
⊆
W
i
+
1, for i
=
1
,
…
,
p
-
1. Let B
p
≔
{
(
B
i
(
A
)
)
i
=
1
p
:
A
∈
A
} be the p-dimensional set-indexed Brownian sheet with the covariance function C
o
v
(
B
p
(
A
)
,
B
p
(
A
′
)
)
≔
λ
I
d
(
A
∩
A
′
)
I
p, for A
,
A
′
∈
A, where I
p is the p
×
p identity matrix. Suppose {
f
i
1
,
…
,
f
i
d
i
} are in C
(
I
d
)
∩
B
V
H
(
I
d
), where C
(
I
d
) is the space of continuous functions on I
d (see Definition A.4 for the definition of B
V
H
(
I
d
)). Then under H
0 it holds that (17)
Σ
-
1
/
2
V
n
1
⋯
n
d
R
n
1
⋯
n
d
⟶
D
B
p
,
f
H
0
≔
B
p
-
p
r
∏
i
=
1
p
W
i
H
B
∗
B
p
,where (18)
p
r
∏
i
=
1
p
W
i
H
B
∗
B
p
≔
p
r
W
i
H
B
∗
B
i
i
=
1
p
.Thereby p
r
W
i
H
B
∗ is a projector such that, for every u
∈
C
(
A
) and A
∈
A, (19)
p
r
W
i
H
B
∗
u
A
≔
∑
j
=
1
d
i
h
f
i
j
,
u
h
f
i
j
A
,
w
h
e
r
e
h
f
i
j
,
u
≔
∫
I
d
R
f
i
j
t
d
u
t
.For t
≔
(
t
k
)
k
=
1
d
∈
I
d, we set u
(
t
) for u
(
∏
k
=
1
d
[
0
,
t
k
]
), and ∫
(
R
) stands for the integral in the sense of Riemann-Stieltjes. Moreover, the limit process B
p
,
f
H
0 is a centered Gaussian process with the covariance function given by (20)
K
A
,
C
≔
diag
⁡
λ
A
∩
C
-
∑
j
=
1
d
1
h
f
1
j
A
h
f
1
j
C
,
…
,
λ
A
∩
C
-
∑
j
=
1
d
p
h
f
p
j
A
h
f
p
j
C
.Proof.
By applying the linear property ofV
n
1
⋯
n
d and Lemma C.2, we have under H
0,(21)
V
n
1
⋯
n
d
R
n
1
⋯
n
d
=
V
n
1
⋯
n
d
E
n
1
⋯
n
d
-
p
r
∏
i
=
1
p
W
i
,
n
1
⋯
n
d
H
B
V
n
1
⋯
n
d
E
n
1
⋯
n
d
.It can be shown by extending the uniform central limit theorem studied in Pyke [23], Bass and Pyke [24], and Alexander and Pyke [25] to its vectorial analog that the term Σ
-
1
/
2
V
n
1
⋯
n
d
(
E
n
1
⋯
n
d
) on the right-hand side of (21) converges weakly to B
p. Therefore we only need to show that the second term satisfies the weak convergence: (22)
Z
n
1
⋯
n
d
≔
Σ
-
1
/
2
p
r
∏
i
=
1
p
W
i
,
n
1
⋯
n
d
H
B
V
n
1
⋯
n
d
E
n
1
⋯
n
d
⟶
D
U
≔
p
r
∏
i
=
1
p
W
i
H
B
∗
B
p
,where U is a p-dimensional centered Gaussian process with the covariance matrix given by (23)
K
U
A
1
,
A
2
=
diag
⁡
∑
j
=
1
d
1
h
f
1
j
A
1
h
f
1
j
A
2
,
…
,
∑
j
=
1
d
p
h
f
p
j
A
1
h
f
p
j
A
2
,
for
A
1
,
A
2
∈
A
.By Prohorov’s theorem it is sufficient to show that, for any finite collection of convex sets A
1
,
…
,
A
r in A and real numbers c
1
,
…
,
c
r, with r
≥
1, it holds that (24)
∑
k
=
1
r
c
k
Z
n
1
⋯
n
d
A
k
⟶
D
∑
k
=
1
r
c
k
U
A
k
, where the left-hand side has the covariance which can be expressed as (25)
Cov
∑
k
=
1
r
c
k
Z
n
1
⋯
n
d
A
k
=
∑
k
=
1
r
∑
l
=
1
r
c
k
c
l
Σ
-
1
/
2
E
B
i
A
k
B
j
A
l
i
,
j
=
1
p
Σ
-
1
/
2
,where B
i
(
A
k
)
≔
p
r
W
i
,
n
1
⋯
n
d
H
B
V
n
1
⋯
n
d
(
ε
i
,
n
1
⋯
n
d
)
(
A
k
), for i
=
1
,
…
,
p and k
=
1
,
…
,
r. Let k and l be fixed. Then by a simultaneous application of the definition of p
r
W
i
,
n
1
⋯
n
d
H
B (see Lemma C.3), (11), the definition of the Riemann-Stieljes integral (cf. Stroock [26], pp. 7–17), and the independence of the array {
ε
i
,
j
1
⋯
j
d
:
1
≤
j
k
≤
n
k
}, we further get (26)
E
B
i
A
k
B
j
A
l
=
∑
w
=
1
d
i
∑
w
′
=
1
d
j
∫
A
k
s
~
i
w
n
1
⋯
n
d
t
λ
I
d
d
t
∫
A
l
s
~
j
w
′
n
1
⋯
n
d
t
λ
I
d
d
t
E
∫
I
d
R
s
~
i
w
n
1
⋯
n
d
t
d
V
n
1
⋯
n
d
ε
i
,
n
1
⋯
n
d
t
∫
I
d
R
s
~
j
w
′
n
1
⋯
n
d
t
d
V
n
1
⋯
n
d
ε
j
,
n
1
⋯
n
d
t
=
∑
w
=
1
d
i
∑
w
′
=
1
d
j
∫
A
k
s
~
i
w
n
1
⋯
n
d
t
λ
I
d
d
t
∫
A
l
s
~
j
w
′
n
1
⋯
n
d
t
λ
I
d
d
t
σ
i
j
n
1
⋯
n
d
∑
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d
s
~
i
w
n
1
⋯
n
d
j
1
n
1
,
…
,
j
d
n
d
s
~
j
w
′
n
1
⋯
n
d
j
1
n
1
,
…
,
j
d
n
d
.By recalling Lemma C.3 the last expression clearly converges to (27)
σ
i
j
∑
w
=
1
d
i
∑
w
′
=
1
d
j
f
i
w
,
f
j
w
′
L
2
λ
I
d
h
f
i
w
A
k
h
f
j
w
′
A
l
=
σ
i
j
∑
w
=
1
d
i
h
f
i
w
A
k
h
f
j
w
A
l
.Hence the matrix E
(
B
i
(
A
k
)
B
j
(
A
l
)
)
i
,
j
=
1
p converges element wise to the symmetric matrix which can be represented as Σ
⊙
D, for a matrix D defined by (28)
D
≔
κ
11
A
k
,
A
l
κ
12
A
k
,
A
l
⋯
κ
1
p
A
k
,
A
l
κ
21
A
k
,
A
l
κ
22
A
k
,
A
l
⋯
κ
2
p
A
k
,
A
l
⋮
⋮
⋮
⋮
κ
p
1
A
k
,
A
l
κ
p
2
A
k
,
A
l
⋯
κ
p
p
A
k
,
A
lwith κ
i
j
(
A
k
,
A
l
)
≔
∑
w
=
1
d
min
⁡
{
i
,
j
}
h
f
i
w
(
A
k
)
h
f
j
w
(
A
l
), for i
,
j
=
1
,
…
,
p. Thereby ⊙ denotes the Hadarmard product defined, for example, in Magnus and Neudecker [27], pp. 53-54. Since Σ
-
1
/
2
(
Σ
⊙
D
)
Σ
-
1
/
2
=
D
⊙
I
p, we successfully have shown that C
o
v
∑
k
=
1
r
c
k
Z
n
1
⋯
n
d
A
k converges to the following general linear combination:(29)
∑
k
=
1
r
∑
l
=
1
r
c
k
c
l
diag
⁡
∑
w
=
1
d
1
h
f
1
w
A
k
h
f
1
w
A
l
,
…
,
∑
w
=
1
d
p
h
f
p
w
A
k
h
f
p
w
A
l
,which is actually the covariance of ∑
k
=
1
r
c
k
U
(
A
k
). Next we observe that, by applying the definition of V
n
1
⋯
n
d and the definition of the Riemann-Stieltjes integral, we can also write ∑
k
=
1
r
c
k
Z
n
1
⋯
n
d
(
A
k
) as follows: (30)
∑
k
=
1
r
c
k
Z
n
1
⋯
n
d
A
k
=
∑
j
1
=
1
n
1
⋯
∑
j
d
=
1
n
d
Σ
-
1
/
2
∑
k
=
1
r
Λ
i
j
1
n
1
,
…
,
j
d
n
d
;
A
k
i
=
1
p
,where (31)
Λ
i
j
1
n
1
,
…
,
j
d
n
d
;
A
k
≔
∑
w
=
1
d
i
c
k
n
1
⋯
n
d
s
~
i
w
n
1
⋯
n
d
j
1
n
1
,
…
,
j
d
n
d
ε
i
,
j
1
⋯
j
d
h
s
~
i
w
n
1
⋯
n
d
A
k
.Let Γ
f
≔
m
a
x
1
≤
w
≤
d
i
;
1
≤
i
≤
p
f
i
w
∞, M
≔
max
1
≤
k
≤
r
⁡
c
k, and Σ
-
1
/
2 be the Euclidean norm of Σ
-
1
/
2. Then by considering the stochastic independence of the array of the p-vector of the random errors, it holds that (32)
0
≤
L
Z
n
1
⋯
n
d
ϵ
≔
∑
j
1
=
1
n
1
⋯
∑
j
d
=
1
n
d
E
∑
k
=
1
r
Λ
i
j
1
n
1
,
…
,
j
d
n
d
;
A
k
i
=
1
p
R
p
2
1
∑
k
=
1
r
Λ
i
j
1
n
1
,
…
,
j
d
n
d
;
A
k
i
=
1
p
R
p
2
≥
ϵ
≤
r
M
2
d
p
Γ
f
2
2
Σ
-
1
2
2
E
E
j
1
⋯
j
d
R
p
2
1
E
j
1
⋯
j
d
R
p
≥
ϵ
n
1
⋯
n
d
r
M
d
p
Γ
f
2
∀
ϵ
>
0
,in which by the well-known bounded convergence theorem (cf. Athreya and Lahiri [28], pp. 57-58) the last term converges to zero. Thus the Lindeberg condition is satisfied. Therefore by the Lindeberg-Levy multivariate central limit theorem studied, for example, in van der Vaart [29], pp. 16, it can be concluded that ∑
k
=
1
r
c
k
Z
n
1
⋯
n
d
(
A
k
) converges in distribution to ∑
k
=
1
r
c
k
U
(
A
k
), where ∑
k
=
1
r
c
k
U
(
A
k
) has the p-variate normal distribution with mean zero and the covariance given by (29).
The tightness ofZ
n
1
⋯
n
d can be shown as follows. By the definition Z
n
1
⋯
n
d can also be expressed as (33)
Z
n
1
⋯
n
d
=
Σ
-
1
/
2
p
r
∏
i
=
1
p
W
i
,
n
1
⋯
n
d
H
B
Σ
1
/
2
Σ
-
1
/
2
V
n
1
⋯
n
d
E
n
1
⋯
n
d
.Since Σ
-
1
/
2
V
n
1
⋯
n
d
(
E
n
1
⋯
n
d
) is tight, hence by recalling Lemma 1 in Billingsley [19], pp. 38–40, it is sufficient to show that the mapping Σ
-
1
/
2
p
r
∏
i
=
1
p
W
i
,
n
1
⋯
n
d
H
B
Σ
1
/
2 is continuous on C
p
(
A
) for every n
1
≥
1
,
…
,
n
d
≥
1. Proposition C.4 finishes the proof.Corollary 2.
By Theorem1 and the well-known continuous mapping theorem (cf. Theorem 5.1 in Billingsley [19]) the distribution of the statistics K
S
n
1
⋯
n
d
,
A and C
v
M
n
1
⋯
n
d
,
A can be approximated, respectively, by those of (34)
K
S
≔
s
u
p
A
∈
A
B
p
-
p
r
∏
i
=
1
p
W
i
H
B
∗
B
p
A
R
p
C
v
M
≔
∫
I
d
B
p
-
p
r
∏
i
=
1
p
W
i
H
B
∗
B
p
A
R
p
2
d
A
.Let t
~
1
-
α and c
~
1
-
α be the (
1
-
α
)th quantile of the distributions of K
S and C
v
M, respectively. When K
S
n
1
⋯
n
d
,
A is used, H
0 will be rejected at level α if and only if K
S
n
1
⋯
n
d
,
A
≥
t
~
1
-
α. Likewise if C
v
M
n
1
⋯
n
d
,
A is used, then H
0 will be rejected at level α if and only if C
v
M
n
1
⋯
n
d
,
A
≥
c
~
1
-
α.


## 3. The Limit of Vn
1
⋯
n
d
(Rn
1
⋯
n
d
) under H
1

The test procedures derived above are consistent in the sense of Definition11.1
.
3 in Lehmann and Romano [30]. That is, the probability of rejection of H
0 under the competing alternative converges to 1. As an immediate consequence we cannot observe the performance of the tests when the model moves away from H
0. Therefore, to be able to investigate the behavior of the tests, we consider the localized model defined as follows:(35)
Z
n
1
⋯
n
d
=
1
n
1
⋯
n
d
g
Ξ
n
1
⋯
n
d
+
E
n
1
⋯
n
d
.When H
0 is true we get the similar least squares residuals as given in Section 1. Therefore, observing Model (35), the test problem will not be altered.In the following theorem, we present the limit process of thep-dimensional set-indexed partial sums process of the residuals under H
1 associated with Model (35).Theorem 3.
Fori
=
1
,
…
,
p, let {
f
i
1
,
…
,
f
i
d
i
} be an ONB of W
i with W
i
⊆
W
i
+
1, for i
=
1
,
…
,
p
-
1. If g
=
(
g
i
)
i
=
1
p
∈
∏
i
=
1
p
B
V
V
(
I
d
) and {
f
i
1
,
…
,
f
i
d
i
}
∈
C
(
I
d
)
∩
B
V
H
(
I
d
) (see Definition A.3 for the notion of B
V
V
(
I
d
)), then, observing (35), we have under H
1 that (36)
Σ
-
1
/
2
V
n
1
⋯
n
d
R
n
1
⋯
n
d
⟶
D
Σ
-
1
/
2
I
p
-
p
r
∏
i
=
1
p
W
i
H
B
∗
h
g
+
B
p
,
f
H
0
,where h
g
≔
(
h
g
i
)
i
=
1
p
:
A
→
C
p
(
A
), with h
g
i
(
A
)
≔
∫
A
g
i
(
t
)
λ
I
d
(
d
t
).Proof.
Considering the linearity ofV
n
1
⋯
n
d and Lemma C.2, when H
1 is true we have (37)
Σ
-
1
/
2
V
n
1
⋯
n
d
R
n
1
⋯
n
d
=
Σ
-
1
/
2
1
n
1
⋯
n
d
V
n
1
⋯
n
d
g
Ξ
n
1
⋯
n
d
-
Σ
-
1
/
2
p
r
∏
i
=
1
p
W
i
,
n
1
⋯
n
d
H
B
1
n
1
⋯
n
d
V
n
1
⋯
n
d
g
Ξ
n
1
⋯
n
d
+
Σ
-
1
/
2
V
n
1
⋯
n
d
E
n
1
⋯
n
d
-
Σ
-
1
/
2
p
r
∏
i
=
1
p
W
i
,
n
1
⋯
n
d
H
B
V
n
1
⋯
n
d
E
n
1
⋯
n
d
.Since, for i
=
1
,
…
,
p, g
i is in B
V
V
(
I
d
), it can be shown that (
1
/
n
1
⋯
n
d
)
V
n
1
⋯
n
d
(
g
(
Ξ
n
1
⋯
n
d
)
)
(
B
) converges uniformly to h
g
B, for every B
∈
A. Also the last two terms on the right-hand side of the preceding equation converge in distribution to B
p
,
f
H
0 by Theorem 1. Thus to the rest we only need to show that p
r
∏
i
=
1
p
W
i
,
n
1
⋯
n
d
H
B
(
1
/
n
1
⋯
n
d
)
V
n
1
⋯
n
d
(
g
n
1
⋯
n
d
) converges to p
r
∏
i
=
1
p
W
i
H
B
B
V
(
I
d
)
∗
h
g. By the definition of component-wise projection we have (38)
p
r
∏
i
=
1
p
W
i
,
n
1
⋯
n
d
H
B
1
n
1
⋯
n
d
V
n
1
⋯
n
d
g
Ξ
n
1
⋯
n
d
=
p
r
W
i
,
n
1
⋯
n
d
H
B
1
n
1
⋯
n
d
V
n
1
⋯
n
d
i
g
i
Ξ
n
1
⋯
n
d
i
=
1
p
=
∑
j
=
1
d
i
∫
I
d
R
s
~
i
j
n
1
⋯
n
d
t
d
1
n
1
⋯
n
d
V
n
1
⋯
n
d
i
g
i
Ξ
n
1
⋯
n
d
t
h
s
~
i
j
n
1
⋯
n
d
i
=
1
p
=
∑
j
=
1
d
i
1
n
1
⋯
n
d
∑
j
1
=
1
n
1
⋯
∑
j
d
=
1
n
d
s
~
i
j
n
1
⋯
n
d
j
1
n
1
,
…
,
j
d
n
d
g
i
j
1
n
1
,
…
,
j
d
n
d
h
s
~
i
j
n
1
⋯
n
d
i
=
1
p
.The right-hand side of the last expression is obtained directly from the definition of the ordinary partial sums (11) and the definition of the Riemann-Stieltjes integral on I
d. Since s
~
i
j
(
n
1
⋯
n
d
) converges uniformly to f
i
j and g
i has bounded variation on I
d for all i and j, then the last expression clearly converges component-wise to ∑
j
=
1
d
i
∫
I
d
(
R
)
f
i
j
(
t
)
d
h
g
i
(
t
)
h
f
i
j
i
=
1
p, which can be written as p
r
∏
i
=
1
p
W
i
H
B
∗
h
g. We are done.Corollary 4.
By Theorem3 the power function of the KS and CvM tests at a level α can now be approximated by computing the probabilities of the form (39)
P
s
u
p
A
∈
A
Σ
-
1
/
2
h
g
-
p
r
∏
i
=
1
p
W
i
H
B
B
V
I
d
∗
h
g
+
B
p
,
f
H
0
R
p
≥
t
~
1
-
α
P
∫
I
Σ
-
1
/
2
h
g
-
p
r
∏
i
=
1
p
W
i
H
B
B
V
I
d
∗
h
g
+
B
p
,
f
H
0
R
p
2
d
A
≥
c
~
1
-
α
,for a fixed g
∈
∏
i
=
1
p
B
V
V
(
I
d
), respectively. In Section 5 we investigate the empirical power functions of the KS and CvM tests by simulation.


## 4. Estimating the Population Covariance Matrix

If the covariance matrixΣ is unknown, as it usually is, it is impossible to use K
S
n
1
⋯
n
d
;
A and C
v
M
n
1
⋯
n
d
;
A in practice. What we propose to do is to employ a consistent estimate of Σ. We need some further notations for expressing the residuals of the model. For i
=
1
,
…
,
p, let Z
i
(
n
1
⋯
n
d
), g
i
(
n
1
⋯
n
d
), and ε
i
(
n
1
⋯
n
d
) be (
n
1
⋯
n
d
)-dimensional column vectors defined by (40)
Z
i
n
1
⋯
n
d
≔
Z
i
j
1
n
1
,
…
,
j
d
n
d
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d
∈
R
n
1
⋯
n
d
g
i
n
1
⋯
n
d
≔
g
i
j
1
n
1
,
…
,
j
d
n
d
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d
∈
R
n
1
⋯
n
d
ε
i
n
1
⋯
n
d
≔
ε
i
j
1
n
1
,
…
,
j
d
n
d
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d
∈
R
n
1
⋯
n
d
.Furthermore, let Z
(
n
1
⋯
n
d
), g
(
n
1
⋯
n
d
), and E
(
n
1
⋯
n
d
) be (
n
1
⋯
n
d
)
×
p-dimensional matrices whose ith column is given, respectively, by the column vectors Z
i
(
n
1
⋯
n
d
), g
i
(
n
1
⋯
n
d
), and ε
i
(
n
1
⋯
n
d
), i
=
1
,
…
,
p. Then Model (35) can also be represented as follows:(41)
Z
n
1
⋯
n
d
=
1
n
1
⋯
n
d
g
n
1
⋯
n
d
+
E
n
1
⋯
n
d
,where, for u
,
v
=
1
,
…
,
p, C
o
v
(
ε
u
(
n
1
⋯
n
d
)
,
ε
v
(
n
1
⋯
n
d
)
)
=
σ
u
v
I
n
1
⋯
n
d, with I
n
1
⋯
n
d being the identity matrix in R
(
n
1
⋯
n
d
)
×
(
n
1
⋯
n
d
).Associated with the subspaceW
i
,
n
1
⋯
n
d we define the design matrix X
i
(
n
1
⋯
n
d
) as an element of R
(
n
1
⋯
n
d
)
×
d
i whose uth column is given by the (
n
1
⋯
n
d
)-dimensional column vector: (42)
f
i
u
j
1
n
1
,
…
,
j
d
n
d
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d
∈
R
n
1
⋯
n
d
,
u
=
1
,
…
,
d
i
,
i
=
1
,
…
,
p
.We denote the column space of X
i
(
n
1
⋯
n
d
) by C
(
X
i
(
n
1
⋯
n
d
)
)
⊂
R
(
n
1
⋯
n
d
)
×
d
i for the sake of brevity. We also define the column-wise projection of any matrix U
(
n
1
⋯
n
d
)
≔
(
U
1
(
n
1
⋯
n
d
)
,
…
,
U
p
(
n
1
⋯
n
d
)
) in R
n
1
⋯
n
d
×
p into the product space ∏
i
=
1
p
C
(
X
i
(
n
1
⋯
n
d
)
) by (43)
p
r
∏
i
=
1
p
C
X
i
n
1
⋯
n
d
U
n
1
⋯
n
d
≔
p
r
C
X
1
n
1
⋯
n
d
U
1
n
1
⋯
n
d
,
…
,
p
r
C
X
p
n
1
⋯
n
d
U
p
n
1
⋯
n
d
.A reasonable estimator of the covariance matrix Σ is denoted by Σ
^
n
1
⋯
n
d, defined by (44)
Σ
^
n
1
⋯
n
d
≔
1
n
1
⋯
n
d
p
r
∏
i
=
1
p
C
X
i
n
1
⋯
n
d
⊥
Z
n
1
⋯
n
d
⊤
p
r
∏
i
=
1
p
C
X
i
n
1
⋯
n
d
⊥
Z
n
1
⋯
n
d
,where (45)
p
r
∏
i
=
1
p
C
X
i
n
1
⋯
n
d
⊥
≔
I
n
1
⋯
n
d
-
p
r
∏
i
=
1
p
C
X
i
n
1
⋯
n
dconstitutes the component-wise orthogonal projector into the orthogonal complement of the product space ∏
i
=
1
p
C
(
X
i
n
1
⋯
n
d
).Zellner [1] and Arnold [22] investigated the consistency of Σ
^
n
1
⋯
n
d toward Σ in the case of the multivariate regression model with W
1
=
⋯
=
W
p. Some difficulties appear when the situation is extended to the case of W
1
≠
⋯
≠
W
p, since it involves the problem of finding the limit of matrices with the components given by inner products of two vectors.Theorem 5.
Suppose the localized model (41) is observed. If H
0 is true, then, under the conditions of Theorem 1, we have Σ
^
n
1
⋯
n
d
→
P
Σ.Proof.
IfH
0 is true, it can be easily shown that (46)
n
1
⋯
n
d
Σ
^
n
1
⋯
n
d
=
p
r
∏
i
=
1
p
C
X
i
n
1
⋯
n
d
⊥
E
n
1
⋯
n
d
⊤
p
r
∏
i
=
1
p
C
X
i
n
1
⋯
n
d
⊥
E
n
1
⋯
n
d
.For technical reason we assume without lost of generality that X
i
(
n
1
⋯
n
d
) is an orthogonal matrix, for i
=
1
,
…
,
p. Hence we further get the representation (47)
n
1
⋯
n
d
Σ
^
n
1
⋯
n
d
=
ε
u
n
1
⋯
n
d
⊤
ε
v
n
1
⋯
n
d
u
,
v
=
1
p
,
p
-
X
v
n
1
⋯
n
d
⊤
ε
u
n
1
⋯
n
d
⊤
X
v
n
1
⋯
n
d
⊤
ε
v
n
1
⋯
n
d
u
,
v
=
1
p
,
p
-
X
u
n
1
⋯
n
d
⊤
ε
u
n
1
⋯
n
d
⊤
X
u
n
1
⋯
n
d
⊤
ε
v
n
1
⋯
n
d
u
,
v
=
1
p
,
p
+
X
u
n
1
⋯
n
d
⊤
ε
u
n
1
⋯
n
d
⊤
X
v
n
1
⋯
n
d
⊤
ε
v
n
1
⋯
n
d
u
,
v
=
1
p
,
p
.Since ε
u
(
j
1
/
n
1
,
…
,
j
d
/
n
d
)
ε
v
(
j
1
/
n
1
,
…
,
j
d
/
n
d
)
u
,
v
=
1
p
,
p are independent and identically distributed random matrices with mean Σ, by the well-known weak law of large numbers, we get (48)
1
n
1
⋯
n
d
∑
j
1
=
1
n
1
⋯
∑
j
d
=
1
n
d
ε
u
j
1
n
1
,
…
,
j
d
n
d
ε
v
j
1
n
1
,
…
,
j
d
n
d
u
,
v
=
1
p
,
p
⟶
P
Σ
.Note that in the practice we consider the polynomial regression model. Hence, for every v
=
1
,
…
,
p, the design matrix X
v
(
n
1
⋯
n
d
) satisfies the so-called Huber condition (cf. Pruscha [31], pp. 115–117). By this reason, for the rest of the terms, we can immediately apply the technique proposed in Arnold [22] to show X
v
(
n
1
⋯
n
d
)
⊤
ε
u
(
n
1
⋯
n
d
)
→
D
N
d
v
(
0
,
σ
u
u
I
d
v
), for all u
,
v
=
1
,
…
,
p. Therefore, we finally get the following component-wise convergence: (49)
1
n
1
⋯
n
d
X
v
n
1
⋯
n
d
⊤
ε
u
n
1
⋯
n
d
⊤
X
v
n
1
⋯
n
d
⊤
ε
v
n
1
⋯
n
d
u
,
v
=
1
p
,
p
⟶
P
O
p
×
p
1
n
1
⋯
n
d
X
u
n
1
⋯
n
d
⊤
ε
u
n
1
⋯
n
d
⊤
X
u
n
1
⋯
n
d
⊤
ε
v
n
1
⋯
n
d
u
,
v
=
1
p
,
p
⟶
P
O
p
×
p
1
n
1
⋯
n
d
X
u
n
1
⋯
n
d
⊤
ε
u
n
1
⋯
n
d
⊤
X
v
n
1
⋯
n
d
⊤
ε
v
n
1
⋯
n
d
u
,
v
=
1
p
,
p
⟶
P
O
p
×
p
,where O
p
×
p is the p
×
p-zero matrix.Remark 6.
SinceΣ
^
n
1
⋯
n
d
-
1
/
2
V
n
1
⋯
n
d
(
R
n
1
⋯
n
d
)
=
Σ
^
n
1
⋯
n
d
-
1
/
2
Σ
1
/
2
Σ
-
1
/
2
V
n
1
⋯
n
d
(
R
n
1
⋯
n
d
), without altering the convergence result presented in Theorem 1, the population variance-covariance matrix Σ can be directly replaced by the consistence estimator Σ
^
n
1
⋯
n
d.


## 5. Calibration of the Tests

The limits of the test statistics are not distribution-free and we need therefore calibration for the distribution of the statistical tests. For the calibration we adapted the idea of residual based bootstrap for multivariate regression studied in Shao and Tu [32] for approximating the distributions of K
S
n
1
⋯
n
d
,
A and C
v
M
n
1
⋯
n
d
,
A.For fixedn
1
,
…
,
n
d, let F
n
1
⋯
n
d be the empirical distribution function of the vectors of least squares residuals {
r
j
1
⋯
j
d
-
R
¯
n
1
⋯
n
d
:
1
≤
j
k
≤
n
k
,
k
=
1
,
…
,
d
} centered at zero vector, where R
¯
n
1
⋯
n
d
≔
(
1
/
(
n
1
⋯
n
d
)
)
∑
j
1
=
1
n
1
⋯
∑
j
d
=
1
n
d
r
j
1
⋯
j
d. Let E
n
1
⋯
n
d
∗
≔
(
E
j
1
⋯
j
d
∗
)
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d be an array of independent and identically distributed random vectors sampled from F
n
1
⋯
n
d and let g
^
(
Ξ
n
1
⋯
n
d
) be the ordinary LSE of g
(
Ξ
n
1
⋯
n
d
), where g
^
(
Ξ
n
1
⋯
n
d
)
=
p
r
∏
i
=
1
p
W
i
Z
n
1
⋯
n
d. Then we generate the array of p-dimensional bootstrap observations which is denoted in this paper by Z
n
1
⋯
n
d
∗
≔
(
Z
j
1
⋯
j
d
∗
)
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d through the model: (50)
Z
n
1
⋯
n
d
∗
=
p
r
∏
i
=
1
p
W
i
Z
n
1
⋯
n
d
+
E
n
1
⋯
n
d
∗
.Based on this model we get the array of p-dimensional bootstrap least squares residuals which is given by the component-wise projection of the bootstrap observations: (51)
R
n
1
⋯
n
d
∗
≔
r
j
1
⋯
j
d
∗
j
1
=
1
,
…
,
j
d
=
1
n
1
,
…
,
n
d
=
Z
n
1
⋯
n
d
∗
-
p
r
∏
i
=
1
p
W
i
Z
n
1
⋯
n
d
∗
. Hence, the bootstrap analog of K
S
n
1
⋯
n
d
,
A and C
v
M
n
1
⋯
n
d
,
A is (52)
K
S
n
1
⋯
n
d
,
A
∗
≔
sup
A
∈
A
⁡
Σ
^
n
1
⋯
n
d
∗
-
1
/
2
V
n
1
⋯
n
d
R
n
1
⋯
n
d
∗
A
R
p
C
v
M
n
1
⋯
n
d
,
A
∗
≔
1
n
1
⋯
n
d
∑
A
∈
A
Σ
^
n
1
⋯
n
d
∗
-
1
/
2
V
n
1
⋯
n
d
R
n
1
⋯
n
d
∗
A
R
p
2
, where (53)
Σ
^
n
1
⋯
n
d
∗
≔
1
n
1
⋯
n
d
p
r
∏
i
=
1
p
C
X
i
n
1
⋯
n
d
⊥
Z
∗
n
1
⋯
n
d
⊤
p
r
∏
i
=
1
p
C
X
i
n
1
⋯
n
d
⊥
Z
∗
n
1
⋯
n
d
.The question regarding the consistency of the bootstrap approximation of thep-dimensional processes Σ
^
n
1
⋯
n
d
∗
-
1
/
2
V
n
1
⋯
n
d
(
R
n
1
⋯
n
d
∗
)
(
A
) for Σ
^
n
1
⋯
n
d
-
1
/
2
V
n
1
⋯
n
d
(
R
n
1
⋯
n
d
)
(
A
) is summarized in the following theorem.Theorem 7.
Let{
f
i
1
,
…
,
f
i
d
i
} be an ONB of  W
i, for i
=
1
,
…
,
p. Suppose the conditions of Theorem 1 are fulfilled. Then under H
0 it holds that (54)
Σ
∗
-
1
/
2
V
n
1
⋯
n
d
R
n
1
⋯
n
d
∗
⟶
D
B
p
-
p
r
∏
i
=
1
p
W
i
H
B
∗
B
p
,where p
r
∏
i
=
1
p
W
i
H
B
∗
B
p is defined in Theorem 1.Proof.
We notice that{
E
j
1
⋯
j
d
∗
:
1
≤
j
k
≤
n
k
,
k
=
1
,
…
,
d
} are independent and identically distributed with E
∗
(
E
1
⋯
1
∗
)
=
0 and Cov
∗
E
1
⋯
1
∗
=
Σ
^
n
1
⋯
n
d
-
R
¯
n
1
⋯
n
d
R
¯
n
1
⋯
n
d
⊤. Hence, the invariance principle implies that (
Σ
^
n
1
⋯
n
d
-
R
¯
n
1
⋯
n
d
R
¯
n
1
⋯
n
d
⊤
)
-
1
/
2
V
n
1
⋯
n
d
(
E
n
1
⋯
n
d
∗
) converges in distribution to B
p. Hence under H
0 we have R
n
1
⋯
n
d
∗
=
E
n
1
⋯
n
d
∗
-
p
r
∏
i
=
1
p
W
i
E
n
1
⋯
n
d
∗ and Σ
n
1
⋯
n
d
∗
-
1
/
2 can be written as (55)
Σ
n
1
⋯
n
d
∗
-
1
/
2
=
Σ
n
1
⋯
n
d
∗
-
1
/
2
Σ
^
n
1
⋯
n
d
-
R
¯
n
1
⋯
n
d
R
¯
n
1
⋯
n
d
⊤
1
/
2
Σ
^
n
1
⋯
n
d
-
R
¯
n
1
⋯
n
d
R
¯
n
1
⋯
n
d
⊤
-
1
/
2
,where it can be shown easily that (56)
R
¯
n
1
⋯
n
d
R
¯
n
1
⋯
n
d
⊤
=
o
P
1
,
E
∗
Σ
^
n
1
⋯
n
d
∗
=
Σ
^
n
1
⋯
n
d
+
o
P
1
,with o
P
(
1
) being the collection of the terms converging in probability to O
p
×
p. Then by recalling Theorem 5 and the linearity of V
n
1
⋯
n
d we only need to show that (57)
Σ
^
n
1
⋯
n
d
-
R
¯
n
1
⋯
n
d
R
¯
n
1
⋯
n
d
⊤
-
1
/
2
V
n
1
⋯
n
d
p
r
∏
i
=
1
p
W
i
E
n
1
⋯
n
d
∗
⟶
D
p
r
∏
i
=
1
p
W
i
H
B
∗
B
p
.The proof is established by imitating the steps of proving convergence result of Theorem 1. We are done.


## 6. Simulation Study

In this section, we report on a simulation study designed to investigate the finite sample size behavior of the KS and CvM tests. We simulate a multivariate model with four components defined on the unit rectangleI
2. The hypothesis under study is H
0
:
g
∈
∏
i
=
1
4
W
i against H
1
:
g
∉
∏
i
=
1
4
W
i, where W
1
≔
[
f
1
], W
2
≔
[
f
1
,
f
2
,
f
3
], W
3
≔
[
f
1
,
f
2
,
f
3
,
f
4
,
f
5
,
f
6
], and W
4
≔
[
f
1
,
f
2
,
f
3
,
f
4
,
f
5
,
f
6
,
f
7
]. Thereby we define f
1
(
t
,
s
)
≔
1, f
2
(
t
,
s
)
=
t, f
3
(
t
,
s
)
=
s, f
4
(
t
,
s
)
=
t
2, f
5
(
t
,
s
)
=
s
2, f
6
(
t
,
s
)
=
t
s, and f
7
(
t
,
s
)
=
t
3, for (
t
,
s
)
∈
I
2. The samples are generated from a localized model Y
(
l
/
n
1
,
k
/
n
2
)
=
(
1
/
n
1
n
2
)
g
(
l
/
n
1
,
k
/
n
2
)
+
E
(
l
/
n
1
,
k
/
n
2
) under the experimental design given by Ξ
n
1
×
n
2
=
{
(
l
/
n
1
,
k
/
n
2
)
∈
I
2
:
1
≤
l
≤
n
1
,
1
≤
k
≤
n
2
}, where g
(
l
/
n
1
,
k
/
n
2
) is defined as (58)
g
l
n
1
,
k
n
2
≔
5
+
ρ
exp
⁡
l
n
1
k
n
2
10
-
5
l
n
1
+
10
k
n
2
+
γ
exp
⁡
l
n
1
k
n
2
10
+
20
l
n
1
-
25
k
n
2
+
10
l
2
n
1
2
-
5
k
2
n
2
2
+
20
l
n
1
k
n
2
+
δ
exp
⁡
l
n
1
k
n
2
30
-
30
l
n
1
-
5
k
n
2
+
20
l
2
n
1
2
-
15
k
2
n
2
2
+
10
l
n
1
k
n
2
+
10
l
3
n
1
3
+
κ
exp
⁡
l
n
1
k
n
2
,for constants ρ, γ, δ, and κ determined prior to the generation of the samples. For fixed n
1 and n
2 and 1
≤
l
≤
n
1 and 1
≤
k
≤
n
2 the vector of random errors E
(
l
/
n
1
,
k
/
n
2
) is generated independently from the 4-variate normal distribution with mean zero and variance-covariance matrix given by (59)
Σ
4
×
4
≔
9
3
-
6
12
3
26
-
7
-
11
-
6
-
7
9
7
12
-
11
7
65
;however, we assume in the computation that Σ
4
×
4 is unknown. It is therefore estimated using Σ
^
n
1
n
2 defined in Section 4. It is important to note that for computational reason we restricted the index set to the Vapnick Chervonenkis Classes (VCC) of subsets of I
2 which is given by the family of closed rectangle with the point (
0,0
) as the essential point. That is, the family {
[
0
,
t
]
×
[
0
,
s
]
:
0
≤
t
,
s
≤
1
}.Whenρ, γ, δ, and κ are set simultaneously to zero then we get the samples which coincide to the model specified under H
0. Conversely, when at least one of them takes a nonzero value then we obtain the samples which can be regarded as from the alternative whose corresponding samples are generated by assigning nonzero values to either one of the constants, ρ, γ, δ, and κ, or the combinations of them.Table1 presents the empirical probabilities of rejection of H
0 for α
=
0.05 and some selected values of ρ, γ, δ, and κ. The empirical powers of the K
S and C
v
M tests are denoted by α
^
K
S and α
^
C
v
M, respectively. The notations σ
^
K
S and σ
^
C
v
M stand, respectively, for the standard deviation of the samples. The critical values of the statistics K
S
n
1
,
n
2
;
A and C
v
M
n
1
,
n
2
;
A are 6.0742 and 7.9910 which are approximated by simulation. For ρ
=
γ
=
δ
=
κ
=
0 the values of α
^
K
S and α
^
C
v
M fluctuate around 0.05 as it should be. This means that, independent of the selected number of the lattice points, both tests attain the specified level of significance.Table 1
The empirical probabilities of rejection of the KS and CvM tests for several selected values ofρ, γ, δ, and κ. The sample sizes are selected to 35
×
40 and 50
×
60. The simulation results are based on 10000 runs.
Sample size
ρ
γ
δ
κ
α
^
K
S
σ
^
K
S
α
^
C
v
M
σ
^
C
v
M
35
×
40
0
0
0
0
4.90%
1.31299
4.60%
2.52127
10
0
0
0
10.60%
1.52902
11.28%
4.11538
20
0
0
0
32.80%
1.71455
33.40%
6.26419
40
0
0
0
88.64%
1.86015
87.22%
10.37546
0
50
0
0
5.03%
1.28540
5.33%
2.63926
0
200
0
0
4.43%
1.22630
5.90%
2.57881
0
300
0
0
4.60%
1.18783
7.03%
2.52180
0
500
0
0
4.63%
0.90516
11.00%
2.32439
0
0
50
0
4.83%
1.27888
5.03%
2.59347
0
0
100
0
5.17%
1.29476
5.20%
2.68005
0
0
200
0
6.03%
1.28627
5.23%
2.53931
0
0
500
0
8.23%
1.27694
7.73%
2.71275
0
0
0
50
5.00%
1.27241
5.03%
2.69462
0
0
0
100
4.77%
1.23645
4.63%
2.56912
0
0
0
300
4.77%
1.27318
5.07%
2.53470
0
0
0
500
4.97%
1.26322
4.97%
2.55223
10
50
0
0
11.34%
1.55144
12.22%
3.82300
10
100
0
0
11.12%
1.49017
12.70%
3.81009
10
200
0
0
11.20%
1.52118
12.06%
3.71303
10
500
0
0
12.46%
1.10222
22.44%
2.87591
20
50
50
50
32.78%
1.76263
33.62%
5.75683
20
50
100
50
31.80%
1.71323
32.52%
5.98971
20
50
200
50
28.16%
1.68594
29.16%
5.79669
20
50
500
50
27.08%
1.52512
28.86%
4.92925
50
50
50
20
98.60%
1.86506
97.60%
13.43552
50
50
50
40
98.60%
1.83424
97.95%
13.10968
50
50
50
60
98.55%
1.91241
98.25%
13.49794
50
×
60
0
0
0
0
4.75%
1.26388
4.55%
2.52770
10
0
0
0
11.35%
1.48925
11.40%
3.70839
20
0
0
0
33.65%
1.75109
33.40%
5.82381
30
0
0
0
65.70%
1.82802
64.90%
8.31756
0
50
0
0
4.70%
1.26907
4.95%
2.60854
0
200
0
0
4.75%
1.23471
5.90%
2.53725
0
300
0
0
5.35%
1.19808
7.35%
2.62291
0
0
50
0
0.0485
1.27867
5.15%
2.57206
0
0
100
0
4.35%
1.23893
4.45%
2.40576
0
0
200
0
5.80%
1.27524
5.50%
2.62313
0
0
500
0
8.45%
1.28708
6.50%
2.64556
0
0
0
50
5.00%
1.28525
4.85%
2.71581
0
0
0
100
4.50%
1.25700
4.70%
2.57227
0
0
0
500
4.70%
1.25102
5.10%
2.51076
20
50
0
0
32.10%
1.70264
32.30%
5.72887
20
100
0
0
33.55%
1.74878
32.50%
5.76105
20
200
0
0
33.30%
1.72065
34.35%
5.71531
20
500
0
0
33.90%
1.61848
46.30%
5.77714
30
50
100
0
62.55%
1.78912
61.80%
8.08037
30
50
200
0
60.50%
1.79929
58.60%
7.76765
30
50
500
0
56.35%
1.78203
54.00%
7.76498
20
20
20
20
32.55%
1.70930
33.05%
5.60510
30
50
50
50
63.10%
1.86549
61.65%
8.63149
50
50
50
50
97.75%
1.88146
97.10%
13.27029
50
100
50
50
97.60%
1.89668
96.90%
13.30456Furthermore, Figure1 exhibits the graphs of the empirical power function of the KS and CvM tests for α
=
0.05 associated with hypothesis H
0 specified above against H
1
:
g
i
∉
W
i, for i
=
1,2
,
3,4. For the four cases we generate the error vectors independently from the same 4-variate normal distribution mentioned above. In the clockwise direction the left-top panel presents the graphs of the power function for testing H
0 against H
1
:
g
1
∉
W
1, the right-top panel is for H
0 against H
1
:
g
2
∉
W
2, the right-bottom is for H
0 against H
1
:
g
3
∉
W
3, and left-bottom is for H
0 against H
1
:
g
4
∉
W
4. The common characteristic of the tests is that the power gets larger as the the model moves away from H
0. The KS tests represented by smooth line tend to have slightly larger power. However, somewhat unexpectedly, in the second case, the CvM test has much larger power.Figure 1
The empirical power functions of the KS test (smooth line) and CvM test (dotted line) for the 4-variate model using50
×
60-point regular lattice. The simulation is under 10000 runs.


## 7. Example of Application

In this example, the proposed method is applied to a mining data studied in Tahir [20]. As introduced in Section 1, the data consist of a simultaneous measurement of the percentages of Nickel (Ni), Cobalt (Co), Ferum (Fe), and other substances like Calcium-Monoxide (CaO), Silicon-Dioxide (S
i
O
2), and Magnesium-Monoxide (MgO). The sample was obtained by drilling bores set according to a three-dimensional lattice of size 7
×
14
×
10 with 7 equidistance rows running west to east, 14 equidistance columns running south to north, and 10 equidistance depths from the surface of the earth to the bottom. To simplify the computation of the test statistics we consider the experimental design as a two-dimensional lattice of size 7
×
14 by taking the average value of the samples measured in the same position. We further assume that the exploration region is given by a closed rectangle so that by suitable rescaling it can be transformed into a closed unit rectangle I
2. Table 2 exhibits, respectively, the pairs scatter plot and Pearson’s correlation coefficient among the percentages of Ni, CaO, Co, the logarithm of the percentages of S
i
O
2 (l
o
g
S
i
O
2), MgO (logMgO), and Fe (logFe). By this reason a multivariate analysis must be conducted in the statistical modelling taking into account the unknown covariance matrix of the vector of the variables. Furthermore, based on the individual scatter plot of the samples which are not presented in this work, it can be inferred that polynomials of lower order seem to be adequate to approximate the population model. More precisely, let Y
≔
(
Y
1
,
Y
2
,
Y
3
,
Y
4
,
Y
5
,
Y
6
)
⊤ be the vector of observations representing the observed percentages of CaO, l
o
g
S
i
O
2, logMgO, Co, Ni, and logFe, respectively. We aim to test the hypothesis (60)
H
0
:
E
Y
=
β
11
β
21
+
β
22
t
+
β
23
s
β
31
+
β
32
t
+
β
33
s
β
41
+
β
42
t
+
β
43
s
β
51
+
β
52
t
+
β
53
s
+
β
54
t
2
+
β
55
t
s
+
β
56
s
2
β
61
+
β
62
t
+
β
63
s
+
β
64
t
2
+
β
65
t
s
+
β
66
s
2
,
0
≤
t
,
s
≤
1
,for some unknown constants β
i
j, i
=
1,2
,
3,4
,
5,6 and j
=
1
,
…
,
d
i, with d
1
=
1, d
2
=
d
3
=
d
4
=
3 and d
5
=
d
6
=
6. For this case we have W
1
≔
[
f
1
], W
2
=
W
3
=
W
4
≔
[
f
1
,
f
2
,
f
3
], and W
5
=
W
6
≔
[
f
1
,
f
2
,
f
3
,
f
4
,
f
5
,
f
6
], with f
1
(
t
,
s
)
≔
1, f
2
(
t
,
s
)
=
t, f
3
(
t
,
s
)
=
s, f
4
(
t
,
s
)
=
t
2, f
5
(
t
,
s
)
=
t
s, and f
6
(
t
,
s
)
=
s
2.Table 2
Pearson’s correlation matrix of the percentages of CaO,log
⁡
S
i
O
2, log
⁡
M
g
O, Ni, log
⁡
F
e, and Co observed over a regular lattice of size 7
×
14. Source of data: Tahir [20].
CaO
log
⁡
S
i
O
2
log
⁡
M
g
O
Ni
log
⁡
F
e
Co
C
a
O
1.0000
0.3949
0.4045
-
0.1285
-
0.1167
-
0.0665
log
⁡
S
i
O
2
0.3949
1.0000
0.8459
-
0.0003
-
0.5414
-
0.4556
log
⁡
M
g
O
0.4045
0.8459
1.0000
-
0.1331
-
0.4968
-
0.3134
Ni
−0.1285
-
0.0003
-
0.1331
1.0000
0.1652
0.1068
log
⁡
F
e
−0.1166
-
0.5414
-
0.4968
0.1652
1.0000
0.5937
Co
−0.0665
-
0.4556
-
0.3134
0.1068
0.5937
1.0000We obtained the valuesK
S
7
×
14
;
A
=
1.44109 and C
v
M
7
×
14
;
A
=
0.936687 with the associated simulated p values of 0.98350 and 0.93670, respectively. We notice that in the computation we consider the VCC {
[
0
,
t
]
×
[
0
,
s
]
:
0
≤
t
,
s
≤
1
} as the index sets instead of A. Hence when using the KS test as well as CvM test the hypothesis will be not rejected for almost all commonly used values of α. There exists a significant evidence that the assumed model is appropriate for describing the functional relationship between the experimental conditional and the percentages of those elements.In the practice some computational difficulties appear for testing using our proposed method. First, to the knowledge of the authors, the analytical formula for computing the critical andp values of the tests have been not yet available in the literatures; therefore we need to approximate them by simulation using computer. Second, although the test procedures are established for a much larger family of sets A, in the application the computation is always restricted to the VCC of subsets of I
d like that of {
∏
i
=
1
d
[
0
,
t
i
]
⊂
I
d
:
0
<
t
i
≤
1
,
i
=
1
,
…
,
d
} or {
∏
i
=
1
d
[
s
i
,
t
i
]
⊂
I
d
:
0
≤
s
i
<
t
i
≤
1
,
i
=
1
,
…
,
d
}.


## 8. Concluding Remark

In this article we have developed an asymptotic method for checking the validity of a general multivariate spatial regression model by considering the multidimensional set-indexed partial sums of the residuals. For the calibration of the distribution of the test statistics we propose the residual based bootstrap for multivariate regression. It is shown by applying imitation technique that the residual bootstrap resampling technique is consistent. In a simulation study the finite sample size behavior of the KS and CvM statistics is investigated in greater detail. For the first-order model CvM test has much larger power, whereas for constant, second-order, and third-order models the powers of the two tests are almost the same.Other possibilities of tests for multidimensional case can be obtained by incorporating a sampling technique according to an arbitrary experimental design. Sometimes because of technical, economic, or ecological reason, practitioners will not or cannot sample the observations equidistantly. One possible approach is to sample according to a continuous probability measure; see, for example, the sampling method proposed in Bischoff [11]. Under this concern we get the so-called weighted KS and CvM tests which can be viewed as generalization of the KS and CvM tests studied in this paper.Instead of considering the least squares residuals of the observations we can also define a test by directly investigating the partial sums of the observations. The limit process will be given by a type of signal plus a noise which is given by the multidimensional set-indexed Brownian sheet. Observing the limit process we can formulate likelihood ratio test based on the Cameron-Matrin-Girsanov density formula of the limit process. Establishing such type of test will be of our concern in our future research project.



---

*Source: 2601601-2016-09-07.xml*