Estimation of Multiresponse Multipredictor Nonparametric Regression Model Using Mixed Estimator

Abstract

In data analysis using a nonparametric regression approach, we are often faced with the problem of analyzing a set of data that has mixed patterns, namely, some of the data have a certain pattern and the rest of the data have a different pattern. To handle this kind of datum, we propose the use of a mixed estimator. In this study, we theoretically discuss a developed estimation method for a nonparametric regression model with two or more response variables and predictor variables, and there is a correlation between the response variables using a mixed estimator. The model is called the multiresponse multipredictor nonparametric regression (MMNR) model. The mixed estimator used for estimating the MMNR model is a mixed estimator of smoothing spline and Fourier series that is suitable for analyzing data with patterns that partly change at certain subintervals, and some others that follow a recurring pattern in a certain trend. Since in the MMNR model there is a correlation between responses, a symmetric weight matrix is involved in the estimation process of the MMNR model. To estimate the MMNR model, we apply the reproducing kernel Hilbert space (RKHS) method to penalized weighted least square (PWLS) optimization for estimating the regression function of the MMNR model, which consists of a smoothing spline component and a Fourier series component. A simulation study to show the performance of proposed method is also given. The obtained results are estimations of the smoothing spline component, Fourier series component, MMNR model, weight matrix, and consistency of estimated regression function. In conclusion, the estimation of the MMNR model using the mixed estimator is a combination of smoothing spline component and Fourier series component estimators. It depends on smoothing and oscillation parameters, and it has linear in observation and consistent properties.

1. Introduction

Often, we have to carry out statistical analyses that involve relationships between several variables in the sense of functional relationships between response variables and predictor variables. The appropriate statistical analysis to use in cases like this is regression analysis. In this regression analysis, basically, we build a mathematical model that is usually called a regression model, where the functional relationship between variables in the model is represented by a regression function, and for analysis purposes, we must estimate this regression function. In general, there are two basic types of regression model, namely, a parametric regression model and nonparametric regression model. Meanwhile, if we combine these two basic types of regression model, we will obtain a regression model called a semiparametric regression model. The parametric regression model is suitable for use in cases where the pattern of the relationship between the response variables and predictor variables is known in the sense of indicating a certain curve shape, such as linear, quadratic, cubic, etc., either through the results of initial investigations based on scatter diagrams or from past experience regarding the relationship between these variables. Hereinafter, the nonparametric regression model is suitable for use in cases where the pattern of the relationship between the response variables and predictor variables is unknown in the sense that it does not indicate the shape of a particular curve. The shape of the curve is only assumed to be smooth, that is, it is contained in a Sobolev space [1].

To estimate the regression functions of the regression models such as a parametric regression model, nonparametric regression model, and semiparametric regression model, several smoothing techniques are used, which are represented by estimators. Until now, several estimators have been discussed by several previous researchers both theoretically and in application in several cases. Some of these estimators are local linear estimators, which are used to determine the boundary correction of regression function of the nonparametric regression model [2], to determine the bias reduction in the estimating regression function [3], and to design a standard growth chart used for assessing toddlers’ nutritional status [4]; as a local polynomial estimator used to estimate regression functions in cases of errors-in-variable [5], in case of correlated errors [6], and to estimate the regression function of a functional data regression model [7] and finite population regression model [8]; as a kernel estimator used for estimating nonparametric regression function [9], for estimating and investigating the consistency property of a regression function [10], and for estimating regression function in case of correlated errors [11]. However, these estimators (i.e., local linear, local polynomial, and kernel) are less flexible because they are highly dependent on the neighborhood of the target point, called bandwidth, so that we need a small bandwidth to estimate a fluctuating data model, and this will cause a curve of estimation result that is too rough. Also, these estimators do not include penalty function as smoothness factor; these estimators only consider the goodness of fit factor. This means that these estimators are not so good for estimating models of data that fluctuate in the subintervals, because these estimators will return the estimation results with large value of mean squared errors (MSEs). On the other hand, spline estimators, especially the smoothing spline estimator, have the ability to handle these problems because the splines consider not only the goodness of fit factor but also the smoothness factor [1,12]. Further, the spline estimators such as the smoothing spline, M-type spline, truncated spline, penalized spline, least square spline, linear spline, and B-spline estimators are more flexible than other estimators to use for estimating the nonparametric regression functions, especially for prediction and interpretation purposes [12,13]. These splines have been used and developed widely in several cases by many researchers. For example, Liu et al. [14] and Gao and Shi [15] used M-type splines for analyzing the variance in correlated data, and for estimating regression functions of nonparametric and semiparametric regression models, respectively; Chamidah et al. [16] used truncated splines to estimate mean arterial pressure for prediction purposes, Chamidah et al. [17] and Lestari et al. [18] developed truncated spline and smoothing spline estimators, respectively, for estimating semiparametric regression models and determining the asymptotic properties of the estimator; Tirosh et al. [19], Irizarry [20], Adams et al. [21,22], Lee [23], and Maharani and Saputro [24] discussed smoothing spline for problems of analyzing fractal-like signals, minimizing risk estimate, modeling ARMA observations and estimating smoothing parameter, selection smoothing parameter using simulation data, and determining GCV criterion, respectively; Wang [13], Wang and Ke [25], Gu [26], and Sun et al. [27] discussed smoothing splines in ANOVA models; Wang et al. [28] applied a bivariate smoothing spline to data of cortisol and ACTH hormones; Lu et al. [29] used a penalized spline for analyzing current status data; Berry and Helwig [30] compared tuning methods for penalized splines; Islamiyati et al. [31,32] developed a least square spline for estimating two responses of nonparametric regression models and discussed linear spline in the modeling of blood sugar; and Kirkby et al. [33] estimated nonparametric density using B-Spline. Additionally, Osmani et al. [34] estimated the coefficient of a rates model using kernel and spline estimators.

In statistical modeling such as regression modeling, we often have to analyze the functional relationship between response variable and predictor variable where there are two or more response variables and the between responses are correlated with each other. The regression models that draw a functional relationship between response variables and predictor variables with correlated responses are the multiresponse nonparametric regression (MNR) model and multiresponse semiparametric regression (MSR) model. Due to there being correlation between responses, in the estimation process of a regression function, we need to include a matrix called a weight matrix. The inclusion of a weight matrix in the estimation process is what differentiates the MNR model from a uniresponse nonparametric regression (UNR) model, where there is no correlation between responses. So, the process of estimating the regression functions of the MNR and MSR models requires a symmetric weight matrix, namely, a diagonal matrix. We can also use several estimators to estimate the regression functions of these MNR and MSR models. One of these estimators is the smoothing spline estimator. Currently, many researchers are paying attention to developing and applying the smoothing spline estimator in many areas of research. For example, Chamidah et al. [17], Lestari et al. [18] discussed smoothing splines in MSR models; Lestari et al. [35] and Wang et al. [28] discussed RKHS in an MNR model and applied an MNR model to determine associations between hormones, respectively. Because of its powerful and flexible properties, the smoothing spline estimator is one of the most popular estimators used for estimating the regression functions of UNR and MNR models.

The previous description shows that the smoothing spline estimator has been used in statistical analysis using regression model approaches such as UNR and MNR models. The smoothing spline estimator provides good fitting results because it can overcome the curve of data, which has a sharp decrease and increase pattern that results a relatively smooth curve. Also, using the smoothing spline provides several advantages in that it has unique statistical properties, enables visual interpretation, can handle smooth data and functions, and can readily handle data that change at certain subintervals [1,12,13]. Meanwhile, the Fourier series estimator is a popular smoothing technique in nonparametric regression modeling. The Fourier series also provides good statistical and visual interpretation among nonparametric regression models. Using the Fourier series estimator provides advantages such as being able to handle data that have a repeating pattern at a certain trend interval and to provide good statistical interpretation [36]. Suparti et al. [37], Mardianto et al. [38], and Amato et al. [39] proposed Fourier series methods for modeling inflation in Indonesia, for modeling longitudinal data, and for approximating separable models, respectively.

The research discussed in the previous description relating to the estimation of nonparametric regression models is still limited to the use of only one type of estimator, whereas in data analysis using a nonparametric regression approach, we are often faced with the problem of analyzing a set of data that has mixed patterns, namely, some of the data have a certain pattern and the rest of the data have a different pattern. To handle these kinds of data, we should use a mixed estimator, namely, a combination of more than one estimator. Mariati et al. [40] discussed a uniresponse multivariable nonparametric regression model using a mixed estimator comprising a smoothing spline and Fourier series that was applied to data of poor households in Bali Province. However, previous researchers, namely, Mariati et al. [40], discussed the use of a mixed estimator to estimate a nonparametric regression model, but they only discussed estimating a UNR model in the sense that there is no correlation between the response variables. In this study, therefore, we theoretically discuss a new estimation method using a mixed estimator for a nonparametric regression model with two or more response variables and predictor variables, and there is a correlation between the response variables. The model is called the multiresponse multipredictor nonparametric regression (MMNR) model. The mixed estimator used for estimating the MMNR model is a combination of smoothing spline and Fourier series suitable for analyzing data wherein some patterns partly change at certain subintervals, and others follow a recurring pattern in a certain trend. Since there is a correlation between responses in the MMNR model, a weight matrix is involved in its estimation process. Also, we apply the reproducing kernel Hilbert space (RKHS) method to penalized weighted least square (PWLS) optimization for estimating the regression function of our MMNR model. In this study, therefore, we discuss theoretically about determining the smoothing spline component and Fourier series component of the MMNR model, determining the goodness of fit and penalty components of PWLS optimization, determining the MMNR model estimate, determining the weight matrix, selecting optimal smoothing and oscillation parameters in the MMNR model, and investigating the consistency of the regression function estimator of the MMNR model. Additionally, the results of this study contribute to the development of the field of theoretical statistics, namely, statistical inference theory, regarding theories of estimation and hypothesis testing in multiresponse multipredictor nonparametric regression based on a mixed smoothing spline and Fourier series estimator.

2. Materials and Methods

To achieve the objectives of this research, in this section, we briefly present materials and methods such as the MMNR model, mixed smoothing spline and Fourier series estimator, reproducing kernel Hilbert space (RKHS), and penalized weighted least square (PWLS) optimization.

2.1. The MMNR Model

We first consider a paired dataset $\left(x_{r1i},x_{r2i},\dots ,x_{rpi},t_{r1i},t_{r2i},\dots ,t_{rqi},y_{ri}\right)$ for $r=\mathrm{1,2},\dots ,R$ and $i=\mathrm{1,2},\dots ,n$. We say that the paired dataset follows an MMNR model if the relationship between $y_{ri}$ and ($x_{r1i},x_{r2i},\dots ,x_{rpi},t_{r1i},t_{r2i},\dots ,t_{rqi}$) satisfies the following functional equation:

(1)$y_{ri}=m_r\left(x_{r1i},x_{r2i},\dots ,x_{rpi},t_{r1i},t_{r2i},\dots ,t_{rqi}\right)+{\epsilon }_{ri}$

