Julie Bertrand
email: julie.bertrand@inserm.fr

Emmanuelle Comets

Marylore Chenel

France Mentré

Some Alternatives to Asymptotic Tests for the Analysis of Pharmacogenetic Data using Nonlinear Mixed Effects Models

Keywords: F -distribution based approach, First-order approximation, Nonlinear mixed effects models, Permutation tests, Pharmacogenetics

Nonlinear mixed effects models allow investigating individual differences in drug concentration profiles (pharmacokinetics) and responses. Pharmacogenetics focusses on the genetic component of this variability. Two tests often used to detect a gene effect on a pharmacokinetic parameter are i) the Wald test, assessing whether estimates for the gene effect are significantly different from 0 and ii) the likelihood ratio test comparing models with and without the genetic effect.

Because those asymptotic tests show inflated type I error on small sample size and/or with unevenly distributed genotypes, we develop two alternatives and evaluate them by means of a simulation study.

First, we assess the performance of the permutation test using the Wald and the likelihood ratio statistics. Second, for the Wald test we propose the use of the F -distribution with four different values for the denominator degrees of freedom. We also explore the influence of the estimation algorithm using both the First-Order Conditional Estimation with interaction linearisation-based algorithm and the Stochastic Approximation Expectation Maximization algorithm. We apply these methods to the analysis of the pharmacogenetics of indinavir in HIV patients recruited in the COPHAR2-ANRS 111 trial. Results of the simulation study show that the permutation test seems appropriate but at the cost of an additional computational burden. One of the four F -distribution based approaches provides a correct type I error estimate for the Wald test and should be further investigated.

Introduction

Pharmacokinetics (PK) studies the time course of a drug in the body. Nonlinear mixed effects models (NLMEM) in the analysis of PK data allow integrating the knowledge accumulated on drug absorption, distribution, metabolism and elimination (ADME), to quantify the inter-individual variability with fewer samples per patient than standard non compartmental approaches. Pharmacogenetics (PG) studies the relationship between this inter-individual variability and variations in the DNA sequence of proteins involved in the ADME mechanisms of the drug.

In previous simulation studies, we have shown that a correction for type I error inflation was required for the two most commonly used tests for differences between groups in NLMEM, the Wald and the likelihood ratio tests (LRT), in data sets of small sample size and/or with unevenly distributed genotypes [START_REF] Bertrand | Comparison of model-based tests and selection strategies to detect genetic polymorphisms influencing pharmacokinetic parameters[END_REF](Bertrand et al., , 2009b)). Also in the context of testing parametric against semi-parametric nonlinear mixed models, [START_REF] Wu | The study of long-term HIV dynamics using semi-parametric non-linear mixed-effects models[END_REF] have shown that asymptotic distributions cannot be used for the Wald test and the LRT when the normality assumption of the random effects is violated, though the LRT was much less affected.

The aim of the present work is to investigate two alternatives to these asymptotic tests to detect a gene effect in PG studies: permutation and F -distribution based tests. Permutation testing is a way of determining whether the null hypothesis of randomness is reasonable, i.e. whether the pattern present in the data could have happen by chance alone [START_REF] Good | Permutation tests: a practical guide to resampling methods for testing hypotheses[END_REF][START_REF] Manly | Randomization, Bootstrap and Monte Carlo Methods in Biology[END_REF]. This alternative requires fewer assumptions than a correction based on simulations from the model under the null hypothesis, but is seldom used in NLMEM [START_REF] Ding | Assessing antiviral potency of anti-HIV therapies in vivo by comparing viral decay rates in viral dynamic models[END_REF]. The second alternative is to correct the Wald test for the underestimation of the variance of parameter estimates in small sample size. We investigate four such corrections based on F -distributions with various values of denominator degrees of freedom (df) which have been proposed in different settings.

In the present work, we study the type I error and power of these two alternatives through simulations. The permutation tests are evaluated for both the Wald and the LRT statistics whereas the F -distributions are used to correct the Wald test only. The same simulation setting as in the simulation studies mentioned before [START_REF] Bertrand | Comparison of model-based tests and selection strategies to detect genetic polymorphisms influencing pharmacokinetic parameters[END_REF](Bertrand et al., , 2009b) ) was used here to enable a comparison with the results of the simulation-based correction studied in these works. As an illustration, we present an analysis of the indinavir concentration-time profiles collected in the COPHAR2 trial.

As in NLMEM the integral in the likelihood has no analytical form, estimation algorithms use model linearisation, such as the first order (FO) and first order conditional estimation with interaction (FOCE-I) [START_REF] Lindstrom | Nonlinear mixed effects models for repeated measures data[END_REF] methods implemented in the NON-MEM and the Phoenix NLME software. Others consider numerical approximations of the likelihood such as the Laplacian or Adaptive Gaussian quadrature algorithms (AGQ) or more recently stochastic approximation of the EM algorithm (SAEM) [START_REF] Deylon | Convergence of a stochastic approximation version of EM algorithm[END_REF] as implemented in the MONOLIX software. In [START_REF] Bertrand | Comparison of model-based tests and selection strategies to detect genetic polymorphisms influencing pharmacokinetic parameters[END_REF]Bertrand et al. ( , 2009b)), we have shown that the bias in the FO algorithm leads to a very large inflation of both asymptotic tests type I error. Thus, to account for the influence of the estimation algorithm in the present work we use FOCE-I and SAEM.

In Section 2, we present the model and how the likelihood and the estimation variance matrix are obtained. Then, we introduce the usual asymptotic tests and both investigated alternatives in Section 3. In Sections 4 and 5 we describe the real data, the simulation study and the evaluation protocol. The results of the evaluation and the illustration are presented in Section 6. We finally discuss our findings and conclusions in Sections 7 and 8.

Evaluation of the likelihood and the estimation variance matrix in NLMEM

Models and notation

To describe the n i -dimensional vector of concentrations y i of a subject i = 1, ..., N , we use a pharmacokinetic function f which is nonlinear in its parameters φ i :

y i = f (X i ; φ i ) + ǫ i , (1)
where X i is the within-subject design vector and

φ i = µ + A i β + Bη i is the subject p-dimensional vector of parameters, A i is the (p × k)-dimensional covariate matrix
modelling the relationship between the covariates and φ i with β the corresponding kcoefficient effects vector. Because the covariate coefficients enter the model nonlinearly through the structural parameter associated to the covariate, the matrix-based notation peculiar to linear mixed effect models can no longer be used and therefore θ, the (p + k)-

dimensional fixed effects vector is decomposed into [µ ′ β ′ ] ′ . B is a (p × q)-dimensional
design matrix, permitting some components of φ i to have no associated random effect when p > q. η i is the random effect q-dimensional vector which follows a Gaussian distribution with null mean and variance-covariance (q×q)-dimensional matrix Ω. ǫ i is the residual error n i -dimensional vector which follows a Gaussian distribution with null mean and variancecovariance (n i × n i )-dimensional matrix Σ i (X i ; φ i , γ) whose diagonal terms are equal to g(X i ; φ i , γ) 2 and off-diagonal terms are equal to 0.

To obtain a combination of constant and proportional error models, g(X i ; φ i , γ) can be set equal to a + bf (X i ; φ i ) with γ = [a b] ′ the vector of the error model parameters. Define the l-dimensional vector of variance parameters as λ = [Vech(Ω) ′ γ ′ ] ′ where the operator Vech(.) creates a column vector from the matrix Ω by stacking its lower diagonal elements below one another. Finally, define the (p+k +l)-dimensional vector of all model parameters

as Ψ = [θ ′ λ ′ ] ′ .

Evaluation of the likelihood

Because of the nonlinearity of the regression function in the random effects, the likelihood of NLMEM cannot be expressed in a closed form. Indeed, for the subject i, the marginal loglikelihood L i (y i ; Ψ) of Ψ for the data y i is given by

