Validating the association between brain activity, as measured in functional MRI, with a combination or a contrast of tasks is usually performed by replicating an experiment in a small group of subjects, and by assessing the presence of a statistically significant average effect across subjects (random effects analyses). While many efforts have been made to control the rate of false detections, statistical characteristics of the data have rarely been studied, and the reliability of the results (supra-thresholds areas that are considered as activated regions) has rarely been assessed. In this work, we take advantage of the large cohort of subjects who underwent the Localizer experiment to study the statistical nature of group data, propose some measures of the reliability of group studies, and address simple methodological questions as : is there, from the point of view of reliability, an optimal statistical threshold for activity maps ? How many subjects should be included in group studies ? What method should be preferred for inference ? Our results suggest that i) optimal thresholds can indeed be found, and are rather lower than usual corrected for multiple comparison thresholds ii) 20 subjects or more should be included in functional neuroimaging studies in order to have sufficient reliability, iii) non-parametric significance assessment should be preferred to parametric methods iv) cluster-level thresholding is more reliable than voxel-based thresholding v) mixed effects tests are much more reliable than random effects tests. Moreover, our study shows that inter-subject variability plays a prominent role in the relatively low sensitivity and reliability of group studies.

Introduction 1.Inter-subject variability in neuroimaging

Many scientists spend lots of efforts in order to obtain statistically significant results in neuroimaging studies, in order to validate a prior hypothesis on brain function, and it is certainly true that one of the greatest difficulties that they have to face is the high variability that is present in their datasets across subjects. In particular, this high degree of variability dramatically reduces the sensitivity of random effects studies. For instance, in [START_REF] Wei | Functional MRI of auditory verbal working memory: long-term reproducibility analysis[END_REF], it has been shown empirically that inter-session variability was much higher than inter-subject variability.

Several techniques have been employed to quantify inter-subject differences [START_REF] Kherif | Group analysis in functional neuroimaging: selecting subjects using similarity measures[END_REF]. The use of such tools shows that inhomogeneous populations can be met quite frequently in neuroimaging studies. This problem is dramatic in general, given the small number of subjects in many studies (10-15 subjects, sometimes less).

Analysing the variability that arises between different datasets is not straightforward, since it is probably of a mixture of random and deterministic or structured factors. The first one is inter-session variability, which is caused by different amount of motion, attention or tiredness of the subjects. However, several specific factors can be proposed in the inter-subject setting :

• Spatial mismatch between datasets. It is known that perfect correspondences between two anatomical images cannot be found in general, and that correspondences should generally be considered as approximative, even after rigid or non-rigid spatial normalization. The magnitude order of such local shifts is probably as large as 1cm in many instances (this can be observed for functional regions like the the motor cortex or the visual areas [Thirion et al., ress, Stiers et al., 2006] or the position of anatomical landmarks [START_REF] Collins | Non-linear cerebral registration with sulcal constraints[END_REF], Hellier et al., 2003]). Across subjects, this effect typically yields a structured but variable pattern.

• Next, the magnitude of the BOLD signal recorded at the same location for several tasks is variable across subjects, and sometimes across sessions [START_REF] Smith | Variability in fMRI: a re-examination of inter-session differences[END_REF], and the precise nature of this variability is not clear. It should be reminded that fMRI is not a quantitative neuroimaging modality, and that the standard use of reporting percent of signal increase is also problematic, due to the ambiguous definition (voxel-based or global average) of the baseline reference.

• Finally, there could be global differences in the cascade of neural activation elicited by a given task, related to genetic or epigenetic differences across subjects, or to different cognitive strategies (for non trivial tasks). This is of course interesting in itself, but at the same time very difficult to account for in traditional studies, and in many instances very difficult to prove.

This is a very sketchy account of possible sources of variability across subjects, but one should remember that all these effects equally treated as confounds and globally modeled as second-level variability [START_REF] Friston | Classical and bayesian inference in neuroimaging: Theory[END_REF] terms in current linear models. It follows that voxel-based random effects analyses, that assess the significance of an effect by comparing its mean value to its variability across subjects are typically very weakly sensitive [Friston et al., 1999, McNamee andLazar, 2004]. Moreover the lack of reliability makes brain mapping procedures [START_REF] Jernigan | More "mapping" in brain mapping: statistical comparison of effects[END_REF] problematic. We illustrate this issue in Fig. 1 with an example taken from the dataset presented next.

Methodology in group analyses

Voxel-based random effects analysis is the standard way to analyse data from group studies (although the extraction of discrete local maxima [Worsley, 2005] presents an attractive alternative). It is in fact a particular instance of mixed effects models, which recently gained in popularity [START_REF] Worsley | A general statistical analysis for fMRI data[END_REF], Friston et al., 2002, Beckmann et al., 2003, Neumann and Lohmann, 2003, Woolrich et al., 2004]. However, there remain some alternatives, i.e. non-systematic parts, in the group data processing strategy.

• Choice of the threshold. Usually, statistical parametric maps (SPMs) are thresholded in such a way the the number or the rate of false positive voxels is controlled, at least asymptotically. In particular, the chosen threshold does not reflect a trade-off between the necessity of controlling both the number of false positive and the number of false negative. One straightforward reason is that it is relatively simple to model the statistical distribution under the null hypothesis, but not under the converse hypothesis.

• It is possible to threshold SPMs at the voxel or at the cluster level. In the latter case, a double thresholding procedure is used : first the map is thresholded at a (relatively lenient) significance level ; second, the size of the resulting connected components is assessed against its distribution under the null hypothesis. This is usually considered as a safe procedure, but fully neglects the possibility of small yet significant activation foci. One of the reasons is that intensive smoothing of the data simply removes such eventuality. • Another alternative consists in the choice of parametric versus non-parametric tests. While parametric tests are particularly efficient and computationally cheap, they are based on possibly unrealistic hypotheses that may reduce their sensitivity. For instance, assessing a Student statistic obtained in a group analysis with respect the Student distribution with the appropriate number of degrees of freedom assumes that the data is normally distributed under the null hypothesis. This hypothesis cannot be checked in the usual, small datasets. Non-parametric tests may avoid these issues [START_REF] Holmes | Nonparametric analysis of statistic images from functional mapping experiments[END_REF], Brammer et al., 1997, Bullmore et al., 1999, Nichols and Holmes, 2002, Hayasaka and Nichols, 2003, Mériaux et al., 2006a], but at a higher computational cost.

• Finally, many other alternatives concern the different pre-processing steps that are standard in typical studies. For simplicity, let us only retain the choice of a smoothing kernel (typically 8 to 12 mm FWHM). It is has been shown [START_REF] Shaw | Evaluating subject specific preprocessing choices in multisubject fMRI data sets using data-driven performance metrics[END_REF], LaConte et al., 2003] that cross-validation schemes could help to optimize these pre-processing choices.

Some alternatives may be proposed to assess the significance of activity in group studies. For instance, parcellation, and thus parcel-based random effects maps can be used [Thirion et al., ress], with a possibly double advantage : if parcels adapt to individual anatomy, they can cope with some parts of the inter-subject variability ; second, this considerably alleviates the multiple comparison problem. Care must however be taken when controlling the test specificity. In general, making statistical inference at the voxel level is a particular, and not necessary optimal, choice.

However, another fundamental question has not received as much attention so far : it concerns the number of subjects that should be included in the study. This number typically represents a trade-off between (a) the cost of making neuroimaging experiments on large cohort of subjects and (b) the necessity to have enough subject for the significance of statistical tests [START_REF] Desmond | [END_REF]Glover, 2002, Murphy and[START_REF] Murphy | [END_REF]. However, the true question could or should me : how many subjects are enough to make the analysis reliable, in terms of avoiding false positives as well as false negatives ? This questions has been raised in several studies, but reliability analysis has to our knowledge, not been used in this context.