where $r=\mathrm{1,2},\dots ,R$ represents responses; $i=\mathrm{1,2},\dots ,n$ represents observations; $y_{ri}$ is the value of the response variable at the $i\mathrm{th}$ observation and $r\mathrm{th}$ response; $m_r$ is the unknown regression function of the $r\mathrm{th}$ response; $x_{r1i},x_{r2i},\dots ,x_{rpi}$ are values of $p$ smoothing spline predictor variables at the $i\mathrm{th}$ observation and $r\mathrm{th}$ response; $t_{r1i},t_{r2i},\dots ,t_{rqi}$ are values of $q$ Fourier series predictor variables at the $i\mathrm{th}$ observation and $r\mathrm{th}$ response; and ${\epsilon }_{ri}$ is the random error at the $i\mathrm{th}$ observation and $r\mathrm{th}$ response, which has mean of zero and variance of ${\sigma }_{ri}^2$. Here, we assume that the correlation between responses is ${\rho }_{ls}=\left\{\begin{array}{c}\begin{array}{ccc}{\rho }_l& for& l=s\end{array}\\ \begin{array}{ccc}0& for& l\ne s\end{array}\end{array}\right.$.

Since the shape of the regression function $m_r$ in the MMNR model is assumed to be unknown and additive in nature, the regression function $m_r$ in Equation (1) can be presented as follows:

(2)$m_r\left(x_{r1i},x_{r2i},\dots ,x_{rpi},t_{r1i},t_{r2i},\dots ,t_{rqi}\right)={\sum }_{j=1}^p{f}_{rj}\left(x_{rji}\right)+{\sum }_{k=1}^q{g}_{rk}\left(t_{rki}\right)$

where $r=\mathrm{1,2},\dots ,R$; $i=\mathrm{1,2},\dots ,n$; ${\sum }_{j=1}^p{f}_{rj}\left(x_{rji}\right)$ is a smoothing spline component of regression function $m_r$ in which $f_{rj}\left(x_{rji}\right)$, $j=\mathrm{1,2},\dots ,p$ are contained in a Sobolev space $W_2^m\left[a_j,b_j\right]$, namely, $f_{rj}\in W_2^m\left[a_j,b_j\right]$; and ${\sum }_{k=1}^q{g}_{rk}\left(t_{rki}\right)$ is a Fourier series component of regression function $m_r$ that is approximated by a Fourier series function. Hereinafter, based on Equations (1) and (2), we have the following MMNR model:

(3)$y_{ri}={\sum }_{j=1}^p{f}_{rj}\left(x_{rji}\right)+{\sum }_{k=1}^q{g}_{rk}\left(t_{rki}\right)+{\epsilon }_{ri}; r=\mathrm{1,2},\dots ,R; i=\mathrm{1,2},\dots ,n.$

Furthermore, the regression function of the MMNR model presented in (3) is estimated using a mixed smoothing spline and Fourier series estimator by applying the reproducing kernel Hilbert space (RKHS) method to the developed penalized weighted least square.

2.2. Mixed Smoothing Spline and Fourier Series Estimator

To estimate the regression function of the MMNR model presented in (3), we use a mixed smoothing spline and Fourier series estimator, which is suitable for analyzing data with some patterns that partly change at certain subintervals, and others that follow a recurring pattern in a certain trend. The regression function $m_r$ in the MMNR model in (3) consists of the component ${\sum }_{j=1}^p{f}_{rj}\left(x_{rji}\right)$, which is approached by smoothing spline, and the component ${\sum }_{k=1}^q{g}_{rk}\left(t_{rki}\right)$, which is approximated by the following Fourier series function:

(4)$g_{rk}\left(t_{rki}\right)=b_{rk}{t}_{rki}+\left(\frac{1}{2}\right){\alpha }_{0rk}+{\sum }_{v=1}^K{\alpha }_{vrk}\cos v{t}_{rki}$

where $r=\mathrm{1,2},\dots ,R$; $k=\mathrm{1,2},\dots ,p$; and $i=\mathrm{1,2},\dots ,n$.

Since we involve a smoothing spline estimator in the process of estimating the regression function, according to Wahba [12] and Wang [13], we should use the reproducing kernel Hilbert space (RKHS) method approach, which is applied to the developed penalized weighted least square (PWLS) optimization.

2.3. Reproducing Kernel Hilbert Space (RKHS)

The reproducing kernel Hilbert space (RKHS) method was first introduced by Aronszajn [41]. The RKHS method was developed for estimating nonparametric regression models and semiparametric regression models by several researchers such as Wang [13], Chamidah et al. [17], Lestari et al. [18], and Kimeldorf and Wahba [42]. We call a Hilbert space $\mathcal{H}$ an RKHS on a set X over field $\mathbb{F}$ if it meets certain conditions, as shown by Aronszajn [41], Berlinet and Thomas-Agnan [43], Paulsen [44], and Yuan and Cai [45]:

• (i). If $\mathcal{H}$ is a vector, then $\mathcal{H}\subset (\mathrm{X},\mathbb{F})$, namely, $\mathcal{H}$ is the subspace of a vector space over $\mathbb{F}$, which is notated by $\mathcal{F}$(X,$\mathbb{F}$);

• (ii). If $\mathcal{H}$ is equipped with an inner product, $⟨ ,⟩$, then it will be a Hilbert space;

• (iii). If $E_y : \mathcal{H}\to \mathbb{F}$ is a linear evaluation functional that is defined by $E_y\left(f\right)=f\left(y\right)$ for every $y\in$ X, then the linear evaluation functional is bounded.

Furthermore, if $\mathcal{H}\subset \mathrm{X}$, namely, $\mathcal{H}$ is an RKHS on $\mathrm{X}$, then there exists a unique vector, $k_y\in \mathcal{H}$, for every $y\in \mathrm{X}$, such that $f\left(y\right)=⟨f , k_y⟩$ for every $f\in \mathcal{H}$. Note that every linear functional that is bounded will be given by the inner product with a unique vector in $\mathcal{H}$, and under this condition, we call the function $k_y$ as a reproducing kernel (RK) for the point $y$. This means that the reproducing kernel (RK) for $\mathcal{H}$ is a two-variable function, which is defined by $K\left(x,y\right)=k_y\left(x\right)$. It implies that $K\left(x,y\right)=k_y\left(x\right)=⟨k_y,k_x⟩$ and ${‖E_y‖}^2={‖k_y‖}^2=⟨k_y,k_y⟩=K(y,y)$ [41,42,43,44,45].

In this study, the RKHS method is employed to the developed penalized weighted least square (PWLS) optimization for obtaining the estimation of the regression function of the MMNR model presented in (3). In the following section, we provide the penalized weighted least square (PWLS) optimization that we developed according to the case of our interest.

2.4. Penalized Weighted Least Square (PWLS) Optimization

Since the regression function of the multiresponse multipredictor nonparametric regression (MMNR) model in (3) consists of two components, namely, the smoothing spline component ${\sum }_{j=1}^p{f}_{rj}\left(x_{rji}\right)$ and the Fourier series component ${\sum }_{k=1}^q{g}_{rk}\left(t_{rki}\right)$, for estimating the regression function, we developed the PWLS optimization for a single-estimator case, namely, the smoothing spline estimator, only as proposed by Wang et al. [28], Lestari et al. [18], and Islamiyati et al. [31], to the two mixed estimators, namely, the smoothing spline and Fourier series estimators. In the process of estimating the regression function of the multiresponse multipredictor nonparametric regression (MMNR) model by using the developed penalized weighted least square (PWLS) optimization, we assume that $g_{rk}\left(t_{rki}\right)$ is a fixed function. Hence, we can estimate the smoothing spline component in the multiresponse multipredictor nonparametric regression (MMNR) model presented in (3) for every response $r=\mathrm{1,2},\dots ,R$ by using the following developed penalized weighted least square (PWLS) optimization:

(5)$\begin{array}{c}\underset{f_{\mathit{rj}}\in W_2^m[a_j,b_j]}{\mathrm{Min}}\left\{n^{-1}{\sum }_{i=1}^n{w}_{ri}\left(y_{ri}-\left[{\sum }_{j=1}^p{f}_{rj}\left(x_{rji}\right)\right.\right.\right.+{\sum }_{k=1}^q\left(b_{rk}{t}_{rki}+\frac{1}{2}{\alpha }_{0rk}+\right.\\ {\left.\left.\left.{\sum }_{v=1}^K{\alpha }_{vrk}\cos v{t}_{rki}\right)\right]\right)}^2+\left.{\sum }_{j=1}^p{\lambda }_{rj}{\int }_{a_j}^{b_j}{\left[f_{rj}\left(x_{rj}\right)\right]}^{2}d{x}_{rj}\right\} \end{array}$

where $0<{\lambda }_{rj}<\infty$; $r=\mathrm{1,2},\dots ,R$; $j=\mathrm{1,2},\dots ,p$; and $w_{ri}$ represents a weight that is an inverse of variances. Note that we include a weight $w_{ri}$ in the developed penalized weighted least square (PWLS) presented in (5) because in this case there is a correlation between responses.

3. Results and Discussions

In this section, we theoretically discuss the estimating process of the multiresponse multipredictor nonparametric regression (MMNR) model. Firstly, suppose that we have a paired dataset such as $\left(x_{r1i},x_{r2i},\dots ,x_{rpi},t_{r1i},t_{r2i},\dots ,t_{rqi},y_{ri}\right)$ for $r=\mathrm{1,2},\dots ,R$ and $i=\mathrm{1,2},\dots ,n$, where the functional relationship between predictor variables $\left(x_{r1i},x_{r2i},\dots ,x_{rpi},t_{r1i},t_{r2i},\dots ,t_{rqi}\right)$ and response variables $y_{ri}$ meets the MMNR model presented in Equation (3), namely:

$\begin{array}{ll}{y}_{ri}& =m_r\left(x_{r1i},x_{r2i},\dots ,x_{rpi},t_{r1i},t_{r2i},\dots ,t_{rqi}\right)+{\epsilon }_{ri}\\ & ={\sum }_{j=1}^p{f}_{rj}\left(x_{rji}\right)+{\sum }_{k=1}^q{g}_{rk}\left(t_{rki}\right)+{\epsilon }_{ri}; r=\mathrm{1,2},\dots ,R ; i=\mathrm{1,2},\dots ,n\end{array}$

where information and assumptions about the multiresponse multipredictor nonparametric regression (MMNR) model are presented in Section 2.1.

Since this multiresponse multipredictor nonparametric regression (MMNR) model has more than one response variable and predictor variable, and has a regression function that is composed of two components, namely, a smoothing spline component and a Fourier series component, to simplify the process of estimating the multiresponse multipredictor nonparametric regression (MMNR) model by using the reproducing kernel Hilbert space (RKHS) method approach that is employed to the developed penalized weighted least square (PWLS) optimization, we should express the multiresponse multipredictor nonparametric regression (MMNR) model presented in Equation (3) in the form of matrix notation, as follows:

(6)$\underset{\left(Rn\times 1\right)}{\underbrace{\left[\begin{array}{c}{y}_{11}\\ y_{12}\\ \begin{array}{c}⋮\\ y_{1n}\\ \begin{array}{c}{y}_{21}\\ y_{22}\\ \begin{array}{c}⋮\\ y_{2n}\\ \begin{array}{c}⋮\\ y_{R1}\\ \begin{array}{c}{y}_{R2}\\ ⋮\\ y_{Rn}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]}}=\underset{\left(Rn\times 1\right)}{\underbrace{\left[\begin{array}{c}{\sum }_{j=1}^p{f}_{1j}\left(x_{1j1}\right)\\ {\sum }_{j=1}^p{f}_{1j}\left(x_{1j2}\right)\\ \begin{array}{c}⋮\\ {\sum }_{j=1}^p{f}_{1j}\left(x_{1jn}\right)\\ \begin{array}{c}{\sum }_{j=1}^p{f}_{2j}\left(x_{2j1}\right)\\ {\sum }_{j=1}^p{f}_{2j}\left(x_{2j2}\right)\\ \begin{array}{c}⋮\\ {\sum }_{j=1}^p{f}_{2j}\left(x_{2jn}\right)\\ \begin{array}{c}⋮\\ {\sum }_{j=1}^p{f}_{Rj}\left(x_{Rj1}\right)\\ \begin{array}{c}{\sum }_{j=1}^p{f}_{Rj}\left(x_{Rj2}\right)\\ ⋮\\ {\sum }_{j=1}^p{f}_{Rj}\left(x_{Rjn}\right)\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]}}+\underset{\left(Rn\times 1\right)}{\underbrace{\left[\begin{array}{c}{\sum }_{k=1}^q{g}_{1k}\left(t_{1k1}\right)\\ {\sum }_{k=1}^q{g}_{1k}\left(t_{1k2}\right)\\ \begin{array}{c}⋮\\ {\sum }_{k=1}^q{g}_{1k}\left(t_{1kn}\right)\\ \begin{array}{c}{\sum }_{k=1}^q{g}_{1k}\left(t_{1kn}\right)\\ {\sum }_{k=1}^q{g}_{2k}\left(t_{2k2}\right)\\ \begin{array}{c}⋮\\ {\sum }_{k=1}^q{g}_{2k}\left(t_{2kn}\right)\\ \begin{array}{c}⋮\\ {\sum }_{k=1}^q{g}_{Rk}\left(t_{Rk1}\right)\\ \begin{array}{c}{\sum }_{k=1}^q{g}_{Rk}\left(t_{Rk2}\right)\\ ⋮\\ {\sum }_{k=1}^q{g}_{Rk}\left(t_{Rkn}\right)\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]}}+\underset{\left(Rn\times 1\right)}{\underbrace{\left[\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{12}\\ \begin{array}{c}⋮\\ {\epsilon }_{1n}\\ \begin{array}{c}{\epsilon }_{21}\\ {\epsilon }_{22}\\ \begin{array}{c}⋮\\ {\epsilon }_{2n}\\ \begin{array}{c}⋮\\ {\epsilon }_{R1}\\ \begin{array}{c}{\epsilon }_{R2}\\ ⋮\\ {\epsilon }_{Rn}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]}}.$

Furthermore, we may write the MMNR model presented in (6) as follows:

(7)$\underset{\left(Rn\times 1\right)}{\underbrace{\left[\begin{array}{c}{y}_{11}\\ y_{12}\\ \begin{array}{c}⋮\\ y_{1n}\\ \begin{array}{c}{y}_{21}\\ y_{22}\\ \begin{array}{c}⋮\\ y_{2n}\\ \begin{array}{c}⋮\\ y_{R1}\\ \begin{array}{c}{y}_{R2}\\ ⋮\\ y_{Rn}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]}}=\underset{\left(Rn\times 1\right)}{\underbrace{\left[\begin{array}{c}{f}_1\left(x_{11}\right)\\ f_1\left(x_{12}\right)\\ \begin{array}{c}⋮\\ f_1\left(x_{1n}\right)\\ \begin{array}{c}{f}_2\left(x_{21}\right)\\ f_2\left(x_{22}\right)\\ \begin{array}{c}⋮\\ f_2\left(x_{2n}\right)\\ \begin{array}{c}⋮\\ f_R\left(x_{R1}\right)\\ \begin{array}{c}{f}_R\left(x_{R2}\right)\\ ⋮\\ f_R\left(x_{Rn}\right)\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]}}+\underset{\left(Rn\times 1\right)}{\underbrace{\left[\begin{array}{c}g\left(t_{11}\right)\\ g_1\left(t_{12}\right)\\ \begin{array}{c}⋮\\ g_1\left(t_{1n}\right)\\ \begin{array}{c}{g}_2\left(t_{21}\right)\\ g_2\left(t_{22}\right)\\ \begin{array}{c}⋮\\ g_2\left(t_{2n}\right)\\ \begin{array}{c}⋮\\ g_R\left(t_{R1}\right)\\ \begin{array}{c}{g}_R\left(t_{R2}\right)\\ ⋮\\ g_R\left(t_{Rn}\right)\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]}}+\underset{\left(Rn\times 1\right)}{\underbrace{\left[\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{12}\\ \begin{array}{c}⋮\\ {\epsilon }_{1n}\\ \begin{array}{c}{\epsilon }_{21}\\ {\epsilon }_{22}\\ \begin{array}{c}⋮\\ {\epsilon }_{2n}\\ \begin{array}{c}⋮\\ {\epsilon }_{R1}\\ \begin{array}{c}{\epsilon }_{R2}\\ ⋮\\ {\epsilon }_{Rn}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]}}.$

Hence, the MMNR model given in (7) can be written as follows:

(8)$\mathbf{y}=\mathbf{f}+\mathbf{g}+\mathsf{\epsilon }$

where $\mathbf{y}={\left[{\mathbf{y}}_{\mathbf{1}},{\mathbf{y}}_{\mathbf{2}},\dots ,{\mathbf{y}}_{\mathbf{R}}\right]}^{\mathrm{T}}$; $\mathbf{f}={\left[{\mathbf{f}}_{\mathbf{1}}\left({\mathbf{x}}_{\mathbf{1}}\right),{\mathbf{f}}_{\mathbf{2}}\left({\mathbf{x}}_{\mathbf{2}}\right),\dots ,{\mathbf{f}}_{\mathbf{R}}\left({\mathbf{x}}_{\mathbf{R}}\right)\right]}^{\mathbf{T}}$; $\mathsf{\epsilon }={\left[{\mathsf{\epsilon }}_{\mathbf{1}},{\mathsf{\epsilon }}_{\mathbf{2}},\dots ,{\mathsf{\epsilon }}_{\mathbf{R}}\right]}^{\mathrm{T}}$;

• $\mathbf{g}={\left[{\mathbf{g}}_{\mathbf{1}}\left({\mathbf{t}}_{\mathbf{1}}\right),{\mathbf{g}}_{\mathbf{2}}\left({\mathbf{t}}_{\mathbf{2}}\right),\dots ,{\mathbf{g}}_{\mathbf{R}}\left({\mathbf{t}}_{\mathbf{R}}\right)\right]}^{\mathbf{T}}$; ${\mathbf{y}}_{\mathbf{1}}={\left[y_{11},y_{12},\dots ,y_{1n}\right]}^{\mathrm{T}}$; ${\mathbf{y}}_{\mathbf{2}}={\left[y_{21},y_{22},\dots ,y_{2n}\right]}^{\mathrm{T}}$;…; ${\mathbf{y}}_{\mathit{R}}={\left[y_{R1},y_{R2},\dots ,y_{Rn}\right]}^{\mathrm{T}}$; ${\mathbf{f}}_{\mathbf{1}}\left({\mathbf{x}}_{\mathbf{1}}\right)={\left[f_1\left(x_{11}\right),f_1\left(x_{12}\right),\dots ,f_1\left(x_{1n}\right)\right]}^{\mathbf{T}}$;

• ${\mathbf{f}}_{\mathbf{2}}\left({\mathbf{x}}_{\mathbf{2}}\right)={\left[f_2\left(x_{21}\right),f_2\left(x_{22}\right),\dots ,f_2\left(x_{2n}\right)\right]}^{\mathbf{T}}$;…; ${\mathbf{f}}_{\mathit{R}}\left({\mathbf{x}}_{\mathbf{R}}\right)={\left[f_R\left(x_{R1}\right),f_R\left(x_{R2}\right),\dots ,f_R\left(x_{Rn}\right)\right]}^{\mathbf{T}}$;

• ${\mathbf{g}}_{\mathbf{1}}\left({\mathbf{t}}_{\mathbf{1}}\right)={\left[g_1\left(t_{11}\right),g_1\left(t_{12}\right),\dots ,g_1\left(t_{1n}\right)\right]}^{\mathbf{T}}$; ${\mathbf{g}}_{\mathbf{2}}\left({\mathbf{t}}_{\mathbf{2}}\right)={\left[g_2\left(t_{21}\right),g_2\left(t_{22}\right),\dots ,g_2\left(t_{2n}\right)\right]}^{\mathbf{T}}$;…;

• ${\mathbf{g}}_{\mathit{R}}\left({\mathbf{t}}_{\mathbf{R}}\right)={\left[g_R\left(t_{R1}\right),g_R\left(t_{R2}\right),\dots ,g_R\left(t_{Rn}\right)\right]}^{\mathbf{T}}$; ${\mathsf{\epsilon }}_{\mathbf{1}}={\left[{\epsilon }_{11},{\epsilon }_{12},\dots ,{\epsilon }_{1n}\right]}^{\mathrm{T}}$;

• ${\mathsf{\epsilon }}_{\mathbf{2}}={\left[{\epsilon }_{21},{\epsilon }_{22},\dots ,{\epsilon }_{2n}\right]}^{\mathrm{T}}$; …; ${\mathsf{\epsilon }}_{\mathbf{R}}={\left[{\epsilon }_{R1},{\epsilon }_{R2},\dots ,{\epsilon }_{Rn}\right]}^{\mathrm{T}}$.

Next, by letting $\mathbf{z}=\mathbf{y}-\mathbf{g}$, we may write the MMNR model given in (8) as follows:

(9)$\mathbf{z}=\mathbf{f}+\mathsf{\epsilon }$

In the following sections, we discuss about determining smoothing spline component and Fourier series component of MMNR model, determining goodness of fit and penalty components of PWLS optimization, estimating the MMNR model, estimating the weight matrix, selecting optimal smoothing and oscillation parameters in the MMNR model, and investigating the consistency of regression function estimator of MMNR model.

3.1. Determining Smoothing Spline Component of MMNR Model

To determine the function of the smoothing spline component of the regression function of multiresponse multipredictor nonparametric regression (MMNR) model presented in Equation (1), we use the RKHS method, which is applied to PWLS optimization. Due to ${\sum }_{j=1}^p{f}_{rj}\left(x_{rji}\right)$ being a smoothing spline component in the MMNR model presented in (3), we use the RKHS method to estimate it with the following steps: Firstly, suppose that we decompose a Hilbert space $\mathcal{H}$ into a direct sum of Hilbert subspaces ${\mathcal{H}}_0$ with basis $\left\{{\gamma }_{rj1},{\gamma }_{rj2},\dots ,{\gamma }_{rjm}\right\}$ where $m$ represents the order of polynomial spline, ${\mathcal{H}}_1$ with basis $\left\{{\beta }_{rj1},{\beta }_{rj2},\dots ,{\beta }_{rjn}\right\}$ where $n$ represents the number of observations, and ${\mathcal{H}}_0\perp {\mathcal{H}}_1$ (i.e., ${\mathcal{H}}_0$ is perpendicular to ${\mathcal{H}}_1$), as follows:

$\mathcal{H}={\mathcal{H}}_0\oplus {\mathcal{H}}_1.$

Hence, we express every function where $f_{rj}\in \mathcal{H}$ as follows:

(10)$\begin{array}{ll}{f}_{rj}& =u_{rj}+v_{rj}={\sum }_{k=1}^m{c}_{rjk}{\gamma }_{rjk}+{\sum }_{s=1}^n{d}_{rjs}{\beta }_{rjs}\\ & ={\mathsf{\gamma }}_{\mathit{r}\mathit{j}}^{\mathit{T}}{\mathbf{c}}_{\mathit{r}\mathit{j}}+{\mathsf{\beta }}_{\mathit{r}\mathit{j}}^{\mathit{T}}{\mathbf{d}}_{\mathit{r}\mathit{j}}\end{array}$

where $u_{rj}\in {\mathcal{H}}_0$; $v_{rj}\in {\mathcal{H}}_1$; $c_{rj}$ and $d_{rj}$ are constans; ${\mathsf{\gamma }}_{\mathit{r}\mathit{j}}={\left({\gamma }_{rj1},{\gamma }_{rj2},\dots ,{\gamma }_{rjm}\right)}^T$; ${\mathsf{\beta }}_{\mathit{r}\mathit{j}}={\left({\beta }_{rj1},{\beta }_{rj2},\dots ,{\beta }_{rjn}\right)}^T$; ${\mathbf{c}}_{\mathit{r}\mathit{j}}={\left(c_{rj1},c_{rj2},\dots ,c_{rjm}\right)}^T$; ${\mathbf{d}}_{\mathit{r}\mathit{j}}={\left(d_{r1j},d_{r2j},\dots ,d_{rmj}\right)}^T$.

Secondly, suppose ${\mathcal{L}}_x\in \mathcal{H}$ is a bounded linear functional, and $f_{rj}\in \mathcal{H}$, then we have the following relationship:

(11)$\begin{array}{ll}{\mathcal{L}}_x{f}_{rj}& ={\mathcal{L}}_x(u_{rj}+v_{rj})\\ & =u_{rj}\left(x_{ri}\right)+v_{rj}\left(x_{ri}\right)\\ & =f_{rj}\left(x_{ri}\right).\end{array}$

Since ${\mathcal{L}}_x\in \mathcal{H}$ is a bounded linear functional, then according to Yuan and Cai [45] and Lestari et al. [35], there exists a representer ${\xi }_i\in \mathcal{H}$ that meets the following equation:

(12)${\mathcal{L}}_x{f}_{rj}=⟨{\xi }_{rji},f_{rj}⟩=f_{rj}\left(x_{ri}\right).$

Based on Equations (10) and (11), we may write Equation (12) as follows:

(13)$f_{rj}\left(x_{ri}\right)=⟨{\xi }_{rji},f_{rj}⟩=⟨{\xi }_{rji},{\mathsf{\gamma }}_{\mathit{r}\mathit{j}}^{\mathit{T}}{\mathbf{c}}_{\mathit{r}\mathit{j}}⟩+⟨{\xi }_{rji},{\mathsf{\beta }}_{\mathit{r}\mathit{j}}^{\mathit{T}}{\mathbf{d}}_{\mathit{r}\mathit{j}}⟩.$

Hence, for $j=1$, we may write Equation (13) as follows:

(14)$f_{r1}\left(x_{ri}\right)=⟨{\xi }_{r1i},{\mathsf{\gamma }}_{\mathit{r}\mathbf{1}}^{\mathit{T}}{\mathbf{c}}_{\mathit{r}\mathbf{1}}⟩+⟨{\xi }_{r1i},{\mathsf{\beta }}_{\mathit{r}\mathbf{1}}^{\mathit{T}}{\mathbf{d}}_{\mathit{r}\mathbf{1}}⟩, i=\mathrm{1,2},\dots ,n.$

Next, based on Equation (14) for $j=1$ and $i=1$, we have $f_{r1}\left(x_{r1}\right)$ as follows:

(15)$$\begin{gathered} \begin{array}{ll}{f}_{r1}\left(x_{r1}\right)& =⟨{\xi }_{r11},{\mathsf{\gamma }}_{\mathit{r}\mathbf{1}}^{\mathit{T}}{\mathbf{c}}_{\mathit{r}\mathbf{1}}⟩+⟨{\xi }_{r11},{\mathsf{\beta }}_{\mathit{r}\mathbf{1}}^{\mathit{T}}{\mathbf{d}}_{\mathit{r}\mathbf{1}}⟩\\ & =c_{r11}⟨{\xi }_{r11},{\gamma }_{r11}⟩+c_{r12}⟨{\xi }_{r11},{\gamma }_{r12}⟩+\dots +c_{r1n}⟨{\xi }_{r11},{\gamma }_{r1m}⟩+ \\ {d}_{r11}⟨{\xi }_{r11},{\beta }_{r11}⟩+d_{r12}⟨{\xi }_{r11},{\beta }_{r12}⟩+\dots +d_{r1n}⟨{\xi }_{r11},{\beta }_{r1n}⟩.\end{array} \end{gathered}$$

Similarly, based on Equation (14) for $j=1$ and $i=2$, we have $f_{r1}\left(x_{r2}\right)$ as follows:

(16)$$\begin{gathered} \begin{array}{ll}{f}_{r1}\left(x_{r2}\right)& =⟨{\xi }_{r12},{\mathsf{\gamma }}_{\mathit{r}\mathbf{1}}^{\mathit{T}}{\mathbf{c}}_{\mathit{r}\mathbf{1}}⟩+⟨{\xi }_{r12},{\mathsf{\beta }}_{\mathit{r}\mathbf{1}}^{\mathit{T}}{\mathbf{d}}_{\mathit{r}\mathbf{1}}⟩\\ & =c_{r11}⟨{\xi }_{r12},{\gamma }_{r11}⟩+c_{r12}⟨{\xi }_{r12},{\gamma }_{r12}⟩+\dots +c_{r1n}⟨{\xi }_{r12},{\gamma }_{r1m}⟩+ \\ {d}_{r11}⟨{\xi }_{r12},{\beta }_{r11}⟩+d_{r12}⟨{\xi }_{r12},{\beta }_{r12}⟩+\dots +d_{r1n}⟨{\xi }_{r12},{\beta }_{r1n}⟩.\end{array} \end{gathered}$$

If we continue this process for $j=1$ and $i=\mathrm{3,4},\dots ,n$ in a similar way, then for $i=n$, we have $f_{r1}\left(x_{rn}\right)$ as follows:

(17)$$\begin{gathered} \begin{array}{ll}{f}_{r1}\left(x_{rn}\right)& =⟨{\xi }_{r1n},{\mathsf{\gamma }}_{\mathit{r}\mathbf{1}}^{\mathit{T}}{\mathbf{c}}_{\mathit{r}\mathbf{1}}⟩+⟨{\xi }_{r1n},{\mathsf{\beta }}_{\mathit{r}\mathbf{1}}^{\mathit{T}}{\mathbf{d}}_{\mathit{r}\mathbf{1}}⟩\\ & =c_{r11}⟨{\xi }_{r1n},{\gamma }_{r11}⟩+c_{r12}⟨{\xi }_{r1n},{\gamma }_{r12}⟩+\dots +c_{r1n}⟨{\xi }_{r1n},{\gamma }_{r1m}⟩+ \\ {d}_{r11}⟨{\xi }_{r1n},{\beta }_{r11}⟩+d_{r12}⟨{\xi }_{r1n},{\beta }_{r12}⟩+\dots +d_{r1n}⟨{\xi }_{r1n},{\beta }_{r1n}⟩.\end{array} \end{gathered}$$

Furthermore, based on Equation (13) for $j=2$, we have $f_{r2}\left(x_{ri}\right)$ as follows:

(18)$f_{r2}\left(x_{ri}\right)=⟨{\xi }_{r2i},{\mathsf{\gamma }}_{\mathit{r}\mathbf{2}}^{\mathit{T}}{\mathbf{c}}_{\mathit{r}\mathbf{2}}⟩+⟨{\xi }_{r2i},{\mathsf{\beta }}_{\mathit{r}\mathbf{2}}^{\mathit{T}}{\mathbf{d}}_{\mathit{r}\mathbf{2}}⟩, i=\mathrm{1,2},\dots ,n.$

Hence, based on Equation (17) for $j=2$ and $i=1$, we have $f_{r2}\left(x_{r1}\right)$ as follows:

(19)$$\begin{gathered} \begin{array}{ll}{f}_{r2}\left(x_{r1}\right)& =⟨{\xi }_{r21},{\mathsf{\gamma }}_{\mathit{r}\mathbf{2}}^{\mathit{T}}{\mathbf{c}}_{\mathit{r}\mathbf{2}}⟩+⟨{\xi }_{r21},{\mathsf{\beta }}_{\mathit{r}\mathbf{2}}^{\mathit{T}}{\mathbf{d}}_{\mathit{r}\mathbf{2}}⟩\\ & =c_{r21}⟨{\xi }_{r21},{\gamma }_{r21}⟩+c_{r22}⟨{\xi }_{r21},{\gamma }_{r22}⟩+\dots +c_{r2n}⟨{\xi }_{r21},{\gamma }_{r2m}⟩+ \\ {d}_{r21}⟨{\xi }_{r21},{\beta }_{r21}⟩+d_{r22}⟨{\xi }_{r21},{\beta }_{r22}⟩+\dots +d_{r2n}⟨{\xi }_{r21},{\beta }_{r2n}⟩.\end{array} \end{gathered}$$

In similar way, based on Equation (17) for $j=2$ and $i=2$, we have $f_{r2}\left(x_{r2}\right)$ as follows:

(20)$$\begin{gathered} \begin{array}{ll}{f}_{r2}\left(x_{r2}\right)& =⟨{\xi }_{r22},{\mathsf{\gamma }}_{\mathit{r}\mathbf{2}}^{\mathit{T}}{\mathbf{c}}_{\mathit{r}\mathbf{2}}⟩+⟨{\xi }_{r22},{\mathsf{\beta }}_{\mathit{r}\mathbf{2}}^{\mathit{T}}{\mathbf{d}}_{\mathit{r}\mathbf{2}}⟩\\ & =c_{r21}⟨{\xi }_{r22},{\gamma }_{r21}⟩+c_{r22}⟨{\xi }_{r22},{\gamma }_{r22}⟩+\dots +c_{r2n}⟨{\xi }_{r22},{\gamma }_{r2m}⟩+ \\ {d}_{r21}⟨{\xi }_{r22},{\beta }_{r21}⟩+d_{r22}⟨{\xi }_{r22},{\beta }_{r22}⟩+\dots +d_{r2n}⟨{\xi }_{r22},{\beta }_{r2n}⟩.\end{array} \end{gathered}$$

We continue this process for $j=2$ and $i=\mathrm{3,4},\dots ,n$ in a similar way such that for $i=n$ we have $f_{r2}\left(x_{rn}\right)$ as follows:

(21)$$\begin{gathered} \begin{array}{ll}{f}_{r2}\left(x_{rn}\right)& =⟨{\xi }_{r2n},{\mathsf{\gamma }}_{\mathit{r}\mathbf{2}}^{\mathit{T}}{\mathbf{c}}_{\mathit{r}\mathbf{2}}⟩+⟨{\xi }_{r2n},{\mathsf{\beta }}_{\mathit{r}\mathbf{2}}^{\mathit{T}}{\mathbf{d}}_{\mathit{r}\mathbf{2}}⟩\\ & =c_{r21}⟨{\xi }_{r2n},{\gamma }_{r21}⟩+c_{r22}⟨{\xi }_{r2n},{\gamma }_{r22}⟩+\dots +c_{r2n}⟨{\xi }_{r2n},{\gamma }_{r2m}⟩+ \\ {d}_{r21}⟨{\xi }_{r2n},{\beta }_{r21}⟩+d_{r22}⟨{\xi }_{r2n},{\beta }_{r22}⟩+\dots +d_{r2n}⟨{\xi }_{r2n},{\beta }_{r2n}⟩.\end{array} \end{gathered}$$

Based on Equations (15)–(17), we obtain the smoothing spline component of the regression function $m_r$ of the MMNR model for $j=1$ and $i=\mathrm{1,2},\dots ,n$ as follows:

(22)$$\begin{gathered} \begin{array}{ll}{\mathbf{f}}_{\mathit{r}\mathbf{1}}=\left[\begin{array}{c}{f}_{r1}\left(x_{r1}\right)\\ f_{r1}\left(x_{r2}\right)\\ \begin{array}{c}⋮\\ f_{r1}\left(x_{rn}\right)\end{array}\end{array}\right]=& \left[\begin{array}{ccc}⟨{\xi }_{r11},{\gamma }_{r11}⟩& ⟨{\xi }_{r11},{\gamma }_{r12}⟩& \begin{array}{cc}\cdots & ⟨{\xi }_{r11},{\gamma }_{r1m}⟩\end{array}\\ ⟨{\xi }_{r12},{\gamma }_{r11}⟩& ⟨{\xi }_{r12},{\gamma }_{r12}⟩& \begin{array}{cc}\cdots & ⟨{\xi }_{r12},{\gamma }_{r1m}⟩\end{array}\\ \begin{array}{c}⋮\\ ⟨{\xi }_{r1n},{\gamma }_{r11}⟩\end{array}& \begin{array}{c}⋮\\ ⟨{\xi }_{r1n},{\gamma }_{r12}⟩\end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}⋮\\ ⟨{\xi }_{r1n},{\gamma }_{r1m}⟩\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{c}_{r11}\\ c_{r12}\\ \begin{array}{c}⋮\\ c_{r1n}\end{array}\end{array}\right]+ \\ \left[\begin{array}{ccc}⟨{\xi }_{r11},{\beta }_{r11}⟩& ⟨{\xi }_{r11},{\beta }_{r12}⟩& \begin{array}{cc}\cdots & ⟨{\xi }_{r11},{\beta }_{r1n}⟩\end{array}\\ ⟨{\xi }_{r12},{\beta }_{r11}⟩& ⟨{\xi }_{r12},{\beta }_{r12}⟩& \begin{array}{cc}\cdots & ⟨{\xi }_{r12},{\beta }_{r1n}⟩\end{array}\\ \begin{array}{c}⋮\\ ⟨{\xi }_{r1n},{\beta }_{r11}⟩\end{array}& \begin{array}{c}⋮\\ ⟨{\xi }_{r1n},{\beta }_{r12}⟩\end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}⋮\\ ⟨{\xi }_{r1n},{\beta }_{r1n}⟩\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{d}_{r11}\\ d_{r12}\\ \begin{array}{c}⋮\\ d_{r1n}\end{array}\end{array}\right]\\ {\mathbf{f}}_{\mathit{r}\mathbf{1}}={\mathbf{U}}_{\mathit{r}\mathbf{1}}{\mathbf{c}}_{\mathit{r}\mathbf{1}}+{\mathbf{V}}_{\mathit{r}\mathbf{1}}{\mathbf{d}}_{\mathit{r}\mathbf{1}}.\end{array} \end{gathered}$$

