Functional and Parametric Estimation in a Semi- and Nonparametric Model with Application to Mass-Spectrometry Data

1 Introduction

We are concerned with the following semi- and nonparametric regression model

where yit is the observed response from i-th individual (i=1,…,n) at time t for (t=1,…,T), xit is the corresponding explanatory variable, αi and βi are individual-specific location and scale parameters and m(⋅) is a baseline intensity function. Here, E(εit)=0, Var(εit)=1, and εit and xit are independent. Of interest here is the simultaneous estimation of αi, βi and m(⋅). We shall assume throughout the paper that εit(i=1,…,n; t=1,…,T) are independent and identically distributed (i.i.d.) with an unknown distribution function, though most results only require that the errors be independent with zero mean.

Model (1) is motivated by analyzing the data generated from mass spectrometer (MS), which is a powerful tool for the separation and large-scale detection of proteins present in a complex biological mixture. Figure 1 is an illustration of MS spectra, which can reveal proteomic patterns or features that might be related to specific characteristic of biological samples. They can also be used for prognosis and for monitoring disease progression, evaluating treatment or suggesting intervention. Two popular mass spectrometers are SELDI-TOF (surface enhanced laser desorption/ionization time-of-fight) and MALDI-TOF (matrix assisted laser desorption and ionization time-of-flight). The abundance of the protein fragments from a biological sample (such as serum) and their time of flight through a tunnel under certain electrical pressure can be measured by this procedure. The y-axis of a spectrum is the intensity (relative abundance) of protein/peptide, and the x-axis is the mass-to-charge ratio (m/z value) which can be calculated using time, length of flight, and the voltage applied. It is known that the SELDI intensity measures have errors up to 50% and that the m/z may shift its value by up to 0.1–0.2% [1]. Generally speaking, many pre-processing steps need to be done before the MS data can be analyzed. Some of the most important steps are noise filtering, baseline correction, alignment, normalization, etc. See, e.g., Guilhaus [2], Banks and Petricoin [3], Baggerly et al. [4], Baggerly et al. [5], Diamandis [6], Feng et al. [7]. We refer readers to Roy et al. [8] for an extensive review about the recent advances in mass-spectrometry data analysis. Here, we assume all the pre-processing steps have already been taken.

Figure 1:

Illustration of MS spectra.

In model (1), m(⋅) represents the common shape for all individuals while αi and βi represents the location and scale parameters for the i-th individual, respectively. Because m(⋅) is unspecified, model (1) may be viewed as a semiparametric model. However, it differs from the usual semi-parametric models in that for model (1), both the parametric and nonparametric components are of primary interest, while in a typical semiparametric setting, the nonparametric component is often viewed as a nuisance parameter. Model (1) contains many commonly encountered regression models as special cases. If all the parametric coefficients αi and βi are known, model (1) reduces to the classical nonparametric regression. On the other hand, if the function m(⋅) is known, then it reduces to the classical linear regression model with each subject having its own regression line. For the present case of αi, βi and function m(⋅) being unknown, the parameters are identifiable only up to a common location-scale change. Thus we assume, without loss of generality, that α1=0 and β1=1. It is also clear that for αi, βi and m(⋅) to be consistently estimable, we need to require that both n and T go to ∞.

There is an extensive literature on semiparametric and nonparametric regression. For semiparametric regression, Begun et al. [9] derived semiparametric information bound while Robinson [10] developed a general approach to constructing n-consistent estimation for the parametric component. We refer to Bickel et al. [11] and Ruppert et al. [12] for detailed discussions on the subject. For nonparametric regression, kernel and local polynomial smoothing methods are commonly used [13–16]. In particular, local polynomial smoothing has many attractive properties including the automatic boundary correction. We refer to Fan and Gijbels [17] and Hardle et al. [18] for comprehensive treatment of the subject.

The existing methods for dealing with nonparametric and semiparametric problems are not directly applicable to model (1). This is due to the mixing of the “finite dimensional” parameters and the nonparametric component. A natural way to handle such a situation is to de-link the two aspects of the estimation through a two-step approach. In this paper, we propose an efficient iterative procedure, alternating between estimation of the parametric component and the nonparametric component. We show that the proposed approach leads to consistent estimators for both the finite-dimensional parameter and the nonparametric function. We also establish the asymptotic normality for the parametric estimators, and the convergence rate for the nonparametric estimate that is then used to derive the optimal bandwidth selection.

2 Main results

In this section, we develop a multi-step approach to estimating both the finite-dimensional parameters αi and βi and the nonparametric baseline intensity m(⋅). Our approach is an iterative procedure which alternates between estimation of αi and βi and that of m(⋅). We show that under reasonable conditions, the estimation for the parametric component is consistent and asymptotically normal when the bandwidth selection are done appropriately. The estimation of the nonparametric component can also attain the optimal rate of convergence.

3 Application to mass spectrometry data

We now apply the proposed multi-step method to a SELDI-TOF mass spectrometry data set from a study of 33 (n=33) liver cancer patients conducted at Changzheng Hospital in Shanghai. For each patient, we extract the m/z values in the region 2,000–10,000 Da (T=21,000), which is believed to contain all the useful information. Figure 2 contains the curves of 10 randomly picked patients.

Figure 2:

Curves of 10 observations and the baseline estimate.