L i (y i ; Ψ) =log p(y i ; φ i , Ψ)dφ i =log p(y i |φ i , Ψ)p(φ i ; Ψ)dφ i (2)
with p(y i |φ i , Ψ) the density of the observations conditioned on the random effects and the model parameters, p(φ i ; Ψ) the density of the individual parameters and p(y i ; φ i , Ψ) the likelihood of the complete data, that is,

-2 log p(y i ; φ i , Ψ) =n i log(2π) + log(|Σ i (X i ; φ i , γ)|) +(y i -f (X i ; φ i )) ′ Σ i (X i ; φ i , γ) -1 (y i -f (X i ; φ i )) + q log(2π) +log(|Ω * |)(φ i -(µ + A i β)) ′ Ω * -1 (φ i -(µ + A i β))
where Ω * = B ′ ΩB.

Different approaches have been proposed to estimate L i (y i ; Ψ). Linearisation based methods consist of taking a Taylor expansion of the model function f and approximating (2) by the loglikelihood of a linear mixed effect model:

-2 L i (y i ; Ψ) = n i log(2π) + log(|V i |) + (y i -E i ) ′ V i -1 (y i -E i )
where E i and V i are respectively the marginal expectation and variance of the vector y i .

Another approach is to use an importance sampling procedure [START_REF] Robert | Monte Carlo Statistical Methods[END_REF] to compute an estimate L i (y i ; Ψ) T of the observed loglikelihood, such as

L i (y i ; Ψ) T = log 1 T T t=1 p(y i |φ i (t) , Ψ) p(φ i (t) ; Ψ) h(φ i (t) )
where the φ i (t) are sampled from an instrumental distribution h(φ i (t) ) which is chosen to minimize the variance of the estimate L i (y i ; Ψ) T . In the MONOLIX software where such a procedure is used, the instrumental function is a non-centred Student distribution [START_REF] Samson | The SAEM algorithm for group comparison tests in longitudinal data analysis based on non-linear mixed-effects model[END_REF][START_REF] Lavielle | MONOLIX (MOdèles NOn LInéaires à effets miXtes)[END_REF].

Evaluation of the estimation variance matrix

The estimation variance matrix is composed of

Var(Ψ) =    Var(θ) Var(θ, λ) Var(θ, λ) ′ Var(λ)   
where Var(θ) is the ((p + k) × (p + k))-estimation variance matrix for the fixed effects, Var(λ) is the (l × l)-estimation variance matrix for the variance components, and Var(θ, λ)

is the ((p + k) × l)-estimation covariance matrix between the fixed effects and the variance components. Based on the Cramer-Rao inequality, the inverse of the Fisher estimation matrix M F is the lower bound of the variance covariance matrix of any unbiased estimators of the parameters. In the framework of normal theory maximum likelihood, M F is computed as the Hessian of the loglikelihood in all the model parameters, and thus

M F = N i=1 -∂ 2 L i (y i ; Ψ) ∂Ψ∂Ψ ′
where N is the total number of subjects. Thus, the (asymptotic) estimate of Var(Ψ) is

M F -1 .
Here as well, linearisation of the model has been proposed to derive an approximate expression of M F . The FOCE-I approach linearises the model around the individual random effect estimates ηi [START_REF] Pinheiro | Mixed Effects Models in S and S-Plus[END_REF].

y i = f (X i ; φi ) + ∂f (X i ; φ i ) ∂η i ′ φi (η i -ηi ) + ǫ i with φi = µ + A i β + ηi .
The approximate marginal expectation E i and variance V i of the vector y i are then given by

E i = f (X i ; φi ) - ∂f (X i ; φ i ) ∂η i ′ φi ηi V i = ∂f (X i ; φ i ) ∂η i ′ φi Ω * ∂f (X i ; φ i ) ∂η i φi + Σ i (X i ; φi , γ) (3)
Alternatively, a Taylor expansion of the model function f can also be performed around the individual parameter estimates φi = μ + A i β + ηi , which provides the following approximation. Then, (1) is rewritten as

y i = f (X i ; φi ) + ∂f (X i ; φ i ) ∂φ i ′ φi (φ i -φi ) + ǫ i
which is reformulated for sake of simplicity:

ỹi = ∂f (X i ; φ i ) ∂φ i ′ φi φ i + ǫ i with ỹi = y i -f (X i ; φi ) + ∂f (X i ;φ i ) ∂φ i ′ φi φi .
Then, the approximate marginal expectation E i and variance V i of the vector ỹi are given by

E i = ∂f (X i ; φ i ) ∂φ i ′ φi (µ + A i β) V i = ∂f (X i ; φ i ) ∂φ i ′ φi Ω * ∂f (X i ; φ i ) ∂φ i φi + Σ i (X i , φi , γ) (4)
This latter approach is implemented in the MONOLIX software.

As we need its expression in the following, let us write here the estimation variance for the fixed effects. Using a first order linearisation of the model where E i and V i are given in either ( 3) or ( 4), we have Var(θ, λ) = 0 and thus

Var(θ) = N i=1 ∂E i ∂θ ′ V i -1 ∂E i ∂θ -1
(5)

Asymptotic tests and alternatives

Asymptotic tests

In asymptotic conditions, testing the null hypothesis H 0 : Cθ = 0 based on the estimates of the fixed effects θ can be carried out with the usual Wald test, comparing the Wald statistic to the critical value of a χ 2 distribution. Here, C is a ((p + k) × U )-contrast matrix and U the number of contrasts one wishes to test. If H 0 is true, then:

W = (C θ) ′ (C ′ Var(θ)C) -1 (C θ) ∼ χ 2 U (6)
H 0 can also be tested using the likelihood ratio test (LRT) that compares the models with Cθ = 0 (reduced) with the full model where θ is estimated. If H 0 is true, then:

S = -2 × (L reduced -L f ull ) ∼ χ 2 U
where L reduced and L f ull are the loglikelihood of the two models.

Permutation based alternative

To perform permutation tests, R data sets are generated by permuting the rows of the covariate matrix from the original data set. For a given test, the statistic Q obs is estimated from the original data and the statistic Q perm is estimated from each of the R data sets.

Thus, r = 1, ..., R values Q permr are obtained which constitute a distribution of the statistic under the null hypothesis of no covariate effect. The permutation p-value is the proportion:

(Card(Q permr Q obs ) + 1)/(R + 1), where the operator Card(.) counts the realisations of

Q permr Q obs .

F -distribution based alternative

Comparing W/U to an F -distribution with an infinite denominator df is equivalent to using the classical Wald test. However, in biology studies, sample size is often small to moderate, thus approximate F -distributions with numerator df= U and non infinite denominator df have been proposed to correct for the departure from the asymptotic. In this section,we consider four different F -distribution based approaches.

The first approach derives from decomposition of the degrees of freedom in the ANOVA (Pinheiro and Bates, 2000) (=DF PB ) with a denominator df equal to N i=1 n i -(N +p+k-U ) and is implemented in the nlme function in R. In the second approach which is implemented in the NLMIXED Procedure in SAS, the computation of the degrees of freedom is based on the number of random effects [START_REF] Wolfinger | Fitting nonlinear mixed models with the new NLMIXED procedure[END_REF] (=DF W ) with a denominator df equal to N -q.

The third approach comes from the multivariate nonlinear models (MNLM) framework [START_REF] Gallant | Seemingly unrelated nonlinear regressions[END_REF]) (=DF G ). In his study on nonlinear regressions contemporaneously but not serially correlated, [START_REF] Gallant | Seemingly unrelated nonlinear regressions[END_REF] observed that estimation variances were underestimated and thus recommended the multiplication of Var(θ) by the factor N/(N -p) and using a denominator df equal to N -p. In classical non linear regression, the downward bias of the maximum likelihood variance estimator is corrected using a factor of n/(n-p), with n being the number of observations and p being the number of model parameters [START_REF] Huet | Statistical Tools for Nonlinear Regression: A Practical Guide with S-PLUS Examples[END_REF]. Here, [START_REF] Gallant | Seemingly unrelated nonlinear regressions[END_REF] proposed substituting N and not N i=1 n i for n because MNLM, as a populationaveraged approach, focusses on the marginal expectation of the response variable [START_REF] Vonesh | Linear and Nonlinear Models for the Analysis of Repeated Measurements[END_REF].

