# A New Mixed Estimator in Nonparametric Regression for Longitudinal Data

Authors: Made Ayu Dwi Octavanny; I Nyoman Budiantara; Heri Kuswanto; Dyah Putri Rahmawati
Journal: Journal of Mathematics
 (2021)

---

## Abstract

We introduce a new method for estimating the nonparametric regression curve for longitudinal data. This method combines two estimators: truncated spline and Fourier series. This estimation is completed by minimizing the penalized weighted least squares and weighted least squares. This paper also provides the properties of the new mixed estimator, which are biased and linear in the observations. The best model is selected using the smallest value of generalized cross-validation. The performance of the new method is demonstrated by a simulation study with a variety of time points. Then, the proposed approach is applied to a stroke patient dataset. The results show that simulated data and real data yield consistent findings.


---

## Body

## 1. Introduction

Nonparametric regression is a statistical method used if the data show an unknown regression curve. The strength of this method is its great ﬂexibility since the data are used to ﬁnd the form of its estimated regression curve without being inﬂuenced by subjective judgements [1]. Some estimators used are spline, Fourier series, kernel, and polynomial.The truncated spline is a function where there is a change in the behaviour of the curve in certain subintervals. The spline is one of the popular estimators in nonparametric regression because it has an excellent visual interpretation. Montoya et al. [2] conducted a simulation study to compare knot selection methods in a penalized regression spline model. Next, the Fourier series is much used to describe curves that present sine and cosine waves. This estimator is commonly used when the data have the characteristics of a periodicity. Bilodeau [3] estimated additive components with functions consisting of truncated Fourier cosine series, using penalized least squares to obtain the coefficients. Furthermore, Tripena [4] developed a Fourier series estimator for bi-response nonparametric regression. Gong and Gao [5] discussed a nonparametric kernel estimation of tax policy’s impact on the demand for private health insurance in Australia.A single estimator in nonparametric regression is commonly used but does not limit the possibility of developing into a mixed estimator. There are many cases where each predictor variable has a different pattern. Therefore, applying only a single estimator can make the regression model’s estimation incorrect and produce a large error. Some previous studies that have been mentioned have limitations in that they can only be used for cross-sectional data. To overcome this limitation, a longitudinal data model has been developed. The use of longitudinal data has increased in recent years because it can be applied in various fields. Longitudinal data are data obtained in observations of independent subjects, where each subject is observed repeatedly at a certain period. This has the advantage of being able to observe changes based on time [6].In [7], the mixed estimator is limited in the sample size it can treat. The present study extends the use of the mixed truncated spline and Fourier series (MTSFS) model to larger sample sizes and various time point designs. Some properties of the new mixed estimator will also be provided. We use generalized cross-validation (GCV) to determine the best model from various knots, oscillations, and smoothing parameters. A simulation study and real data are employed to demonstrate the performance of the proposed method. The case study includes the factors that affect the Glasgow Coma Scale (GCS) on stroke patients, i.e., body temperature (BT) and pulse rate (PR).The rest of this paper is divided into several main topics. Section2 introduces the materials and methods used in this study. Section 3 consists of five subsections: the theory of the new mixed estimator, its properties, the selection of the best model, the simulation study, and a case study. The details of the new mixed estimator are presented in Section 3.1, followed by its properties in Section 3.2. Section 3.3 presents how to select the optimum knot point, oscillation parameter, and smoothing parameter to obtain the best model. The results of a simulation study of the proposed method are presented in Section 3.4, followed by an application to real data in Section 3.5. Section 4 concludes the paper.


## 2. Materials and Methods

Supposeyit is the response variable and xit and zit are the corresponding predictor variables with sample size of n subjects i=1,2,…,n, with each subject having T observations t=1,2,…,T. The relations between the response and predictors for the nonparametric regression model for longitudinal data are(1)yit=μx1it,…,xpit,z1it,…,zqit+εit,where μ is the regression curve and εit is a random error. Assume that the form of the regression curve μ is unknown and additive, so that(2)μx1it,…,xpit,z1it,…,zqit=∑j=1pfjixjit+∑k=1qgkizkit,where ∑j=1pfjixjit is the truncated spline component and ∑k=1qgkizkit is the Fourier series component.The functionfji,j=1,2,…p, is an approximation using truncated spline functions and gki,k=1,2,…,q, is that using Fourier series. The estimator μ is obtained through a two-step optimization, i.e., penalized weighted least squares (PWLS) and weighted least squares (WLS).


## 3. Results and Discussion

### 3.1. MTSFS Model for Longitudinal Data

