Model Testing Based on Regression Spline

Abstract

Tests based on regression spline are developed in this chapter for testing nonparametric functions in nonparametric, partial linear and varying-coefficient models, respectively. These models are more flexible than linear regression model. However, one important problem is if it is really necessary to use such complex models which contain nonparametric functions. For this purpose, p-values for testing the linearity and constancy of the nonparametric functions are established based on regression spline and fiducial method. In the application of spline-based method, the determination of knots is difficult but plays an important role in inferring regression curve. In order to infer the nonparametric regression at different smoothing levels (scales) and locations, multi-scale smoothing methods based on regression spline are developed to test the structures of the regression curve and compare multiple regression curves. It could sidestep the determination of knots; meanwhile, it could give a more reliable result in using the spline-based method.

1. Introduction

It is well known that the model which contains nonparametric functions, such as partial linear model and varying-coefficient model, plays an important role in applications due to its flexible structure. However, in practice, investigators often want to know whether it is really necessary to fit the data with such more complex models rather than a simpler model. This amounts to testing the linearity of nonparametric functions in a regression model. In this chapter, we first consider the following three frequently used regression models.

Nonparametric regression model:

$$y=f\left(x\right)+\epsilon .$$E1

Partial linear regression model:

$$y=Z^{\prime }b+f\left(x\right)+\epsilon .$$E2

Varying-coefficient model:

$$y=z_1{f}_1\left(x_1\right)+\dots +z_p{f}_p\left(x_p\right)+\epsilon .$$E3

In models (1)–(3), $y$ is the response variable, $Z=\left(z_1,\cdots ,z_p\right)$ is a p-dimensional regressor, $x$ and $x_1,\cdots ,x_p$ are covariant taking values in a finite interval, $\epsilon$ is the error, b is a parameter vector, and f(x) and $f_j\left(x_j\right),j=1,2,\dots ,p$ are unknown smooth functions. Usually we suppose that $\left(z,x\right)$ and $\mathrm{\epsilon }$ are independent and $\epsilon ˜F\left(\bullet /\sigma \right)$, where $F$ is a known cumulate distribution function (cdf) with mean 0 and variance 1; $\sigma$ is unknown. Without loss of generality, we can suppose that $x$ and $x_1,\cdots ,x_p$ take values in [0, 1]. We try to test the linearity of $f\left(x\right)$ in models (1) and (2) and the constancy of $f_j\left(x_j\right)$ in model (3) for some $j\in \left(1,2,\cdots ,p\right)$.

The hypothesis testing in nonparametric regression model was considered in many papers. Härdle and Mammen [1] developed the visible difference between a parametric and a nonparametric curve estimates. Based on smoothing techniques, many tests were constructed for testing the linearity in regression model; see Hart [2], Cox et al. [3], and Cox and Koh [4] for a review. Recently, Fan et al. [5] studied a generalized likelihood ratio statistic, which behaves well in large sample case. Tests based on penalized criterion were developed by Eubank and Hart [6] and Baraud [7].

The linearity of partial linear regression model (2) was studied by Bianco and Boente [8], Liang et al. [9], and Fan and Huang [10]. There are also many other papers concerning such testing problems (see [11, 12, 13, 14, 15, 16], among others). The constancy of the functional coefficient $f_j\left(x_j\right)$ in varying-coefficient model (3) was studied in Fan and Zhang [17], Cai et al. [18], Fan and Huang [19], You and Zhou [20], and Tang and Cheng [21]. Local polynomials and smoothing spline methods to estimate the coefficients in model (3) can be seen in Hoover et al. [22], Wu et al. [23], and so on.

The critical values of most of the previous tests were obtained by Wilks theorem or bootstrap method. So such tests only behave well in the case of relatively large sample size. This chapter would give some testing procedures based on regression spline and the fiducial method [24] in Section 2. It has a good performance even when the sample size is small.

In using the regression spline, the key problem is the determination of knots used in spline interpolation. As we know that, for smoothing methods such as kernel-based method and smoothing spline, the smoothness is controlled by smoothing parameters. For the well-known kernel estimate, the bandwidth that is extremely big or small might leads to over-smoothing or under-smoothing, respectively. In order to avoid the selection of an optimal smoothing parameter, multi-scale smoothing method was introduced by Chaudhuri and Marron [25, 26] based on kernel estimation for exploring structures in data. This multi-scale method is known as significant zero crossings of derivatives (SiZer) methodology. The basic idea of SiZer is to infer a nonparametric model by using a wide range of smoothing parameter (bandwidth) values rather than only using one “optimal” value in some sense.

There have been many versions of SiZer for various applications, such as the local likelihood version of SiZer in Li and Marron [27], the robust version of SiZer in Hannig and Lee [28], and the quantile version of SiZer in Park et al. [29]. In addition, Marron and deUñaÁlvarez [30] applied SiZer to estimate length biased, censored density, and hazard functions; Kim and Marron [31] utilized SiZer for jump detection and Park and Kang [32] applied SiZer to compare regression curves. The smoothing spline version of SiZer was proposed by Marron and Zhang [33]. It used the tuning parameter (penalty parameter) that controls the size of penalty as the smoothing parameter is.

Comparing with bandwidth for kernel-based method and tuning parameter for smoothing spline, it is more difficult to determine the number of knots and their positions. For this reason a multi-scale smoothing method based on regression spline is proposed in Section 3 to test the structures of nonparametric regression model. The proposed multi-scale method does not involve the determination of the “best” number of knots and can be extended easily to a more general case.

2. Tests for nonparametric function based on regression spline

In this section, the linearity of function $f\left(x\right)$ in model (1) is tested based on regression spline and fiducial method. Then, the proposed test procedure for model (1) is extended to test the linearity of model (2) and the constancy of function coefficient in model (3), respectively.

2.1. Test the linearity of nonparametric regression model

Without loss of generality, we suppose that $x$ in model (1) takes values in [0, 1] and the set of knots is T = {0 $=t_1<t_2,\cdots ,<t_m=1\}$. In order to estimate model (1), nonparametric function f(x) is fitted by $k$th order splines with knots T. This means that

$$f\left(x\right)\approx \sum _{j=1}^{m+k-1}{\beta }_j{g}_j\left(x\right),$$E4

where ${\beta }_j$ is coefficient and $g_j\left(x\right),$ $j=1,2,\cdots ,m+k-1,$ is basis function for order k splines, over the knots $t_1,t_2,\cdots ,t_m.$

With $n$-independent observations $Y=\left(y_1,y_2,\cdots ,y_n\right)\in {ℝ}^n$, the basis matrix $G_{n\times \left(m+k-1\right)}$ is defined by $G=\left\{g_j\left(x_i\right)\right\},x_i\text{is the designed point},i=1,2,\cdots ,n;j=1,2,\cdots ,m+k-1.$ Hence, model (1) can be approximated as $Y\approx \mathit{G\beta }+\epsilon$. The least squares estimator of coefficients is

$$\widehat{\beta }={\left(G^{T}G\right)}^{-1}{G}^{T}Y,$$E5

and the estimator of $f\left(x_i\right)$ can be expressed as

$$\widehat{Y}={\left\{\widehat{f}\left(x_1\right),\widehat{f}\left(x_2\right),\cdots ,\widehat{f}\left(x_n\right)\right\}}^T=G{\left(G^{T}G\right)}^{-1}{G}^{T}Y.$$E6

For testing the linearity of model (1), linear spline is used to approximate $f\left(x\right)$. It means that basis function $g_j\left(x\right)$ is a linear function:

$$g_1\left(x\right)=\frac{-x-t_2}{\left(t_2-t_1\right){\mathbb{l}}_2\left(t\right)},$$