Next, in the same way as in the process of obtaining Equation (22), we obtain ${\mathbf{f}}_{\mathit{r}\mathit{j}}$ for $j=\mathrm{2,3},\dots ,p$ as follows:

(23)$\left.\begin{array}{c}{\mathbf{f}}_{\mathit{r}\mathbf{2}}={\mathbf{U}}_{\mathit{r}\mathbf{2}}{\mathbf{c}}_{\mathit{r}\mathbf{2}}+{\mathbf{V}}_{\mathit{r}\mathbf{2}}{\mathbf{d}}_{\mathit{r}\mathbf{2}}\\ {\mathbf{f}}_{\mathit{r}\mathbf{3}}={\mathbf{U}}_{\mathit{r}\mathbf{3}}{\mathbf{c}}_{\mathit{r}\mathbf{3}}+{\mathbf{V}}_{\mathit{r}\mathbf{3}}{\mathbf{d}}_{\mathit{r}\mathbf{3}}\\ \begin{array}{c}⋮\\ {\mathbf{f}}_{\mathit{r}\mathit{p}}={\mathbf{U}}_{\mathit{r}\mathit{p}}{\mathbf{c}}_{\mathit{r}\mathit{p}}+{\mathbf{V}}_{\mathit{r}\mathit{p}}{\mathbf{d}}_{\mathit{r}\mathit{p}}\end{array}\end{array}\right\}$

Since $⟨{\xi }_{rji},{\beta }_{rji}⟩=⟨{\beta }_{rji},{\beta }_{rji}⟩$, we have matrix ${\mathbf{V}}_{\mathit{r}\mathit{j}}$ as follows:

(24)${\mathbf{V}}_{\mathit{r}\mathit{j}}=\left[\begin{array}{ccc}⟨{\beta }_{rj1},{\beta }_{rj1}⟩& ⟨{\beta }_{rj1},{\beta }_{rj2}⟩& \begin{array}{cc}\cdots & ⟨{\beta }_{rj1},{\beta }_{rjn}⟩\end{array}\\ ⟨{\beta }_{rj2},{\beta }_{rj1}⟩& ⟨{\beta }_{rj2},{\beta }_{rj2}⟩& \begin{array}{cc}\cdots & ⟨{\beta }_{rj2},{\beta }_{rjn}⟩\end{array}\\ \begin{array}{c}⋮\\ ⟨{\beta }_{rjn},{\beta }_{rj1}⟩\end{array}& \begin{array}{c}⋮\\ ⟨{\beta }_{rjn},{\beta }_{rj2}⟩\end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}⋮\\ ⟨{\beta }_{rjn},{\beta }_{rjn}⟩\end{array}\end{array}\end{array}\right].$

Based on Equations (22) and (23), we obtain the smoothing spline component for the $r\mathrm{th}$ response, namely, ${\mathbf{f}}_{\mathit{r}}$, as follows:

(25)$\begin{array}{ll}{\mathbf{f}}_{\mathit{r}}=\left[\begin{array}{c}{\mathbf{f}}_{\mathit{r}\mathbf{1}}\\ {\mathbf{f}}_{\mathit{r}\mathbf{2}}\\ \begin{array}{c}⋮\\ {\mathbf{f}}_{\mathit{r}\mathit{p}}\end{array}\end{array}\right]& =\left[\begin{array}{c}{\mathbf{U}}_{\mathit{r}\mathbf{1}}{\mathbf{c}}_{\mathit{r}\mathbf{1}}+{\mathbf{V}}_{\mathit{r}\mathbf{1}}{\mathbf{d}}_{\mathit{r}\mathbf{1}}\\ {\mathbf{U}}_{\mathit{r}\mathbf{2}}{\mathbf{c}}_{\mathit{r}\mathbf{2}}+{\mathbf{V}}_{\mathit{r}\mathbf{2}}{\mathbf{d}}_{\mathit{r}\mathbf{2}}\\ \begin{array}{c}⋮\\ {\mathbf{U}}_{\mathit{r}\mathit{p}}{\mathbf{c}}_{\mathit{r}\mathit{p}}+{\mathbf{V}}_{\mathit{r}\mathit{p}}{\mathbf{d}}_{\mathit{r}\mathit{p}}\end{array}\end{array}\right]\\ & =\left[\begin{array}{ccc}{\mathbf{U}}_{\mathit{r}\mathbf{1}}& 0& \begin{array}{cc}\cdots & 0\end{array}\\ 0& {\mathbf{U}}_{\mathit{r}\mathbf{2}}& \begin{array}{cc}\cdots & 0\end{array}\\ \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}⋮\\ {\mathbf{U}}_{\mathit{r}\mathit{p}}\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{\mathbf{c}}_{\mathit{r}\mathbf{1}}\\ {\mathbf{c}}_{\mathit{r}\mathbf{2}}\\ \begin{array}{c}⋮\\ {\mathbf{c}}_{\mathit{r}\mathit{p}}\end{array}\end{array}\right]+\left[\begin{array}{ccc}{\mathbf{V}}_{\mathit{r}\mathbf{1}}& 0& \begin{array}{cc}\cdots & 0\end{array}\\ 0& {\mathbf{V}}_{\mathit{r}\mathbf{2}}& \begin{array}{cc}\cdots & 0\end{array}\\ \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}⋮\\ {\mathbf{V}}_{\mathit{r}\mathbf{2}}\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{\mathbf{d}}_{\mathit{r}\mathbf{1}}\\ {\mathbf{d}}_{\mathit{r}\mathbf{2}}\\ \begin{array}{c}⋮\\ {\mathbf{d}}_{\mathit{r}\mathit{p}}\end{array}\end{array}\right]\\ & ={\mathbf{U}}_{\mathit{r}}{\mathbf{c}}_{\mathit{r}}+{\mathbf{V}}_{\mathit{r}}{\mathbf{d}}_{\mathit{r}}.\end{array}$

Finally, based on Equation (25), we obtain the smoothing spline component for all responses $r=\mathrm{1,2},\dots , R$, namely, $\mathbf{f}$, as follows:

(26)$$\begin{gathered} \begin{array}{ll}\underset{(Rn\times 1)}{\underbrace{\mathbf{f}}}& =\underset{(Rn\times 1)}{\underbrace{{\left(f_{11},f_{12},\dots ,f_{1n},f_{21},f_{22},\dots ,f_{2n},\dots ,f_{R1},f_{R2},\dots ,f_{Rn}\right)}^T}}\\ & =\underset{(Rn\times Rp)}{\underbrace{diag\left(u_{11},u_{12},\dots ,u_{1p},u_{21},u_{22},\dots ,u_{2p},\dots ,u_{R1},u_{R2},\dots ,u_{Rp}\right)}}\underset{(Rp\times 1)}{\underbrace{{\left(c_{11},c_{12},\dots ,c_{1p},c_{21},c_{22},\dots ,c_{2p},\dots ,c_{R1},c_{R2},\dots ,c_{Rp}\right)}^T}}+ \\ \underset{(Rn\times Rp)}{\underbrace{diag\left(v_{11},v_{12},\dots ,v_{1p},v_{21},v_{22},\dots ,v_{2p},\dots ,v_{R1},v_{R2},\dots ,v_{Rp}\right)}}\underset{(Rp\times 1)}{\underbrace{{\left(d_{11},d_{12},\dots ,d_{1p},d_{21},d_{22},\dots ,d_{2p},\dots ,d_{R1},d_{R2},\dots ,d_{Rp}\right)}^T}}\\ \mathbf{f}=\mathbf{U}\mathbf{c}+\mathbf{V}\mathbf{d}\end{array} \end{gathered}$$

where

• $\mathbf{f}={\left(f_{11},f_{12},\dots ,f_{1n},f_{21},f_{22},\dots ,f_{2n},\dots ,f_{R1},f_{R2},\dots ,f_{Rn}\right)}^T$;

• $\mathbf{U}=diag\left(u_{11},u_{12},\dots ,u_{1p},u_{21},u_{22},\dots ,u_{2p},\dots ,u_{R1},u_{R2},\dots ,u_{Rp}\right)$;

• $\mathbf{c}={\left(c_{11},c_{12},\dots ,c_{1p},c_{21},c_{22},\dots ,c_{2p},\dots ,c_{R1},c_{R2},\dots ,c_{Rp}\right)}^T$;

• $\mathbf{V}=diag\left(v_{11},v_{12},\dots ,v_{1p},v_{21},v_{22},\dots ,v_{2p},\dots ,v_{R1},v_{R2},\dots ,v_{Rp}\right)$;

• $\mathbf{d}={\left(d_{11},d_{12},\dots ,d_{1p},d_{21},d_{22},\dots ,d_{2p},\dots ,d_{R1},d_{R2},\dots ,d_{Rp}\right)}^T$.

Thus, the smoothing spline component of the MMNR model can be presented in matrix notation as given in Equation (26), that is, $\mathbf{f}=\mathbf{U}\mathbf{c}+\mathbf{V}\mathbf{d}$.

3.2. Determining Fourier Series Component of MMNR Model

We consider the Fourier series component of the regression function of the MMNR model presented by Equation (3), namely, ${\sum }_{k=1}^q{g}_{rk}\left(t_{rki}\right)$. Note that the curve shape of $g_{rk}\left(t_{rki}\right)$ is unknown and only assumed to be contained in a continuous space $C\left(0,\pi \right)$. Based on Equation (4), this component is approximated by Fourier series function as follows:

(27)$\begin{array}{ll}{g}_{rk}\left(t_{rki}\right)& =b_{rk}{t}_{rki}+{\frac{1}{2}\alpha }_{0rk}+{\sum }_{v=1}^K{\alpha }_{vrk}\cos v{t}_{rki}; k=\mathrm{1,2},\dots ,q\\ & =b_{rk}{t}_{rki}+{\frac{1}{2}\alpha }_{0rk}+{\alpha }_{1rk}\cos t_{rki}+{\alpha }_{2rk}\cos 2t_{rki}+\dots +{\alpha }_{Krk}\cos K{t}_{rki}\end{array}$

Next, for $i=\mathrm{1,2},\dots ,n$, we may express the Fourier series function given in Equation (27) as follows:

(28)$\begin{array}{ll}{g}_{rk}\left(t_{rki}\right)& =\left[\begin{array}{ccc}{t}_{rk1}& 1/2& \begin{array}{ccc}\cos t_{rk1}& \cos 2t_{rk1}& \begin{array}{cc}\cdots & \cos K{t}_{rk1}\end{array}\end{array}\\ t_{rk2}& 1/2& \begin{array}{ccc}\cos t_{rk2}& \cos 2t_{rk2}& \begin{array}{cc}\cdots & \cos K{t}_{rk2}\end{array}\end{array}\\ \begin{array}{c}⋮\\ t_{rkn}\end{array}& \begin{array}{c}⋮\\ 1/2\end{array}& \begin{array}{c}\begin{array}{ccc}⋮ & ⋮& \begin{array}{cc}⋮ & ⋮\end{array}\end{array}\\ \begin{array}{ccc}\cos t_{rkn}& \cos 2t_{rkn}& \begin{array}{cc}\cdots & \cos K{t}_{rkn}\end{array}\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{b}_{rk}\\ {\alpha }_{0rk}\\ \begin{array}{c}{\alpha }_{1rk}\\ {\alpha }_{2rk}\\ \begin{array}{c}⋮\\ {\alpha }_{Krk}\end{array}\end{array}\end{array}\right]; k=\mathrm{1,2},\dots ,q\\ & ={\mathbf{T}}_{\mathit{r}\mathit{k}}{\mathsf{\alpha }}_{\mathit{r}\mathit{k}}\end{array}$

Hence, based on Equation (28), we can write the Fourier series component of the regression function of the MMNR model for $i=\mathrm{1,2},\dots ,n$ as follows:

(29)$\begin{array}{ll}{\mathbf{g}}_{\mathit{r}}& ={\sum }_{k=1}^q{g}_{rk}\left(t_{rki}\right)={\sum }_{k=1}^q{\mathbf{T}}_{\mathit{r}\mathit{k}}{\mathsf{\alpha }}_{\mathit{r}\mathit{k}}={\mathbf{T}}_{\mathit{r}\mathbf{1}}{\mathsf{\alpha }}_{\mathit{r}\mathbf{1}}+{\mathbf{T}}_{\mathit{r}\mathbf{2}}{\mathsf{\alpha }}_{\mathit{r}\mathbf{2}}+\dots +{\mathbf{T}}_{\mathit{r}\mathit{q}}{\mathsf{\alpha }}_{\mathit{r}\mathit{q}}\\ & ={\mathbf{T}}_{\mathit{r}}{\mathsf{\alpha }}_{\mathit{r}}\end{array}$

where ${\mathsf{\alpha }}_{\mathit{r}}={\left[b_{r1},{\alpha }_{0r1},{\alpha }_{1r1},\dots ,{\alpha }_{Kr1},\dots ,b_{rq},{\alpha }_{0rq},{\alpha }_{1rq},\dots ,{\alpha }_{Krq}\right]}^T$, $r=\mathrm{1,2},\dots ,R$, and

${\mathbf{T}}_{\mathit{r}}=\left[\begin{array}{ccc}{t}_{r11}& 1/2& \begin{array}{ccc}\cos t_{r11}& \cdots & \begin{array}{ccc}\cos K{t}_{r11}& \cdots & \begin{array}{ccc}{t}_{rq1}& 1/2& \begin{array}{ccc}\cos t_{rq1}& \dots & \cos K{t}_{rq1}\end{array}\end{array}\end{array}\end{array}\\ t_{r12}& 1/2& \begin{array}{ccc}\cos t_{r12}& \cdots & \begin{array}{ccc}\cos K{t}_{r12}& \cdots & \begin{array}{ccc}{t}_{rq2}& 1/2& \begin{array}{ccc}\cos t_{rq2}& \dots & \cos K{t}_{rq2}\end{array}\end{array}\end{array}\end{array}\\ \begin{array}{c}⋮\\ t_{r1n}\end{array}& \begin{array}{c}⋮\\ 1/2\end{array}& \begin{array}{c}\begin{array}{ccc}⋮ & ⋮ & \begin{array}{ccc}⋮& ⋮ & \begin{array}{ccc}⋮ & ⋮ & \begin{array}{ccc}⋮ & ⋮ & ⋮\end{array}\end{array}\end{array}\end{array}\\ \begin{array}{ccc}\cos t_{r1n}& \cdots & \begin{array}{ccc}\cos K{t}_{r1n}& \cdots & \begin{array}{ccc}{t}_{rqn}& 1/2& \begin{array}{ccc}\cos t_{rqn}& \cdots & \cos K{t}_{rqn}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]$

Furthermore, based on Equation (29) and in the same way as in the process of obtaining Equation (26), we obtain the Fourier series component for all responses $r=\mathrm{1,2},\dots , R$, namely, $\mathbf{g}$, as follows:

(30)$\begin{array}{ll}\underset{(Rn\times 1)}{\underbrace{\mathbf{g}}}& =\underset{(Rn\times 1)}{\underbrace{{\left(g_{11},g_{12},\dots ,g_{1n},g_{21},g_{22},\dots ,g_{2n},\dots ,g_{R1},g_{R2},\dots ,g_{Rn}\right)}^T}}\\ & =\underset{(Rn\times Rq)}{\underbrace{diag\left({\mathbf{T}}_{11},{\mathbf{T}}_{12},\dots ,{\mathbf{T}}_{1q},{\mathbf{T}}_{21},{\mathbf{T}}_{22},\dots ,{\mathbf{T}}_{2q},\dots ,{\mathbf{T}}_{R1}, {\mathbf{T}}_{R2},\dots ,{\mathbf{T}}_{Rq}\right)}}\underset{(Rq\times 1)}{\underbrace{{\left({\mathsf{\alpha }}_{11}^T,{\mathsf{\alpha }}_{12}^T,\dots ,{\mathsf{\alpha }}_{1q}^T,{\mathsf{\alpha }}_{21}^T,{\mathsf{\alpha }}_{22}^T,\dots ,{\mathsf{\alpha }}_{2q}^T,\dots ,{\mathsf{\alpha }}_{R1}^T, {\mathsf{\alpha }}_{R2}^T,\dots ,{\mathsf{\alpha }}_{Rq}^T\right)}^T}}\\ & =\mathbf{T}\mathsf{\alpha }\end{array}$

where

• $\mathbf{g}={\left(g_{11},g_{12},\dots ,g_{1n},g_{21},g_{22},\dots ,g_{2n},\dots ,g_{R1},g_{R2},\dots ,g_{Rn}\right)}^T$;

• $\mathbf{T}=diag\left({\mathbf{T}}_{11},{\mathbf{T}}_{12},\dots ,{\mathbf{T}}_{1q},{\mathbf{T}}_{21},{\mathbf{T}}_{22},\dots ,{\mathbf{T}}_{2q},\dots ,{\mathbf{T}}_{R1}, {\mathbf{T}}_{R2},\dots ,{\mathbf{T}}_{Rq}\right)$;

• $\mathsf{\alpha }={\left({\mathsf{\alpha }}_{11}^T,{\mathsf{\alpha }}_{12}^T,\dots ,{\mathsf{\alpha }}_{1q}^T,{\mathsf{\alpha }}_{21}^T,{\mathsf{\alpha }}_{22}^T,\dots ,{\mathsf{\alpha }}_{2q}^T,\dots ,{\mathsf{\alpha }}_{R1}^T, {\mathsf{\alpha }}_{R2}^T,\dots ,{\mathsf{\alpha }}_{Rq}^T\right)}^T$.

3.3. Determining Goodness of Fit and Penalty Components of PWLS Optimization

To estimate the regression function of the MMNR model, we use the PWLS as presented in Equation (5), namely, by taking the solution to the following PWLS optimization:

$$\begin{gathered} \underset{f_{\mathit{rj}}\in W_2^m[a_j,b_j]}{\mathrm{Min}}\left\{n^{-1}{\sum }_{i=1}^n{w}_{ri}\left(y_{ri}-\left[{\sum }_{j=1}^p{f}_{rj}\left(x_{rji}\right)\right.\right.\right.+{\sum }_{k=1}^q\left(b_{rk}{t}_{rki}+\frac{1}{2}{\alpha }_{0rk}+\right. \\ {\left.\left.\left.{\sum }_{v=1}^K{\alpha }_{vrk}\cos v{t}_{rki}\right)\right]\right)}^2+\left.{\sum }_{j=1}^p{\lambda }_{rj}{\int }_{a_j}^{b_j}{\left[f_{rj}\left(x_{rj}\right)\right]}^{2}d{x}_{rj}\right\} \end{gathered}$$

where $0<{\lambda }_{rj}<\infty$; $r=\mathrm{1,2},\dots ,R$; $j=\mathrm{1,2},\dots ,p$; and $w_{ri}$ represents a weight that is an inverse of random error variances.

Based on Equations (8) and (9), we have the MMNR model, which can be expressed as $\mathbf{y}=\mathbf{f}+\mathbf{g}+\mathsf{\epsilon }$ or $\mathbf{z}=\mathbf{f}+\mathsf{\epsilon }$ where z $=\mathbf{y}-\mathbf{g}$. Therefore, by considering Equation (26), the PWLS optimization can be written in matrix notation as follows:

(31)$\mathbf{Q}(\mathbf{c},\mathbf{d})=n^{-1}{\left(\mathbf{z}-\mathbf{f}\right)}^T\mathbf{W}\left(\mathbf{z}-\mathbf{f}\right)+{\sum }_{j=1}^p{\lambda }_{rj}{\int }_{a_j}^{b_j}{\left[f_{rj}\left(x_{rj}\right)\right]}^{2}d{x}_{rj}$

where $\mathbf{W}$ is a symmetrical weight matrix that is the inverse of a covariance matrix of random errors. Details of the weight matrix can be found in Lestari et al. [18,35]. In addition, the PWLS optimization function presented in Equation (31) shows that the goodness of fit component of PWLS optimization is $n^{-1}{\left(\mathbf{z}-\mathbf{f}\right)}^T\mathbf{W}\left(\mathbf{z}-\mathbf{f}\right)$.

Next, we determine the penalty component of PWLS optimization. In the PWLS optimization function presented in Equation (31), ${\sum }_{j=1}^p{\lambda }_{rj}{\int }_{a_j}^{b_j}{\left[f_{rj}\left(x_{rj}\right)\right]}^{2}d{x}_{rj}$ is a penalty component of PWLS optimization. In this step, we also express ${\int }_{a_j}^{b_j}{\left[f_{rj}\left(x_{rj}\right)\right]}^{2}d{x}_{rj}$ in matrix notation as follows:

(32)$\begin{array}{ll}{\int }_{a_j}^{b_j}{\left[f_{rj}\left(x_{rj}\right)\right]}^{2}d{x}_{rj}& ={‖{\mathbf{P}}_{\mathbf{1}}{\mathbf{f}}_{\mathit{r}\mathit{j}}‖}^2=⟨{\mathbf{P}}_{\mathbf{1}}{\mathbf{f}}_{\mathit{r}\mathit{j}},{\mathbf{P}}_{\mathbf{1}}{\mathbf{f}}_{\mathit{r}\mathit{j}}⟩\\ & =⟨{\mathbf{P}}_{\mathbf{1}}\left({\mathsf{\gamma }}_{\mathit{r}\mathit{j}}^{\mathit{T}}{\mathbf{c}}_{\mathit{r}\mathit{j}}+{\mathsf{\beta }}_{\mathit{r}\mathit{j}}^{\mathit{T}}{\mathbf{d}}_{\mathit{r}\mathit{j}}\right),{\mathbf{P}}_{\mathbf{1}}\left({\mathsf{\gamma }}_{\mathit{r}\mathit{j}}^{\mathit{T}}{\mathbf{c}}_{\mathit{r}\mathit{j}}+{\mathsf{\beta }}_{\mathit{r}\mathit{j}}^{\mathit{T}}{\mathbf{d}}_{\mathit{r}\mathit{j}}\right)⟩\\ & ={\mathbf{d}}_{\mathit{r}\mathit{j}}^{\mathit{T}}\left({\mathsf{\beta }}_{\mathit{r}\mathit{j}}^{\mathit{T}}{\mathsf{\beta }}_{\mathit{r}\mathit{j}}^{\mathit{T}}\right){\mathbf{d}}_{\mathit{r}\mathit{j}}\\ & =\left[\begin{array}{ccc}{d}_{rj1}& d_{rj2}& \begin{array}{cc}\cdots & d_{rjn}\end{array}\end{array}\right]\left[\begin{array}{ccc}⟨{\beta }_{rj1},{\beta }_{rj1}⟩& ⟨{\beta }_{rj1},{\beta }_{rj2}⟩& \begin{array}{cc}\cdots & ⟨{\beta }_{rj1},{\beta }_{rjn}⟩\end{array}\\ ⟨{\beta }_{rj2},{\beta }_{rj1}⟩& ⟨{\beta }_{rj2},{\beta }_{rj2}⟩& \begin{array}{cc}\cdots & ⟨{\beta }_{rj2},{\beta }_{rjn}⟩\end{array}\\ \begin{array}{c}⋮\\ ⟨{\beta }_{rjn},{\beta }_{rj1}⟩\end{array}& \begin{array}{c}⋮\\ ⟨{\beta }_{rjn},{\beta }_{rj2}⟩\end{array}& \begin{array}{cc}\begin{array}{c}⋮\\ \cdots \end{array}& \begin{array}{c}⋮\\ ⟨{\beta }_{rjn},{\beta }_{rjn}⟩\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{d}_{rj1}\\ d_{rj2}\\ \begin{array}{c}⋮\\ d_{rjn}\end{array}\end{array}\right]\\ & ={\mathbf{d}}_{\mathit{r}\mathit{j}}^{\mathit{T}} {\mathbf{V}}_{\mathit{r}\mathit{j}} {\mathbf{d}}_{\mathit{r}\mathit{j}}.\end{array}$

Hence, based on Equation (32), we obtain the penalty component of PWLS optimization, which can be expressed in matrix notation for the $r\mathrm{th}$ response as follows:

(33)$\begin{array}{ll}{\sum }_{j=1}^p{\lambda }_{rj}{\int }_{a_j}^{b_j}{\left[f_{rj}\left(x_{rj}\right)\right]}^{2}d{x}_{rj}& ={\sum }_{j=1}^p{\mathit{\lambda }}_{\mathit{r}\mathit{j}} {\mathbf{d}}_{\mathit{r}\mathit{j}}^{\mathit{T}} {\mathbf{V}}_{\mathit{r}\mathit{j}} {\mathbf{d}}_{\mathit{r}\mathit{j}}\\ & =\left[\begin{array}{ccc}{\mathit{d}}_{\mathit{r}\mathbf{1}}& {\mathit{d}}_{\mathit{r}\mathbf{2}}& \begin{array}{cc}\cdots & {\mathit{d}}_{\mathit{r}\mathit{p}}\end{array}\end{array}\right]\left[\begin{array}{ccc}{\mathbf{V}}_{\mathit{r}\mathbf{1}}{\mathit{\lambda }}_{\mathit{r}\mathbf{1}}& 0& \begin{array}{cc}\cdots & 0\end{array}\\ 0& {\mathbf{V}}_{\mathit{r}\mathbf{2}}{\mathit{\lambda }}_{\mathit{r}\mathbf{2}}& \begin{array}{cc}\cdots & 0\end{array}\\ \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}⋮\\ {\mathbf{V}}_{\mathit{r}\mathit{p}}{\mathit{\lambda }}_{\mathit{r}\mathit{p}}\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{\mathit{d}}_{\mathit{r}\mathbf{1}}\\ {\mathit{d}}_{\mathit{r}\mathbf{2}}\\ \begin{array}{c}⋮\\ {\mathit{d}}_{\mathit{r}\mathit{p}}\end{array}\end{array}\right]\\ & =\left[\begin{array}{ccc}{\mathit{d}}_{\mathit{r}\mathbf{1}}& {\mathit{d}}_{\mathit{r}\mathbf{2}}& \begin{array}{cc}\cdots & {\mathit{d}}_{\mathit{r}\mathit{p}}\end{array}\end{array}\right]\left[\begin{array}{ccc}{\mathit{\lambda }}_{\mathit{r}\mathbf{1}}\mathbf{I}& 0& \begin{array}{cc}\cdots & 0\end{array}\\ 0& {\mathit{\lambda }}_{\mathit{r}\mathbf{2}}\mathbf{I}& \begin{array}{cc}\cdots & 0\end{array}\\ \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}⋮\\ {\mathit{\lambda }}_{\mathit{r}\mathit{p}}\mathbf{I}\end{array}\end{array}\end{array}\right]\left[\begin{array}{ccc}{\mathbf{V}}_{\mathit{r}\mathbf{1}}& 0& \begin{array}{cc}\cdots & 0\end{array}\\ 0& {\mathbf{V}}_{\mathit{r}\mathbf{2}}& \begin{array}{cc}\cdots & 0\end{array}\\ \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{c}⋮\\ 0\end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}⋮\\ {\mathbf{V}}_{\mathit{r}\mathit{p}}\end{array}\end{array}\end{array}\right]\left[\begin{array}{c}{\mathit{d}}_{\mathit{r}\mathbf{1}}\\ {\mathit{d}}_{\mathit{r}\mathbf{2}}\\ \begin{array}{c}⋮\\ {\mathit{d}}_{\mathit{r}\mathit{p}}\end{array}\end{array}\right]\\ & ={\mathbf{d}}_{\mathit{r}}^T {\mathbf{L}}_{\mathbf{r}} {\mathbf{V}}_{\mathit{r}}{\mathbf{d}}_{\mathit{r}}.\end{array}$

Furthermore, based on Equation (33) and in the same way as in the process of obtaining Equation (26), the penalty component for all responses $r=\mathrm{1,2},\dots ,R$ can be obtained as follows:

Let ${\mathbf{P}}_{\mathit{r}}={\sum }_{j=1}^p{\lambda }_{rj}{\int }_{a_j}^{b_j}{\left[f_{rj}\left(x_{rj}\right)\right]}^{2}d{x}_{rj}={\sum }_{j=1}^p{\mathit{\lambda }}_{\mathit{r}\mathit{j}}{\mathbf{d}}_{\mathit{r}\mathit{j}}^{\mathit{T}}{\mathbf{V}}_{\mathit{r}\mathit{j}}{\mathbf{d}}_{\mathit{r}\mathit{j}}$ be the $r\mathrm{th}$ response penalty, and the penalty component for all responses $r=\mathrm{1,2},\dots ,R$ is presented by $\mathbf{P}$ as follows:

(34)$\begin{array}{ll}\underset{(Rn\times 1)}{\underbrace{\mathbf{P}}}=& \underset{(Rn\times Rp)}{\underbrace{{\left({\mathbf{d}}_{11},{\mathbf{d}}_{21},\dots ,{\mathbf{d}}_{\mathit{R}\mathbf{1}},{\mathbf{d}}_{12},{\mathbf{d}}_{22},\dots ,{\mathbf{d}}_{\mathit{R}\mathbf{2}},\dots ,{\mathbf{d}}_{\mathbf{1}\mathit{p}},{\mathbf{d}}_{\mathbf{2}\mathit{p}},\dots ,{\mathbf{d}}_{\mathit{R}\mathit{p}}\right)}^T}}\times \\ & \underset{(Rp\times Rp)}{\underbrace{diag\left({\mathsf{\lambda }}_{11}\mathbf{I},{\lambda }_{21}\mathbf{I},\dots ,{\lambda }_{\mathit{R}\mathbf{1}}\mathbf{I},{\mathsf{\lambda }}_{12}\mathbf{I},{\mathsf{\lambda }}_{22}\mathbf{I},\dots ,{\lambda }_{\mathit{R}\mathbf{2}}\mathbf{I},\dots ,{\lambda }_{\mathbf{1}\mathit{p}}\mathbf{I},{\lambda }_{\mathbf{2}\mathit{p}}\mathbf{I},\dots ,{\lambda }_{\mathit{R}\mathit{p}}\mathbf{I}\right)}}\times \\ & \underset{(Rp\times Rp)}{\underbrace{diag\left({\mathbf{V}}_{11},{\mathbf{V}}_{21},\dots ,{\mathbf{V}}_{\mathit{R}\mathbf{1}},{\mathbf{V}}_{12},{\mathbf{V}}_{22},\dots ,{\mathbf{V}}_{\mathit{R}\mathbf{2}},\dots ,{\mathbf{V}}_{\mathbf{1}\mathit{p}},{\mathbf{V}}_{\mathbf{2}\mathit{p}},\dots ,{\mathbf{V}}_{\mathit{R}\mathit{p}}\right)\times }}\\ & \underset{(Rp\times 1)}{\underbrace{\left({\mathbf{d}}_{11},{\mathbf{d}}_{21},\dots ,{\mathbf{d}}_{\mathit{R}\mathbf{1}},{\mathbf{d}}_{12},{\mathbf{d}}_{22},\dots ,{\mathbf{d}}_{\mathit{R}\mathbf{2}},\dots ,{\mathbf{d}}_{\mathbf{1}\mathit{p}},{\mathbf{d}}_{\mathbf{2}\mathit{p}},\dots ,{\mathbf{d}}_{\mathit{R}\mathit{p}}\right)}}\\ & ={\mathbf{d}}^{\mathit{T}}\mathbf{L}\mathbf{V}\mathbf{d}\end{array}$

where

• $\mathbf{d}=\left({\mathbf{d}}_{11},{\mathbf{d}}_{21},\dots ,{\mathbf{d}}_{\mathit{R}\mathbf{1}},{\mathbf{d}}_{12},{\mathbf{d}}_{22},\dots ,{\mathbf{d}}_{\mathit{R}\mathbf{2}},\dots ,{\mathbf{d}}_{\mathbf{1}\mathit{p}},{\mathbf{d}}_{\mathbf{2}\mathit{p}},\dots ,{\mathbf{d}}_{\mathit{R}\mathit{p}}\right)$;

• $\mathbf{L}=diag\left({\mathsf{\lambda }}_{11}\mathbf{I},{\lambda }_{21}\mathbf{I},\dots ,{\lambda }_{\mathit{R}\mathbf{1}}\mathbf{I},{\mathsf{\lambda }}_{12}\mathbf{I},{\mathsf{\lambda }}_{22}\mathbf{I},\dots ,{\lambda }_{\mathit{R}\mathbf{2}}\mathbf{I},\dots ,{\lambda }_{\mathbf{1}\mathit{p}}\mathbf{I},{\lambda }_{\mathbf{2}\mathit{p}}\mathbf{I},\dots ,{\lambda }_{\mathit{R}\mathit{p}}\mathbf{I}\right)$;

• $\mathbf{V}=diag\left({\mathbf{V}}_{11},{\mathbf{V}}_{21},\dots ,{\mathbf{V}}_{\mathit{R}1},{\mathbf{V}}_{12},{\mathbf{V}}_{22},\dots ,{\mathbf{V}}_{\mathit{R}2},\dots ,{\mathbf{V}}_{1\mathit{p}},{\mathbf{V}}_{2\mathit{p}},\dots ,{\mathbf{V}}_{\mathit{R}\mathit{p}}\right)$.

Thus, the goodness of fit component of the PWLS optimization is presented as follows:

$\begin{array}{l}{n}^{-1}{\sum }_{i=1}^n{w}_{ri}\left(y_{ri}-\left[{\sum }_{j=1}^p{f}_{rj}\left(x_{rji}\right)+{\sum }_{k=1}^q\left(b_{rk}{t}_{rki}+\frac{1}{2}{\alpha }_{0rk}+\right.\right.\right.{\left.\left.\left.{\sum }_{v=1}^K{\alpha }_{vrk}\cos v{t}_{rki}\right)\right]\right)}^2\\ \hspace{1em}=n^{-1}{\left(\mathbf{z}-\mathbf{f}\right)}^T\mathbf{W}\left(\mathbf{z}-\mathbf{f}\right)\end{array}$

and the penalty component of the PWLS optimization is presented as follows:

${\sum }_{j=1}^p{\lambda }_{rj}{\int }_{a_j}^{b_j}{\left[f_{rj}\left(x_{rj}\right)\right]}^{2}d{x}_{rj}={\mathbf{d}}^{\mathit{T}}\mathbf{L}\mathbf{V}\mathbf{d}.$

3.4. Estimating the MMNR Model

The PWLS optimization presented in Equation (31) is used to obtain the estimation of the MMNR model. Also, we definitely consider Equations (26), (30) and (34). In this section, we provide a complete explanation of the estimation process of the MMNR model. In Equation (31), we have the PWLS as follows:

$\begin{array}{ll}\mathbf{Q}(\mathbf{c},\mathbf{d})& =n^{-1}{\left(\mathbf{z}-\mathbf{f}\right)}^{\mathbf{T}}\mathbf{W}\left(\mathbf{z}-\mathbf{f}\right)+{\sum }_{j=1}^p{\lambda }_{rj}{\int }_{a_j}^{b_j}{\left[f_{rj}\left(x_{rj}\right)\right]}^{2}d{x}_{rj}\\ & =n^{-1}{\left(\mathbf{z}-\mathbf{f}\right)}^{\mathbf{T}}\mathbf{W}\left(\mathbf{z}-\mathbf{f}\right)+{\mathbf{d}}^{\mathbf{T}}\mathbf{L}\mathbf{V}\mathbf{d}\\ & =n^{-1}\left[{\left(\mathbf{z}-\left(\mathbf{U}\mathbf{c}+\mathbf{V}\mathbf{d}\right)\right)}^{\mathbf{T}}\mathbf{W}\left(\mathbf{z}-\left(\mathbf{U}\mathbf{c}+\mathbf{V}\mathbf{d}\right)\right)\right]+{\mathbf{d}}^{\mathbf{T}}\mathbf{L}\mathbf{V}\mathbf{d}\\ & =n^{-1}\left[{\left(\mathbf{z}-\mathbf{U}\mathbf{c}-\mathbf{V}\mathbf{d}\right)}^{\mathbf{T}}\mathbf{W}\left(\mathbf{z}-\mathbf{U}\mathbf{c}-\mathbf{V}\mathbf{d}\right)\right]+{\mathbf{d}}^{\mathbf{T}}\mathbf{L}\mathbf{V}\mathbf{d}\\ & =n^{-1}\left\{\left[{\left(\mathbf{z}-\mathbf{U}\mathbf{c}-\mathbf{V}\mathbf{d}\right)}^{\mathbf{T}}\mathbf{W}\left(\mathbf{z}-\mathbf{U}\mathbf{c}-\mathbf{V}\mathbf{d}\right)\right]+n \mathbf{L}{\mathbf{d}}^{\mathbf{T}}\mathbf{V}\mathbf{d}\right\}\\ & =n^{-1}\{{\mathbf{z}}^{\mathbf{T}}\mathbf{W}\mathbf{z}-{\mathbf{z}}^{\mathbf{T}}\mathbf{W}\mathbf{U}\mathbf{c}-{\mathbf{z}}^{\mathbf{T}}\mathbf{W}\mathbf{V}\mathbf{d}-{\mathbf{c}}^{\mathbf{T}}{\mathbf{U}}^{\mathbf{T}}\mathbf{W}\mathbf{z}+{\mathbf{c}}^{\mathbf{T}}{\mathbf{U}}^{\mathbf{T}}\mathbf{W}\mathbf{U}\mathbf{c}+\\ & {\mathbf{c}}^{\mathbf{T}}{\mathbf{U}}^{\mathbf{T}}\mathbf{W}\mathbf{V}\mathbf{d}-{\mathbf{d}}^{\mathbf{T}}{\mathbf{V}}^{\mathbf{T}}\mathbf{W}\mathbf{z}+{\mathbf{d}}^{\mathbf{T}}{\mathbf{V}}^{\mathbf{T}}\mathbf{W}\mathbf{U}\mathbf{c}+{\mathbf{d}}^{\mathbf{T}}{\mathbf{V}}^{\mathbf{T}}\mathbf{W}\mathbf{V}\mathbf{d}+n\mathbf{L}{\mathbf{d}}^{\mathbf{T}}\mathbf{V}\mathbf{d}\}.\end{array}$

In this step, we first determine the estimator for $\mathbf{d}$, namely, $\widehat{\mathbf{d}}$, as follows:

(35)$\begin{array}{ll}\frac{\partial \mathbf{Q}(\mathbf{c},\mathbf{d})}{\partial \mathbf{d}}{|}_{\mathbf{d}=\widehat{\mathbf{d}}}=\mathbf{0}& \iff -2{\mathbf{V}}^{\mathbf{T}}\mathbf{W}\mathbf{z}+2{\mathbf{V}}^{\mathbf{T}}\mathbf{W}\mathbf{U}\mathbf{c}+2{\mathbf{V}}^{\mathbf{T}}\mathbf{W}\mathbf{V}\mathbf{d}+\mathbf{2}n\mathbf{L}\mathbf{V}\mathbf{d}=\mathbf{0}\\ & \iff -2{\mathbf{V}}^{\mathbf{T}}\left[-\mathbf{z}\mathbf{W}+\mathbf{U}\mathbf{W}\mathbf{c}+\left(\mathbf{W}\mathbf{V}+n\mathbf{L}\right)\widehat{\mathbf{d}}\right]=\mathbf{0}\\ & \iff -\mathbf{z}\mathbf{W}+\mathbf{U}\mathbf{W}\mathbf{c}+\left(\mathbf{W}\mathbf{V}+n\mathbf{L}\right)\widehat{\mathbf{d}}=\mathbf{0}.\end{array}$

Let $\mathbf{H}=(\mathbf{W}\mathbf{V}+n\mathbf{L})$ such that Equation (35) can be written as follows:

(36)$-\mathbf{z}\mathbf{W}+\mathbf{U}\mathbf{W}\mathbf{c}+\mathbf{H}\widehat{\mathbf{d}}=\mathbf{0}\iff \widehat{\mathbf{d}}={\mathbf{H}}^{-\mathbf{1}}\left(\mathbf{z}\mathbf{W}-\mathbf{U}\mathbf{W}\mathbf{c}\right).$

Next, we determine the estimator for $\mathbf{c}$, namely, $\widehat{\mathbf{c}}$, as follows:

(37)$\frac{\partial \mathbf{Q}(\mathbf{c},\mathbf{d})}{\partial \mathbf{c}}{|}_{\mathbf{c}=\widehat{\mathbf{c}}}=\mathbf{0} \iff -\mathbf{2}{\mathbf{U}}^{\mathbf{T}}\mathbf{W}\mathbf{z}+\mathbf{2}{\mathbf{U}}^{\mathbf{T}}\mathbf{W}\mathbf{U}\widehat{\mathbf{c}}+\mathbf{2}{\mathbf{U}}^{\mathbf{T}}\mathbf{W}\mathbf{V}\mathbf{d}=\mathbf{0}.$

By substituting Equation (36) into Equation (37), we obtain the following equation:

(38)$-\mathbf{2}{\mathbf{U}}^{\mathbf{T}}\mathbf{W}\mathbf{z}+\mathbf{2}{\mathbf{U}}^{\mathbf{T}}\mathbf{W}\mathbf{U}\widehat{\mathbf{c}}+\mathbf{2}{\mathbf{U}}^{\mathbf{T}}\mathbf{W}\mathbf{V}\left[{\mathbf{H}}^{-\mathbf{1}}\left(\mathbf{z}\mathbf{W}-\mathbf{U}\mathbf{W}\widehat{\mathbf{c}}\right)\right]=\mathbf{0}.$

Since $\mathbf{H}=(\mathbf{W}\mathbf{V}+n\mathbf{L})$, then we have $\mathbf{V}{\mathbf{H}}^{-\mathbf{1}}=(\mathbf{I}-n\mathbf{L}{\mathbf{H}}^{-\mathbf{1}})$, and we may therefore write Equation (38) as follows:

(39)$-\mathbf{2}{\mathbf{U}}^{\mathbf{T}}\mathbf{W}\mathbf{z}+\mathbf{2}{\mathbf{U}}^{\mathbf{T}}\mathbf{W}\mathbf{U}\widehat{\mathbf{c}}+\mathbf{2}{\mathbf{U}}^{\mathbf{T}}\mathbf{W}\left(\mathbf{I}-n\mathbf{L}{\mathbf{H}}^{-\mathbf{1}}\right)\left(\mathbf{z}\mathbf{W}-\mathbf{U}\mathbf{W}\widehat{\mathbf{c}}\right)=\mathbf{0}.$

Equation (39) returns the following result:

(40)$\widehat{\mathbf{c}}={\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{z}.$

By substituting Equation (40) into Equation (36), we obtain the following result:

(41)$\widehat{\mathbf{d}}={\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\left[\mathbf{I}-\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\right]\mathbf{z}.$

Hence, based on Equations (26), (40) and (41), we obtain the following estimated smoothing spline component:

$\begin{array}{ll}\widehat{\mathbf{f}}& =\mathbf{U}\widehat{\mathbf{c}}+\mathbf{V}\widehat{\mathbf{d}}\\ & =\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{z}+\mathbf{V}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\left(\mathbf{I}\mathbf{-}\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\right)\mathbf{z}\\ & =\left[\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}+\mathbf{V}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\left(\mathbf{I}-\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\right)\right]\mathbf{z}\\ & =\left[\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}+\mathbf{V}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\left(\mathbf{I}-\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\right)\right]\left(\mathbf{y}-\mathbf{g}\right)\\ & =\mathbf{A}(\mathsf{\lambda },\mathbf{K})\left(\mathbf{y}-\mathbf{g}\right).\end{array}$

Thus, the estimation of the smoothing spline component of the MMNR model notated as ${\widehat{\mathbf{f}}}_{\left(\mathit{\lambda },\mathit{K}\right)}(\mathbf{x},\mathbf{t})$ is given by

(42)${\widehat{\mathbf{f}}}_{\left(\mathit{\lambda },\mathit{K}\right)}\left(\mathbf{x},\mathbf{t}\right)=\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\left(\mathbf{y}-\mathbf{g}\right)=\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{g}.$

where $\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)=\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}+\mathbf{V}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\left(\mathbf{I}-\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\right)$, and $\mathbf{W}$ is a special symmetric weight matrix, namely, a diagonal matrix, which is an inverse of a covariance matrix of random errors. In this case, based on the assumption of random errors of the MMNR model, the weight matrix $\mathbf{W}$ is given as follows [18,35]:

$\mathbf{W}={\left[\mathit{d}\mathit{i}\mathit{a}\mathit{g}\left({\mathbf{W}}_{\mathbf{1}},{\mathbf{W}}_{\mathbf{2}},\dots ,{\mathbf{W}}_{\mathit{R}}\right)\right]}^{-\mathbf{1}}$

where

${\mathbf{W}}_{\mathbf{r}}=\left(\begin{array}{ccc}{\sigma }_{r1}^2& {\sigma }_{r\left(\mathrm{1,2}\right)}& \begin{array}{cc}\cdots & {\sigma }_{r(1,n)}\end{array}\\ {\sigma }_{r\left(\mathrm{2,1}\right)}& {\sigma }_{r2}^2& \begin{array}{cc}\cdots & {\sigma }_{r(2,n)}\end{array}\\ \begin{array}{c}⋮\\ {\sigma }_{r(n,1)}\end{array}& \begin{array}{c}⋮\\ {\sigma }_{r(n,2)}\end{array}& \begin{array}{cc}\begin{array}{c}\ddots \\ \cdots \end{array}& \begin{array}{c}⋮\\ {\sigma }_{rn}^2\end{array}\end{array}\end{array}\right), r=\mathrm{1,2},\dots ,R.$

In the previous explanation, we obtained ${\widehat{\mathbf{f}}}_{\left(\lambda ,K\right)}(\mathbf{x},\mathbf{t})$, which is presented in (42). The next step is determining the Fourier series component of the regression function of the MMNR model. Let us look again at the MMNR model presented in Equation (8):

$\mathbf{y}=\mathbf{f}+\mathbf{g}+\mathsf{\epsilon }$

We substitute Equation (42) into Equation (8) such that we obtain the following equation:

(43)$\begin{array}{cc}\hfill \mathbf{y}& =\mathbf{f}+\mathbf{g}+\mathsf{\epsilon }\hfill \\ & =\mathbf{A}(\mathsf{\lambda },\mathbf{K})\left(\mathbf{y}-\mathbf{g}\right)+\mathbf{g}+\mathsf{\epsilon }\hfill \\ \hfill \iff \mathsf{\epsilon }& =\mathbf{y}-\left[\mathbf{A}\right(\mathsf{\lambda },\mathbf{K})\left(\mathbf{y}-\mathbf{g}\right)+\mathbf{g}]\hfill \\ & =\mathbf{y}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}+\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{g}-\mathbf{g}.\hfill \end{array}$

Next, we substitute Equation (30), namely, $\mathbf{g}=\mathbf{T}\mathsf{\alpha }$, into Equation (43) such that we obtain the equation for random error as follows:

$\mathsf{\epsilon }=\mathbf{y}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}+\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}\mathsf{\alpha }-\mathbf{T}\mathsf{\alpha }.$

Hence, we have

(44)$\begin{array}{ll}{\mathsf{\epsilon }}^{\mathbf{T}}\mathsf{\epsilon }& ={\left[\mathbf{y}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}+\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}\mathsf{\alpha }-\mathbf{T}\mathsf{\alpha }\right]}^{\mathbf{T}}\left[\mathbf{y}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}+\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}\mathsf{\alpha }-\mathbf{T}\mathsf{\alpha }\right]\\ & ={\mathbf{y}}^{\mathbf{T}}\mathbf{y}-2{\mathbf{y}}^{\mathbf{T}}\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}+\mathbf{2}{\mathsf{\alpha }}^{\mathbf{T}}{\mathbf{T}}^{\mathbf{T}}\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}-\mathbf{2}{\mathsf{\alpha }}^{\mathbf{T}}{\mathbf{T}}^{\mathbf{T}}\mathbf{y}-\mathbf{2}{\mathsf{\alpha }}^{\mathbf{T}}{\mathbf{T}}^{\mathbf{T}}{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}+\\ & \hspace{1em}\mathbf{2}{\mathsf{\alpha }}^{\mathbf{T}}{\mathbf{T}}^{\mathbf{T}}\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}-\mathbf{2}{\mathsf{\alpha }}^{\mathbf{T}}{\mathbf{T}}^{\mathbf{T}}\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}\mathsf{\alpha }+{\mathbf{y}}^{\mathbf{T}}{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}+\\ & \hspace{1em}{\mathsf{\alpha }}^{\mathbf{T}}{\mathbf{T}}^{\mathbf{T}}{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}\mathsf{\alpha }+{\mathsf{\alpha }}^{\mathbf{T}}{\mathbf{T}}^{\mathbf{T}}\mathbf{T}\mathsf{\alpha }.\end{array}$

Furthermore, based on Equation (44), we determine the estimator for $\mathsf{\alpha }$, namely, $\widehat{\mathsf{\alpha }}$, as follows:

(45)$\begin{array}{l}\frac{\partial \left({\mathsf{\epsilon }}^{\mathbf{T}}\mathsf{\epsilon }\right)}{\partial \mathsf{\alpha }}{|}_{\mathsf{\alpha }=\widehat{\mathsf{\alpha }}}=\mathbf{0}\\ \\ \iff \mathbf{2}{\mathbf{T}}^{\mathbf{T}}\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}-\mathbf{2}{\mathbf{T}}^{\mathbf{T}}\mathbf{y}-\mathbf{2}{\mathbf{T}}^{\mathbf{T}}{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}+\mathbf{2}{\mathbf{T}}^{\mathbf{T}}\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}+\\ -\mathbf{4}{\mathbf{T}}^{\mathbf{T}}\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}\widehat{\mathsf{\alpha }}+\mathbf{2}{\mathbf{T}}^{\mathbf{T}}{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}\widehat{\mathsf{\alpha }}+\mathbf{2}{\mathbf{T}}^{\mathbf{T}}\mathbf{T}\widehat{\mathsf{\alpha }}=\mathbf{0}\\ \iff \widehat{\mathsf{\alpha }}={\left[\left(\mathbf{2}-{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\right)\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}-\mathbf{T}\right]}^{-\mathbf{1}}\left({\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}-\mathbf{I}\right)\left(\mathbf{I}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\mathbf{y}.\end{array}$

Additionally, by substituting Equation (45) into Equation (30), we obtain the estimator for the Fourier series component of the regression function of the MMNR model, namely, ${\widehat{\mathbf{g}}}_{\left(\lambda ,K\right)}(\mathbf{x},\mathbf{t})$, as follows:

(46)$\begin{array}{ll}{\widehat{\mathbf{g}}}_{\left(\mathit{\lambda },\mathit{K}\right)}\left(\mathbf{x},\mathbf{t}\right)& =\mathbf{T}{\left[\left(\mathbf{2}-{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\right)\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}-\mathbf{T}\right]}^{-\mathbf{1}}\left({\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}-\mathbf{I}\right)\left(\mathbf{I}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\mathbf{y}\\ & =\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}\end{array}$

where $\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)=\mathbf{T}{\left[\left(\mathbf{2}-{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\right)\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}-\mathbf{T}\right]}^{-\mathbf{1}}\left({\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}-\mathbf{I}\right)\left(\mathbf{I}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\right)$.

Hereinafter, by substituting Equation (46) into Equation (42) we obtain the estimator for the smoothing spline component of the regression function of the MMNR model as follows:

(47)$\begin{array}{ll}{\widehat{\mathbf{f}}}_{\left(\mathit{\lambda },\mathit{K}\right)}\left(\mathbf{x},\mathbf{t}\right)& =\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\left[\mathbf{I}-\mathbf{T}{\left(\left(\mathbf{2}-{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\right)\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}-\mathbf{T}\right)}^{-\mathbf{1}}\left({\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}-\mathbf{I}\right)\left(\mathbf{I}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\right]\mathbf{y}\\ & =\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}\end{array}$

where $\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)=\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\left[\mathbf{I}-\mathbf{T}{\left(\left(\mathbf{2}-{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\right)\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}-\mathbf{T}\right)}^{-\mathbf{1}}\left({\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}-\mathbf{I}\right)\left(\mathbf{I}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\right]$.

Finally, based on Equations (2), (8), (46) and (47), we obtain the estimation of the regression function of the MMNR model, which consists of estimations of the smoothing spline component and Fourier series component as follows:

${\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})={\widehat{\mathbf{f}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})+{\widehat{\mathbf{g}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t}).$

Thus, the estimation result of the MMNR model by using the mixed smoothing spline and Fourier series estimator is as follows:

(48)$\begin{array}{ll}\widehat{\mathbf{y}}={\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})& ={\widehat{\mathbf{f}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})+{\widehat{\mathbf{g}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})\\ & =\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{y}\\ & =\left(\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\mathbf{y}\end{array}$

where λ represents the smoothing parameter of the smoothing spline estimator, $\mathbf{K}$ represents the oscillation parameter of the Fourier series estimator,

$\begin{array}{l}\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)=\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\left[\mathbf{I}-\mathbf{T}{\left(\left(\mathbf{2}-{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\right)\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}-\mathbf{T}\right)}^{-\mathbf{1}}\left({\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}-\mathbf{I}\right)\left(\mathbf{I}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\right],\\ \mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)=\mathbf{T}{\left[\left(\mathbf{2}-{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\right)\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}-\mathbf{T}\right]}^{-\mathbf{1}}\left({\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}-\mathbf{I}\right)\left(\mathbf{I}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\right),\mathrm{and}\\ \mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)=\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}+\mathbf{V}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\left(\mathbf{I}-\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\right).\end{array}$

Next, in the following section, we discuss the estimation of weight matrix W.

3.5. Estimating Weight Matrix W

Suppose that we have a paired dataset $\left(x_{r1i},x_{r2i},\dots ,x_{rpi},t_{r1i},t_{r2i},\dots ,t_{rqi},y_{ri}\right)$ for $r=\mathrm{1,2},\dots ,R$ and $i=\mathrm{1,2},\dots ,n$, which follows the MMNR model, and we assume that $\mathbf{y}={\left[{\mathbf{y}}_{\mathbf{1}},{\mathbf{y}}_{\mathbf{2}},\dots ,{\mathbf{y}}_{\mathbf{R}}\right]}^{\mathrm{T}}$ is a normally distributed random sample with a mean of m and covariance of ${\mathbf{W}}^{-\mathbf{1}}$ such that we have a likelihood function as follows:

$L\left(\mathbf{m},\mathbf{W}|\mathbf{y}\right)={\prod }_{j=1}^n\left\{\frac{1}{{\left(2\pi \right)}^{\frac{n}{2}}{\left({\mathbf{W}}_r^{-1}\right)}^{\frac{1}{2}}}\exp \left(-\frac{1}{2}{\left({\mathbf{y}}_{rj}-{\mathbf{m}}_{rj}\right)}^{\mathbf{T}}{\mathbf{W}}_{\mathit{r}}({\mathbf{y}}_{rj}-{\mathbf{m}}_{rj})\right)\right\}.$

Since $\mathbf{W}=diag({\mathbf{W}}_{\mathbf{1}},{\mathbf{W}}_{\mathbf{2}},\dots ,{\mathbf{W}}_{\mathbf{R}})$, then we may write the likelihood function as follows:

$\begin{array}{ll}L\left(\mathbf{m},\mathbf{W}|\mathbf{y}\right)=& \frac{1}{{\left(2\pi \right)}^{\frac{n^2}{2}}{\left({\mathbf{W}}_1^{-1}\right)}^{\frac{n}{2}}}\mathit{exp}\left(-\frac{1}{2}{\sum }_{j=1}^n{\left({\mathbf{y}}_{\mathbf{1}\mathit{j}}-{\mathbf{m}}_{\mathbf{1}\mathbf{j}}\right)}^{\mathbf{T}}{\mathbf{W}}_{\mathbf{1}}({\mathbf{y}}_{\mathbf{1}\mathit{j}}-{\mathbf{m}}_{\mathbf{1}\mathbf{j}})\right)\times \\ & \frac{1}{{\left(2\pi \right)}^{\frac{n^2}{2}}{\left({\mathbf{W}}_2^{-1}\right)}^{\frac{n}{2}}}\mathit{exp}\left(-\frac{1}{2}{\sum }_{j=1}^n{\left({\mathbf{y}}_{\mathbf{2}\mathit{j}}-{\mathbf{m}}_{\mathbf{2}\mathbf{j}}\right)}^{\mathbf{T}}{\mathbf{W}}_{\mathbf{2}}({\mathbf{y}}_{\mathbf{2}\mathit{j}}-{\mathbf{m}}_{\mathbf{2}\mathbf{j}})\right)\times \dots \times \\ & \frac{1}{{\left(2\pi \right)}^{\frac{n^2}{2}}{\left({\mathbf{W}}_R^{-1}\right)}^{\frac{n}{2}}}\exp \left(-\frac{1}{2}{\sum }_{j=1}^n{\left({\mathbf{y}}_{\mathit{R}\mathit{j}}-{\mathbf{m}}_{\mathbf{r}\mathbf{j}}\right)}^{\mathbf{T}}{\mathbf{W}}_{\mathbf{R}}\left({\mathbf{y}}_{\mathit{R}\mathit{j}}-{\mathbf{m}}_{\mathbf{R}\mathbf{j}}\right)\right).\end{array}$

Hence, the estimated weight matrix can be obtained by determining the solution to the following optimization:

$\begin{array}{ll}\underset{\mathbf{W}}{\mathbf{Q}=\mathrm{Max}}L\left(\mathbf{m},\mathbf{W}|\mathbf{y}\right)=& \underset{{\mathbf{W}}_{\mathbf{1}}}{\mathrm{Max}}\left\{\frac{1}{{\left(2\pi \right)}^{\frac{n^2}{2}}{\left({\mathbf{W}}_1^{-1}\right)}^{\frac{n}{2}}}\mathit{exp}\left(-\frac{1}{2}{\sum }_{j=1}^n{\left({\mathbf{y}}_{\mathbf{1}\mathit{j}}-{\mathbf{m}}_{\mathbf{1}\mathbf{j}}\right)}^{\mathbf{T}}{\mathbf{W}}_{\mathbf{1}}({\mathbf{y}}_{\mathbf{1}\mathit{j}}-{\mathbf{m}}_{\mathbf{1}\mathbf{j}})\right)\right\}\times \\ & \underset{{\mathbf{W}}_{\mathbf{2}}}{\mathrm{Max}}\left\{\frac{1}{{\left(2\pi \right)}^{\frac{n^2}{2}}{\left({\mathbf{W}}_2^{-1}\right)}^{\frac{n}{2}}}\mathit{exp}\left(-\frac{1}{2}{\sum }_{j=1}^n{\left({\mathbf{y}}_{\mathbf{2}\mathit{j}}-{\mathbf{m}}_{\mathbf{2}\mathbf{j}}\right)}^{\mathbf{T}}{\mathbf{W}}_{\mathbf{2}}({\mathbf{y}}_{\mathbf{2}\mathit{j}}-{\mathbf{m}}_{\mathbf{2}\mathbf{j}})\right)\right\}\times \cdots \times \\ & \underset{{\mathbf{W}}_{\mathit{R}}}{\mathrm{Max}}\left\{\frac{1}{{\left(2\pi \right)}^{\frac{n^2}{2}}{\left({\mathbf{W}}_R^{-1}\right)}^{\frac{n}{2}}}\mathit{exp}\left(-\frac{1}{2}{\sum }_{j=1}^n{\left({\mathbf{y}}_{\mathit{R}\mathit{j}}-{\mathbf{m}}_{\mathbf{R}\mathbf{j}}\right)}^{\mathbf{T}}{\mathbf{W}}_{\mathbf{R}}\left({\mathbf{y}}_{\mathit{R}\mathit{j}}-{\mathbf{m}}_{\mathbf{R}\mathbf{j}}\right)\right)\right\}.\end{array}$

Hereinafter, according to Johnson and Wichern [46], we can determine the maximum value of each component of the likelihood function by taking the following equations:

${\widehat{\mathbf{W}}}_{\mathbf{1}}=\frac{{\widehat{\mathsf{\epsilon }}}_{\mathbf{1}}{\widehat{\mathsf{\epsilon }}}_{\mathbf{1}}^{\mathbf{T}}}{Rn}=\frac{\left({\mathbf{y}}_{\mathbf{1}}-{\widehat{\mathbf{m}}}_{\mathbf{1}}\right){\left({\mathbf{y}}_{\mathbf{1}}-{\widehat{\mathbf{m}}}_{\mathbf{1}}\right)}^{\mathbf{T}}}{Rn} , {\widehat{\mathbf{W}}}_{\mathbf{2}}=\frac{{\widehat{\mathsf{\epsilon }}}_{\mathbf{2}}{\widehat{\mathsf{\epsilon }}}_{\mathbf{2}}^{\mathbf{T}}}{Rn}=\frac{\left({\mathbf{y}}_{\mathbf{2}}-{\widehat{\mathbf{m}}}_{\mathbf{2}}\right){\left({\mathbf{y}}_{\mathbf{2}}-{\widehat{\mathbf{m}}}_{\mathbf{2}}\right)}^{\mathbf{T}}}{Rn} , \dots , {\widehat{\mathbf{W}}}_{\mathit{R}}=\frac{{\widehat{\mathsf{\epsilon }}}_{\mathbf{R}}{\widehat{\mathsf{\epsilon }}}_{\mathbf{R}}^{\mathbf{T}}}{Rn}=\frac{\left({\mathbf{y}}_{\mathit{R}}-{\widehat{\mathbf{m}}}_{\mathbf{R}}\right){\left({\mathbf{y}}_{\mathbf{R}}-{\widehat{\mathbf{m}}}_{\mathit{R}}\right)}^{\mathbf{T}}}{Rn}.$

Furthermore, based on Equation (48), we have

$\widehat{\mathbf{y}}={\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}\left(\mathbf{x},\mathbf{t}\right)=\mathbf{M}(\mathsf{\lambda },\mathbf{K})\mathbf{y}$

where $\mathbf{M}\left(\mathsf{\lambda },\mathbf{K}\right)=\left(\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\right)$. Hence, we obtain the maximum likelihood estimator for the weight matrix $\mathbf{W}$ as follows:

$\widehat{\mathbf{W}}=diag\left({\widehat{\mathbf{W}}}_{\mathbf{1}},{\widehat{\mathbf{W}}}_{\mathbf{2}},\dots ,{\widehat{\mathbf{W}}}_{\mathit{R}}\right)$

where ${\widehat{\mathbf{W}}}_{\mathbf{1}}=\frac{\left(\mathbf{I}-{\mathbf{M}}_{\mathbf{1}}(\mathsf{\lambda },\mathbf{K},{\widehat{\mathit{\sigma }}}_{\mathbf{1}}^{\mathbf{2}})\right){\mathbf{y}}_{\mathbf{1}}{\mathbf{y}}_{\mathbf{1}}^{\mathit{T}}{\left(\mathbf{I}-{\mathbf{M}}_{\mathbf{1}}(\mathsf{\lambda },\mathbf{K},{\widehat{\mathit{\sigma }}}_{\mathbf{1}}^{\mathbf{2}})\right)}^{\mathit{T}}}{Rn}$, ${\widehat{\mathbf{W}}}_{\mathbf{2}}=\frac{\left(\mathbf{I}-{\mathbf{M}}_{\mathbf{2}}(\mathsf{\lambda },\mathbf{K},{\widehat{\mathit{\sigma }}}_{\mathbf{2}}^{\mathbf{2}})\right){\mathbf{y}}_{\mathbf{2}}{\mathbf{y}}_{\mathbf{2}}^{\mathit{T}}{\left(\mathbf{I}-{\mathbf{M}}_{\mathbf{2}}(\mathsf{\lambda },\mathbf{K},{\widehat{\mathit{\sigma }}}_{\mathbf{2}}^{\mathbf{2}})\right)}^{\mathit{T}}}{\mathit{R}\mathit{n}}$, …, ${\widehat{\mathbf{W}}}_{\mathit{R}}=\frac{\left(\mathbf{I}-{\mathbf{M}}_{\mathit{R}}(\mathsf{\lambda },\mathbf{K},{\widehat{\mathit{\sigma }}}_{\mathit{R}}^{\mathbf{2}})\right){\mathbf{y}}_{\mathit{R}}{\mathbf{y}}_{\mathit{R}}^{\mathit{T}}{\left(\mathbf{I}-{\mathbf{M}}_{\mathit{R}}(\mathsf{\lambda },\mathbf{K},{\widehat{\mathit{\sigma }}}_{\mathit{R}}^{\mathbf{2}})\right)}^{\mathit{T}}}{\mathit{R}\mathit{n}}$.

This shows that the estimated weight matrix obtained is a symmetric matrix, especially a diagonal matrix whose main diagonal components are the estimated weight matrices of the first response, second response, and so on, up to the R-th response.

Next, in the following section, we provide discussion on selecting optimal smoothing and oscillation parameters in the MMNR model.

3.6. Selecting Optimal Smoothing and Oscillation Parameters in the MMNR Model

In statistical analysis using an MMNR model approach based on a mixed smoothing spline and Fourier series estimator, selecting the optimal smoothing parameter and oscillation parameter, namely, (${\mathit{\lambda }}_{Opt},{\mathbf{K}}_{opt}$), cannot be neglected, and it is crucial to good regression function fitting of the MMNR model based on the mixed smoothing spline and Fourier series estimator. In other words, the mixed smoothing spline and Fourier series estimator highly depends on both the smoothing parameter and oscillation parameter. Therefore, in this study, we also need to discuss selecting the optimal smoothing and oscillation parameters. There are some criteria used for selecting the smoothing parameter, including minimizing CV (cross-validation), GCV (generalized cross-validation), Mallows’ ${\mathrm{C}}_p$, and AIC (Akaike’s information criterion) [1,12]. However, according to Ruppert and Carroll [47], for good fitting of a regression function based on a spline estimator, Mallows’ ${\mathrm{C}}_p$ and GCV were very satisfactory.

In this study, we use the GCV criterion, which has been developed for multiresponse cases to determine the optimal smoothing parameter and oscillation parameter values for good regression function fitting of the multiresponse multipredictor nonparametric regression (MMNR) model. Firstly, we determine the mean squared errors (MSEs) of the regression function presented in (48) as follows:

(49)$\begin{array}{ll}MSE\left(\mathsf{\lambda },\mathrm{K}\right)& =\frac{{\left[\mathbf{y}-\widehat{\mathbf{y}}\right]}^{\mathbf{T}}\mathbf{W}\left[\mathbf{y}-\widehat{\mathbf{y}}\right]}{n}\\ & =n^{-1}\left\{{\left[\mathbf{y}-{\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})\right]}^{\mathbf{T}}\mathbf{W}\left[\mathbf{y}-{\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})\right]\right\}\\ & =n^{-1}\left\{{\left[\mathbf{y}-\left(\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\mathbf{y}\right]}^{\mathbf{T}}\mathbf{W}\left[\mathbf{y}-\left(\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\mathbf{y}\right]\right\}\\ & =n^{-1}{\left[\left(\mathbf{I}-\left(\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\right)\mathbf{y}\right]}^{\mathbf{T}}\mathbf{W}\left[\left(\mathbf{I}-\left(\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\right)\mathbf{y}\right]\\ & =n^{-1}{‖{\mathbf{W}}^{\frac{1}{2}}\left(\mathbf{I}-\left(\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\right)\mathbf{y}‖}^2\end{array}$

where λ represents the smoothing parameter of the smoothing spline estimator, $\mathbf{K}$ represents the oscillation parameter of the Fourier series estimator,

• $\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)=\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\left[\mathbf{I}-\mathbf{T}{\left(\left(\mathbf{2}-{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\right)\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}-\mathbf{T}\right)}^{-\mathbf{1}}\left({\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}-\mathbf{I}\right)\left(\mathbf{I}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\right]$,

• $\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)=\mathbf{T}{\left[\left(\mathbf{2}-{\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}\right)\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\mathbf{T}-\mathbf{T}\right]}^{-\mathbf{1}}\left({\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)}^{\mathbf{T}}-\mathbf{I}\right)\left(\mathbf{I}-\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)\right)$, and

• $\mathbf{A}\left(\mathsf{\lambda },\mathbf{K}\right)=\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}+\mathbf{V}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\left(\mathbf{I}-\mathbf{U}{\left({\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\mathbf{U}\right)}^{-\mathbf{1}}{\mathbf{U}}^{\mathbf{T}}{\mathbf{H}}^{-\mathbf{1}}\mathbf{W}\right)$.

Secondly, based on the $MSE\left(\mathsf{\lambda },\mathrm{K}\right)$ given in Equation (49), we have the GCV (generalized cross-validation) for the MMNR model as follows:

(50)$\begin{array}{ll}GCV\left(\mathsf{\lambda },\mathrm{K}\right)& =\frac{MSE\left(\mathsf{\lambda },\mathrm{K}\right)}{{\left[n^{-1}trace\left(\mathbf{I}-\left(\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\right)\right]}^2}\\ & =\frac{n^{-1}{‖{\mathbf{W}}^{\frac{1}{2}}\left(\mathbf{I}-\left(\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\right)\mathbf{y}‖}^2}{{\left[n^{-1}trace\left(\mathbf{I}-\left(\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\right)\right]}^2}.\end{array}$

Hence, based on the GCV function presented in (50), we can obtain the optimal smoothing parameter ($\mathsf{\lambda }$) and optimal oscillation parameter (K) values, namely, (${\mathit{\lambda }}_{Opt},{\mathbf{K}}_{opt}$), by taking the solution to the following optimization:

(51)$\begin{array}{ll}{G}_{opt}\left({\mathit{\lambda }}_{Opt},{\mathbf{K}}_{opt}\right)& =\underset{\lambda \in {\mathbb{R}}^{+};K\in {\mathbb{R}}^{+}}{\mathrm{Min}}\left\{GCV\left(\mathsf{\lambda },\mathrm{K}\right)\right\}\\ & =\underset{\lambda \in {\mathbb{R}}^{+};K\in {\mathbb{R}}^{+}}{\mathrm{Min}}\left\{\frac{n^{-1}{‖{\mathbf{W}}^{\frac{1}{2}}\left(\mathbf{I}-\left(\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\right)\mathbf{y}‖}^2}{{\left[n^{-1}trace\left(\mathbf{I}-\left(\mathbf{F}\left(\mathsf{\lambda },\mathbf{K}\right)+\mathbf{G}\left(\mathsf{\lambda },\mathbf{K}\right)\right)\right)\right]}^2}\right\}\end{array}$

where ${\mathbb{R}}^{+}$ represents a positive real number set. Thus, the optimal values of both smoothing parameter and oscillation parameter, namely, (${\mathit{\lambda }}_{Opt},{\mathbf{K}}_{opt}$), can be obtained by taking the solution to the optimization function given in Equation (51).

Furthermore, in the following section, we determine one of the statistically sound estimator properties, namely, the consistency property.

3.7. Consistency of Regression Function Estimator of MMNR Model

To investigate the consistency of the regression function estimator of the MMNR model, ${\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})$, firstly, we should investigate asymptotic properties of the mixed smoothing spline and Fourier series estimator, ${\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})$, based on an integrated mean squared error (IMSE) criterion. In this study, we develop the IMSE from the uniresponse with one predictor case proposed by Wand and Jones [48] and multiresponse with one predictor case proposed by Lestari et al. [35] into a multiresponse with multipredictor case. Next, by considering Equation (48), the IMSE can be decomposed into two components, namely, the ${Biased}^2\left(\mathit{\delta }\right)$ component and $Var\left(\mathit{\delta }\right)$ component, as follows:

(52)$\begin{array}{ll}IMSE\left(\mathit{\delta }\right)& =E{\int }_a^b\left[{\left({\mathit{m}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})-{\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})\right)}^{\mathbf{T}}\mathbf{W}\left({\mathit{m}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})-{\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})\right)\right]d\mathbf{t}\\ & ={Biased}^2\left(\mathit{\delta }\right)+Var\left(\mathit{\delta }\right)\end{array}$

where ${Biased}^2\left(\mathit{\delta }\right)={\int }_a^{b}E[{\left({\mathit{m}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})-E{\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})\right)}^{\mathbf{T}}\mathbf{W}\left({\mathit{m}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})-E{\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})\right)]d\mathbf{t},$ $Var\left(\mathit{\delta }\right)={\int }_a^{b}E\left[{\left(E{\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})-{\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})\right)}^{\mathbf{T}}\mathbf{W}\left(E{\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})-{\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})\right)\right]d\mathbf{t}$, $\mathit{\delta }=(\mathsf{\lambda },\mathbf{K})$, $\mathbf{W}$ is a weight matrix, and ${\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})$ as given in Equation (48) is a regression function of the MMNR model.

Next, by using Theorem 2 in Lestari et al. [35], we obtain ${Biased}^2\left(\mathit{\delta }\right)\le O\left(\mathit{\delta }\right)$ as $n\to \infty$, where $O(·)$ represents the “big oh”. We can find the details of the “big oh” in Wand and Jones [48], and Sen and Singer [49]. Also, by using Theorem 3 in Lestari et al. [35], we obtain $Var({\mathit{\delta }}_r)\le O\left(\frac{1}{n{\left({\mathit{\delta }}_r\right)}^{1/2m}}\right)$ as $n\to \infty$ and for every response $r=\mathrm{1,2},\dots ,R$. Hence, we obtain the asymptotic property of the mixed regression function estimator of the MMNR model based on the IMSE criterion as follows:

(53)$IMSE\left(\mathit{\delta }\right)={Bias}^2\left(\mathit{\delta }\right)+Var\left(\mathit{\delta }\right)\le O\left(\mathit{\delta }\right)+O\left(\mathsf{\theta }\right) \mathrm{as} n\to \infty$

where $\mathit{\delta }={\left({\mathit{\delta }}_{\mathbf{1}},{\mathit{\delta }}_{\mathbf{2}},\dots ,{\mathit{\delta }}_{\mathit{R}}\right)}^{\mathrm{T}}$, $\mathsf{\theta }={\left(\begin{array}{ccc}\frac{\mathbf{1}}{\mathit{n}{\mathit{\delta }}_{\mathbf{1}}^{\frac{\mathbf{1}}{\mathbf{2}\mathit{m}}}}& \frac{\mathbf{1}}{\mathit{n}{\mathit{\delta }}_{\mathbf{2}}^{\frac{\mathbf{1}}{\mathbf{2}\mathit{m}}}}& \begin{array}{cc}\cdots & \frac{\mathbf{1}}{\mathit{n}{\mathit{\delta }}_{\mathit{R}}^{\frac{\mathbf{1}}{\mathbf{2}\mathit{m}}}}\end{array}\end{array}\right)}^{\mathit{T}}$, ${\mathit{\delta }}_{\mathit{r}}=({\mathsf{\lambda }}_{\mathit{r}},{\mathbf{K}}_{\mathit{r}})$, and $r=\mathrm{1,2},\dots ,R$.

Furthermore, based on Equation (53), we obtain an inequality as follows:

$IMSE\left(\mathit{\delta }\right)={Bias}^2\left(\mathit{\delta }\right)+Var\left(\mathit{\delta }\right)\le O\left(\mathit{\delta }\right)+O\left(\mathsf{\theta }\right)\le O\left(n\mathit{\delta }\right) \mathrm{as} n\to \infty .$

Hence, according to Sen and Singer [49] and Serfling [50], for any small positive number, $\vartheta >0$, we obtain the following relationship:

(54)$\underset{n\to \infty }{\mathrm{Lim}}P\left(\left|\frac{IMSE\left(\mathit{\delta }\right)}{n\mathit{\delta }}\right|>\vartheta \right)\le \underset{n\to \infty }{\mathrm{Lim}}P\left(\left|IMSE\left(\mathit{\delta }\right)\right|>\vartheta \right)=0.$

Since $P\left(\left|IMSE\left(\mathit{\delta }\right)\right|>\vartheta \right)=1-P(\left|IMSE\left(\mathit{\delta }\right)\right|\le \vartheta )$ then by considering Equation (54) and by applying the properties of probability, we obtain the following equation:

(55)$\begin{array}{r}\underset{n\to \infty }{\mathrm{Lim}}\left(1-P(\left|IMSE\left(\mathit{\delta }\right)\right|\le \vartheta \right))=0 \iff 1-\underset{n\to \infty }{\mathrm{Lim}}P\left(\left|IMSE\left(\mathit{\delta }\right)\right|\le \vartheta \right)=0\\ \iff \underset{n\to \infty }{\mathrm{Lim}}P\left(\left|IMSE\left(\mathit{\delta }\right)\right|\le \vartheta \right)=1.\end{array}$

Equation (55) shows that the regression function estimator, ${\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})$, of the multiresponse multipredictor nonparametric regression (MMNR) model obtained by using the mixed smoothing spline and Fourier series estimator is a consistent estimator based on the integrated mean squared error (IMSE) criterion.

Next, in the following section, we provide a simulation study.

3.8. Simulation Study

In this section, we provide a simulation study to show the performance of the mixed smoothing spline and Fourier series estimator in estimating the MMNR model. In other words, the simulation study provided is intended to show the sensitivity of parameter selection and its impact on model performance. In this simulation study, we used generation data of $n=50$ with a correlation of 0.95 and a variance of 0.25. Here, we use a three-response MMNR model as follows:

(56)$\left.\begin{array}{c}{y}_{1i}=17{\left(x_{1i}-0.8\right)}^3-1.5t_{1i}-\cos \left(6\pi t_{1i}\right)+{\epsilon }_{1i} , i=\mathrm{1,2},\dots ,50.\\ y_{2i}=18{\left(x_{2i}-0.8\right)}^3-1.75t_{2i}-1.2\cos \left(6\pi t_{2i}\right)+{\epsilon }_{2i} , i=\mathrm{1,2},\dots ,50.\\ y_{3i}=19{\left(x_{3i}-0.8\right)}^3-2t_{3i}-1.5\cos \left(6\pi t_{3i}\right)+{\epsilon }_{3i} , i=\mathrm{1,2},\dots ,50.\end{array}\right\}$

Next, to show the influence of smoothing parameters (λ) and oscillation parameters (K), we use different smoothing parameters by generating lambda (λ) randomly, and different oscillation parameters by moving K from 1 to 4. We carry out this process 15 times until the optimal values of smoothing parameters (λ) and oscillation parameters (K) are obtained.

Simulation results including MSE (mean squared error) values, minimum GCV (generalized cross-validation) values, coefficient of determination ($R^2$) values, and lambda (λ) values for every oscillation parameter value $(\mathbf{K}=\mathrm{1,2},\mathrm{3,4})$ are given in Table 1.

Table 1 presents the optimal lambda values based on the minimum GCV values from 15 repetitions of the MMNR model presented in Equation (56) for every oscillation parameter (K) where $\mathbf{K}=\mathrm{1,2},\mathrm{3,4}$. From Table 1, we can observe that every different value of both the smoothing parameter and oscillation parameter will give different performance to the estimation performance. These are indicated by the different coefficient of determination ($R^2$) and MSE (mean squared errors) values as the values change due to oscillation parameters (K) and smoothing parameters (λ).

