Antoine Barbieri
email: antoine.barbieri@umontpellier.fr

Jean Peyhardi

Thierry Conroy

Sophie Gourgou

Christian Lavergne

Caroline Mollevi

Item response models for the longitudinal analysis of health-related quality of life in cancer clinical trials

Keywords: Item response theory, Mixed models, Ordinal categorical data, Longitudinal analysis, Health-related quality of life

Introduction

Endpoints refer to biological and clinical measures to assess the efficiency of new therapeutic strategies. The overall survival endpoint is the gold standard to show a clinical benefit of these strategies and treatments. Therapeutic treatments being more efficient and increasing the patients' lifetime, the overall survival endpoint may become insufficient to show a significant difference between two treatments. It is then necessary to consider a longer follow-up or a larger cohort of patients to have a sufficient number of events and a good statistical power [START_REF] Fiteni | Endpoints in cancer clinical trials[END_REF], both representing considerable costs. Thus, to conclude to the benefit of a new treatment, other endpoints have emerged and the health-related quality of life (HRQoL) is currently one of the most important. In cancer clinical trials, the patient-reported outcomes are increasingly used to analyze a clinical benefit for medical decision-making [START_REF] Fiteni | Endpoints in cancer clinical trials[END_REF]. The HRQoL endpoint may seem more pertinent to show the interest of a new therapy in some cases such as the palliative or geriatric situations. However, there are conceptual and methodological brakes underlying to the concept and the assessment of HRQoL. Indeed, HRQoL is a multidimensional concept regarding the physical, psychological and social functions as well as symptoms associated with the disease and treatments. Another conceptual brake is the subjectivity of its measurement. Indeed, patients report their feelings about their HRQoL thanks to self-reported questionnaires. Both arguments preclude the use of HRQoL as sole primary endpoint in clinical trials.

In oncology, HRQoL is assessed using a general questionnaire for a set of different cancers, and an additional specific questionnaire associated with each type of cancer [START_REF] Aaronson | The European Organization for Research and Treatment of Cancer QLQ-C30: a quality-of-life instrument for use in international clinical trials in oncology[END_REF][START_REF] Cella | The Functional Assessment of Cancer Therapy scale: development and validation of the general measure[END_REF]. Each questionnaire decomposes the HRQoL to measure several under-concepts (dimensions of HRQoL) which themselves comprise one or several items. The items are built on the Likert scales in which the response variable is ordinal. Thus, considering several items for a given dimension, HRQoL data are composed of multiple ordinal responses. Also, the questionnaires are filled by the subjects themselves, and collected at different times defined in the trial protocol (usually at inclusion, during treatment and follow-up). These repeated measures are used to assess the evolution of the subject's HRQoL over time. In Europe, these questionnaires are developed and validated by the European organization for research and treatment of cancer (EORTC). The standard questionnaire currently used in oncology is the EORTC Quality of Life Questionnaire -Core 30 (EORTC QLQ-C30) [START_REF] Aaronson | The European Organization for Research and Treatment of Cancer QLQ-C30: a quality-of-life instrument for use in international clinical trials in oncology[END_REF], together with the scoring procedure proposed by the EORTC [START_REF] Fayers | EORTC QLQ-C30 Scoring Manual[END_REF]. The score is then calculated for each dimension and for each subject, corresponding to the average of the item responses for a single dimension, and expressed on a scale ranging from 0 to 100. The interpretation is such that high functional scores reflect good functional capacities and a good HRQoL level, and conversely, high symptomatic scores represent strong symptoms and point out difficulties. The use of scoring procedures is common in practice because the statistical methods for quantitative variables are more powerful and easier to implement and interpret [START_REF] Grilli | Multilevel Models for Ordinal Data[END_REF]. But, in a Likert scale, the gap which separates each adjacent category of response ("not at all", "a little", "quite a bit" and "very much") may not be the same, and the HRQoL score calculation does not take into account this characteristic. Another drawback in the score use is that subjects could have different item outcomes and obtain the same score. In this situation, the score does not make a distinction between these subjects [START_REF] Gorter | Why item response theory should be used for longitudinal questionnaire data analysis in medical research[END_REF].

The longitudinal statistical models classically used in oncology are performed on the summary score through using the linear mixed models (LMM) or time-to-event models [START_REF] Anota | Comparison of three longitudinal analysis models for the health-related quality of life in oncology: a simulation study[END_REF]. In the LMM, the variable associated with the HRQoL score is considered as a Gaussian variable while it presents the characteristics of an ordinal variable, being non-continuous and bounded. These models allow to take into account the correlation introduced by repeated measurements on the same patient (collection of the HRQoL questionnaires over time) and different covariates such as time, treatment group, age... However, the use of the LMM for HRQoL analysis is scientifically questionable given the characteristics of the score. Furthermore, many symptomatic dimensions are composed of only one item, the HRQoL score has exactly the same properties than ordinal categorical data, and using the LMM is not appropriated. Thus, if a ceiling or floor effect is observed, the categorical feature is even more marked when one of the two extreme categories is over-represented. The second approach for the longitudinal analysis of HRQoL is based on the time-to-event models: the time-to-deterioration (TTD) and the time-untildefinitive-deterioration (TUDD) [START_REF] Hamidou | Time to Deterioration in Quality of Life Score as a Modality of Longitudinal Analysis in Patients with Breast Cancer[END_REF]. Survival approaches are often used and thus well-known in the oncology, and are appreciated for their easiness to interpret result and their good understanding by clinicians. In these models, an event is classically defined by the (definitive or not) deterioration of the HRQoL score between baseline and a follow-up time, given a minimally clinically important difference [START_REF] Anota | Time to health-related quality of life score deterioration as a modality of longitudinal analysis for health-related quality of life studies in oncology: do we need RECIST for quality of life to achieve standardization?[END_REF]. The lack of homogeneity of the methods used for the HRQoL data analyses in different oncology clinical trials is also a real obstacle to the comparison of results. Indeed, the LMM and TUDD approaches show results which may sometimes seem contradictory, and with different interpretations, but they may also be complementary. An example can be taken comparing two similar cancer clinical trials investigating the effect of bevacizumab. In the first trial [START_REF] Chinot | Bevacizumab plus Radiotherapy-Temozolomide for Newly Diagnosed Glioblastoma[END_REF], HRQoL analysis through TUDD showed that the bevacizumab group had a later deterioration of HRQoL compared with patients in the standard group. Conversely, in the second trial [START_REF] Gilbert | A Randomized Trial of Bevacizumab for Newly Diagnosed Glioblastoma[END_REF], HRQoL analysis using the LMM showed a worse HRQoL overtime in the bevacizumab group.

Interest in the HRQoL endpoint is growing rapidly in cancer clinical trials and it is essential to use a suitable methodology to analyze HRQoL data, taking into account the data properties (repeated measurements of the multiple ordinal responses). In our study, we first focused on the different and most adapted models to analyze HRQoL from raw data, i.e. directly on the item outcomes. Studies on psychometric properties from questionnaires such as the one used for HRQoL have been ongoing for a long time [START_REF] Edelen | Applying item response theory (IRT) modeling to questionnaire development, evaluation, and refinement[END_REF][START_REF] Jafari | Using Rasch rating scale model to reassess the psychometric properties of the Persian version of the PedsQLTM 4.0 Generic Core Scales in school children[END_REF], known as the item response theory (IRT). The IRT models link the individual's item responses and the latent variable which represents the studied HRQoL concept. They can be seen as generalized linear mixed models (GLMM) for ordinal responses with a particular parameterization of the linear predictor. The interest for this kind of model to analyze the data, including the longitudinal analyzes, is increasing [START_REF] Titman | Item response theory and structural equation modelling for ordinal data: Describing the relationship between KIDSCREEN and Life-H. Statistical Methods in Medical Research[END_REF][START_REF] Hardouin | Power and sample size determination for group comparison of patient-reported outcomes using polytomous Rasch models[END_REF][START_REF] Gorter | Why item response theory should be used for longitudinal questionnaire data analysis in medical research[END_REF][START_REF] Santos | Multidimensional and longitudinal item response models for non-ignorable data[END_REF]. However, to our knowledge, there is no work that discusses of the choice of one of the different IRT models over the others, and specially for HRQoL longitudinal analysis. First, we propose in section 2 a conceptual selection of these models through practical and methodological arguments. In this selection, we replace IRT models in the GLMM framework using the new specification of generalized linear models (GLM) for categorical responses, proposed by [START_REF] Peyhardi | A new specification of generalized linear models for categorical responses[END_REF]. Then, to complete the comparisons already done between the IRT models and the LMM on their capacity to detect fixed effects, we focus in section 3 on the sensibility of these models to detect the random effects through a simulation study. Finally, section 4 presents an application of the chosen IRT model on real data from a multicenter randomized phase III clinical trial in first-line metastatic pancreatic cancer patients [START_REF] Conroy | FOLFIRINOX versus gemcitabine for metastatic pancreatic cancer[END_REF].

Conceptual selection of IRT models

