Simultaneous confidence bands for contrasts between
several nonlinear regression curves
Xiaolei Lu∗ and Satoshi Kuriki†
Abstract — We propose simultaneous confidence bands of the hyperbolic-type for the contrasts between several nonlinear(curvilinear) regression curves. The critical value of a confidence band is determined from the distribution of the maximumof a chi-square random process defined on the domain of explanatory variables. We use the volume-of-tube method toderive an upper tail probability formula of the maximum of a chi-square random process, which is suﬃciently accurate
in commonly used tail regions. Moreover, we prove that the formula obtained is equivalent to the expectation of theEuler-Poincaré characteristic of the excursion set of the chi-square random process, and hence conservative. This result istherefore a generalization of Naiman’s inequality for Gaussian random processes. As an illustrative example, growth curvesof consomic mice are analyzed.

Keywords: Chi-square random process; expected Euler-characteristic heuristic; growth curve; Naiman’s inequality; volume-of-tube method.
1 Introduction
This paper concerns multiple comparisons of k (≥ 3) nonlinear (curvilinear) regression curves estimatedfrom independent k groups. For each group i = 1, . . . , k, let
βTi f (x), x ∈ X,(1.1)

be the true regression curve, where f (x) = ( f1 (x), . . . , fp (x))T is a known regression basis vector function,βi = (βi1 , . . . , βip )T is an unknown parameter vector, and X ⊆ R is the domain of explanatory variables.
Suppose that an estimator b
βi of βi is available whose distribution is the normal distribution Np (βi , r−1Σ),iwhere ri is a known constant representing the size of the group i. The p × p matrix Σ is supposed tobe known, or at least known up to a constant Σ = σ2 Σ0 . In the case of the latter, we suppose that anindependent estimator b
σ2 of σ2 is available.
∑Let C denote the set of vectors c = (c1 , . . . , ck )T such that ki=1 ci = 0. The focus of this paper is the∑construction of 1 − α simultaneous confidence bands for all the contrasts ki=1 ci βTi f (x) between the kregression curves for all x ∈ X and c ∈ C, where X is a given finite interval [a, b]. Specifically, according
to the traditional form of the point estimate plus or minus a probability point times the estimatedstandard error, we construct a 1 − α simultaneous confidence band of the formvu
t kkk∑∑∑ c2i  f (x)T Σ f (x),ci βTi f (x) ∈ci b
βTi f (x) ± bα (1.2)
ri i=1i=1

i=1where bβTi f (x) is the estimator of βTi f (x) in (1.1). The critical value bα is determined such that the eventin (1.2) for all x ∈ X and c ∈ C holds with a probability of at least 1 − α.
Our problem typically arises from the following experimental design. Suppose that for each groupi, and for each explanatory variable x j , j = 1, . . . , n, we have observations yi j1 , . . . , yijri as objective
variables with ri replications, which are assumed to follow the modelyi jh = βTi f (x j ) + εijh ,i = 1, . . . , k, j = 1, . . . , n, h = 1, . . . , ri .

(1.3)

∗ Department of Statistical Sciences, Graduate University for Advanced Studies, 10-3 Midoricho, Tachikawa, Tokyo 190-8562,
Japan† The Institute of Statistical Mathematics, 10-3 Midoricho, Tachikawa, Tokyo 190-8562, Japan, Email: kuriki@ism.ac.jp, Tel:
+81-50-5533-8423

Here, random errors εijh are assumed to be independently distributed as the normal distribution
N(0, σ(x j )2 ). Then, the least squares estimator b
βi of βi has the multivariate normal distribution Np (βi , r−1Σ),iwhere−1n∑1TΣ = f(x)
f(x)
jj

σ(x j )2j=1