Some lemmas and theorems are provided to obtain an MTSFS model for longitudinal data. Lemma1 presents the goodness of fit of the Fourier series component, Lemma 2 presents the penalty component, and Lemma 3 gives the solution for the truncated spline component. Theorem 1 presents the PWLS optimization, and Theorem 2 presents WLS optimization. The results are summarized as follows.Lemma 1.
If the Fourier series component in equation (2) is given by ∑k=1qgkizkit, then the goodness of fit can be formulated as follows:(3)N−1y∗−Zc′Wy∗−Zc,where N=n×T, W is the nT×nT weighting matrix, y∗=y−f, and(4)y∗=y1∗′y2∗′…yn∗′′,yi∗=yi1∗yi2∗…yiT∗′,f=f1′f2′…fn′′,fi=∑j=1pfjixji1,∑j=1pfjixji2,…,∑j=1pfjixjiT′,Z=Z10…00Z2…0⋮⋮⋱⋮00…Zn,c=c1c2⋮cn,Zi=z1i112cosz1i1cos2z1i1…cosHz1i1…zqi112coszqi1cos2zqi1…cosHzqi1z1i212cosz1i2cos2z1i2…cosHz1i2…zqi212coszqi2cos2zqi2…cosHzqi2⋮⋮⋮⋮⋰⋮⋰⋮⋮⋮⋮⋰⋮z1iT12cosz1iTcos2z1iT…cosHz1iT…zqiT12coszqiTcos2zqiT…cosHzqiT,ci=d1ic01ic11ic21i…cH1i…dqic0qic1qic2qi…cHqi′.Proof.
The functiongki is assumed to be unknown and contained in the space C0,π. The function is approximated using Fourier series with a trend, modified from Bilodeau [3]:(5)gkizkit=dkizkit+12c0ki+∑h=1Hchkicoshzkit.
Equation (5) can be written in matrix form(6)g=g1′g2′⋮gn′=Z1c1Z2c2⋮Zncn=Z10…00Z2…0⋮⋮⋱⋮00…Znc1c2⋮cn=Zc,where Z is a nT×2+Hnq matrix and c is a 2+Hnq vector.
The nonparametric model in equation (2) can be rewritten as(7)yit−∑j=1pfjixjit=∑k=1qgkizkit+εit,yit∗=∑k=1qgkizkit+εit,i=1,2,…,n,t=1,2,…,T.
Then, the goodness of fit for equation (7) can be written as(8)N−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tgkizkit2=N−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tdkizkit+12c0ki+∑h=1Hchkicoshzkit2.
The model in equation (8) can be written in matrix form:(9)N−1Wy∗−Zc2=N−1y∗−Zc′Wy∗−Zc.Lemma 2.
If the penalty component is given,(10)∑k=1qλk∫0π2πgki″zk2dzk,then(11)∑k=1qλk∫0π2πgki″zk2dzk=c′Dλc,where(12)c=c1c2…cn′,ci=d1ic01ic11ic21i…cH1i…dqic0qic1qic2qi…cHqi′,Dλ=D1λ0…00D2λ…0⋮⋮⋱⋮00…Dnλ,Diλ=d10…00d2…0⋮⋮⋱⋮00…dq,Ik=λk14λk24…λkH4,dk=diag0,0,Ik,k=1,2,…,q.Proof.
Regarding equation (5), we define(13)gki″zk=ddzkddzkdkizkit+12c0ki+∑h=1Hchkicoshzkit=−∑h=1Hh2chkicoshzkit.
Consequently,(14)Pkgk=∫0π2π∑h=1Hh2chkicoshzkit2dzk=2π∫0π∑h=1Hh2chkicoshzkit2+2∑h<mHh2chkicoshzkitm2cmkicosmzkitdzk.
Let(15)A=2π∑h=1H∫0πh2chkicoshzkit2dzk,B=2π2∑h<lH∫0πh2chkicoshzkitl2clkicoslzkitdzk.
Thus,(16)Pkgk=A+B=∑h=1Hh4chki2.
The penalty in equation (16) can be rewritten in matrix form(17)c′Dλc.
Regarding the goodness of fit in Lemma1 and the penalty component in Lemma 2, we obtain the PWLS optimization as(18)N−1∑i=1n∑t=1Twityit∗−∑k=1qgkizkit2+∑k=1qλk∫0π2πgki″zk2dzk,0<λk<∞.
Equation (18) can be rewritten in the form(19)N−1y∗−Zc′Wy∗−Zc+c′Dλc.Theorem 1.
If the goodness of fit is as in Lemma1 and the penalty component is as in Lemma 2, then the Fourier series component obtained by minimizing the PWLS in equation (18) is(20)g^k,h,λx,z=Ly∗,where y∗=y−f and L=ZZ′WZ+NDλ−1Z′W.Proof.
The optimization in equation (18) can be written as(21)Mingk∈C0,πN−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tgkizkit2+∑k=1qλk∫0π2πgk″zk2dzk.
Equation (21) can be rewritten in the form(22)Minc∈C0,πN−1y∗−Zc′Wy∗−Zc+c′Dλc=Minc∈C0,πQc.
We obtain(23)Qc=N−1y∗′Wy∗−2N−1c′Z′Wy∗+N−1c′Z′WZc+c′Dλc.
The completion of the optimization in (23) is obtained by taking the partial derivative of Qc with respect to c and setting it equal to zero, i.e.,(24)∂Qc∂c=0,giving the result(25)c^=Z′WZ+NDλ−1Z′Wy∗.
By substituting (25) into (6), we obtain(26)g^k,h,λx,z=Zc^=ZZ′WZ+NDλ−1Z′Wy∗=Ly∗.
The nonparametric regression model in equation (7) can be written as(27)y∗=g^k,h,λx,z=Ly∗,where y∗=y−f and L=ZZ′WZ+NDλ−1Z′W.Lemma 3.
If the truncated spline component in equation (2) is ∑j=1pfjixjit, then the WLS is(28)I−Ly−I−LMγ′WI−Ly−I−LMγ,where L=ZZ′WZ+NDλ−1Z′W.Proof.
The functionfji is a linear truncated spline function with s knot for each xj,j=1,2,…,p:(29)fjixjit=αjixjit+∑u=1sβujixjit−Kuji+,where xjit−Kuji+=xjit−Kuji,xjit≥Kuji,0,xjit<Kuji.
Equation (29) can be rewritten in matrix form(30)f=X|Sα¯β=Mγ,where(31)X=X10⋯00X2⋯0⋮⋮⋱⋮00⋯Xn,S=S10⋯00S2⋯0⋮⋮⋱⋮00⋯Sn,α=α1α2⋮αn,β=β1β2⋮βn,Xi=x1i1x2i1…xpi1x1i2x2i2…xpi2⋮⋮⋱⋮x1iTx2iT…xpiT,αi=α1iα2i⋮αpi,Si=x1i1−K11i+1…x1i1−Ks1i+1…xpi1−K1pi+1…xpi1−Kspi+1x1i2−K11i+1…x1i2−Ks1i+1…xpi2−K1pi+1…xpi2−Kspi+1⋮⋱⋮⋱⋮⋱⋮x1iT−K11i+1…x1iT−Ks1i+1…xpiT−K1pi+1…xpiT−Kspi+1,βi=β11i…βs1i…β1pi…βspiT.
The MTSFS model for longitudinal data in equation (1) can be written in the form(32)y=f+g+ε.
By substituting (26) into (32), we obtain(33)y=f+Ly∗+ε.
To obtain the estimator of truncated spline component, equation (33) can be written as(34)y−Ly∗=f+ε,y−Ly−f=f+ε,I−Ly=I−Lf+ε.
Equation (34) can be rewritten as(35)I−Ly=I−LMγ+ε.
Thus,(36)ε=I−Ly−I−LMγ.
As a consequence, the WLS is given by(37)ε′ε=I−Ly−I−LMγ′WI−Ly−I−LMγ.Theorem 2.
Suppose the WLS is as in Lemma3. Then, the mixed estimator obtained by minimizing WLS is(38)f^k,h,λx,z=MJ−1Ky,where J=2−L′WLM−WM and K=L′−IWI−L.Proof.
The WLS optimization in Lemma3 can be written as(39)MinγI−Ly−I−LMγ′WI−Ly−I−LMγ=MinγQγ.
We obtain(40)Qγ=y′Wy+y′L′WLy+γ′M′L′WLMγ+γ′M′WMγ−2y′L′Wy+2γ′M′L′Wy+−2N−1γ′M′Wy−2N−1γ′M′L′WLy+2N−1γ′M′WLy−2N−1γ′M′WLMγ.
The complete optimization is obtained by setting equal to zero the partial derivative ofQγ with respect to γ, that is,(41)∂Qγ∂γ=0,giving the result(42)γ^=2−L′WLM−WM−1L′−IWI−Ly=J−1Ky.
By substitutingγ^ into equation (30), we obtain(43)f^k,h,λx,z=Mγ^=MJ−1Ky=Ak,h,λy.
We obtainc^ by substituting (43) into (25):(44)c^=Z′WZ+NDλ−1Z′Wy∗=Z′WZ+NDλ−1Z′Wy−f^=Z′WZ+NDλ−1Z′Wy−Mγ^=Z′WZ+NDλ−1Z′Wy−MJ−1Ky=Z′WZ+NDλ−1Z′WI−MJ−1Ky.
As a consequence,(45)g^k,h,λx,z=Zc^=ZZ′WZ+NDλ−1Z′WI−MJ−1Ky=LI−MJ−1Ky=Bk,h,λy.
By substitutingγ^ and c^ into the MTSFS model for longitudinal data, we obtain(46)y^=μ^k,h,λx,z=f^k,h,λx,z+g^k,h,λx,z=Mγ^+Zc^=MJ−1Ky+LI−MJ−1Ky=MJ−1K+LI−MJ−1Ky=Ck,h,λy,where J=2−L′WLM−WM, K=L′−IWI−L, and L=ZZ′WZ+NDλ−1Z′W.