Next, based on all the minimum GCV values presented in Table 1, we chose the minimum GCV value from the four oscillation parameter (K) values. In Table 1, the minimum GCV value is produced by $\mathbf{K}=1$, namely, 0.5786323, with an MSE value of 1.02363794, $R^2$ value of 0.95311712, and ${\lambda }_1=0.84886039$; ${\lambda }_2=0.0021173$; ${\lambda }_3=0.001$. In Table 1, these values are shown in bold. These values are the optimal values and the best MMNR estimate is obtained based on these optimal values. Therefore, these values are then used to determine the estimated values of the MMNR model and to produce the plots of the estimated values of the MMNR model presented in Equation (56).

Hereinafter, we obtained the plots of estimated model values versus predictor (observation) values for every response together with the plots of their actual values. The plots of estimated model values versus predictor (observation) values for the first response together with the plots of actual values of the first response are presented in Figure 1. The plots of estimated model values versus predictor (observation) values for the second response together with the plots of actual values of the second response are presented in Figure 2. The plots of estimated model values versus predictor (observation) values for the third response together with the plots of actual values of the third response are presented in Figure 3.

From Figure 1, Figure 2 and Figure 3, we can observe that the values of the estimated MMNR models using the mixed smoothing spline and Fourier series estimator are close to the actual values. This is shown by the performance curves of the values of the estimated MMNR models, which tend to be similar to the actual values curves. Hence, based on the results of the simulation study presented in Table 1, and Figure 1, Figure 2 and Figure 3, it shows that the estimation value of the MMNR model in (56) depends on both the smoothing spline parameter and oscillation parameter. In other words, it shows that the selection of the smoothing spline parameter and oscillation parameter is sensitive and impacts the MMNR model’s performance.

Hereinafter, based on the results of the simulation study, we obtained the surface plots of the estimation results of the MMNR model presented in Equation (56). These surface plots of the estimated MMNR model are presented in three-dimensional plots (see Figure 4, Figure 5 and Figure 6).

4. Conclusions

The estimation result of the MMNR model by using the mixed smoothing spline and Fourier series estimator, presented by ${\widehat{\mathit{y}}=\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})$, is a combination of the smoothing spline component and Fourier series component estimators presented by ${\widehat{\mathbf{f}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})$ and ${\widehat{\mathbf{g}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})$, respectively. The estimator of the MMNR model is highly dependent on the smoothing parameter ($\lambda$) and oscillation parameter (K). These optimal parameter values, namely, (${\mathit{\lambda }}_{Opt},{\mathbf{K}}_{opt}$), are determined by using the GCV criterion to determine the values that minimize the function $GCV\left(\mathsf{\lambda },\mathrm{K}\right)$. Also, the estimator of the MMNR model by using the mixed smoothing spline and Fourier series estimator that we have obtained is linear with respect to the observations given in Equation (48). Additionally, the regression function estimator, ${\widehat{\mathit{m}}}_{\left(\mathsf{\lambda },\mathbf{K}\right)}(\mathbf{x},\mathbf{t})$, of the MMNR model obtained by using the mixed smoothing spline and Fourier series estimator is a consistent estimator based on the IMSE criterion. This means that the mixed smoothing spline and Fourier series estimator is statistically a good estimator because it meets one of the criteria for the goodness of an estimator, namely, being consistent. It is therefore suitable for analyzing data with patterns that partly change at certain subintervals, and others that follow a recurring pattern in a certain trend. In addition, the estimated weight matrix that we obtained is a diagonal matrix where the main diagonal elements are the estimated weight matrix of the first response, the estimated weight matrix of second response, and so on, until the estimated weight matrix of the R-th response. The results of this study contribute to developing statistical inference theories such as estimation and hypothesis testing theories, especially the development of statistical inference theory of nonparametric regression.

Author Contributions

All authors contributed to this research article, namely, conceptualization, N.C., B.L. and I.N.B.; methodology, N.C., B.L. and I.N.B.; software, B.L., N.C. and D.A.; validation, B.L., N.C., I.N.B. and D.A.; formal analysis, B.L., N.C., I.N.B. and D.A.; investigation, resource and data curation, B.L., N.C., D.A. and I.N.B.; writing—original draft preparation, B.L. and N.C.; writing—review and editing, B.L. and N.C.; visualization, B.L. and N.C.; supervision, B.L., N.C., I.N.B. and D.A.; project administration, N.C. and B.L.; funding acquisition, N.C. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

All data are contained within this article.

Acknowledgments

The authors thank Airlangga University for technical support. The authors are also grateful to the editors and anonymous peer-reviewers of the Symmetry journal for providing comments, corrections, criticisms, and suggestions that were useful for improving the quality of this article.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Plots of estimated model values versus predictor (observation) values for the first response ([Forumla omitted. See PDF.]) together with the plots of actual values of the first response.

Figure 2. Plots of estimated model values versus predictor (observation) values for the second response ([Forumla omitted. See PDF.]) together with the plots of actual values of the second response.

Figure 3. Plots of estimated model values versus predictor (observation) values for the third response ([Forumla omitted. See PDF.]) together with the plots of actual values of the third response.

Figure 4. Surface plot of the estimated model versus predictors for the first response ([Forumla omitted. See PDF.]).

Figure 5. Surface plot of the estimated model versus predictors for the second response ([Forumla omitted. See PDF.]).

Figure 6. Surface plot of the estimated model versus predictors for the third response ([Forumla omitted. See PDF.]).

Table 1

Values of MSE, minimum GCV, $R^2$, and λ.

K	MSE	Minimum GCV	${\mathit{R}}^2$	λ
1	1.02363794	0.5786323	0.95311712	${\mathit{\lambda }}_{\mathbf{1}}=\mathbf{0.84886039}\mathbf{;}$ ${\mathit{\lambda }}_{\mathbf{2}}=\mathbf{0.00211737}$;${\mathit{\lambda }}_{\mathbf{3}}=\mathbf{0.001}$.
2	2.20132482	2.0945904	0.90018631	${\lambda }_1=0.55180324;$ ${\lambda }_2=0.89946915$;${\lambda }_3=0.54313339$.
3	2.09512311	2.17049788	0.90527677	${\lambda }_1=0.74451376;$ ${\lambda }_2=0.62104441$;${\lambda }_3=0.27310874$.
4	2.03215324	2.10132858	0.90769321	${\lambda }_1=0.56110414;$ ${\lambda }_2=0.26679346$; ${\lambda }_3=0.4535796$.

References

1. Eubank, R.L. Nonparametric Regression and Spline Smoothing; 2nd ed. Marcel Dekker: New York, NY, USA, 1999.

2. Cheruiyot, L.R. Local linear regression estimator on the boundary correction in nonparametric regression estimation. J. Stat. Theory Appl.; 2020; 19, pp. 460-471. [DOI: https://dx.doi.org/10.2991/jsta.d.201016.001]

3. Cheng, M.-Y.; Huang, T.; Liu, P.; Peng, H. Bias reduction for nonparametric and semiparametric regression models. Stat. Sin.; 2018; 28, pp. 2749-2770. [DOI: https://dx.doi.org/10.5705/ss.202017.0058]

4. Chamidah, N.; Zaman, B.; Muniroh, L.; Lestari, B. Designing local standard growth charts of children in East Java province using a local linear estimator. Int. J. Innov. Creat. Change; 2020; 13, pp. 45-67.

5. Delaigle, A.; Fan, J.; Carroll, R.J. A design-adaptive local polynomial estimator for the errors-in-variables problem. J. Am. Stat. Assoc.; 2009; 104, pp. 348-359. [DOI: https://dx.doi.org/10.1198/jasa.2009.0114]

6. Francisco-Fernandez, M.; Vilar-Fernandez, J.M. Local polynomial regression estimation with correlated errors. Comm. Stat. Theory Methods; 2001; 30, pp. 1271-1293. [DOI: https://dx.doi.org/10.1081/STA-100104745]

7. Benhenni, K.; Degras, D. Local polynomial estimation of the mean function and its derivatives based on functional data and regular designs. ESAIM Probab. Stat.; 2014; 18, pp. 881-899. [DOI: https://dx.doi.org/10.1051/ps/2014009]

8. Kikechi, C.B. On local polynomial regression estimators in finite populations. Int. J. Stats. Appl. Math.; 2020; 5, pp. 58-63.

9. Wand, M.P.; Jones, M.C. Kernel Smoothing; 1st ed. Chapman and Hall/CRC: New York, NY, USA, 1995.

10. Cui, W.; Wei, M. Strong consistency of kernel regression estimate. Open J. Stats.; 2013; 3, pp. 179-182. [DOI: https://dx.doi.org/10.4236/ojs.2013.33020]

11. De Brabanter, K.; De Brabanter, J.; Suykens, J.A.K.; De Moor, B. Kernel regression in the presence of correlated errors. J. Mach. Learn. Res.; 2011; 12, pp. 1955-1976.

12. Wahba, G. Spline Models for Observational Data; SIAM: Philadelphia, PA, USA, 1990.

13. Wang, Y. Smoothing Splines: Methods and Applications; Taylor & Francis Group: Boca Raton, FL, USA, 2011.

14. Liu, A.; Qin, L.; Staudenmayer, J. M-type smoothing spline ANOVA for correlated data. J. Multivar. Anal.; 2010; 101, pp. 2282-2296. [DOI: https://dx.doi.org/10.1016/j.jmva.2010.06.001]

15. Gao, J.; Shi, P. M-Type smoothing splines in nonparametric and semiparametric regression models. Stat. Sin.; 1997; 7, pp. 1155-1169.

16. Chamidah, N.; Lestari, B.; Massaid, A.; Saifudin, T. Estimating mean arterial pressure affected by stress scores using spline nonparametric regression model approach. Commun. Math. Biol. Neurosci.; 2020; 2020, 72.

17. Chamidah, N.; Lestari, B.; Budiantara, I.N.; Saifudin, T.; Rulaningtyas, R.; Aryati, A.; Wardani, P.; Aydin, D. Consistency and asymptotic normality of estimator for parameters in multiresponse multipredictor semiparametric regression model. Symmetry; 2022; 14, 336. [DOI: https://dx.doi.org/10.3390/sym14020336]

18. Lestari, B.; Chamidah, N.; Budiantara, I.N.; Aydin, D. Determining confidence interval and asymptotic distribution for parameters of multiresponse semiparametric regression model using smoothing spline estimator. J. King Saud Univ.-Sci.; 2023; 35, 102664. [DOI: https://dx.doi.org/10.1016/j.jksus.2023.102664]

19. Tirosh, S.; De Ville, D.V.; Unser, M. Polyharmonic smoothing splines and the multidimensional Wiener filtering of fractal-like signals. IEEE Trans. Image Process.; 2006; 15, pp. 2616-2630. [DOI: https://dx.doi.org/10.1109/TIP.2006.877390]

20. Irizarry, R.A. Choosing Smoothness Parameters for Smoothing Splines by Minimizing an Estimate of Risk. Available online: https://www.biostat.jhsph.edu/~ririzarr/papers/react-splines.pdf (accessed on 3 February 2024).

21. Adams, S.O.; Ipinyomi, R.A.; Yahaya, H.U. Smoothing spline of ARMA observations in the presence of autocorrelation error. Eur. J. Stats. Prob.; 2017; 5, pp. 1-8.

22. Adams, S.O.; Yahaya, H.U.; Nasiru, O.M. Smoothing parameter estimation of the generalized cross-validation and generalized maximum likelihood. IOSR J. Math.; 2017; 13, pp. 41-44.

23. Lee, T.C.M. Smoothing parameter selection for smoothing splines: A simulation study. Comput. Stats. Data Anal.; 2003; 42, pp. 139-148. [DOI: https://dx.doi.org/10.1016/S0167-9473(02)00159-7]

24. Maharani, M.; Saputro, D.R.S. Generalized cross-validation (GCV) in smoothing spline nonparametric regression models. J. Phys. Conf. Ser.; 2021; 1808, 12053. [DOI: https://dx.doi.org/10.1088/1742-6596/1808/1/012053]

25. Wang, Y.; Ke, C. Smoothing spline semiparametric nonlinear regression models. J. Comp. Graph. Stats.; 2009; 18, pp. 165-183. [DOI: https://dx.doi.org/10.1198/jcgs.2009.0010]

26. Gu, C. Smoothing Spline ANOVA Models; Springer: New York, NY, USA, 2002.

27. Sun, X.; Zhong, W.; Ma, P. An asymptotic and empirical smoothing parameters selection method for smoothing spline ANOVA models in large samples. Biometrika; 2021; 108, pp. 149-166. [DOI: https://dx.doi.org/10.1093/biomet/asaa047] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34294943]

28. Wang, Y.; Guo, W.; Brown, M.B. Spline smoothing for bivariate data with applications to association between hormones. Stat. Sin.; 2000; 10, pp. 377-397.

29. Lu, M.; Liu, Y.; Li, C.-S. Efficient estimation of a linear transformation model for current status data via penalized splines. Stat. Methods Med. Res.; 2020; 29, pp. 3-14. [DOI: https://dx.doi.org/10.1177/0962280218820406] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30592240]

30. Berry, L.N.; Helwig, N.E. Cross-validation, information theory, or maximum likelihood? A comparison of tuning methods for penalized splines. Stats; 2021; 4, pp. 701-724. [DOI: https://dx.doi.org/10.3390/stats4030042]

31. Islamiyati, A.; Zakir, M.; Sirajang, N.; Sari, U.; Affan, F.; Usrah, M.J. The use of penalized weighted least square to overcome correlations between two responses. BAREKENG J. Ilmu Mat. Dan Terap.; 2022; 16, pp. 1497-1504. [DOI: https://dx.doi.org/10.30598/barekengvol16iss4pp1497-1504]

32. Islamiyati, A.; Raupong,; Kalondeng, A.; Sari, U. Estimating the confidence interval of the regression coefficient of the blood sugar model through a multivariable linear spline with known variance. Stat. Transit. New Ser.; 2022; 23, pp. 201-212. [DOI: https://dx.doi.org/10.2478/stattrans-2022-0012]

33. Kirkby, J.L.; Leitao, A.; Nguyen, D. Nonparametric density estimation and bandwidth selection with B-spline basis: A novel Galerkin method. Comput. Stats. Data Anal.; 2021; 159, 107202. [DOI: https://dx.doi.org/10.1016/j.csda.2021.107202]

34. Osmani, F.; Hajizadeh, E.; Mansouri, P. Kernel and regression spline smoothing techniques to estimate coefficient in rates model and its application in psoriasis. Med. J. Islam. Repub. Iran; 2019; 33, 90. [DOI: https://dx.doi.org/10.47176/mjiri.33.90] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31696084]

35. Lestari, B.; Chamidah, N.; Aydin, D.; Yilmaz, E. Reproducing kernel Hilbert space approach to multiresponse smoothing spline regression function. Symmetry; 2022; 14, 2227. [DOI: https://dx.doi.org/10.3390/sym14112227]

36. Bilodeau, M. Fourier smoother and additive models. Can. J. Stat.; 1992; 20, pp. 257-269. [DOI: https://dx.doi.org/10.2307/3315313]

37. Suparti, S.; Prahutama, A.; Santoso, R.; Devi, A.R. Spline-Fourier’s Method for Modelling Inflation in Indonesia. E3S Web Conf.; 2018; 73, 13003. [DOI: https://dx.doi.org/10.1051/e3sconf/20187313003]

38. Mardianto, M.F.F.; Gunardi,; Utami, H. An analysis about Fourier series estimator in nonparametric regression for longitudinal data. Math. Stats.; 2021; 9, pp. 501-510. [DOI: https://dx.doi.org/10.13189/ms.2021.090409]

39. Amato, U.; Antoniadis, A.; De Feis, I. Fourier series approximation of separable models. J. Comput. Appl. Math.; 2002; 146, pp. 459-479. [DOI: https://dx.doi.org/10.1016/S0377-0427(02)00398-9]

40. Mariati, M.P.A.M.; Budiantara, I.N.; Ratnasari, V. The application of mixed smoothing spline and Fourier series model in nonparametric regression. Symmetry; 2021; 13, 2094. [DOI: https://dx.doi.org/10.3390/sym13112094]

41. Aronszajn, N. Theory of reproducing kernels. Trans. Am. Math. Soc.; 1950; 68, pp. 337-404. [DOI: https://dx.doi.org/10.1090/S0002-9947-1950-0051437-7]

42. Kimeldorf, G.; Wahba, G. Some results on Tchebycheffian spline functions. J. Math. Anal. Appl.; 1971; 33, pp. 82-95. [DOI: https://dx.doi.org/10.1016/0022-247X(71)90184-3]

43. Berlinet, A.; Thomas-Agnan, C. Reproducing Kernel Hilbert Spaces in Probability and Statistics; Kluwer Academic: Norwell, MA, USA, 2004.

44. Paulsen, V.I. An Introduction to the Theory of Reproducing Kernel Hilbert Space. Research Report. 2009; Available online: https://www.researchgate.net/publication/255635687_AN_INTRODUCTION_TO_THE_THEORY_OF_REPRODUCING_KERNEL_HILBERT_SPACES (accessed on 24 March 2022).

45. Yuan, M.; Cai, T.T. A reproducing kernel Hilbert space approach to functional linear regression. Ann. Stat.; 2010; 38, pp. 3412-3444. [DOI: https://dx.doi.org/10.1214/09-AOS772]

46. Johnson, R.A.; Wichern, D.W. Applied Multivariate Statistical Analysis; Prentice Hall: New York, NY, USA, 1982.

47. Ruppert, D.; Carroll, R. Penalized Regression Splines, Working Paper; School of Operation Research and Industrial Engineering, Cornell University: Ithaca, NY, USA, 1997.

48. Wand, M.P.; Jones, M.C. Kernel Smoothing; Chapman & Hall: London, UK, 1995.

49. Sen, P.K.; Singer, J.M. Large Sample in Statistics: An Introduction with Applications; Chapman & Hall: London, UK, 1993.

50. Serfling, R.J. Approximation Theorems of Mathematical Statistics; John Wiley: New York, NY, USA, 1980.