is the inverse of the p × p information matrix. Data with this structure often appear in growth curve
analysis and longitudinal data analysis.
2 Main result
Recall that our aim is to obtain the threshold bα such that the event in (1.2) for all x ∈ X and c ∈ Cholds with a probability of at least 1 − α. For this purpose, we adopt the volume-of-tube method(Takemura and Kuriki, 2002). This is an integral geometric approach originating from Naiman (1986),and is equivalent to the expected Euler-characteristic heuristic (Adler and Taylor, 2007) in our setting.
In the unit sphere Sp−1 of the p-dimensional Euclidean space, define a trajectory Γ of a normalizedbasis vector function by{}1

Γ = ψ(x) | x ∈ X ⊂ S

p−1withψ(x) =

Σ 2 f (x)1∥Σ 2 f (x)∥.

Then, bα can be obtained from geometric quantities of Γ.
Theorem (Lu and Kuriki, 2016). Let{ ())}()(Γ( k )bP(b) = √ 2 k−1 |Γ| P χ2k > b2 − P χ2k−2 > b2 + P χ2k−1 > b2 ,πΓ( 2 )where |Γ| is the length of the curve Γ. Let b = bα be the solution of bP(b) = α. Then, bα is a conservative andasymptotically exact (when α → ∞) threshold for the 1 − α simultaneous confidence bands (1.2).
Note that when the map ψ : X → Sp−1 is one-to-one,∫d|Γ| =
ψ(x) dx.
dx
X(2.1)

3 Growth curve analysis
Mice are one of the most popular model organisms, and are often used in genomic research. Figure 3.1depicts the average body weights of male mice from four diﬀerent strains measured from 2 to 20 weeks
after birth. The four strains are C57BL/6 (referred to as B6), MSM/Ms (MSM), B6-Chr17MSM (B6-17),and B6-ChrXTMSM (B6-XT). Among these, B6 is the most common laboratory strain and serves as thestandard. MSM is a wild-derived strain that has contrasting properties to B6 such as non-black color,small size, and aggressive behavior. B6-17 and B6-XT are artificial strains known as consomic micemade from B6 and MSM. B6-17 has all the chromosomes from B6, but only chromosome 17 from MSM;
B6-XT has all the chromosomes from B6, but only half of the X chromosome from MSM. By comparing
the consomic strains with B6, we expect to reveal the role of each chromosome.
The dataset we used is publicly available as Supplemental Table S1 of Takada, et al. (2008). In theirexperiments, the weight (unit: gram) yijh of the hth individual from strain i was measured at time pointx j . The measurement time points were {x1 , . . . , x10 } = {2, 4, . . . , 20} (n = 10). This dataset includes theaverage body weight yij of strain i at time x j , and its standard error
1∑
yij =yijh ,rivtri

h=1c i j) =s.e.(y1 ∑(yijh − yij )2 ,r2i h=1ri

35
30
25
20weight (g)

15
10
5

5

10

15

20weekFigure 3.1: Average body weights of mice from four strains.
(sample mean: ◦ (B6), + (B6-17), ⋄ (B6-XT), △ (MSM);
fitted curve: — (B6), · · · (B6-17), – – (B6-XT), – · – (MSM))

as well as the number ri of individuals of strain i.
In the following analysis, we use k = 3 groups (strains) B6 (i = 1), B6-17 (i = 2), and B6-XT (i = 3).
The number of individuals are r1 = 12, r2 = 24, and r3 = 12.
We fit the model (1.3) to these data. We estimate the variance asbσ(x j )2 = ∑k1

rik ∑∑

i=1 (ri − 1) i=1 h=1(yijh − yij )2 = ∑k1k∑

i=1 (ri − 1) i=1c i j )2 ,r2i s.e.(ywhich is used as the true value σ(x j )2 hereafter. One particular feature of this dataset is that theexperiment is well controlled and the measurement errors are quite small.
As the basis function f (x), we consider a family of basis functions( ())
x−2f (x) = fd,m (x) = Bd(m − d) − (i − d − 1),
20 − 21≤i≤mwhere
Bd (x) =(d+1∑(−1)
r=0d+1−r d + 1r)(r − x)d+d!

(de Boor, 1978, page 89). fd,m (x) consists of m B-spline bases with equally-spaced knots at intervals of(20 − 2)/(m − d). Note that fd,m (x) is piecewise of class Cd .
In the range d = 2, 3, 4 and m = d+1, d+2, . . . , n (= 10), we searched for the best model that minimizes
AIC and BIC. In both criteria, the minimizer was (d, m) = (2, 5), which we use as the true value hereafter.
Suppose that we are interested in the period X = [a, b] = [2, 20]. An approximate value of the lengthof Γ in (2.1) is given byN∑
|Γ| ≈ψ(xt ) − ψ(xt−1 ) ,t=1

where xt = a+t(b−a)/N, t = 0, 1, . . . , N. When N = 10, 000, the approximate value of |Γ| is 6.989 = 2.225π.
Using this value, the critical value is bα = 3.258 (α = 0.05).

To compare k groups, various types of contrasts are used. For a pairwise comparison between groupi and group j, we choose c = (. . . , 0, 1 , 0, . . . , 0, −1, 0, . . .). For the comparison of groups {i, j} and groupithk, we use(c = . . . , 0,jth)rjri, 0, . . . , 0,, 0, . . . , 0, −1, 0, . . . .
kthri + r j
ri + r j
ithjthFigure 3.2 depicts the diﬀerence curves of strains B6-17 vs. B6 (left) and B6-XT vs. B6 (right), and their
95% simultaneous confidence bands. In the left panel, the horizontal line representing zero diﬀerence isalmost between the confidence bands. This means that there is no significant diﬀerence between B6-17
and B6. In contrast, in the right panel, after about week 14, the horizontal line is outside the confidence
bands, thereby indicating that B6-XT and B6 are diﬀerent during this period.

difference (g)
−1 0
1
2

B6−XT vs. B6

difference (g)
−1 0
1
2

B6−17 vs. B6510week

15

20

510week

15

20

Figure 3.2: Diﬀerences of body weights and 95% confidence bands.
For a fixed x, the test statistic for the null hypothesis H0,x : βT1 f (x) = · · · = βTk f (x) is
2∑kkTb
∑


r
βf(x)
i
1i=1i
2T .
bχ (x) =rβf(x)−i∑ikf (x)T Σ f (x) i=1 i=1 ri
For a fixed x, the null distribution is the chi-square distribution with k −1 degrees of freedom. However,for the overall null hypothesis H0 : βT1 f (x) = · · · = βTk f (x) for all x ∈ X, the distribution of the maximumof the chi-square random process should be used. The critical value b2α can be used for this purpose. In
this dataset, as already shown in Figure 3.2, after about week 14, the hypothesis of equality is rejected.
References[1] Adler, R. J. and Taylor, J. E. (2007). Random Fields and Geometry, Springer.
[2] de Boor, C. (1978). A Practical Guide to Splines, Springer.
[3] Lu, X. and Kuriki, S. (2016). Simultaneous confidence bands for contrasts between several nonlinear
regression curves, arXiv:1510.05077.
[4] Naiman, D. Q. (1986). Conservative confidence bands in curvilinear regression. The Annals of Statistics, 14 (3), 896–906.
[5] Takada, T., Mita, A., Maeno, A., Sakai, T., Shitara, H., Kikkawa, Y., Moriwaki, K, Yonekawa, H., andShiroishi, T. (2008). Mouse inter-subspecific consomic strains for genetic dissection of quantitativecomplex traits, Genome Research, 18 (3), 500–508. http://genome.cshlp.org/content/18/3/500[6] Takemura, A. and Kuriki, S. (2002). On the equivalence of the tube and Euler characteristic methods
for the distribution of the maximum of Gaussian fields over piecewise smooth domains. The Annalsof Applied Probability, 12 (2), 768–796.