Home | About | Journals | Submit | Contact Us | Français |

**|**Comput Math Methods Med**|**v.2013; 2013**|**PMC3610364

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Data Set
- 3. Feature Description
- 4. Classification
- 5. Classification Results
- 6. Discussion
- 7. Conclusion
- References

Authors

Related links

Comput Math Methods Med. 2013; 2013: 904267.

Published online 2013 March 12. doi: 10.1155/2013/904267

PMCID: PMC3610364

Suxian Cai,^{
1
} Shanshan Yang,^{
1
} Fang Zheng,^{
1
} Meng Lu,^{
1
} Yunfeng Wu,^{
1
,}^{*} and Sridhar Krishnan^{
2
}

*Yunfeng Wu: Email: gro.eeei@uw.y

Academic Editor: Thierry Busso

Received 2012 October 25; Revised 2013 January 31; Accepted 2013 February 11.

Copyright © 2013 Suxian Cai et al.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Analysis of knee joint vibration (VAG) signals can provide quantitative indices for detection of knee joint pathology at an early stage. In addition to the statistical features developed in the related previous studies, we extracted two separable features, that is, the number of atoms derived from the wavelet matching pursuit decomposition and the number of significant signal turns detected with the fixed threshold in the time domain. To perform a better classification over the data set of 89 VAG signals, we applied a novel classifier fusion system based on the dynamic weighted fusion (DWF) method to ameliorate the classification performance. For comparison, a single leastsquares support vector machine (LS-SVM) and the Bagging ensemble were used for the classification task as well. The results in terms of overall accuracy in percentage and area under the receiver operating characteristic curve obtained with the DWF-based classifier fusion method reached 88.76% and 0.9515, respectively, which demonstrated the effectiveness and superiority of the DWF method with two distinct features for the VAG signal analysis.

The knee is the largest and most complex joint in the human body [1]. A pair of knees support nearly the entire weight of the human body and help the body perform different locomotion functions. Knee osteoarthritis is a common form of rheumatic disorder caused by the degeneration or damage of articular cartilage in the knee joint. Detection of knee joint pathology at an early stage can help clinicians apply appropriate therapeutical or surgical procedures to retard the degenerative process in the affected knee joint [2, 3].

Arthroscopy is usually performed as a semi-invasive surgical procedure for knee joint disorder detection [1]. Physicians can inspect the interior of a joint with an arthroscope fiber that is inserted into the knee through a small incision. Although arthroscopy is often used as the gold standard for relatively low-risk assessment of joint surfaces [4], it cannot be applied to patients whose knees are in a highly degenerated state due to arthritis, ligamentous instability, meniscectomy, or patellectomy. Furthermore, arthroscopy is not suited for routine examinations of the articular cartilage because the same incision may not undergo repeated arthroscope fiber invasions due to bacterial infection.

Magnetic resonance imaging can assist in the characterization of *in vivo* orthopaedic condition of articular cartilage and is also a popular noninvasive method for the assessment of knee joint degeneration [5]. The weakness of the magnetic resonance imaging is that such a technique is only able to display anatomical morphology [6]. The imaging protocol cannot support functional condition detection of the knee joint during leg movement, because the subject has to lay down throughout the magnetic resonance scanning procedure.

Vibration arthrometry is an alternative technology for noninvasive detection of knee pathologies [4, 7]. The knee joint vibration arthrographic (VAG) signal can be recorded by accelerometer or electrostethoscope sensors attached on the surface of the knee cap [8–10]. For the healthy adults, their articular surfaces in the knee are smooth and slippery without any cartilage friction or collision. However, the vibrations generated due to the friction between the degenerative articular cartilages are expected to present anomalous patterns in the amplitude and frequency scales [2]. The vibration arthrometry not only provides a clinical option of the noninvasive and low-cost routine detection for knee joint disorders but also supports the functional study on the knee joint during leg movements [1].

Computer-aided analysis of knee joint VAG signals is very useful for screening and monitoring of articular cartilage disorders at an early stage [11–13]. Based on the noninvasive detection results, the computational algorithms may effectively help the medical experts make an accurate decision, so that the frequency of the diagnostic open surgery with arthroscope can be reduced [8, 14–16]. Adaptive filtering techniques based on the least-mean-square (LMS) and recursive least-squares lattice (RLSL) algorithms were used to remove muscle contraction interference present in VAG signals [17]. Tavathia et al. [18] and Moussavi et al. [19] proposed different linear prediction models and adaptive segmentation methods for parameterization and screening of VAG signals. Jiang et al. [20] extended the application of vibration arthrometry to artificial knee joints *in vitro* and analyzed the VAG signal with the root-mean-squared (RMS) value and the parameters of an autoregressive (AR) model. Matching pursuit (MP) time-frequency distributions (TFDs) of VAG signals were studied by Krishnan et al. [9], and a modified local discriminant bases algorithm was proposed by Umapathy and Krishnan [14] to distinguish the abnormal VAG signals from the normal ones. In order to simplify the procedures of signal processing and decision making, Rangayyan and Wu [11, 12] proposed the statistical parameters including form factors, skewness, kurtosis, probability density function entropy, variance of mean-squared values, and turns count with adaptive threshold, for the screening of VAG signals based on the radial-basis function network (RBFN). Mu et al. [21] used the linear and nonlinear strict 2-surface proximal classifiers to test a subset of the aforementioned statistical parameters. In this paper, the number of atoms derived from the wavelet MP decomposition and the turns count detected with the fixed threshold in the waveform variability analysis are extracted as features, and a classifier fusion system based on the dynamic weighted fusion (DWF) method is proposed for the classification of the VAG signals.

The VAG data were collected by the research group of Rangayyan, University of Calgary, Canada [9]. The experimental protocol was approved by the Conjoint Health Research Ethics Board of the University of Calgary. Each subject was requested to sit on a rigid table in a relaxed position with the leg being tested freely and suspended in midair. The knee joint vibration was measured by placing a miniature accelerometer at the middle position of the patella [8]. Minute electrical charges were generated by the accelerometer based on the acceleration and deceleration of the knee movement when the subject swung the leg over an approximate angle range of 135° to 0° and back to 135° in the duration *T* = 4s. Each VAG signal was conditioned by an isolation preamplifiers to prevent the aliasing effects. The signal was then amplified and digitized with a data acquisition board and the National Instruments LabVIEW software at the sampling rate *f*_{s} = 2kHz with 12-bit resolution per sample. Auscultation of the knee joint using a stethoscope was also performed, and a qualitative description of the sound intensity and type was recorded.

