Determining confidence interval and asymptotic distribution for parameters of multiresponse semiparametric regression model using smoothing spline estimator

3.1 Regression function estimation

We may present the MSR model (1) as follows:

We can rewrite model (2) as follows:

Suppose ${\widehat{\mathit{\beta }}}_k$ is the true WLS estimate of ${\mathit{\beta }}_k$ . Hence, we can express model (3) as follows:

Next, let ${\mathbf{y}}^{\mathbf{\ast }}={({\mathbf{y}}_{1i}^{\ast },{\mathbf{y}}_{2i}^{\ast },\cdots ,{\mathbf{y}}_{\mathit{pi}}^{\ast })}^T$ ; $\mathbf{g}={({\mathbf{g}}_1,{\mathbf{g}}_2,\cdots ,{\mathbf{g}}_p)}^T$ ; $\epsilon ={({\epsilon }_{1i},{\epsilon }_{2i},\cdots ,{\epsilon }_{\mathit{pi}})}^T$ ; and $\mathbf{t}={({\mathbf{t}}_{1i},{\mathbf{t}}_{2i},\cdots ,{\mathbf{t}}_{\mathit{pi}})}^T$ where ${\mathbf{y}}_{\mathit{ki}}^{\ast }={(y_{k1}^{\ast },y_{k2}^{\ast },\cdots ,y_{k{n}_k}^{\ast })}^T$ ; ${\epsilon }_{\mathit{ki}}={({\epsilon }_{\mathrm{k}1},{\epsilon }_{k2},\cdots ,{\epsilon }_{k{n}_k})}^T$ ;

${\mathbf{t}}_{\mathit{ki}}={(t_{k1},t_{k2},\cdots ,t_{k{n}_k})}^T$ ; and $k=1,2,\cdots ,p$ .

Hence, we can present the MSR model (4) in the following matrix equation:

where $E\left(\epsilon \right)=0$ , $Cov\left(\epsilon \right)={\mathbf{W}}^{-1}$ (namely).

