# Nonparametric Regression Model for Longitudinal Data with Mixed Truncated Spline and Fourier Series

Authors: Made Ayu Dwi Octavanny; I. Nyoman Budiantara; Heri Kuswanto; Dyah Putri Rahmawati
Journal: Abstract and Applied Analysis
 (2020)

---

## Abstract

Existing literature in nonparametric regression has established a model that only applies one estimator to all predictors. This study is aimed at developing a mixed truncated spline and Fourier series model in nonparametric regression for longitudinal data. The mixed estimator is obtained by solving the two-stage estimation, consisting of a penalized weighted least square (PWLS) and weighted least square (WLS) optimization. To demonstrate the performance of the proposed method, simulation and real data are provided. The results of the simulated data and case study show a consistent finding.


---

## Body

## 1. Introduction

Regression analysis is aimed at modeling the association between the predictor and the response. If the data pattern shows an unknown regression curve, nonparametric regression is used [1]. However, if the form of the regression curve is known, parametric regression can be applied [2]. Additionally, nonparametric regression has high flexibility because the data is expected to find its regression curve estimation form without being influenced by the researcher’s subjectivity [3]. In this study, we have analyzed several models such as kernel, spline [4–7], and Fourier series [8].A spline estimator, which has an excellent ability to handle data with changes at subspecified intervals [9], was obtained using penalized least square optimization [10] and the Bayesian approach [11]. A spline estimator can be applied for cross-sectional data as well as longitudinal data. Additionally, several studies on nonparametric regression for longitudinal data have been addressed using kernel estimator [12, 13], generalized spline regression [14], and mixed-effects model [7]. Fourier series, which is useful to explain curves that show sine and cosine waves, is generally used if the data pattern is unknown and there is a tendency to iterate.A considerable amount of research has used only one estimator for each predictor. However, because each predictor can have a different pattern, it was proposed to develop a mixed estimator. Recently, Sudiarsa et al. [15] discussed a study of the mixed estimator of the truncated spline and Fourier series. The study, which only discussed cross-sectional data, did not obtain a model for each subject as it did not include longitudinal data. Consequently, this study cannot be used to investigate response behavior based on the time change.Although some research has been carried out on a mixed estimator, no studies have explored multisubject data so far. This paper proposes a new methodology for a mixed estimator of the truncated spline and Fourier series in the nonparametric regression for longitudinal data. This study addresses the gap in previous research by obtaining a mixed estimator of the truncated spline and Fourier series in the nonparametric regression for longitudinal data and applying it to simulated data and a case study.This study is organized as follows. We briefly explain the materials and methods used in our study in Section2. Section 3 consists of three subsections: the developed theory, simulation study, and case study. We present the developed nonparametric regression theory for longitudinal data with a mixed estimator of the truncated spline and Fourier series with two-stage estimation in Section 3.1. In Section 3.2, we conduct a simulation study based on the developed theory to assess the proposed estimator’s behavior. To illustrate the applicability of the model, we use a dataset of patients with pulmonary tuberculosis in Section 3.3. Section 4 presents the conclusion.


## 2. Materials and Methods

Longitudinal data hasn independent subjects and T observations for each subject. Given paired data x1it,⋯,xpit,z1it,⋯,zqit,yit, which consist of p and q predictors with n subjects, each subject has T observations. The relationship between x1it,⋯,xpit,z1it,⋯,zqit and yit, which followed a nonparametric regression model for longitudinal data, is as follows:
(1)yit=μitx1it,⋯,xpit,z1it,⋯,zqit+εit,i=1,2,⋯,n,t=1,2,⋯,T.Each regression curve is additive so that the model can be expressed as
(2)μitx1it,⋯,xpit,z1it,⋯,zqit=∑j=1pfjixjit+∑k=1qgkizkit,i=1,2,⋯,n,t=1,2,⋯,T.∑j=1pfjixjit is the truncated spline component and ∑k=1qgkizkit is the Fourier series component.This study’s first objective is to obtain the mixed estimator of the truncated spline and Fourier series in nonparametric regression for longitudinal data. To achieve this, we propose a two-stage estimation method. The first stage is estimating the components of the Fourier series using the penalized weighted least square (PWLS) method. The second stage is estimating the truncated spline component using the weighted least square (WLS) method. For the second goal, that is, a simulation study, we generate functions that meet the truncated spline and Fourier series characteristics. In the third step, we apply the developed theory to a dataset of patients with pulmonary tuberculosis.


## 3. Results and Discussion

### 3.1. Mixed Model of Truncated Spline and Fourier Series with Two-Stage Estimation