In the present study, we used a total of 89 VAG signals (recorded from 51 healthy volunteers and 38 subjects with knee joint pathologies), the same as investigated in a few previous related studies [11, 12, 14]. The normal signals were recorded from the healthy subjects identified by physical examinations. The abnormal signals, an example of which is shown in Figure 1(b), were collected from symptomatic patients scheduled to undergo arthroscopic examinations independent of the VAG studies. The knee joint disorders were associated with chondromalacia of different grades, meniscal tear, tibial chondromalacia, and anterior cruciate ligament injuries, which were confirmed in arthroscopy. Compared with the normal VAG signal displayed in Figure 1(a), we can observe that the abnormal signal exhibits a higher degree of variability in the time domain, as illustrated in Figure 1(b).

VAG signals are nonstationary in nature, that is, such signals exhibit time-varying spectral characteristics, and they cannot be accurately represented by common signal processing techniques such as Fourier transform and AR modeling [9]. It is therefore better to use a joint time-frequency approach for VAG analysis. For time-frequency representation of a signal, the methodology of MP introduced by Mallat and Zhang [22] is able to decompose the signal using basis functions with good time-frequency properties which are referred to as atoms.

The MP method is a so-called “greedy” algorithm that successively approximates a signal *x*(*t*) of *N* samples with orthogonal projections onto elements from a waveform dictionary = {*d*_{r}(*t*)}_{rΓ} of *P* vectors, in which $||{d}_{r}||=\sqrt{[\int {d}_{r}^{\mathrm{2}}(t)dt]}=\mathrm{1}$. Gabor function, local cosine trees, and wavelet packets are often applied to build up dictionaries for MP applications. In this investigation, we implemented the Daubechies wavelet MP decomposition [23], because the Daubechies wavelets are a family of orthogonal wavelets that have a support of minimum size for a given number of vanishing moments [24], and such wavelets can be used for signal decomposition with excellent time and scale properties [25]. The projection of VAG signal *x*(*t*) using the dictionary of wavelet packet bases, *d*_{rm}(*t*), calculated with a Daubechies 8 (db8) filter can be formulated as

$$\begin{array}{c}x\left(t\right)=\sum _{m=\mathrm{0}}^{M-\mathrm{1}}{a}_{m}{d}_{{r}_{m}}\left(t\right),\end{array}$$

(1)

where *a*_{m} are the expansion coefficients and *M* denotes the iterations of decomposition. And the wavelet MP decomposition can be implemented as follows. In the beginning, *x*(*t*) is projected on a vector *d*_{r0}(*t*) , and the residue *R*^{1}*x*(*t*) is computed, that is,

$$\begin{array}{c}x\left(t\right)=x,{d}_{{r}_{\mathrm{0}}}{d}_{{r}_{\mathrm{0}}}\left(t\right)+{R}^{\mathrm{1}}x\left(t\right),\end{array}$$

(2)

where *x*, *d*_{r0} denotes the inner product (projection). Since the first atom *d*_{r0}(*t*) is orthogonal to *R*^{1}*x*(*t*), we have

$$\begin{array}{c}{||x||}^{\mathrm{2}}={|x,{d}_{{r}_{\mathrm{0}}}|\mathrm{2}+{||{R}^{\mathrm{1}}x||}^{\mathrm{2}}.}^{}\end{array}$$

(3)

In order to minimize ||*R*^{1}*x*||, *r*_{0} Γ is chosen such that |*x*, *d*_{r0}| is maximum, that is,

$$\begin{array}{c}|x,{d}_{{r}_{\mathrm{0}}}|\ge \underset{r\mathrm{\Gamma}}{\mathrm{sup}}\end{array}$$

(4)

The MP iterates this procedure by subdecomposing the residue. And the VAG signal *x*(*t*) after *M* iterations of decomposition is then expressed as

$$\begin{array}{c}x\left(t\right)=\sum _{m=\mathrm{1}}^{M-\mathrm{1}}{R}^{m}x,{d}_{{r}_{m}}{d}_{{r}_{m}}\left(t\right)+{R}^{M}x\left(t\right),\end{array}$$

(5)

where |*x*, *d*_{rm}| ≥ sup _{rΓ} |*x*, *d*_{r}| and *R*^{0}*x*(*t*) *x*(*t*). Since the residue term *R*^{M}*x*(*t*) can be regarded as noise after sufficient iterations, one common approach to stop the iterative process depends on the convergence of residual energy, ||*R*^{m}*x*||^{2}. In the previous related studies [9, 26], Krishnan et al. used the decay parameter defined by Mallat and Zhang [22] to end the MP decomposition based on the Gabor function dictionary. Although the decay parameter as an iterative indicator is well devised in theory, for a wavelet packet dictionary that contains *P* = *N*log _{2}*N* vectors, each MP iteration then requires *O*(*N*log _{2}*N*) operations. For the VAG data analyzed, each signal consists of *N* = *f*_{s} × *T* = 8000 samples, such that the wavelet MP decomposition is computationally expensive.

In the present study, we propose using a signal-to-noise ratio (SNR) as the alternative indicator to determine the iterations of the wavelet MP decomposition. According to the VAG signals in Figure 1, it can be inferred that, with a given SNR, the number of the wavelet MP decomposition iterations for a normal signal will be fewer than that for an abnormal signal, because the abnormal signal is much more noisy and also contaminated by a larger amount of artifacts such as muscle contraction interference. Thus, the number of MP iterations can be considered as a potential feature for classification applications. On the other hand, the MP decomposition with many great iterations can provide an excellent value of SNR and is also suited for denoising applications [26], but such an implementation is time consuming as mentioned above. It is therefore necessary to search for an appropriate value of SNR that makes a tradeoff between efficiency and effectiveness of the wavelet MP decomposition. After testing the SNR with different values, we found that the SNR of 15dB could be an excellent indicator to determine the wavelet MP iterations. The *P* value (0.0002) of the number of atoms (Natom) obtained with the Student's *t*-test [27] indicated that the Natom values (numerically equal to the number of decomposition iterations *M*) were significantly different between the normal and abnormal VAG signals. Such results also confirm our assumption mentioned above.