$$g_{k-1}\left(x\right)=\frac{x-t_{k-2}}{\left(t_{k-1}-t_{k-2}\right){\mathbb{l}}_{k-1}\left(t\right)}-\frac{x-t_k}{\left(t_k-t_{k-1}\right){\mathbb{l}}_k\left(t\right)},3\le k\le m,$$E7

$$g_m\left(x\right)=\frac{x-t_{m-1}}{\left(t_m-t_{m-1}\right){\mathbb{l}}_m\left(t\right)}.$$

In this case, the approximated function in (4) is a linear interpolation with $k$ =1. The true value is ${\beta }_j=f\left(t_j\right),j=1,2,\cdots ,m$. The linearity of function $f\left(x\right)$ can be written as

$${\mathrm{H}}_0:\frac{{\beta }_2-{\beta }_1}{t_2-t_1}=\frac{{\beta }_3-{\beta }_2}{t_3-t_2}=\cdots =\frac{{\beta }_m-{\beta }_{m-1}}{t_m-t_{m-1}}.$$

Null hypothesis ${\mathrm{H}}_0$ can be expressed in matrix as $L\prime \beta =0$,

where

$$L^{\prime }=\left[\begin{array}{lll}{h}_2-h_1-h_2& h_{1}0\cdots & 000\\ \cdots \cdots & \cdots \cdots \cdots & \cdots \cdots \cdots \\ 00& 00\cdots & h_{m-1}-h_{m-1}-h_{m-2}{h}_{m-2}\end{array}\right],$$

where $h_j=t_{j+1}-t_j,j=1,2,\cdots ,m-2$. Null hypothesis ${\mathrm{H}}_0$ is equivalent to the following one:

$${\mathrm{H}}_0^{\ast }:L^{\prime }\beta =0.$$E8

The p-value for testing hypothesis ${\mathrm{H}}_0^{\ast }$ will be derived by the fiducial method in the following context. Assume that matrix $G$ has full rank, and let $\epsilon ˜\mathit{\sigma N}\left(0,1\right)$. In model $Y=\mathit{G\beta }+\epsilon$, the sufficient statistic of $\left(\beta ,{\sigma }^2\right)$ is $\left(\widehat{\beta },S^2\right),$ where $\widehat{\beta }$ is defined in (5) and

$$S^2=Y^{\prime }\left(I-P_G\right)Y,P_G=G{\left(G^{\prime }G\right)}^{-1}{G}^{\prime }.$$

By Dawid and Stone [34], the sufficient statistic can be represented as a functional model:

$$\widehat{\beta }=\beta +\sigma {\left(G^{\prime }G\right)}^{-\frac{1}{2}}{E}_1,S=\sigma E_2,E=(E_1,E_2)\sim Q,$$E9

where $Q$ is the probability measure of $E=(E_1$, $E_2)$ and $E_1\sim N\left(0,I_m\right)$ and, independently, $E_2^2\sim {\chi }^2\left(n-m\right).$ From linear regression model, the fiducial model of $\beta$ can be obtained:

$$\widehat{\beta }=\beta +\frac{S}{E_2}{\left(G^{\prime }G\right)}^{-\frac{1}{2}}{E}_1,E=(E_1,E_2)\sim Q.$$E10

Given $\left(\widehat{\beta },S^2\right),$ the distribution of the right side in fiducial model is the fiducial distribution of $\beta$. That is, the fiducial distribution of $\beta$ is the conditional distribution of $R\left(E;\widehat{\beta },S^2\right)$ when $\left(\widehat{\beta },S^2\right)$ is given, where

$$R\left(E;\widehat{\beta },S^2\right)=\widehat{\beta }-\frac{S}{E_2}{\left(G^{\prime }G\right)}^{-\frac{1}{2}}{E}_1.$$E11

For testing hypothesis ${\mathrm{H}}_0^{\ast }$, the p-value is defined as

$$p\left(\widehat{\beta },S^2\right)=Q\left({‖L^{\prime }\left[R\left(E;\widehat{\beta },S^2\right)-{\mathrm{E}}_{Q}R\left(E;\widehat{\beta },S^2\right)\right]‖}_{\mathrm{\Sigma }}^2\ge {‖L^{\prime }{\mathrm{E}}_{Q}R\left(E;\widehat{\beta },S^2\right)‖}_{\mathrm{\Sigma }}^2\right),$$E12

where $Q$(·) and ${\mathrm{E}}_Q$ express the probability for an event and the expectation of a random variable under $Q$, respectively, and $\mathrm{\Sigma }$ is the conditional covariance matrix of $L^{\prime }{\mathrm{E}}_{Q}R\left(E;\widehat{\beta },S^2\right)$ given $\widehat{\beta },S^2$ and ${‖v‖}_{\mathrm{\Sigma }}^2=v^{\prime }{\mathrm{\Sigma }}^{-1}v$ for a vector $v.$

According to the definition of generalized pivotal quantity in [35], $R\left(E;\widehat{\beta },S^2\right)$ is a generalized pivotal quantity and also a fiducial pivotal quantity about $\beta$. Naturally, $L^{\prime }R\left(E;\widehat{\beta },S^2\right)$ is the fiducial pivotal quantity about $L^{\prime }\beta$. With the definition of $Q$ in Eq. (10), we have that

$$p\left(\widehat{\beta },S^2\right)=1-F_{m-2,n-m}\left(\frac{\left(n-m\right){\widehat{\beta }}^{\prime }L{\left(L^{\prime }{\left(G^{\prime }G\right)}^{-1}L\right)}^{-1}{L}^{\prime }\widehat{\beta }}{\left(m-2\right)S^2}\right),$$E13

where $F_{m-2,n-m}$ is the cdf of $F$-distribution with degrees of freedom $m-2\text{and}n-m$.

Under model (1) and the hypothesis that $f\left(x\right)$ is a linear function, null hypothesis ${\mathrm{H}}_0^{\ast }$ given in (8) is true. Suppose that the error is normally distributed, then the p-value given in Eq. (12) distributes as uniform distribution on interval (0, 1). On the other hand, under some mild condition, the test procedure based on $p\left(\widehat{\beta },S^2\right)$ is consistent. Which means that $p\left(\widehat{\beta },S^2\right)$ tends to be zero in probability 1 if ${\mathrm{H}}_0^{\ast }$ is false. The corresponding theoretical proof of the large sample properties and finite sample properties of $p\left(\widehat{\beta },S^2\right)$ is the same as the proof given in Li et al. [36].

In applications, we need to check some hypotheses as follows:

$${\mathrm{H}}_{01}:f\left(x\right)=C⟺{\beta }_1={\beta }_2=\cdots ={\beta }_m,$$

$${\mathrm{H}}_{02}:f\left(x\right)=\mathit{Cx}⟺\frac{{\beta }_2-{\beta }_1}{t_2-t_1}=\frac{{\beta }_3-{\beta }_2}{t_3-t_2}=\cdots =\frac{{\beta }_m-{\beta }_{m-1}}{t_m-t_{m-1}},\text{and},{\beta }_1=0.$$

The p-values for testing ${\mathrm{H}}_{01}$ and ${\mathrm{H}}_{02}$ can be obtained by replacing $L$ in (12) by $L_{01}$ and $L_{02}$, respectively, where $L_{02}=\left(e_1,L\right),e_1={\left(1,0,0,\cdots ,0\right)}^{\prime }$ and

$$L_{01}=\left[\begin{array}{ccc}{h}_2-h_1& 0\cdots & 00\\ \cdots \cdots & \cdots & \cdots \cdots \\ 00& 0\cdots & h_m-h_{m-1}\end{array}\right].$$E14