Lemmas and theorems are used to obtain a nonparametric regression model for longitudinal data with a mixed estimator. The regression curve component that is approximated by the Fourier series estimator is presented in Lemma1 and the penalty component for the Fourier series function is presented in Lemma 2. Following the PWLS form in Lemma 3, we estimate the Fourier series component by using PWLS in Theorem 4. The regression curve component that is approximated by the truncated spline estimator is presented in Lemma 5, and we estimate the truncated spline component using the WLS method in Theorem 6. The results are summarized as follows:Lemma 1.
Ifgkizkit is approached by the Fourier series function, then the goodness of fit is
(3)N−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tgkizkit2=N−1y∗−Zc′Wy∗−Zc,where N=n×T and W is the nT×nT weighting matrix.Proof.
The regression curvegkizkit is a regression curve of an unknown shape and is contained in continuous space C0,π. The component of gkizkit,k=1,2,⋯,q, in Equation (2) is approximated by the Fourier series function with the trend line as follows:
(4)gkizkit=dkizkit+12c0ki+∑h=1HchkicoshzkitIf the regression curveg in Equation (4) involves only one predictor, then it can be written as follows:
(5)g=g1′g2′⋯gn′′,(6)gi=gkizki1gkizki2⋯gkizkiT′,i=1,2,⋯,n.By using Equation (4), Equation (6) can be written in the form of a matrix as follows:
(7)gi=gkizki1gkizki2⋮gkizkiT=dkizki1+12c0ki+c1kicos1zki1+c2kicos2zki1+⋯+cHkicosHzki1dkizki2+12c0ki+c1kicos1zki2+c2kicos2zki2+⋯+cHkicosHzki2⋮dkizkiT+12c0ki+c1kicos1zkiT+c2kicos2zkiT+⋯+cHkicosHzkiT=zki112cos1zki1cos2zki1⋯cosHzki1zki212cos1zki2cos2zki2⋯cosHzki2⋮⋮⋮⋮⋱⋮zkiT12cos1zkiTcos2zkiT⋯cosHzkiTdkic0kic1kic2ki⋮cHki=Zkicki.The Fourier series function in the nonparametric regression component for longitudinal data with predictorq can be expressed in the following form:
(8)gi=∑k=1qZkicki=Z1ic1i+Z2ic2i+⋯+Zqicqi=Zici.So,(9)g=g1′g2′⋮gn′=Z1c1Z2c2⋮Zncn=Z10⋯00Z2⋯0⋮⋮⋱⋮00⋯Znc1c2⋮cn=Zc.Z is a nT×2+Hnq matrix as follows:
(10)Z=Z10⋯00Z2⋯0⋮⋮⋱⋮00⋯Zn,with
(11)Zi=z1i112cosz1i1cos2z1i1⋯cosHz1i1z1i212cosz1i2cos2z1i2⋯cosHz1i2⋮⋮⋮⋮⋱⋮z1iT12cosz1iTcos2z1iT⋯cosHz1iT⋯⋯⋯⋯zqi112coszqi1cos2zqi1⋯cosHzqi1zqi212coszqi2cos2zqi2⋯cosHzqi2⋮⋮⋮⋮⋱⋮zqiT12coszqiTcos2zqiT⋯cosHzqiT.Whereas,c is a 2+Hnq×1 vector given by
(12)c=c1c2⋮cn,ci=d1ic01ic11ic21i⋯cH1i⋯dqic0qic1qic2qi⋯cHqi′.To estimate the Fourier series component, the nonparametric regression model in Equation (1) can be written as
(13)yit−∑j=1pfjixjit=∑k=1qgkizkit+εit,(14)yit∗=∑k=1qgkizkit+εit,i=1,2,⋯,n,t=1,2,⋯,T.The model in Equation (14) can be written in matrix form:
(15)y∗=g+ε.Then, a goodness of fit for the model is formed as follows:(16)N−1∑i=1n∑t=1Twityit∗−∑k=1qgkizkit2,with N=n×T. If the function gkizkit is approached by a Fourier series function as in Equation (4), then the goodness of fit can be presented in the form
(17)N−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tgkizkit2=N−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tdkizkit+12c0ki+∑h=1Hchkicoshzkit2=N−1y∗−Zc′Wy∗−Zc,with W as a weighting matrix for the regression of longitudinal data.Lemma 2.
If the Fourier series is given, then the penalty component is(18)∑k=1qλk∫0π2πgki′′zk2dzk=∑k=1qλk∫0π2πgki′′zk2dzk=c′Dλc.Proof.
The penalty component∑k=1qλk∫0π2/πgki′′zk2dzk in the PWLS optimization based on Equation (4) can be obtained as follows:
(19)gki′′zk=ddzkddzkdkizkit+12c0ki+∑h=1Hchkicoshzkit=−∑h=1Hh2chkicoshzkit.As a result,(20)Pkgk=∫0π2π∑h=1Hh2chkicoshzkit2dzk=2π∫0π∑h=1Hh2chkicoshzkit2+2∑h<mHh2chkicoshzkitm2cmkicosmzkitdzk.To simplify, we defined(21)A=2π∑h=1H∫0πh2chkicoshzkit2dzk,B=2π2∑h<lH∫0πh2chkicoshzkitl2clkicoslzkitdzk.The value ofA will be obtained as follows.
(22)A=2π∑h=1H∫0πh2chkicoshzkit2dzk=1π∑h=1Hh4chki2zkit+2hsinhzkit0π=∑h=1Hh4chki2.Furthermore, the value of B is given by
(23)B=2π2∑h<lH∫0πh2chkicoshzkitl2clkicoslzkitdzk=4π∑h<lHhl2chkiclki∫0πcoshzkitcoslzkitdzk=0.Based on Equations (22) and (23), it can be written as follows:
(24)Pkgk=A+B=∑h=1Hh4chki2.Fori=1,2,⋯,n, we obtained
(25)∑k=1qλkPkgk=∑i=1nci′Dλci=c1′Dλc1+c2′Dλc2+⋯+cn′Dλcn=c′Dλc,where
(26)Dλ=Dλ10⋯00Dλ2⋯0⋮⋮⋱⋮00⋯Dλq,Dλi=00λi14λi24⋯λiH400λi14λi24⋯λiH4⋮⋮⋮⋮⋱⋮00λi14λi24⋯λiH4..Thus, the penalty component can be expressed in a matrix form as follows:(27)∑k=1qλk∫0π2πgki′′zk2dzk=c′Dλc.Lemma 3.
If the goodness of fit component is presented in Lemma1 and the penalty component is given by Lemma 2, then the PWLS is
(28)N−1y∗−Zc′Wy∗−Zc+c′Dλc.In general, PWLS is defined as follows:(29)N−1∑i=1n∑t=1Twityit∗−∑k=1qgkizkit2+∑k=1qλk∫0π2πgki′′zk2dzk,0<λk<∞.Besides, PWLS can be presented in the form of a matrix as follows:(30)N−1y∗−Zc′Wy∗−Zc+c′Dλc.The next step is to obtain a Fourier series estimator in nonparametric regression for longitudinal data derived in Theorem4.Theorem 4.
If paired data are given, which follows the nonparametric regression model for longitudinal data, then the mixed estimator that minimizes PWLS in Lemma3,
(31)mingk∈C0,πN−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tgkizkit2+∑k=1qλk∫0π2πgk′′zk2dzk,is g^k,h,λx,z=Ly∗ with y∗=y−f and L=ZZ′WZ+NDλ−1Z′W.Proof.
The first estimation step in the mixed estimator of the truncated spline and Fourier series in the nonparametric regression model for longitudinal data is performed by estimating the form of the Fourier series estimator by using the PWLS method. The PWLS in Equation (29) can be written in the form of a matrix as follows:
(32)Qc=N−1y∗−Zc′Wy∗−Zc+c′Dλc=N−1y∗′Wy∗−2N−1c′Z′Wy∗+N−1c′Z′WZc+c′Dλc.Next, we complete PWLS optimization using the following steps:(33)minc∈R2+HnqQc=minc∈R2+HnqN−1y∗′Wy∗−2N−1c′Z′Wy∗+N−1c′Z′WZc+c′Dλc.To complete the optimization, the estimators are obtained by performing a partial derivative ofQc concerning c and the results are equaled to zero. The given results are
(34)c^=Z′WZ+NDλ−1Z′Wy∗.By substitutingc^ into Equation (9), we get
(35)g^k,h,λx,z=Zc^=ZZ′WZ+NDλ−1Z′Wy∗=Ly∗.So, the model in Equation (15) can be written as
(36)y∗=g^k,h,λx,z=Ly∗,with y∗=y−f and L=ZZ′WZ+NDλ−1Z′W.Lemma 5.
Iffjixjit is approached with the truncated spline function, then the WLS is
(37)I−Ly−I−LMγ′WI−Ly−I−LMγ.W is a nT×nT weighting matrix.Proof.
fjixjit,j=1,2,⋯,p, is a truncated spline estimator component. The component regression curve f is a linear truncated spline function defined as follows:
(38)fjixjit=αjixjit+∑u=1sβujixjit−Kuji+,with truncated function
(39)xjit−Kuji+=xjit−Kuji,xjit≥Kuji,0,xjit<Kuji.If the component of the regression curvef involves one predictor, then it can be written as follows:
(40)f=f1′f2′⋯fn′′,(41)fi=∑j=1pfjixji1∑j=1pfjixji2⋯∑j=1pfjixjiT′,i=1,2,⋯,n.By using Equation (41), it can be described in the form of a matrix as follows:
(42)fi=fjixji1fjixji2⋮fjixjiT=xji1xji2⋮xjiTαji+xji1−K1ji+1xji1−K2ji+1⋯xji1−Ksji+1xji2−K1ji+1xji2−K2ji+1⋯xji2−Ksji+1⋮⋮⋱⋮xjiT−K1ji+1xjiT−K2ji+1⋯xjiT−Ksji+1β1jiβ2ji⋮βsji=Xjiαji+Sjiβji.The truncated spline function in the nonparametric regression component for longitudinal data withp predictors can be expressed in the following form:
(43)fi=∑j=1pXjiαji+Sjiβji=Xiαi+Siβi.So,(44)f=f1′f2′⋮fn′=X1α1X2α2⋮Xnαn+S1β1S2β2⋮Snβn=X10⋯00X2⋯0⋮⋮⋱⋮00⋯Xnα1α2⋮αn+S10⋯00S2⋯0⋮⋮⋱⋮00⋯Snβ1β2⋮βn=Xα+Sβ.Furthermore, it can be written as follows:(45)f=X10⋯00X2⋯0⋮⋮⋱⋮00⋯XnS10⋯00S2⋯0⋮⋮⋱⋮00⋯Snα1α2⋮αn¯β1β2⋮βn=XSα¯β=Mγ.M is a nT×1+snp matrix and γ is a 1+snp×1 vector.The mixed model of nonparametric regression for longitudinal data in Equation (1) can be written in a matrix as follows:
(46)y=f+g+ε.By substituting Equation (36) with Equation (46), we get
(47)y=f+Ly∗+ε.The truncated spline component can be written as(48)y−Ly∗=f+ε,y−Ly−f=f+ε,y−Ly+Lf=f+ε,y−Ly=f−Lf+ε,I−Ly=I−Lf+ε.Substituting Equation (45) with Equation (48), we obtain
(49)I−Ly=I−LMγ+ε.If the functionfjixjit is approximated by the truncated spline function as in Equation (38), then
(50)ε=I−Ly−I−LMγ.Therefore, we obtain the WLS by(51)ε′ε=I−Ly−I−LMγ′WI−Ly−I−LMγ,where W is a weighting matrix for the regression of longitudinal data. Next, the truncated spline estimator in the nonparametric regression for longitudinal data is derived in Theorem 6.Theorem 6.
If paired data which follows the nonparametric regression model for longitudinal data is given, then the mixed estimator that minimizes WLS in Lemma5,
(52)minγ∈R1+snpI−Ly−I−LMγ′WI−Ly−I−LMγ,is f^k,h,λx,z=MJ−1Ky with J=2−L′WLM−WM and K=L′−IWI−L..Proof.
The second estimation stage of the mixed estimator of the truncated spline and Fourier series in the nonparametric regression model for longitudinal data is performed using the WLS method. The estimator can be obtained by completing the WLS optimization as follows:(53)minγ∈R1+snpI−Ly−I−LMγ′WI−Ly−I−LMγ=minγ∈R1+snpQγ.The estimators are obtained by performing a partial derivative ofQγ and the results are equaled to zero. The partial derivative results are as follows:
(54)γ^=2−L′WLM−WM−1L′−IWI−Ly=J−1Ky,where J=2−L′WLM−WM and K=L′−IWI−L.By substitutingγ^ into the form of a truncated spline estimator component as in Equation (45), we obtain
(55)f^k,h,λx,z=Mγ^=MJ−1Ky,(56)f^k,h,λx,z=Ak,h,λy.After obtainingγ^, we obtain c^ by substituting Equation (56) with Equation (34).
(57)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,c^=Z′WZ+NDλ−1Z′WI−MJ−1Ky.So,(58)g^k,h,λx,z=Zc^=ZZ′WZ+NDλ−1Z′WI−MJ−1Ky=LI−MJ−1Ky,g^k,h,λx,z=Bk,h,λy.By substitutingγ^ and c^ into the mixed estimator of the truncated spline and the Fourier series in the nonparametric regression for longitudinal data, the following estimation results are obtained.
(59)μ^k,h,λ=f^k,h,λ+g^k,h,λ=Mγ^+Zc^=MJ−1Ky+LI−MJ−1Ky=MJ−1K+LI−MJ−1Ky,μ^k,h,λ=Ck,h,λy,with L=ZZ′WZ+NDλ−1Z′W, J=2−L′WLM−WM, and K=L′−IWI−L.The mixed estimator depends on the optimum knot point, oscillation parameter, and smoothing parameter. To obtain the best model, it is essential to select the optimum parameter. One of the criteria to select the optimum parameter is the generalized cross-validation (GCV) method [11]. The GCV function of the nonparametric regression model for longitudinal data is as follows:
(60)GCVKopt,Hopt,λopt=N−1I−Ck,h,λ2N−1traceI−Ck,h,λ2.The optimum knot point, oscillation parameter, and smoothing parameter are obtained by solving the minimum optimization, as presented in Equation (60).