From Figure 2, we can observe that the much noise has been reduced in the normal and abnormal VAG signals, reconstructed with 145 and 784 MP atoms, respectively. The results of noise removal in the VAG signals, as depicted in Figure 2, are even better than those in the previous study using the Gabor function dictionary and the energy decay parameter [26].

Besides the time-frequency MP decomposition, the waveform variability analysis in the time domain may be useful for classification as well. According to Rangayyan [28], a signal sample can be identified as a “turn” if it satisfies the following two conditions at the same time: (1) it represents a change in direction in the signal, that is, an alteration in the sign of the derivative (from positive to negative or vice versa); (2) the difference between its amplitude and that of the preceding sample is over a certain threshold. Willison [29] used the turns count method to analyze the electromyographic signal. The experiments showed that the electromyographic signal recorded from a patient with myopathy usually possesses more turns than the signal of a healthy subject at a comparable level of volitional effort. In the present work, we first normalized each VAG signal in the amplitude range from zero to unity, the same as in the recent studies [11, 12]. In each signal, the amplitude of all samples was amplified with the same scale, so that the variability information of the signal can be preserved. Before applying the turns count method in the signal, we implemented a filtering procedure using a 10th-order lowpass Butterworth filter (cutoff frequency: 50Hz) with unit gain at direct current (DC) [23]. This lowpass Butterworth filter causes a delay of 100 samples (or 0.05s), which was calibrated after the filtering procedure in our experiments. The reason that we used the lowpass Butterworth filter instead of the signal reconstructed with the MP atoms, as shown in Figure 2, is that the MP method is unable to eliminate the artifacts such as the interference caused by muscle contractions or 50 or 60Hz power-supply lines.

In the past work of Rangayyan and Wu [12], the threshold to determine the significance of a turn was adaptively set to be 0.5*σ*_{v}, where *σ*_{v} denotes the standard deviation of the VAG signal analyzed. Although the turns detected with the adaptive threshold provide good discriminant information for classification ofVAG signals, the number of turns counted from a normal signal is larger than that of an abnormal one, because the standard deviation of the normal VAG signal is usually smaller than that of the abnormal signal [12]. Such a result, however, somewhat deviates from our expectation that the turns associated with an abnormal VAG signal would be larger in number due to a higher degree of variability. In the present study, we fixed the amplitude threshold at 0.2 to compute the turns over the normalized and filtered VAG signals.

Figure 3 shows the results of the turns count with the fixed threshold (TCFT) method for the VAG signals in Figure 1, which were normalized and processed with the Butterworth filter. We can observe that more significant turns have been identified in the abnormal signal, as marked in Figure 3(b), in comparison with the normal signal shown in Figure 3(a). In addition, the *P* value of the TCFT obtained with the Student's *t*-test is 0.0013 (significance level: *P* < 0.01), which indicates a significant difference between the normal and abnormal signals.

In addition to two aforementioned features, we also considered the other five features extracted from the same VAG data set in our previous work [11–13, 30]. These features included the form factors computed for the first half (FF1) and the second half (FF2) of each VAG signal [11]; the variance of the mean-squared value (VMS) of the each signal [12]; the mean value (*μ*) of the Parzen-window probability density function of each signal [13]; the fractal dimension (FD) estimated by the power spectral analysis [30]. Total seven features were combined in the vector form for the following pattern analysis task.

To perform the signal classifications, we applied a single least-squares support vector machine (LS-SVM) and the ensembles of several component LS-SVM classifiers, the details of which are presented as follows.

The support vector machine (SVM) proposed by Cortes and Vapnik [31] is a type of universal approximator, the learning of which follows the structural risk minimization criterion [32]. To optimize the SVM model parameters, a subset of the representative training data is selected to be the support vectors, which are considered to be the most informative for the classification task. By choosing the nonlinear inner-product kernels in the network, the SVM is able to perform the same function as the polynomial learning machine, radial basis function network, or multilayer perceptron with a single hidden layer [33, 34]. The LS-SVM was proposed by Suykens et al. [35] as a reformulation to the standard SVM, with an improvement of the moderate complexity. The learning of the LS-SVM is implemented by minimizing a regularized least-squares cost function with equality constraints, under the Kuhn-Tucker condition [36]. Recently, the LS-SVM has also been widely used in a number of biomedical applications [37–39].

To determine the kernel function suited for the VAG signal classification, we implemented the LS-SVM using the linear, polynomial, sigmoid, and Gaussian kernels, one by one specified by different model parameters, and then evaluated each LS-SVM with the leave-one-out (LOO) method. As a type of cross-validation approach, the LOO partition procedure repeatedly used each signal once for the validation and the remaining signals for training [34]. By checking the accuracy and the optimal separating hyperplane provided by each LS-SVM, we chose the polynomial kernel function, the degree and intercept parameters of which equal to 2 and 1, respectively, and set the regularization parameter of the LS-SVM to be 5.

As an emerging machine learning methodology based on the “divide and conquer” principle [33], ensemble of classifiers was widely used in the literature [40–44], with the aim to achieve a better performance versus a single classifier. By combing a finite number of component neural networks (CNNs) with a well-devised combination rule [41–43] or fusion strategy [44–46], a neural network ensemble is expected to provide an informative overall decision that is supposedly superior to that attained by any one of the CNNs acting solely [33].

The most popular ensemble algorithms are AdaBoost [47] and Bagging [48]. The AdaBoost works by repeatedly training a given type of weak-learning machine from different distributed training data sets and then combining their outputs. The distribution of training data for the current CNN is boosted depending on the performance of previous CNNs, that is, the training data that are incorrectly predicted by previous CNNs will be chosen with priority to train the current CNN. In spite of the effectiveness, the AdaBoost is very sensitive to outliers and sometimes results in overfitting [49]. On the other hand, the Bagging algorithm introduces the bootstrap approach [50] into the training data resampling procedure [48] and aggregates the CNNs with the simple average strategy [33]. In the bootstrap procedure, each data sample was selected separately at random from the original data set such that a particular data sample could appear multiple times in a bootstrap-generated data set. The bias of the Bagging ensemble would converge by averaging, while the variance falls much smaller than that of each CNN.