The fourth approach is an extension to NLMEM of the method developed by [START_REF] Fai | Approximate F-tests of multiple degree of freedom hypotheses in generalized least squares analyses of unbalanced split-plot experiments[END_REF] (=DF FC ) that is implemented in the DDFM=SATTERTH option of the MIXED Procedure in SAS. [START_REF] Fai | Approximate F-tests of multiple degree of freedom hypotheses in generalized least squares analyses of unbalanced split-plot experiments[END_REF] proposed using the spectral decompo-

sition of C ′ Var(θ)C in (6) to decompose W into a sum of U squared Student statistics with of df υ u = 2(c ′ u Var(θ)c u ) 2 Var(c ′ u Var(θ)c u ) with u = 1, .
., U the decomposition index and where Var(c ′ u Var(θ)c u ) can be approximated by the delta method using the estimates of Var(λ) and the estimates of the model parameters:

Var(c ′ u Var(θ)c u ) ≈ ∂c ′ u Var(θ)c u ∂λ ′ Var(λ) ∂c ′ u Var(θ)c u ∂λ (7)
When U > 1, [START_REF] Fai | Approximate F-tests of multiple degree of freedom hypotheses in generalized least squares analyses of unbalanced split-plot experiments[END_REF] showed that W/U follows an approximate Fdistribution with denominator df:

υ = 2 U u=1 υ u υ u -2 U u=1 υ u υ u -2 -U
As, in NLMEM, Var(θ) includes derivatives, we use the following property of matrix

derivatives dA -1 dx = -A -1 dA dx A -1 to take the derivative of c ′ u Var(θ)c u in (7) with respect
to each element l of the vector λ:

∂c ′ u Var(θ)c u ∂λ l = c ′ u   -Var(θ) ∂ N i=1 ∂E i ∂θ ′ V i -1 ∂E i ∂θ ∂λ l Var(θ)   c u (8)
We also use this property to take the derivative of N i=1 ∂E i ∂θ

′ V i -1 ∂E i
∂θ in (8) with respect to each element l of λ, so we have, for each element q of V ech(Ω),

∂ N i=1 ∂E i ∂θ ′ V i -1 ∂E i ∂θ ∂V ech(Ω) q = N i=1 ∂E i ∂θ ′ -V i -1 ∂f (X i ; φ i ) ∂φ i ′ φi ∂f (X i ; φ i ) ∂φ i φi V i -1 ∂E i ∂θ
and for the error model parameters,

∂ N i=1 ∂E i ∂θ ′ V i -1 ∂E i ∂θ ∂a 2 = N i=1 ∂E i ∂θ ′ -V i -1 diag 1 + bf (X i ; φi ) a V i -1 ∂E i ∂θ ∂ N i=1 ∂E i ∂θ ′ V i -1 ∂E i ∂θ ∂b 2 = N i=1 ∂E i ∂θ ′ -V i -1 diag af (X i ; φi ) b + f (X i ; φi ) 2 V i -1 ∂E i ∂θ
In these calculations, we linearise the model around the individual parameters estimates as described in Section 2.3.

Real Data and Simulation Study

Real data

As in our previous papers [START_REF] Bertrand | Comparison of model-based tests and selection strategies to detect genetic polymorphisms influencing pharmacokinetic parameters[END_REF](Bertrand et al., , 2009b)), we illustrate the different approaches with data from a PK substudy of the COPHAR2-ANRS 111 study, a multicentre noncomparative pilot trial of early therapeutic drug monitoring in HIV-positive patients naïve of treatment. We focus on the PK sub-study from the group of patients receiving indinavir doses of 400, 600 and 800 mg boosted with ritonavir at a dose of 100 mg twice a day. Patients were genotyped for the exons 21 and 26 of the ABCB1 gene which code for the P-glycoprotein, and for polymorphisms on gene coding for proteins involved in the metabolism of indinavir;

the CYP 3A4*1B, CYP 3A5*3 and *6 polymorphisms.

The indinavir pharmacokinetic profiles were determined at 1, 3, 6, and 12h following drug administration, 2 weeks after the treatment onset (Fig. 1).

[Figure 1 about here.]

Simulation setting

An extended description of the simulations can be found in [START_REF] Bertrand | Comparison of model-based tests and selection strategies to detect genetic polymorphisms influencing pharmacokinetic parameters[END_REF] and is briefly summarised below.

Parameters from the simulated NLMEM were set based on a preliminary analysis of the indinavir data without covariates. The concentrations at time t were simulated using a one compartment model at steady state (τ = 12 h) with first order absorption (k a = 1.4 h -1 ), first order elimination (k = 0.2 h -1 ) and apparent volume of distribution (V /F = 102 L).

f (t; k a , k, V /F ) = D V /F k a k a -k e -kt 1 -e -kτ - e -kat 1 -e -kaτ
We use a diagonal Ω matrix and a proportional error model (g(t; φ i , γ) = bf (t, φ i )) setting the dose (D) to 400 mg. The model parameters are expressed in term of natural logarithms to achieve the positivity requirement of pharmacokinetic parameters, µ = [log(k a ) log(k) log(V /F )] ′ with random effect standard deviations set to 1.13, 0.41 and 0.26, respectively and the coefficient of variation for the residual error set to 20% (b=0.2).

We simulate a diplotype of SN P 1 and SN P 2 , with a distribution mimicking that of exons 26 and 21 of the ABCB1 gene [START_REF] Sakaeda | MDR1 genotype-related pharmacokinetics and pharmacodynamics[END_REF]; we obtained frequencies for the rare homozygotes, heterozygotes and common homozygotes of respectively 24%, 48% and 28% for SN P 1 and 29%, 44% and 27% for SN P 2 . As, in the intestine, the P-glycoprotein restricts drug entry into the body, we consider an effect on the drug bioavailability (F) through the volume of distribution V/F.

The genetic coefficient values are chosen to be consistent with results found in the literature for ABCB1 polymorphisms on drugs disposition [START_REF] Marzolini | Polymorphisms in human MDR1 (P-glycoprotein): recent advances and clinical relevance[END_REF] and provide clinically relevant effect, with V/F and CL/F (=k × V /F ) increasing from 105.4 to 200.5 L and 21.1 to 40.1 L/h respectively between common and rare homozygotes for SN P 1 .

In this work, the tests under study assess only the effect of SN P 1 on V /F even if we simulated diplotypes to allow for population genetic factors such as linkage disequilibrium.

In the model including the covariate, the effect coefficient vector is β = [β 1 β 2 ] ′ . As we test for the effect of the SNP as a whole, we have U = 2 and C =    0 0 0 1 0

0 0 0 0 1    ′ .
We simulated 200 data sets with N = 40 subjects and n = 4 samples as in the COPHAR2 trial under both the null hypothesis of no genetic effect (H 0 ) and the alternative hypothesis of a genetic effect (H 1 ) (Fig. 2).

[Figure 2 about here.]

Evaluation protocol

Simulation study

In the present research we used the FOCE-I and SAEM algorithms implemented in the softwares NONMEM version 7 and MONOLIX version 2.1 respectively to estimate the parameters, the likelihood and the estimation variance matrix. In NONMEM, the starting values were set to the simulated ones + 10% except for the standard deviation of k a which was set to 0.7. In MONOLIX, the starting values were similar for the fixed effects but the random effect standard deviations were set to 1.

For the likelihood estimation in MONOLIX, the number of iterations for the importance sampling T was set to 10000.