### 3.2. The Properties of MTSFS Model for Longitudinal Data

This section provides the MTSFS model’s properties for longitudinal data, i.e., it is biased and linear in observations. The mixed estimator is biased, as proved by(47)Eμ^k,h,λx,z=Ef^k,h,λx,z+g^k,h,λx,z=EAk,h,λy+Bk,h,λy=ECk,h,λy=Ck,h,λEy.Assuming thatCk,h,λ≠I,(48)Eμ^k,h,λx,z=Ck,h,λEy=Ck,h,λfk,h,λx,z+gk,h,λx,z=Ck,h,λμk,h,λx,z.The result in equation (48) showed that the mixed estimator is biased because Eμ^k,h,λx,z≠μk,h,λx,z. Even though the mixed estimator is biased, the mixed estimator is linear in observations proved by equation (49) below.(49)μ^k,h,λx,z=f^k,h,λx,z+g^k,h,λx,z=Ak,h,λy+Bk,h,λy=Ck,h,λy.


### 3.3. The Selection of the Optimal Number of Knot Points, Oscillation Parameter, and Smoothing Parameter

The MTSFS model is very dependent on the number of the knot points, oscillation parameter, and smoothing parameter. The best model is obtained by using the optimal values of these parameters. In semiparametric and nonparametric regression, there are several methods to obtain the best regression model. One of the popular methods is generalized cross-validation (GCV).Craven and Wahba [8] introduced a modified cross-validation (CV) method called GCV. This method was developed to overcome the shortcomings of complex CV calculations. GCV has several advantages: it is simple and efficient in calculation, invariant to transformation, and does not require information about variant. In addition, GCV has better asymptotic properties than other methods [9, 10].In this study, the value of GCV is a criterion that can be used to determine the best model from variety of knots, oscillations, and smoothing parameters. The criterion for selecting the best model includes taking the model with the lowest GCV value. The modified GCV method of the MTSFS model for longitudinal data is stated as follows.(50)GCVk,h,λ=N−1I−Ck,h,λy2N−1traceI−Ck,h,λ2.The optimal number of knot points, oscillation, and smoothing parameters is obtained by minimizingGCVk,h,λ.


