Benjamin Mahieu
email: benjamin.mahieu@inrae.fr

Pascal Schlich
email: pascal.schlich@inrae.fr

Michel Visalli

Hervé Cardot

A multiple-response chi-square framework for the analysis of Free-Comment and Check-All-That-Apply data Authors

Keywords:

Highlights

-The usual chi-square framework is not well suited to analyze FC and CATA data -A modified multiple-response chi-square framework is introduced -This new framework takes into account the specificities of FC and CATA data -R package (MultiResponseR) available upon request to the authors and on GitHub

Introduction

Free-Comment (FC) [START_REF] Ten Kleij | Text analysis of open-ended survey responses: a complementary method to preference mapping[END_REF] and Check-All-That-Apply (CATA) [START_REF] Adams | Advantages and uses of check-all-that-apply response compared to traditional scaling of attributes for salty snacks[END_REF] are word citation occurrence-based methods that aim at collecting product descriptions from consumers using either their own words or a mutual predefined list of descriptors. These descriptions are collected without any quantification or product comparison. At the panel level, the collected data constitute count data that are usually stored in a contingency table that contains the number of times each descriptor (in columns) was cited for each product (in rows).

The analysis of these data starts by testing whether overall differences exist between the products. Two approaches can be distinguished to do so. The first one consists of performing a chi-square test while the second one is based on a combination of Cochran's Q statistics [START_REF] Meyners | Existing and new approaches for the analysis of CATA data[END_REF]. Pursuing the analyses further is only recommended if the existence of overall differences between products is established. In this case, these differences can be visualized using Correspondence Analysis (CA). CA enables to represent the structure of the dependence between products and descriptors on a factorial map that decomposes the whole dependence into axes of maximal and decreasing dependence. As a final step of the analysis, it is important to determine which descriptors are significantly associated with which product. Again, two approaches can be distinguished to do so.

The first one is multidimensional alignment [START_REF] Meyners | Existing and new approaches for the analysis of CATA data[END_REF] that consists of considering a descriptor significantly positively (resp. negatively) associated to a product when their vectors in the sensory space depicted by the CA form an angle lower than or equal to 45° (resp. higher than or equal to 135°). The second approach consists of testing each cell of the contingency table against the null hypothesis of independence using a chi-square test or a Fisher's exact test (Mahieu, Visalli, & Schlich, 2020a;[START_REF] Symoneaux | Comment analysis of consumer's likes and dislikes as an alternative tool to preference mapping. A case study on apples[END_REF].

All of these approaches but the combination of Cochran's Q statistics are based on the chi-square statistic. The chi-square statistic can be directly used to test for overall differences between the products before performing the CA. The total inertia of CA is the chi-square statistic divided by the grand sum of the contingency table, also called phi-square index. Since multidimensional alignment relies on the CA, it depends also on the chi-square statistic. Finally, the tests per cell approach directly rely on the chisquare statistic since Fisher's exact test can, roughly speaking, be seen as an exact chi-square test.

These common practices assume that all citations are independent experimental units within an evaluation, which is not the case since citations of descriptors by a given subject for a given product are not independent. Instead, one evaluation, i.e. the entire set of descriptors cited by one subject for one product, should be considered as an experimental unit [START_REF] Loughin | Testing for Association in Contingency Tables with Multiple Column Responses[END_REF]. Indeed, considering citations as experimental units implies computing incorrect expected values under the null hypothesis of independence between products and descriptors [START_REF] Loughin | Testing for Association in Contingency Tables with Multiple Column Responses[END_REF], resulting in an incorrect chi-square statistic. Subsequent analyses of FC and CATA data based on this chi-square statistic are thus also incorrect and can sometimes lead to wrong interpretations.

The present paper aims to overcome the previous limitations by introducing the multiple-response chi-square framework based on the multiple-response chi-square statistic of [START_REF] Loughin | Testing for Association in Contingency Tables with Multiple Column Responses[END_REF]. This new framework considers experimental units as being the evaluations rather than the citations. First, some notations are introduced and the multiple-response chi-square test of [START_REF] Loughin | Testing for Association in Contingency Tables with Multiple Column Responses[END_REF] is presented and adapted to the context of FC and CATA data. Second, the multipleresponse Correspondence Analysis (MR-CA) is introduced. Third, the transposition of the methodologies presented in Mahieu et al. (2020a) to the multiple-response chisquare framework is established. Fourth, examples of the benefits of the new framework are given on real CATA data. Finally, an overall discussion and a conclusion are given. Let us denote by the probability of descriptor to be cited for product . What is under investigation is whether differs from one product to another. Using the above notations, the following hypotheses are considered:

Material and methods

: = ⋯ = = , ∀ ∈ 1, … , : !"#$ $ ∈ 1, … , % , & ∈ 1, … , '# ℎ ≠ & $*+ℎ %$ ≠ ,
Note that this does not correspond to a classical test of homogeneity since, for each product , multiple descriptors can be selected. Under the null hypothesis, the expected number of citations of descriptors for product , denoted by -. /, is equal to × and can be estimated by × 0 . The following test statistic, called multiple-response chi-square statistic, is thus introduced:

1 23 4 = 5 5 6 - × 0 8 4 × 0 9
As × 0 = × : 0 × 0 ;, 1 23 4 can also be expressed as:

1 23 4 = 5 5 < -× : 0 × 0 ;= 4 × : 0 × 0 ; 9
As in [START_REF] Loughin | Testing for Association in Contingency Tables with Multiple Column Responses[END_REF], it can be shown that the asymptotic distribution of this test statistic under the null hypothesis is complicated because descriptors might not be selected independently. A reasonable option for estimating the distribution of In usual CA, the products are compared to each other according to their profile. The profile of each product is defined as the proportion of citations of each descriptor for this product relatively to the total number of citations (all descriptors combined) elicited by this same product. Thus, in the context of FC and CATA data, when products elicit different average citation rates (all descriptors combined) then absolute differences in descriptors' citation rates between products are distorted due to this "citation rescaling". The degree of distortion depends on the degree of differences in citation rates between products. For more details on the usual CA, one can refer e.g. to [START_REF] Greenacre | Correspondence Analysis in Practice[END_REF]. The previous assertions are also applicable to Hellinger-distance-based CA [START_REF] Rao | A Review of Canonical Coordinates and an Alternative to Correspondence Analysis Using Hellinger Distance[END_REF][START_REF] Vidal | Comparison of Correspondence Analysis based on Hellinger and chi-square distances to obtain sensory spaces from check-all-that-apply (CATA) questions[END_REF] because this latter is also based on the products' profiles.

MR-CA overcomes the above limitation by scaling products according to their number of evaluations instead of their number of received citations. It results in comparing products based on their average proportions of citations for each descriptor. This "evaluation scaling" only has importance in the case of unbalanced design. Indeed, products that are more evaluated are likely to elicit more citations of all descriptors and it is necessary to put products on an equal footing before comparing them. To summarize, the propensity of some products to elicit more citations than others does not affect MR-CA while it does with usual CA.

When applied to FC and CATA data, MR-CA can be seen as standing at the frontier between the usual CA of the descriptor by product contingency table and the PCA of the products' average profiles depicted by the descriptors' proportions of citations.

MR-CA performs the PCA of the products' average proportions of citations but weighting the descriptors proportionally to their citation rate as in usual CA.

Definition

Similarly, to the usual CA based on the singular value decomposition of the matrix of standardized residuals defined by the usual chi-square statistic, the MR-CA is based on the singular value decomposition of the matrix of standardized residuals defined by the multiple-response chi-square statistic. Using the notations defined in the previous section, let us consider:

-> a column matrix of size × 1 whose elements equals 0 , ∈ 1, … , -? a column matrix of size × 1 whose elements equals 0 , ∈ 1, … , -@ 3 a diagonal matrix of size × whose diagonal elements equal 0 , ∈ 1, … , -@ A a diagonal matrix of size × whose diagonal elements equal 0 , ∈ 1, … , -B a matrix of size × whose general term equal 0 , ∈ 1, … , , ∈ 1, … , Using these notations, the MR-CA is based on the singular value decomposition of the matrix C defined as:

C = @ 3 D 4E B ->? F G@ A D 4
Let us denote by H the matrix of left singular vectors of C, I the diagonal matrix of singular values of C and J the matrix of right singular vectors of C such that C = HIJ F . Similarly to the usual CA, the principal coordinates of the products are defined as @ 3 D K L HI and the so-called contribution coordinates [START_REF] Greenacre | Contribution Biplots[END_REF] of the descriptors are defined as J. Note that since this system of coordinates defines a strict biplot as defined in [START_REF] Gabriel | The biplot graphic display of matrices with application to principal component analysis[END_REF], it is suggested to use arrows rather than points to display the descriptors' coordinates. This could help practitioners to remember to interpret relations between products and descriptors as scalar products (orthogonal projection) and not "proximities". Different systems of coordinates could be used for displaying results of MR-CA similarly to usual CA [START_REF] Greenacre | Tying up the loose ends in simple, multiple, joint correspondence analysis[END_REF].

However, the one proposed here has two benefits: it enables interpreting maps similarly to Principal Component Analysis (PCA) biplots and the coordinates of the columns (descriptors) reflect their respective contribution to the inertia and to the distances between rows (products) [START_REF] Greenacre | Tying up the loose ends in simple, multiple, joint correspondence analysis[END_REF].

Equivalently, the MR-CA can be defined as the PCA of the matrix @ 3 D B@ A D K L . This latter definition of MR-CA better highlights that the distance between two products ≠ ′ ∈ 1, … , in the sensory space depicted by MR-CA called multiple-response chi-square distance is equal to:

N OP L E , & G = Q 5 < - , , = 4 9
From the definition of the multiple-response chi-square distance, one can see that the weight given to each product is proportional to its number of evaluations rather than its number of received citations as it is in usual CA. Finally, it should be noted that the number of axes obtained by MR-CA is equal to the minimum between -1 and , as in a PCA in which descriptors act as variables and products as individuals, while in usual CA it is equal to the minimum between -1 and -1. This difference in the number of axes is because usual CA centers both rows (products) and columns (descriptors) while MR-CA centers only rows.

Statistical inference for multiple-response chi-square framework

This section transposes the methodologies from Mahieu et al. (2020a) to the multiple chi-square framework.

The dimensionality test of the dependence

Conceptual aims for Free-Comment and Check-All-That-Apply data

The aim of this test is twofold. First, it investigates if at least one axis of the MR-CA is significant, that is if some overall differences exist between the products. If no axis is significant, interpreting subsequent analyses including the outputs from MR-CA might lead to over-interpretations. If at least one axis is significant, the second aim of the test is to determine the number of axes that can be considered significant and thus interpreted. Because drawing sensory conclusions based on more than three or four axes can be difficult visually, the number of significant axes is taken into account in subsequent proposed analyses, which are simpler to interpret from a sensory point of view.

Technical aspects

It is possible to test if the dependence of each MR-CA axis is significant with a stepwise procedure similarly as for the usual CA (Mahieu et al., 2020a). The idea is to test, at each step R ER > 1G, whether the hypothesis of independence between products and descriptors is still rejected while the dependence captured by the axes 1 to R -1 was removed. In other words, it is tested if the strength of the dependence is still large enough to be considered significant.

As seen in the previous section, the total number of MR-CA axes, denoted T, is W U = :@ 3 4 C U @ A 4 + >? F ; ×

The multiple-response chi-square test associated with the test statistic 1 23 4 U enables testing if the k-th axis of the MR-CA captures a significant dependence between products and descriptors. Note that if R = 1 then this test corresponds to the multipleresponse chi-square test defined in section 2.1.

The multiple-response chi-square statistic of the products by descriptors contingency table is related to the eigenvalues of the MR-CA by the following equation:

1 23 4 = × 5 X Y Z Y
where 1 23 4 is the multiple-response chi-square statistic of the contingency table, is the total number of evaluations and X Y is the i-th eigenvalue of the MR-CA. This relation enables to compute each 1 23 4 U as:

∀R, 1 23 4 U = × 5 X Y Z Y U
To estimate the distribution of each 1 23 4 U under the null hypothesis, it is proposed to randomly permute the response vectors along products within each subject (Mahieu et al., 2020a;[START_REF] Meyners | Existing and new approaches for the analysis of CATA data[END_REF][START_REF] Meyners | Statistical inference for temporal dominance of sensations data using randomization tests[END_REF][START_REF] Wakeling | A New Significance Test for Consensus in Generalized Procrustes Analysis[END_REF][START_REF] Winkler | Multi-level block permutation[END_REF], a response vector referring to all citations given for one product by one subject.

To summarize, the dependence between products and descriptors captured by each MR-CA axis can be tested following these steps:

Confidence ellipses and discrimination of the products

In MR-CA, as well as in every multivariate analysis providing a product map, superimposing confidence ellipses on product coordinates is crucial to estimate if products are well discriminated. A total bootstrap procedure [START_REF] Cadoret | Construction and evaluation of confidence ellipses applied at sensory data[END_REF] is proposed to achieve this objective. This procedure consists of generating virtual panels by randomly resampling with replacement the subjects of the actual panel. Then, the product configurations of the virtual panels are rotated on the product configuration of the actual panel thanks to Procrustes rotations. A confidence ellipse is then constructed for each product based on the coordinates of its rotated bootstrap replicates. It is proposed to rely on the significant axes, indicated by the test of dependence presented in section 2.3.1, to determine the number of axes to account for the Procrustes rotations in the total bootstrap procedure.

For each pair of products, to determine if the two products are significantly different, it is proposed to rely on the total bootstrap test [START_REF] Mahieu | Free-comment outperformed check-all-that-apply in the sensory characterisation of wines with consumers at home[END_REF] considering the null hypothesis that the two products are not different. For each pair of products, a canonical discriminant analysis based on the rotated bootstrap replicates of the two products is performed. The rotated bootstrap replicates of the two products are then projected on the axis resulting from the canonical discriminant analysis. The distribution of the paired differences of the projected bootstrap replicates is estimated. Finally, the probability of zero to belong to this distribution is estimated and used as a p-value of the test. It is proposed to perform the total bootstrap tests on the significant axes.

2.3.3. Determination of the significant associations between products and descriptors: multiple-response tests per cell

Conceptual aims for Free-Comment and Check-All-That-Apply data

These tests aim to investigate the relations between descriptors and products. In particular, they investigate for a given descriptor and a given product if this descriptor is cited for this product in a proportion that significantly differs from the overall average citation proportion of this descriptor all products combined. The tests can be one-sided (positive differences) or two-sided (both positive and negative differences):

this choice is up to the discretion of the practitioner. A discussion is given about this choice in Mahieu et al. (2020a).

Technical aspects

It is proposed to define a multiple-response test per cell to test the following hypotheses for a given ∈ 1, … , and a given ∈ 1, … , :

: = : ≠
The multiple-response test per cell is based on a Monte-Carlo procedure. In this procedure, for each product ∈ 1, … , , evaluations are randomly drawn among the subjects having evaluated and only one evaluation is randomly drawn among each of these subjects. This enables constructing a virtual contingency table under the null hypothesis accounting for both the subject structure of the data and the nonindependence of the citations. Indeed, one evaluation is randomly drawn from each subject having evaluated and one randomly drawn evaluation (that respect the joint distributions of citations of the descriptors) contributes to several cells in the virtual contingency table .    A large number of virtual contingency tables under the null hypothesis can be generated by repeating this procedure. Then, for each cell, the proportion of

E * G
under the null hypothesis having an equal or a more extreme value than the observed constitute a p-value of the test. The multiple-response tests per cell can be performed with a two-sided alternative hypothesis or a one-sided greater alternative hypothesis.

Finally, it is proposed to perform the multiple-response tests per cell on the derived contingency table corresponding to the significant axes (Mahieu et al., 2020a), denoted W aYb , and defined following the reconstitution formula as:

W aYb = :@ 3 4 C aYb @ A 4 + >? F ; ×

Where C aYb = H aYb I aYb J aYb F with H aYb the matrix of left singular vectors of C corresponding to the significant axes, I aYb the diagonal matrix of singular values of C corresponding to the significant axes and J aYb the matrix of right singular vectors of C corresponding to the significant axes.

To perform the multiple-response tests per cell on W aYb rather than on the observed contingency table results in a gain of power without any inflation of the type I error as suggested by the simulation results presented in the Appendix. The simulation results also suggest that the smaller the number of significant axes and the intensity of the dependence between products and descriptors, the higher the gain of power.

Examples

These examples from two CATA datasets aim to compare outputs obtained from analyses belonging to the usual chi-square framework to those obtained from analyses belonging to the multiple-response chi-square framework. Although these examples deal with CATA datasets, note that the multiple-response chi-square framework is also appropriate to analyze FC data.

Datasets

The datasets are the same from [START_REF] Mahieu | An investigation of the stability of Free-Comment and Check-All-That-Apply in two consumer studies on red wines and milk chocolates[END_REF].

The study took place at the Barry Callebaut© Company, Belgium. Seventy regular consumers of milk chocolates (at least once every two weeks) were recruited among the employees of the Barry Callebaut© Company (not implied in sensory and consumer research). They performed a CATA task on four milk chocolates having different recipes: a standard Belgian milk chocolate, a Swiss milk chocolate, a milk compound chocolate, and a protein base milk chocolate. The four products were presented according to a Williams Latin square design. For each product, the CATA task was carried out according to two sensory modalities: texture in the mouth followed by flavor in the mouth. All the CATA descriptors were selected thanks to the expertise of sensory experts from the Barry Callebaut© Company. The collected data were then stored in two contingency tables, one per sensory modality, by cross tabulating the citation counts of the descriptors (columns) by the products (rows).

Since sensory interpretation is out of the scope of this paper dedicated to the comparison of the two chi-square frameworks, the descriptors were renamed D1, D2, etc. and the products were renamed P1, P2, P3, and P4. Finally, for the texture dataset, an additional product called P5 was artificially created. This product is exactly P4 except that for P5 the number of received citations for every descriptor has been divided by two as compared to P4. This was made to illustrate the differences between the multiple-response chi-square framework and the usual chisquare framework.

Analyses

All analyses were performed using R 4.0.2 (R Core Team, 2020). The analyses belonging to the multiple-response chi-square framework were performed using the R package "MultiResponseR" developed for this purpose by the authors.

The two contingency tables were analyzed using the following procedure. An alpha risk (Type I error) of 10% was considered as the significance level.

The dimensionality of the dependence between products and descriptors was determined within each chi-square framework using the dimensionality test (2000 simulations) presented in Mahieu et al. (2020a) for the usual chi-square framework and using the dimensionality test (2000 simulations) presented in section 2.3.1 for the multiple-response chi-square framework.

When at least one axis was significant, the corresponding CA (usual or multipleresponse) was performed on the contingency table. Outputs of each CA were displayed using a standard biplot [START_REF] Greenacre | Contribution Biplots[END_REF]. For each CA, confidence ellipses for the products' coordinates in the sensory space were computed with a total bootstrap procedure using 2000 bootstrap samples. The Procrustes rotations were performed on the significant axes. For each pair of products, a total bootstrap test was performed on the significant axes for assessing the significance of product difference.

For each pair of product and descriptor (cell), a Fisher's exact test was performed for the usual chi-square framework and a multiple-response test per cell as described in section 2.3.3 (2000 simulations) was performed for the multiple-response chi-square framework. All tests per cell were performed with a one-sided greater alternative hypothesis and conducted on the derived contingency table corresponding to the significant axes. 1 shows that whatever the sensory modality and the axis considered, the eigenvalues of the CA are higher in the multiple-response framework than in the usual one. This suggests that the usual framework underestimates the dependence between products and descriptors. This line of reasoning is reinforced by the example treated by [START_REF] Loughin | Testing for Association in Contingency Tables with Multiple Column Responses[END_REF] as they obtained a lower p-value (which is partly a function of the effect size) for their chi-square test in the multipleresponse framework than in the usual one. On the dimensionality of the dependence, Table 1 shows that similar conclusions are provided between products and descriptors by the two chi-square frameworks concerning the flavor dataset: three axes capture significant dependence. However, the dependence on the third axis appears more certain (p=0.0054) in the multiple-response chi-square framework than in the usual one (p=0.0914). Concerning the texture dataset, only two axes capture significant dependence within the usual chi-square framework while three axes capture significant dependence within the multiple-response chi-square framework.

Results

Sensory modality

Chi

Fig. 1 shows that for the texture dataset, the maps depicted by the two first axes of the usual CA (Fig. 1(a)) and the MR-CA (Fig. 1(b)) are very similar: all the products except P5 and all the descriptors have the same position on the two maps. The only difference between these maps is the location of P5 being different from P4 and closer to the origin in MR-CA (Fig. 1(b)) as compared to usual CA (Fig. 1(a)). The reason for this difference lies in the fact that P4 and P5 have the same profile (repartition of citations) in the usual CA. On the contrary, the MR-CA captures that P5 received fewer citations than P4 for all the descriptors but still following the same pattern of association with the descriptors. This explains the position of P5 relative to P4: P5 deviates from independence in the same direction that P4 (same pattern of association with the descriptors) but P5 is closer to the origin of the coordinates system than P4 (received fewer citations). Concerning the third significant axis obtained with the multiple-response chi-square framework on the texture dataset (Fig. 1(c)), it mainly traduces that P5 received fewer citations than P4 for all descriptors, which is logical. Note that the usual CA is unable to capture this difference between P4 and P5, which explains the non-significance of the third axis for this CA.

For the flavor dataset, Fig. 2 shows that the spaces provided by the usual CA and the MR-CA exhibit different configurations for both products and descriptors. For every descriptor, there is at least one other product that received more citations than P3.

Thus, in MR-CA, it is associated with no descriptor, which explains its position: P3 lies at the opposite of every descriptor loadings (Fig. 2(c) & Fig. 2(d)). On the contrary, in usual CA, P3 seems to be associated with D1, D5, and D6 and slightly with D2 (Fig. 2(a) & Fig. 2(b)). Indeed, in usual CA, the number of citations received by P3 for every descriptor is rescaled according to its total number of received citations. Thus, the fact that for every descriptor there is at least one other product that received more citations than P3 is erased in the usual CA. These features of P3 are the principal explanation of the differences between the spaces provided by MR-CA and usual CA, both applied on the flavor dataset. Another notable difference is between the maps depicted by the two first axes of the usual CA and of the MR-CA (Fig. 2(a) & Fig. 2(c)). On these maps, P1 and P2 appear to be more associated with D1, D4, and D6 in MR-CA as compared to the usual CA. This difference is due to the opposite phenomenon that occurred with P3: P1 and P2 received much more citations than P3 and P4 for these descriptors and the usual CA distorts this difference while the MR-CA does not.

Concerning the total bootstrap tests, whatever the considered sensory modality and whatever the considered chi-square framework, the conclusions they provided were the same except when considering the pair P4/P5 and the texture dataset. In the usual chi-square framework, P4 is for sure not different from P5 (p = 1) while P4 and P5 are significantly different in the multiple-response chi-square framework (p < 0.001). Of course, this is perfectly in line with Fig. 1 For texture, Fig. 3 shows that differences in the significant associations concern the pairs: P2-D2, P5-D3, and P5-D8. The pair P2-D2 is significant in the multipleresponse and not in the usual chi-square framework because P2 received more citations of this descriptor than the other products except P3. Concerning the product P5, it is noticeable that in the usual framework, it is significantly associated with the same descriptors as P4 (D3, D5, and D8), which was expected since P5 has the same profile that P4 in this framework. The pair P5-D8 being significant in the usual framework with a percentage of citations (25.71%) lower than the one of product P3 (28.57%) which is not significant nicely illustrates the issue of the "citation rescaling" due to considering the citations as experimental units. Since P3 and P5 were evaluated the same number of times, it is counterintuitive to have the one with the lowest proportion of citations significant and not the other. However, in the multipleresponse framework, both P3 and P5 are not significantly associated with D8, which is consistent. Regarding the pair P5-D3, the association is not significant in the multiple-response chi-square framework while it is in the usual chi-square framework. This difference is due to the "citation rescaling" that occurs in the usual chi-square framework and not in the multiple-response one.

Concerning the flavor dataset, several differences are shown in Fig. 3 between the conclusions provided by the two chi-square frameworks on descriptor by product significant associations. As was suggested by Fig. 2, P1 and P2 are significantly associated with D1, D4, and D6 in the multiple-response chi-square framework while only P2 is significantly associated with only D4 and D6 in the usual chi-square framework. This difference is because P1 and P2 received much more citations than P3 and P4 for these descriptors. On the contrary, without the "citation rescaling", since P3 and P4 received fewer citations, they got less significance in the multiple response framework; precisely, P3-D5 and P4-D3 are no longer significant in this framework. Finally, it is noticeable that the counterintuitive conclusion in the usual chi-square framework on the significant association of D3 with P4 and not with P1 and P2 while these received a higher percentage of D3 citations than P4, no longer holds in the multiple-response chi-square framework.

Discussion

To the best of our knowledge, it is the first time that a chi-square framework properly taking into account multiple-response data is introduced. The proposed analyses including the test of dimensionality, the product confidence ellipses, the pairwise product comparisons, and the product by descriptor association tests, the three of them being conducted on the significant axes, are all originals. This multipleresponse chi-square framework fits perfectly to FC and CATA data. However, this multiple-response chi-square framework is not restricted to be used only in sensory and consumer science and can be used to analyze any multiple-response data whatever the field they come from.

The examples presented in this paper showed that the multiple-response chi-square framework is better suited than the usual chi-square framework to analyze FC and CATA data. A major benefit of using the multiple-response chi-square framework is that when the experimental design is balanced, every product is equally weighted. This is more appropriate and leads to logical outputs as opposed to the usual chisquare framework that can lead to counterintuitive outputs. Indeed, it sounds more logical to weight the products equally and not rescale them according to their number of received citations when they have been evaluated the same number of times. Note that an equivalent weighting of the products using the usual chi-square framework is almost impossible since products are very unlikely to receive the same number of citations at the panel level. The multiple-response tests per cell introduced in this paper take into account all the specific aspects of FC and CATA data, especially the non-independence of citations between descriptors.

The conclusions provided by the two chi-square frameworks are not always necessarily different. For example, they would have been almost the same on the texture dataset if P5 had not been artificially added to the dataset. The more different the citation rates (all descriptors combined) between products are, the more the conclusions drawn from the usual chi-square framework will differ from the multipleresponse one. The products likely receive different numbers of citations when some products have few sensory characteristics while some others have a lot or when some products present obvious characteristics while the characteristics of the other products are more subtle; these kinds of situations are likely to occur in sensory evaluation.

Since the multiple-response chi-square framework relies heavily on Monte-Carlo and bootstrap simulations, the results of the proposed analyses are not instantaneous.

For the datasets used as examples, it took around 30 seconds by dataset to obtain the results of all analyses. However, this computation time increases with the number of evaluations and thus with the number of subjects and products. For large datasets (e.g. 3000 evaluations), it takes around 5 minutes to obtain the results using the settings of this paper.

Conclusion

For the analysis of Free-Comment and Check-All-That-Apply data, the paper proposes to replace the usual chi-square framework with a new multiple-response chi-square framework taking into account dependence among citations within an evaluation. It is thus statistically valid while the former was not. The new framework includes a test of dimensionality, a Correspondence Analysis with confidence ellipses, a test for pairwise product comparison, and a test of significance of product by descriptor associations. Note that ellipses, tests of product comparisons, and tests of association with descriptors are the three of them computed on the significant axes of dependence. The basic difference introduced by this new framework is not to longer consider citations (one descriptor for one product by one subject) as experimental units, but to rely on evaluations (vector of citations for one product by one subject) as being the experimental units. Simulations showed that testing the significance of product by descriptor associations on the significant axes of dependence increased power in detecting product by descriptor associations without any inflation of the type I error. The new approaches are supported by an R package called "MultiResponseR" and available upon request to the authors and on GitHub.

Appendix: Simulations

To investigate the benefits and/or the downsides of performing the multiple-response tests per cell on the derived contingency table corresponding to the significant axes, simulations of sensory data were performed. For every simulation, 60 subjects, 5 products, and 10 descriptors were considered. The 5 products were considered as being evaluated by the 60 consumers, as it is common in sensory evaluation. The descriptors marginal probabilities were randomly chosen and were the following: 0.20, 0.56, 0.26, 0.23, 0.21, 0.30, 0.20, 0.42, 0.52, 0.75. From these marginal probabilities, the matrix of expected probabilities under the null hypothesis of independence between products and descriptors was computed. This matrix contained 50 cells (5 products × 10 descriptors). Some deviation from independence was then added iteratively to these expected probabilities such that at each step, one axis of dependence was added orthogonally to the previous axes. On the first added axis, two products were differentiated on six descriptors. On the second added axis, two products were differentiated on four descriptors. On the third added axis, two products were differentiated on two descriptors. On the fourth added axis, four products were differentiated on four descriptors. This enabled to control the true dimensionality of the dependence between products and descriptors. The cells that deviated from the null hypothesis did with equal intensity but with opposite direction to keep the marginal probabilities fixed. Two levels of deviation intensity were considered: 0.1 and 0.2. 8 matrices (4 levels of dimensionality × 2 levels of deviation intensity) of probabilities were thus generated. Each of the 8 matrices contains 50 cells (5 products × 10 descriptors).

For each of these 8 matrices, 1000 datasets were simulated. Each of these datasets was generated by adding 60 individual data (the subjects). Each individual data was generated by performing a random Bernoulli draw for each of the 50 cells according to the specified probability given in the matrix.

For each of the 8000 datasets (8 matrices of probabilities × 1000 generated datasets), the number of significant axes was considered unknown and was determined using the dimensionality test presented in section 2.3.1. The multipleresponse tests per cell were then performed on either the observed table or the derived contingency table corresponding to the significant axes returned by the test.

The p-values of the multiple-response tests per cell were stored.

For each combination of the factors deviation intensity (0.1 or 0.2), dimensionality (one axis, two axes, etc.), and table (observed or derived) and for each of the 50 cells, the proportion of test (among the 1000 datasets) rejecting the null hypothesis was computed at the following nominal alpha risks: 5%, 7.5%, and 10%. Then, the results from a given cell were assigned either to the group H0 if its probability was not modified or to the group H1 otherwise. Finally, the average proportion of rejection of the null hypothesis was computed within each group (H0 or H1), number of dimensions, and deviation intensity. The results are presented in Table 3. Table 3: Average proportion of rejection of the null hypothesis among the 1000 simulations depending on the deviation intensity, 562 the dimensionality, the nominal alpha risk, the table considered, and the deviation from the null hypothesis or not.

Table 3 shows that the empirical type I error never exceed the nominal alpha risk in group H0 for both approaches, which suggests that both approaches are valid. It can be seen that the empirical type I error in the H0 group was even slightly lower when considering the derived table which is a nice feature.

The percentage of rejections in group H1 (estimating test power) was higher when considering the derived table as compared to the observed table whatever the combination of factors considered except with a dimensionality of 4. Therefore, performing the multiple-response tests per cell on the derived contingency table corresponding to the significant axes enables gaining power without increasing type I error. It should also be noted that the smaller the dimensionality of the dependence, the higher the gain of power. It is logical because a low dimensionality maximizes the difference between the derived table and the observed one. Finally, it should also be noted that the gain in power is higher with the lower independence deviation (0.1 vs 0.2), that is with the more complex/subtle situation. This is a nice feature arguing in favor of performing the multiple-response tests per cell on the derived contingency table corresponding to the significant axes.

2. 1 .

 1 Notations and multiple-response chi-square test of homogeneity Let us consider an FC or a CATA experiment where subjects evaluated products on descriptors. Each product ∈ 1, … , has been evaluated times and the total number of evaluations is equal to = ∑ . Note that in the particular case of balanced experimental design, i.e. when all subjects evaluated all products, then = × . Let us denote by the number of citations of descriptor ∈ 1, … , for product during the evaluations and by the number of citations of descriptor during all the evaluations.

  equal to the minimum between -1 and . Let us consider H U the matrix of the T -R + 1 last left singular vectors of C, I U the diagonal matrix of the T -R + 1 last singular values of C and J U the matrix of the T -R + 1 last right singular vectors of C such that C U = H U I U J U F . Let us denote by 1 23 4 U the multiple-response chi-square statistic of the derived contingency table corresponding to the T -R + 1 last axes of the MR-CA denoted W U and defined following the reconstitution formula as:

  number of contingency tables by randomly permuting the response vectors along products within each subject (ii) Perform MR-CA on each of the simulated contingency tables under permutation having an equal or a larger value than the observed 1 23 4 U .

  (a), Fig. 1(b), and Fig. 1(c).

Fig. 1 :

 1 Fig. 1: Biplot from Correspondence Analysis of the texture dataset: (a) usual CA (axes 1-2), (b) MR-CA (axes 1-2) and (c) MR-CA (axes 3-2).

Fig. 2 :

 2 Fig. 2: Biplot from Correspondence Analysis of the flavor dataset: (a) usual CA (axes 1-2), (b) usual CA (axes 3-2), (c) MR-CA (axes 1-2) and (d) MR-CA (axes 3-2).

Fig. 3 :

 3 Fig. 3: Descriptors by product percentages of citations across the panel. Highlighted cells denote a significant (α = 10%) Fisher exact tests per cell in the usual chi-square framework or a significant (α = 10%) multiple-response test per cell (2000 simulations) in the MR chi-square framework

Table 1 :

 1 Eigenvalues of Correspondence Analysis and corresponding p-values (in brackets) for testing the number of significant axes in the usual and multipleresponse frameworks for the two datasets Table

		-square framework	Axis 1	Axis 2	Axis 3	Axis 4
	Texture	Usual Multiple-response 0.907 (<0.001) 0.323 (<0.001) 0.079 (<0.001) 0.002 (0.6146) 0.447 (<0.001) 0.162 (<0.001) 0.001 (0.9970) 0 (1)
	Flavor	Usual Multiple-response 0.557 (<0.001) 0.089 (<0.001) 0.013 (0.0054) 0.243 (<0.001) 0.012 (0.0154) 0.003 (0.0914)	/ /

Acknowledgments

This study is part of a Ph.D. financed by the Region Bourgogne-Franche-Comté and the SensoStat Company.

The authors would like to thank Barry Callebaut© for providing the chocolate samples.

Calculations in the Appendix were performed using HPC resources from PSIUN CCUB (Centre de Calcul de l'Université de Bourgogne, France). The Institut de Mathématiques de Bourgogne is supported by the EIPHI Graduate School (contract ANR-17-EURE-0002).


M O Naumenko-Dézes

I M Villa

Y Rolland

S Gallet

P Lanari

reveals intragrain heterogeneity

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Three main issues are currently limiting the accuracy and precision of the 40 Ar/ 39 Ar dating technique and consequently its reliability: (1) the systematic uncertainty of the flux monitor ages (e.g. the Fish Canyon Tuff sanidine: [START_REF] Cebula | The Fish Canyon tuff, a potential standard for the 40 Ar/ 39 Ar and fission-track methods[END_REF][START_REF] Rivera | A refined astronomically calibrated 40 Ar/ 39 Ar age for Fish Canyon sanidine[END_REF], which can be due at least partly due to mineralogical inhomogeneities (e.g. [START_REF] Bachmann | The Fish Canyon magma body, San Juan volcanic field, Colorado: rejuvenation and eruption of an upper-crustal batholith[END_REF]; (2) the lack of understanding of the intra-grain 40 Ar/ 39 Ar age distribution;

(3) the uncertainty on the 40 K decay constant and its branching ratio [START_REF] Steiger | Subcommission on geochronology: convention on the use of decay constants in geo-and cosmochronology[END_REF][START_REF] Naumenko-Dèzes | 40 K-40 Ca -87 Rb-87 Sr age comparison: constraints on the 40 K decay constant[END_REF].

Issue (3) could, at least in principle, be solved by geological intercalibration: the age determination of a suitable calibrator sample by two or more geochronometric systems yields the unknown decay constant(s) as a function of one decay constant taken as a reference (usually that of 238 U), provided the calibrator records a "pointlike" geological history [START_REF] Begemann | Call for an improved set of decay constants for geochronological use[END_REF] and is homogeneous within the analytical resolution (JCGM, 2008, entry F.2.6.2). This is a long-standing problem inherent in intercalibration experiments [START_REF] Villa | IUPAC-IUGS status report on the half-lives of 238 U, 235 U and 234 U[END_REF][START_REF] Villa | IUGS-IUPAC recommendation on the half-lives of 147 Sm and 146 Sm[END_REF].

The metrological requirement reveals a contradiction with the natural heterogeneity of many natural samples and their un-point-like extensive retrograde reaction history. For the calibrators of the 40 K decay half-life the scarcity of chemically homogeneous samples had been addressed by Naumenko- [START_REF] Naumenko-Dèzes | 40 K-40 Ca -87 Rb-87 Sr age comparison: constraints on the 40 K decay constant[END_REF]. The two least ill-preserved samples analyzed in that study were two museum-grade mica megacrysts: the Phalaborwa phlogopite (ca. 2 Ga old), and the Rubikon lepidolite (ca. 0.5 Ga old), which had been dated by Rb-Sr and K-Ca. The Rb-Sr age of both samples (calculated with the 87 Rb decay constant recommended by Villa et al., 2015) agreed with the independently known U-Pb ages. This is strong evidence that the Rb-Sr system records no significant post-magmatic perturbation.

In contrast, the K-Ca ages, when calculated using the [START_REF] Steiger | Subcommission on geochronology: convention on the use of decay constants in geo-and cosmochronology[END_REF] decay constant (λ40) and branching ratio (B), were apparently younger than the U-Pb age. However, by analogy with the Rb-Sr system, the K-Ca system is expected to record the same lack of detectable chemical retrogression, as both systems feature alkali parents and earth alkali daughters. Therefore, Naumenko- [START_REF] Naumenko-Dèzes | 40 K-40 Ca -87 Rb-87 Sr age comparison: constraints on the 40 K decay constant[END_REF] had interpreted the apparent age discordance as a constraint for a (λ40, B) pair different from that proposed by [START_REF] Steiger | Subcommission on geochronology: convention on the use of decay constants in geo-and cosmochronology[END_REF]. An additional complication was that [START_REF] Dèzes | Refined isotopic compositions of K, Ca and a complementary comparison of the 40K-40Ca, 40K-40Ar and[END_REF] had obtained preliminary apparent K-Ar ages in disagreement with all of the above ages. One reason for the discordance could be if diffusive loss of Ar* had occurred and had been faster than that of Sr* and Ca*. Such an effect for the Madagascar sanidine had been hypothesized by [START_REF] Nägler | In pursuit of the 40 K branching ratio: K-Ca and 39 Ar-40 Ar dating of gem silicates[END_REF], and documented by [START_REF] Flude | Observation of centimetre-scale argon diffusion in alkali feldspars: implications for 40Ar/39Ar thermochronology[END_REF]. Another reason could have been compositional heterogeneities (near the resolution of EPMA analyses, approx. 1 %) resulting from dissolution/reprecipitation. The spatial distribution of chemical differences resulting from these two resetting mechanisms is distinctly different [START_REF] Villa | Diffusion in mineral geochronometers: present and absent[END_REF]: diffusion necessarily produces a bell-shaped concentration profile, whereas recrystallization usually results in patchy zonations.

This study endeavours to assess the above-mentioned limiting factors by exploring the spatial distribution of the radiogenic 40 Ar ( 40 Ar*) and 40 K concentrations using high precision 40 Ar/ 39 Ar dating on subgrain aliquots of the Phalaborwa and Rubikon micas studied by Naumenko- [START_REF] Naumenko-Dèzes | 40 K-40 Ca -87 Rb-87 Sr age comparison: constraints on the 40 K decay constant[END_REF]. Because both samples appear fresh, it was expected that 40 Ar/ 39 Ar ages would be regularly distributed in the grains and would thus make it possible to obtain constraints on the K-Ar and K-Ca double decay rate.

Analytical techniques

Phlogopite (Fig. 1) from the Phalaborwa layered intrusion, South Africa, was obtained from the Phalaborwa Mine Collection, courtesy K. Mezger. The Phalaborwa intrusion consists of carbonatites, phoscorite and pyroxenites. Each layer of this intrusion was dated with U-Pb; the concordant U-Pb ages indicate a rapid emplacement sequence [START_REF] Wu | In situ determination of U-Pb ages and Sr-Nd-Hf isotopic constraints on the petrogenesis of the Phalaborwa carbonatite Complex, South Africa[END_REF]. Its fast cooling was suggested as a recommendation to choose it as a 40 K half-life calibrator (Y. Amelin, pers. comm., 2014), following the assertion by [START_REF] Steiger | Subcommission on geochronology: convention on the use of decay constants in geo-and cosmochronology[END_REF] that intercalibrations are reliable if samples look fresh and cooled quickly. Five phlogopite sheets up to a cm size were hand-picked from phoscorite. The sheets were first split along cleavage planes with a steel blade so as to obtain a thin, transparent, visually inclusion-free layer, photographed, then cut into long strips ca. 1-2 mm wide. Strips were then cut into 0.2-0.4 mm wide rectangular chips with a mass of the order of 0.1 mg, recording and photographing the position of each of them within the original grain. The mica rectangles were individually packed into thin, high purity aluminum foil disks that were stacked in an aluminum tube interspaced with interlaboratory age monitors and with the Rubikon lepidolite (see below) for irradiation. The mass of the thin rectangular pieces was about 10 times larger than the typical UV Laser microprobe sample size; in order to have 37 Ar and 38 Ar beams large enough for a sufficiently precise measurement, some samples consisted of larger chips or several rectangles. Rb-Sr dating on another fragment of the same hand specimen yielded an age of 2056±5.2 Ma (Naumenko- [START_REF] Naumenko-Dèzes | 40 K-40 Ca -87 Rb-87 Sr age comparison: constraints on the 40 K decay constant[END_REF], concordant with the U-Pb and Pb-Pb ages of baddeleyite and zircon (2060.3±0.4 Ma, [START_REF] Kumar | Precise Pb baddeleyite geochronology by the thermal extractionthermal ionization mass spectrometry method[END_REF].

Lepidolite from the Rubikon complex, Namibia, was collected from a Li and Rb bearing pegmatite; a ca. 3 cm lepidolite foil was obtained from Schweizerisches Naturhistorisches Museum, Bern, courtesy B.A. Hofmann.

Existing age data for Rubikon are a baddeleyite U-Pb age of 505.5±2.7 Ma [START_REF] Melcher | Tantalum-(niobium-tin) mineralisation in African pegmatites and rare metal granites: Constraints from Ta-Nb oxide mineralogy, geochemistry and U-Pb geochronology[END_REF], coincident with the Rb-Sr age, 504.7±4.2 Ma (Naumenko- [START_REF] Naumenko-Dèzes | 40 K-40 Ca -87 Rb-87 Sr age comparison: constraints on the 40 K decay constant[END_REF]. The lepidolite foil was first cleaved and then cut into rectangular chips in the same way as the Phalaborwa phlogopite. The chips were photographed, wrapped and irradiated in the same package as the latter. The 270 MWh irradiation was performed at McMaster University, Canada, avoiding Cd shielding in order to preserve the petrological information deriving from 38 ArCl. For reasons not anticipated, many samples were only measured after 37 ArCa had decayed.

For our purpose of charting the geometrical distribution of 40 Ar* in the megacryst it is sufficient to use selfmonitoring, i.e., calculating the irradiation parameter J for one arbitrarily chosen chip (P97 and P67, respectively: Table 1) assumed to have a K-Ar age equal to the Rb-Sr age previously obtained by us on a different fragment of the same megacryst (Rubikon) or hand specimen (Phalaborwa). Self-monitoring circumvents the acute problem that the correct half-life of 40 K is unknown. By referring the relative variations of the 40 Ar*/ 39 Ar ratio to the correct Rb-Sr age, instead of referring them to the 40 Ar*/ 39 Ar ratio of an irradiation monitor having a different age, the problem of taking into account the unknown curvature of the 40 Ar*/ 39 Ar increase with time is reduced to approximating the growth curve with a linear increase for the time interval between 2000 and 2060 Ma, or between 490 and 505 Ma. Thus, we calculated the K-Ar ages using the [START_REF] Steiger | Subcommission on geochronology: convention on the use of decay constants in geo-and cosmochronology[END_REF] constants, but the relative ages among the chips remain the same (to less than 0.1 % bias) if we use any of the (λ40, B) pairs discussed by Naumenko- [START_REF] Naumenko-Dèzes | 40 K-40 Ca -87 Rb-87 Sr age comparison: constraints on the 40 K decay constant[END_REF]. If the mica were chronologically homogeneous, the apparent ages of all chips should be equal, after correcting for the smooth irradiation gradient. If the apparent ages of the pieces show an additional difference after correcting for the irradiation gradient, and this difference is due to diffusive Ar loss, then the position of the piece within the original crystal should define a bell-shaped diffusion profile. If the relationship between apparent age and position is patchy, the rejuvenation mechanism controlling the 40 Ar* and K distributions can only be due to an additional process faster than diffusion, such as retrograde reactions.

Both Rubikon and Phalaborwa were measured in two different laboratories, on two different mass spectrometers and using two different degassing protocols. In Milano Bicocca, samples were pre-heated to ca. 250 °C for 30 minutes under vacuum and then step-heated in a double-vacuum furnace and analyzed in a NuInstruments Noblesse mass-spectrometer equipped with one Faraday collector and two ion counters [START_REF] Bosio | Tephrochronology and chronostratigraphy of the Miocene Chilcatay and Pisco formations (East Pisco Basin, Peru)[END_REF]. The intercalibration of the three collectors was guaranteed in each measurement, for the instantaneous source conditions, by measuring at least two isotopes on two different collectors. In Géoazur, Nice, individual chips were mounted in a drilled copper plate in a steel chamber, maintained at 10 -8 -10 -9 mbar, not preheated, and stepheated with a 100 kHz IR-CO2 laser used at 5-15 % power output during 60 s. Argon isotopes were measured in the multi-collector ARGUS VI spectrometer with five Faraday detectors (10 11 Ω for 40 Ar, 10 12 Ω resistors for other masses). Collector gain calibration was performed by the computer-controlled application of predetermined voltages to each collector. Collector sensitivity was controlled every week by peak jumping measurements of air pipettes and corresponding correction was applied. Over the course of one year, the relative gain drift between collectors did not exceed 19 ppm. System blanks were measured before and after each measurement. The average of these two measurements was then subtracted from the measured sample value. The blanks typically accounted for less than 0.06% of the 40 Ar, 0.6-3.81% of the 38 Ar, and 0.12-0.71% of the 36 Ar. Mass discrimination for both mass spectrometers was monitored by daily analyses of air pipette volumes. Géoazur used ArArCalc© software v2.5.2 for data processing [START_REF] Koppers | ArArCALC-software for 40 Ar/ 39 Ar age calculations[END_REF], Milano Bicocca used a dedicated calculation spreadsheet. Spreadsheets with the complete analytical data from both laboratories can be found in A_Supplementary Materials, in which 40 Ar/ 39 Ar ages are calculated with fixed J (see above), but not corrected for the flux gradient. Ages in Table 1 are corrected for flux gradient.

Electron probe microanalyses (EPMA) were obtained on a JEOL 8200 electron microprobe at Institut für Geologie, Bern University, by mounting other, nearby cleavage fragments of the same megacryst parallel and perpendicularly to the (001) plane, the latter exposing the interior of the TOT phyllosilicate structure [START_REF] Heri | Petrology and geochronology of the "muscovite standard" B4M[END_REF]. The data of the EPMA analyses were calibrated with XMapTools (Lanari et al. 2014a;2019), which allows the precise determination of mineral species at the ca. 10 μm scale (e.g., Lanari et al., 2014b). This provides petrological constraints for the interpretation of the 40 Ar/ 39 Ar data. All analyzed elements (Mg, Mn, Al, Fe, Ti, Si, K) show resolvable concentration variations (Figs. 3,4,and B_Supplementary Materials).

In addition to the mica megacrysts, several interlaboratory monitors (Fish Canyon sanidine, Alder Creek sanidine, McClure Mountains amphibole) were also analyzed; mg-sized aliquots were wrapped in flat Al foil disks, whose positions were also recorded. Their chronological heterogeneity, as already documented by [START_REF] Jicha | Re-evaluation of the ages of 40 Ar/ 39 Ar sanidine standards and supereruptions in the western US using a Noblesse multi-collector mass spectrometer[END_REF] and Naumenko- [START_REF] Naumenko-Dèzes | 40 K-40 Ca -87 Rb-87 Sr age comparison: constraints on the 40 K decay constant[END_REF], was reproduced here. The role of the interlaboratory monitors in the present experiment was minor, other than to confirm the irradiation gradient of 0.6%/cm already inferred from the selfmonitoring.

Results

The apparent 40 Ar/ 39 Ar ages of the lepidolite chips span a 6 % interval from 504±10 Ma down to 473±5 Ma (Fig. 5, Table 1). It is important to point out that individual chip ages of Rubikon and Phalaborwa are not calibrated absolute ages; instead, they represent intra-crystal variations of the 40 Ar*/ 39 Ar ratio normalized to one arbitrary chip (P97 and P67, respectively). The principal result is the scatter of the apparent chip ages, not their absolute value. It could be that some chips show calibrated K-Ar ages higher than the Rb-Sr and U-Pb age (due to removal of K), but more probably the highest K-Ar ages should be equal, within uncertainty, to the Rb-Sr and U-Pb age. The age variations are neither a function of the position in the irradiation can nor of the position inside the crystal.

The apparent 40 Ar/ 39 Ar ages of the phlogopite chips vary from 2081±15 Ma down to 1957±13 Ma (Fig. 6, Table 1). The age variations neither correlate with the position of the chips in the crystal nor with the position of the chips in the irradiation can (Fig. 7). The darker colored sheets showed ages reset below 2010 Ma. Two of the most clean and thin sheets yielded ages close to the 2060 Ma age of the intrusion. The spatial age variation is random, a 2050 Ma chip may lie next to a 2005 Ma one. As the samples were analyzed in two laboratories, the comparison confirms that the observed age distribution is not an analytical artifact but a primary effect.

Both micas in both laboratories gave invariably "plateau" ages (Fig. 8 shows one example of two neighboring chips irradiated in adjacent positions). "Plateaus" had often been interpreted in the past as a sign of ideality. The "plateau age" variations largely exceeding the 2 sigma level imply instead that at most one of the "plateaus" can reflect ideal behaviour, whereas all others do not. This reproduces results obtained by [START_REF] Foland | 40 Ar/ 39 Ar incremental heating plateaus for biotites with excess argon[END_REF] and [START_REF] Popov | Removing a mask of alteration: Geochemistry and age of the Karadag volcanic sequence in SE Crimea[END_REF], calling into question the identification of "plateaus" with sample ideality. Therefore, the questions should be: what brings about the observed variations, what diagnostic tools does the 40 Ar/ 39 Ar systematics provide, what correction (if any) can be applied, and what do our findings imply for the viability of micas as precise half-life calibrators?

The observed chronological patchiness of the two museum megacrysts required further petrological investigations. Electron probe microanalysis of the mica mounts showed clear variations in major element concentrations (Figs. 34): a more-than-doubling of TiO2 (from 0.6 to 1.6 %) and relative increases around 20 % of FeO (between 7.5 and 9.5 %) and MnO (between 1.8 and 2.1%). Some schlieren-like variations of K, Al and Mn concentrations were observed in the (001) plane of the lepidolite (Fig. 3), unrelated to visible microfractures. The variations in the plane perpendicular to (001) were even more revealing. Cleavage gaps (possibly enhanced by the polishing) are obvious, gaps being between 3 and 25 µm thick; no gradual transition between the crystal and the gaps (where all element concentrations are zero) can be recognized. In addition, near the upper margin of the element map the lepidolite shows a prominent curved vertical feature rich in Si and poor in Al and Mn. It also shows an entire layer poorer in K and Al and slightly richer in F and Mn. The phlogopite (Fig. 4) shows the same (001)-parallel cleavages with lower element concentrations in the gaps; however, unlike lepidolite, Ti, Fe and Al (more insoluble elements) do not fall to zero like Mg and Si do. It appears that at least some cleavages were partly filled. This is more evident near the bottom of the element map, where in a 20 µm-thick layer Mg reaches the double of the concentration it has in the rest of the map. In the plane parallel to (001), two microfractures are visible, one of them being totally unrelated to the 150 % variation of the Ti concentration.

We interpret these observations as evidence that micas were affected by two processes: (1) fluid-assisted, open-system chemical reactions causing cation transport and element redistribution; (2) retrograde mineral formation. Some micas show signs of thin layers of vermiculite and/or pyrophyllite and/or talc developing in between mica sheets, even if initially a visual inspection by optical microscopy failed to reveal it. After EPMA element mapping had demonstrated mineral replacement, a dedicated re-investigation under the optical microscope did reveal some fuzzy darker bands that can reach 100 µm width. As the element maps show that alteration is patchy at the mm scale, it is easy to explain why 100 µm sized chips that were randomly cut across alteration patches showed haphazard ages. In the brighter and cleaner-looking sheets the alteration bands are < 5 µm wide. Alteration fronts are more likely to start along cracks, which favor the ingress of metasomatic fluids, which would cause younger ages near cracks. However, the presence of narrow (5 µm wide) bands with anomalous element concentrations in the middle of mica sheets shows that the chemically open-system alteration followed three-dimensional pathways and did affect even chips that had seemed unaltered under an optical microscope. At a scale ≈ 1 mm, the heterogeneities are averaged away and only a slightly younger whole-grain plateau is recorded. here. This agrees with the common observation that newly-formed crystals in presence of fluids do not retain any excess or inherited argon (e.g., [START_REF] Sanchez | Dating low-temperature deformation by 40 Ar/ 39 Ar on white mica, insights from the Argentera-Mercantour Massif (SW Alps)[END_REF].

We sought the causes of the age discordance amongst the samples for those chips whose 38 ArCl signal was sufficiently well-resolved from the background, i.e., those nine with a sample mass > 0.18 mg. We also obtained a detectable 37 Ar signal on three samples measured within two months of irradiation. The Ca/K ratios are highest in the low temperature steps, which we attribute to heterochemical retrogression phases, and decrease by 2-3 orders of magnitude in the steps representing the degassing of the micas sensu stricto. Two Phalaborwa samples had identical Ca/K ratios of 0.0018 ± 0.0003 (P14) and 0.0017 ± 0.0009 (P32) in the isochemical steps (the steps in which the Cl/K ratio (and the Ca/K ratio, where available) shows a constant chemical signature); the isochemical steps of one Rubikon sample averaged Ca/K = 0.0031 ± 0.0008. In the three studied instances, the Ca/K ratios of the isochemical steps show no correlation with the corresponding ages.

Where the 38 ArCl and 37 ArCa are well resolved from background we can investigate isochemical ages (as defined by [START_REF] Villa | Geochronology and isotope transport systematics in a subsurface granite from the Larderello-Travale geothermal system (Italy)[END_REF], i.e., the average of ages in the isochemical steps; normally isochemical steps yield indistinguishable ages [START_REF] Villa | Age discordance and mineralogy[END_REF], as they do here. The steps with a constant Cl/K ratio in the Phalaborwa and Rubikon samples are those with the lowest Cl/K. In contrast, the first step degassed at low temperature always showed high Cl/K, 1-2 orders of magnitude higher than the isochemical steps. The age of the first step is generally younger. The negative correlation between age and Cl/K ratio is statistically significant at the 95 % confidence level only in few cases (C_Supplementary Material), even if it appears ubiquitous. Because of this correlation trend, in principle the most reliable estimate of the K-Ar age is not given by the "plateau" age but by the leftmost (Cl-poorest) end-point of the Cl/K-age alteration trend, even if in practice the two coincide. However, clear as the petrochronological evidence may be, neither the Phalaborwa phlogopite nor the Rubikon lepidolite is suitable as decay constant calibrators for the K-Ar system, as the extrapolated correction of their measured apparent ages to model the "pre-alteration" K-Ar ages entails possible systematic errors exceeding 1 %.

By inference, this conclusion should be applied to all cases where no conclusive proof of the absence of retrogression was provided. The stratigraphic and geodynamic modeling at the km scale should always take into account the physical and chemical phenomena occurring at the sub-μm scale [START_REF] Villa | From the nm to the Mm: isotopes, atomic-scale processes, and continent-scale tectonic models[END_REF][START_REF] Hochella | Nanogeoscience: From origins to cutting-edge applications[END_REF]Ju et al., 2020, and references therein).

Conclusions

This study demonstrates a previously neglected phenomenon: even "rapidly cooled", clear, visually unaltered, museum-grade micas can be incompletely pristine geochronometers. In particular, small (5 µm wide) alteration layers between mica sheets already can affect K-Ar ages by several percent if the analyzed chip is small enough. Large variations of the Ti and Mg concentrations are a sign of retrograde reactions that correspond to the large K-Ar age variations.

Step heating "plateaus", as were obtained for all the samples in this study, do not guarantee that the K-Ar system was not partially reset. The spatial distribution of the apparent ages of subgrain chips does not support the idea that K-Ar ages only represent cooling ages. The present study documents subtle, patchy open-system compositional changes, which were detected by electron microprobe. This chemical heterogeneity requires one or more episodes of mass transfer mediated by post-magmatic aqueous fluid circulation.

The diagnostic tools that the 40 Ar/ 39 Ar systematics provide are a judicious use of the isotope correlation plots, especially the Cl/K-age and Ca/K-age diagrams. Identifying trends and extrapolating them to an end-member composition can provide a correction for the effects of metasomatism and alteration, but the precision of the extrapolation is limited by the precision of the electron probe microanalysis that constrains the pristine chemical composition of the unaltered patches. The potential to use the intercomparison of mineral ages as a calibration of half-lives suffers from the same limitation on precision.

Figure captions

  The observation that Rb-Sr ages are close to the U-Pb ages was used by Naumenko-[START_REF] Naumenko-Dèzes | 40 K-40 Ca -87 Rb-87 Sr age comparison: constraints on the 40 K decay constant[END_REF] to hypothesize that the geochemically similar K-Ca parent-daughter pair likely also should give an age concordant with the U-Pb age. From the assumed K-Ca age and the measured40 Ca*/ 40 K ratio they constrained the half-life and branching ratio of the 40 K decay. The closed behaviour of the Rb-Sr system despite the introduction of Mg and Ti (which requires at least local recrystallization of the T-O-T mica framework hosting Mg and Ti) suggests that the earth-alkali daughter isotopes 87 Sr* and 40 Ca* are reprecipitated into the recrystallizing vermiculite/talc pseudomorph, whereas incompatible 40 Ar* daughter isotopes are lost, giving the rejuvenated K-Ar ages presented

Figure 1 .

 1 Figure 1. Thin-section pictures of the Phalaborwa sample. Left, plane polarized light, right, cross polarized light.

Figure 2 .

 2 Figure 2. Scheme of the irradiation package. Green, irradiation monitors; PhB, phlogopite; Lep, lepidolite; P1-P102, individually packed rectangular mica chips. Individual chip positions in the stack increased by approx. 0.5 mm relative to the previous one.

Figure 3 .

 3 Figure3. Two examples of electron probe element maps of the Rubikon lepidolite, parallel and perpendicular to the (001) plane. Data treated with XMapTools, versions 2018 and 2020(Lanari et al. 2014a; 2019). More element maps can be found in the B_Supplementary Materials.

Figure 4 .

 4 Figure 4. Two examples of electron probe element maps of the Phalaborwa phlogopite, parallel and perpendicular to the (001) plane. Data treated with XMapTools(Lanari et al. 2014a; 2019). More element maps can be found in the B_Supplementary Materials.

Figure 5 .

 5 Figure 5. Apparent 40 Ar/ 39 Ar ages of Rubikon lepidolite chips, normalized to chip P97 assumed to have an age equal to the Rb-Sr age.

Figure 6 .

 6 Figure 6. Apparent 40 Ar/ 39 Ar ages of Phalaborwa phlogopite chips, normalized to chip P67 assumed to have an age equal to the Rb-Sr age.

Figure 7 .

 7 Figure 7. Normalized (self-monitored) age vs. position of rectangular chips in the irradiation canister.

Figure 8 .

 8 Figure 8. "Plateau" ages for two neighboring chips of the Phalaborwa phlogopite. The height of the boxes is the 2 sigma uncertainty of each step.

Acknowledgements

Samples were donated by K. Mezger and B.A. Hofmann. Funding for M.O.N.-D. were provided by the Schweizerische Nationalfonds, grants 200021-131916 and 200020-153126. Analyses in Nice Sophia Antipolis were supported by a Géoazur internal fund. Valentina Barberini (Milano Bicocca) is thanked for laboratory maintenance. Reviews by Anonymous, D.V. Popov and B.S. Kamber greatly improved the manuscript.