### 3.2. Simulation Study

To demonstrate the performance of the proposed method, we created one sample sizen=10 with t=13. For the simulation study, we considered ten models for each subject. The models are generated from the formula that contains two different functions to represent the truncated spline and Fourier series pattern. A polynomial function is used to present the truncated spline, while a trigonometry function is used to present the Fourier series. Additionally, xit and zit are generated from U0,1 distribution, and random errors εit are generated from a multivariate normal distribution.Using ten subjects and two predictors, the formula for generated data is stated as follows:(61)yit=fitxit+gitzit+εit,i=1,2,⋯,10,t=1,2,⋯,13.The simulation study is applied based on these models, as shown in Table1.Table 1
Model for the simulation study.
SubjectModel1y1t=8.25x1t−11−x1t2+sin2πz1t+cos2πz1t+ε1t2y2t=9.55x2t−11−x2t2+sin2πz2t+cos2πz2t+ε2t3y3t=6.75x3t−11−x3t2+sin2πz3t+cos2πz3t+ε3t4y4t=3.25x4t−11−x4t2+sin2πz4t+cos2πz4t+ε4t5y5t=5.55x5t−11−x5t2+sin2πz5t+cos2πz5t+ε5t6y6t=7.75x6t−11−x6t2+sin2πz6t+cos2πz6t+ε6t7y7t=4.25x7t−11−x7t2+sin2πz7t+cos2πz7t+ε7t8y8t=8.55x8t−11−x8t2+sin2πz8t+cos2πz8t+ε8t9y9t=9.75x9t−11−x9t2+sin2πz9t+cos2πz9t+ε9t10y10t=4.75x10t−11−x10t2+sin2πz10t+cos2πz10t+ε10tFigure1 illustrates the partial relationship between the response and each predictor variable. It can be seen that the relationship between predictor x and the response for each subject tends to change at certain subintervals, which is suitable for the truncated spline estimator. The relationship between z and the response for each subject has a repetitive pattern with a particular trend line, which is suitable for the Fourier series estimator.Figure 1
Scatterplot of simulated data for each subject and predictor: (a) truncated spline; (b) Fourier series.
(a)(b)Wu and Zhang [7] stated that a regression’s performance strongly depends on good knot locations and a good choice of the number of knots. In general, the number of knots is smaller than the sample size n. Considering the scatterplot of simulated data and computational convenience, our study uses three knots K=1,2,3 and three oscillation parameters H=1,2,3. To choose the optimum parameter, we use the minimum GCV criteria. Table 2 provides a summary of the GCV for varying knots and oscillations. What is remarkable is that using one knot and one oscillation with λ=0.033, we obtain the best model with the lowest GCV, 3.112. This model yields satisfactory results with a root mean square error (RMSE) of 0.837.Table 2
Summary of simulation results.
KnotsOscillationλGCV110.0333.11221.7153.14230.9673.116210.0333.19721.8003.20930.9503.158310.0013.45220.1003.47130.0013.441