Reliability can be thought of as the confidence level that can be given to a certain result. Following [START_REF] Liou | Bridging functional MR images and scientific inference: Reproducibility maps[END_REF] we define it practically as the agreement between independent measurements of a certain effect.

Taking advantage of a large dataset

In this work, we use the large dataset collected in a localizer experiment (see the companion paper) to address empirically some of the questions listed above. First, we perform a relatively elementary statistical description of the dataset. In particular we focus on the deviation from normality in the distribution of effects for a given contrast, and on the spatial distribution of second-level variance in the dataset.

Then we concentrate on the reliability of the results when varying different features in the statistical procedure. Reliability analysis is based on binary i.e. thresholded maps obtained from distinct subgroups of subjects. One measurement consists in modelling the activated/inactivated status of the voxels in the image as a mixture of binary distributions, and to assess the reproducibility of the activated/inactivated status using this model. The second measure is based on the distance between clusters of a certain size in the thresholded binary maps. A large distance means that no correspondence can be found between supra-threshold clusters across groups of subjects, hence that the maps are not very reliable. We use this framework to study the effect of the threshold value, then of the size of the group of subjects, and finally the effect of the method itself on the reliability of the maps. All this is possible due to the large sample size of the population of subjects. Importantly also, we show that our cross-validation procedure is not biased by possible confounds like a correlation between reliability and sensitivity, and that both reliability measures yield similar results.

An important point is that this procedure can be performed for different cognitive contrasts which have different characteristics in terms of contrast-to-noise ratio or spatial variability.

Materials and Methods

Dataset

The experimental paradigm is described with more detail in the companion paper. Briefly, we used an event-related experimental paradigm consisting of ten conditions. Subjects underwent a series of stimuli or were engaged in task such as passive viewing of horizontal or a vertical checkerboards, left or right click after audio or video instruction, computation (subtraction) after video or audio instruction, sentence listening, from audio or visual modality. Events were randomly occurring in time (mean inter stimulus interval: 3s), with ten occurrences per event type (except motor button clicks for which there are only five trials per session).

Eighty-one right-handed subjects participated in the study. The subjects gave informed consent and the protocol was approved by local ethics committee. Functional images were acquired on a 3T Bruker scanner using an EPI sequence (TR = 2400ms, TE = 60ms, matrix size=64 × 64, FOV = 24cm × 24cm). Each volume consisted of n a 4mm-thick axial slices without gap, where n a varied from 26 to 40 according to the session. A session comprised 130 scans. The first four functional scans were discarded to allow the MR signal to reach steady state. Anatomical T1 images were acquired on the same scanner, with a spatial resolution of 1 × 1 × 1.2 mm 3 . fMRI data processing consisted in 1) temporal Fourier interpolation to correct for between-slice timing, 2) motion estimation. For all subjects, motion estimates were smaller than 1mm and 1 degree, 3) spatial normalization of the functional images, re-interpolation to 3×3×3 mm 3 , and 4) smoothing (5mm FWHM). This pre-processing was performed with the SPM2 software (see e.g. [START_REF] Ashburner | Human Brain Function[END_REF]). Datasets were also analyzed using the SPM2 software, using standard high-pass filtering and AR(1) whitening. For further analysis, the voxel-based estimated effects for several contrasts of interest were retained.

We determined a global brain mask for the group by considering all the voxels that belong to at least half of the individual brain masks defined with SPM2. It comprises approximately 60000 voxels. Note that considering merely the intersection of the individual masks yields about 30000 voxels only.

Elementary Statistical description of the dataset

In this section, we select a few contrasts of interest, and study the statistical distribution of the corresponding parameters in each voxel. Using a first level (subject-specific) GLM, one can obtain parametric estimates of the BOLD activity at each voxel in each subject: For each subject s ∈ [1..S] and each voxel v ∈ [1..V ], we have a parameter estimate β(s, v), and a variance estimate σ(s, v) 2 .

The first question that may arise is whether the effects β(s, v) are normally distributed or not, since this is a key assumption in classical (Random Effects) group studies. We have used the test of D'Agostino-Pearson [Zar, 1999], based on the computation of the skew and the kurtosis (third and forth order cumulants) of the values { β(s, v)}, s = 1..S in each voxel v. This provides us with a P-value of the D'Agostino-Pearson statistic under the null (normal) hypothesis. For the sake of visualization, we convert the P-value into a z-value. We have then repeated the procedure based on the normalized effects {τ

(s, v) = β(s,v) σ(s,v) }, s ∈ [1..S], v ∈ [1.
.V ] which removes a potential variability in signal scaling across the population. Note that at the group level, the normalization through the residual magnitude has a much greater impact than the deviation from normality due to the fact that σ(s, v) is estimated with a finite (ν = 100) number of degrees of freedom.

Then, assuming a two-level normal model of the data

β(s, v) = β(s, v) + ε(s, v), with ε(s, v) ∼ N (0, σ(s, v) 2 ) (1) β(s, v) = β(v) + ζ(s, v), with η(s, v) ∼ N (0, v g (v) 2 ) (2)
where β(s, v) is the true effect for subject s, β(s, v) is the estimated effect for subject s, and β(v) is the average effect in the population at voxel v; ε(s, v) and η(s, v) are first-level (estimation) and second-level (inter-subject) normal residual terms. The first equation represents thus the subjectspecific estimation of the signal and the second one the group-level model. We have estimated the second level variance v g in each voxel, since it plays a central role in many group-level statistics.

In particular, an interesting question is whether β(v) and v g (v) are independent or not. Note that v g is estimated by maximizing the likelihood of the data given v g . Newton or EM estimation schemes can be used [START_REF] Worsley | A general statistical analysis for fMRI data[END_REF][START_REF] Mériaux | Robust statistics for nonparametric group analysis in fMRI[END_REF]. In this work, we use a Newton estimation scheme.

Group-level analysis methods: Voxel-based statistics

We review here different techniques used for voxel-based inter-subject activation detection. We consider a given contrast of interest. For each subject s ∈ [1..S] and each voxel v ∈ [1..V ], we have a parameter estimate β(s, v), and a variance estimate σ(s, v) 2 .

A random Effects (RFX) statistic is based on model ( 1)-( 2), in which the first level variance is neglected. It is defined as

ρ(v) = mean s∈[1..S] β(s, v) (S -1)var s∈[1..S] β(s, v) (3)
Under the null hypothesis, assuming a normal distribution for { β(s, v)}, s = 1..S, ρ(v) is Studentdistributed with (S -1) degrees of freedom, and the P-value under the null hypothesis can be assessed with or without correction for multiple comparisons1 . Alternatively, a non-parametric scheme can be used to estimate the distribution of ρ(v) under the null hypotheses, based on milder assumptions [Hayasaka andNichols, 2003, Mériaux et al., 2006a]. In this work, we use the analytical threshold.

A Mixed Effects (MFX) statistic takes into account the first-level variance: assuming a group or (second-level) variance v g (v) at each voxel v, the MFX statistic writes:

µ(v) = S s=1 β(s, v) σ(s, v) 2 + v g (v) S s=1 1 σ(s, v) 2 + v g (v) -1 2 (4)
Intuitively, MFX may perform better than RFX since it down-weights the observations with high first-level variance. The distribution of the quantity µ(v) under the null hypothesis is difficult to assess [START_REF] Woolrich | Multi-level linear modelling for fMRI group analysis using Bayesian inference[END_REF]. We rely on an non-parametric scheme as in [Mériaux et al., 2006a[START_REF] Mériaux | Robust statistics for nonparametric group analysis in fMRI[END_REF]: we tabulate the values of µ(v) for different sign swaps of each subject's dataset in order to generate a distribution under the null hypothesis, and compare the actual values with their null distribution. A quicker but very conservative approximation (µ ∼ t S-1 , t S-1 being the Student law with S -1 degrees of freedom) is also possible.

A possible alternative consists in neglecting the group variance v g in the previous formula. This yields a pseudo-MFX statistic, which we denote henceforth as ΨFX:

Ψ(v) = S s=1 β(s, v) σ(s, v) 2 S s=1 1 σ(s, v) 2 -1 2 (5)
Note that this is the statistic proposed in [START_REF] Neumann | Bayesian second-level analysis of functional magnetic resonance images[END_REF]. The difference is that we perform a frequentist test by estimating the distribution of Ψ(v) under the null hypothesis by random sign swaps ( which we refer to as non-parametric approach), exactly as for the MFX test.

As an alternative, we also use Wilkoxon's signed rank statistic (WKX) [START_REF] Hollander | Nonparametric statistical inference[END_REF], which starts with sorting the absolute effects in ascending order, then sums up the ranks modulated by the corresponding effect's sign, yielding:

W (v) = S s=1 sign( β(s, v))rank( β(s, v)) (6)
This statistic is data-independent, thus its significance is assessed very easily. Unlike the previous statistics, it does not assess the positivity of the average effect, but the asymmetry of the estimated effects β(s, v) with respect to 0, the null hypothesis being that β(s, v) are distributed symmetrically about 0. The main interest of this statistic is that it is not based on the hypothesis that the ( β(s, v)), s = 1..S are normally distributed.

Group-level analysis methods: Higher-level statistics

Higher-level or non voxel-based analyses statistical inference methods include cluster-based inference and parcel-based inference.

Cluster-based inference [Hayasaka andNichols, 2003, Mériaux et al., 2006a] is simply an extension of the voxel-based procedures, based on a double thresholding of a statistic map: first, a threshold is performed at the voxel level, then supra-threshold clusters are kept whenever their size is statistically large enough. In our implementation, we measure connectivity using the 18nearest neighbours of each voxel in 3D, and estimate the P-values at the cluster-level using the non-parametric framework.

Parcel-based inference is a different scheme in which parcels are defined across subjects using anatomical and/or functional information. Two possible schemes have been presented in [START_REF] Flandin | Improved detection sensitivity of functional MRI data using a brain parcellation technique[END_REF] and [Thirion et al., ress], based on a Gaussian Mixture Model (GMM) and a hierarchical approach respectively. A key issue of both techniques is to obtain functionally, but also spatially connected parcels. Statistics, e.g. RFX, can then be computed at the parcel level by working on parcel-based signal average instead of voxel-based signal (PRFX statistic). The advantage is that some spatial relaxation is possible in the definition of the parcels, allowing for a better spatial registration of functional information. A related difficulty arises when the same functional data is used to build the parcels and perform the test, since false positives may arise from the spatial relaxation procedure. The present work therefore aims at checking the reliability of PRFX maps.

Assessing the reliability of activation maps

In this work, we propose two measures to assess the reliability of the activation maps derived from group analysis. This first one, based on mixture of binomial distribution, characterizes the stability of the status (active/inactive) of each voxel of the dataset. The second one measures how frequently clusters of voxels are found at similar locations in the normalized MNI/Talairach space across subjects. We use these measures in a bootstrap framework that enable us to characterize the reproducibility of activation maps obtained at the group level.

Reliability measure at the voxel level

In order to estimate the reliability of a statistical model, we need a method to compare statistical maps issued from the same technique, but sampled from different groups of subjects. We use the reliability indexes elaborated in [START_REF] Genovese | Estimating testretest reliability in functional MR imaging. I: Statistical methodology[END_REF], Liou et al., 2003, Liou et al., 2005]. Assume that a statistical procedure (e.g. thresholding) yields binary maps g 1 , .., g R for different groups of subjects. At each voxel v, [g 1 (v), .., g R (v)] is thus a R-dimensional binary vector. At the image level, the distribution of G(v) = R r=1 g r (v) is modelled by a mixture of two binary distributions, one for the null hypothesis, one for the converse hypothesis: Let π 1 A be the probability that a truly active voxel is declared active, π 0 A = 1π 0 A the probability that a truly active voxel is declared inactive, π 1 I , the probability that a truly inactive voxel is declared active, π 0 I = 1π 1 I the probability that an truly inactive voxel is declared inactive, and λ the proportion of truly activated voxels. Then, using a spatial independence assumption, the log-likelihood of the data 7), with R = 8. The Y axis is in log-coordinates for the sake of readability.

writes

log(P (G)|λ, π 0 A , π 0 I ) = cst + V v=1 log λ(π 0 A ) R-G(v) (π 1 A ) G(v) + (1 -λ)(π 0 I ) R-G(v) (π 1 I ) G(v) (7)
Assuming R ≥ 3 the three free parameters, π 0 A , π 0 I , λ can be estimated using EM or Newton's methods. Note that optimizing the model over its different parameters sequentially, and using an adequate initialization, we could run the model for R = 2, though with higher variability in the estimation. An example of mixture model is given in Fig 2.

Given these estimates, the coherence index κ, known as Cohen's kappa is computed to measure the concordance of the different observations. Let p 0 = λπ 1 A + (1λ)π 0 I be the proportion of agreement between the observations and the model. p 0 should be compared to the proportion of agreement by chance p C = λπ 0 + (1λ)(1π 0 ), where π 0 = λπ 0 A + (1λ)π 0 I is the proportion of voxels declared inactive. The proportion of agreement corrected for chance is thus

κ = p 0 -p C 1 -p C (8)
In this setting, 0 ≤ κ ≤ 1 measures the accordance of the bimodal model with the data, which in turns reflects the coherence of the binary maps given as input to the model ( 7). If κ is close to 0, there is a very little agreement on which voxels are active, while there is a very good agreement if κ is close to 1. λ can also be kept as an index of the test sensitivity. Note that more complex -and realistic-models have been proposed in the literature [START_REF] Maitra | Test-retest reliability estimation of functional MRI data[END_REF], in which the parameter λ is allowed to vary spatially . However, our main purpose is not activation detection, but to obtain a global reliability measurement; for this reason, we keep the basic setting.

Reliability measure at the cluster level

Another way to assess the reliability of the results is to compare the positions of the clusters of supra-threshold voxels that arise through any group analysis. Assuming that the binary maps g 1 , .., g R are obtained from different groups of subjects through a thresholding procedure, one can post-process them in order to yield connected components. The connected components with a size greater than a given threshold η are then retained, and their centre of mass (cm) is computed: let x r i , i = 1..I(r) be the spatial coordinates of the cms derived from map g r , we propose the following average distance between any two maps:

Φ = 1 R(R -1) R r=1 s∈(1..R)-{r} 1 I(r) I(r) i=1 min j∈(1..I(s)) φ( x r i -x s j ), (9)
where

φ(x) = 1 -exp -x 2 2δ 2
is a penalty function that is close to zero when the cluster centroids are properly matched and close to 1 otherwise. Φ represents the average mismatch between the cm of a supra-threshold component in a given map and the closest cm of any supra-threshold cluster obtained from another map. Appropriate penalty terms are used to handle the case I(r) = 0. We have performed some experiments using η = 10 voxels or η = 30 voxels, and use δ = 6mm.

Procedure for the assessment of reliability

The procedure consists in dividing the population of 81 subjects in R = 2, 3, 4, 5, 6 or 8 disjoint groups of S = 40, 27, 20, 16, 13 and 10 subjects respectively. The computation of different statistics, the derivation of an adequate threshold and the thresholding are performed in the different subgroups, and global reliability measures are derived from the ensuing binary maps. This procedure is repeated 100 times for each instance, yielding a distribution of the indexes κ, λ and Φ for each possible technique/parameter.