The alternatives proposed in the present work are evaluated in terms of type I error and power. The type I error estimate is the percentage of the 200 data sets simulated under H 0 where the test is significant and the power estimate is the same percentage but for the 200 data sets simulated under H 1 . The prediction interval around 5% for 200 data sets is [2%; 8%]. The estimates of type I error and power obtained with permutations or F -distributions are compared to values obtained using the asymptotic tests (theoretical threshold of p-value=0.05) and using a simulation based correction (threshold built from simulations under the model with Cθ = 0).

As described in Section 3, the type I error of the asymptotic tests is obtained using a threshold of 5.99 which is the 95 th percentile of a χ 2 distribution with 2 degrees of freedom.

The type I error for the simulation based correction is obtained using a threshold of 10.18 and 7.05 with FOCE-I and 7.80 and 6.96 with SAEM for the Wald test and the LRT, respectively. These values are the 95 th percentiles of the distributions obtained for both the Wald test and the LRT statistics on 1000 data sets simulated under the null hypothesis in [START_REF] Bertrand | Comparison of model-based tests and selection strategies to detect genetic polymorphisms influencing pharmacokinetic parameters[END_REF] using the same setting as described in Section 4 of the present paper.

We also evaluate the p-value distribution for both tests on the 200 data sets simulated under H 0 . If the test statistic does follow its reference distribution, these p-values should follow an uniform distribution on the interval [0, 1]. Thus, we computed the Kolmogorov distance between the uniform distribution on the interval [0, 1] and the p-value distributions observed for the Wald test and the LRT under their asymptotic form and using the proposed alternatives with both estimation algorithms [START_REF] Ruckdeschel | S4 classes for distributions[END_REF]. For a sample of 200 values, the Kolmogorov-Smirnov (K-S) test of equality of one-dimensional probability distributions is significant at a level of 5% when the Kolmogorov distance is below 1.358 √ 200 = 0.096.

For the permutation test, the number R of permutations is set to 1000.

The fourth approach based on F -distribution (DF F C ) is evaluated only with SAEM in MONOLIX. Our implementation of this denominator df calculation is nested within the MONOLIX code, following the computation of the estimation variance matrix. We did not program the extension with the FOCE-I algorithm as we had no easy access to the NONMEM code dedicated to the estimation variance matrix.

Application to real data

The indinavir concentrations were analysed with the same PK model as in the simulation study, using the SAEM algorithm. The details of the covariate model building strategy are described in Bertrand et al. (2009a). Briefly, we tested for demographic and biological covariates along with the genetic polymorphisms. First, a screening on the estimated individual parameters was performed, followed by a forward selection based on LRT. Then, the p-value of the covariates remaining in the final model were assessed using the asymptotic LRT as well as the four F -distribution based Wald tests for illustrative purpose.

Results

Simulation study

As shown on Table 1, SAEM achieved convergence on all data sets and always provided the estimation variance matrix, whereas FOCE-I failed to converge and/or compute the estimation variance matrix on few data sets.

[Table 1 about here.]

The thresholds used for the asymptotic test, the simulation based correction, and the first three F -distribution based approaches, are the same for the 200 data sets under both hypotheses. For the permutation test and the DF FC approach, the threshold comes from, respectively, the permutation based distribution and the estimated denominator df, which are both specific to the data set.

With SAEM, the type I error for the asymptotic Wald test and LRT are significantly inflated. With FOCE-I, the type I error for the asymptotic Wald test is significantly inflated whereas the asymptotic LRT estimate is at the upper limit of the prediction interval around 5%. Of note, the Wald test and LRT asymptotic estimates on the 1000 data sets used to derive the simulation-based threshold are 8.9% and 7.6% with SAEM and 12.2% and 7.4% The range of denominator df obtained with the DF FC method is quite narrow and happens to include that of the DF W approach.

Fig. 4 shows that for the Wald test with SAEM, one should use a denominator df of 10 to obtain a type I error of exactly 5%, and between 6 and 27 for an estimate non significantly different from 5%.

[Figure 4 about here.]

Real data

As described in Bertrand et al. (2009a), the indinavir concentration time data are adequately fitted by a one compartment model parametrised in first-order absorption rate, apparent volume of distribution and apparent elimination clearance (Cl/F). The model with no co-variate has an absorption constant of 1.3 h -1 with an important inter-individual variability of 118%, an elimination clearance of 21.9 L.h -1 and a volume of distribution of 93.9 L with inter-individual variabilities of 34.4% and 19.3%, respectively. The standard deviations of the random effect are expressed as the coefficient of variation of the pharmacokinetic parameters, as in the model the latter are coded in natural logarithms. estimates are in accordance with the literature on indinavir given in combination with ritonavir. The coefficient of variation for the residual error is 44.5%. All of the estimation relative standard errors (RSE) were below 25% with the exception of k a and V /F (around 30% and 60%, respectively).

After ascending selection using the asymptotic LRT, only an effect of the CYP 3A4*1B polymorphism effect on k a (p-value=0.02) and an age effect on Cl/F (p-value=0.03) remain in the model. The p-values of the permutation tests are 0.04 and 0.1 for the CYP 3A4*1B polymorphism and age effects, respectively. The corresponding p-values are 0.014 and 0.043 using the DF P B approach, 0.018 and 0.048 using the DF W approach, 0.023 and 0.057 using the DF G approach, and 0.02 and 0.047 using the DF F C approach.

Only the effect of the CYP 3A4*1B polymorphism remains in the final model, as the effect of age is discarded based on the p-value estimates from the permutation test and the DF G approach.

Discussion