### 3.3. Case Study

After conducting the simulation, we applied the proposed model to the case to confirm the results of the previous simulation. The data for this research was obtained from a study conducted by Fernandes and Solimun [16], that is, patients with pulmonary tuberculosis disease. Pulmonary tuberculosis is a contagious disease caused by Mycobacterium tuberculosis, which can attack various organs, particularly the lungs. This disease is typical among women in their productive years (ages 15-50 years). The World Health Organization (WHO) declared tuberculosis a global emergency in 1992 [17]. WHO report in 2013 stated that there were 8.6 million tuberculosis cases in 2012, of which 40% of the cases were in Southeast Asia. In a further report, Indonesia was noted as the country with the second-largest number of cases, 2.8 million, in 2015.This study’s dataset consists of four patientsn=4 that represent radiological images of the thorax (stadium), which are minimal lesion, mod advance, far advance, and KP Miller. suPAR level as a response y with mL units in several observation periods every two weeks for six months of treatment t=13 is observed. The predictor variables are the erythrocyte sedimentation rate x with mm/hour units and body mass index z with kg/m3 units.The partial relationship between the suPAR level and each predictor variable for each subject is presented in Figure2. There were changes in data patterns in the four subjects observed for six months with measurements taken every two weeks. The plot in Figure 2 shows a different pattern for each predictor. For this reason, we propose a nonparametric regression approach based on a mixed estimator for longitudinal data. The erythrocyte sedimentation rate x will be approached by a truncated spline estimator, while the body mass index z will be approached by a Fourier series estimator.Figure 2
Scatterplot between suPAR levels with each predictor for the subject: (a)minimal lesion; (b) mod advance; (c) far advance; (d) KP Miller.
(a)(b)(c)(d)In this case study, similar to the simulation study, only three knot points and three oscillation parameters were used. From the various knots and oscillation parameter results, we obtained the GCV values listed in Table3. Interestingly, the data in this table is that the minimum GCV achieved one oscillation parameter under the same conditions as the simulation study, that is, one knot point. However, it has a different smoothing parameter, λ=1.25. This model provides a GCV value of 0.239 with a RMSE of 0.197. The knot point location and the results of the parameter estimation for each subject of patients with pulmonary tuberculosis are presented in Tables 4 and 5, respectively.Table 3
Summary of case study results.
KnotsOscillationλGCV111.2500.23920.0300.25030.1000.321210.0650.25920.0250.26130.9750.266310.0250.25920.0250.26131.1500.266Table 4
Knot point location for each subject.
SubjectKnot point locationMinimal lesion6.184Mod advance12.449Far advance18.816KP Miller57.714Table 5
The results of parameter estimation.
SubjectParameterEstimationMinimal lesionα01-0.481α110.105β111-0.119d110.004c011-0.849c1115.6×10−6Mod advanceα020.441α12-0.038β1120.041d12-0.069c0124.642c1121.1×10−7Far advanceα03-0.073α130.004β113-0.004d130.008c013-1.903c1131.4×10−6KP Millerα04-0.407α140.073β114-0.009d140.058c014-0.513c1141.6×10−6Based on the optimal knot points in Table3 and the parameter estimation for each subject in Table 4, the nonparametric regression model based on a mixed estimator for longitudinal data can be written as follows:
(1)
Model estimation for subjectminimal lesion:(62)y^1t=−0.481+0.105x11t−0.119x11t−6.184++0.004z11t−0.849+5.6×10−6cosz11t=−1.33+0.105x11t−0.119x11t−6.184++0.004z11t+5.6×10−6cosz11t(2)
Model estimation for subjectmod advance:(63)y^2t=0.441−0.038x12t+0.041x12t−12.449+−0.069z12t+4.642+1.1×10−7cosz12t=5.083−0.038x12t+0.041x12t−12.449+−0.069z12t+1.1×10−7cosz12t(3)
Model estimation for subjectfar advance:(64)y^3t=−0.073+0.004x13t−0.004x13t−18.816++0.008z13t−1.903+1.4×10−6cosz13t=−1.976+0.004x13t−0.004x13t−18.816++0.008z13t+1.4×10−6cosz13t(4)
Model estimation for subjectKP Miller:(65)y^4t=−0.407+0.073x14t−0.009x14t−57.714++0.058z14t−0.513+1.6×10−6cosz14t=−0.92+0.073x14t−0.009x14t−57.714++0.058z14t+1.6×10−6cosz14t