First, we choose the traditional RFX analysis procedure and adequate parameters and evaluate the distribution of the different indexes for three contrasts of interest. This is important to understand how well the indexes are characteristic of the amount, the spread and the variability of supra-threshold activity. In particular, it is important that the estimated reliability indexes are e.g. less variable for a given contrast than across contrasts.

Second, we evaluate the choice of the threshold on the different indexes, in the case of the voxelbased RFX test. While the sensitivity index certainly decreases while the threshold increases, the behaviour of the reliability may be more complex, due to the trade-off between false positive and false negatives rates.

Third, we study the behavior of the different measurements when the number of subjects in the group varies; while it is obvious that reliability increases with the group size, it is not clear whether there exists a plateau and at which level. Previous studies [START_REF] Desmond | [END_REF]Glover, 2002, Murphy and[START_REF] Murphy | [END_REF] suggest a steady increase of sensitivity with the group size.

Finally, we choose the following statistics: RFX, RFX on smoothed (12mm FWHM instead of 5mm) effect maps (SRFX), MFX, Wilkoxon(WKX), Cluster-level RFX (CRFX), Parcel-based RFX (PRFX) and ΨFX. RFX, SRFX, MFX, ΨFX and PRFX maps are thresholded at the p < 0.001 level, uncorrected for multiple comparisons. CRFX is thresholded at p < 0.01, uncorrected level at the voxel level, then at p < 0.01, at the cluster level.

PRFX maps are computed for a number Q = 500 of parcels. Since the parcel centres are defined at the group level in Talairach space, the voxels in the group result map are assigned to the parcel with the closest center in Talairach space. This results in a piecewise constant map, the pieces resulting from a Voronoi parcellation of the group mask into parcels. Note that in our bootstrap procedure, such boundaries are defined independently in each subgroup of subject. For parcellation, we use the hierarchical procedure presented in [Thirion et al., ress]; from our experiments (not shown), it yields slightly better results than the GMM procedure [START_REF] Flandin | Improved detection sensitivity of functional MRI data using a brain parcellation technique[END_REF].

Results

Statistical model of the inter-subject data

We have performed the D'Agostino-Pearson test on the effects β(v) of all the voxels, as well as the normalized effects β σ (v). This yields a map for each contrast. We present three of them, for the contrasts left click-right click, audio instructions-video instructions and computation-reading, thresholded at the p < 0.001 uncorrected level. We also present the inter-subject variability maps v g (v) computed in the mixed-effect model (2). We present these maps together with the RFX map based on 81 subjects in Fig. 345.

In either case, the regions with highest group variance are found in the regions with highest RFX statistics in absolute values; some of them are absent in the maps 3-5, where signed statistics are presented.

The observation of these maps suggests that

• High variance areas tend to co-localize with the activated areas. This implies that the parameters v g (v) and β(v) are certainly not independent, and that statistics that are penalized by the group variance may not be very efficient in general.

• Non-normality is very significant in wide regions of the brain.

• Deviation from the normality hypothesis is much lower for the normalized effects β σ than for the raw effects β. It means dimensionless first-level statistics yield more homogeneous quantities across subjects than effects expressed in percents of baseline signal increase. Note that all the z-maps are thresholded at z = 8 for numerical reasons. The color scale of the variance image has been chosen arbitrarily in order to have supra-threshold areas that are comparable with the other maps. The variance is expressed in squared percentage of the BOLD mean signal.

• Deviation from normality of the effects does not specifically co-localize with activated areas, but, in several cases it coincides with the boundaries of activated areas.

Reliability measurements for different cognitive contrasts

We have applied an RFX analysis for different cognitive contrasts, using R = 5 groups of S = 16 and an RFX threshold θ = 3.1. The contrasts are left click-right click, audio instructions-video instructions and computation-reading. The reliability index κ, the proportion of putative true positives λ and the inter-cluster distance penalty Φ are given in Fig. 6. It shows that κ and λ have a different behaviour and are strongly dependent on the cognitive contrast under study. For instance, the left motor contrast activates relatively small regions with a relatively low reliability; the auditory-selective contrast activates larger regions, with high reproducibility ; the computationselective contrast activates larger regions, but with low reliability. The inter-cluster distance penalty Φ does not discriminate between the different contrasts as strongly as κ. As could have been expected, it has the opposite behaviour (maximal for the computation contrast, minimal for the auditory contrast).

How does the threshold affect the reliability of the analysis

Here we study the behaviour of our reliability measures when applied to a thresholded RFX map, when we let the threshold vary. The reliability measure is computed for 100 different splits of the population of subjects into R = 5 groups of S = 16 subjects, in the case of the left click-right click contrast. The threshold varies from θ = 2.2 (p < 0.015, uncorrected) to θ = 4 (p < 3.2.10 -5 ,uncorrected) in steps of 0.2. Logically, the sensitivity parameter λ decreases when θ increases (see Fig. 7(b)). More interestingly, κ reaches a maximum for θ ⋆ ∼ 2.7, but the index remains close at least for θ < 3.5 as can be seen in seen in Fig. 7(a). Accordingly, the inter-cluster distance penalty Φ is minimized for a threshold θ ⋆ ∼ 3. The correspondence of this results is interesting, given that these two similarity measures are obtained independently, and based on different considerations. Note that we have obtained similar results when studying the other contrasts with slightly higher (auditory contrast) or lower (computation contrast) threshold values. Thereafter, we retain the threshold θ = 3.1 (p < 0.001, uncorrected for multiple comparisons) for RFX analyses.

How many subjects are necessary to obtain a reliable RFX map

We study the dependence of κ, λ and Φ when we let the size S of the group vary. We base our investigation on the left click-right click contrast, with group maps thresholded at the θ = 3.1 level. The results are presented in Fig. 8. It shows that the reliability increases with the group size, which was not unexpected. The sensitivity also increases with the group size. Interestingly, the reliability reaches a plateau only for S > 20. The inter-cluster distance penalty Φ has a similar behaviour, with a plateau for S > 20 subjects when η = 10, while lower values are reached when using η = 30. : Dependence of the reproducibility and of the sensitivity of the group RFX analysis on the functional contrast under consideration. These results are obtained by drawing 5 disjoint groups of S = 16 subjects in the population of 81 subjects, and applying the procedure described in section 2.5.1. The threshold is θ = 3.1 (a) The results are more reliable for the contrast that shows auditory regions than for a contrast that shows motor activity or a contrast that shows the regions involved in the computation task. (b) By contrast the size of the putatively activated areas is greater for the contrast that shows regions involved in computation, and smaller for the contrast that shows the regions involved in motor activity. (c-d) The cluster variability penalty Φ is presented for clusters of more than η = 10 (c) or η = 30 (d) voxels. The behaviour is as expected, with the smallest value for the auditory-specific contrast.

Comparison of different group analysis methods

Now we study how the reliability index behaves for different statistical methods: The voxel-based RFX, the RFX test after 12 mm smoothing of the data -instead of 5mm-(SRFX), the MFX test, the parcel-based RFX test (PRFX), the cluster-level thresholded RFX (CRFX), The Wilkoxon test (WKX), and pseudo-MFX test ΨFX. RFX, SRFX, MFX, WKX, PRFX and ΨFX maps are thresholded at the p < 0.001 level, uncorrected for multiple comparisons. The CRFX map is first thresholded at the p < 0.01, uncorrected level, then at the p < 0.01 cluster-level. The results are obtained by bootstrapping in groups of size S = 10. The results are presented in Fig. 9.