This section concerns a conceptual selection of IRT models for the longitudinal analysis of HRQoL in cancer clinical trials. HRQoL raw data are repeated measurements of ordinal multiple responses. The GLMM for ordinal responses seem well suitable to analyze this kind of data. The use of random effects takes into account the inter-patient variability and the correlation between the repeated measurements for each single patient. The IRT models have been increasingly used to analyze health data deriving from self-questionnaires made of polytomous responses [START_REF] Hardouin | Spatio-temporal Rasch analysis of quality of life outcomes in the French general population. Measurement invariance and group comparisons[END_REF][START_REF] Anota | Comparison of three longitudinal analysis models for the health-related quality of life in oncology: a simulation study[END_REF][START_REF] Barbieri | Applying the Longitudinal Model from Item Response Theory to Assess the Health-Related Quality of Life in the PRODIGE 4/ACCORD 11 Randomized Trial[END_REF]. These models turn out to be GLMM for polytomous data with a specific parameterization of the linear predictor taking into account the multiple outcomes. For ordinal responses, three families of regression models are described: the families of adjacent models [START_REF] Masters | A rasch model for partial credit scoring[END_REF][START_REF] Agresti | Analysis of Ordinal Categorical Data[END_REF], cumulative models [START_REF] Samejima | Estimation of Latent Ability Using a Response Pattern of Graded Scores1[END_REF][START_REF] Mccullagh | Regression models for ordinal data (with discussion)[END_REF], and sequential models [START_REF] Tutz | Sequential item response models with an ordered response[END_REF][START_REF] Fahrmeir | Multivariate Statistical Modelling Based on Generalized Linear Models[END_REF]. Many IRT models are proposed for the analysis of this kind of data, often with no explanation regarding the choice of one model over another.

In this section, we use the new specification of the GLM for categorical data proposed by [START_REF] Peyhardi | A new specification of generalized linear models for categorical responses[END_REF] to talk about the relevance of the models adopted in the context of the HRQoL longitudinal analysis in cancer clinical trials. Each model is defined according to three components (r,F,Z): the ratio of probabilities (r), the cumulative distribution function (F ), and the parameterization of the linear predictor determined by the design matrix (Z). For the GLMM framework, we extended this new specification to the quadruplet (r,F,Z,U) with Z the design matrix of fixed effects and U the design matrix of random effects. The relationship between these components is determined by R = F(Zβ + U ξ). Given the linear predictor η = Zβ + U ξ and π

(j) iv = (π (j) iv0 , . . . , π (j)
iv,M -1 ) the truncated vector 1 of conditional probabilities with π (j) ivm = Pr Y (j) iv = m|ξ i the conditional probability that subject i(i = 1, . . . , n) selects the category m ∈ {0, . . . , M j } for item j(j = 1, . . . , J) at visit v(v = 1, . . . , V ) given individual random effect, we defined:

• R = r m π (j) iv i,j,v,m ; • F η (j) ivm i,j,v,m = F η (j) ivm i,j,v,m
where η = η

(j) ivm i,j,v,m
.

After a presentation of the IRT parameterization used concerning the linear predictor, we compare different polytomous IRT models on the basis of the link function (ratio of probabilities and the cumulative distribution function (CdF)) using methodological and practical arguments.

The IRT parameterization of the linear predictor

The IRT probabilistic models emerged following the works of Georg Rasch [START_REF] Rasch | Probabilistic Models for Some Intelligence and Attainment Tests[END_REF] on dichotomous responses, and were then extended to ordinal responses. Considering the three families of adjacent, cumulative and sequential models, there are three associated famous IRT models [START_REF] Boeck | Explanatory Item Response Models: A Generalized Linear and Nonlinear Approach[END_REF][START_REF] Bacci | A Class of Multidimensional Latent Class IRT Models for Ordinal Polytomous Item Responses[END_REF], respectively the graded response model [START_REF] Samejima | Estimation of Latent Ability Using a Response Pattern of Graded Scores1[END_REF], the (generalized) partial credit model [START_REF] Masters | A rasch model for partial credit scoring[END_REF][START_REF] Muraki | A Generalized Partial Credit Model: Application of an EM Algorithm[END_REF], and the sequential model [START_REF] Tutz | Sequential item response models with an ordered response[END_REF]. These models link the individual's item responses to the unidimensional latent variable which represents a concept not directly measurable. In oncology framework, the concept is HRQoL relatively to one specific HRQoL dimension.

From the IRT, the specific parameterization of the linear predictor η

(j)
im is built into two parts: the individual part and the item part. The best-known way is to consider the following decomposition:

η (j) im = α j (θ i -δ jm ) , (1)
where θ i is associated with an unidimensional random variable (currently assuming to be standard normal for identifiability) representing the latent value for the subject i, δ jm and α j being the item parameters.

1 This truncated vector is sufficient to characterize the categorical distribution since π

(j) ivM = 1 -M -1 c=0 π (j) ivc ;
Generally called difficulty parameter, δ jm is the intercept (or threshold) associated with the item j for the category m ∈ {1, . . . , M j }. The parameter α j is called the discrimination parameter of item j, and represents the sensitivity of each response probability according to the value of the latent trait. Indeed, the more the discrimination parameter value is high, the more the item allows well discriminating two individuals with a close latent trait value. However, the predictor is no longer linear for IRT models using discrimination parameters because it includes a product of parameters. Thus, these models do not belong to the class of GLMM [START_REF] Liu | A mixed-effects regression model for longitudinal multivariate ordinal data[END_REF].

In oncology, HRQoL analysis is carried out using IRT models which do not include the discrimination parameters (fixed to one for all items). Thus, these IRT models are within the class of GLMM. Concerning the longitudinal analysis, several studies proposed to extend some IRT models using the linear decomposition of the latent variable θ with fixed and random effects [START_REF] Hardouin | Spatio-temporal Rasch analysis of quality of life outcomes in the French general population. Measurement invariance and group comparisons[END_REF][START_REF] Verhagen | Longitudinal measurement in health-related surveys. A Bayesian joint growth model for multivariate ordinal responses[END_REF][START_REF] Huber | Mathematical Methods in Survival Analysis, Reliability and Quality of Life[END_REF]:

θ iv = x iv β + u v ξ i , (2)
with the vector β is associated with the fixed effects, the vector ξ i with the subject-specific random effects and the index v the current visit. In the equation ( 2), θ iv is thus the estimation of latent process at the visit v.

The probability ratio: structure of the models

The linear predictor is not directly related to the response probability but to a particular transformation ratio. The choice of the ratio is related to the nature of response from the ordering assumption among categories. Thus, reference ratio [START_REF] Peyhardi | A new specification of generalized linear models for categorical responses[END_REF] for nominal response is excluded in this work because the HRQoL responses are ordinal. Let's consider the simple situation from GLM with one item with (M + 1) response categories. The three model families for ordinal data are distinguished by the choice of the ratio of probabilities r (π) = (r 0 (π) , . . . , r M -1 (π)). Each model is summarized by M equations {r m (π) = F (η m )} m=0,...,M -1 with η m = δ m -θ, highlighting the decomposition of the link function which is determined through the ratio of probabilities and the CdF. Indeed, we may distinguish different ratios of probabilities for these different families, respectively, for the cumulative models,

r m (π) = π 0 + . . . + π m , m = 0, . . . , M -1; (3)
for the adjacent models,

r m (π) = π m π m + π m+1 , m = 0, . . . , M -1; (4)
and, for the sequential models,

r m (π) = π m π m + . . . + π M , m = 0, . . . , M -1.
In the IRT, adjacent and cumulative models are usually presented given the reverse permutation [START_REF] Samejima | Estimation of Latent Ability Using a Response Pattern of Graded Scores1[END_REF][START_REF] Masters | A rasch model for partial credit scoring[END_REF][START_REF] Bacci | A Class of Multidimensional Latent Class IRT Models for Ordinal Polytomous Item Responses[END_REF]. This permutation is defined as the reversal of category order [START_REF] Mccullagh | Regression models for ordinal data (with discussion)[END_REF]. Assuming that the considered CdF is symmetric (i.e. the corresponding probability density function is symmetric about the y-axis), these models are invariant under this permutation [START_REF] Peyhardi | A new specification of generalized linear models for categorical responses[END_REF]. For our application context, this is as an advantage for result interpretation. A lower item-response category reflects a lower level of the symptomatic dimensions whereas it represents a higher level of capacity for the functional dimensions. The reverse permutation for the functional dimensions makes it easier and intuitive for clinicians to present their results. This allows the homogenization of the result interpretation as it is done in the scoring procedure proposed by the EORTC (for functional dimensions, the score scale is reversed compared with the item responses categories order) [START_REF] Fayers | EORTC QLQ-C30 Scoring Manual[END_REF]. However, sequential models correspond to process ordering and reversing the process may change its nature. These models are not reversible (i.e. no invariant under the reverse permutation). Thus, sequential models will not be used and only the adjacent and cumulative models which correspond to scale ordering as used for HRQoL measurements, will be consider.