## 3.1. Mixed Model of Truncated Spline and Fourier Series with Two-Stage Estimation

Lemmas and theorems are used to obtain a nonparametric regression model for longitudinal data with a mixed estimator. The regression curve component that is approximated by the Fourier series estimator is presented in Lemma1 and the penalty component for the Fourier series function is presented in Lemma 2. Following the PWLS form in Lemma 3, we estimate the Fourier series component by using PWLS in Theorem 4. The regression curve component that is approximated by the truncated spline estimator is presented in Lemma 5, and we estimate the truncated spline component using the WLS method in Theorem 6. The results are summarized as follows:Lemma 1.
Ifgkizkit is approached by the Fourier series function, then the goodness of fit is
(3)N−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tgkizkit2=N−1y∗−Zc′Wy∗−Zc,where N=n×T and W is the nT×nT weighting matrix.Proof.
The regression curvegkizkit is a regression curve of an unknown shape and is contained in continuous space C0,π. The component of gkizkit,k=1,2,⋯,q, in Equation (2) is approximated by the Fourier series function with the trend line as follows:
(4)gkizkit=dkizkit+12c0ki+∑h=1HchkicoshzkitIf the regression curveg in Equation (4) involves only one predictor, then it can be written as follows:
(5)g=g1′g2′⋯gn′′,(6)gi=gkizki1gkizki2⋯gkizkiT′,i=1,2,⋯,n.By using Equation (4), Equation (6) can be written in the form of a matrix as follows:
(7)gi=gkizki1gkizki2⋮gkizkiT=dkizki1+12c0ki+c1kicos1zki1+c2kicos2zki1+⋯+cHkicosHzki1dkizki2+12c0ki+c1kicos1zki2+c2kicos2zki2+⋯+cHkicosHzki2⋮dkizkiT+12c0ki+c1kicos1zkiT+c2kicos2zkiT+⋯+cHkicosHzkiT=zki112cos1zki1cos2zki1⋯cosHzki1zki212cos1zki2cos2zki2⋯cosHzki2⋮⋮⋮⋮⋱⋮zkiT12cos1zkiTcos2zkiT⋯cosHzkiTdkic0kic1kic2ki⋮cHki=Zkicki.The Fourier series function in the nonparametric regression component for longitudinal data with predictorq can be expressed in the following form:
(8)gi=∑k=1qZkicki=Z1ic1i+Z2ic2i+⋯+Zqicqi=Zici.So,(9)g=g1′g2′⋮gn′=Z1c1Z2c2⋮Zncn=Z10⋯00Z2⋯0⋮⋮⋱⋮00⋯Znc1c2⋮cn=Zc.Z is a nT×2+Hnq matrix as follows:
(10)Z=Z10⋯00Z2⋯0⋮⋮⋱⋮00⋯Zn,with
(11)Zi=z1i112cosz1i1cos2z1i1⋯cosHz1i1z1i212cosz1i2cos2z1i2⋯cosHz1i2⋮⋮⋮⋮⋱⋮z1iT12cosz1iTcos2z1iT⋯cosHz1iT⋯⋯⋯⋯zqi112coszqi1cos2zqi1⋯cosHzqi1zqi212coszqi2cos2zqi2⋯cosHzqi2⋮⋮⋮⋮⋱⋮zqiT12coszqiTcos2zqiT⋯cosHzqiT.Whereas,c is a 2+Hnq×1 vector given by
(12)c=c1c2⋮cn,ci=d1ic01ic11ic21i⋯cH1i⋯dqic0qic1qic2qi⋯cHqi′.To estimate the Fourier series component, the nonparametric regression model in Equation (1) can be written as
(13)yit−∑j=1pfjixjit=∑k=1qgkizkit+εit,(14)yit∗=∑k=1qgkizkit+εit,i=1,2,⋯,n,t=1,2,⋯,T.The model in Equation (14) can be written in matrix form:
(15)y∗=g+ε.Then, a goodness of fit for the model is formed as follows:(16)N−1∑i=1n∑t=1Twityit∗−∑k=1qgkizkit2,with N=n×T. If the function gkizkit is approached by a Fourier series function as in Equation (4), then the goodness of fit can be presented in the form
(17)N−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tgkizkit2=N−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tdkizkit+12c0ki+∑h=1Hchkicoshzkit2=N−1y∗−Zc′Wy∗−Zc,with W as a weighting matrix for the regression of longitudinal data.Lemma 2.
If the Fourier series is given, then the penalty component is(18)∑k=1qλk∫0π2πgki′′zk2dzk=∑k=1qλk∫0π2πgki′′zk2dzk=c′Dλc.Proof.
The penalty component∑k=1qλk∫0π2/πgki′′zk2dzk in the PWLS optimization based on Equation (4) can be obtained as follows:
(19)gki′′zk=ddzkddzkdkizkit+12c0ki+∑h=1Hchkicoshzkit=−∑h=1Hh2chkicoshzkit.As a result,(20)Pkgk=∫0π2π∑h=1Hh2chkicoshzkit2dzk=2π∫0π∑h=1Hh2chkicoshzkit2+2∑h<mHh2chkicoshzkitm2cmkicosmzkitdzk.To simplify, we defined(21)A=2π∑h=1H∫0πh2chkicoshzkit2dzk,B=2π2∑h<lH∫0πh2chkicoshzkitl2clkicoslzkitdzk.The value ofA will be obtained as follows.
(22)A=2π∑h=1H∫0πh2chkicoshzkit2dzk=1π∑h=1Hh4chki2zkit+2hsinhzkit0π=∑h=1Hh4chki2.Furthermore, the value of B is given by
(23)B=2π2∑h<lH∫0πh2chkicoshzkitl2clkicoslzkitdzk=4π∑h<lHhl2chkiclki∫0πcoshzkitcoslzkitdzk=0.Based on Equations (22) and (23), it can be written as follows:
(24)Pkgk=A+B=∑h=1Hh4chki2.Fori=1,2,⋯,n, we obtained
(25)∑k=1qλkPkgk=∑i=1nci′Dλci=c1′Dλc1+c2′Dλc2+⋯+cn′Dλcn=c′Dλc,where
(26)Dλ=Dλ10⋯00Dλ2⋯0⋮⋮⋱⋮00⋯Dλq,Dλi=00λi14λi24⋯λiH400λi14λi24⋯λiH4⋮⋮⋮⋮⋱⋮00λi14λi24⋯λiH4..Thus, the penalty component can be expressed in a matrix form as follows:(27)∑k=1qλk∫0π2πgki′′zk2dzk=c′Dλc.Lemma 3.
If the goodness of fit component is presented in Lemma1 and the penalty component is given by Lemma 2, then the PWLS is
(28)N−1y∗−Zc′Wy∗−Zc+c′Dλc.In general, PWLS is defined as follows:(29)N−1∑i=1n∑t=1Twityit∗−∑k=1qgkizkit2+∑k=1qλk∫0π2πgki′′zk2dzk,0<λk<∞.Besides, PWLS can be presented in the form of a matrix as follows:(30)N−1y∗−Zc′Wy∗−Zc+c′Dλc.The next step is to obtain a Fourier series estimator in nonparametric regression for longitudinal data derived in Theorem4.Theorem 4.
If paired data are given, which follows the nonparametric regression model for longitudinal data, then the mixed estimator that minimizes PWLS in Lemma3,
(31)mingk∈C0,πN−1∑i=1n∑t=1Twityit∗−∑k=1q∑t=1Tgkizkit2+∑k=1qλk∫0π2πgk′′zk2dzk,is g^k,h,λx,z=Ly∗ with y∗=y−f and L=ZZ′WZ+NDλ−1Z′W.Proof.
The first estimation step in the mixed estimator of the truncated spline and Fourier series in the nonparametric regression model for longitudinal data is performed by estimating the form of the Fourier series estimator by using the PWLS method. The PWLS in Equation (29) can be written in the form of a matrix as follows:
(32)Qc=N−1y∗−Zc′Wy∗−Zc+c′Dλc=N−1y∗′Wy∗−2N−1c′Z′Wy∗+N−1c′Z′WZc+c′Dλc.Next, we complete PWLS optimization using the following steps:(33)minc∈R2+HnqQc=minc∈R2+HnqN−1y∗′Wy∗−2N−1c′Z′Wy∗+N−1c′Z′WZc+c′Dλc.To complete the optimization, the estimators are obtained by performing a partial derivative ofQc concerning c and the results are equaled to zero. The given results are
(34)c^=Z′WZ+NDλ−1Z′Wy∗.By substitutingc^ into Equation (9), we get
(35)g^k,h,λx,z=Zc^=ZZ′WZ+NDλ−1Z′Wy∗=Ly∗.So, the model in Equation (15) can be written as
(36)y∗=g^k,h,λx,z=Ly∗,with y∗=y−f and L=ZZ′WZ+NDλ−1Z′W.Lemma 5.
Iffjixjit is approached with the truncated spline function, then the WLS is
(37)I−Ly−I−LMγ′WI−Ly−I−LMγ.W is a nT×nT weighting matrix.Proof.
fjixjit,j=1,2,⋯,p, is a truncated spline estimator component. The component regression curve f is a linear truncated spline function defined as follows:
(38)fjixjit=αjixjit+∑u=1sβujixjit−Kuji+,with truncated function
(39)xjit−Kuji+=xjit−Kuji,xjit≥Kuji,0,xjit<Kuji.If the component of the regression curvef involves one predictor, then it can be written as follows:
(40)f=f1′f2′⋯fn′′,(41)fi=∑j=1pfjixji1∑j=1pfjixji2⋯∑j=1pfjixjiT′,i=1,2,⋯,n.By using Equation (41), it can be described in the form of a matrix as follows:
(42)fi=fjixji1fjixji2⋮fjixjiT=xji1xji2⋮xjiTαji+xji1−K1ji+1xji1−K2ji+1⋯xji1−Ksji+1xji2−K1ji+1xji2−K2ji+1⋯xji2−Ksji+1⋮⋮⋱⋮xjiT−K1ji+1xjiT−K2ji+1⋯xjiT−Ksji+1β1jiβ2ji⋮βsji=Xjiαji+Sjiβji.The truncated spline function in the nonparametric regression component for longitudinal data withp predictors can be expressed in the following form:
(43)fi=∑j=1pXjiαji+Sjiβji=Xiαi+Siβi.So,(44)f=f1′f2′⋮fn′=X1α1X2α2⋮Xnαn+S1β1S2β2⋮Snβn=X10⋯00X2⋯0⋮⋮⋱⋮00⋯Xnα1α2⋮αn+S10⋯00S2⋯0⋮⋮⋱⋮00⋯Snβ1β2⋮βn=Xα+Sβ.Furthermore, it can be written as follows:(45)f=X10⋯00X2⋯0⋮⋮⋱⋮00⋯XnS10⋯00S2⋯0⋮⋮⋱⋮00⋯Snα1α2⋮αn¯β1β2⋮βn=XSα¯β=Mγ.M is a nT×1+snp matrix and γ is a 1+snp×1 vector.The mixed model of nonparametric regression for longitudinal data in Equation (1) can be written in a matrix as follows:
(46)y=f+g+ε.By substituting Equation (36) with Equation (46), we get
(47)y=f+Ly∗+ε.The truncated spline component can be written as(48)y−Ly∗=f+ε,y−Ly−f=f+ε,y−Ly+Lf=f+ε,y−Ly=f−Lf+ε,I−Ly=I−Lf+ε.Substituting Equation (45) with Equation (48), we obtain
(49)I−Ly=I−LMγ+ε.If the functionfjixjit is approximated by the truncated spline function as in Equation (38), then
(50)ε=I−Ly−I−LMγ.Therefore, we obtain the WLS by(51)ε′ε=I−Ly−I−LMγ′WI−Ly−I−LMγ,where W is a weighting matrix for the regression of longitudinal data. Next, the truncated spline estimator in the nonparametric regression for longitudinal data is derived in Theorem 6.Theorem 6.
If paired data which follows the nonparametric regression model for longitudinal data is given, then the mixed estimator that minimizes WLS in Lemma5,
(52)minγ∈R1+snpI−Ly−I−LMγ′WI−Ly−I−LMγ,is f^k,h,λx,z=MJ−1Ky with J=2−L′WLM−WM and K=L′−IWI−L..Proof.
The second estimation stage of the mixed estimator of the truncated spline and Fourier series in the nonparametric regression model for longitudinal data is performed using the WLS method. The estimator can be obtained by completing the WLS optimization as follows:(53)minγ∈R1+snpI−Ly−I−LMγ′WI−Ly−I−LMγ=minγ∈R1+snpQγ.The estimators are obtained by performing a partial derivative ofQγ and the results are equaled to zero. The partial derivative results are as follows:
(54)γ^=2−L′WLM−WM−1L′−IWI−Ly=J−1Ky,where J=2−L′WLM−WM and K=L′−IWI−L.By substitutingγ^ into the form of a truncated spline estimator component as in Equation (45), we obtain
(55)f^k,h,λx,z=Mγ^=MJ−1Ky,(56)f^k,h,λx,z=Ak,h,λy.After obtainingγ^, we obtain c^ by substituting Equation (56) with Equation (34).
(57)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,c^=Z′WZ+NDλ−1Z′WI−MJ−1Ky.So,(58)g^k,h,λx,z=Zc^=ZZ′WZ+NDλ−1Z′WI−MJ−1Ky=LI−MJ−1Ky,g^k,h,λx,z=Bk,h,λy.By substitutingγ^ and c^ into the mixed estimator of the truncated spline and the Fourier series in the nonparametric regression for longitudinal data, the following estimation results are obtained.
(59)μ^k,h,λ=f^k,h,λ+g^k,h,λ=Mγ^+Zc^=MJ−1Ky+LI−MJ−1Ky=MJ−1K+LI−MJ−1Ky,μ^k,h,λ=Ck,h,λy,with L=ZZ′WZ+NDλ−1Z′W, J=2−L′WLM−WM, and K=L′−IWI−L.The mixed estimator depends on the optimum knot point, oscillation parameter, and smoothing parameter. To obtain the best model, it is essential to select the optimum parameter. One of the criteria to select the optimum parameter is the generalized cross-validation (GCV) method [11]. The GCV function of the nonparametric regression model for longitudinal data is as follows:
(60)GCVKopt,Hopt,λopt=N−1I−Ck,h,λ2N−1traceI−Ck,h,λ2.The optimum knot point, oscillation parameter, and smoothing parameter are obtained by solving the minimum optimization, as presented in Equation (60).