From the point of view of reliability, the WKX and RFX tests have the worst performance overall, while the SRFX performs slightly better. CRFX, PRFX and MFX techniques yield higher reliability, but ΨFX yield the highest values. The results are more variable with PRFX than with other techniques; this reflects the fact that PRFX is based on a smaller number of volume elements, so that statistical tests have a more crispy behaviour.

CRFX, MFX, and to a lesser extent, PRFX tests are more sensitive, i.e. have a larger fraction of generally activated voxels, than voxel-based tests. Finally, the average supra-threshold cluster distance Φ is minimal for ΨFX, and relatively low for MFX. It is approximately similar for the other techniques.

Discussion

Normality and Second-level variance

From Figures 345, one of the most striking effects is the co-localization of high second-level variance areas with high RFX statistic areas. Numerically, such an effect is not expected since the RFX is defined as the quotient of the mean effect by its standard deviation. Its interpretation could be that 1) the contrast-to-noise ratio (CNR) of the BOLD effect is highly variable across subjects, and by definition this effect does not appear in non-activated areas 2) spatial mis-registration2 implies that at a given voxel, i.e. a given position in MNI space, some subjects have activity while other subjects have no activity there, thus widening the signal distribution. For simple contrasts such as those used (left or right click, sentence listening), different cognitive strategies should be ruled out. This inflated variance effect certainly deserves more investigation, given its prominent effect on statistics (sensitivity and reliability): for instance, the ΨFX statistic -that does not take into account the group variance, hence is simply a weighted average of the subject-based effects-seems more reliable than the MFX statistic, which is itself much more reliable than the RFX statistic (see Fig. 9). The effect of group variance is also an argument in favour of Bayesian analysis of fMRI data, where the reference signal level is not 0 [START_REF] Friston | Posterior probability maps and SPMs[END_REF].

Non-normality is another important factor. To our knowledge, it has not been investigated earlier, since it requires a high number of subjects. Interestingly, the importance of non-normality is minored when considering normalized effects τ (s, v) instead of raw effects β (s, v). This confirms that first-level statistics can play an import role on group statistics. Interestingly, several areas with significant non-normality are found at the periphery of activation maxima, confirming the impact of spatial shifts in group statistics. once again, further investigations on non-normality may be performed, e.g. searching different groups of subject in the population or outlier subjects (see [START_REF] Kherif | Group analysis in functional neuroimaging: selecting subjects using similarity measures[END_REF]). The effect of non-normality can be evaluated by comparing the theoretical thresholds, based on normality assumption, with those derived non-parametrically, by random sign swap of the effects across subjects. In general, it is advisable to use non-parametric assessment to obtain reliable thresholds [Mériaux et al., 2006a]. However, the choice of the statistic should not necessarily be based on the avoidance of normal noise model: for instance, the Wilkoxon statistic, that adapts to non-normal data, did not perform better than other statistics in our experiments (see Fig. 9).

Measuring the reliability of group studies

The reliability of an activation pattern measures how systematically a given voxel or region will be found when performing a group study in one or other group of subjects. Taking advantage of the great number of subjects, we have used a bootstrap procedure, and two measures for assessing the reliability of the group studies: one models the activated/non-activated state of locations (voxel) as a mixture of binomial distributions, and quantifies the difference between the null and the active mode, while the other defines how well clusters of supra-threshold activity match across groups.

The first criterion has already been proposed in the literature ; it has the advantage of yielding very stable results across splits (see Figs. 6,9); one reason is that all the R groups are used in each single computation of the parameters, while the cluster-based measure is based on pairwise comparisons. However, care must be taken because the estimation may be trapped in local minima (although we have never observed convergence problems in our experiments), or because the joint estimation of the different parameters may imply some non-trivial interaction between the parameters (e.g. the sensitivity λ might not be independent from κ). More importantly, results at the voxel level are not as important as the presence of a strong local maximum or a significant cluster, which deserve being reported.

This has incited use to develop a second measure (see Eq. ( 9)), which takes into account only extended clusters and compares the position of their centres. Note that the penalty function Φ stabilizes to Φ ≃ 1 as soon as the distance exceeds 12mm.This is because clusters distant from e.g. 20 mm cannot be considered as homologous, as well as clusters distant from e.g. 50 mm (this is true because we are reporting group results; when reporting individual results, greater variability might be allowed). Averaging across supra-threshold clusters yields an idea of how frequently close clusters will be obtained across groups of subjects. This pairwise measure is somewhat more variable than the voxel-based indexes, but it yields an independent confirmation of possible differences in reliability.

As we have noticed, the dependence on the contrast studied is high with respect to the bootstrap dispersion, confirming the worthiness of these measures (Fig. 6). In general terms, κ and Φ have a similar behaviour (which corresponds to opposite numerical fluctuations: κ is high when Φ is low and vice versa). This was not obvious, given that the two measures are independent and based on completely different approaches. It suggests that our results are not artefactual, but really intrinsic to the data.

Our setting for the study of the reliability may also be used to compare competing pre-processing techniques or analysis frameworks, in addition to previous contributions based on cross-validation [START_REF] Strother | The quantitative evaluation of functional neuroimaging experiments: the NPAIRS data analysis framework[END_REF] and information theory [START_REF] Kjems | The quantitative evaluation of functional neuroimaging experiments: mutual information learning curves[END_REF].

Is there any best threshold ?

The fundamental question of finding an optimal threshold to label areas as activated has rarely been tackled, since it requires the modelling at the voxel level of both the null and the converse hypothesis to control both the false positive and false negative rates. This is possible here, thanks to the high number of subjects. Interestingly, we find a relatively low value for the optimal threshold (θ ⋆ ∼ 2.7 when considering κ, θ ⋆ ∼ 3 when considering Φ ; note that these two measures are independent). The corresponding P-values (0.0035 -0.001) are not conservative, so that such thresholds do not allow a very strict control of the rate of false positives. Family-wise error control procedures such as Bonferroni, Random Field Theory [START_REF] Ashburner | Human Brain Function[END_REF], and, to a lesser extent, False Discovery Rate [START_REF] Genovese | Thresholding of statistical maps in functional neuroimaging using the false discovery rate[END_REF], typically imply the use of much higher thresholds. In this study, we have chosen a relatively lenient threshold p < 10 -3 uncorrected because specificity control was not our main point. However, a good compromise between the control of false positives and the reliability may be the use of cluster-level or parcel-level inference.

We obtained very similar results with functional contrasts, like the auditory contrast, that are stronger in terms of signal, the only difference being that the optimal threshold was slightly higher, between 3 and 3.5. However, it is not obvious that our results generalize to datasets with different structure and our point is certainly not justify lenient thresholding procedures. Nevertheless, the question of an optimized threshold should be addressed more systematically in neuroimaging studies.

How does the sample size affect the reliability of the results ?

Another fundamental question concerns the number of subjects that should be included in a study. The point here is not only the sensitivity [START_REF] Desmond | [END_REF]Glover, 2002, McNamee andLazar, 2004], but also the reliability [START_REF] Murphy | An empirical investigation into the number of subjects required for an event-related fMRI study[END_REF] of the results. Our results clearly indicate that S = 20 is it a minimum if one wants to have acceptable reliability. As far as we know, most studies currently do not have this number of subjects, and one might be concerned with the reliability of many findings from neuroimaging studies.

One might object that in our case, only one session was available for each subject, and that the quick event-related design might yield poor results in terms of detection. However, the results that we describe are related to very basic contrasts (auditory and motor activity) for which we could check that most of the subjects (motor contrast) or even every subject (auditory contrast) had significant functional activity in expected regions, which has to be compared with the subtle functional contrasts that are currently topics of investigation. Moreover, our finding confirms earlier simulations and studies [START_REF] Desmond | [END_REF]Glover, 2002, Murphy and[START_REF] Murphy | [END_REF]. For these reasons, we think that this paper should be considered as an incentive to include more subjects in neuroimaging studies.