Several studies have evaluated by simulation the performance of tests for discrete covariate on continuous responses using NLMEM with various designs and estimation methods (see the review in the discussion of Bertrand et al. (2009b)). Permutation tests and F -distribution based tests have been compared in the linear mixed effects framework (Lin and Heagerty, 2004;[START_REF] Routledge | P-values from permutation and F-tests[END_REF]. In particular [START_REF] Routledge | P-values from permutation and F-tests[END_REF] Because the covariates enter the model in a linear fashion, we also consider using a F -distribution for the Wald test, a widely used method in LMEM [START_REF] Verbeke | Linear Mixed Models in Practice. A SAS-Oriented Approach[END_REF]. The df of the DF PB approach which is derived from the multivariate ANOVA differs greatly from the df used within the mixed effect framework. Indeed, in the multivariate ANOVA, the subject concentrations are considered as independent variables under the influence of the genotype and the multivariate statistic is based on the comparison of the between and within variance matrices whose covariance terms inform on the correlation between the variables. Thus in this approach, n variables are considered in N subjects leading to an N × n sample size. Then, one degree of freedom is subtracted per parameter estimates minus the number of contrast tested. On the other hand in the mixed effects model, the independent variable is the subject with N the sample size. For DF W and DF G , one degree of freedom is subtracted per random effect variance and fixed effect respectively. For DF FC , the information available in the data for the parameter estimation is quantified through the variance of parameter estimates rather than by subtracting degrees of freedom to the sample size. Here, all these methods however fail to correct for the Wald type I error inflation, as the denominator df estimates for the 200 data sets are too large. Actually, the optimal denominator df is even smaller than the number of independent observations to compensate for the downward bias of the variance of the estimates, as shown by the performance of the DF G method.

Indeed, the only method based on an F -distribution that corrects for the type I error inflation of the Wald test is the DF G method with SAEM. With the FOCE-I algorithm, the Kolmogorov distance suggested an improvement using the DF G method despite the fact that we observed an inflation of the type I error, indicating that the correction may be useful even in this case. To go past the simple design under study here, we considered two other designs from a previous study evaluating the Wald test and the LRT (Bertrand et al., 2009b): a design optimised using the PFIM software [START_REF] Retout | PFIM Interface 2.1. UMR[END_REF] including N = 80 subjects sorted in 4 groups with n = 2 samples and a combined design with N = 20 subjects having n = 4 samples plus N = 80 subjects with only one trough concentration (n = 1). The type I error estimates using the classical Wald test and the DF G method are 8.7% and 5.7%, respectively for the optimal design and 8.4% and 4.5%, respectively for the combined design.

These results show that the DF G method also corrects the type I error inflation observed in these designs. However, it does so not only through approximating the distribution for the Wald test but essentially through inflating the estimation variance for the effect coefficient by the factor N/(N -p). This factor, which brings Bessel's correction for the unbiased estimator of the variance to mind, was recommended by [START_REF] Gallant | Seemingly unrelated nonlinear regressions[END_REF] based only on a Monte Carlo study in the multivariate nonlinear model framework. Thus, caution should be used and further investigation is needed as to designs involving different number of samples per subject and in other simulation settings.

One limitation of the present research is that we consider the method developed by [START_REF] Fai | Approximate F-tests of multiple degree of freedom hypotheses in generalized least squares analyses of unbalanced split-plot experiments[END_REF] and not the method developed by [START_REF] Kenward | Small sample inference for fixed effects from restricted maximum likelihood[END_REF]. The latter method is indeed based on restricted maximum likelihood (REML) which is not commonly used in NLMEM. Further, [START_REF] Faes | The effective sample size and an alternative small-sample degrees-of-freedom method[END_REF] have recently revisited the effective sample size in correlated data and provided an approximate degrees of freedom method. They used REML estimation, thus we did not assess their approach. [START_REF] Meza | REML estimation of variance parameters in nonlinear mixed effects models using the SAEM algorithm[END_REF] have developed a REML extension of the SAEM algorithm not implemented in MONOLIX yet. A perspective of the present work would therefore be to assess the [START_REF] Kenward | Small sample inference for fixed effects from restricted maximum likelihood[END_REF] and [START_REF] Faes | The effective sample size and an alternative small-sample degrees-of-freedom method[END_REF] methods using this extension. Of note, a Bayesian framework would be appealing, not the least because it would circumvent the issue of the type I error by summarizing the available information in a posterior distribution instead of a p-value.

Conclusion

While this research confirms the feasibility of permutation tests in pharmacogenetic studies for both the LRT and the Wald test, further investigations in other simulation settings are required to recommend the F -distribution based approach proposed in multivariate nonlinear models for the Wald test. As permutation does only slightly better than the asymptotic alternatives for an additional computing burden, a sensible course of action could be to use it only when decisions based respectively on the asymptotic test and the Gallant alternative are discordant. * Threshold for the distribution of W, not W/2. † The threshold for the permutation test and both the threshold and the denominator df for the DFFC approach are displayed as median [range] over the T data sets simulated under H0.

  with FOCE-I. For the permutation test, the type I error estimates of the Wald test and the LRT are non-significantly different from the nominal level of 5%, with both algorithms. Of note with FOCE-I, the maximal value observed over the 195 thresholds is four times larger than that obtained with SAEM for the Wald test. Using the simulation-based correction or the permutation test, the corrected power estimates for the Wald test are much lower using FOCE-I than SAEM [Figure 3 about here.] As displayed on Fig. 3, all the Kolmogorov distance estimates are below the K-S threshold for the permutation test and the simulation based correction. With SAEM and FOCE-I,the Kolmogorov distance estimate of the DF G approach is similar to that obtained using permutation test or simulations based correction with both estimation algorithms, although the type I error estimate is about 10% with FOCE-I and not significantly different from the nominal level of 5% with SAEM. For the asymptotic test and the DF PB approach, the Kolmogorov distance and type I error estimates are close and significantly inflated. Indeed, the large denominator df of the DF PB approach leads to a threshold quite close to that of a χ 2 divided by 2. The DF FC and the DF W methods also have close and significantly inflated estimates of Kolmogorov distance and type I error, though lower than the asymptotic test.

  concluded that when unable to guarantee, in advance, the reliability of p-values based on the F -distribution one should use a permutation test. However, as far as we know, such an evaluation has not yet been published in the specific context of NLMEM. Both simulation and permutation based-approaches are relatively time consuming with computing time estimates which are about 3.4 h and 24 min compared to 19.2 s and 12.4 s for the asymptotic Wald test with SAEM and FOCE-I, respectively. The permutation test is found to be a robust alternative to correction based on simulations under H 0 from a model, with type I error close to 5%. With regard to the main assumption of permutation tests, our simulation setting ensures exchangeability of the observations as the genetic covariate only affected the fixed effects. The external validity of a permutation test is often questioned as it is by construction 'sample dependent'. Manly (1998) argues however that it is equivalent to question the representativeness of the study sample which also conditions the external validity of classical tests.

Figure 1 .Figure 2 .Figure 3 .Figure 4 .

 1234 Figure 1. Indinavir pharmacokinetic concentrations versus time profiles for the 40 patients included in the COPHAR2-ANRS 111 study sorted by dose and genotype for the CYP 3A4*1B polymorphism. The plain, dashed and dotted lines represent patients with a 400, 600 and 800 mg dose, respectively.

  of data sets on which the test could be performed.Predicted interval for a nominal level of 5% = [2.0 -8.0].

Table 1

 1 Threshold, denominator degrees of freedom (df ), type I error, and power estimates (%) of the Wald test and the LRT for different approaches, using both FOCE-I and SAEM

	Algorithm Test Approach	denominator	Threshold	T Type I T Power
			df		error
	FOCE	Wald asymptotic	-	5.99 *	199 13.1 199 53.8
		simulation	-	10.18 *	199 4.5 199 18.6
		permutation	-	9.09[6.9;34.65] * † 195 4.6 197 27.4
		DF P B	117	3.07	199 12.1 199 52.3
		DF W	37	3.25	199 11.1 199 48.2
		DF G	37	3.25	199 10.1 199 44.2
		LRT asymptotic	-	5.99 *	200 8.0 199 76.4
		simulation	-	7.05 *	200 5.5 199 69.8
		permutation	-	6.57[5.75;7.63] * † 198 6.5 199 72.4
	SAEM	Wald asymptotic	-	5.99 *	200 10.0 200 78.0
		simulation	-	7.80 *	200 5.5 200 71.0
		permutation	-	7.33[6.33;8.46] * † 200 6.0 200 73.0
		DF P B

Acknowledgements

The authors would like to thank IFR02 and Hervé Le Nagard for the use of the Centre de Biomodélisation as well as Prof. Marc Lavielle for the valuable help he provided in using MONOLIX.

During this research, Julie Bertrand was supported by a grant from the Institut de Recherches Internationales Servier (France).


Mathieu Lagrange
email: mathieu.lagrange@ircam.fr

Alexey Ozerov
email: alexey.ozerov@technicolor.com

Emmanuel Vincent
email: emmanuel.vincent@inria.fr

ROBUST SINGER IDENTIFICATION IN POLYPHONIC MUSIC USING MELODY ENHANCEMENT AND UNCERTAINTY-BASED LEARNING

Enhancing specific parts of a polyphonic music signal is believed to be a promising way of breaking the glass ceiling that most Music Information Retrieval (MIR) systems are now facing. The use of signal enhancement as a pre-processing step has led to limited improvement though, because distortions inevitably remain in the enhanced signals that may propagate to the subsequent feature extraction and classification stages. Previous studies attempting to reduce the impact of these distortions have relied on the use of feature weighting or missing feature theory. Based on advances in the field of noise-robust speech recognition, we represent the uncertainty about the enhanced signals via a Gaussian distribution instead that is subsequently propagated to the features and to the classifier. We introduce new methods to estimate the uncertainty from the signal in a fully automatic manner and to learn the classifier directly from polyphonic data. We illustrate the results by considering the task of identifying, from a given set of singers, which one is singing at a given time in a given song. Experimental results demonstrate the relevance of our approach.

INTRODUCTION

Being able to focus on specific parts of a polyphonic musical signal is believed to be a promising way of breaking the glass ceiling that most Music Information Retrieval (MIR) tasks are now facing [START_REF] Aucouturier | Improving timbre similarity: How high is the sky?[END_REF]. Many approaches were recently proposed to enhance specific signals (e.g., vocals, drums, bass) by means of source separation methods [START_REF] Durrieu | Source/filter model for unsupervised main melody extraction from polyphonic audio signals[END_REF][START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF].

The benefit of signal enhancement has already been proven for several MIR classification tasks, such as singer identification [START_REF] Fujihara | A modeling of singing voice robust to accompaniment sounds and its application to singer music information retrieval[END_REF][START_REF] Mesaros | Singer identification in polyphonic music using vocal separation and pattern recognition methods[END_REF], instrument recognition [START_REF] Heittola | Musical instrument recognition in polyphonic audio using sourcefilter model for sound separation[END_REF], tempo estimation [START_REF] Chordia | Using source separation to improve tempo detection[END_REF], and chord recognition [START_REF] Reed | Minimum classification error training to improve isolated chord recognition[END_REF]. In most of those works, signal enhancement was used as a pre-processing step. Since the enhancement process must operate with limited prior knowledge about the properties of the specific parts to be enhanced, distortions inevitably remain in the enhanced signals that propagate to the subsequent fea-Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. c 2012 International Society for Music Information Retrieval. ture extraction and classification stages resulting in limited improvement or even degradation of the classification accuracy.

A few studies have attempted to reduce the impact of these distortions on the classification accuracy. In [START_REF] Fujihara | A modeling of singing voice robust to accompaniment sounds and its application to singer music information retrieval[END_REF][START_REF] Kitahara | Instrument identification in polyphonic music: Feature weighting to minimize influence of sound overlaps[END_REF], feature weighting and frame selection techniques were proposed that associate a constant reliability weight to each feature over all time frames or to all features in each time frame. In practice, however, distortions affect different features in different time frames so that the assumption of constant reliability does not hold. A more powerful approach consists of estimating and exploiting the reliability of each feature within each time frame. A first step in this direction was taken in [START_REF] Eggink | Application of missing feature theory to the recognition of musical instruments in polyphonic audio[END_REF], where recognition of musical instruments in polyphonic audio was achieved using the missing feature theory. This theory adopted from noise-robust speech recognition assumes that only certain features are observed in each time frame while other features are missing and thus discarded from the classification process [START_REF] Cooke | Robust automatic speech recognition with missing and unreliable acoustic data[END_REF].

Nevertheless, the approach in [START_REF] Eggink | Application of missing feature theory to the recognition of musical instruments in polyphonic audio[END_REF] has the following three limitations. First, such binary uncertainty (either observed or missing) does not account for partially distorted features nor for correlations between the distortions affecting different features. To avoid this limitation, it was proposed in the speech recognition field to use the so-called Gaussian uncertainty [START_REF] Deng | Dynamic compensation of HMM variances using the feature enhancement uncertainty computed from a parametric model of speech distortion[END_REF], where the distortions over a feature vector are modeled as a zero-mean multivariate Gaussian with possibly non-diagonal covariance matrix. Second, this approach necessitates clean data to train the classifiers, while for some tasks, e.g., singer identification, collecting such clean data may be impossible. Third, the approach in [START_REF] Eggink | Application of missing feature theory to the recognition of musical instruments in polyphonic audio[END_REF] relies on manual f0 annotation and its use in a fully automatic system has not been demonstrated.

The contribution of this paper is threefold: (1) promoting the use of Gaussian uncertainty instead of binary uncertainty for robust classification in the field of MIR, (2) using a fully automatic procedure for Gaussian uncertainty estimation, (3) learning classifiers directly from noisy data with Gaussian uncertainty.

To illustrate the potential of the proposed approach we consider in this paper the task of singer identification in popular music and address it, in line with [START_REF] Fujihara | A modeling of singing voice robust to accompaniment sounds and its application to singer music information retrieval[END_REF][START_REF] Mesaros | Singer identification in polyphonic music using vocal separation and pattern recognition methods[END_REF], using Gaussian Mixture Model (GMM)-based classifiers and Mel-frequency cepstral coefficients (MFCCs) as features. We consider this task since it is one of the MIR classification tasks for which the benefit of signal enhancement is most obvious. Indeed, the information about singer identity is mostly concentrated in the singing voice signal.

The remainder of this paper is organized as follows. Some background about singer identification and baseline approaches is provided in Section 2. The proposed approach based on Gaussian uncertainty is detailed in Section 3. Experiments are presented in Section 4 and a discussion is provided in Section 5.

SINGER IDENTIFICATION

Background

When it comes to characterizing a song from its content, identifying the singer that is performing at a given time in a given song is arguably an interesting and useful piece of information. Indeed, most listeners have a strong commitment to the singer while listening to a given song. However, the literature about automatic singer identification is relatively scarce, compared for example with musical genre detection. This may be explained by several difficulties that pose interesting challenges for research in machine listening.

First, the human voice is a very flexible and versatile instrument and very small changes in its properties have noticable effects on human perception. Second, the musical accompaniment that forms the background is very diverse and operates at about the same loudness as the singing voice. Hence, very little can be assumed on both sides and the influence of the background cannot be neglected.

For humans, though, it is relatively easy to focus on the melody sung by the singer as our hearing system is highly skilled at segregating human vocalizations within cluttered acoustical environments. This segregation is also made possible by compositional choices. For example, most of the time in pop music, only one singer is singing at a time, and if not, the others are background vocals that are usually more easily predictable and sung at a relatively low volume.

From an application perspective, singing voice enhancement is expected to be useful for the identification of singers which have sung with different bands or with different instrumentations, such as unplugged versions. More on the so-called album effect can be found in [START_REF] Kim | Towards quantifying the album effect in artist identification[END_REF]. In this case, classifying the mixture signal will induce variability in the singer models due to occlusion, while classifying the singing voice signal alone should provide better identification. The same remark applies to the case where a song features multiple singers and one needs to identify which singer is singing at a given time. For some other repertoires where the notions of singer and artist/band are very tightly linked, it is questionable whether the singing voice signal suffices for classification, because the musical background can also provide discriminative cues. Nevertheless, singing voice enhancement is likely to remain beneficial by enabling the computation of separate features over the singing voice and over the background and their fusion in the classification process. In this paper, for simplicity, we illustrate the potential of our approach by considering the singing voice signal only unless otherwise stated.

Baseline Approaches

More formally, let us assume that each recording x f n (also called mixture), represented here directly in the Short Term Fourier Transform (STFT) domain, f = 1, . . . , F and n = 1, . . . , N being respectively frequency and time indices, is the sum of two contributions: the main melody (here the singing voice) v f n and the accompaniment a f n . This can be written in the following vector form:

x n = v n + a n , (1)
where

x n = [x 1n , . . . , x F n ] T , v n = [v 1n , . . . , v F n ] T and a n = [a 1n , . . . , a F n ] T .
We assume that there are K singers to be recognized, and for each singer there is a sufficient amount of training and testing mixtures. In line with [START_REF] Fujihara | A modeling of singing voice robust to accompaniment sounds and its application to singer music information retrieval[END_REF][START_REF] Mesaros | Singer identification in polyphonic music using vocal separation and pattern recognition methods[END_REF], we adopt a singer identification approach based on MFCC features and GMMs.

Without any melody enhancement such an approach consists in the following two steps [START_REF] Kim | Singer identification in popular music recordings using voice coding features[END_REF] (Fig. 1):

1. Learning: For each singer k = 1, . . . , K, the corresponding GMM model is estimated in the maximum likelihood (ML) sense from the features (here MFCCs) ȳ computed directly from the training mixtures of that singer.

Decoding:

A testing mixture x is assigned to the singer k for which the likelihood of model θ k evaluated on the features extracted in the same way is maximum 1 . In order to gain invariance with respect to the accompaniment, one needs to separate the contribution of the accompaniment and the singer within the mixture. This separation may be embedded within the classifier, as in [START_REF] Tsai | Automatic singer recognition of popular music recordings via estimation and modeling of solo vocal signals[END_REF]. In this case, the separation has to be performed in the feature domain, usually the log Mel spectrum.

Alternatively, melody enhancement can be applied as a pre-processing step [START_REF] Fujihara | A modeling of singing voice robust to accompaniment sounds and its application to singer music information retrieval[END_REF][START_REF] Mesaros | Singer identification in polyphonic music using vocal separation and pattern recognition methods[END_REF] over the spectrogram of the mixture. since the spectrogram have better spectral resolution than the log Mel spectrum, this approach can potentially achieve better discrimination, as in that case, the features (MFCCs) are no longer computed from the audio mixture, but from the corresponding melody estimate v (Fig. 2).

PROPOSED APPROACH

Inspired by some approaches in speech processing [START_REF] Deng | Dynamic compensation of HMM variances using the feature enhancement uncertainty computed from a parametric model of speech distortion[END_REF], we propose to consider Gaussian uncertainty by augmenting the melody estimates v by a set of covariance matrices Σ v representing the errors about these estimates. This Gaussian uncertainty is first estimated in the STFT domain, then propagated through MFCC computation, and finally exploited for GMM learning and decoding steps (Fig. 3).

Melody Enhancement

Given the mixture, we assume that each STFT frame v n of the melody is distributed as

v n |x n ∼ N (v n , Σv,n ), (2)
and we are looking for an estimate of vn and Σv,n . In this study, we have chosen the melody enhancement method 2 proposed by Durrieu et al. [START_REF] Durrieu | Source/filter model for unsupervised main melody extraction from polyphonic audio signals[END_REF]. This method has shown very promising results for vocals enhancement task within the 2011 Signal Separation Evaluation Campaign (SiSEC 2011) [START_REF] Araki | The 2011 signal separation evaluation campaign (SiSEC2011): -audio source separation[END_REF] and its underlying probabilistic model facilitates STFT domain uncertainty computation.

The main melody v, usually a singer, is modeled thanks to a source/filter model, and the accompaniment a is modeled using Non-negative Matrix Factorization (NMF) model. The leading voice is assumed to be harmonic and monophonic. The separation system mainly tracks the leading voice following two cues: first its energy, and second the smoothness of the melody line. Therefore, the resulting separated leading voice is usually the instrument or voice that is the most salient in the mixture, over certain durations of the signal. Overall this modeling falls into the framework of constrained hierarchical NMF with Itakura-Saito divergence [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF], which allows a probabilistic Gaussian interpretation [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF].

More precisely the method is designed for stereo mixtures. Let mixing equation

x j,f n = v j,f n + a j,f n (3)
be a stereophonic version of the monophonic mixing equation [START_REF] Adiloglu | An uncertainty estimation approach for the extraction of individual source features in multisource recordings[END_REF], where j = 1, 2 is the channel index and equations ( 1) and ( 3) are related for any signal s j,f n as A probabilistic Gaussian interpretation of modeling in [START_REF] Durrieu | Source/filter model for unsupervised main melody extraction from polyphonic audio signals[END_REF] assumes v j,f n and a j,f n are zero-mean Gaussians that are mutually independent and independent over channel j, frequency f and time n. The corresponding constrained hierarchical NMF structured modeling allows the estimation of their respective variances σ 2 v,j,f n and σ 2 a,j,f n from the multichannel mixture. With these assumptions the posterior distribution of v j,f n given x j,f n can be shown to be Gaussian with mean

s f n = (s 1,f n + s 2,f n )/2. (4)
vj,fn = σ 2 v,j,f n σ 2 v,j,f n + σ 2 a,j,f n x j,f n (5)
obtained by Wiener filtering, as in [START_REF] Durrieu | Source/filter model for unsupervised main melody extraction from polyphonic audio signals[END_REF], and the variance [START_REF] Rose | Integrated models of signal and background with application to speaker identification in noise[END_REF]

σ2 v,j,f n = σ 2 v,j,f n σ 2 a,j,f n σ 2 v,j,f n + σ 2 a,j,f n . (6)
Finally, thanks to the posterior between-channel independence of v j,f n and the down-mixing (4), vn and Σv,n in (2) are computed as

vn = (v 1,f n + v2,fn )/2 f , (7)
Σv,n = diag (σ 2 v,1,f n + σ2 v,2,f n )/2 f . (8)
Note that any Gaussian model-based signal enhancement method, e.g., one of the methods implementable via the general source separation framework in [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF], is suitable to compute this kind of uncertainty in the time-frequency domain.

Uncertainty Propagation during MFCC Computation

Let M(•) be the nonlinear transform used to compute an M -dimensional MFCC feature vector y n ∈ R M . It can be expressed as [START_REF] Adiloglu | An uncertainty estimation approach for the extraction of individual source features in multisource recordings[END_REF]

y n = M(v n ) = D log(M|v n |), (9)
where D is the M ×M DCT matrix, M is the M ×F matrix containing the Mel filter coefficients, and | • | and log(•) are both element-wise operations.

In line with (2), we assume that the clean (missing) feature y n = M(v n ) is distributed as

y n |x n ∼ N (ȳ n , Σy,n ), (10)
which is an approximation because of the Gaussian assumption ( 2) and the nonlinear nature of M(•).

To compute the feature estimate ȳn and its Gaussian uncertainty covariance Σy,n we propose to use the Vector Taylor Series (VTS) method [START_REF] Moreno | A vector Taylor series approach for environment-independent speech recognition[END_REF] that consists in linearizing the transform M(•) by its first-order vector Taylor expansion in the neighborhood of the voice estimate vn : [START_REF] Goto | RWC music database: popular, classical, and jazz music databases[END_REF] where J M (v n ) is the Jacobian matrix of M(v n ) computed in v n = vn . This leads to the following estimates of the noisy feature value ȳn and its uncertainty covariance Σy,n [START_REF] Fujihara | A modeling of singing voice robust to accompaniment sounds and its application to singer music information retrieval[END_REF], as propagated through this (now linear) transform:

y n = M(v n ) ≈ M(v n ) + J M (v n ) (v n -vn ),
ȳn = M(v n ), ( 12
)
Σy,n = D M M|v n |1 1×F Σv,n D M M|v n |1 1×F T , ( 13
)
where 1 1×F is an 1 × F vector of ones and the magnitude | • | and the division are both element-wise operations.

GMM Decoding and Learning with Uncertainty

Each singer is modeled by a GMM θ = {µ i , Σ i , ω i } I i=1 , where i = 1, . . . , I are mixture component indices, and µ i , Σ i and ω i ( i ω i = 1) are respectively the mean, the covariance matrix and the weight of the i-th component. In other words, each clean feature vector y n is modeled as follows:

p(y n |θ) = I i=1 ω i N (y n |µ i , Σ i ), (14)
where

N (y n |µ i , Σ i ) 1 (2π) M |Σ i | - (y n -µ i ) T Σ -1 i (y n -µ i ) 2 . ( 15
)
Since the clean feature sequence y = {y n } n is not observed, its likelihood, given model θ, cannot be computed using [START_REF] Kim | Towards quantifying the album effect in artist identification[END_REF]. Thus in the "likelihood computation" step (Fig. 3), we rather compute the likelihood of the noisy features ȳ given the uncertainty and the model, that can be shown to be equal to [START_REF] Deng | Dynamic compensation of HMM variances using the feature enhancement uncertainty computed from a parametric model of speech distortion[END_REF]: [START_REF] Mesaros | Singer identification in polyphonic music using vocal separation and pattern recognition methods[END_REF] We see that in this likelihood Gaussian uncertainty covariance Σy,n adds to the prior GMM covariance Σ i , thus adaptively decreasing the effect of signal distortion.

p(ȳ| Σy , θ) = N n=1 I i=1 ω i N (ȳ n |µ i , Σ i + Σy,n ).
In the "model learning" step (Fig. 3), we propose to estimate the GMM parameters θ by maximizing the likelihood [START_REF] Mesaros | Singer identification in polyphonic music using vocal separation and pattern recognition methods[END_REF]. This can be achieved via the iterative Expectation-Maximization (EM) algorithm introduced in [START_REF] Ozerov | GMM-based classification from noisy features[END_REF] and summarized in Algorithm 1. The derivation of this algorithm is omitted here due to lack of space and the Matlab source code for GMM decoding and learning is available at http: //bass-db.gforge.inria.fr/amulet. Algorithm 1 One iteration of the EM algorithm for the likelihood integration-based GMM learning from noisy data. E step. Conditional expectations of natural statistics:

γ i,n ∝ ω i N (ȳ n |µ i , Σ i + Σy,n ), and i γ i,n = 1, (17) ŷi,n = W i,n (ȳ n -µ i ) + µ i , (18)
R yy,i,n = ŷi,n ŷT i,n + (I -W i,n ) Σ i , (19)
where

W i,n = Σ i Σ i + Σy,n -1 . (20)
M step. Update GMM parameters:

ω i = 1 N N n=1 γ i,n , (21)
µ i = 1 N n=1 γ i,n N n=1 γ i,n ŷi,n , (22)
Σ i = 1 N n=1 γ i,n N n=1 γ i,n R yy,i,n -µ i µ T i . (23) 4. EXPERIMENTS 4.1 Database
For our evaluation, we consider a subset of the RWC Popular Music Database [START_REF] Goto | RWC music database: popular, classical, and jazz music databases[END_REF] which has previously been considered in [START_REF] Fujihara | A modeling of singing voice robust to accompaniment sounds and its application to singer music information retrieval[END_REF] for the same task. It consists of 40 songs sung by 10 singers, five of which were male (denoted by a to e) and the five others female (denoted by f to j). This set is then divided into the four groups of songs considered in [START_REF] Fujihara | A modeling of singing voice robust to accompaniment sounds and its application to singer music information retrieval[END_REF], each containing one song by each singer . Each of those songs is then split into segments of 10 seconds duration. Among those segments, only the ones where a singing voice is present (not necessarily during the whole duration of the segment) are kept unless otherwise stated.

Considering short duration segments instead of the whole song is done for two reasons. First, it makes the task more generic in the sense that multiple singers can also potentially be tracked within a same song. Second, it allows us to gain statistical relevance during the cross validation by enlarging the number of tests.

Methods

For each of those segments, features are computed and classified using the three methods depicted in Figures 1 to 3. The first one, named mix, consists in computing the features directly from the mixture, and serves as a baseline. The second method, termed v-sep, consider melody enhancement as a pre-processing step. The main melody enhancement system considered in this study is available under two versions: a version focusing on the voiced part of the singing voice and another version attempting to jointly enhance the voiced and the unvoiced parts of the singing voice (see [START_REF] Durrieu | Source/filter model for unsupervised main melody extraction from polyphonic audio signals[END_REF] for details). In the following, only the results of the former are reported since the latter led to much smaller classification accuracy. When the estimated vocals signal has zero power in a given time frame, the resulting MFCCs may be undefined. Such frames are discarded. The last method, termed v-sep-uncrt, consists in exploiting the estimated uncertainty about the enhancement process. For all those methods, we considered MFCC features and dropped the first coefficient, thus discarding energy information. Mixtures of 32 Gaussians are then trained using 50 iterations of the EM algorithm for each singer. For testing, the likelihood of each singer model is computed for each segment and the one with the highest likelihood is selected as the estimate.

Results

The aforementioned four groups of songs are considered for a 4-fold cross validation. For each fold, the selected group is used for testing and the data of the three remaining ones are used for training the models. The average detection accuracy are shown in Table 1. Compared to the baseline, v-sep and v-sep-uncrt achieve better performance while considering segments, indicating that focusing on the main harmonic source within the segment is beneficial for identifying the singer. That is, the level of feature invariance gained by the separation process more than compensates for the distortions it induces.

Considering the uncertainty estimate adds a significant level of improvement in the v -sep case. We assume that this gain of performance is obtained because the use of uncertainty allows us to focus on the energy within the spectrogram that effectively belongs to the voice and that the use of the uncertainty allows us to robustly consider standard features (MFCCs).

Performing a majority vote over the all the segments (in this case the likelihood of each singer is taken into account even if no singing voice is present) of each song gives an accuracy of 85% and restricting the vote to only the sung segments gives a 94% accuracy. These numbers can respectively be considered as worst and best cases. It is therefore likely that a complete system that would incorporate a music model to discard segments with only music would achieve an accuracy that is between those bounds. Although a more formal comparison would be needed, we believe that those results compare favorably with the performances obtained in [START_REF] Fujihara | A modeling of singing voice robust to accompaniment sounds and its application to singer music information retrieval[END_REF] using specialized features on the same dataset while standard MFCC features were used here. It is also interesting to notice that in this case of songlevel decisions, considering the separation without uncertainty does not give any improvement compared to the mix baseline.

DISCUSSION

We have presented in this paper a computational scheme for extracting meaningful information in order to tackle a music retrieval task: singer identification. This is done by considering an enhanced version of the main melody that is more or less reliable in specific regions of the time/frequency plane. Instead of blindly making use of this estimate, we propose in this paper to consider how uncertain the separation estimate is during the modeling phase. This allows us to give more or less importance to the features depending on how reliable they are in different time frames, both during the training and the testing phases. For that purpose, we adopted the Gaussian uncertainty framework and introduced new methods to estimate the uncertainty in a fully automatic manner and to learn GMM classifiers directly from polyphonic data.

One should notice that the proposed scheme is not tied to the task considered in this paper. It is in fact completely generic and may be easily applied to other GMM-based MIR classification tasks where the prior isolation of a specific part of the music signal could be beneficial. The only part that would require adaptation is the derivation of VTS uncertainty propagation equations for other features than MFCCs. Uncertainty handling for other classifiers than GMM has also received some interest recently in the speech processing community.

The experiments reported in this paper provide us with encouraging results. Concerning this specific task of singer identification, we intend to exploit both the enhanced singing voice and accompaniment signals and to experiment on other datasets with a wider range of musical styles. In particular, we believe that the hip-hop/rap musicals genres would be an excellent testbed both from a methodological and application point of view, as many songs feature several singers: knowing which singer is performing at a given time is a useful piece of information. Finally, we would like to consider other content based retrieval tasks in order to study the relevance of this scheme for a wider range of applications.

Figure 1 .

 1 Figure 1. The standard classification scheme.

Figure 2 .

 2 Figure 2. Considering melody enhancement as a preprocessing step.

Figure 3 .

 3 Figure 3. Proposed approach with melody enhancement and Gaussian uncertainty.

Table 1 .

 1 Average accuracy of the tested methods per singing segment and per song, considering either all the segments or only those segments where singing voice is present in the latter case.

	Accuracy (%)		per 10 sec. singing segment			per song
	input	Fold 1 Fold 2 Fold 3 Fold 4 Total all seg. sung seg.
	mix	51	53	55	38	49	57	64
	v-sep	60	63	53	43	55	57	64
	v-sep-uncrt	71	72	84	83	77	85	94

In order not to overload the notations, the singer index k is omitted hereafter, where applicable.

The Python source code is available at http://www.durrieu. ch/research/jstsp2010.html

ACKNOWLEDGMENTS

This work was partly supported by the Quaero project funded by Oseo and ANR-11-JS03-005-01. The authors wish to thank Jean-Louis Durrieu for kindly providing his melody estimation source code to the community and Mathias Rossignol for his useful comments.