## 3.2. Simulation Study

To demonstrate the performance of the proposed method, we created one sample sizen=10 with t=13. For the simulation study, we considered ten models for each subject. The models are generated from the formula that contains two different functions to represent the truncated spline and Fourier series pattern. A polynomial function is used to present the truncated spline, while a trigonometry function is used to present the Fourier series. Additionally, xit and zit are generated from U0,1 distribution, and random errors εit are generated from a multivariate normal distribution.Using ten subjects and two predictors, the formula for generated data is stated as follows:(61)yit=fitxit+gitzit+εit,i=1,2,⋯,10,t=1,2,⋯,13.The simulation study is applied based on these models, as shown in Table1.Table 1
Model for the simulation study.
SubjectModel1y1t=8.25x1t−11−x1t2+sin2πz1t+cos2πz1t+ε1t2y2t=9.55x2t−11−x2t2+sin2πz2t+cos2πz2t+ε2t3y3t=6.75x3t−11−x3t2+sin2πz3t+cos2πz3t+ε3t4y4t=3.25x4t−11−x4t2+sin2πz4t+cos2πz4t+ε4t5y5t=5.55x5t−11−x5t2+sin2πz5t+cos2πz5t+ε5t6y6t=7.75x6t−11−x6t2+sin2πz6t+cos2πz6t+ε6t7y7t=4.25x7t−11−x7t2+sin2πz7t+cos2πz7t+ε7t8y8t=8.55x8t−11−x8t2+sin2πz8t+cos2πz8t+ε8t9y9t=9.75x9t−11−x9t2+sin2πz9t+cos2πz9t+ε9t10y10t=4.75x10t−11−x10t2+sin2πz10t+cos2πz10t+ε10tFigure1 illustrates the partial relationship between the response and each predictor variable. It can be seen that the relationship between predictor x and the response for each subject tends to change at certain subintervals, which is suitable for the truncated spline estimator. The relationship between z and the response for each subject has a repetitive pattern with a particular trend line, which is suitable for the Fourier series estimator.Figure 1
Scatterplot of simulated data for each subject and predictor: (a) truncated spline; (b) Fourier series.
(a)(b)Wu and Zhang [7] stated that a regression’s performance strongly depends on good knot locations and a good choice of the number of knots. In general, the number of knots is smaller than the sample size n. Considering the scatterplot of simulated data and computational convenience, our study uses three knots K=1,2,3 and three oscillation parameters H=1,2,3. To choose the optimum parameter, we use the minimum GCV criteria. Table 2 provides a summary of the GCV for varying knots and oscillations. What is remarkable is that using one knot and one oscillation with λ=0.033, we obtain the best model with the lowest GCV, 3.112. This model yields satisfactory results with a root mean square error (RMSE) of 0.837.Table 2
Summary of simulation results.
KnotsOscillationλGCV110.0333.11221.7153.14230.9673.116210.0333.19721.8003.20930.9503.158310.0013.45220.1003.47130.0013.441