Reliability of the different statistical tests

One of the most important practical questions is to describe or design the most efficient ways to perform group studies in neuroimaging. Our experiments are not exhaustive. In particular, some other techniques could have been introduced, as well as combination of techniques used in this work. But our first aim is to find general guidelines.

First of all, we think that non-parametric assessment of functional activity should be preferred to analytical tests, which rely on unvalidated and sometimes wrong hypotheses. This can be done using adapted toolboxes e.g. SnPM [START_REF] Hayasaka | Validating Cluster Size Inference: Random Field and Permutation Methods[END_REF] or distance [Mériaux et al., 2006a]. It is worthwhile to note that the implementation of the tests in C reduces computation time to very reasonable time (cluster-level P-values can e.g. be computed in less than one minute on a dataset of ten subjects). Non-parametric estimation of the significance benefits both to the sensitivity and to the reproducibility of the studies.

Second, Mixed-Effects models should systematically be preferred to mere random effects analyses: there is some information in the first level of the data that improves the estimation of the group effects/variance and statistic.

Third, cluster-and parcel-based inference should be preferred to voxel-based thresholding. Cluster-level inference is of frequent use, which benefits to the sensitivity and the reliability of group analyses. However, it is based on the assumption that activated regions are large, which is not necessarily true. Parcel-based inference may thus be a interesting alternative, since it further allows some spatial relaxation in the subject-to-subject correspondence. The price to pay is a larger variability of the results due to a more crispy decision function (activated vs no activated). We recommend the combination of one of these techniques together with MFX. By contrast, stronger smoothing did not increase significantly the reliability of the results.

Fourth, to our surprise, ΨFX was found to be the most reliable technique. Although the statistic function does not take into account the group variance -as argued earlier, this is probably the reason of its higher performance -its distribution under the null hypothesis is tabulated by random swaps of the effects signs, so that it is indeed a valid group inference technique. However, care should be taken when using it: the thresholds have to be computed voxel per voxel, i.e. are not spatially stationary, second, the statistic value itself has no obvious interpretation, by contrast with the RFX or MFX statistics.

Conclusion

This analysis is a also start point to develop new strategies for brain mapping. Several directions may be addressed in the future:

• First trying to relate inter-subject variability to behavioral differences and individual or psychological characteristics of the subjects. Once again, such investigation may be undertaken only on large databases of subjects, and we the data basis used in this experiment might and probably will be used in such a framework.

• Second, efforts will further be made to relate spatial functional variability to anatomical variability. While some cortex-based analysis reports have indicated a greater sensitivity than standard volume-based mappings [START_REF] Fischl | Highresolution intersubject averaging and a coordinate system for the cortical surface[END_REF], statistical evidence is still lacking, and it is not clear at all how much can be gained when taking into account macroanatomical features, e.g. sulco-gyral anatomy. Similarly, diffusion-based imaging may add useful information to improve cross-subject brain cartography [START_REF] Behrens | A consistent relationship between local white matter architecture and functional specialisation in medial frontal cortex[END_REF].

• Third, at a statistical level, we think that intermediate levels of descriptions could be used more systematically between the subjects and the group level. Looking for outlier subjects, possible subgroups and so on can be investigated [START_REF] Kherif | Group analysis in functional neuroimaging: selecting subjects using similarity measures[END_REF], Thirion et al., 2005, Thirion et al., 2006], though finding a meaningful distance and separation criteria is not straightforward. For instance, it would be interesting to know what proportion of subjects had a significant activity in a given region; such a simple question requires to solve issues in across-subjects correspondences and in statistical thresholding (how can one be sure that two foci of activity in two subjects are homologous ?).

Finally, we hope that these procedures and results will provide useful guidelines in the experimental procedures and analysis of group functional neuroimaging datasets.

Figure 1 :

 1 Figure1: Illustration of the low sensitivity and weak reliability of supra-threshold patterns in standard group studies. (a) For a functional contrast that shows regions involved in a computation task, we show activity maps thresholded at a p < 0.001 level after a random effect analysis on 6 disjoint groups of 13 subjects; the position of the view is z=37mm in the MNI normalized space.(b) In the same plane, here is the same map computed from all the subjects together. Note the low sensitivity and weak reliability of the maps in (a).

Figure 2 :

 2 Figure 2: Example of mixture of binomial distribution. The empirical histogram of G(v) is modelled by the model in Eq. (7), with R = 8. The Y axis is in log-coordinates for the sake of readability.

Figure 3 :

 3 Figure 3: Statistical model of the effects for the left click-right click contrast, on S = 81 subjects. (a) z-value associated withe the RFX test; (b) group variance estimate; (c) z-value of the D'Agostino-Pearson test for normality of the effects β ; (d) z-value of the D'Agostino-Pearson test applied to the normalized effects β σ .Note that all the z-maps are thresholded at z = 8 for numerical reasons. The color scale of the variance image has been chosen arbitrarily in order to have supra-threshold areas that are comparable with the other maps. The variance is expressed in squared percentage of the BOLD mean signal.

Figure 4 :

 4 Figure 4: Statistical model of the effects for the audio instructions-video instructions contrast, on S = 81 subjects. (a) z-value associated withe the RFX test; (b) group variance estimate; (c) z-value of the D'Agostino-Pearson test for normality of the effects β; (d) z-value of the D'Agostino-Pearson test applied to the normalized effects β σ .Note that all the z-maps are thresholded at z = 8 for numerical reasons. The color scale of the variance image has been chosen arbitrarily in order to have supra-threshold areas that are comparable with the other maps. The variance is expressed in squared percentage of the BOLD mean signal.

Figure 5 :

 5 Figure 5: Statistical model of the effects for the computation-reading contrast, on S = 81 subjects. (a) z-value associated withe the RFX test; (b) group variance estimate; (c) z-value of the D'Agostino-Pearson test for normality of the effects β; (d) z-value of the D'Agostino-Pearson test applied to the normalized effects β σ .Note that all the z-maps are thresholded at z = 8 for numerical reasons. The color scale of the variance image has been chosen arbitrarily in order to have supra-threshold areas that are comparable with the other maps. The variance is expressed in squared percentage of the BOLD mean signal.

  Figure6: Dependence of the reproducibility and of the sensitivity of the group RFX analysis on the functional contrast under consideration. These results are obtained by drawing 5 disjoint groups of S = 16 subjects in the population of 81 subjects, and applying the procedure described in section 2.5.1. The threshold is θ = 3.1 (a) The results are more reliable for the contrast that shows auditory regions than for a contrast that shows motor activity or a contrast that shows the regions involved in the computation task. (b) By contrast the size of the putatively activated areas is greater for the contrast that shows regions involved in computation, and smaller for the contrast that shows the regions involved in motor activity. (c-d) The cluster variability penalty Φ is presented for clusters of more than η = 10 (c) or η = 30 (d) voxels. The behaviour is as expected, with the smallest value for the auditory-specific contrast.

Figure 7 :

 7 Figure7: Dependence of the reproducibility, the sensitivity, and the distance between suprathreshold clusters of the group RFX analysis on the threshold chosen to binarize the statistic maps. These results are obtained by drawing 5 disjoint groups of S = 16 subjects in the population of 81 subjects, and applying the procedure described in section 2.5.1. This is performed on the images of the left click-right click contrast. (a) The reproducibility index κ shows a maximum for θ ∼ 2.7. (b) The sensitivity decreases when θ increases.(c,d) The average distance between supra-threshold clusters of more than 10(c) or 30(d) voxels across groups has a minimum around θ ∼ 3.