The smoothing spline estimator of function $\mathbf{g}$ in model (6) can be determined by solving the PWLS: $$\underset{{\mathbf{g}}_1,\cdots {,\mathbf{g}}_p\in W_2^m}{\mathrm{Min}}\left\{N^{-1}\right[{\left({\mathbf{y}}_1^{\ast }-{\mathbf{g}}_1\right)}^T{\mathbf{W}}_1\left({\mathbf{y}}_1^{\ast }-{\mathbf{g}}_1\right)+{{\left({\mathbf{y}}_2^{\ast }-{\mathbf{g}}_2\right)}^T\mathbf{W}}_2\left({\mathbf{y}}_2^{\ast }-{\mathbf{g}}_2\right)+\cdots +$$ $${\left({\mathbf{y}}_p^{\ast }-{\mathbf{g}}_p\right)}^T{\mathbf{W}}_p\left({\mathbf{y}}_p^{\ast }-{\mathbf{g}}_p\right)]+{\lambda }_1\underset{a_1}{\stackrel{b_1}{\int }}{\left({\mathbf{g}}_1^{\left(m\right)}\left({\mathbf{t}}_1\right)\right)}^{2}d{\mathbf{t}}_1+$$

where $N={\sum }_{k=1}^p{n}_k$ ; ${\mathbf{W}}_1,{\mathbf{W}}_2,\cdots ,{\mathbf{W}}_p$ are weight matrices that are inverse of covariance matrix, $\mathbf{g}\in W_2^{\mathit{m}}[a,b]$ , and ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_p$ are smoothing parameters that set the balance between good fit and smoothness of estimation.

Based on Eq.(6), it is easy to show that the covariance matrix of random errors in MSR model (1) is:

where ${\mathbf{W}}_k^{-1}=\left(\begin{array}{ccc}{\sigma }_{k1}^2& {\sigma }_{k\left(1,2\right)}& \begin{array}{cc}\cdots & {\sigma }_{k(1,n_1)}\end{array}\\ {\sigma }_{k\left(2,1\right)}& {\sigma }_{k2}^2& \begin{array}{cc}\cdots & {\sigma }_{k(2,n_2)}\end{array}\\ \begin{array}{c}⋮\\ {\sigma }_{k(n_k,1)}\end{array}& \begin{array}{c}⋮\\ {\sigma }_{k(n_k,2)}\end{array}& \begin{array}{cc}\ddots & ⋮\\ \cdots & {\sigma }_{k{n}_k}^2\end{array}\end{array}\right)$ , and $k=1,2,\cdots ,p$ .

Solution to optimization PWLS in (7) is obtained by using RKHS method. We can read details of RKHS in Aronszajn (1950), Eubank (1988), Wahba (1990), Gu (2002), and Wang (2011). Firstly, we express the model (4) into general smoothing spline regression model (Wang, 2011):

where $i=1,2,\cdots ,n_k$ ; $k=1,2,\cdots ,p$ ; $g_k\in {\mathcal{G}}_k$ is a function which unknown and smooth contained in Hilbert space ${\mathcal{G}}_k$ ; and ${\mathcal{L}}_{t_{\mathit{ki}}}\in {\mathcal{G}}_k$ is a linear function and bounded.

Suppose we may decompose the Hilbert space ${\mathcal{G}}_k$ into direct sum of two subspaces ${\mathcal{F}}_k$ and ${\mathcal{H}}_k$ such that we have:

where ${\mathcal{F}}_k$ is orthogonal to ${\mathcal{H}}_k$ . Hence, for every function $g_k\in {\mathcal{G}}_k$ , $k=1,2,\cdots ,p$ can be expressed as follows:

$g_k=f_k+h_k$ ; $f_k\in {\mathcal{F}}_k$ ; $h_k\in {\mathcal{H}}_k$ .

Next, if $\left\{{\delta }_{k1},{\delta }_{k2},\cdots ,{\delta }_{{\mathit{km}}_k}\right\}$ is basis of space ${\mathcal{F}}_k$ and $\left\{{\xi }_{k1},{\xi }_{k2},\cdots ,{\xi }_{{\mathit{kn}}_k}\right\}$ is basis of space ${\mathcal{H}}_k$ , then we can express every function $g_k\in {\mathcal{G}}_k$ , $k=1,2,\cdots ,p$ as follows:

where ${\delta }_k={({\delta }_{k1},{\delta }_{k2},\cdots ,{\delta }_{{\mathit{km}}_k})}^T$ ; ${\mathbf{c}}_k={(c_{k1},c_{k2},\cdots ,c_{{\mathit{km}}_k})}^T$ ; ${\xi }_k={({\xi }_{k1},{\xi }_{k2},\cdots ,{\xi }_{{\mathit{km}}_k})}^T$ ; and $${\mathbf{d}}_k={(d_{k1},d_{k2},\cdots ,d_{{\mathit{km}}_k})}^T$$

Hereinafter, since ${\mathcal{L}}_{t_{\mathit{ki}}}\in {\mathcal{G}}_k$ is bounded linear function and $g_k\in {\mathcal{G}}_k$ , $k=1,2,\cdots ,p$ then we have: $${\mathcal{L}}_{t_{\mathit{ki}}}{g}_k={\mathcal{L}}_{t_{\mathit{ki}}}\left(f_k+h_k\right)={\mathcal{L}}_{t_{\mathit{ki}}}\left(f_k\right)+{\mathcal{L}}_{t_{\mathit{ki}}}\left(h_k\right)$$ $$=f_k\left(t_{\mathit{ki}}\right)+h_k\left(t_{\mathit{ki}}\right)$$

Based on Eq. (12) and Riesz representation theorem (Wang, 2011), there is a representer ${\omega }_{\mathit{ki}}\in {\mathcal{G}}_k$ of ${\mathcal{L}}_{t_{\mathit{ki}}}$ such that:

${\mathcal{L}}_{t_{\mathit{ki}}}{g}_k=\langle {\omega }_{\mathit{ki}},g_k\rangle =g_k\left(t_{\mathit{ki}}\right)$ ;. $g_k\in {\mathcal{G}}_k$

where $\langle \bullet ,\bullet \rangle$ notates a product of inner. By considering Eq.(11) and inner-product properties, the following equation is obtained:

Next, by using Eq. (13) for $k=1$ we get:

Hence, based on Eq. (14) for $i=1,2,\cdots ,n_1$ we have: $${\mathbf{g}}_1\left({\mathbf{t}}_1\right)={\left(g_1\left(t_{11}\right),\cdots ,g_1\left(t_{{1n}_1}\right)\right)}^T={\mathbf{A}}_1{\mathbf{c}}_1+{\mathbf{B}}_1{\mathbf{d}}_1$$

where ${\mathbf{c}}_1={\left(c_{11},c_{12},\cdots ,c_{{1m}_1}\right)}^T$ ; ${\mathbf{d}}_1={\left(d_{11},d_{12},\cdots ,d_{{1n}_1}\right)}^T$ ;

${\mathbf{A}}_1=\left(\begin{array}{ccc}\langle {\omega }_{11},{\delta }_{11}\rangle & \cdots & \langle {\omega }_{11},{\delta }_{{1m}_1}\rangle \\ ⋮& \ddots & ⋮\\ \langle {\omega }_{{1n}_1},{\delta }_{11}\rangle & \cdots & \langle {\omega }_{{1n}_1},{\delta }_{{1m}_1}\rangle \end{array}\right)$ ; and ${\mathbf{B}}_{\left(1\right)}=\left(\begin{array}{ccc}\langle {\omega }_{11},{\xi }_{11}\rangle & \cdots & \langle {\omega }_{11},{\xi }_{{1n}_1}\rangle \\ ⋮& \ddots & ⋮\\ \langle {\omega }_{{1n}_1},{\xi }_{11}\rangle & \cdots & \langle {\omega }_{{1n}_1},{\xi }_{{1n}_1}\rangle \end{array}\right)$ .

Similarly, we get:

${\mathbf{g}}_2\left({\mathbf{t}}_2\right)={\mathbf{A}}_2{\mathbf{c}}_2+{\mathbf{B}}_2{\mathbf{d}}_2$ , ${\mathbf{g}}_3\left({\mathbf{t}}_3\right)={\mathbf{A}}_3{\mathbf{c}}_3+{\mathbf{B}}_3{\mathbf{d}}_3$ … ${\mathbf{g}}_p\left({\mathbf{t}}_p\right)={\mathbf{A}}_p{\mathbf{c}}_p+{\mathbf{B}}_p{\mathbf{d}}_p$ .

Therefore, generally, the following expression of $\mathbf{g}\left(\mathbf{t}\right)$ is obtained: $$\mathbf{g}\left(\mathit{t}\right)={\left({\mathbf{g}}_1\left({\mathbf{t}}_1\right),{\mathbf{g}}_2\left({\mathbf{t}}_2\right),\cdots ,{\mathbf{g}}_p\left({\mathbf{t}}_p\right)\right)}^{\mathit{T}}={\left({\mathbf{A}}_1{\mathbf{c}}_1,{\mathbf{A}}_2{\mathbf{c}}_2,\cdots ,{\mathbf{A}}_p{\mathbf{c}}_p\right)}^T+{\left({\mathbf{B}}_1{\mathbf{d}}_1,{\mathbf{B}}_2{\mathbf{d}}_2,\cdots ,{\mathbf{B}}_p{\mathbf{d}}_p\right)}^T$$ $$=diag\left({\mathbf{A}}_1,{\mathbf{A}}_2,\cdots ,{\mathbf{A}}_p\right){\left({\mathbf{c}}_1^{\mathrm{T}},{\mathbf{c}}_2^{\mathrm{T}},\cdots ,{\mathbf{c}}_p^T\right)}^T+diag\left({\mathbf{B}}_1,{\mathbf{B}}_2,\cdots ,{\mathbf{B}}_p\right){\left({\mathbf{d}}_1^T,{\mathbf{d}}_2^T,\cdots ,{\mathbf{d}}_p^T\right)}^T$$

where $\mathrm{N}={\sum }_{k=1}^p{n}_k$ ; $\mathrm{M}={\sum }_{k=1}^p{m}_k$ ; A is a matrix with dimension $\mathrm{N}\times \mathrm{M}$ ; c is a vector with dimension $\mathrm{M}\times 1$ ; B is a matrix with dimension $\mathrm{N}\times \mathrm{N}$ ; and d is a vector with dimension $\mathrm{N}\times 1$ . Generally, based on Eq.(15), the MSR model (6) can be written as follows:

Hereafter, to obtain regression function estimation of MSR model (16), we determine the solution to PWLS (7) which can be presented as follows: $$\underset{g_{k\in {\mathcal{G}}_k}}{\mathrm{Min}}\left\{{\Vert {\mathrm{W}}^{\frac{1}{2}}\epsilon \Vert }^2\right\}=\underset{g_{k\in {\mathcal{G}}_k}}{\mathrm{Min}}\left\{{\Vert {\mathrm{W}}^{\frac{1}{2}}\left({\mathbf{y}}^{\ast }-\mathbf{g}\right)\Vert }^2\right\}$$

with constraint $\underset{a_k}{\stackrel{b_k}{\int }}{\left(g_k^{\left(m\right)}\left(t_k\right)\right)}^2{\mathit{dt}}_k<{\gamma }_k$ , ${\gamma }_k\ge 0$ . Solution to the PWLS optimization is same as the solution to the following PWLS optimization:

where ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_p$ are smoothing parameters. These smoothing parameters set the balance between ${\mathrm{N}}^{-1}{\left({\mathbf{y}}^{\ast }-\mathbf{g}\right)}^T\mathbf{W}\left({\mathbf{y}}^{\ast }-\mathbf{g}\right)$ , as goodness of fit, and ${\sum }_{k=1}^p{\lambda }_k{\int }_{a_k}^{b_k}{\left(g_k^{\left(m\right)}\left(t_k\right)\right)}^2{\mathit{dt}}_k$ as the smoothness. To solve PWLS optimization (17), we decompose the penalty in (17) such that we get: $$\sum _{k=1}^p{\lambda }_k\underset{a_k}{\stackrel{b_k}{\int }}{\left(g_k^{\left(m\right)}\left(t_k\right)\right)}^{2}d{t}_k=d^T\varphi Bd$$

where $\varphi =diag({\lambda }_1{I}_{n_1},{\lambda }_2{I}_{n_2},\cdots ,{\lambda }_p{I}_{n_p})$ . Also, we get the goodness of fit: $${\mathrm{N}}^{-1}{\left({\mathbf{y}}^{\ast }-\mathbf{g}\right)}^T\mathbf{W}\left({\mathbf{y}}^{\ast }-\mathbf{g}\right)={\mathrm{N}}^{-1}{\left({\mathbf{y}}^{\ast }-\mathbf{A}\mathbf{c}-\mathbf{B}\mathbf{d}\right)}^T\mathbf{W}\left({\mathbf{y}}^{\ast }-\mathbf{A}\mathbf{c}-\mathbf{B}\mathbf{d}\right)$$

Hence, by combining penalty and goodness of fit, we obtain PWLS optimization whose solutions are:

$\widehat{\mathbf{c}}={\left({\mathbf{A}}^T{\mathbf{D}}^{-1}\mathbf{W}\mathbf{A}\right)}^{-1}{\mathbf{A}}^T{\mathbf{D}}^{-1}\mathbf{W}{\mathbf{y}}^{\mathbf{\ast }}$ and. $\widehat{\mathbf{d}}={\mathbf{D}}^{-1}\mathbf{W}\left[\mathbf{I}-\mathbf{A}{\left({\mathbf{A}}^T{\mathbf{D}}^{-1}\mathbf{W}\mathbf{A}\right)}^{-1}{\mathbf{A}}^T{\mathbf{D}}^{-1}\mathbf{W}\right]{\mathbf{y}}^{\mathbf{\ast }}$

where $D=WB+N\varphi I$ . Therefore, the estimated regression function in nonparametric component of MSR model (1) or (6) is:

Based on Eq. (3), we can express Eq. (18) as:

Hence, the sum of squared errors (SSE) is given by:

Next, by minimizing the SSE, we obtain the estimation of parameter $\beta$ namely $\widehat{\beta }$ as follows:

where $\widehat{\beta }$ is a WLS estimator for parameters in parametric component of MSR model (1). Furthermore, by substituting Eq. (22) into Eq. (18), we get estimator of $\mathbf{g}$ as follows:

where $\widehat{\mathbf{g}}$ is smoothing spline estimator for regression function $\mathbf{g}$ in nonparametric component of MSR model (1).

Finally, by considering MSR model (1) and based on estimation results given by equations (22) and (23), we obtain MSR model estimation based on smoothing spline as follows:

where $\mathbf{H}={\mathbf{H}}_{\mathit{par}}+{\mathbf{H}}_{\mathit{nonpar}}$ ; ${\mathbf{H}}_{\mathit{par}}=\mathbf{X}{\left[{\mathbf{X}}^T{\left(\mathbf{I}-{\mathbf{H}}_{\mathit{\lambda }}\right)}^T\left(\mathbf{I}-{\mathbf{H}}_{\lambda }\right)\mathbf{X}\right]}^{-1}{\mathbf{X}}^T{\left(\mathbf{I}-{\mathbf{H}}_{\lambda }\right)}^T\left(\mathbf{I}-{\mathbf{H}}_{\lambda }\right)$ ; and ${\mathbf{H}}_{\mathit{nonpar}}={\mathbf{H}}_{\lambda }[\mathbf{I}-\mathbf{X}{\left({\mathbf{X}}^{\mathit{T}}{\left(\mathbf{I}-{\mathbf{H}}_{\lambda }\right)}^{\mathit{T}}\left(\mathbf{I}-{\mathbf{H}}_{\lambda }\right)\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\mathit{T}}{\left(\mathbf{I}-{\mathbf{H}}_{\lambda }\right)}^{\mathit{T}}\left(\mathbf{I}-{\mathbf{H}}_{\lambda }\right)]$ .

Based smoothing spline in MSR model, estimator of $\mathbf{g}$ given in (23) is called weighted partial smoothing spline estimator of regression function of MSR model (1).

3.2 Determining confidence interval of β

To determine a CI, we use pivotal quantity (Sahoo, 2013). We assume that ${\epsilon }_{\mathit{ki}}$ in (1) follows Normal distribution that independent and identic with mean zero and variance ${\sigma }_{\mathit{ki}}^2$ or we write ${\epsilon }_{\mathit{ki}}{}_{i.i.d}N(0,{\sigma }_{\mathit{ki}}^2)$ where ${\sigma }_{\mathit{ki}}^2$ is unknown. Next, the $100\left(1-\alpha \right)\%$ CI for ${\beta }_{\mathit{ki}}$ , $k=1,2,\cdots ,p$ ; $i=1,2,\cdots ,n_k$ is designed such that we have a pivotal quantity of parameter ${\beta }_{\mathit{ki}}$ :

where $\mathrm{\Omega }={(\mathbf{I}-{\mathbf{H}}_{\mathit{\lambda }})}^{\mathrm{T}}\left(\mathbf{I}-{\mathbf{H}}_{\mathit{\lambda }}\right)\mathbf{X}$ , $\mathrm{M}\mathrm{S}\mathrm{E}\left(\mathrm{\lambda }\right)=\frac{{(\mathbf{y}-\mathrm{\Omega }{\left({\mathrm{\Omega }}^{\mathrm{T}}\mathrm{\Omega }\right)}^{-1}{\mathrm{\Omega }}^{\mathrm{T}}\mathbf{y})}^{\mathrm{T}}(\mathbf{y}-\mathrm{\Omega }{\left({\mathrm{\Omega }}^{\mathrm{T}}\mathrm{\Omega }\right)}^{-1}{\mathrm{\Omega }}^{\mathrm{T}}\mathbf{y})}{\left({\sum }_{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)}$ , ${\beta }_{\mathit{ki}}$ is the $i^{\mathit{th}}$ element for $k^{\mathit{th}}$ response of parameters vector $\beta$ , and ${\left({\mathrm{\Omega }}^{\mathrm{T}}\mathrm{\Omega }\right)}_{\mathit{ii}}^{-1}$ is diagonal element of ${\left({\mathrm{\Omega }}^{\mathrm{T}}\mathrm{\Omega }\right)}^{-1}$ . We can use GCV or CV instead of MSE to overcome over fitting (Amini & Roozbeh, 2015; Roozbeh, 2018; and Roozbeh et al., 2020).

Hereinafter, if $\mathbf{Z}=\mathrm{\Omega }{\left({\mathrm{\Omega }}^{\mathrm{T}}\mathrm{\Omega }\right)}^{-1}{\mathrm{\Omega }}^{\mathrm{T}}$ , then MSE(λ) in (25) is given by:

where $\mathbf{P}=(\mathbf{I}-\mathbf{Z})$ . Hence, the pivotal quantity (25) can be expressed as follows:

The pivotal quantity (27) follows a $t$ -student distribution with $\left({\sum }_{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)$ degree of freedom.

Furthermore, to determine the 100(1 – $\alpha$ )% CI for ${\beta }_{\mathit{ki}}$ , $k=1,2,\cdots ,p;$ $i=1,2,\cdots ,n_k$ , we must take the solution to probability equation:

where $L_{\mathit{ki}}$ is lower limit of CI and $U_{\mathit{ki}}$ is upper limit of CI, and $(1-\alpha )$ is level of confidence. Next, we substitute Eq.(27) into Eq.(28) so that we get:

We can write Eq.(29) as:

where $=L_{\mathit{ki}}\sqrt{\left(\frac{{\mathbf{y}}^{\mathrm{T}}\mathbf{P}\mathbf{y}}{\left({\sum }_{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)}\mathbf{V}\right)}$ ; $U=U_{\mathit{ki}}\sqrt{\left(\frac{{\mathbf{y}}^{\mathrm{T}}\mathbf{P}\mathbf{y}}{\left({\sum }_{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)}\mathbf{V}\right)}$ ;

$\mathbf{V}={\left({\mathrm{\Omega }}^{\mathrm{T}}\mathrm{\Omega }\right)}_{\mathit{ii}}^{-1}$ and $=\mathbf{I}-{\mathrm{\Omega }\left({\mathrm{\Omega }}^{\mathrm{T}}\mathrm{\Omega }\right)}^{-1}{\mathrm{\Omega }}^{\mathrm{T}}=\mathbf{I}-\mathbf{Z}$ .

If interval length of CI is shortest then the CI is good. Therefore, we find values of $L_{\mathit{ki}}\in \mathbb{R}$ and $U_{\mathit{ki}}\in \mathbb{R}$ that results length of CI in (30) is the shortest. If $\mathit{length}(L_{\mathit{ki}},U_{\mathit{ki}})$ is length of CI in (30), then we have: $$length\left(L_{\mathit{ki}},U_{\mathit{ki}}\right)=$$ $$\left({\widehat{\beta }}_{\mathit{ki}}-L_{\mathit{ki}}\sqrt{\left(\frac{{\mathbf{y}}^{\mathrm{T}}\mathbf{P}\mathbf{y}}{\left(\sum _{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)}\mathbf{V}\right)}\right)-\left({\widehat{\beta }}_{\mathit{ki}}-U_{\mathit{ki}}\sqrt{\left(\frac{{\mathbf{y}}^{\mathrm{T}}\mathbf{P}\mathbf{y}}{\left(\sum _{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)}\mathbf{V}\right)}\right)$$ $$=(U_{\mathit{ki}}-L_{\mathit{ki}})\sqrt{\left(\frac{{\mathbf{y}}^{\mathrm{T}}\mathbf{P}\mathbf{y}}{\left(\sum _{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)}\mathbf{V}\right)}$$

Hence, the shortest length of CI for ${\beta }_{\mathit{ki}}$ is determined by taking the solution to optimization:

that meets the condition:

${\int }_{L_{\mathit{ki}}}^{U_{\mathit{ki}}}T\left(s\right)ds=1-\alpha$ or $K\left(U_{\mathit{ki}}\right)-K\left(L_{\mathit{ki}}\right)-\left(1-\alpha \right)=0$ (32).

where $T(\bullet )$ represents distribution of probability of $t_{\left({\sum }_{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)}$ and $K(\bullet )$ represents distribution of cumulative probability of $t_{\left({\sum }_{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)}$ . Next, by applying Lagrange method, it results equation as follows: $$R\left(L_{\mathit{ki}},U_{\mathit{ki}},\gamma \right)=\left(U_{\mathit{ki}}-L_{\mathit{ki}}\right)\sqrt{\left(\frac{{\mathbf{y}}^{\mathrm{T}}\mathbf{P}\mathbf{y}}{\left(\sum _{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)}\mathbf{V}\right)}+$$

where $\gamma$ is constant of Lagrange. Hereafter, the following equations are obtained:

From equations (34) and (35), we obtain the following relationship:

The Eq.(37) implies $L_{\mathit{ki}}=U_{\mathit{ki}}$ or $L_{\mathit{ki}}=-U_{\mathit{ki}}$ . Since, $L_{\mathit{ki}}=U_{\mathit{ki}}$ is not satisfied, then the shortest CI can be determined from the $L_{\mathit{ki}}$ and $U_{\mathit{ki}}$ values which fulfill:

By using $\left(1-\alpha \right)$ level of confidence $,$ the $L_{\mathit{ki}}$ and $U_{\mathit{ki}}$ values which fulfill condition (38) can be obtained from the $t_{\left({\sum }_{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)}$ distribution table.

Consequently, the shortest smoothing spline CI for parameters of MSR model fulfills the following probability: $$P\left[{\widehat{\beta }}_{\mathit{ki}}-U_{\mathit{ki}}\sqrt{\left(\frac{{\mathbf{y}}^{\mathrm{T}}\mathbf{P}\mathbf{y}}{\left(\sum _{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)}\mathbf{V}\right)}\le {\beta }_{\mathit{ki}}\le {\widehat{\beta }}_{\mathit{ki}}+U_{\mathit{ki}}\sqrt{\left(\frac{{\mathbf{y}}^{\mathrm{T}}\mathbf{P}\mathbf{y}}{\left(\sum _{k=1}^p{n}_k\right)-n_k(1+q_k+r_k)}\mathbf{V}\right)}\right]=1-\alpha$$

where value of $U_{\mathit{ki}}$ can be determined from Eq.(38) which is ${\int }_{U_{\mathit{ki}}}^{\infty }T\left(s\right)ds=\frac{\alpha }{2}$ . Hence, we have: $$P\left[{\widehat{\beta }}_{\mathit{ki}}-t_{(\frac{\alpha }{2};N-n_k(1+q_k+r_k))}\sqrt{\left(\frac{{\mathbf{y}}^{\mathrm{T}}\mathbf{P}\mathbf{y}}{N-n_k(1+q_k+r_k)}\mathbf{V}\right)}\le {\beta }_{\mathit{ki}}\le {\widehat{\beta }}_{\mathit{ki}}+t_{(\frac{\alpha }{2};N-n_k(1+q_k+r_k))}\sqrt{\left(\frac{{\mathbf{y}}^{\mathrm{T}}\mathbf{P}\mathbf{y}}{N-n_k(1+q_k+r_k)}\mathbf{V}\right)}\right]=1-\alpha$$ $$N=\sum _{k=1}^p{n}_k$$

Finally, by using distribution of $t$ -student, the $100\left(1-\alpha \right)\%$ CIs parameters ${\beta }_{\mathit{ki}}$ , $k=1,2,\cdots ,p$ ; $i=1,2,\cdots ,n_k$ of MSR model (1) are:

where $N={\sum }_{k=1}^p{n}_k$ ; $\mathbf{V}={\left({\mathrm{\Omega }}^{\mathrm{T}}\mathrm{\Omega }\right)}_{\mathit{ii}}^{-1}$ ; $\mathbf{P}=\mathbf{I}-{\mathrm{\Omega }\left({\mathrm{\Omega }}^{\mathrm{T}}\mathrm{\Omega }\right)}^{-1}{\mathrm{\Omega }}^{\mathrm{T}}$ ; $\mathrm{\Omega }={(\mathbf{I}-{\mathbf{H}}_{\mathit{\lambda }})}^{\mathrm{T}}\left(\mathbf{I}-{\mathbf{H}}_{\mathit{\lambda }}\right)\mathbf{X}$ ; and ${\mathbf{H}}_{\mathit{\lambda }}$ is given in (19). The asymptotic distribution of ${\widehat{\beta }}_{\mathit{ki}}$ is Normal as presented by Theorem 2 in section 3.3.

3.3 Determining asymptotic distribution

For investigating asymptotic distribution of $\widehat{\beta }$ , we consider the following lemmas and theorem.