2.2. Test the linearity of partial linear model

To test the linearity of model (2), p-value can be established analogously. With $n$-independent observations $Y=\left(y_1,y_2,\cdots ,y_n\right)\in {ℝ}^n$, model (2) can be represented as.

$$y_i=Z_i^{\prime }b+f\left(x_i\right)+{\epsilon }_i,i=1,2,\cdots ,n,$$

where $Z_i^{\prime }={\left(z_{i1},\cdots ,z_{\mathit{ip}}\right)}^{\prime },b={\left(b_1,\cdots ,b_p\right)}^{\prime },x_i,i=1,2,\cdots ,n\mathrm{are}\text{fixed designed points}$. With the approximation of $f\left(x\right)$ given in (4), model (2) can be approximated by $Y\approx \mathit{X\theta }+\epsilon$, where $X={\left(\mathbb{Z},G\right)}_{n\times \left(p+m+1\right)};$ $\mathbb{Z}=$ ($z_{\mathit{ij}}$)$;i=1,2,\cdots ,n;j=1,2,\cdots ,p;G$ is the same as above; and $\theta ={\left(b^{\prime },{\beta }^{\prime }\right)}^{\prime }$. Then p-value for testing the linearity of model (2) can be defined by replacing $G$ in (12) by $X$, $\beta$ by $\theta ,$ and $L$ by $L_{03}$, respectively,$L_{03}={\left(0_{\left(m-2\right)\times p},L^{\prime }\right)}^{\prime }$.

The large sample and finite sample properties of the testing procedure for model (2) are the same as the test procedure for model (1).

2.3. Test the constancy of functional coefficient in varying-coefficient model

For model (3), investigators often want to know whether the coefficients are really varying; this means to test the constancy of the coefficient functions, that is, testing hypothesis:

$${\mathrm{H}}_{31}:f_j\left(x\right)=C_j\text{for}j=1,2,\cdots ,p\text{and some constant}{C}_j,$$E15

$${\mathrm{H}}_{32}:f_{j0}\left(x\right)=C_{j0}\text{for some}j=j_0\text{and some constant}{C}_{j0}.$$E16

With the set of knots T = {0 $=t_1<t_2,\cdots ,<t_m=1\},$ coefficient $f_j\left(x\right)$ can also be approximated by

$$f_j\left(x\right)=\sum _{k=1}^m{\beta }_{\mathit{js}}{g}_j\left(x\right),j=1,2,\cdots ,p,$$

where the true value of ${\beta }_{\mathit{js}}=f_j\left(t_k\right)$. Basic functions $g_j,j=1,2,\cdots ,m+1$ were defined in (7). The varying-coefficient model (3) is approximately represented as

$$Y=\mathit{X\beta }+\epsilon ,$$E17

where $X=\left(F_1,\cdots ,F_p\right)$ is $n\times \mathit{mp}$ matrix and $F_j=\left(z_{\mathit{ji}}{f}_k\left(x_i\right)\right),k=1,2,\cdots ,m,i=1,2,\cdots ,n,j=1,2,\cdots ,p$. $\beta ={\left({\beta }_1^{\prime },\cdots ,{{\beta }_p}^{\prime }\right)}^{\prime }$ is $\mathit{mp}$-dimensional parametric vector, ${\beta }_j={\left(f_j\left(t_1\right),\cdots ,f_p\left(t_m\right)\right)}^{\prime }$.

It is worth noting that under null hypothesis ${\mathrm{H}}_{31}$ defined in (15), regression model (3) is equivalent to model (17). However, this equivalence does not hold under null hypothesis ${\mathrm{H}}_{32}$ defined in (16). Null hypotheses ${\mathrm{H}}_{31}$ and ${\mathrm{H}}_{32}$ can be expressed in matrix as the following two, respectively:

$${\mathrm{H}}_{31}^{\ast }:L_1^{\prime }\beta =0,$$E18

$${\mathrm{H}}_{32}^{\ast }:L_2^{\prime }\beta =0,$$E19

where $L_1^{\prime }$ is $p\left(m-1\right)\times \mathit{mp}$ matrix.

$$L_1^{\prime }=\left(\begin{array}{ccc}{L_{01}}^{\prime }& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {L_{01}}^{\prime }\end{array}\right),L_{01}=\left[\begin{array}{ccc}1-1& 0\cdots & 00\\ \cdots \cdots & \cdots & \cdots \cdots \\ 00& 0\cdots & 1-1\end{array}\right],$$

$$L_2^{\prime }={\left(0_{\left(m-1\right)\times \left(m{j}_0-m\right)},L^{\prime },0_{\left(m-1\right)\times \left(\mathit{mp}-m{j}_0\right)}\right)}_{\left(m-1\right)\times \mathit{mp}}.$$

In the same way as the p-value in (13) is defined, p-value to test hypotheses ${\mathrm{H}}_{31}^{\ast }$ and ${\mathrm{H}}_{32}^{\ast }$ can be defined as below if the error $\epsilon$ distributes as normal distribution:

$$p_{31}\left(\widehat{\beta },S^2\right)=1-F_{p\left(m-1\right),n-\mathit{mp}}\left(\frac{\left(n-\mathit{mp}\right){\widehat{\beta }}^{\prime }{L}_1{\left(L_1^{\prime }{\left(X^{\prime }X\right)}^{-1}{L}_1\right)}^{-1}{L}_1^{\prime }\widehat{\beta }}{p\left(m-1\right)S^2}\right),$$E20

$$p_{32}\left(\widehat{\beta },S^2\right)=1-F_{m-1,n-\mathit{mp}}\left(\frac{\left(n-\mathit{mp}\right){\widehat{\beta }}^{\prime }{L}_2{\left(L_2^{\prime }{\left(X^{\prime }X\right)}^{-1}{L}_2\right)}^{-1}{L}_2^{\prime }\widehat{\beta }}{\left(m-1\right)S^2}\right).$$E21

According to the above discussion, it can be seen that $p_{31}\left(\widehat{\beta },S^2\right)$ is uniformly distributed over (0, 1) under hypothesis $H_{31}^{\ast }$. However, under null hypothesis $H_{32}^{\ast }$, varying-coefficient model (2) is not linear. Hence, there is a difference between the distribution function of $p_{32}\left(\widehat{\beta },S^2\right)$ under $H_{32}^{\ast }$ and uniform distribution. This difference has an accurate expression, which can be seen in Li et al. [37] (Theorem 3). On the other hand, $p_{31}\left(\widehat{\beta },S^2\right)$ and $p_{32}\left(\widehat{\beta },S^2\right)$ both tend to be zero in probability if null hypotheses are false when sample size tends to be infinity under some mild conditions. The corresponding proof was provided also in Li et al. [37].

3. Multi-scale method based on regression spline

For regression spline, the number of knots controls the smoothness of the estimator. The determination of knots is important and plays a large influence on the inference results. The GCV method is usually used to choose an optimal number of knots. While, but after the number of knots is given, the determination of the optimal positions of knots is difficult. Shi and Li [38] chose knots by placing an additional new knot to reduce the value of GCV, until it could not be reduced by placing any additional knots. Hence, once a knot was selected, it cannot be removed from the knot set. Mao and Zhao [39] determined the locations of knots conditioned on the number of knots $m$ first and chose $m$ later by GCV criterion. In fact, the locations of knots can be considered as parameters which can be estimated from data. This is the free-knot spline; see DiMatteo et al. [40] and Sonderegger and Hannig [41]. However, the estimation of the optimal locations is computationally intractable, and replicate knots might appear in the estimated knot vectors [42].