Figure 8 :

 8 Figure 8: Dependence of the reproducibility κ (a), the sensitivity λ (b) and the average distance between supra-threshold cluster centroids Φ (c) of the group RFX analysis on the group size. The reliability is assessed considering disjoint groups of size S = 10, 13, 16, 20, 27, 40 within the population of 81 subjects. This is performed on the images of the left click-right click contrast. (a) The reproducibility increases with S and reaches a plateau for S > 20. (b) The size of putatively activated areas steadily increases with S. (c-d) The average intra-cluster distance decreases with S; it reaches a plateau for S > 20 when η = 10 (c), whereas it further decreases when η = 30 (d).

Figure 9 :

 9 Figure 9: Dependence of the reliability κ (a), the sensitivity λ (b), inter-supra-threshold cluster distance penalty Φ (c-d) of the statistical analysis on the group statistic used. Φ is based clusters of size greater than η = 10 (c) or η = 30 (d). These quantities assessed considering R = 8 disjoint groups of size S = 10 within the population of 81 subjects, using the left click-right click contrast.

Given our definition of the group mask, it may occur that functional data is available in a subsample of the population of size n, with S

≤ n ≤ S. In such a case, S is replaced by n in the formulas. In the present work, we systematically apply such corrections, even though we do not mention them.

Spatial mis-registration may be artefactual (incorrect normalization) or not (intrinsically different functional anatomy.)


Mihai Telescu

Riccardo Trinchero

Igor Simone Stievano

Noël Tanguy

Worst-Case Optimization of a Digital Link for Wearable Electronics in a Stochastic Framework

Keywords: stochastic optimization, stochastic analysis, metamodeling, textile transmission line, wearable electronics

  L'archive ouverte pluridisciplinaire

I. INTRODUCTION

In engineering, the term "uncertainties" designates tolerances and imperfections but also random environment variables such as temperature or humidity. Accounting for such uncertainties is a very important aspect of any design flow. For decades, Monte Carlo analysis has been a common strategy for a posteriori validation of a design [START_REF] Zhang | Development of robust interconnect model based on design of experiments and multiobjective optimization[END_REF]. Monte Carlo is the simplest form of stochastic analysis: system parameters vary within specific tolerance margins and the designer evaluates the adequacy of the system under test according to a chosen metric. Monte Carlo is computationally costly and in recent years more efficient techniques have been developed (see [START_REF] Manfredi | Stochastic transmission line analysis via polynomial chaos methods: an overview[END_REF][START_REF] Ding | How Affine Arithmetic Helps Beat Uncertainties in Electrical Systems[END_REF] and the references therein).

However, the common strategy is still to carry out design optimizations for a nominal case and use stochastic analysis as a posteriori validation. In some cases, yield optimization may take place, again, a posteriori, using the nominal design as a starting point.

In this paper we propose an alternative and less conventional approach. Uncertainties are embedded in the design flow from the early stages and parameter optimization algorithms become uncertainty-aware.

The feasibility and strength of the approach is demonstrated on a flexible differential link including a passive equalizer described in Fig. 1.

II. SCENARIO DESCRIPTION

Wearable electronics, sometimes referred to as "smart textiles", have been attracting a lot of interest in the past two decades [START_REF] Cottet | Electrical characterization of textile transmission lines[END_REF][START_REF] Marculescu | Electronic textiles: A platform for pervasive computing[END_REF][START_REF] Shi | Smart Textile-Integrated Microelectronic Systems for Wearable Applications[END_REF][START_REF] Nayak | Electronic textiles for military personnel[END_REF] There is particular interest for the field in the defense industry [START_REF] Nayak | Electronic textiles for military personnel[END_REF] but many civilian applications also exist. The idea of creating apparel featuring sensor hubs, processors and various interfaces brings new challenges in terms of signal and power integrity. Transmission lines are characterized by much higher fabrication imperfections than in the case of regular PCBs and suffer greater impact from environment parameters such as temperature and humidity. Bending may also be an issue.

In the framework described above, the example in Fig. 1 is considered. It consists of a differential serial link implemented by conductive copper wires weaved in a textile structure. The link features a differential driver at the near end and its paired receiver at the far end. A passive equalizer is inserted between the driver and the line. It is assumed that the equalizer and drivers are implemented either using conventional electronics or semi-flexible substrates and are affected by uncertainties that are negligible with respect to those affecting the textile transmission line. This assumption is reasonable if one assumes the link to be part of a military body armor, or specialized garment, featuring mixed materials. The objective is to optimize the equalizer by properly choosing the values of the passive components , , accounting for various uncertainties of the channel by observing the eye diagram on the receiver side.

III. STOCHASTIC OPTIMIZATION STRATEGY

A. General Mathematical Framework

Consider a cost function , . If is a constant and is a deterministic variable, finding the minimum of , is a classic problem. Local minima may in some cases pose a challenge, but it is an issue which has been addressed successfully in many ways. Methods such as genetic algorithms, particle swarm, ant colony and many others are now mainstream. Now assume that is a "hidden" random variable representing an uncertain parameter the designer cannot control directly; remains deterministic at this stage. Along the axis, multiple realizations of , exist as varies according to a probabilistic law. When conducting optimization in this framework, the very objective of the operation needs to be redefined.

Two scenarios appear obvious to the authors, although arguably many other interesting ones can be imagined. The first scenario is an average optimization which can be expressed as finding such that: Argmin ∈ , (i.e., the expected values computed for all the possible realizations of the cost function due to the random variable is minimal). The second is worst-case optimization which can be defined as follows.

Let

!"

be a new deterministic cost function such that:

!" #$ % ∶ ' , < ≥ α+, ∈ (1)
where ' denotes the probability computed on the stochastic variable for any values of the deterministic variable ∈ and , ∈ 0,1 is a parameter defining the level of confidence.

Naively put, if one runs a very large number of Monte Carlo simulations for every value of and collects the maxima, then these maxima form

!"
for , → 1 . The worst-case optimum is such that:

Argmin ∈ 0 !" 1
Note that the optimization process does not require analytical knowledge of

!"
but simply a reduced number of pointwise estimations.

Once the criterion has been defined, let us assume that the designer attempts to solve the problem with an iterative algorithm of their choice. As the algorithm progresses, for every value of , the probability distribution of needs to be explored. This can quickly become computationally challenging.

Moreover, in many practical cases one will have to deal with multivariable cost functions, 2 , . . , 4 , 2 , . . , 5 , yielding a problem that is practically unsolvable by mere brute force and requires the use of computationally efficient methods. A solution is provided in paragraph D. Also, for practical reasons, confidence level , will be equated to 0.99 in this implementation.

Throughout the paper random variables 2 , 7 , . . , 4 will be referred to as subject variables. The designer has no direct control over them. Variables 2 , 7 , . . , 5 will be referred to as object variables. We consider object variables deterministic at this stage although paths to relax this assumption are being investigated.

B. Particular case of textile interconnect line

For the test case considered in this paper we have two object variables characterizing the equalizer ( 2 , 7

) and three subject variables characterizing the link, 2,7,8 . Coefficient 2 translates the variation of the characteristic impedance according to the tolerances available in [START_REF] Cottet | Electrical characterization of textile transmission lines[END_REF]. It is caused by imperfections in the fabrication process. Coefficient 7 represents temperature and 8 the variation of the line length.

A metric derived from the eye diagram, known as eye SNR and subsequently denoted as 94: is used to evaluate the quality of the link, thus playing the role of the cost function 2 , 7 , 8 , , in the previous section. It is defined as:

94: ; < =; > ? < @? > (2)
where A 2 and A B are the logical "1" and respectively "0" NRZ levels and C 2 and C B are the variances of the distributions of all samples interpreted as either "1" or "0". Basically 94: is a measure of the eye vertical opening weighted by the "smudginess" of the eye. In our implementation the cost function will be -94: since the common convention in optimization is to seek minima rather than maxima.