### 3.4. Simulation Study

This section presents the use of the MTSFS model on simulation data to see the performance of the estimators obtained. The simulation was done with two predictor variables, sample sizen = 20, and varied number of time points, T=5,10,15. We considered 20 models for each subject generated from a formula that contains two different functions representing the truncated spline and Fourier series pattern. A polynomial function is used to present the truncated spline, while trigonometric functions are used to present the Fourier series. The predictors are generated from U0,1, and the random errors εit are generated from a multivariate normal distribution.The weight matrix is user-specified [6], and in this paper, we use W=N−1I so that each of the measurements is treated equally. In this study, we use two numbers of knots (K=1 and K=2) and several oscillation parameters (H=1,2,3,4). These simulation studies based on GCV are presented in Table 1.Table 1
Comparison of the values of GCV based on variations in number of time points.
Number of time pointsNumber of knotsNumber of oscillationsλGCV51110.00004.932628.77554.934930.81634.0640410.00004.92112110.00005.905120.20415.912830.81634.7601410.00005.901310110.20413.8454210.00003.805130.81632.7655410.00003.8482210.20414.1491210.00004.106030.81632.9306410.00004.151915110.20413.813327.95923.820730.81632.5196410.00003.8124210.20413.966325.10203.971730.81632.5693410.00003.9660Based on Table1, for the model with T=5, it appears that the smallest GCV value occurs when the oscillation H=3 for K=1 and K=2. The same is also seen at T=10 and T=15. The smallest GCV occurs when the oscillation H=3 for both knot points.In general, the larger the number of time points, the smaller the resulting value of the GCV. Furthermore, the greater the number of knot points, the larger the GCV. Other results show that a larger oscillation is not guaranteed to produce large or small values of the GCV. So, it is necessary to choose the optimum oscillation that produces the smallest GCV.


### 3.5. Application to Real Data

The MTSFS model for longitudinal data obtained has been applied to the stroke patient dataset. These data were taken after an initial study of GCS on stroke patients and the factors that influenced it. The pattern of the relations between the predictors and response followed the characteristics of a truncated spline and a Fourier series. There is a predictor with the form of a truncated spline, which is changing in certain subintervals. The other predictor has the form of a Fourier series, with a repeating pattern.Stroke is a noncommunicable disease. The number of those who suffer from it continues to increase in the world. In 2013, stroke was the second leading cause of death globally (11.8% of all deaths) after ischemic heart disease (14.8% of all deaths). Besides, stroke is the third cause of disability, namely, 4.5% of all causes of disability [11]. Based on the Global Burden of Disease (GBD) Study 2016, the estimated global lifetime risk of stroke for those aged 25 years or older almost reached 25% [12]. The global prevalence of stroke in 2017 was 104.2 million people. In addition, age-standardized stroke prevalence rates were highest in Eastern Europe, North Africa, the Middle East, and Central and East Asia. Several countries in Europe, North Africa, and Central Asia have the highest rates of stroke mortality. Indonesia is one of the countries with the highest death rate due to stroke. According to the Indonesia Basic Health Survey 2007, stroke was the highest cause of death (15.4%) [13]. The prevalence of stroke in Indonesia in 2013 was 7% and increased to 10.9%, according to the Indonesia Basic Health Survey 2018. Furthermore, the stroke prevalence for those aged 15 years or older was 10.85%.Stroke patients often experience head injuries due to falls. Trauma or head injury requires supervision to ensure further medical treatment. GCS was initially used to assess consciousness level after head injury and is now used in the medical field of both acute and trauma patients. According to Champion [14], trauma severity assessment is used to measure the severity of the injury and describe the patient’s case's severity. An injury to the head will activate the immune response and release of interleukins by activating white blood cells and increasing the temperature. Injury to the hypothalamus also increases interference with body temperature regulation [15]. Fever is generally defined as an increase in body temperature above average and has been identified as one of the causes of worsening head injury [16]. Based on previous studies, the pulse was found to be dysfunctional in patients with severe brain injury. Changes in pulse rate are associated with increased mortality in brain injury patients. Therefore, changes in patients’ pulse rates with severe brain injury should be carefully observed [17].This article uses GCS as a response variable with 20 stroke patients (n=20) with 14-day (T=14) measurements for each subject. Simultaneously, the predictor variables are body temperature (BT) and pulse rate (PR). The partial relationship between GCS and each predictor variable is illustrated in Figure 1. The plot in Figure 1 shows a diﬀerent pattern for each predictor. For this reason, we propose the MTSFS model for longitudinal data. Body temperature will be approximated by a truncated spline estimator, while the pulse rate will be approximated by a Fourier series estimator.Figure 1
Scatterplot of 20 stroke patients. (a) GCS and body temperature. (b) GCS and pulse rate.
(a)(b)Based on Table2, some scenarios were performed to compare the effectiveness of the proposed model, i.e., using a single estimator and the mixed estimator. By using the GCV formula in Section 3.3, we selected the best model according to the minimum GCV criterion. As shown in Table 2, modeling GCS in stroke patients produces the smallest GCV value of 45.3517, achieved at when the number of knots is K=1 and number of oscillation is H=3 with λ=10. This study revealed that the MTSFS model is better than using a single estimator. This model yields a RMSE = 2.43. The important result was that the number of knots and the number of the oscillation parameters of the best model in the simulation study are in line with those for the case study. Similar to the simulation study, the best model had one knot and the oscillation parameter H equal to three.Table 2
Summary of case study results.
ModelNumber of knotsNumber of oscillationsλGCVTruncated spline estimator1——391858.1Fourier series estimator—1—323174.42161621.03101785.9475480.3Mixed truncated spline and Fourier series110.204148.5434210.000047.4535310.000045.351749.387845.3519210.204151.6948210.000050.687239.387848.734547.142948.7347Figure2 presents a comparison between the response variable (red line) and fitted values (blue line) using the MTSFS model. From the graph, we can see that some of the fitted values resemble the pattern of the real data; although some do not, there are deviations. Despite its limitations, the present study certainly adds to our understanding of the theory of the new mixed estimator in nonparametric regression for longitudinal data.Figure 2
Comparison between actual and fitted values for each subject.