The cumulative models also have additional properties [START_REF] Mccullagh | Regression models for ordinal data (with discussion)[END_REF], including that they are invariant when successive categories are gathered. Thus, if one category is not observed, it can be combined with its successive categories without changing the model. Another advantage of the cumulative models is their interpretation through a continuous latent variable. Indeed, the continuous latent variable Ỹ underlying the model exists and a direct link with the response variable Y through the thresholds presumed to be strictly increasing (-

∞ = δ 0 < δ 1 < . . . < δ M < δ M +1 = +∞) is such as: {Y = m} if δ m < Ỹ ≤ δ m+1 , m = 0, . . . , M ,
where Ỹ = θ + ε and ε is the error term distributed following the CdF. Here, the latent variable Ỹ represents HRQoL and its interpretation is then equivalent to the interpretation of the response variable using a LMM. However, an advantage of the adjacent models is that there is no constraint affecting the model estimation. Nonetheless, the cumulative models have to respect constraints, which can make difficult the model estimation, particularly in the case of non-proportional design of linear predictor [START_REF] Peyhardi | A new specification of generalized linear models for categorical responses[END_REF]. For the proportional design, a common slope (θ) is considered for all categories, else the slope is dependent of the category (θ m ). Let the simple parameterization of the linear predictor η m = θ m -δ m for m ∈ {1, . . . , M } where θ m and δ m are the slope and the intercept associated with the category m, respectively. Considering proportional design (θ = θ 1 = . . . = θ M ), the cumulative models refer to the principle of thresholds [START_REF] Mccullagh | Regression models for ordinal data (with discussion)[END_REF][START_REF] Hedeker | A random-effects ordinal regression model for multilevel analysis[END_REF] with the constraint they have to be strictly increasing such as -∞ < δ 1 < . . . < δ M < +∞. Considering the non-odd proportional models, the constraint then becomes -∞ < η M < . . . < η 1 < +∞ which is more difficult to verify. For the longitudinal analysis of HRQoL data in oncology, the proportional design is considered and to verify the constraint only on the threshold is easier.

Table 1 summarizes some characteristics of these three families of models which are important for the longitudinal analysis of HRQoL in cancer clinical trials. In this context, a proportional design of the linear predictor is preferred. Under this parameterization, there is few difficulties to respect the cumulative models constraints and to estimate them. The adjacent models seem more flexible than cumulative models because they are always define for all linear predictor. But, their interpretation of the results using the cumulative model is more intuitive than adjacent models.

The cumulative distribution function

The latest component discussed in the IRT model selection is the CdF. Each model probability can be defined with any CdF. As commonly seen in the IRT models, and thank to the reversibility property, the adjacent and cumulative models are used in descending order. For the cumulative model, the probabilities are defined from the equation ( 3) and given the F as:

   π 0 = 1 -F (η 1 ) π m = F (η m ) -F (η m+1 ) , m = 1, M -1 π M = F (η M )
.

(5)

Then, general expression of the sequential model whatever the CdF used is such that

π m = F (η m ) m-1 k=1 {1 -F (η k )} ,
where m = 1, . . . , M and 0 k=1 {.} = 1 [START_REF] Fahrmeir | Multivariate Statistical Modelling Based on Generalized Linear Models[END_REF]. Such general equations have never been presented for adjacent models, only described with the logistic CdF [START_REF] Masters | A rasch model for partial credit scoring[END_REF][START_REF] Muraki | A Generalized Partial Credit Model: Application of an EM Algorithm[END_REF][START_REF] Fahrmeir | Multivariate Statistical Modelling Based on Generalized Linear Models[END_REF][START_REF] Agresti | Analysis of Ordinal Categorical Data[END_REF][START_REF] Hardouin | Spatio-temporal Rasch analysis of quality of life outcomes in the French general population. Measurement invariance and group comparisons[END_REF][START_REF] Anota | Comparison of three longitudinal analysis models for the health-related quality of life in oncology: a simulation study[END_REF]. However, the different response probabilities can be presented from the adjacent ratio and according to a general CdF (F ):

               π 0 = 1 1 + M m=1 m k=1 F (η k ) 1-F (η k ) π m = m k=1 F (η k ) 1-F (η k ) 1 + M m=1 m k=1 F (η k ) 1-F (η k ) , m = 1, . . . , M (6)
The CdF choice is especially used to best fit the data. Let's four CdF from two different kinds: the most commonly used symmetric distributions, the logistic and Gaussian distributions, and the two asymmetric distributions, the Gumbel min and Gumbel max distributions. The two later distributions are respectively defined by F (η) = exp(-exp(-η)) for the Gumbel max distribution and by F (η) = 1 -exp(-exp(η)) for the Gumbel min distribution.

Figure 1a shows different slopes depending on the particular CdF. The CdF allows to take into account the influence of linear predictor (η) change on the response probability evolution. In general IRT parameterization (equation 1), the slope adjustment is managed by the discrimination parameter. Depending on different discrimination parameter values, Figure 1b presents the CdF logistic according to the individual latent variable. This item parameter has the task of fitting the CdF for each considered item to distinguish more the different response variable. F (α j (θ -δ j )). Figure 1a presents the different CdF for one item j given α j = 1 and δ j = 0. Figure 1b presents the logistic CdF adjustment for three items (j = 1, 2, 3) with different α j and δ j = 0, given the linear predictor defined in equation (1).

In the literature, the cumulative model is presented according to the use of several of the previously mentioned CdF [START_REF] Samejima | Estimation of Latent Ability Using a Response Pattern of Graded Scores1[END_REF][START_REF] Fahrmeir | Multivariate Statistical Modelling Based on Generalized Linear Models[END_REF][START_REF] Liu | A mixed-effects regression model for longitudinal multivariate ordinal data[END_REF], while the adjacent models are most often presented with the logistic CdF. However, for both, we recommend the use of a symmetric CdF. As seen in Figure 2b, the IRT parametrization is a subtle and effective way to take into account the multiple item outcome characteristics in GLMM for categorical responses. In the context of the HRQoL in clinical trial, the HRQoL dimension considers a small set of items which are correlated and measure an unique latent variable. The discrimination parameter is routinely not use in this kind of analysis.

Relatively to the literature, Table 2 shows the specification of the famous polytomous IRT models following the different components. For IRT models being within the class of GLMM, we proposed to define them with the four components (r,F,Z q ,U r ). The kind of considered location item parameters can be indicated by the index q where q = 1 for including only difficulty parameters. Let q = 2 for considering the rating scale model [START_REF] Andrich | A rating formulation for ordered response categories[END_REF] parameterization where the difficulty parameters are common for all items and one shift parameter is considered for each item. Regarding the random part, the number of random effects can be then indicate by the index r. For the classical IRT parametrization presented in Table 2, only one random effect (r = 1) is taken into account : the capacity parameter θ. For IRT models including discrimination parameters for each item, we proposed to replace the components Z and U by a component specifying that the predictor is no longer linear (nl), such as (r,F,nl).

Table 2: Specification of the famous IRT model following the components : (r,F,Z p ,U r ) for the GLMM and (r,F,nl) for IRT model with no longer linear predictor. Index p denotes the number of kind of item parameters considered in the IRT model and r the number of random effect.

IRT models η

(j) im (r,F,Z p ,U r ) Rating scale model θ i -(δ m + τ j ) (adjacent,logistic,Z 2 ,U 1 ) Partial credit model θ i -δ jm (adjacent,logistic,Z 1 ,U 1 ) Sequential Rasch model θ i -δ jm (sequential,logistic,Z 1 ,U 1 ) Graded response model α j (θ i -δ jm ) (cumulative,logistic,nl) Generalized partial credit model α j (θ i -δ jm ) (adjacent,logistic,nl)
3 Simulation study

In the previous section, we focused on the use of the mixed models for ordinal data analysis and their relevance in the HRQoL analysis in oncology was discussed. Some comparisons studies exist between these different approaches [START_REF] Blanchin | Comparison of CTT and Rasch-based approaches for the analysis of longitudinal Patient Reported Outcomes[END_REF][START_REF] Anota | Comparison of three longitudinal analysis models for the health-related quality of life in oncology: a simulation study[END_REF][START_REF] Barbieri | Applying the Longitudinal Model from Item Response Theory to Assess the Health-Related Quality of Life in the PRODIGE 4/ACCORD 11 Randomized Trial[END_REF], mainly on the fixed part of the mixed models to identify the trend of latent trait. [START_REF] Anota | Comparison of three longitudinal analysis models for the health-related quality of life in oncology: a simulation study[END_REF] had shown an equivalent capacity to detect a fixed effect for the LMM and for one of the IRT models. Indeed, even if the LMM take into account the HRQoL score, which is a summary variable, this approach is at least equivalent to the IRT models in terms of power.

In this simulation study, the adjacent and cumulative models with the same parameterization of the linear predictor and the logistic CdF were used (as usually in the IRT models). The aim of the following section is to reinforce these comparisons between the LMM and the IRT models on the random part of the mixed models. The datasets were simulated from an IRT model (adjacent and cumulative models). Regarding the parameterization, two subject-specific random effects ξ i0 and ξ i1 were considered, respectively associated with the intercept and the slope. Of course, the usefulness of the random effect introduction in the model is strongly associated with the observed data. HRQoL is a subjective endpoint, and the individual random effect ξ i0 is thus entirely justified. Indeed, it is easy to imagine that each patient has a different level of HRQoL at baseline. The random slope is more questionable, indeed, the assumption that the specific HRQoL evolution of one single patient diverges from the average evolution for the whole population, is less obvious than the previous one. In this section, the capacity of the mixed models to detect the slope random effect was thus studied. No group effect was considered in this simulation study.

Design

We want to study the capacity of each model to detect the random effect ξ i1 associated with time (random slope). The two subject-specific random effects are considered independent where ξ i0 ∼ N (0, σ 2 0 ) and ξ i1 ∼ N (0, σ 2 1 ). The following model choice study is performed on the basis of the Bayesian information criteria (BIC) where two models were considered: M 2 with the two random effects (r,F,Z 1 , U 2 ) and M 1 excluding the random slope (r,F,Z 1 , U 1 ). For the IRT models, the linear decomposition of the latent trait θ iv only took into account the time as a fixed effect. The two considered models with proportional design are:

M 2 : θ iv = (t v -t 0 ) β 1 + ξ i0 + (t v -t 0 ) ξ i1 M 1 : θ iv = (t v -t 0 ) β 1 + ξ i0 (7)
In order to best reflect the EORTC QLQ-C30 questionnaire, the most frequent HRQoL dimension with two items (j = 1, 2) comprising four response categories (m ∈ {0, . . . , M } with M = 3), was used to design the simulation study. A sample size of three hundred subjects (i = 1, . . . , n with n = 300) and eight follow-up time (v = 0, . . . , 7), as for the trial presented in the previous section, were considered. The datasets were simulated from a multinomial distribution. The different response probabilities π

(j) ivm = Pr Y (j)
iv = m|θ iv , δ j concerning the subject i for item j were determined by equation ( 6) for the adjacent model and by equation ( 5) for the cumulative model, given:

• the item parameters (δ j1 , δ j2 , δ j3 ) j=1,2 ;

• the latent trait (θ iv ) deduced in accordance with equation ( 7);

• the logistic CdF, F η (j) ivm = exp η (j) ivm 1 + exp η (j) ivm , where η (j) ivm = θ iv -δ jm .
The values of the parameters used were deduced from the pain symptom data of the clinical trial presented in the previous section. We considered two kinds of difficulty parameters: near δ ne = (δ ne 1 , δ ne 2 ) and far δ f a = (δ f a 1 , δ f a 2 ). These parameter values were chosen in order to illustrate several scenarios described in Table 3. The different scenarios were due with the different associations between the model used to simulate the data, (adjacent,logistic,Z 1 , U r ) r=1,2 or (cumulative,logistic,Z 1 , U r ) r=1,2 , and the different considered values of the difficulty parameters. Table 3 shows the simulated responses expected at baseline (t = 0). The responses simulated across time depended of the considered coefficient β 1 . Each scenario was simulated N = 500 times.

Table 3: Values of difficulty parameters used to simulate the data and expected responses at t 0 under each studied scenarios.

Difficulty parameters Models

δ ne 1 = (-1.6, 1, 1.45)

δ f a 1 = (-2.1, 1, 2.75) (r, F, Z 1 , U r ) r=1,2
δ ne 2 = (-0.8, 1.15, 1.9)

δ f a 2 = (-1.25, 1.4, 3.3) (adjacent,logistic,Z 1 , U r ) r=1,2
balanced responses focus on center categories (1 and 2) (cumulative,logistic,Z 1 , U r ) r=1,2 focus on extreme categories (0,1 and 3) balanced responses

Concerning the LMM, the scoring procedure proposed by the EORTC was considered [START_REF] Fayers | EORTC QLQ-C30 Scoring Manual[END_REF], and the score associated with a symptomatic dimension was first calculated using the simulated data. Let the two simulated ordinal outcomes y

(1) iv and y

(2) iv concerning the individual i at the visit v, the related score was:

S iv = J=2 j=1 y (j) iv 2 100 M
Similarly to the parameterization in equation ( 7), we took into account the related choice model with:

M 2 : S iv = β l 0 + (t v -t 0 ) β l 1 + ξ i0 + (t v -t 0 ) ξ i1 + ε iv M 1 : S iv = β l 0 + (t v -t 0 ) β l 1 + ξ i0 + ε iv
where β l 0 is the fixed parameter associated with the intercept, the ξ l are the random effects normally distributed with the mean equals to zero and ε iv ∼ N (0, σ 2 ε ) the error term.

Results

Table 4 shows the capacity of the three models (adjacent model, cumulative model and LMM) to detect the random slope given different scenarios (Table 1). When we simulated the data under M 2 according to the random effect variances estimated from real data, each model detected the random slope (ξ i1 ) in 100% of cases whatever the different situations. On the contrary, under M 1 , the simulated model M 1 was correctly chosen in most cases. For all simulations under M 1 , the cumulative model seemed to detect the random slope although it was not included in the simulation step. Moreover, the IRT model which was not use to simulate the data, wrongly detected this random effect given a negative value of β 1 and the difficulty parameter coefficients δ ne . This could be explained by the fact that the difficulty parameters were not uniformly separated around zero and also because they were too close. Indeed, given β 1 < 0, the probabilities to observe the upper categories were very small over time and under-represented in comparison with the lower categories (as illustrated in Figure 2a). In the specific case where β 1 = -0.3, the IRT model which did not simulate the data could not explain the different outcomes only with the fixed effect and the random intercept, and it compensated the lack of information with the random slope.

We then could expect symmetric results from β 1 (positive values) considering the opposite sign of the difficulty parameters because of the reversibility property of the IRT models.

On the contrary, the LMM was stable and thus allowed making the good choice of model whatever the β 1 values and the IRT model used to simulate the data. Concerning the IRT models where only one model out of the two detected the random effect ξ i1 , the most suitable model seemed the one not detecting this random effect.

The capacity of the different models to detect the random slope when its variance value changes, is presented in Table 5. All models were sensitive to the signal-to-noise ratio. Indeed, the more β 1 increased, the less the random effect provided information. This was well characterized as the capacity to detect the random effect for greater variances when the signal was strong. In this case, the signal provided the essential information explaining the different responses. In the model comparison, the LMM was less sensitive than the IRT models. Indeed, the LMM detected the random slope for a greater variance of this one whatever the β 1 value. This result was expected because the LMM is based on the HRQoL score which is a summary variable with less information than the raw data. Thus, the IRT models are more sensitive in all cases. Comparing the two IRT models, there is a tendency for the random slope model to be preferred under the cumulative model regardless of whether it is the true model model or not. On the contrary, in the specific case where β 1 = -0.3, the IRT model used to simulate the data was less efficient than the other IRT model which detected a random slope to remedy the lack of information. This was coherent with our previous results shown in Table 4. Finally, the more β 1 was close to zero, the more the models detected the random slope for a low variance.

Table 4: Frequency (on N = 500 datasets) of the M 1 selection according to the BIC, given t v = (0, 1, 2, 4, 6, 8, 10, 12) and σ 2 0 = 1.5. For r = 1, 2, the (adjacent,logistic,Z 1 , U r ) models and the (cumulative,logistic,Z 1 , U r ) models are denoted respectively by AM and CM. For the random component,

U 1 if σ 2 1 = 0 and U 2 if σ 2 1 > 0.
Simulated Scenarios Model AM using δ ne CM using δ f a CM using δ ne AM using Table 5:

δ f a σ 2 1 β 1 LMM AM CM LMM AM CM LMM AM CM LMM AM CM 0.2 -0.3 0 0 0 0 0 0 0 0 0 0 0 0 0.2 0.3 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.
Frequency (on N = 500 datasets) of the M 2 selection according to the BIC, given t v = (0, 1, 2, 4, 6, 8, 10, 12) and σ 2 0 = 1.5.

For the (adjacent,logistic,Z 1 , U r ) models and the (cumulative,logistic,Z 1 , U r ) models are denoted respectively by AM and CM. For the random component,

U 1 if σ 2 1 = 0 and U 2 if σ 2 1 > 0.
Simulated Scenarios Model AM using δ ne CM using δ f a CM using δ ne AM using δ

Application to HRQoL data

In this section, we performed a longitudinal analysis of the HRQoL data from a multicenter randomized phase III clinical trial in first-line metastatic pancreatic cancer patients: PRODIGE4/ACCORD11 [START_REF] Conroy | FOLFIRINOX versus gemcitabine for metastatic pancreatic cancer[END_REF]. Three hundred and forty-two patients were randomly assigned to Folfirinox (experimental arm) versus Gemcitabine (control arm) regimens. The detailed inclusion and exclusion criteria, the study design and protocol, the treatment, the compliance to the questionnaires, and the HRQoL analyses have previously been published [START_REF] Conroy | FOLFIRINOX versus gemcitabine for metastatic pancreatic cancer[END_REF][START_REF] Gourgou-Bourgade | Impact of FOLFIRINOX Compared With Gemcitabine on Quality of Life in Patients With Metastatic Pancreatic Cancer: Results From the PRODIGE 4/ACCORD 11 Randomized Trial[END_REF][START_REF] Barbieri | Applying the Longitudinal Model from Item Response Theory to Assess the Health-Related Quality of Life in the PRODIGE 4/ACCORD 11 Randomized Trial[END_REF]. The patients filled the EORTC QLQ-C30 questionnaire themselves at different follow-up times defined in the protocol: at baseline, day 15, day 30, and at months 2, 4, 6, 8, and 10. The different measuring times reflected the longitudinal aspect of the HRQoL and allowed the assessment of the change of HRQoL for each dimension. Analyses were performed using the SAS software (version 9.3) [START_REF] Institute | SAS/STAT 9.3 User's Guide: Mixed Modeling[END_REF][START_REF] Boeck | Explanatory Item Response Models: A Generalized Linear and Nonlinear Approach[END_REF]. For all previous arguments, the cumulative models are preferred for the longitudinal analysis of the HRQoL. There, we considered the (cumulative,logistic,Z 1 , U 2 ) model with to analyze the data. In HRQoL study in oncology, the analysis is carried out for each HRQoL dimension. Given one HRQoL dimension with few correlated items, the discrimination parameters could be considered equals to one for each item. The distinction between the multiple-item responses is only achieved though difficulty parameters (thresholds) [START_REF] Anota | Comparison of three longitudinal analysis models for the health-related quality of life in oncology: a simulation study[END_REF][START_REF] Barbieri | Applying the Longitudinal Model from Item Response Theory to Assess the Health-Related Quality of Life in the PRODIGE 4/ACCORD 11 Randomized Trial[END_REF]. For the HRQoL longitudinal analysis with the subject i (i = 1, . . . , 342), the visit v (v = 1, . . . , 8), the item j with M j response categories, the (cumulative,logistic,Z 1 , U 2 ) model is defined by:

Pr Y (j) iv ≥ m|θ i = exp η (j) ivm 1 + exp η (j) ivm , (8)
with the following linear predictor is considered in the analyses:

     η (j) ivm =θ iv -δ jm θ iv =g i β 1 + (t v -t 0 ) β 2 + g i (t v -t 0 ) β 3 + ξ i0 + (t v -t 0 ) ξ i1 (9)
where:

• δ jm is the difficulty parameter (threshold) associated with the category m of item j;

• t v is the date of the visit v, and t 0 is the date of baseline;

• g i = 1 if the patient i belongs the experimental group (Folfirinox), g i = 0 if the patient i belongs the control group (Gemcitabine);

• β 1 is the effect difference at baseline between Folfirinox and Gemcitabine groups;

• β 2 is the slope (evolution) of health-related quality-of-life perception for the Gemcitabine group;

• β 2 +β 3 is the slope (evolution) of health-related quality-of-life perception for the Folfirinox group;

• ξ i0 and ξ i1 are respectively the subject-specific random effects associated with the intercept and the slope such as (ξ i0 , ξ i1 ) ∼ N (0, Σ), Σ being the unstructured covariance matrix.

These HRQoL data have been already analyzed with different approaches. Gourgou-Bourgade et al. ( 2012) have presented the results using time-to-event models. They concluded for a better HRQoL in the Folfirinox arm than the Gemcitabine arm. Then, [START_REF] Barbieri | Applying the Longitudinal Model from Item Response Theory to Assess the Health-Related Quality of Life in the PRODIGE 4/ACCORD 11 Randomized Trial[END_REF] have presented the results through the LMM and the partial credit model extended for the longitudinal analysis (adjacent,logistic,Z 1 , U 2 ). The conclusions of both mixed models are similar.

For the (cumulative,logistic,Z 1 , U 2 ) model, Table 6 shows the estimations of fixed parameters, their standard deviation and the associated P-value of the Wald test. Concerning the functional dimension, we performed a reverse permutation on the functional scales for an intuitive interpretation. This allows to consider that an increase of the latent trait θ is associated with an increase of the functional capacity (improvement of HRQoL) or increase of the symptoms (deterioration of HRQoL). For all HRQoL dimensions, there should be no difference at baseline (β 1 = 0) in a randomized clinical trial. However, we observed a significantly difference concerning the diarrhea symptom between the two groups at baseline (p = 0.007 * * ). It referred to a difference between the two arms at day 15, day 30 (during treatment period) but no necessarily at baseline. Then, the perception of diarrhea symptom remained higher in the Folfirinox arm over time. This result is expected because the Folfirinox is more toxic than Gemcitabine and is known to cause more diarrhea symptom.

Regarding the others dimensions, the HRQoL changed over time for several dimensions (emotional functioning, pain, insomnia, constipation and appetite loss) returning a significant improvement in terms of HRQoL perception. Only the pain showed a significantly different evolution between the two arms (p = 0.04). Indeed, the patients receiving the Folfirinox had a perception of pain which decreased significantly more over time than the patients receiving the Gemcitabine.

Its interpretability and intuitive illustration is one of the many advantages of the cumulative models. The constraints on the item parameter in these models allow an interpretation through the latent variable (e.g. comparing the proportions of the response categories for one specific item over time or between different groups during a fixed time). Figure 2 presents the evolution concerning the probability of response either over time (Figure 2a) or between group (Figure 2b). This example is illustrated through the first item of the pain symptom of the clinical trial previously presented. The probability (π m ) for a patient to response the category m corresponds to the area under the curve delimited by the horizontal lines. For both groups, Figure 2a shows that the probability to choose the categories 2 or 3 decreased over time while the probability to choose the category 0 increased. At baseline, the response proportion for the categories 0, 1, 2 and 3 were respectively π 0 = 0.10, π 1 = 0.62, π 2 = 0.22 and π 3 = 0.06 for each group. Then, the evolution of the proportions showed a decrease of the level of pain between the baseline and the 4-month visit, and, finally, a decrease of the latent trait over time. Likewise, Figure 2b shows the different response proportions between the two groups at four months. In the control group, the proportions were π 0 = 0.29, π 1 = 0.61, π 2 = 0.08 and π 3 = 0.02 for the categories 0, 1, 2 and 3. In the experimental group, they were π 0 = 0.47, π 1 = 0.48, π 2 = 0.04 and π 3 = 0.01. The probability to response category 3 was the lowest whatever the group, but was even less likely for patients in experimental group than in control group. In contrary, the probability to response category 0 was more likely in experimental group than in control group. The observed gab corresponded to the difference between the two linear predictors associated with each group only four months after the baseline. One of the interests of this illustration concerns the clinical interpretation. The IRT models thus offer a complete analysis: the general analysis of a HRQoL dimension and the specific analysis for each item [START_REF] Edelen | Applying item response theory (IRT) modeling to questionnaire development, evaluation, and refinement[END_REF]. (9) between the two groups four month after the Baseline (t 4 ).

Discussion

We have explored the suitable mixed models for the longitudinal analysis of the HRQoL in oncology. This data coming from questionnaires through Likert scales, we focused on regression models for ordinal data. These models have been specified with three components, the linear predictor parameterization, the ratio of probabilities and the CdF [START_REF] Peyhardi | A new specification of generalized linear models for categorical responses[END_REF]. In oncology, the analysis being performed on multiple-item measurements associated with one HRQoL dimension [START_REF] Fayers | EORTC QLQ-C30 Scoring Manual[END_REF], the specific IRT parameterization of the linear predictor is thus used. The item parameters allow to distinguish the outcomes from different items which measure an unique unidimensional latent variable. This latent variable was decomposed linearly to take into account the different covariates in the fixed part of the model and to incorporate subject-specific random effects. The analysis with IRT models is the richer because they are based on raw data [START_REF] Gorter | Why item response theory should be used for longitudinal questionnaire data analysis in medical research[END_REF]. The analysis can be made on one specific item through the item parameters or on the whole HRQoL dimension [START_REF] Edelen | Applying item response theory (IRT) modeling to questionnaire development, evaluation, and refinement[END_REF]. Indeed, these models take into consideration all available information from the data, it is why the use of this kind of model is more and more studied [START_REF] Gorter | Why item response theory should be used for longitudinal questionnaire data analysis in medical research[END_REF]. Then, concerning the choice of the model family, the cumulative and adjacent models are preferred. From the ratios of probabilities which characterize them and the symmetric CdF, the practical property of the invariant seems important to interpret the results. The cumulative models also assume an underlying continuous latent variable that is associated with a linear mixed regression model [START_REF] Mccullagh | Regression models for ordinal data (with discussion)[END_REF][START_REF] Hedeker | A random-effects ordinal regression model for multilevel analysis[END_REF]. This allows a better interpretation and illustration of the results such as the easy analysis of the evolution of the response proportions of the different categories over time or between groups, given one item. The adjacent models show the advantage not to have any constraint for the model estimation. These models can thus be preferred when the regression is performed on the item part of the linear predictor, given non-proportional design. Finally, the choice of the CdF essentially depends on the observed data and properties which interest the users. These IRT models are reversible only if the CdF is symmetric. Thu, the use a commonly symmetric CdF is preferred (the logistic and the Gaussian distributions). From the conceptual IRT model selection, the cumulative model seems the most suitable given its advantages for the longitudinal analysis of HRQoL in cancer clinical trials.

The simulation study showed that the IRT model capacity to detect the random effect was better than the LMM currently used. This result seems natural because the LMM is based on the study of a summary variable with less information. Thus, the variability from data is also reduced. Of course, the usefulness of the random effect introduction in the model is strongly associated with the observed data. Moreover, the more the difficulty parameters were distinct and the influence of covariates was stronger, the less the random effect provided information. All these results confirmed that the IRT models allow a more detailed analysis to interpret the results from a specific dimension or item. Whatever the IRT model used to generate the data, the LMM remained competitive through these simulations. However, the IRT model that did not generate data, seemed more sensitive to the random slope than the other IRT model used to simulate the dataset. Indeed, in some cases, it tended to detect the random slope while it did not exist. In case where one of the two models detects the random slope, the use of the model not detecting the effect as it is, seems the most appropriated concerning a data-driven choice. However, we recommend to use only one kind of model (with same components discussed previously) allowing to make the results comparable across HRQoL dimensions.

An aspect that remains to be discussed is the multidimensional aspect of HRQoL. Nowadays in oncology, the different dimensions are analyzed independently of one another, and this causes the problem of multiple tests. An approach to consider the all HRQoL dimensions would be the use of structural equation modeling. This would allow to show the influence of each HRQoL dimension through some factors to explain the global HRQoL and potential structural links between the latent variables.

Figure 1 :

 1 Figure 1: Relationship between the CdF and the IRT parameterization for dichotomous items where F (η) =

Figure 2 :

 2 Figure 2: Interpretation of the (cumulative,logistic,Z 1 , U 2 ) model through its underlying latent variable concerning the pain symptom (including the items 9 and 19). The estimated difficulty parameter for the item 9 are δ 9,1 = -2.1, δ 9,2 = 1 and δ 9,3 = 2.75, and for the item 19 : δ 19,1 = -1.28, δ 19,2 = 1.40 and δ 19,3 = 3.34.Figure 2a: the different HRQoL evolution of the latent variable and response variable (Y (j) ) j=9,19 between the two groups.Figure 2b: the different proportions (π m ) of different response categories of Y(9) between the two groups four month after the Baseline (t 4 ).

  Figure 2a: the different HRQoL evolution of the latent variable and response variable (Y (j) ) j=9,19 between the two groups.Figure 2b: the different proportions (π m ) of different response categories of Y