In this conference paper the authors focus on worst-case optimization essentially seeking an equalizer defined by { F, F} which guarantees a level of 94: in the most unfavorable conditions defined by { 2 , 7 , 8 }.

C. Algorithm

Modern optimization algorithms typically use multiple evaluations of the cost function per iteration in an attempt to converge on the global minimum and avoid being trapped in local minima. The authors have successfully experimented with both particle swarm and genetic algorithms for a variety of examples in the stochastic framework. The results presented in this paper are based on particle swarm optimization (PSO).

PSO is a heuristic optimization method inspired by the behavior of birds [START_REF] Kennedy | Particle swarm optimization[END_REF]. Particles are "virtual agents" meant to explore a manifold space and described by their position and their velocity. Let point H defined by 2 2I , . . , 5

5I

be the position of a particle. As the algorithm iterates H follows the law:

H J + 1 H J + L J + 1 (3)
where L denotes the "velocity" of the particle (a displacement in the manifold space in fact) and J denotes the iteration. [START_REF] Kennedy | Particle swarm optimization[END_REF][START_REF] Chatterjee | Nonlinear inertia weight variation for dynamic adaptation in particle swarm optimization[END_REF]. Once the cost function is evaluated in point H , the velocity is updated according to an equation that takes into account the best position occupied by the particle, the overall best position ever occupied by any particle in the swarm, but also includes an "inertia" -a tendency to maintain a trajectory. The particle thus "learns" from its own experience but also from the collective experience of the swarm. In our implementation velocity and inertia are updated according to [START_REF] Chatterjee | Nonlinear inertia weight variation for dynamic adaptation in particle swarm optimization[END_REF].

The following steps provide a general description of the optimization strategy, without specific implementation details which are beyond the scope of this conference paper.

1. Select the number of particles in the swarm; initialize positions, velocities, inertias. For our test case the term "positions" designates various pairs { F, F} (following a Latin hypercube method for example).

2. For every pair { F, F} compute a closed form of the probability density function of 2 , 7 , 8 , , using the method described in paragraph D.

Use probability density functions to evaluate metric (in our case

!"

).

4. Select best case among the explored pairs { F, F}.

5. Update inertia, velocities, then positions.

6. If maximum number of iterations or convergence tolerance is reached end. Else go to step 2 and repeat.

Note that what makes the optimization stochastic is the presence of step 2 and the nature of the metric in step 3. If step 2 is suppressed and a deterministic metric is used in step 3 the algorithm reverts to traditional, deterministic PSO.

D. Kernel Density Estimation

Let M I denote the probability density function of 2 , . . , 4 , 2 , . . , 5 in a given point H defined by 2 2I , . . , 5 5I . The shape of this probability density function depends on random variables 2 , . . , 4 . It is possible to derive a non-parametric, unsupervised, model of M I via what is known as kernel density estimation [START_REF] Bowman | Applied Smoothing Techniques for Data Analysis[END_REF] according to the expression:

M I ≅ 2 !O ∑ Q 0 " R =" RS O 1 ! T2
where I is a realization of I 2 , . . , 4 , 2 , . . , 5 for 2 2 , . . , 4 4 with 2 2I , . . , 5 5I . Q is the socalled kernel; the Gaussian function is a common choice of kernel. Finally, ℎ > 0 is a smoothing parameter called window or bandwidth.

The main advantage of this approach is that it may be used even when W, the number of subject variables, is very large. Also, it is readily available in Matlab as the ksdensity function.

For the results obtained in the present paper the authors used this standard implementation, albeit with a narrow windowing more adequate for capturing the "tail" of the probability density function.

IV. RESULTS

The context is the one of the link of Fig. 1. An NRZ encoding is assumed at a data rate of 2 Gbs. For the sake of simplicity linear models (similar to the ones implemented in IBIS-AMI simulations were used to model the buffers). If needed enhanced models of such devices can be adopted [START_REF] Signorini | Present and future of I/O-buffer behavioral macromodels[END_REF][START_REF] Siviero | Constructive Signal Approximations for Fast Transient Simulation of Coupled Channels[END_REF][START_REF] Stievano | Reliable eyediagram analysis of data links via device macromodels[END_REF]. All computations of 94: used a pre-recorded random stream of 3000 bits and were performed using Matlab function comm.EyeDiagram. Future implementations might take advantage of recent advances in fast eye diagram scoring (see [START_REF] Li | An eye diagram improvement method using simulation annealing algorithm[END_REF] and the references therein).

The algorithm presented in section III.C is used to compute an equalizer optimized for a worst-case scenario. The results highlight the difference between deterministic (nominal) optimization and stochastic optimization. A large simulation involving a systematic sweep of parameters and was also conducted in order to validate the optimization results. It consisted in a grid search investigation of 100 { F, F} points and 100 Monte Carlo runs per point. A worstcase optimum was in { F, F} {90 , 2.5 `a} (see Fig. 3), which turns out to be different from the nominal case optimum illustrated in Fig. 2, { F, F} {77 , 3.1 `a}.

Figures 2 and 3 illustrate the third iteration of a PSO process with 15 particles in a deterministic case (figure 2) and, respectively, in a stochastic case (figure 3). Yellow contour lines allow the identification of the zones where the value of 94: is large. The optima obtained by systematic inspection are indicated by the crossing of dashed lines. Crosses mark the position of the particles. The global best obtained by the swarm is indicated by a black dot. Note that the simulations are distinct (i.e. the results of the nominal optimization are not used to initialize the stochastic one).

The figures illustrate that in one case the algorithm is simply faced with a deterministic problem (a single contour map) whereas in the second case multiple realizations of the contour map are possible (only a selection was actually plotted).

As far as the efficiency is concerned, it is important to point out that the stochastic optimization process reaches the optimum after only 3 iterations with 15 particles, each iteration requiring 25 Monte Carlo runs per particle.

Fig. 1 .

 1 Fig. 1. System under test: digital link for wearable applications.

  The nominal per unit line parameters of the textile interconnect are A ! X 0.479 [/] , ! X 25.5 `a/] and ! X 5.63 /]. The nominal line length is d ! X 30 e]. Uncertainties are taken into account according to the following relations 912, 1.074 , 7 ∈ -5, 35 and 8 ∈ 0.85, 1.15 . One may note that variable 2 impacts line parameters in correlated manner. It translates the variation in impedance due to imperfect manufacturing and matches tolerances indicated in[START_REF] Cottet | Electrical characterization of textile transmission lines[END_REF]. Variable 7 is simply the ambient temperature and affects the resistance of the copper wire. Variable 8 translates the variation in line length due to garment sizing as well as fabrication (cutting and weaving).

Fig. 2 .

 2 Fig. 2. Traditional optimization (nominal case). Third iteration of PSO.

Fig. 3 .

 3 Fig. 3. Stochastic optimization. Third iteration of PSO. Multiple realizations of the cost function are displayed.

Without equalization, the value of 94: is 2.3 in a worstcase scenario. Introducing the stochastically optimized shunt equalizer increases 94: to 4.1. The respective eyes are compared in Fig. 4.

V. CONCLUDING REMARKS

This paper investigates the feasibility of an alternative worst-case optimization scheme for textile interconnects. Specifically, an equalizer was designed by an optimization process which considers system uncertainties from the very beginning. The proposed approach combines an unsupervised methodology, which is used to deal with the stochastic nature of the cost function in an efficient way, and a custom implementation of a particle swarm optimizer. The results clearly demonstrate the effectiveness of the proposed strategy. Uncertainty-aware optimization is an interesting approach which can be a very powerful tool for designers who need to meet ambitious performance boundaries while having to work with high tolerances or with models where deterministic assumptions no longer hold.