## 3.1. MTSFS Model for Longitudinal Data

Some lemmas and theorems are provided to obtain an MTSFS model for longitudinal data. Lemma1 presents the goodness of fit of the Fourier series component, Lemma 2 presents the penalty component, and Lemma 3 gives the solution for the truncated spline component. Theorem 1 presents the PWLS optimization, and Theorem 2 presents WLS optimization. The results are summarized as follows.Lemma 1.
If the Fourier series component in equation (2) is given by ∑k=1qgkizkit, then the goodness of fit can be formulated as follows:(3)N−1y∗−Zc′Wy∗−Zc,where N=n×T, W is the nT×nT weighting matrix, y∗=y−f, and(4)y∗=y1∗′y2∗′…yn∗′′,yi∗=yi1∗yi2∗…yiT∗′,f=f1′f2′…fn′′,fi=∑j=1pfjixji1,∑j=1pfjixji2,…,∑j=1pfjixjiT′,Z=Z10…00Z2…0⋮⋮⋱⋮00…Zn,c=c1c2⋮cn,Zi=z1i112cosz1i1cos2z1i1…cosHz1i1…zqi112coszqi1cos2zqi1…cosHzqi1z1i212cosz1i2cos2z1i2…cosHz1i2…zqi212coszqi2cos2zqi2…cosHzqi2⋮⋮⋮⋮⋰⋮⋰⋮⋮⋮⋮⋰⋮z1iT12cosz1iTcos2z1iT…cosHz1iT…zqiT12coszqiTcos2zqiT…cosHzqiT,ci=d1ic01ic11ic21i…cH1i…dqic0qic1qic2qi…cHqi′.Proof.
The functiongki is assumed to be unknown and contained in the space C0,π. The function is approximated using Fourier series with a trend, modified from Bilodeau [3]:(5)gkizkit=dkizkit+12c0ki+∑h=1Hchkicoshzkit.
Equation (5) can be written in matrix form(6)g=g1′g2′⋮gn′=Z1c1Z2c2⋮Zncn=Z10…00Z2…0⋮⋮⋱⋮00…Znc1c2⋮cn=Zc,where Z is a nT×2+Hnq matrix and c is a 2+Hnq vector.
The nonparametric model in equation (2) can be rewritten as(7)yit−∑j=1pfjixjit=∑k=1qgkizkit+εit,yit∗=∑k=1qgkizkit+εit,i=1,2,…,n,t=1,2,…,T.
Then, the goodness of fit for equation (7) can be written as(8)N−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tgkizkit2=N−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tdkizkit+12c0ki+∑h=1Hchkicoshzkit2.
The model in equation (8) can be written in matrix form:(9)N−1Wy∗−Zc2=N−1y∗−Zc′Wy∗−Zc.Lemma 2.
If the penalty component is given,(10)∑k=1qλk∫0π2πgki″zk2dzk,then(11)∑k=1qλk∫0π2πgki″zk2dzk=c′Dλc,where(12)c=c1c2…cn′,ci=d1ic01ic11ic21i…cH1i…dqic0qic1qic2qi…cHqi′,Dλ=D1λ0…00D2λ…0⋮⋮⋱⋮00…Dnλ,Diλ=d10…00d2…0⋮⋮⋱⋮00…dq,Ik=λk14λk24…λkH4,dk=diag0,0,Ik,k=1,2,…,q.Proof.
Regarding equation (5), we define(13)gki″zk=ddzkddzkdkizkit+12c0ki+∑h=1Hchkicoshzkit=−∑h=1Hh2chkicoshzkit.
Consequently,(14)Pkgk=∫0π2π∑h=1Hh2chkicoshzkit2dzk=2π∫0π∑h=1Hh2chkicoshzkit2+2∑h<mHh2chkicoshzkitm2cmkicosmzkitdzk.
Let(15)A=2π∑h=1H∫0πh2chkicoshzkit2dzk,B=2π2∑h<lH∫0πh2chkicoshzkitl2clkicoslzkitdzk.
Thus,(16)Pkgk=A+B=∑h=1Hh4chki2.
The penalty in equation (16) can be rewritten in matrix form(17)c′Dλc.
Regarding the goodness of fit in Lemma1 and the penalty component in Lemma 2, we obtain the PWLS optimization as(18)N−1∑i=1n∑t=1Twityit∗−∑k=1qgkizkit2+∑k=1qλk∫0π2πgki″zk2dzk,0<λk<∞.
Equation (18) can be rewritten in the form(19)N−1y∗−Zc′Wy∗−Zc+c′Dλc.Theorem 1.
If the goodness of fit is as in Lemma1 and the penalty component is as in Lemma 2, then the Fourier series component obtained by minimizing the PWLS in equation (18) is(20)g^k,h,λx,z=Ly∗,where y∗=y−f and L=ZZ′WZ+NDλ−1Z′W.Proof.
The optimization in equation (18) can be written as(21)Mingk∈C0,πN−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tgkizkit2+∑k=1qλk∫0π2πgk″zk2dzk.
Equation (21) can be rewritten in the form(22)Minc∈C0,πN−1y∗−Zc′Wy∗−Zc+c′Dλc=Minc∈C0,πQc.
We obtain(23)Qc=N−1y∗′Wy∗−2N−1c′Z′Wy∗+N−1c′Z′WZc+c′Dλc.
The completion of the optimization in (23) is obtained by taking the partial derivative of Qc with respect to c and setting it equal to zero, i.e.,(24)∂Qc∂c=0,giving the result(25)c^=Z′WZ+NDλ−1Z′Wy∗.
By substituting (25) into (6), we obtain(26)g^k,h,λx,z=Zc^=ZZ′WZ+NDλ−1Z′Wy∗=Ly∗.
The nonparametric regression model in equation (7) can be written as(27)y∗=g^k,h,λx,z=Ly∗,where y∗=y−f and L=ZZ′WZ+NDλ−1Z′W.Lemma 3.
If the truncated spline component in equation (2) is ∑j=1pfjixjit, then the WLS is(28)I−Ly−I−LMγ′WI−Ly−I−LMγ,where L=ZZ′WZ+NDλ−1Z′W.Proof.
The functionfji is a linear truncated spline function with s knot for each xj,j=1,2,…,p:(29)fjixjit=αjixjit+∑u=1sβujixjit−Kuji+,where xjit−Kuji+=xjit−Kuji,xjit≥Kuji,0,xjit<Kuji.
Equation (29) can be rewritten in matrix form(30)f=X|Sα¯β=Mγ,where(31)X=X10⋯00X2⋯0⋮⋮⋱⋮00⋯Xn,S=S10⋯00S2⋯0⋮⋮⋱⋮00⋯Sn,α=α1α2⋮αn,β=β1β2⋮βn,Xi=x1i1x2i1…xpi1x1i2x2i2…xpi2⋮⋮⋱⋮x1iTx2iT…xpiT,αi=α1iα2i⋮αpi,Si=x1i1−K11i+1…x1i1−Ks1i+1…xpi1−K1pi+1…xpi1−Kspi+1x1i2−K11i+1…x1i2−Ks1i+1…xpi2−K1pi+1…xpi2−Kspi+1⋮⋱⋮⋱⋮⋱⋮x1iT−K11i+1…x1iT−Ks1i+1…xpiT−K1pi+1…xpiT−Kspi+1,βi=β11i…βs1i…β1pi…βspiT.
The MTSFS model for longitudinal data in equation (1) can be written in the form(32)y=f+g+ε.
By substituting (26) into (32), we obtain(33)y=f+Ly∗+ε.
To obtain the estimator of truncated spline component, equation (33) can be written as(34)y−Ly∗=f+ε,y−Ly−f=f+ε,I−Ly=I−Lf+ε.
Equation (34) can be rewritten as(35)I−Ly=I−LMγ+ε.
Thus,(36)ε=I−Ly−I−LMγ.
As a consequence, the WLS is given by(37)ε′ε=I−Ly−I−LMγ′WI−Ly−I−LMγ.Theorem 2.
Suppose the WLS is as in Lemma3. Then, the mixed estimator obtained by minimizing WLS is(38)f^k,h,λx,z=MJ−1Ky,where J=2−L′WLM−WM and K=L′−IWI−L.Proof.
The WLS optimization in Lemma3 can be written as(39)MinγI−Ly−I−LMγ′WI−Ly−I−LMγ=MinγQγ.
We obtain(40)Qγ=y′Wy+y′L′WLy+γ′M′L′WLMγ+γ′M′WMγ−2y′L′Wy+2γ′M′L′Wy+−2N−1γ′M′Wy−2N−1γ′M′L′WLy+2N−1γ′M′WLy−2N−1γ′M′WLMγ.
The complete optimization is obtained by setting equal to zero the partial derivative ofQγ with respect to γ, that is,(41)∂Qγ∂γ=0,giving the result(42)γ^=2−L′WLM−WM−1L′−IWI−Ly=J−1Ky.
By substitutingγ^ into equation (30), we obtain(43)f^k,h,λx,z=Mγ^=MJ−1Ky=Ak,h,λy.
We obtainc^ by substituting (43) into (25):(44)c^=Z′WZ+NDλ−1Z′Wy∗=Z′WZ+NDλ−1Z′Wy−f^=Z′WZ+NDλ−1Z′Wy−Mγ^=Z′WZ+NDλ−1Z′Wy−MJ−1Ky=Z′WZ+NDλ−1Z′WI−MJ−1Ky.
As a consequence,(45)g^k,h,λx,z=Zc^=ZZ′WZ+NDλ−1Z′WI−MJ−1Ky=LI−MJ−1Ky=Bk,h,λy.
By substitutingγ^ and c^ into the MTSFS model for longitudinal data, we obtain(46)y^=μ^k,h,λx,z=f^k,h,λx,z+g^k,h,λx,z=Mγ^+Zc^=MJ−1Ky+LI−MJ−1Ky=MJ−1K+LI−MJ−1Ky=Ck,h,λy,where J=2−L′WLM−WM, K=L′−IWI−L, and L=ZZ′WZ+NDλ−1Z′W.