There are some noticeable features in the data. All curves appear to be continuous. They peak simultaneously around certain locations; at each location, curves have the same shape but with different heights. All those features are captured well by our model. Since the observed values of m/z for each person may fluctuate, we need to perform registration to make the analysis easier. Here, we use the observations from the first individual and set his/her m/z values as the reference. Then we use the linear interpolation method to compute the intensities of all the other individuals at the reference m/z locations. After that we get the preprocessed data which has the same m/z values for each observation.

We use the cross-validation method described in Section 2.3 to select the optimal bandwidths with K=33, i.e., leave-one-out cross validation. We compute the MSPE at the grid of h=2,4,6,…,40 and h∗=2,4,6,…,40. Table 1 contains a portion of the result with h=30,32,…,40 and h∗=2,4,6,…,20.

As we can see in Table 1, the minimum MSPE occurs at the location of h=34 and h∗=4, which agrees with our theory that h and h∗ should not be chosen with the same rate for the purpose of estimating the nonparametric component.

Finally, we use the selected bandwidths to estimate the location and scale parameters as well as the nonparametric curve for the whole data set. The estimated parameters are reported in Table 2 and the baseline nonparametric curve estimation is shown in the lower part of Figure 2. From Table 2, we can see that each individual has very different regression coefficients, which was also verified by looking at Figure 2. In addition, comparing the estimated curve for the baseline intensities with the real curves of 10 observations, it is clear that the majority of the peaks and shapes are captured by the nonparametric estimate with appropriate degree of smoothing. A graphical representation of the raw curve of a single individual (16th subject) is also illustrated in Figure 3 along with estimates derived from m˜(⋅) and mˆ(⋅). We can see from the figure that the estimate from mˆ(⋅) is notably better than that from m˜(⋅), which shows that multi-step procedure is effective in improving the estimates for the baseline curve. We observed similar phenomenon for all the other subjects.

Figure 3:

Nonparametric estimates of the curve from m˜(⋅) and mˆ(⋅) for the 16th object.

4 Simulation studies

We conduct simulations to assess the performance of the proposed method for parameter and curve estimation. The true curve m(⋅) is chosen from a moving average smoother of the cross-sectional mean of a fraction of real Mass Spectrometry data in Section 3 after log transformation. We set 10000 m/z values equally-spaced from 1 to 10,000 (T=10,000) and the number of individuals n=30. The true values of the parameters αi,βi,i=1,2,…,n for each individual are shown in Table 3. And the error terms εit are sampled independently from N(0,σ2), where we set σ=0.25.

We apply our multi-step procedure to the simulated data with different choices of the bandwidth. The number of runs is 100. The estimated parameters αˆi and βˆi are shown in Table 3 along with the standard errors. We set h=35, which leads to the smallest SSE of m˜ shown in Table 4. It is evident that the estimation is very accurate for all the location and scale parameters. In addition, for each h and h∗, the computation time is less than two minutes on a laptop.

From Table 4, we can see that the global optimal bandwidths are h=25,h∗=36. It is interesting to note that the optimal bandwidth for m˜(⋅) is h=35, which is different from the optimal bandwidth for the final estimator.

To evaluate the quality of the our multi-step estimation method for the nonparametric baseline function, we consider a classical nonparametric estimation on another set of data where the same true function m(⋅) is used but αi=0,βi=1 for all i=1,…,n. We applied the same local linear estimation with different bandwidths. The mean SSE of the estimated m(⋅) from 100 repetitions for different hs are given in Table 5. When we applied the multi-step estimation procedure, the best mean SSE we achieved in Table 4 is very close to the minimal mean SSE 0.4442 for the oracle estimator. This comparison confirms that there is little loss of information in the proposed method when both parametric and nonparametric components are estimated simultaneously.

We use cross-validation to get a data-driven choice of the bandwidths. Here, we set K=5 to get a mean MSPE of every different bandwidth choices of both steps over 100 runs, and the optimal bandwidths are those with the minimum mean MSPE. The mean MSPE values are shown in Table 6, from which we can see that the smallest value is located at h=25,h∗=36, which is quite close to the optimal bandwidths h=25 and h∗=38 in Table 4. Therefore, the cross-validation idea appears to work well in terms of selecting the best bandwidths.

5 Discussion

This paper proposes a semi- and nonparametric model suitable for analyzing the mass spectra data. The model is flexible and intuitive, capturing the main feature in the MS data. Both the parametric and nonparametric components have natural interpretation. A multi-step iterative algorithm is proposed for estimating both the individual location and scale regression coefficients and the nonparametric function. The algorithm combines local linear fitting and the least squares method, both of which are easy to implement and computationally efficient. Both simulation studies and real data analysis demonstrate that the proposed multi-step procedure works well.

The local linear fitting for the nonparametric function estimation maybe replaced with other nonparametric estimation techniques. Because the location and scale parameters are subject specific, the empirical Bayes method [19] may be used. In addition, nonparametric Bayes may also be applicable with the nonparametric function being modeling as a realization of Gaussian process.

The proposed model and the associated iterative estimation method do not account for the random error in the measurement of X. It is desirable to incorporate the measurement error into the model [20].

Many studies involving MS data are aimed at classifying patients of different disease types. The information of peaks are usually applied as the basis of the classifier. The proposed method provides a natural way to identify the peaks for different group of patients by using the multi-step estimation procedure on each group and finding out the corresponding nonparametric baseline function. From the estimated baseline function, the information of peaks can be easily extracted, which can then be used for classification.