On the other hand, many statisticians think that the statistical inference based on one smoothing level is not reliable although it is the optimal one. Therefore, multi-scale method is developed to estimate and test nonparametric regression curves. Chaudhuri and Marron [25, 26] proposed a multi-scale method to explore the significant structures (local minima and maxima or global trend) in data, which is known as SiZer. Significant zero crossings of derivatives (SiZer) is a powerful visualization technique for exploratory data analysis. It applies a large range of smoothing parameter values to do statistical inference simultaneously and use a 2D colored map (SiZer map) to summarize all of the results inferred at different smoothing levels (scales) and locations.

In this section, a regression spline version of SiZer is proposed for exploring structures of curve and comparing multiple regression curves, respectively. The proposed SiZer employs the number of knots as smoothing parameter (scales). For the sake of simplicity, linear spline is employed first to construct SiZer, which is denoted as SiZerLS. In addition, another version of SiZer—SiZerSS—is introduced, which is proposed in Marron and Zhang [33]. In SiZerSS, smoothing spline is used to infer the monotonicity of $f\left(x\right)$, and the tuning parameter (penalty parameter) that controls the size of penalty is chosen to be as the smoothing parameter. Finally, SiZer-RS, a version of SiZer based on higher-order spline interpolation, is constructed to compare multiple regression curves at different scales and locations simultaneously.

In order to understand SiZerLS clearly, we first present an example in which SiZerLS are simulated. This example is modified from Hannig and Lee [28] with the same regression function:

$$f\left(x\right)=5+4.2\left(1+\frac{\left|x-0.3\right|}{0.03}\right)-4+5.1\left(1+\frac{\left|x–0.7\right|}{0.01}\right)-4.$$

The observations generated from model (1) with $200$ equally spaced design points from (0, 1) and $\sigma \sim N\left(\mathrm{0,0.5}\right)$ are plotted in Figure 1. Estimator ${\widehat{f}}_m\left(x\right)$ denotes the linear spline smoother obtained from (6) using $m$ equally spaced knots chosen from (0, 1). The curves of ${\widehat{f}}_m\left(x\right)$ with different values of $m$ are plotted in Figure 1 too. The simulated SiZerLS map and SiZerSS map are shown in Figure 2, respectively.

Figure 1.

200 observations and the estimated curves based on different knot sets.

Figure 2.

BYP SiZerLS and SiZerSS for detecting peaks of data.

In Figure 2, BYP SiZerLS is SiZerLS map based on multiple testing procedures, BYP, where BYP denotes the multiple testing procedure proposed in Benjamini and Yekutieli [43]. SiZerSS is the smoothing spline version of SiZer. The two SiZers are simulated under the same range of scales and nominal level 0.05. There are four colors in SiZer maps: red indicates that the estimated regression curve is significantly decreasing; blue indicates that the estimated regression curve is significantly increasing; purple indicates that the curve is neither significantly increasing nor decreasing; gray shows that there are no sufficient data for conducting reasonable statistical inference. Figure 1 preliminarily shows that SiZer maps can locate peaks well. The theoretical foundation of SiZerLS and SiZerSS will be discussed in more detail at a later stage.

3.1. Construction of SiZerLS map for exploring features of regression curve

The proposed SiZerLS map will be constructed on the basis of the p-values with multiple testing adjustment. The p-value for testing the monotonicity of the smoothed curve is defined first based on linear spline approximation and fiducial method in the same way as p-values in Section 2. Consequently, multiple testing adjustment is discussed detailedly to control the row-wise false discovery rate (FDR) of SiZerLS.

In the view of SiZer, all of the useful information is included in the smoothed curve, which is defined below. Suppose we have observations ${\left\{x_i,y_i\right\}}_{i=1}^n$ from regression model (1). By linear spline estimation, estimator ${\widehat{f}}_m\left(x\right)$ can be obtained:

$${\widehat{f}}_m\left(x\right)=\mathit{g}{\left(x\right)}^{\prime }{\left(G^{T}G\right)}^{-1}{G}^{T}Y,$$E22

where $\mathit{g}\left(x\right)\mathit{=}{\left(g_1\left(x\right),g_2\left(x\right),\cdots ,g_m\left(x\right)\right)}^{\prime }$; $g_j\left(x\right),j=1,\cdots ,m$ are the basis functions defined in (7) on the basis of $m$ knots; and $G$ is the matrix defined in Section 2. The smoothed curve at smoothing level $m$ is denoted as.

$$f_m\left(x\right)=\mathrm{E}\left({\widehat{f}}_m\left(x\right)\right)=\mathit{g}{\left(x\right)}^{\prime }{\left(G^{T}G\right)}^{-1}{G}^T\mathit{f},$$

where $\mathit{f}\mathit{=}{\left\{f\left(x_1\right),f\left(x_2\right),\cdots ,f\left(x_n\right)\right\}}^{\prime }$. SiZer focuses on $f_m\left(x\right).$ Its monotonicity is determined totally by ${\left(G^{T}G\right)}^{-1}{G}^T\mathit{f}$. Hence, it is enough to test the following $m-1$ pairs of null hypotheses:

$$H_{\mathit{Ik}}=f_m\left(t_k\right)=e_k^{\prime }{\left(G^{\prime }G\right)}^{-1}{G}^{\prime }\mathit{f}\le e_{k+1}^{\prime }{\left(G^{\prime }G\right)}^{-1}{G}^{\prime }\mathit{f}=f_m\left(t_{k+1}\right)\text{(and)}$$

$$H_{\mathit{Dk}}=f_m\left(t_k\right)=e_k^{\prime }{\left(G^{\prime }G\right)}^{-1}{G}^{\prime }\mathit{f}\ge e_{k+1}^{\prime }{\left(G^{\prime }G\right)}^{-1}{G}^{\prime }\mathit{f}=f_m\left(t_{k+1}\right),k=1,2,\cdots ,m-1,$$E23

where $e_k$ is an $m$-dimensional column vector having 1 in the $k$th entry and zero elsewhere. Let $b$ denote ${\left(G^{\prime }G\right)}^{-1}{G}^{\prime }\mathit{f}$. Then, $H_{\mathit{Ik}}$ and $H_{\mathit{Dk}}$ can be written as

$$H_{\mathit{Ik}}^{\ast }=L_{k}b\le 0,k=1,2,\cdots ,m-1;H_{\mathit{Dk}}^{\ast }=L_{k}b\ge 0,k=1,2,\cdots ,m-1,$$E24

where $L_k\triangleq (e_k^{\prime }-e_{k+1}^{\prime }$). The p-values to test hypotheses in (24) under linear model $Y=\mathit{Gb}+\epsilon$ can be defined using pivotal quantity about $b$. This pivotal quantity is $R\left(E;\widehat{\beta },S^2\right)$, which is defined in (11). The p-value for testing $H_{\mathit{Ik}}^{\ast }$ is the fiducial probability that null hypothesis holds:

$$P_{\mathit{Ik}}^{\ast }\left(\widehat{\beta },S\right)=P\left\{L_{k}R\left(E;\widehat{\beta },S\right)\le 0\right\}=P\left\{L_k\widehat{\beta }-\frac{S}{E_2}{\left(G^{\prime }G\right)}^{-\frac{1}{2}}{E}_1\le 0\right\}$$

$$=P\left\{\frac{\sqrt{n-m}{L}_k{\left(G^{\prime }G\right)}^{-1}{G}^{\prime }{E}_1}{{\left(L_k{\left(G^{\prime }G\right)}^{-1}{L}_k^{\prime }\right)}^{\frac{1}{2}}{E}_2}\ge \frac{\sqrt{n-m}\widehat{\beta }}{S{\left(L_k{\left(G^{\prime }G\right)}^{-1}{L}_k^{\prime }\right)}^{\frac{1}{2}}}\right\},$$E25