Since the LS-SVM is not a weak-learning machine [35], we used the Bagging algorithm rather than the AdaBoost for the ensemble of 5 component LS-SVMs (CSVMs) which were labeled from CSVM1 to CSVM5 in numerical sequence. The CSVMs combined with the Bagging algorithm were trained by different bootstrap-generated data sets. The number of signals in each bootstrap-generated data set was of equal size to the original VAG data set. The testing data set for each CSVM was the same as the original VAG data set. Because the training data for each CSVM were generated using the bootstrap approach, it is not necessary to apply the LOO method to the ensemble system any more.

According to the linear combination rule of the Bagging, the CSVMs are simply averaged in the ensemble, so that the effectiveness of the ensemble would be affected by some of the CSVMs with poor performance, because the simple average strategy treats all the CSVMs equally. With the aim to utilize the diverse knowledge generated by the CSVMs, we applied a dynamic weighted fusion (DWF) rule to adaptively combine the CSVMs in the classifier fusion system [51].

Suppose that a total of *K* CSVMs are linearly combined in the classifier fusion system. The local decision generated by the *k*th CSVM is denoted as *g*_{k}(**f**^{n}), with regard to the feature vector of the *n*th VAG signal, **f**^{n}. The classifier fusion system then provides the overall classification decision, *g*_{DWF}(**f**^{n}), by linearly combining the CSVMs with the weights *w*_{k}(**f**^{n}) that are varied from one signal to another. Thus, the DWF-based ensemble output can be formulated as

$$\begin{array}{c}{g}_{\text{DWF}}\left({\mathbf{f}}^{n}\right)=\sum _{k=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right){g}_{k}\left({\mathbf{f}}^{n}\right).\end{array}$$

(6)

The nonnegative and normalization constraints on the fusion weights, as widely accepted in the literature [41, 43, 45], can be written as

$$\begin{array}{c}\sum _{k=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right)=\mathrm{1},{w}_{k}\left({\mathbf{f}}^{n}\right)\ge \mathrm{0}.\end{array}$$

(7)

The task of the DWF is to determine the fusion weights that help the ensemble system provide an overall classification decision with higher accuracy. To achieve this goal, let us study the error term of the CSVMs and the DWF-based ensemble. Concerning the *k*th CSVM, the squared error that characterizes the difference between the local decision and the desired class label, *l*(**f**^{n}), in relation to the *n*th VAG signal is

$$\begin{array}{c}{e}_{k}^{\mathrm{2}}\left({\mathbf{f}}^{n}\right)={\left[l\left({\mathbf{f}}^{n}\right)-{g}_{k}\left({\mathbf{f}}^{n}\right)\right]}^{\mathrm{2}}.\end{array}$$

(8)

Then, the squared error of the ensemble, *e*_{DWF}^{2}(**f**^{n}), is estimated in an analogous manner to that of each CSVM. Consider that the fusion weights are normalized, as presented in (7), the class label *l*(**f**^{n}) can be split by multiplying the fusion weights, so that *e*_{DWF}^{2}(**f**^{n}) is derived as follows:

$$\begin{array}{c}{e}_{\text{DWF}}^{\mathrm{2}}\left({\mathbf{f}}^{n}\right)={\left[l\left({\mathbf{f}}^{n}\right)-\sum _{k=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right){g}_{k}\left({\mathbf{f}}^{n}\right)\right]}^{\mathrm{2}}\\ ={\left[\sum _{k=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right)l\left({\mathbf{f}}^{n}\right)-\sum _{k=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right){g}_{k}\left({\mathbf{f}}^{n}\right)\right]}^{\mathrm{2}}\\ =\left\{\sum _{k=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right)\left[l\left({\mathbf{f}}^{n}\right)-{g}_{k}\left({\mathbf{f}}^{n}\right)\right]\right\}\\ \times \left\{\sum _{j=\mathrm{1}}^{K}{w}_{j}\left({\mathbf{f}}^{n}\right)\left[l\left({\mathbf{f}}^{n}\right)-{g}_{j}\left({\mathbf{f}}^{n}\right)\right]\right\}& =\sum _{k=\mathrm{1}K}^{}\sum _{j=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right){w}_{j}\left({\mathbf{f}}^{n}\right)\left[l\left({\mathbf{f}}^{n}\right)-{g}_{k}\left({\mathbf{f}}^{n}\right)\right]& \times \left[l\left({\mathbf{f}}^{n}\right)-{g}_{j}\left({\mathbf{f}}^{n}\right)\right]=\sum _{k=\mathrm{1}K}^{}\sum _{j=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right){w}_{j}\left({\mathbf{f}}^{n}\right){e}_{k}\left({\mathbf{f}}^{n}\right){e}_{j}\left({\mathbf{f}}^{n}\right),\end{array}$$

(9)

where *e*_{k}(**f**^{n}) = [*l*(**f**^{n}) − *g*_{k}(**f**^{n})] represents the instantaneous error of the *k*th CSVM.

Considering (7) and (9), the minimization of the squared error of the ensemble is equivalent to the constrained quadratic programming (CQP) problem specified as follows:

$$\begin{array}{c}\text{minimize}{e}_{\text{DWF}}^{\mathrm{2}}\left({\mathbf{f}}^{n}\right)=\sum _{k=\mathrm{1}K}^{}\sum _{j=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right){w}_{j}\left({\mathbf{f}}^{n}\right){e}_{k}\left({\mathbf{f}}^{n}\right){e}_{j}\left({\mathbf{f}}^{n}\right),\text{subject to}\sum _{k=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right)=\mathrm{1},{w}_{k}\left({\mathbf{f}}^{n}\right)\ge \mathrm{0}.\end{array}$$

(10)

In order to solve the CQP problem presented in (10), we applied the Lagrange multiplier method [36] and defined the cost function as