Table 1 :

 1 Summary of the characteristics for the three model families

	Models

1 : for some non proportional design models

Table 6 :

 6 Estimations of fixed effect parameters (β p ) p=1,2,3 of the (cumulative,logistic,Z 1 , U 2 ) model. All HRQoL dimensions of the EORTC QLQ-C30 are considered.

	HRQoL
	Dimensions	Coefficient Standard error	Pvalue
	Global Health Status
	β 2	0.098	0.070	0.166
	β 3	0.130	0.085	0.128
	Physical functioning
	β 2	-0.150	0.077	0.051
	β 3	0.122	0.098	0.212
	Role functioning
	β 2	-0.011	0.081	0.892
	β 3	0.157	0.103	0.131
	Emotional functioning
	β 2	0.335	0.070	< .001 * * *
	β 3	0.001	0.086	0.992
	Cognitive functioning
	β 2	-0.002	0.054	0.972
	β 3	0.088	0.067	0.189
	Social functioning
	β 2	0.010	0.073	0.888
	β 3	0.116	0.093	0.211
	Fatigue
	β 2	-0.087	0.085	0.308
	β 3	-0.033	0.107	0.761
	Nausea and vomiting
	β 2	-0.052	0.060	0.393
	β 3	-0.069	0.072	0.336
	Pain
	β 2	-0.330	0.076	< .001 * * *
	β 3	-0.188	0.092	0.040 *
	Dyspnea
	β 2	-0.060	0.075	0.420
	β 3	-0.093	0.088	0.295
	Insomnia
	β 2	-0.359	0.080	< .001 * * *
	β 3	0.046	0.083	0.627
	Appetite loss
	β 2	-0.354	0.072	< .001 * * *
	β 3	-0.026	0.080	0.747
	Constipation
	β 2	-0.325	0.077	< .001 * * *
	β 3	0.003	0.083	0.974
	Diarrhea
	β 1	0.739	0.272	0.007 * *
	β 2	0.018	0.067	0.792
	β 3	-0.026	0.076	0.786
	Financial difficulties
	β 2	-0.522	0.282	0.066
	β 3	0.302	0.208	0.146

Acknowledgements

We thank Dr. Hélène de Forges for her editorial assistance and UNICANCER for the data from PRODIGE4 / ACCORD11 clinical trial which is used in this paper.

Funding

This study was supported by a grant from the French Public Health Research Institute (www.iresp.net) under the 2012 call for projects as part of the 2009-2013 Cancer Plan.


Pierre De Micheli

Amico Settefrati

Stéphane Marie

Julien Barlier

Patrice Lasne

Benjamin Scholtes

Marc Bernacki

François Bay

Dr.-Ing P De Micheli

Dr -Ing A Settefrati

Dr -Ing S Marie

Dipl.-Ing J Barlier

Dipl.-Ing P Lasne

Prof Dipl.-Ing B M Scholtes

Prof F Bernacki

Bay

Towards the simulation of the whole manufacturing chain processes with FORGE®

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

Introduction

It is quite obvious that final properties of forged parts depend strongly on the main forming processes, but they also depend significantly on all the other processes of the manufacturing chain. As a consequence, the core of a metal forming Finite Element (FE) software is and remains the thermo mechanical computation of the forming processes, with a lot of efforts and innovations delivered to enhance robustness, accuracy and efficiency in this domain (classical and complex rheological behaviours modelling, anisotropic mesh adaptation, bi-mesh methods, multi material simulation, smoothed tools, steady state simulations) [START_REF] Revil-Baudard | A detailed description of the Gurson-Tvergaard-Needleman model within a mixed velocity-pressure finite element formulation[END_REF][2][3][START_REF] Perchat | Multi criteria anisotropic adaptive remeshing for sheet metal forming in FORGE® In[END_REF][START_REF] De Micheli | Dramatic speed-up in FEM simulations of various incremental forming processes thanks to multi-mesh implementation in FORGE® In[END_REF][START_REF] Chenot | Finite element simulation of multi material metal forming In[END_REF][START_REF] Hachani | A smoothing procedure based on quasi-C1 interpolation for 3D contact mechanics with applications to metal forming In[END_REF][START_REF] Ripert | Efficient numerical simulation of steady-state metal forming processes by an iterative surface calculation In[END_REF]. However, in order to tackle their daily challenges and better communicate with suppliers and customers, the blacksmith industry will also need a more global software package able to simulate the whole manufacturing chain. Transvalor falls within this approach with the development of FORGE® and our other softwares in order to be able to provide this global metal forming simulation tool. The knowledge in mechanical and numerical engineering in this domain has been improved a lot in the last decades, but does not permit to fully simulate the manufacturing chain yet. Experimental mechanical engineers are still working hard to enhance our comprehension of the physical phenomena and to describe them in models, as well as numerical mechanical engineers, working to deliver accurate, stable and efficient numerical platforms to run those new models on a wide range of industrial processes simulations. A brief bibliographic review on the state of the art concerning numerical simulation of the whole manufacturing chain will be presented in the first part of this article. It will underline both physical challenges, where experimental studies are needed, and numerical challenges, where numerical improvements are needed. Then, as a software editor, we will focus the scope of this article on some points where physical phenomena are at least partially known and already described in generic-enough models, but where the numerical methods are not completely satisfying and could be improved. It will naturally lead us to the following parts of this article, focusing on innovative methods proposed by Transvalor to better address those numerical challenges. Microstructure prediction, induction hardening simulation and nitriding simulation will be detailed and discussed.

State of the art

For economic reasons, the metallurgical industry has always looked for reducing the number of steps in the manufacturing chain and in the same time still improving the final parts properties. As the mechanical behaviour is closely linked to the material structure at the micro and nanoscale, a quantitatively correct prediction of the structural state from casting until the final thermal and thermochemical treatments becomes crucial in order to meet both requirements.

Billet history

The first step of a classic metal forming simulation consists in providing the initial billet properties. Most of the time, the billet history is neglected, and the initial state of the metal is described empirically, using in general constant parameters in the whole part. However, when defects may exist, when initial properties are not uniform, measuring the properties distribution is not straightforward. Simulating the billet history from the mill to predict those parameters would be more convenient than having to perform a complex bunch of experiments.

The numerical answer to this problem consists in connecting an existing casting simulation tool with a metal forming simulation tool. Main metal forming FE software packages have made collaboration with other software editors [9,10]. TRANSVALOR, unlike its competitors, has interconnected two in-house software packages, THERCAST® and FORGE®, to ensure a maximal flexibility and compatibility. Available physical models allow predicting yet grain structure evolution during solidification, composition of phases in equilibrium, cracks, porosities, macrosegregations distributions and for most of them their evolutions during the whole manufacturing chain [11,12] (see figure 1). No specific part will be dedicated to this issue in this article.

Microstructural evolution during forging

In response to the strong demand from the metallurgical industry to reduce characterization campaigns and to introduce more realistic physical behaviour in their simulations, numerous models are proposed today to link chemical composition, microstructural evolution and mechanical behaviour.

Mean field models are generally associated with low computation times making them suitable for full industrial computations. In these approaches, the microstructure is macroscopically described selecting representative material parameters (grain size, inclusions, phase percentage, precipitates, etc.) considered in an average way. The identification of the physical laws governing the evolution of these parameters allows predicting the microstructural evolution and their influence on the mechanical behaviour.

Concerning the modelling of recrystallization phenomena, mean field models are becoming years after years more physic, more generic and more accurate [13,14]. Even if, to our knowledge, no model is generic enough to be easily adapted to a reasonable range of alloys and processes, software editors are also working hard in order to integrate them in their codes. The main point for the software editor is not to deliver some default models, but to make them easy to modify, and to allow an advanced user to implement new ones.

Heat treatments

The control of the structural evolutions occurring during heat treatments of metallic alloys is one of the best ways to confer at the material particular properties adapted to its future industrial application.

Physical models used to predict metallurgical structure evolutions in steels specifically on cooling and quenching are very similar for all main metal forming FE software packages.

These models are generally based on the description of the transformation kinetics of austenite during isothermal treatments (represented on TTT diagrams). The anisothermal kinetics can then be calculated using the Johnson-Mehl-Avrami-Kolmogorov (JMAK)

formalism [21][22][23][24][START_REF] Kolmogorov | Statistical theory of crystallization of metals In[END_REF] and the additivity principle [26]. Consideration of the properties (i) of each phase (rheology, thermophysical parameters) and (ii) of each phase transformation (latent heat, crystal structure change) allows predicting parts distorsions. More basic semi-empirical models also exist for steel austenitization and annealing simulations.

If the thermomechanical computation and the phase transformations simulation are fully coupled one to the other, thermal computation in case of heat treatments remains challenging on the numerical point of view.