## 3.2. The Properties of MTSFS Model for Longitudinal Data

This section provides the MTSFS model’s properties for longitudinal data, i.e., it is biased and linear in observations. The mixed estimator is biased, as proved by(47)Eμ^k,h,λx,z=Ef^k,h,λx,z+g^k,h,λx,z=EAk,h,λy+Bk,h,λy=ECk,h,λy=Ck,h,λEy.Assuming thatCk,h,λ≠I,(48)Eμ^k,h,λx,z=Ck,h,λEy=Ck,h,λfk,h,λx,z+gk,h,λx,z=Ck,h,λμk,h,λx,z.The result in equation (48) showed that the mixed estimator is biased because Eμ^k,h,λx,z≠μk,h,λx,z. Even though the mixed estimator is biased, the mixed estimator is linear in observations proved by equation (49) below.(49)μ^k,h,λx,z=f^k,h,λx,z+g^k,h,λx,z=Ak,h,λy+Bk,h,λy=Ck,h,λy.


## 3.3. The Selection of the Optimal Number of Knot Points, Oscillation Parameter, and Smoothing Parameter

The MTSFS model is very dependent on the number of the knot points, oscillation parameter, and smoothing parameter. The best model is obtained by using the optimal values of these parameters. In semiparametric and nonparametric regression, there are several methods to obtain the best regression model. One of the popular methods is generalized cross-validation (GCV).Craven and Wahba [8] introduced a modified cross-validation (CV) method called GCV. This method was developed to overcome the shortcomings of complex CV calculations. GCV has several advantages: it is simple and efficient in calculation, invariant to transformation, and does not require information about variant. In addition, GCV has better asymptotic properties than other methods [9, 10].In this study, the value of GCV is a criterion that can be used to determine the best model from variety of knots, oscillations, and smoothing parameters. The criterion for selecting the best model includes taking the model with the lowest GCV value. The modified GCV method of the MTSFS model for longitudinal data is stated as follows.(50)GCVk,h,λ=N−1I−Ck,h,λy2N−1traceI−Ck,h,λ2.The optimal number of knot points, oscillation, and smoothing parameters is obtained by minimizingGCVk,h,λ.