$$\begin{array}{c}C({w}_{\mathrm{1}}\left({\mathbf{f}}^{n}\right),\dots ,{w}_{K}\left({\mathbf{f}}^{n}\right),\lambda \left({\mathbf{f}}^{n}\right))\\ =\sum _{k=\mathrm{1}K}^{}\sum _{j=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right){w}_{j}\left({\mathbf{f}}^{n}\right){e}_{k}\left({\mathbf{f}}^{n}\right){e}_{j}\left({\mathbf{f}}^{n}\right)-\lambda \left({\mathbf{f}}^{n}\right)\left[\sum _{k=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right)-\mathrm{1}\right],\end{array}$$

(11)

where the nonnegative coefficient *λ*(**f**^{n}) represents the Lagrange multiplier, which varies from one signal to another.

According to the weak Lagrangian principle [36], the optimum solution to the CQP problem, {**w***(**f**^{n}), *λ**(**f**^{n})}, is the stationary point of the cost function presented in (11) and satisfies the following unique equations:

$$\begin{array}{c}\frac{C({w}_{\mathrm{1}}\left({\mathbf{f}}^{n}\right),\dots ,{w}_{K}\left({\mathbf{f}}^{n}\right),\lambda \left({\mathbf{f}}^{n}\right)){w}_{k}\left({\mathbf{f}}^{n}\right)}{}=\mathrm{2}\sum _{j=\mathrm{1}}^{K}{w}_{j}\left({\mathbf{f}}^{n}\right){e}_{k}\left({\mathbf{f}}^{n}\right){e}_{j}\left({\mathbf{f}}^{n}\right)-\lambda \left({\mathbf{f}}^{n}\right)=\mathrm{0},\frac{C({w}_{\mathrm{1}}\left({\mathbf{f}}^{n}\right),\dots ,{w}_{K}\left({\mathbf{f}}^{n}\right),\lambda \left({\mathbf{f}}^{n}\right))\lambda \left({\mathbf{f}}^{i}\right)}{}=\sum _{k=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right)-\mathrm{1}=\mathrm{0}.\end{array}$$

(12)

The optimal weights *w*_{k}(**f**^{n}) of the DWF that minimize the squared error of the ensemble system can be obtained by solving (12), that is,

$$\begin{array}{c}{w}_{k}^{}\end{array}$$

(13)

Because the error term of the CSVM can be estimated when the *n*th VAG class label *l*(**f**^{n}) is given and the CSVM model parameters are specified, the optimal fusion weights can be directly computed according to (13).

Now let us divert our attention to the DWF-based ensemble error term. Considering (9) and (13), we have

$$\begin{array}{c}{e}_{\text{DWF}}^{\mathrm{2}}\left({\mathbf{f}}^{n}\right)=\sum _{k=\mathrm{1}K}^{}\sum _{j=\mathrm{1}}^{K}{w}_{k}\left({\mathbf{f}}^{n}\right){w}_{j}\left({\mathbf{f}}^{n}\right){e}_{k}\left({\mathbf{f}}^{n}\right){e}_{j}\left({\mathbf{f}}^{n}\right)& =\frac{\mathrm{1}}{{\sum}_{i=\mathrm{1}}^{K}{\sum}_{j=\mathrm{1}}^{K}{e}_{i}^{-\mathrm{1}}\left({\mathbf{f}}^{n}\right){e}_{j}^{-\mathrm{1}}\left({\mathbf{f}}^{n}\right)}.\end{array}$$

(14)

It is clear that both of *e*_{k}^{2}(**f**^{n}) and *e*_{DWF}^{2}(**f**^{n}) are nonnegative, that is, *e*_{k}^{2}(**f**^{n}) ≥ 0 and ∑_{i=1}^{K}∑_{j=1}^{K}*e*_{i}^{−1}(**f**^{n})*e*_{j}^{−1}(**f**^{n}) ≥ 0. To compare with the squared error of the CSVM, we may employ the division operator such that

$$\begin{array}{c}\frac{{e}_{k}^{\mathrm{2}}\left({\mathbf{f}}^{n}\right)}{{e}_{\text{DWF}}^{\mathrm{2}}\left({\mathbf{f}}^{n}\right)}=\mathrm{1}+\sum _{\begin{array}{c}i=\mathrm{1}\\ i\ne k\end{array}K}^{}\sum _{\begin{array}{c}j=\mathrm{1}\\ j\ne k\end{array}}^{K}\frac{{e}_{k}^{\mathrm{2}}\left({\mathbf{f}}^{n}\right)}{{e}_{i}\left({\mathbf{f}}^{n}\right){e}_{j}\left({\mathbf{f}}^{n}\right)}\ge \mathrm{1}.\end{array}$$

(15)

Therefore, the optimal fusion weights derived in (15) guarantee that the DWF-based ensemble system more or less likely outperforms a single CSVM.

For a fair performance comparison, the CSVMs combined by the DWF-based ensemble were with the same bootstrap-generated training data and model parameters as those CSVMs combined by the Bagging algorithm. The testing data for the CSVMs were the entire data set of the VAG signals.

In addition, we also selected some subsets of fixed size 15 to train the CSVMs, and would like to evaluate the generalization capability of the proposed DWF-based ensemble in the case of small-size training data. The subset for each CSVM was actively selected according to the quadratic Renyi entropy maximization criterion [52] as follows. In the experiments, we first randomly partitioned the entire VAG data set into two subsets: the working subset containing 15 VAG signals and the candidate subset containing the remaining 74 signals. In each iteration step, we randomly selected one VAG signal **f**^{c} from the candidate subset to replace an arbitrary signal **f**^{w} in the working subset. If the quadratic Renyi entropy of the new working subset increases, then the working signal **f**^{w} should be replaced by the candidate signal **f**^{c}; otherwise, **f**^{w} remains in the working subset, and the candidate signal **f**^{c} should be rejected and returned to the candidate subset [35]. After a few iteration steps, the quadratic Renyi entropy would become stable when reaching the maximum value, then the active selection of the fixed-size subset is terminated. The subsets for the 5 CSVMs were selected independently, and the corresponding quadratic Renyi entropy values with respect to the iteration steps are showed in Figure 4. After the model parameters of the CSVMs were optimized by the corresponding subsets of the fixed size, the original VAG data were input to the CSVMs and the DWF-based ensemble for testing and comparison purpose.