First, very local phenomena occur at the part boundaries and an accurate simulation is possible only if the mesh size is adapted to the problem. The use of isotropic or anisotropic mesh adaptation is then critical. Moreover, modelling of heat supply for furnace or blowtorch heating and bath or spray cooling is generally quite well approximated with well-defined boundary conditions. However, in the case of induction heating and hardening, the electromagnetic computation of the eddy current in the part is necessary to deduce the heating rate distribution. The problem is driven by the well known Maxwell equations, and the challenge is now on the numerical side. The common approach, based on connection elements between the inductor and the part, is well suited for 2D computations. However, 3D computations are almost impossible because of the big band width, which means bad convergence, high computation times and bad parallel scalability [27,28]. The other approach, chosen by TRANSVALOR, consists in discretizing the whole domain in one single mesh (parts, inductors, air) [29-35]. This method is perfectly adapted to 3D computations and is easily parallelizable, with a much smaller band width. The work to extand it to parallel computation and to processes with complex relative motions of the part and the inductor will be presented in section 5 of this article.

Thermochemical treatments

The difference between a pure heat treatment and a thermochemical treatment is that in the second one, at least one physical element is diffused toward the metal.

In the case of carburizing, only carbon is diffused in the matrix, which can be simulated by most of the main metal forming simulation software packages. Diffusion properties and boundary conditions are quite well known and not really problematic. As for heat treatment, mesh adaptation is here the key numerical feature to obtain both efficiency and accuracy in the computations. The subsequent quenching can then be simulated taking into account the shift of the phase transformations kinetics due to the modification of the local carbon rate in the material.

Nitriding is a more complex process. First, on the numerical side, two elements are diffused in the metal (Nitrogen and Carbon), which cannot be simulated by all software packages. Secondly, on the physical side, several precipitates can appear and dissolve

Models description

Semi-empirical JMAK approach

The analytical modelling of recrystallization phenomena very often uses modified versions of the JMAK equation. For isothermal and steady state conditions, the different recrystallized volume fractions X and recrystallized grain sizes d evolutions were established as a function of initial grain size d 0 and processing variables such as equivalent strain rate ̅ , temperature T and equivalent strain ̅ . An overview of well known relations is given in Table 1.

Dynamic recrystallization Static recrystallization

( ( ̅ ) ) ̇ ( ) ̇ ( ) ̇ ( ) ( ( ) ) ̅ ̇ ( ) ̅ ̇ ( )
Table 1: Example of equations describing (i) the kinetics of dynamic and static recrystallization and (ii) the evolution of the recrystallized grain sizes.

With an extremely low CPU time, this approach is well-suited for computations at each integration point of a finite element mesh. Implemented in the FORGE® user routines library, the model has undergone several modifications in order to be used with any thermo-mechanical local loading and to take into account multiple recrystallization steps that can be encountered during multi-pass processes [START_REF] Teodorescu | Development of a level set methodology to simulate grain growth in the presence of real secondary phase particles and stored energy-application to a nickel-base superalloy[END_REF]. Comparison between simulation and experimental results for Inconel718 in figure 2 illustrates the good agreement obtained as soon as the parameter values are correctly calibrated.

Physically based approach, mean field model

The JMAK based model requires a significant experimental effort and the parameters are difficult to interpret physically. In order to have predictive capabilities in complex industrial conditions, physically based models have emerged to link the metallurgical principles.

First of all, the microstructure is described by a set of internal averaged variables representative of the material (average dislocation density and equivalent grain size for the recrystallized and non recrystallized grains). The microstructural evolution is then directly given by the modelling of the physical laws governing the evolution of these parameters.

During a deformation step for example, the evolution of the dislocation density depends on two competiting processes namely work-hardening and dynamic recovery softening and is often described by the Kocks-Mecking model:

√ eq. 1
Where the term √ represents the dislocation density storage and is the softening parameter. Dynamic recrystallization modelling necessitates (i) the modelling of the nucleation rate of the recrystallized grains and (ii) the modelling of the grains growth kinetics. The mean dislocation density can then be calculated by:

( ) eq. 2
Where is the dislocation density of the deformed matrix, is the dislocation density of the recrystallized grains and is the recrystallized volume fraction.

Since the metal behaviour depends on the dislocation density, such approach allows a direct coupling with the rheological behavior of the material. As an example, the Taylor equation hereafter defines the flow stress from the dislocation density by:

√ eq. 3
With  the shear modulus, b the Burgers vector, M the Taylor factor, a "dislocation free" yield stress and α a constant [13,39].

The user routines structure proposed in FORGE® allows implementing or modifying easily such kind of models with the possibility to couple the material rheology with the metallurgical evolution (see figure 3a). These mean field approaches are generally associated with low CPU times, making them suitable for macroscale computations over a whole part (see figure 3b).

Full field approach

The idea behind these approaches is that morphology and topology play a non-negligible role in the evolution of the microstructure, which is difficult to capture with other models.

Microstructures are explicitly represented and simulations are performed on a

Representative Volume Element (RVE) at the mesoscopic scale. The polycristal is constructed respecting the topological characteristics of the grains and the metallurgical properties [40,41]. This approach is based on a level set (LS) description of interfaces in a finite-element context: for each individual grain, a signed distance function , gives at any point the distance to the grain boundary. Thus the interface is given by the level 0 of the function . This method has the advantage of avoiding the difficult problem of tracking interfaces.

The model works in 2D and in 3D and computational cost is reduced thanks to an anisotropic mesh adaptation around the interfaces (see figure 4). As the interface moves, periodic remeshing is performed such that the refinement zone always coincides with the interface position. This approach has been introduced to model primary recrystallization, including the nucleation stage and has been extended to take into account the grain growth stage [19,20,41]. Moreover in this LS context, Zener pinning phenomenon can be taken into account in a natural way [42]. Thus this aspect will be considered in the development of the DIGI-µ® software. In primary recrystallization, the kinetics of interface motion is directly linked to the state variables stored in the mesh, related to the stored strain energy.

In the case of pure grain growth, the curvature of the grain boundaries pilots the interface motion (see figure 5 and6).

Conclusion

Three approaches are proposed to predict metallurgical evolutions occurring during thermo-mechanical processes. The two first approaches are based on mean field models.

Well-suited for large-scale simulations, they have been implemented into FORGE® software package through the user routines library. The third approach simulates the microstructural evolution at the mesoscopic scale where the microstructure components are fully modeled. Modelling at this scale can be a help for understanding complex microstructural phenomena but it can also be used to optimize/calibrate higher scale models (like mean field models). These simulations cannot be directly integrated in a macroscale finite element computation, but they are much closer to physical phenomena and allow describing in a natural way the materials in terms of microstructural features.

Induction hardening 4.1 Principle

Induction heating and more precisely induction hardening are nowadays widely used.

Those processes allow a quick and local heating, but predicting the temperature distribution without the use of numerical simulation in real 3D configurations is quite complex. FORGE® proposes a unique formulation for metal forming FE software, enabling the full 3D computation of the electromagnetic field in the whole domain (air, workpiece, surrounding parts, inductors), coupled to the standard FORGE® computation.

Finite element formulation

The problem consists in solving Maxwell equations system in order to deduce the heat power created by the circulations of the eddy currents in the part. Several simplifications can however be applied to Maxwell equations system in the case of induction heating.

First, due to the dimensions of industrial induction heating installations and the current frequencies used in the inductor, from about 50 Hz to 1MHz, the magneto-quasi-static hypothesis can be done.

The problem is then written in a single non-linear equation thanks to the A-V formulation

[29] described in the first line of eq. 4. ⃑ represents the magnetic vector potential, the electric scalar potential, the electrical resistivity and the magnetic permeability. In order to have a unique solution, we have to add the two divergence free gauges corresponding to the second and third line of eq. 4.