## 3.3. Case Study

After conducting the simulation, we applied the proposed model to the case to confirm the results of the previous simulation. The data for this research was obtained from a study conducted by Fernandes and Solimun [16], that is, patients with pulmonary tuberculosis disease. Pulmonary tuberculosis is a contagious disease caused by Mycobacterium tuberculosis, which can attack various organs, particularly the lungs. This disease is typical among women in their productive years (ages 15-50 years). The World Health Organization (WHO) declared tuberculosis a global emergency in 1992 [17]. WHO report in 2013 stated that there were 8.6 million tuberculosis cases in 2012, of which 40% of the cases were in Southeast Asia. In a further report, Indonesia was noted as the country with the second-largest number of cases, 2.8 million, in 2015.This study’s dataset consists of four patientsn=4 that represent radiological images of the thorax (stadium), which are minimal lesion, mod advance, far advance, and KP Miller. suPAR level as a response y with mL units in several observation periods every two weeks for six months of treatment t=13 is observed. The predictor variables are the erythrocyte sedimentation rate x with mm/hour units and body mass index z with kg/m3 units.The partial relationship between the suPAR level and each predictor variable for each subject is presented in Figure2. There were changes in data patterns in the four subjects observed for six months with measurements taken every two weeks. The plot in Figure 2 shows a different pattern for each predictor. For this reason, we propose a nonparametric regression approach based on a mixed estimator for longitudinal data. The erythrocyte sedimentation rate x will be approached by a truncated spline estimator, while the body mass index z will be approached by a Fourier series estimator.Figure 2
Scatterplot between suPAR levels with each predictor for the subject: (a)minimal lesion; (b) mod advance; (c) far advance; (d) KP Miller.
(a)(b)(c)(d)In this case study, similar to the simulation study, only three knot points and three oscillation parameters were used. From the various knots and oscillation parameter results, we obtained the GCV values listed in Table3. Interestingly, the data in this table is that the minimum GCV achieved one oscillation parameter under the same conditions as the simulation study, that is, one knot point. However, it has a different smoothing parameter, λ=1.25. This model provides a GCV value of 0.239 with a RMSE of 0.197. The knot point location and the results of the parameter estimation for each subject of patients with pulmonary tuberculosis are presented in Tables 4 and 5, respectively.Table 3
Summary of case study results.
KnotsOscillationλGCV111.2500.23920.0300.25030.1000.321210.0650.25920.0250.26130.9750.266310.0250.25920.0250.26131.1500.266Table 4
Knot point location for each subject.
SubjectKnot point locationMinimal lesion6.184Mod advance12.449Far advance18.816KP Miller57.714Table 5
The results of parameter estimation.
SubjectParameterEstimationMinimal lesionα01-0.481α110.105β111-0.119d110.004c011-0.849c1115.6×10−6Mod advanceα020.441α12-0.038β1120.041d12-0.069c0124.642c1121.1×10−7Far advanceα03-0.073α130.004β113-0.004d130.008c013-1.903c1131.4×10−6KP Millerα04-0.407α140.073β114-0.009d140.058c014-0.513c1141.6×10−6Based on the optimal knot points in Table3 and the parameter estimation for each subject in Table 4, the nonparametric regression model based on a mixed estimator for longitudinal data can be written as follows:
(1)
Model estimation for subjectminimal lesion:(62)y^1t=−0.481+0.105x11t−0.119x11t−6.184++0.004z11t−0.849+5.6×10−6cosz11t=−1.33+0.105x11t−0.119x11t−6.184++0.004z11t+5.6×10−6cosz11t(2)
Model estimation for subjectmod advance:(63)y^2t=0.441−0.038x12t+0.041x12t−12.449+−0.069z12t+4.642+1.1×10−7cosz12t=5.083−0.038x12t+0.041x12t−12.449+−0.069z12t+1.1×10−7cosz12t(3)
Model estimation for subjectfar advance:(64)y^3t=−0.073+0.004x13t−0.004x13t−18.816++0.008z13t−1.903+1.4×10−6cosz13t=−1.976+0.004x13t−0.004x13t−18.816++0.008z13t+1.4×10−6cosz13t(4)
Model estimation for subjectKP Miller:(65)y^4t=−0.407+0.073x14t−0.009x14t−57.714++0.058z14t−0.513+1.6×10−6cosz14t=−0.92+0.073x14t−0.009x14t−57.714++0.058z14t+1.6×10−6cosz14t


## 4. Conclusions

Based on the simulation study and the case study, we selected the best model by using the minimum GCV. The higher knot point or oscillation parameter does not produce a high GCV and vice versa. Therefore, we tried several combinations of knot points and the oscillation parameter to choose the best model. The result of the case study of patients suffering from pulmonary tuberculosis is similar to the simulation study. This study found that the best model uses a one knot point and one oscillation with differentλ. It can be concluded that the simulation study supports the results of the case study.A limitation of this study is that it does not investigate other sample sizes. Consequently, we cannot compare the performance of the developed theory for different sample sizes. Despite its limitations, the study certainly adds to our understanding of the mixed estimator’s new theory in nonparametric regression for longitudinal data.



---

*Source: 4710745-2020-12-10.xml*