The classification accurate rate in percentage of the single LS-SVM with the LOO method (LS-SVM/LOO) was 87.64%, which was better than the result (accuracy: 61.8%) obtained by the Fisher's linear discriminant analysis using the same feature set. In addition, the receiver operating characteristic (ROC) curve technique was implemented to test the overall diagnostic performance for all classification methods. According to Table 1, the area under the ROC curve (*A*_{z}) obtained with the LS-SVM/LOO was 0.7523 with a standard error (SE) of 0.0536, as illustrated in Figure 5.

ROC curves provided by (a) the five component LS-SVM classifiers in the ensemble; (b) the LS-SVM evaluated by the leave-one-out method, the Bagging, and the dynamic weighted fusion (DWF) method. The ROC results were obtained with the testing data.

Classification results of different classifiers and ensemble systems on the testing data set. LS-SVM/LOO: the least-squares support vector machine evaluated with the leave-one-out (LOO) method. CSVM: the component least-squares support vector machines. **...**

In the ensemble experiments, the highest accuracy provided the CSVMs was 83.15%, whereas the lowest one was only 76.41%. Despite that the *A*_{z} value of the CSVM4 was slightly smaller than that of the LS-SVM/LOO, any of the other four CSVMs outperformed the LS-SVM/LOO with higher *A*_{z} values under the ROC curves (see Figure 5).

Referring to Table 1 and Figure 5, the superiority of the DWF-based ensemble in the diagnostic performance was prominent. The DWF-based ensemble provided the highest overall accuracy of 88.76%, the best *A*_{z} value of 0.9515, and the lowest SE of 0.0244. It can be observed from Figure 5 that the ROC curve obtained with the DWF-based ensemble was consistently over either of that provided by the Bagging or the single LS-SVM/LOO classifier. In comparison with the single LS-SVM/LOO classifier, the Bagging ensemble can improve the diagnostic ROC curve, but its accurate rate was even 8.99% lower than that of the LS-SVM/LOO classifier.

Regarding the experiments with small-size training data input, the CSVMs optimized by the training data of fixed size 15 provided relatively poor diagnostic results on the testing data (the entire VAG data set), as depicted in Figure 6. The highest *A*_{z} value produced by CSVM4 was only 0.692 (SE: 0.0571), and the worst *A*_{z} value produced by CSVM2 was 0.5026 (SE: 0.0637), the latter of which was just slightly better than a wild guess. However, the DWF-based ensemble was not affected by its CSVMs which had difficulty in the case of lacking enough training data. The *A*_{z} value under the ROC curve provided by the DWF-based ensemble still reached as high as 0.9494, with the SE value of 0.0345. Such results indicated that the DWF-based ensemble had a better generalization capability when coping with the training data of smaller size.

The Natom and TCFT features derived in the present study indicated the significance of difference (*P* < 0.01) between the normal and abnormal VAG signals, and helped the classifiers improve their diagnostic performance. In the previous work [12], the signal turns were adaptively determined by the threshold that was a half of the standard deviation of the VAG signal analyzed. As the variance of an abnormal VAG signal is commonly larger than that of a normal signal, the number of signal turns detected from the abnormal signal is smaller, which cannot describe the oscillations in the abnormal signal. To overcome such a drawback, the lowpass Butterworth filter and the fixed threshold method were introduced in the present study. With the fixed threshold, more signal turns would be detected from the abnormal VAG signals, which revealed the degree of higher variability in amplitude resulting from different knee joint disorders. Compared with the adaptive threshold, the fixed threshold could help physicians establish a straightforward but effective discriminant criterion in knee joint VAG signal monitoring.

Regarding the signal classification, the results of the single LS-SVM/LOO are much better than the previous studies using the logistic regression analysis with AR coefficients as features (accuracy: 68.9%) [17]; or with the energy, energy spread, frequency, and frequency spread features derived from the Gabor MP method (accuracy: 68.9%, *A*_{z}: 0.68) [9]. The improvement of classification accuracy was largely contributed by the Natom (*P* = 0.0002) and TCFT (*P* = 0.0013) features developed in the present study.

Although the accurate rate of the Bagging 8.99% was lower than the LS-SVM/LOO, the ensemble provided an *A*_{z} of 0.8483, which was much larger than that of the LS-SVM/LOO (*A*_{z}: 0.7523), which implied that the Bagging performed better on the prediction of the true positive (abnormal) signals associated with the knee joint disorders. On the other hand, it is worth noting that the Bagging ensemble did not outperform the five CSVMs. The Bagging ensemble was worse in accurate rate than either CSVM2 or CSVM5, and its *A*_{z} value was lower than that of CSVM1. The reason why the Bagging did not achieve the better performance may be related to the robustness of the LS-SVM. According to the remarks of Breiman's work [48], “Bagging stable classifiers is not a good idea”, because the Bagging ensemble with stable classifiers can only slightly improve the accuracy, but leading to more computational complexity.

The DWF-based ensemble produced higher classification accuracy and better ROC diagnostic performance than any of the CSVMs, along with the Bagging ensemble and the single LS-SVM/LOO. The ROC curve result of the DWF-based ensemble is consistently better than that of the RBFN classifier with the features of form factors, skewness, kurtosis, and entropy (*A*_{z}: 0.8172) [11] or with the features of variance of mean-squared values and turns count with adaptive threshold (*A*_{z}: 0.9174) [12] in the recent related studies. In addition, the DWF method is also comparable to the nonlinear strict 2-surface proximal classifier (*A*_{z}: 0.95) proposed by Mu et al. [21]. The fusion weights in the DWF-based ensemble were dynamically optimized, which guaranteed the superiority of diagnostic performance.

Analysis of VAG signals using advanced digital signal processing and pattern recognition techniques is able to provide distinct indicators of degenerative articular cartilage surfaces, and has high potential for noninvasive detection of knee joint pathology [53, 54]. In the feature extraction experiments, the features of Natom and TCFT, respectively, derived from the time-frequency wavelet MP decomposition and time-domain signal variability analysis are separable with significant *P* values. Using these features, the classification by means of the LS-SVM/LOO is superior to the logistic regression analysis used in the previous studies. In addition, we utilized the ensembles of classifiers to effectively ameliorate the overall classification performance. Compared with the most popular Bagging algorithm, the DWF-based ensemble used in the present study can significantly improve the classification accuracy and the ROC curve with higher *A*_{z} and lower SE values, over the entire VAG data set.