{ ⃑ ⃑ ⃑⃑ ( ⃑ ⃑⃑ ⃑ ) ⃑ ⃑⃑ ⃑ ⃑⃑ ⃑ ( ⃑ ⃑⃑ )
eq. 4

The magnetic field ⃑⃑ and the electric field ⃑⃑ are linked to ⃑ and by:

{ ⃑⃑ ⃑ ⃑⃑ ⃑ ⃑⃑ ⃑ ⃑ ⃑⃑ eq. 5
This equation is discretized in space on the whole domain (inductors, air, parts, and concentrators) using Nedelec edge elements [45]. This equation is solved over a few electric current periods until the periodic state is reached (in general three periods are reasonable). Then, the mean local heating rate on the last electric period is computed and used as a heat source term in the standard thermo mechanical computation of FORGE®.

The thermo mechanical time step used in FORGE® can of course be several orders of magnitude bigger than the current period in the inductor. It is however reasonable to compute several thermo mechanical time steps with the same heat source term. A new electromagnetic computation is done only when the computation configuration has undergone significant enough changes (temperature changes, making and evolve, or relative motions of the parts and the inductors).

This electromagnetic formulation, in comparison with the boundary element appraoach, presents the advantages to generate symmetric numerical problems with small bandwidth, which ensure a good performance for sequential and parallel resolutions (Single Program Multiple Data paradigm) with iterative solvers [START_REF] Coupez | From a Direct Solver To a Parallel Iterative Solver in 3-d Forming Simulation[END_REF]47].

Results

This formulation enables the computation of full 3D complex geometries for the workpieces and the inductors. We will see in the following example local induction hardening on a crankshaft's crankpin's journal with an adapted inductor (see figure 7), with both motions of the crankshaft and of the inductor (see figure 8). After this local heating, the part is quenched in polymer and water. During the process, the crankshaft turns around its main axis at 60 rpm. The inductor also turns around the main axis of the crankshaft at 60rpm, but also around the axis of the crankpin's journal at -60 rpm.

The austenization rate prediction during heating on a cutting plane in the crankpin's journal can be seen in figure 8.

Figure 8: Austenite rate on a cutting plane in the crankpin's journal during the first second of the induction heating process.

A precise distribution of phases during the heating and the following quenching process can be computed. Final hardness and residual stresses can be predicted by the solver (see figure 9). We can see in this example that the martensite rate is not homogeneous on the crankpin's journal, which is indeed what is observed in the real process. The use of magnetic field concentrators should be studied in order to improve the process.

Conclusion

The technical choices made by FORGE® enable to perform with high accuracy a very wide range of induction heating simulations, with complex 3D geometries and even complex motions between the inductor and the part. Those computations are fully parallelizable, with a good parallel efficiency. Thanks to the parallel version, complex induction heating calculations with relative motions between the part and the inductor can be done in a reasonable time on parallel cluster machine. The features presented here will be available in the next commercial release of FORGE®.

Nitriding

Principle

Nitriding can be applied on various metals under various conditions. It consists most commonly in diffusing nitrogen into the surface of a steel part at about 550°C to harden the surface and improve its mechanical properties. This process has a lot of advantages in terms of fatigue life, wear, friction, corrosion and a relatively good predictiveness of the induced distortions. However, it remains very long and then expensive, with about 100 hours for treating less than 1 mm of depth.

Directly at the surface, a thin white layer from 0 to 30 microns composed of nitride ε and nitride γ' is formed. Then nitrogen will diffuse from this layer in a solid solution of ferrite to form what is called the diffusion layer, which is about 0.1mm to 1mm thick (see figure 10). Some of this freely diffusing nitrogen will precipitates with various elements of the alloy, included carbon. This is the reason why carbon diffusion can play an important role, even without carbon enrichment of the part: the gradient created in the freely diffusing carbon because of nitrogen and carbon precipitates will also engender a diffusion of the dissolved carbon, and then can modify the composition at equilibrium of the alloy. The presence of the nitrogen precipitations is the main cause of hardening and of the creation of compressive stresses.

Finite element formulation

Our task during the DEFINIT project consisted in adapting the dedicated 1D nitriding simulation tool, developed at Arts et Metiers ParisTech in Aix-en-Provence, to full industrial 3D simulations, and confront it to the industrial partners' results [36,37]. In our finite element simulation, the problem is described by adding several variables at each element: the rate of free nitrogen, free carbon (free means here dissolved, free to diffuse), and the rate of phases and precipitates. First, the evolution of free carbon and nitrogen rates is computed at each node by solving two diffusion problems, described by the following equation. C represents here the chemical element concentration free to diffuse, and D the diffusivity.

( ⃑⃑⃑⃑⃑⃑⃑⃑⃑⃑⃑ ( )) eq. 6

Then, the composition of the metal in the element is computed depending on temperature (which is in general constant), carbon rates and nitrogen rates. Precipitates can be formed or dissolved, also modifying the free carbon and nitrogen rate. All data needed to compute the metal evolution can be obtained with the software JMatPro or Matcalc. The composition computed at each element enables local volume variation computation and hardness estimations.

Of course, such a computation only makes sense if the mesh is fine enough at the boundaries to model the phenomenon, which means that the mesh size must at all cost be smaller than 0.2 mm for a 1mm diffusion layer, or even smaller.

Results

We have the privilege to present here some results we obtained on a fully nitrided gear, provided with courtesy of Airbus Helicopters, one of the partners of the DEFINIT project (see figure 11).

Figure 11: Geometry of the gear provided by Airbus Helicopters (the blue zone is not affected by the nitriding process), and the corresponding anisotropic mesh used to the finite element simulation (a fourth of a tooth).

The use of anisotropic mesh is here of first importance in order to describe the phenomenon occurring in the diffusion layer, but on a fourth of a tooth. We can obtain a mesh size equal to 0.05mm at the boundaries with only 22000 nodes in the mesh. In comparison, we would have needed 1.5 millions of nodes with an adapted isotropic mesh to obtain the same precision (see figure 12).

Conclusion

With those simulations, FORGE® has opened the doors of full 3D simulation of nitriding on industrial parts. First results are encouraging. Thanks to anisotropic meshing, accurate predictions can be obtained in the diffusion layer of a full industrial part. Parameters data files can be provided by JMatPro® or ThermoCalc® for first approximation simulations, allowing the user to easily predict phase and precipitates distribution, and then deducing hardness, residual stresses and distortion thanks to the full thermo mechanical simulation.

We will carry on this ongoing work to provide you as soon as possible an advanced industrialized nitriding simulation tool in the commercial version of FORGE®.

Conclusion

In order to match the increasing quality requirements for the mechanical parts, a better understanding of the effect of the manufacturing processes on metal composition and microstructure evolution is needed. In the last decades, a lot of experimental work has been carried out to study the different physical phenomena at work, and build models.

Those models have been implemented in FE software packages and have democratized microstructure computations.

The main issue from the experimental point of view is not only to describe the phenomena occurring in the different processes but also to obtain the material parameters for the models. This long term work requires high technical and scientifical knowledge as well as advanced laboratory devices. Our historic and on-going links with Cemef Mines ParisTech and our collaborations with renowned customers and academic research laboratories allow us enriching our software packages with the most relevant models and corresponding associated data.

From the numerical point of view, this article has shown that the main issues for microstructure and surface treatments simulations were not so different than those encountered in usual metal forming simulations. The key is to provide a performing simulation tool able to run those models on complex industrial workpieces, to catch local phenomena with always more accuracy, and everything in a reasonable computation time.

In order to describe with accuracy the microstructure evolution in a component, with mean field models, during induction hardening or nitriding, or even more to follow grain topology on a representative volume element, the first requirement is to discretize the domain finely enough. Isotropic and anisotropic mesh adaptation is the first answer in order to reduce as much as possible the number of degrees of freedom. However, when it is not anymore possible to reduce it, parallel and massive parallel computations must be used. For that reason, each TRANSVALOR's software package is based on a full parallel approach, which allows obtaining very high parallel efficiency and very high speedups.

To complete that, an open data structure based on user routine concept allows the user to adapt those models to its particular needs and then to improve the predictions provided by the simulation.

In order to make sure that any user will obtain with FORGE® and any other TRANSVALOR's software packages the best accuracy cost ratio for his advanced physical model applied to industrial configurations, those lines of work will remain in the future among Transvalor's highest priorities.

Figure 1 :

 1 Figure 1: (a) Following porosities from casting to final shape with THERCAST® and FORGE®. (b) Grain structure predicted using CAFE modelling in THERCAST® will be soon importable in FORGE® [12].

Figure 2 :

 2 Figure 2: Comparison of simulated and measured average ASTM grain sizes on an Inconel718 part after a rotary forging process (courtesy Tecnalia, ES).

  Figure 3: (a) Simplified FORGE®'s organisation chart. (b) Microstructural evolution calculated with a mean field model during a reducer-rolling simulation

Figure 4 :

 4 Figure 4: Polycrystal description in 3D with finite element mesh (anisotropic meshing near grain boundaries)

Figure 5 :Figure 6 :

 56 Figure 5: Microstructure evolution of a 304L stainless steel undergoing the heat treatment cycle defined on the left (pure grain growth evolution). Initial configuration with a mean grain diameter of 90µm

Figure 7 :

 7 Figure 7: Geometry of the part of the crankshaft and of the inductor used for local induction hardening.

Figure 9 :

 9 Figure 9: Hardness (HV) and equivalent strain distribution in the crankshaft's journal after quenching.

Figure 10 :

 10 Figure 10: Micrograph showing the white layer and the diffusion layer after nitriding (Arts et Metiers ParisTech Aix-en-Provence)

Figure 12 :

 12 Figure 12: Comparison between an anisotropic (22000 nodes) and an isotropic mesh (1.5 millions of nodes) providing the same accuracy for nitriding simulations.

  It is possible in FORGE® thanks to the user routines library [15].

	If this mean field framework is quite convenient, it can be synonymous for a given
	material of a large amount of experiments with advanced laboratory devices. Therefore,
	the idea to replace these heavy and expensive experiments by full field simulations on a
	Representative Volume Element (RVE) is born. Over the last decades, several mesoscale
	numerical methods have been developed to simulate explicitly the microstructure
	evolution due to recrystallization during hot forming processes [16-18]. The Center for
	Material Forming (Cemef Mines ParisTech) has recently proposed a solution based on a
	FE formulation and a level-set description of phase and grain interfaces to model
	recrystallization and grain growth phenomena [19, 20] which will soon end up on the first
	industrial solution DIGI-µ® developed by TRANSVALOR.

  at equilibrium depending on the local chemical composition and temperature for different alloys. Gathering all this available information in user routines library, nitriding simulations have been performed with FORGE®, with results in pretty good agreement with experiments and 1D dedicated academic simulation software, as we will see in section 5 of this article.Metal forming modelling can be predictive only if the strain rate, strain and temperature dependency of the flow behavior is correctly described. The mechanical properties and behavior of metallic materials depending mostly on the dislocation content and structure,

	3 Microstructure evolution modelling during hot forging
	3.1 Principle

themselves, and the boundary conditions through the white layer are not yet well understood. Experimental studies start to obtain encouraging results on some alloys, and have traduced them into 1D dedicated simulation tools

[36, 37]

. Thermodynamical calculations provided for example by JMatPro® or ThermoCalc® can predict phases and precipitates composition this points out the need to incorporate microstructure concepts into those models. The goal is to correctly describe the main physical mechanisms occurring in metals during thermomechanical processes i.e. work-hardening, recovery, grain boundary migration, nucleation and grain growth related to dynamic, static or metadynamic recrystallization and phase transformation.