## 3.4. Simulation Study

This section presents the use of the MTSFS model on simulation data to see the performance of the estimators obtained. The simulation was done with two predictor variables, sample sizen = 20, and varied number of time points, T=5,10,15. We considered 20 models for each subject generated from a formula that contains two different functions representing the truncated spline and Fourier series pattern. A polynomial function is used to present the truncated spline, while trigonometric functions are used to present the Fourier series. The predictors are generated from U0,1, and the random errors εit are generated from a multivariate normal distribution.The weight matrix is user-specified [6], and in this paper, we use W=N−1I so that each of the measurements is treated equally. In this study, we use two numbers of knots (K=1 and K=2) and several oscillation parameters (H=1,2,3,4). These simulation studies based on GCV are presented in Table 1.Table 1
Comparison of the values of GCV based on variations in number of time points.
Number of time pointsNumber of knotsNumber of oscillationsλGCV51110.00004.932628.77554.934930.81634.0640410.00004.92112110.00005.905120.20415.912830.81634.7601410.00005.901310110.20413.8454210.00003.805130.81632.7655410.00003.8482210.20414.1491210.00004.106030.81632.9306410.00004.151915110.20413.813327.95923.820730.81632.5196410.00003.8124210.20413.966325.10203.971730.81632.5693410.00003.9660Based on Table1, for the model with T=5, it appears that the smallest GCV value occurs when the oscillation H=3 for K=1 and K=2. The same is also seen at T=10 and T=15. The smallest GCV occurs when the oscillation H=3 for both knot points.In general, the larger the number of time points, the smaller the resulting value of the GCV. Furthermore, the greater the number of knot points, the larger the GCV. Other results show that a larger oscillation is not guaranteed to produce large or small values of the GCV. So, it is necessary to choose the optimum oscillation that produces the smallest GCV.


## 3.5. Application to Real Data