The authors would like to thank Dr. Ranagaraj M. Rangayyan, Dr. Cyril B. Frank, and Dr. Gordon D. Bell from University of Calgary, for the work of data acquisition. This work was supported in part by the National Natural Science Foundation of China under Grant no. 81101115, the Natural Science Foundation of Fujian Province of China under Grant no. 2011J01371, the Fundamental Research Funds for the Central Universities of China under Grant no. 2010121061, the Natural Sciences and Engineering Research Council of Canada (NSERC), and the Canada Research Chairs Program. The authors are also grateful to the anonymous reviewers for the constructive comments. The preliminary works were presented at the 16th International Conference on Digital Signal Processing and the 2009 International Conference on Computational Intelligence for Measurement Systems and Applications.

1. Wu Y, Krishnan S, Rangayyan RM. Computer-aided diagnosis of knee-joint disorders via vibroarthrographic signal analysis: a review. *Critical Reviews in Biomedical Engineering*. 2010;38(2):201–224. [PubMed]

2. Frank CB, Rangayyan RM, Bell GD. Analysis of knee joint sound signals for non-invasive diagnosis of cartilage pathology. *IEEE Engineering in Medicine and Biology Magazine*. 1990;9(1):65–68. [PubMed]

3. Vigorita VJ. *Orthopaedic Pathology*. Philadelphia, PA, USA: Lippincott Williams and Wilkins; 1999.

4. McCoy GF, McCrea JD, Beverland DE. Vibration arthrography as a diagnostic aid in diseases of the knee. A preliminary report. *Journal of Bone and Joint Surgery B*. 1987;69(2):288–293. [PubMed]

5. Manaster BJ, Crim J, Rosenberg ZS. *Diagnostic and Surgical Imaging Anatomy: Knee, Ankle, Foot*. Philadelphia, PA, USA: Lippincott Williams and Wilkins; 2007.

6. Eckstein F, Adam C, Sittek H, et al. Non-invasive determination of cartilage thickness throughout joint surfaces using magnetic resonance imaging. *Journal of Biomechanics*. 1997;30(3):285–289. [PubMed]

7. Kernohan WG, Barr DA, McCoy GF, Mollan RAB. Vibration arthrometry in assessment of knee disorders: the problem of angular velocity. *Journal of Biomedical Engineering*. 1991;13(1):35–38. [PubMed]

8. Rangayyan RM, Krishnan S, Bell GD, Frank CB, Ladly KO. Parametric representation and screening of knee joint vibroarthrographic signals. *IEEE Transactions on Biomedical Engineering*. 1997;44(11):1068–1074. [PubMed]

9. Krishnan S, Rangayyan RM, Bell GD, Frank CB. Adaptive time-frequency analysis of knee joint vibroarthrographic signals for noninvasive screening of articular cartilage pathology. *IEEE Transactions on Biomedical Engineering*. 2000;47(6):773–783. [PubMed]

10. Krishnan S, Rangayyan RM, Bell GD, Frank CB. Auditory display of knee-joint vibration signals. *Journal of the Acoustical Society of America*. 2001;110(6):3292–3304. [PubMed]

11. Rangayyan RM, Wu YF. Screening of knee-joint vibroarthrographic signals using statistical parameters and radial basis functions. *Medical and Biological Engineering and Computing*. 2008;46(3):223–232. [PubMed]

12. Rangayyan RM, Wu Y. Analysis of vibroarthrographic signals with features related to signal variability and radial-basis functions. *Annals of Biomedical Engineering*. 2009;37(1):156–163. [PubMed]

13. Rangayyan RM, Wu Y. Screening of knee-joint vibroarthrographic signals using probability density functions estimated with Parzen windows. *Biomedical Signal Processing and Control*. 2010;5(1):53–58.

14. Umapathy K, Krishnan S. Modified local discriminant bases algorithm and its application in analysis of human knee joint vibration signals. *IEEE Transactions on Biomedical Engineering*. 2006;53(3):517–523. [PubMed]

15. Kim KS, Seo JH, Kang JU, Song CG. An enhanced algorithm for knee joint sound classification using feature extraction based on time-frequency analysis. *Computer Methods and Programs in Biomedicine*. 2009;94(2):198–206. [PubMed]

16. Wu Y, Krishnan S. Combining least-squares support vector machines for classification of biomedical signals: a case study with knee-joint vibroarthrographic signals. *Journal of Experimental and Theoretical Artificial Intelligence*. 2011;23(1):63–77.

17. Krishnan S, Rangayyan RM, Bell GD, Frank CB, Ladly KO. Adaptive filtering, modelling and classification of knee joint vibroarthrographic signals for non-invasive diagnosis of articular cartilage pathology. *Medical and Biological Engineering and Computing*. 1997;35(6):677–684. [PubMed]

18. Tavathia S, Rangayyan RM, Frank CB, Bell GD, Ladly KO, Zhang YT. Analysis of knee vibration signals using linear prediction. *IEEE Transactions on Biomedical Engineering*. 1992;39(9):959–970. [PubMed]

19. Moussavi ZMK, Rangayyan RM, Bell GD, Frank CB, Ladly KO, Zhang YT. Screening of vibroarthrographic signals via adaptive segmentation and linear prediation modeling. *IEEE Transactions on Biomedical Engineering*. 1996;43(1):15–23. [PubMed]

20. Jiang CC, Lee JH, Yuan TT. Vibration arthrometry in the patients with failed total knee replacement. *IEEE Transactions on Biomedical Engineering*. 2000;47(2):219–227. [PubMed]

21. Mu T, Nandi AK, Rangayyan RM. Screening of knee-joint vibroarthrographic signals using the strict 2-surface proximal classifier and genetic algorithm. *Computers in Biology and Medicine*. 2008;38(10):1103–1111. [PubMed]

22. Mallat SG, Zhang Z. Matching pursuits with time-frequency dictionaries. *IEEE Transactions on Signal Processing*. 1993;41(12):3397–3415.