where the subscript $\mathit{Ik}$ of $P_{\mathit{Ik}}^{\ast }$ represents the interval ($t_k,t_{k+1}$) in which we test monotonicity and $m$ represents the number of knots used in linear interpolation. In addition, p-value $P_{\mathit{Dk}}^{\ast }\left(\widehat{\beta },S\right)$ for testing $H_{\mathit{Dk}}^{\ast }$ satisfies equation $P_{\mathit{Ik}}^{\ast }\left(\widehat{\beta },S\right)+P_{\mathit{Dk}}^{\ast }\left(\widehat{\beta },S\right)=1$.

It is worth noting that p-value $P_{\mathit{Ik}}^{\ast }\left(\widehat{\beta },S\right)$ is uniformly distributed on (0,1) if all of the hypotheses $H_{\mathit{Ik}},H_{\mathit{Dk}},k=1,2,\cdots ,m-1$ are true (regression function is a constant). In applications, p-value $P_{\mathit{Ik}}^{\ast }\left(\widehat{\beta },S\right)$ for testing $H_{\mathit{Ik}}$ can be approximated as below when $n\to \infty .$ This approximation is reasonable (see Theorem 1 in [44]):

$$P_{\mathit{Ik},m}\left(\widehat{\beta },S\right)\triangleq 1-\mathrm{\Phi }\left(\frac{\sqrt{n-m}{L}_k\widehat{\beta }}{S{\left(L_k{\left(G^{\prime }G\right)}^{-1}{L}_k\right)}^{1/2}}\right).$$E26

The proposed SiZerLS map will be constructed on the basis of the above p-values with multiple testing adjustment. In fact, SiZer is a visual method for exploratory data analysis, and it focuses on exploring features that really exist in data instead of testing whether some assumed features are statistically significant in a strict way. FDR is the expected proportion of the false positives among all discoveries, and FDR can be either permissive or conservative according to the number of hypotheses. Considering that different numbers of hypotheses need to be tested for SiZerLS with respect to various smoothing parameters, the multiple testing adjustment to control FDR would be better if used to improve the exploratory property of SiZer. Hence, the well-known multiple testing procedure which was proposed in Benjamini and Yekutieli [43] (denoted as BYP) is applied to control the row-wise FDR of SiZerLS. The BYP was proved to control FDR under α for any dependent test statistics.

3.1.1. Benjamin-Yekutieli procedure to control FDR (BYP)

Suppose that we have obtained p-values $P_{\mathit{Ik},m}\left(\widehat{\beta },S\right)$ for testing hypotheses $H_{\mathit{Ik}}$ in (23), $k=1,2,\cdots ,m-1$:

1. Order p-values $P_{\mathit{Ik},m}^{\ast }$ and get the ordered p-values $P_{I\left(1\right),m},P_{I\left(2\right),m},\cdots ,P_{I\left(m-1\right),m}$.

2. For a given p-value α, find the largest $i$ for $k=1,2,\cdots ,m-1$ for which $P_{I\left(i\right),m}\le \frac{\mathit{k\alpha }}{\left(m-1\right){\sum }_{j=1}^{m-1}\frac{1}{j}}$ and reject all $H_{I\left(k\right),m}$ for $k=1,2,\cdots ,m-1$.

The detailed steps to construct SiZerLS with BYP adjustment are given below:

Step 1. Construct 2D grid map. Without loss of generality, we assume that designed points $x_i,i=1,2,\cdots ,n$ are chosen from [0, 1]. Then the 2D map is a rectangular area [0, 1; ${\text{log}}_{10}\left(1/m_{\max },\right),{\text{log}}_{10}\left(1/m_{\min }\right)]$; see BYP SiZerLS displayed in Figure 2. The value of $m$ is determined by the following rule: $m=\text{round}\left(1/{10}^l\right)$, where function $\text{round}$ (∙) is the nearest integer function and $l$ takes equally spaced values from interval $\left[{\text{log}}_{10}\left(1/m_{\min }\right),{\text{log}}_{10}\left(1/m_{\max }\right)\right]$. For a given $m$, abscissa $x$ takes values at the corresponding knots ${\mathbf{T}}_m=\left\{t_1,t_2,\cdots ,t_m\right\}$. On the basis of different values of $m$ and ${\mathbf{T}}_m$, the 2D map is divided into many pixels.

Step 2. Calculate p-values for each pixel. Each pixel in the 2D map constructed in step 1 is determined by two adjacent knots and a determined $m$. For pixel $\left(t_k,t_{k+1},m=m_0\right)$, we calculate p-value $P_{\mathit{Ik},m_0}$ and $P_{\mathit{Dk},m_0}$ for testing hypotheses $H_{\mathit{Ik},m_0}^{\ast }$ and $H_{\mathit{Dk},m_0}^{\ast }$, respectively, with $m_0$ knots.

Step 3. Multiple testing adjustment. For a given value $m=m_0$, carry out multiple testing procedure BYP using p-values $P_{\mathit{Ik},m_0}$ ($P_{\mathit{Dk},m_0}$), $k=1,2,\cdots ,m_0$, obtained from step 2 to test the fowling family of hypotheses simultaneously:

$$\left\{H_{I1,m_0}^{\ast },H_{I2,m_0}^{\ast },\cdots ,H_{I{m}_0-1,m_0}^{\ast }\right\}\left(H_{D1,m_0}^{\ast },H_{D2,m_0}^{\ast },\cdots ,H_{D{m}_0-1,m_0}^{\ast }\right).$$

Step 4. Color pixels. According to the multiple testing results at smoothing level $m_0$ if $H_{\mathit{Ik}}^{\ast }$ is rejected and $H_{\mathit{Dk}}^{\ast }$ is accepted, pixel $\left(t_k,t_{k+1},m=m_0\right)$ is colored red to indicate significant decreasing. On the contrary, if $H_{\mathit{Ik},,m_0}^{\ast }$ is accepted and $H_{\mathit{Dk},m_0}^{\ast }$ rejected, pixel $\left(t_k,t_{k+1},m=m_0\right)$ is colored blue to show significant increasing; purple is used for no significant trend in other cases.

In SiZer map, gray indicates that no sufficient data can be used to test the monotonicity of regression function at point $x$ with $m$ knots. Such sufficiency is quantified as effective sample size (ESS). Noting that the number of nonzero elements in the $k$th column of $G$ has a demonstrable effect on the inference in interval $\left(t_k,t_{k+1}\right),$ and it is determined directly by how many observations are included in $\left(t_k,t_{k+1}\right)$, we define $\mathit{ESS}\left(t_k,m\right)$ as.

$${\left(\mathit{ESS}\left(t_1,m\right),\mathit{ESS}\left(t_2,m\right),\cdots ,\mathit{ESS}\left(t_m,m\right)\right)}^{\prime }\triangleq G^{\prime }G{\left(1,1,\cdots ,1\right)}^{\prime }.$$

In SiZerLS map, pixel $\left(t_k,t_{k+1},m=m_0\right)$ would be colored gray if.

$$\text{min}\left(\mathit{ESS}\left(t_k,m_0\right),\mathit{ESS}\left(t_{k+1},m_0\right)\right)<5.$$

In order to avoid selecting knots, $m$ equally spaced knots or equal $x$-quantiles are used in interpolation. The smoothing level of regression spline estimate is controlled by $m$ together with the positions of knots. The level of smoothness should be reduced to detect some local fine feature; however, the total number of knots should be limited to avoid excessive under-smoothing in a wide range. In applications of SiZerLS, the range of scales is recommended to include the coarsest smoothing level, $m=2$, and the finest smoothing level, ${\mathrm{avg}}_{x\in {{\mathbf{T}}_m}_{\max }}\mathit{ESS}\left(x,m_{\max }\right)$ < 5.

3.2. Construction of SiZerSS map for exploring features of regression curve

SiZerSS given in Marron and Zhang [33] employed smoothing spline to construct SiZer map for nonparametric model (1). Given ${\left\{x_i,y_i\right\}}_{i=1}^n$ and a smoothing parameter $\lambda$, the smoothing spline estimator is the function ${\widehat{f}}_{\lambda }$ that minimizes the regularization criterion over function $f$:

$$\sum _{i=1}^n{\omega }_i{\left[y_i-f\left(x_i\right)\right]}^2+\lambda \int {\left[f^{\prime \prime }\left(x\right)\right]}^2\mathrm{d}x.$$E27

By simple calculation, we can get the estimator vector:

$${\widehat{\mathbf{f}}}_{\mathbf{\lambda }}\mathbf{=}\left({\widehat{f}}_{\lambda }\left(x_1\right),{\widehat{f}}_{\lambda }\left(x_2\right),\mathbf{\cdots },{\widehat{f}}_{\lambda }\left(x_n\right)\right)\mathbf{=}{\left(W+\mathit{\lambda K}\right)}^{-1}\mathit{WY}=A_{\lambda }Y,$$E28

where weight matrix $W=\text{diag}\left({\omega }_1,{\omega }_2,\cdots ,{\omega }_n\right)$ and the hat matrix $A_{\lambda }={\left(W+\mathit{\lambda K}\right)}^{-1}W$.

In order to construct SiZerSS, the derivative of $f$ at any point $x$ needs to be estimated along with its variance. Let $s_i=x_{i+1}-x_i$ and $n\times \left(n-1\right)$ matrix $Q=\left\{q_{\mathit{ij}}\right\},i=1,2,\cdots ,n,j=2,\cdots ,n-1$, where $q_{j-1,j}=s_{j-1}^{-1},q_{\mathit{jj}}=-s_{j-1}^{-1}-s_j^{-1},q_{j+1,j}=s_j^{-1}$, and $q_{i,j}=0$ for $\left|i-j\right|\ge 2.$ Let $\left({\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_n\right)=\left(f^{\prime \prime }\left(x_1\right),f^{\prime \prime }\left(x_2\right)\cdots ,f^{\prime \prime }\left(x_n\right)\right)$. By the definition of natural cubic spline, $f^{“}\left(x_1\right)=f^{“}\left(x_n\right)=0$. Let $\gamma ={\left({\gamma }_2,\cdots ,{\gamma }_{n-1}\right)}^{\prime }$. According to Theorem 2.1 of Green and Silverman [45], the vectors $\mathbf{f}$ and $\gamma$ specify a natural cubic spline $f$ if and only if $Q^{\prime }\mathbf{f}=\mathit{R\gamma }$,

where $R$ is a ($n-2)\times \left(n-2\right)$ symmetric matrix with elements $r_{\mathit{ij}},i=2,\cdots ,n-1,j=2,\cdots ,n-1,$ which is given by $r_{\mathit{ii}}=\frac{1}{3}\left(s_{i-1}+s_i\right),r_{i,i+1}=r_{i+1,i}=\frac{1}{6}{s}_i$ and $r_{\mathit{ij}}=0$ for $\left|i-j\right|\ge 2$. The estimator $\widehat{\gamma }$ can be obtained from equation $\left(R+\lambda Q^{\prime }Q\right)\gamma =Q^{\prime }Y$. Then estimator $\widehat{f}\left(x\right)$ and ${\widehat{f}}^{\prime }\left(x\right)$ can be written as a linear combination of $\widehat{\mathbf{f}}$ and $\widehat{\gamma }$. Let $h_i\left(x\right)=x-x_i,i=1,2,\cdots ,n.$ When $x<x_{1.}$

$${\widehat{f}}_{\lambda }\left(x\right)={\widehat{f}}_{\lambda }\left(x_1\right)+h_1\left(x\right)\left\{\frac{{\widehat{f}}_{\lambda }\left(x_2\right)-{\widehat{f}}_{\lambda }\left(x_1\right)}{s_1}-\frac{s_1}{6}{\widehat{\gamma }}_2\right\},{\widehat{f}}^{\prime }\left(x\right)=\frac{{\widehat{f}}_{\lambda }\left(x_2\right)-{\widehat{f}}_{\lambda }\left(x_1\right)}{s_1}-\frac{s_1}{6}{\widehat{\gamma }}_2.$$

When $x_i\le x\le x_{i+1}$, let ${\delta }_i\left(x\right)=\left(1+\frac{h_i\left(x\right)}{s_i}\right){\widehat{\gamma }}_{i+1}+\left(1-\frac{h_{i+1}\left(x\right)}{h_i}\right){\widehat{\gamma }}_i$ for $i=1,2,\cdots ,n$,

$${\widehat{f}}_{\lambda }\left(x\right)=\frac{h_i\left(x\right){\widehat{f}}_{\lambda }\left(x_{i+1}\right)-h_{i+1}\left(x\right){\widehat{f}}_{\lambda }\left(x_i\right)}{s_i}+\frac{h_i\left(x\right)h_{i+1}\left(x\right){\delta }_i\left(x\right)}{6},$$

$${{\widehat{f}}_{\lambda }}^{\prime }\left(x\right)=\frac{{\widehat{f}}_{\lambda }\left(x_{i+1}\right)-{\widehat{f}}_{\lambda }\left(x_i\right)}{s_i}+\frac{h_i\left(x\right)h_{i+1}\left(x\right)\left({\widehat{\gamma }}_{i+1}-{\widehat{\gamma }}_i\right)}{6s_i}+\frac{h_i\left(x\right)+h_{i+1}\left(x\right)}{6}{\delta }_i\left(x\right).$$

(When) $x>x_n$

$${\widehat{f}}_{\lambda }\left(x\right)={\widehat{f}}_{\lambda }\left(x_n\right)+\frac{h_n\left(x\right)}{6}\left\{\frac{{\widehat{f}}_{\lambda }\left(x_n\right)-{\widehat{f}}_{\lambda }\left(x_{n-1}\right)}{s_{n-1}}+s_{n-1}{\widehat{\gamma }}_{n-1}\right\},$$

$${{\widehat{f}}_{\lambda }}^{\prime }\left(x\right)=\frac{1}{6}\left\{\frac{{\widehat{f}}_{\lambda }\left(x_n\right)-{\widehat{f}}_{\lambda }\left(x_{n-1}\right)}{s_{n-1}}+s_{n-1}{\widehat{\gamma }}_{n-1}\right\}.$$

The variance of ${{\widehat{f}}_{\lambda }}^{\prime }\left(x\right)$ can be calculated easily if the estimator of ${\sigma }^2,$ the variance of the error in model (1), is obtained. ${\sigma }^2$ can be estimated by the sum of squared residuals $\sum {\left(y_i-{\widehat{f}}_{\lambda }\left(x_i\right)\right)}^2$. If ${\sigma }^2$ is a function of $x$, ${\sigma }^2\left(x\right)$ can be estimated by ${\left(y_i-{\widehat{f}}_{\lambda }\left(x\right)\right)}^2$. The confidence interval of ${f_{\lambda }}^{\prime }\left(x\right)$ are of the form:

$${{\widehat{f}}_{\lambda }}^{\prime }\left(x\right)\pm q.\widehat{\mathrm{SD}}\left({{\widehat{f}}_{\lambda }}^{\prime }\left(x\right)\right),$$E29

where $q$ is based on the nominal level. For details, see Section 3 of Chaudhuri and Marron [25].

SiZerSS can be constructed as SiZerLS. For different values of $x$, if interval (29) contains zero, pixel $\left(x,\lambda \right)$ is colored purple; if confidence interval is on the right side of zero, blue is used to indicate increasing; otherwise, red is used to imply decreasing. Gray is used to indicate that there is no sufficient data to do reliable inference. The sufficiency can be found in Chaudhuri and Marron [25].

The simulated SiZerLS and SiZerSS maps are displayed in Figure 2, where the red and blue regions locate the bumps of regression curve accurately. This simulation illustrates the good behavior of SiZerLS and SiZerSS in exploring features in data.

3.3. Construction of SiZer-RS map for comparing multiple regression curves

The comparison of two or more populations is a common problem and is of great practical interest in statistics. In this subsection, comparison of multiple regression curves in a general regression setting is developed based on regression spline. Suppose we have $n=\sum _{i=1}^k{n}_i$ independent observations from the following $k$ regression models:

$$y_{\mathit{ij}}=f_i\left(x_{\mathit{ij}}\right)+{\sigma }_i\left(x_{\mathit{ij}}\right){\epsilon }_{\mathit{ij}},i=1,2,\cdots ,k,j=1,2,\cdots ,n_i,$$E30

where $x_{\mathit{ij}}$ s are covariates, the errors ${\epsilon }_{\mathit{ij}}\sim N\left(0,1\right)$ s are independent and identically distributed errors, $f_i\left(·\right)$ is the regression function, and ${\sigma }_i^2\left(·\right)$ is the conditional variance function of the $i$th population. We are concerned about whether the $k$ populations in model (30) are equal; if not, what is the difference? To this end, a multi-scale method, SiZer-RS, based on regression spline is proposed to compare $f_i\left(·\right)$ across multiple scales and locations.

As described in Park and Kang [32], the choice of smoothing parameter is also important for comparing regression curves. They developed SiZer for the comparison of regression curves based on local linear smoother. SiZer map for comparing regression curves is a 2D color map, which consists of a large number of pixels. Each pixel is indexed by a scale (smoothing parameter) and a location; the color of a pixel indicates the result for testing the equality of two or more multiple regression curves at the corresponding location and scale. SiZer provides us with more information about the locations of the differences among the regression curves if they do exist. Park et al. [46] developed an ANOVA-type test statistic and conducted it in scale space for testing the equality of more than two regression curves.

The works mentioned above are kernel-based method. Besides it, regression spline is an important smoothing device and is used widely in applications. For a given smoothing parameter $m$ (the number of knots used in regression spline), the p-value for testing the equality of $k$ regression curves at point $x$ is established. Consequently, SiZer-RS is constructed in the same way as SiZerLS for comparing multiple retrogression curves based on higher-order spline interpolation.

For a given smoothing parameter $m$ (the number of knots used in regression spline), the smoothed curve is defined as $f_{i,m}\left(x\right)=\mathrm{E}({\widehat{f}}_{i,m}$ (x)), where ${\widehat{f}}_{i,m}\left(x\right)$ is the regression spline estimator. SiZer-RS for comparing multiple regression curves is based on the testing results for testing null hypothesis:

$$H_{m,x}:f_{1,m}\left(x\right)=f_{2,m}\left(x\right)=\cdots =f_{k,m}\left(x\right),$$E31

at point $x$ with smoothing parameter $m$. Without loss of generality, we still suppose that the explanatory variable $x$ takes value from [0, 1]. On the basis of a knot set ${\mathbf{T}}_m=\left\{0=t_1<t_2,\cdots ,<t_m=1\right\}$, we have the approximation:

$$f_i\left(x\right)\approx \sum _{s=1}^{m+q-1}{\beta }_{i,s}{g}_{m,s}\left(x\right)\triangleq N_m{\left(x\right)}^{\prime }{\beta }_i^m,$$E32

where ${\beta }_i^m={\left({\beta }_{i,1},{\beta }_{i,2},\cdots ,{\beta }_{i,m+p-1}\right)}^{\prime }.$ The estimator of $f_i\left(x\right)$ at smoothing level $m$ can be obtained ${\widehat{f}}_{i,m}\left(x\right)=N_m{\left(x\right)}^{\prime }{\widehat{\beta }}_i^m$, in which, $N_m\left(x\right)=\left\{g_{m,s}\left(x\right),s=1,2,\cdots ,m+q-1\right\}.$ If $q=3$, $N_m{\left(x\right)}^{\prime }$ is defined below:

$$N_m^l\left(x\right)=\left(t_l-t_{l-4}\right)\left[t_{l-4},t_{l-3},t_{l-2},t_{l-1},t_l\right]{\left(t-x\right)}_{+}^3,l=2,3,\cdots ,m+3,$$

where $t_l=t_{\text{min}\left(\text{max}\left(l,1\right),m\right)}$ for $l=-2,-1,\cdots ,m+3$:

$${\left(t-x\right)}_{+}^3=\left\{\begin{array}{c}{\left(t-x\right)}^3,t>x\\ 0,t\le x\end{array}\right..$$

For a function $g\left(·\right),\left[t_{l-4},t_{l-3},t_{l-2},t_{l-1},t_l\right]g\left(·\right)$ denotes the fourth-order divided difference of $g\left(·\right)$, that is:

$$\left\{\begin{array}{l}\left[t_1,t_2\right]g=g^{\prime }\left(t\right),\text{if}{t}_1=t_2=t\hfill \\ \left[t_1,t_2\right]g=\frac{g\left(t_2\right)-g\left(t_1\right)}{t_2-t_1}\text{otherwise},\\ \left[t_1,t_2,\cdots ,t_k\right]g=g^{\left(k-1\right)}\left(t\right),\mathit{\text{if}}{t}_1=\cdots =t_k\\ \left[t_1,t_2,\cdots ,t_k\right]g=\frac{\left[t_2,t_3,\cdots ,t_k\right]g-\left[t_1,t_2,\cdots ,t_{k-1}\right]g}{t_k-t_1},\text{otherwise}.\end{array}\right.$$

Then model (31) can be approximately written as the following linear regression model:

$$Y_i=G_i^m{\beta }_i^m+{\mathrm{\Sigma }}_i{E}_i,$$E33

where

$$Y_i={\left(y_{i1},y_{i2},\cdots ,y_{i{n}_i}\right)}^{\prime },G_i^m={\left(N_m^l\left(x_i\right)\right)}_{n_i\times \left(m+2\right)},{\mathrm{\Sigma }}_i=\text{diag}\left\{{\sigma }_i\left(x_{\mathit{ij}}\right)\right\},E_i={\left({\epsilon }_{i1},{\epsilon }_{i2},\cdots ,{\epsilon }_{i{n}_i}\right)}^{\prime }.$$

At first, we suppose ${\mathrm{\Sigma }}_i$ is known and then replace it by its available estimator.

From regression model (33), we can get the estimator ${\widehat{\beta }}_i^m={\left(G_i^{m\prime }{\mathrm{\Sigma }}_i^{-1}{G}_i^m\right)}^{-1}{G_i^m}^{\prime }{\mathrm{\Sigma }}_i^{-1}{Y}_i.$ Let ${\mathrm{b}}_i^m$ denote the expectation of ${\widehat{\beta }}_i^m$:

$${\mathrm{b}}_i^m=\mathrm{E}\left({\widehat{\beta }}_i^m\right)={\left({G_i^m}^{\prime }{\mathrm{\Sigma }}_i^{-1}{G}_i^m\right)}^{-1}{G_i^m}^{\prime }{\mathrm{\Sigma }}_i^{-1}{\mathbf{f}}_i,$$

where ${\mathbf{f}}_i={\left(f_i\left(x_{i1}\right),\cdots ,f_i\left(x_{i{n}_i}\right)\right)}^{\prime }$. Therefore, the smoothed curve

$$f_{i,m}\left(x\right)=\mathrm{E}\left({\widehat{f}}_{i,m}\left(x\right)\right)=\mathrm{E}\left[N_m{\left(x\right)}^{\prime }{\left(G_i^{m\prime }{\mathrm{\Sigma }}_i^{-1}{G}_i^m\right)}^{-1}{G}_i^{m\prime }{\mathrm{\Sigma }}_i^{-1}{Y}_i\right]=N_m{\left(x\right)}^{\prime }{\mathrm{b}}_i^m.$$E34

Denote ${\mathbf{b}}^m={\left({\mathrm{b}}_1^{m\prime },{\mathrm{b}}_2^{m\prime },\cdots ,{\mathrm{b}}_k^{m\prime }\right)}^{\prime },$ and correspondingly, denote its estimator as ${\widehat{\beta }}^m={\left({\beta }_1^{m\prime },{\beta }_2^{m\prime },\cdots ,{\beta }_k^{m\prime }\right)}^{\prime }$. Hypothesis $H_{m,x}$ can be presented as

$$H_{m,x}:L_m\left(x\right){\mathbf{b}}^m=0_{k-1},$$E35

where

$$L_m\left(x\right)=\left[\begin{array}{lllll}{N}_m\left(x\right)& N_m\left(x\right)& N_m\left(x\right)& & N_m\left(x\right)\\ -N_m\left(x\right)& 0& 0& \cdots & 0\\ 0& -N_m\left(x\right)& 0& & 0\\ \cdots & \cdots & \cdots & & \cdots \\ 0& 0& 0& \cdots & 0\\ 0& 0& 0& & -N_m\left(x\right)\end{array}\right]$$

is a $\left(k-1\right)\times k\left(m+q-1\right)$ matrix.

The p-value for testing hypothesis $H_{m,x}$ in (35) can be defined as

$$p_{m,x}\left({\widehat{\beta }}_i^m,{\widehat{\mathrm{\Sigma }}}_m\right)=P\{T^m{\left(x\right)}^{\prime }{L}_m{\left(x\right)}^{\prime }{\left[L_m\left(x\right){\widehat{\mathrm{\Sigma }}}_m{L}_m{\left(x\right)}^{\prime }\right]}^{-1}{L}_m\left(x\right)T^m\left(x\right)$$

$$\ge {\widehat{\beta }}_i^{m\prime }{L}_m{\left(x\right)}^{\prime }{\left[L_m\left(x\right){\widehat{\mathrm{\Sigma }}}_m{L}_m{\left(x\right)}^{\prime }\right]}^{-1}{L}_m\left(x\right){\widehat{\beta }}_i^m\},$$E36

where $T^m\left(x\right)\triangleq {\left\{{\left(G_i^{m\prime }{\widehat{\mathrm{\Sigma }}}_{i,m}^{-1}{G}_i^m\right)}^{-1}{G}_i^{m\prime }{\widehat{\mathrm{\Sigma }}}_{i,m}^{-\frac{1}{2}}{E}_i\right\}}^{\prime },i=1,2,\cdots ,k;{\widehat{\mathrm{\Sigma }}}_{i,m}=\text{diag}\left\{{\widehat{\sigma }}_i\left(x_{\mathit{ij}}\right),j=1,2,\cdots ,n_i\right\}$ is an estimator of the variance matrix of the $i$th regression model and

$${\widehat{\mathrm{\Sigma }}}_m=\text{diag}\left\{{\left(G_i^{m\prime }{\widehat{\mathrm{\Sigma }}}_{i,m}^{-1}{G}_i^m\right)}^{-1},i=1,2,\cdots ,k\right\}$$

is an estimator of the variance matrix of $T^m\left(x\right)$ given ${\widehat{\beta }}_i^m,{\widehat{\sigma }}_i^{m2},i=1,2,\cdots ,k.$ The estimator of ${\sigma }_i\left(x_{\mathit{ij}}\right)$ can be found in Li and Xu [36], where the smoothing parameter $,m_p$, can be used as a pilot smoothing parameter, which is different from $m$ used in ${\widehat{f}}_{i,m}\left(x\right)$. SiZer-RS map can be constructed based on different values of $m_p$, which represents the different trade-offs between the structure of regression curve and errors.

The two SiZer maps given in Figure 4 are constructed using the data plotted in Figure 3 to compare three regression curves $f_1\left(x\right)=f_x\left(x\right)=0,f_3\left(x\right)=0.5\text{sin}\left(2\mathit{\pi x}\right)$. Since the variance of errors is a constant, it can be estimated by the sum of squares of residues. In this case, pilot smoothing parameter is avoided [47, 48]. The two blue regions in Figure 4 clearly show their difference across interval (0, 1). The gray color indicates that there is no sufficient data that can be used to get credible testing results at $x$ and nearby. The sufficiency is quantized as $\mathit{ESS}\left(x,m\right)$ for SiZer-RS, and pixel $\left(x,m\right)$ is colored gray if $\mathit{ESS}\left(x,m\right)<5$:

$$\mathit{ESS}\left(x,m\right)\triangleq \underset{i=1,2,\cdots ,k}{\text{min}}\left\{N_m\left(x\right)G_i^{m\prime }{G}_i^m{\left(1,1,\cdots ,1\right)}^{\prime }\right\}.$$

Figure 3.

200 observations.

Figure 4.

SiZer-RS map.

Figure 4 shows that SiZer-RS map can explore the differences between regression curves accurately.

It is worth noting that, for SiZer-RS map, the coarsest smoothing level should be $m=q+1$ to ensure the effectiveness of the $q$th regression spline and the finest smoothing level is recommend to be the one such that $\underset{x\in \left[x_1,x_2,\cdots ,x_g\right]}{\mathrm{avg}}\mathit{ESS}\left(x,m\right)<5$, where $x_1,x_2,\cdots ,x_g$ are points at which hypothesis $H_{m,x}$ is tested and pixels are produced by combing different values of $m$. In applications, a wide range of values of $m_p$ can be used to generate a family of SiZer-RS maps. Particularly, $m_p$ and $m$ can both be used as smoothing parameters simultaneously to construct a 3D SiZer-RS map [47, 48].

4. Conclusion

This chapter introduces regression spline method for testing the parametric form of nonparametric regression function in nonparametric, partial linear, and varying-coefficient models, respectively. The corresponded p-values are established based on fiducial method and spline interpolation. The test procedures on the basis of the proposed p-value are accurate in some cases and are consistent under some mild conditions, which means that the p-value tends to be zero when null hypothesis is false as sample size and the number of knots used in spline interpolation tend to be infinity. Hence, the proposed test procedures are performed well especially in small sample size case.

The spline-based method frequently used smoothing method, which can be used easily with other statistical methods. When using the spline-based method, the smoothing level is controlled by the number of knots and their positions. In order to sidestep the determination of knots and obtain more reliable results, multi-scale smoothing methods are proposed based on spline regression to infer structures of regression function. The multi-scale method is a visual method to do inference at different locations and smoothing levels. In addition, the smoothing spline version of multi-scale method is also introduced. The proposed multi-scale method can also be used for comparing multiple regression curves. Some real data examples illustrate the practicability of the proposed multi-scale method.

The MATLAB code of SiZerLL and other versions of SiZer based on kernel smoother is available from the homepage of Professor Marron JS; the MATLAB code of SiZerLS can be downloaded from the following website:

http://www.tandfonline.com/doi/suppl/10.1080/10618600.2014.1001069?scroll=top.