The MTSFS model for longitudinal data obtained has been applied to the stroke patient dataset. These data were taken after an initial study of GCS on stroke patients and the factors that influenced it. The pattern of the relations between the predictors and response followed the characteristics of a truncated spline and a Fourier series. There is a predictor with the form of a truncated spline, which is changing in certain subintervals. The other predictor has the form of a Fourier series, with a repeating pattern.Stroke is a noncommunicable disease. The number of those who suffer from it continues to increase in the world. In 2013, stroke was the second leading cause of death globally (11.8% of all deaths) after ischemic heart disease (14.8% of all deaths). Besides, stroke is the third cause of disability, namely, 4.5% of all causes of disability [11]. Based on the Global Burden of Disease (GBD) Study 2016, the estimated global lifetime risk of stroke for those aged 25 years or older almost reached 25% [12]. The global prevalence of stroke in 2017 was 104.2 million people. In addition, age-standardized stroke prevalence rates were highest in Eastern Europe, North Africa, the Middle East, and Central and East Asia. Several countries in Europe, North Africa, and Central Asia have the highest rates of stroke mortality. Indonesia is one of the countries with the highest death rate due to stroke. According to the Indonesia Basic Health Survey 2007, stroke was the highest cause of death (15.4%) [13]. The prevalence of stroke in Indonesia in 2013 was 7% and increased to 10.9%, according to the Indonesia Basic Health Survey 2018. Furthermore, the stroke prevalence for those aged 15 years or older was 10.85%.Stroke patients often experience head injuries due to falls. Trauma or head injury requires supervision to ensure further medical treatment. GCS was initially used to assess consciousness level after head injury and is now used in the medical field of both acute and trauma patients. According to Champion [14], trauma severity assessment is used to measure the severity of the injury and describe the patient’s case's severity. An injury to the head will activate the immune response and release of interleukins by activating white blood cells and increasing the temperature. Injury to the hypothalamus also increases interference with body temperature regulation [15]. Fever is generally defined as an increase in body temperature above average and has been identified as one of the causes of worsening head injury [16]. Based on previous studies, the pulse was found to be dysfunctional in patients with severe brain injury. Changes in pulse rate are associated with increased mortality in brain injury patients. Therefore, changes in patients’ pulse rates with severe brain injury should be carefully observed [17].This article uses GCS as a response variable with 20 stroke patients (n=20) with 14-day (T=14) measurements for each subject. Simultaneously, the predictor variables are body temperature (BT) and pulse rate (PR). The partial relationship between GCS and each predictor variable is illustrated in Figure 1. The plot in Figure 1 shows a diﬀerent pattern for each predictor. For this reason, we propose the MTSFS model for longitudinal data. Body temperature will be approximated by a truncated spline estimator, while the pulse rate will be approximated by a Fourier series estimator.Figure 1
Scatterplot of 20 stroke patients. (a) GCS and body temperature. (b) GCS and pulse rate.
(a)(b)Based on Table2, some scenarios were performed to compare the effectiveness of the proposed model, i.e., using a single estimator and the mixed estimator. By using the GCV formula in Section 3.3, we selected the best model according to the minimum GCV criterion. As shown in Table 2, modeling GCS in stroke patients produces the smallest GCV value of 45.3517, achieved at when the number of knots is K=1 and number of oscillation is H=3 with λ=10. This study revealed that the MTSFS model is better than using a single estimator. This model yields a RMSE = 2.43. The important result was that the number of knots and the number of the oscillation parameters of the best model in the simulation study are in line with those for the case study. Similar to the simulation study, the best model had one knot and the oscillation parameter H equal to three.Table 2
Summary of case study results.
ModelNumber of knotsNumber of oscillationsλGCVTruncated spline estimator1——391858.1Fourier series estimator—1—323174.42161621.03101785.9475480.3Mixed truncated spline and Fourier series110.204148.5434210.000047.4535310.000045.351749.387845.3519210.204151.6948210.000050.687239.387848.734547.142948.7347Figure2 presents a comparison between the response variable (red line) and fitted values (blue line) using the MTSFS model. From the graph, we can see that some of the fitted values resemble the pattern of the real data; although some do not, there are deviations. Despite its limitations, the present study certainly adds to our understanding of the theory of the new mixed estimator in nonparametric regression for longitudinal data.Figure 2
Comparison between actual and fitted values for each subject.


## 4. Conclusions

This article presented a new mixed estimator for nonparametric regression models for longitudinal data, combining the truncated spline estimator and Fourier series to obtain better estimation results if the predictors have different data patterns. A new two-step method for estimating the parameters used penalized weighted least squares and weighted least squares.In both the simulation study and case study, the best model was selected by using the minimum generalized cross-validation (GCV). A higher number of knots or oscillation parameter does not produce a high GCV. Therefore, several combinations of knots and oscillations had to be tried to determine the best model. Interestingly, the number of knots and oscillations of the best model is the same for the simulation study as for the case study with real data: one knot and three oscillations.The proposed model estimator is a useful alternative for estimating nonparametric regression curves for longitudinal data. The presented results suggest that the next study can use different weighting matrices and more knots, so that researchers can compare the results for improving the model performance. Another possible line of further research could be simulation studies using another function to evaluate the performance of the proposed model.



---

*Source: 3909401-2021-11-15.xml*