23. Wu Y, Krishnan S. Classification of knee-joint vibroarthrographic signals using time-domain and time-frequency domain features and least-squares support vector machine. Proceedings of the 16th International Conference on Digital Signal Processing (DSP '09); July 2009; Santorini, Greece. pp. 361–366.

24. Daubechies I. Wavelet transform, time-frequency localization and signal analysis. *IEEE Transactions on Information Theory*. 1990;36(5):961–1005.

25. Mallat S. *A Wavelet Tour of Signal Processing*. 2nd edition. San Diego, CA, USA: Academic Press; 1999.

26. Krishnan S, Rangayyan RM. Automatic de-noising of knee-joint vibration signals using adaptive time-frequency representations. *Medical and Biological Engineering and Computing*. 2000;38(1):2–8. [PubMed]

27. Hahn GJ, Shapiro SS. *Statistical Models in Engineering*. Hoboken, NJ, USA: Wiley; 1994.

28. Rangayyan RM. *Biomedical Signal Analysis: A Case-Study Approach*. New York, NY, USA: IEEE and Wiley; 2002.

29. Willison RG. Analysis of electrical activity in healthy and dystrophic muscle in man. *Journal of Neurology, Neurosurgery, and Psychiatry*. 1964;27(5):386–394. [PMC free article] [PubMed]

30. Rangayyan RM, Oloumi F, Wu Y, Cai S. Fractal analysis of knee-joint vibroarthrographic signals via power spectral analysis. *Biomedical Signal Processing and Control*. 2013;8(1):23–29.

31. Cortes C, Vapnik V. Support-vector networks. *Machine Learning*. 1995;20(3):273–297.

32. Vapnik VN. *Statistical Learning Theory*. New York, NY, USA: Wiley; 1998.

33. Haykin S. *Neural Networks: A Comprehensive Foundation*. 2nd edition. Englewood Cliffs, NJ, USA: Prentice Hall PTR; 1998.

34. Duda RO, Hart PE, Stork DG. *Pattern Classification*. 2nd ed edition. New York, NY: Wiley; 2001.

35. Suykens JAK, Van Gestel T, De Brabanter J, De Moor B, Vandewalle J. *Least Squares Support Vector Machines*. Singapore: World Scientific Publishing; 2002.

36. Nash SG, Sofer A. *Linear and Nonlinear Programming*. Columbus, OH, USA: McGraw-Hill; 1995.

37. Wu Y, Krishnan S. Computer-aided analysis of gait rhythm fluctuations in amyotrophic lateral sclerosis. *Medical and Biological Engineering and Computing*. 2009;47(11):1165–1171. [PubMed]

38. Wu Y, Krishnan S. Statistical analysis of gait rhythm in patients with Parkinson’s disease. *IEEE Transactions on Neural Systems and Rehabilitation Engineering*. 2010;18(2):150–158. [PubMed]

39. Wu Y, Shi L. Analysis of altered gait cycle duration in amyotrophic lateral sclerosis based on nonparametric probability density function estimation. *Medical Engineering and Physics*. 2011;33(3):347–355. [PubMed]

40. Hansen LK, Salamon P. Neural network ensembles. *IEEE Transactions on Pattern Analysis and Machine Intelligence*. 1990;12(10):993–1001.

41. Hashem S, Schmeiser B. Improving model accuracy using optimal linear combinations of trained neural networks. *IEEE Transactions on Neural Networks*. 1995;6(3):792–794. [PubMed]

42. Kittler J, Hatef M, Duin RPW, Matus J. On combining classifiers. *IEEE Transactions on Pattern Analysis and Machine Intelligence*. 1998;20(3):226–239.

43. Ueda N. Optimal linear combination of neural networks for improving classification performance. *IEEE Transactions on Pattern Analysis and Machine Intelligence*. 2000;22(2):207–215.

44. Kuncheva LI. A theoretical study on six classifier fusion strategies. *IEEE Transactions on Pattern Analysis and Machine Intelligence*. 2002;24(2):281–286.

45. Wu Y, Wang C, Ng SC, Madabhushi A, Zhong YX. Breast cancer diagnosis using neural-based linear fusion strategies. Proceedings of the 13th International Conference on Neural Information Processing (ICONIP '06); October 2006; Hong Kong. pp. 165–175.

46. Sinha A, Chen H, Danu DG, Kirubarajan T, Farooq M. Estimation and decision fusion: a survey. *Neurocomputing*. 2008;71(13-15):2650–2656.

47. Freund Y, Schapire RE. A Decision-Theoretic Generalization of On-Line Learning and an Application to Boosting. *Journal of Computer and System Sciences*. 1997;55(1):119–139.

48. Breiman L. Bagging predictors. *Machine Learning*. 1996;24(2):123–140.

49. Ridgeway G. The State of boosting. *Computing Science and Statistics*. 1999;31:172–181.

50. Efron B, Tibshirani R. *An Introduction to the Bootstrap*. New York, NY, USA: Chapman and Hall; 1993.

51. Wu Y, Krishnan S. An adaptive classifier fusion method for analysis of knee-joint vibroarthrographic signals. Proceedings of the IEEE International Conference on Computational Intelligence for Measurement Systems and Applications (CIMSA '09); May 2009; Hong Kong. pp. 190–193.

52. Girolami M. Orthogonal series density estimation and the kernel eigenvalue problem. *Neural Computation*. 2002;14(3):669–688. [PubMed]

53. Cai S, Wu Y, Xiang N, et al. Detrending knee joint vibration signals with a cascade moving average filter. Proceedings of the 34th Annual International Conference of IEEE Engineering in Medicine and Biology Society (EMBC '12); August 2012; San Diego, CA, USA. pp. 4357–4360. [PubMed]

54. Wu Y, Cai S, Xu F, Shi L, Krishnan S. Chondromalacia patellae detection by analysis of intrinsic mode functions in knee joint vibration signals. IFMBE Proceedings of 2012 World Congress on Medical Physics and Biomedical Engineering (WC '12); May 2012; Beijing, China. pp. 493–496.

Articles from Computational and Mathematical Methods in Medicine are provided here courtesy of **Hindawi Publishing Corporation**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |