Denis Pallez

Christel Dartigues

Parisa Ghodous

Michel Martinez

Data Architecture for Collaborative Conceptual Design

Keywords: -Computer Aided Design, Shape synthesis, Multiple View Points, Product Data Integration, Concurrent engineering

Currently, computer aided systems have concentrated on the capture and representation of geometric shape and technical information, as opposed to providing support for product design. The aim of this paper is to propose a framework for an integrated software that assists designers in the early design phases to work in co-operative and collaborative manner.

I. INTRODUCTION

The technological advances carried out these last years in the field of products development led the researchers to elaborate the approaches that reduce the cost and time of product development, enhance the quality of product and help the designers to be more creative. These objectives are difficult to obtain due to the great number of phases, which should be carried out during the product development (Fig. 1) and the important number of experts of different disciplines that are involved. Concurrent Engineering is considered one of the key concepts that enables the companies to reach these objectives [START_REF] Ghodous | Editors of Advances in Concurrent Engineering[END_REF].

Concurrent engineering is a wide field of research. The researchers who are interested in concurrent engineering, work on different aspects such as [START_REF] Prasad | Concurrent Engineering Fundamentals[END_REF]:

• Philosophical aspect deals with the boundaries of the responsibility and the authority, culture and organization management, • Methodological aspect deals with system thinking, approaches to system complexity, systems engineering, product realization taxonomy and system integration, • Conceptual aspect deals with concurrency and simultaneity, modes of concurrency and cooperation, work flow mapping and • Virtual aspect deals with capturing life cycle intent and information modeling. Our work is situated on the virtual aspect. In the product design, each expert works with his own applications, handles his own data, has knowledge and constraints which are specific to his field of work and has his own point of view on the product [START_REF] Ghodous | Collaborative and standard design and manufacturing model[END_REF]. This heterogeneity implies many problems, in particular the exchange and the sharing of information with the other experts. It is thus necessary to make it possible to each expert to represent adequately its data while facilitating integration and communication of their data with the other experts.

In this context, our objective is to develop a concurrent and collaborative system allowing experts to participate on product development as soon as possible and help them to work together [START_REF] Dartigues | Modélisation collaborative des activités de conception[END_REF]. In the following, firstly, we present the function to form mapping approach as a conceptual design technique. Then we explain why it is more realistic to consider this methodology in a collaborative context and we present enhanced version of the function to form mapping approach presented in [START_REF] Gardan | On shape to specifications adequacy[END_REF][START_REF] Gardan | From Functions to Shapes[END_REF]. Finally, we present the data models that are necessary to achieve a collaborative and conceptual design aided system based on this methodology. Due to the reasons of normalization, these models are described by EXPRESS-G formalism [START_REF]Part 11, Description Methods: The EXPRESS Language Reference Manual[END_REF].

II. THE FUNCTION TO FORM MAPPING AS A TECHNIQUE FOR CONCEPTUAL DESIGN

The aim of any design is to obtain as soon as possible the product so as to be most profitable as possible. In a nonroutine design, it is delicate and extremely complex to obtain the best products answering customer's specifications. Indeed, a great number of different experts of different disciplines participate in the product design. The difficulty lies in the fact that those persons have to collaborate. Shape is today the main representation of a product, even though the current trend is to remove geometry from its central position in order to add high level information. In fact, geometry constitutes the starting point of many activities such as mechanical optimizations, kinematics simulations, and so on. The function to form mapping appears to be one of the most important activities of the design process and up till now happened manually. This activity is very important both for the choices that are made and for the quantity of work that it represents [START_REF] Gardan | A methodology for a function-to-shape translation tool in foundry[END_REF].

That is why we have the desire to assist designers in the first stages of product design. The objective is not to construct the shape automatically but to both automate a certain number of heavy and tiresome tasks, and assist the designer during the first stages of design. In best case, this assistance makes the designers' stimulation possible by presenting them solutions that they had not thought before.

We suppose that the first information we have at our disposal are included in the customer's specifications which is expressed in natural language (Fig. 2a). Due to the functional decomposition (Fig. 2b) using functions of making, maintaining, prevention, control [START_REF] Keuneke | Device Representation. The significance of functional knowledge[END_REF] and allowing [START_REF] Mcdowell | Conceptual design for polymer composite assemblies[END_REF], behaviors of the product can be identified (Fig. 2c) in order to respect the philosophy of FBS [START_REF] Tomiyama | A CAD for functional Design[END_REF][START_REF] Umeda | Supporting conceptual design based on the function-behaviour-state modeller[END_REF][START_REF] Ranta | Integration of Functional and Featurebased product modelling --the IMS/GNOSIS experience[END_REF]. Next, with the help of the same philosophy, structure of the product can be obtained (Fig. 2d). This model can be identified with product's components linked by assembly relation. In some cases, we have noticed that a function can also be decomposed into a set of constraints on physical parameters, called intermediate specifications (Fig. 2e). Parameters used in the intermediate specifications are called intermediate parameters [START_REF] Gardan | On shape to specifications adequacy[END_REF][START_REF] Gardan | From Functions to Shapes[END_REF] and are of a rather high level. Intermediate parameters are defined as quantifiable and measurable entities referring to the physical world. Just as the parameters and the specifications, constraints established in the intermediate specifications are also called intermediate constraints. In this way, it will be possible to specify the internal characteristics of the component or eventually its links with other components. By instance, consider the function "to be handled easily by a human-being" for a water bottle. This function can be decomposed into some constraints such as: • Weight less than 2 kilo, • Volume not greater than 3 litres, • Material used must be light and so on… In [START_REF] Gardan | From Functions to Shapes[END_REF], we have defined an intermediate constraint by a quadruple <IP, R, Exp, W> where IP is an intermediate parameter, R is a relation among {<, >, =, ≠}, Exp is an arithmetic expression, W is the relative weight of constraint in comparison with the other constraints of the intermediate specifications.

Once a great number of intermediate constraints are obtained, we consider that it is sufficient to propose shapes solutions to the designer that achieved those intermediate constraints, and so that satisfied functions, from which the constraints are deduced. We suppose that this technique could be applied in a more general way. From the product's specifications (Fig. 2a) to the intermediate constraints (Fig. 2d), actions achieved by the designer during the design process are currently manual. They might be assisted with difficulties in so far as manipulated information are mainly expressed in natural language. However, as intermediate constraints can be formalized, we can attempt to automate the intermediate constraints to form mapping in order to assist the designer. For this stage, we propose to generate the component shape of the product. Before generating the product shape that is composed of components, the complexity is reduced by generating a component after the other and next by assembling shapes of each component of the product. A component will be produced from its own intermediate parameters. Each component is a rigid, finite and homogenous solid. A set of shapes called "solutions space" (Fig. 2f) is generated for a component from primitive shapes. Naturally, several shapes could be proposed to the user: we must detect the shapes, which are best satisfying the intermediate constraints (Fig. 2g). For that, we have defined a Satisfaction Degree SD ic (s) of an intermediate constraint ic by a given shape s as a real number in the [0,1] interval. This number expresses the quality with which s satisfies the intermediate constraint ic. If SD ic (s) is near zero, the constraint is badly satisfied, if it is near one, the constraint is well satisfied. The method to compute SD ic (s) changes depending on whether s is a primitive solid or a complex object [START_REF] Gardan | On shape to specifications adequacy[END_REF]. Fig. 2 Overall architecture of the system In this section, we will first detail the way to generate solutions space. Then, we will present the advantages and drawbacks of such method, and finally we will investigate how to apply it in collaborative environment.

A. Generation of the solutions space

The solutions space is composed of various topologies that are solutions (like spheres, boxes or more complex shapes...) that we call shapes classes. To model geometry associated with each topology, we associate each shapes class with a set of physical parameters of lower level. They are called terminal parameters for distinguishing them from physical parameters and are strongly related to the geometry (the radius of a sphere, the width, the length and the height of a box...). They must represent the shapes class geometry in a coherent way. For instance, a tetrahedron can be characterized indifferently by the lengths of each side of the base and the height on the one hand, or by three angles and the height on the other hand. Choosing a real value for each terminal parameter determine a solution shape. Finally, the solutions space is generated as automatically as possible by using expert's rules alike algorithms that translate intermediate constraints into terminal constraints. As a single intermediate constraint can generate several terminal constraints, we have to compute the intersection between the all generated terminal constraints into intervals for each terminal parameter.

B. Synthesis of this methodology

First of all, one has to consider this methodology as a means of automatically generating solutions shapes that may give ideas to the designer for the final product shape. We applied such a methodology within the framework of the filling system for foundry mould design [START_REF] Gardan | A methodology for a function-to-shape translation tool in foundry[END_REF]. We validated the methodology by taking the specificity of the foundry (trade features) into account. We showed that this approach was applicable to this case and that we obtained quickly results better than those based on an expert's knowledge. Moreover, the design process of this methodology consists in presenting several solutions shapes that best satisfy the intermediate constraints to the designer. Next, the design process does not stop there in so far as the designer has to modify the intermediate specifications if the resulting shapes are not good enough for him and recursively until to obtain the most promising shapes. In fact, the previous method could be hard to use by only one person because the intermediate constraints' consistency is in charge of the designer, even if there is several thousands of constraints. In addition, the same designer ought to manipulate all the intermediate constraints in order to produce an acceptable solution shape: to do this, he ought to know the signification of all the constraints; that is rather unrealistic in most cases. A more suitable methodology is to distinguish between intermediate constraints from a specific domain of activity and intermediate constraints from other domain of activity. As a consequence, we treat in the next section the way to apply this approach in a collaborative context with the aim of reducing productivity time but also of managing a lot of experts points of view.

C. Methodology applied in a collaborative context

As it seems to be unrealistic to design the product shape only by a designer due to the large number of design domain that intervenes, we propose to put the previous methodology in a collaborative context. First of all, the different stages that consist in translating product's specifications (Fig. 2a) into intermediate specifications (Fig. 2e) are unchanged except the fact that there is at least one expert per design domain. Consequently, we can group the intermediate constraints based on the same design domain together (Fig. 3a): for instance, we have distinguished mechanical, electrical and thermal design. The adaptation of the previous methodology in the collaborative context consists in applying the method for each design domain. Consequently, each expert of one domain will generate a solutions space (Fig. 3b) from his own intermediate specifications and his own primitive shape library. Recursively, each expert can modify its solutions space by modifying its intermediate specifications as the previous method did. Once each expert is satisfied about his solutions space, one has to merge all the solutions space from the different design domain in order to find the most promising solutions that satisfied all the domains (Fig. 3c). The difficulty lies in the fact that each solutions space is composed by shape classes, and each of them is composed by terminal parameters that are totally different in theory because of their origin of different domain. As a consequence, it is unrealistic to merge all the solutions space from all the design domain because few shape class combinations could be done. Currently, techniques that come in our mind are stochastic or morphing methods as studied in [START_REF] Gardan | Towards a specifications-to-shape translation tool[END_REF]. Finally, once a technique will be found, one can visualize the most promising solutions using computation of the Satisfaction Degree.

III. PROPOSED MODEL FOR REPRESENTATION OF EXPERT VIEWPOINT

The analysis of current works on product modeling shows that current single fixed representations are inadequate to model the various concepts present in multidisciplinary product development situation. Consequently, the dynamic representation of multiple views of a product based on functional contexts seems to be necessary. Depending on the view taken, certain properties and descriptions of the object become relevant. A comprehensive model of a product must be able to built depending on the particular need.

Consider the example of the mechanical part in Fig. 3d. The mechanical design model of this product is different of electrical engineer's model or thermal model. Any model should allow a dynamic evolution and must be capable of accommodating multiple concepts unambiguously and consistently so that the elements could not be duplicated. Any inconsistencies between the various models have to be discovered and corrected. This process may go through several iterations. The result is a set of models, one per consulting discipline, where, although each set represents the product using a different point of view, the comprehensive representation is consistent. There is no attempt to integrate the various sets into one set.

The basic description of a product differs from one viewer to another. Each view may represent a product with different elements and different composition hierarchies.

No one model contains a full comprehensive description of the product but each model should be consistent regarding to the object being described. Different descriptions of the same elements and different subsets of these descriptions in different models exist. Some models used for product description are shown in Fig. 4. We detail the model of intermediate specification

View point model represents the information about a product from a particular view and the relationship between the different viewpoints. Several functions are related to each viewpoint.

Purpose explains why an object does what it does and it is related to the human socio-cultural environment concept, The purpose model represents the purpose and the relationship between the purposes.

Function is what an object does. Functional model represents function and the different relationships between the functions.

Behavior is how the object does what is does. Behavior model represents the behavior, the relationships between the behaviors and the relationships with functions.

Structure is what the object is. We use the STEP standard to represent the concepts related to product data (the structure of product and its relations with other models such as shape, materials, tolerances, etc). The definition of a product in the STEP product data model is any physical object, which is produced by either natural or manufacturing processes. Any part or assembly that contributes to a product is also considered to be a product. A car is a product while its wheels and engine assemblies are considered as other products. Furthermore, each of these products can be further decomposed into smaller components or products. The details of these models are described in [START_REF] Ghodous | Collaborative and standard design and manufacturing model[END_REF]. We present the detail of this model in the following.

A. The intermediate specification model

We represent the intermediate specification concepts by EXPRESS-G formalism (Fig. 5). EXPRESS-G is a graphical language developed by ISO 10303 STEP [START_REF]Part 11, Description Methods: The EXPRESS Language Reference Manual[END_REF]. The EXPRESS-G basic notations used in Figures include entities (rectangles); super-type/subtype relationships (thick solid lines); required attributes (normal lines); relationship for optional attributes (dashed lines). Additionally, the direction of an attribute is symbolized by an open circle, where the circle represents the "many" side of a "one to many" relationship.

In Fig. Lastly, due to the fact that each expert handles his own whole of constraints, there are relations between them (IC-Relationship). These relations are of two types. We define the Same_as relation, which expresses the fact that two experts handle the same constraint. Their expression and their total weight will be then identical, but their local weight might be different. We also define the NCU relation (Next Constraint Usage), which means that an intermediate constraint can be compose of one or more other constraints.

IV. CONCLUSION AND FUTURE WORKS

The analysis of current works on product modeling shows that current single fixed representations are inadequate to model the various concepts present in multidisciplinary product development situation. Consequently, the dynamic representation of multiple views of a product based on functional contexts seems to be necessary. In this paper, we have studied the problem of representation of experts' multiple-view in a collaborative conceptual design environment. We have considered the function to form mapping approach as a conceptual design method and we have presented the possibilities to automatically check whether a shape satisfies the specifications. Our opinion is that whether the shape is simple (a primitive solid), or complex (combined object) this operation can not be performed without participation and collaboration of different human experts. However, it seems that the software will be more and more self-sufficient as its experience will grow, particularly for simple shapes.

Our future works concerns how to maintain the consistency of models and how to integrate the different shapes solution spaces. We have defined some coherence rules between different models. Each expert at any time can define his model and may collaborate with the other models. As a consequence, when one model is manipulated, corresponding effects will be made automatically in the other. Future works will focus on the definition and formalization of these rules to improve proposed multiple-view model.

V. REFERENCES

Fig. 3

 3 Fig. 3 New methodology in a collaborative environment

Fig. 4

 4 Fig. 4 Some necessary models and their relationship

Fig. 5 :

 5 Fig. 5: Intermediate specifications model expressed in a EXPRESS-format

  5, the intermediate specifications are the set of different intermediate constraints . The intermediate constraints are elaborated from the functions that characterize the product. As the functions relate to a given view point (mechanical view point, electrical view point, thermal view point and so on), the intermediate constraints also relate to this same view point. On the other hand, each intermediate constraint is defined on a component, i.e.: is related to the structure of the product (entities Product-Definition and Product-Definition-Relationship).

	The intermediate constraints are expressed in the following
	way: first of all, an intermediate constraint described by an
	expression. This expression contains an intermediate
	parameter (which can be the volume, the weight or any
	engineering parameter), a relation (<, >, =, ≠, etc.) and a
	mathematical expression, which can possibly include other
	intermediate parameters. With each intermediate constraint
	are associated two weights: a local weight and a total
	weight. These weights are used to give a list of priority of
	the constraints that must be carried out for the expert of a
	precise field (local weight) and for all of the experts (total
	weight).


Esther D Ellen
email: esther.ellen@wur.nl

Vincent Ducrocq

Bart J Ducro

Roel F Veerkamp

Piter Bijma

Combining survival analysis and a linear animal model to

de documents scientifiques 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.

Background

Mortality due to cannibalism in laying hens is a worldwide economic, health, and welfare problem, occurring in all types of commercial poultry housing systems [START_REF] Blokhuis | Studies of feather pecking in poultry[END_REF]. Due to the likely prohibition of beak-trimming in the European Union in the near future, this problem will increase if no further actions are taken, and, therefore, needs to be solved urgently.

One of the possibilities is to use genetic selection [START_REF] Muir | Group selection for adaptation to multiple-hen cages: Selection program and direct responses[END_REF][START_REF] Jones | Genetic selection for poultry behaviour: Big bad wolf or friend in need?[END_REF]. However, selection for lower mortality has not been very effective in most cases [START_REF] Preisinger | Internationalisation of breeding programmes -breeding egg-type chickens for a global market[END_REF]. First, heritabilities of mortality are low, ranging between 3.2% and 9.9%, leading to low accuracy [START_REF] Mielenz | Mortality of laying hens housed in single and group cages[END_REF][START_REF] Craig | Fearful and associated responses of caged White Leghorn hens: Genetic parameters estimates[END_REF][START_REF] Robertson | The heritability of all-or-none traits: Viability of poultry[END_REF][START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF][START_REF] Ducrocq | Genetic improvement of laying hens viability using survival analysis[END_REF]. Second, censoring is high (animals still alive at the end of the testing period have no record on survival time) [START_REF] Ducrocq | Genetic improvement of laying hens viability using survival analysis[END_REF], leading to low accuracy as well. Third, traditional methods for selection against mortality can lead to unfavourable response to selection, because these methods ignore the social effect an individual has on it's group members (so-called social interactions) [START_REF] Muir | Group selection for adaptation to multiple-hen cages: Selection program and direct responses[END_REF][START_REF] Griffing | Selection in reference to biological groups I. Individual and group selection applied to populations of unordered groups[END_REF][START_REF] Ellen | Genetic improvement of traits affected by interactions among individuals: Sib selection schemes[END_REF][START_REF] Muir | Group selection for adaptation to multiple-hen cages: Selection program and responses[END_REF].

Heritabilities for survival traits are often estimated using a linear animal model [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF][START_REF] Bijma | Multilevel Selection 2: Estimating the genetic parameters determining inheritance and response to selection[END_REF]. However, a linear animal model does not take into account the fact that some animals are still alive at the end of the testing period (so-called censored records), for these animals the true survival days is unknown. Furthermore, linear models do not properly account for the nature of survival data, because survival data are usually heavily skewed [START_REF] Kachman | Applications in survival analysis[END_REF]. Survival analysis [START_REF] Kalbfleisch | The statistical analysis of failure time data[END_REF] appropriately accounts for both censoring and non-normality in the data. Survival analysis is used to examine either the length of time an individual survives or the length of time until an event occurs. Models for survival analysis can be built from a hazard function, which measures the risk of an event to occur, given that the individual has survived up to time t [START_REF] Kachman | Applications in survival analysis[END_REF][START_REF] Kleinbaum | Survival analysis: A self-learning text[END_REF].

Social interactions occur when individuals are kept together in a group. Wolf [START_REF] Wolf | Genetic architecture and evolutionary constraint when the environment contains genes[END_REF] has mentioned that the environment provided by group members is often the most important component of the environment experienced by an individual in that group. There is clear evidence that social interactions contribute to the heritable variation in traits [START_REF] Muir | Group selection for adaptation to multiple-hen cages: Selection program and direct responses[END_REF][START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF][START_REF] Bijma | Multilevel Selection 2: Estimating the genetic parameters determining inheritance and response to selection[END_REF][START_REF] Wolf | Genetic architecture and evolutionary constraint when the environment contains genes[END_REF][START_REF] Wade | Group selection among laboratory populations of Tribolium[END_REF][START_REF] Wade | An experimental study of group selection[END_REF][START_REF] Muir | Incorporation of competitive effects in forest tree or animal breeding programs[END_REF][START_REF] Moore | The inheritance of social dominance, mating behaviour and attractiveness to mates in male Nauphoeta cinerea[END_REF]. For instance, social interactions have a substantial genetic effect on mortality due to cannibalism [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF][START_REF] Bijma | Multilevel Selection 2: Estimating the genetic parameters determining inheritance and response to selection[END_REF][START_REF] Wolf | Genetic architecture and evolutionary constraint when the environment contains genes[END_REF][START_REF] Muir | Incorporation of competitive effects in forest tree or animal breeding programs[END_REF][START_REF] Brichette | A genetic analysis of intraspecific competition for growth in mussel cultures[END_REF][START_REF] Arango | Estimation of variance components including competitive effects of Large White growing gilts[END_REF][START_REF] Van Vleck | Effect of competition on gain in feedlot bulls from Hereford selection lines[END_REF][START_REF] Bergsma | The contribution of social effects to heritable variation in finishing traits of domestic pigs (Sus scrofa)[END_REF][START_REF] Chen | Estimation of genetic parameters for average daily gain using models with competition effects[END_REF]. Bijma et al. [START_REF] Bijma | Multilevel Selection 2: Estimating the genetic parameters determining inheritance and response to selection[END_REF] and Ellen et al. [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF] have found that 1 / 3 to 2 / 3 of the heritable variation in survival days is due to social interactions. To reduce mortality due to cannibalism, the classical model for a given genotype must be extended to consider not only the individuals' direct effect of its own genes, but also the associative genetic effect of the individual on the phenotypes of its group members [START_REF] Griffing | Selection in reference to biological groups I. Individual and group selection applied to populations of unordered groups[END_REF]. Muir [START_REF] Muir | Group selection for adaptation to multiple-hen cages: Selection program and direct responses[END_REF] has clearly shown that selection methods targeting both direct and associative genetic effects (group selection) results in a decrease in mortality due to cannibalism in laying hens, whereas selection based on only the direct genetic effect (individual selection) results in an increase in mortality [START_REF] Craig | Group selection for adaptation to multiple-hen cages: Beak-related mortality, feathering, and body weight responses[END_REF]. Furthermore, Muir [START_REF] Muir | Incorporation of competitive effects in forest tree or animal breeding programs[END_REF] has found that, in Japanese quail, group selection results in decreased mortality and increased bodyweight. However, so far associative genetic effects have not been implemented in existing software for survival analysis. To analyse data on mortality due to cannibalism, a solution might be to combine survival analysis and a linear animal model including associative effects.

Ducrocq et al. [START_REF] Ducrocq | Implementation of an approximate multitrait BLUP evaluation to combine production traits and functional traits into a total merit index[END_REF] have proposed a two-step approach for multiple trait evaluation of longevity and production traits in dairy cattle, which faces similar problems. The two-step approach is a combination of survival analysis and a linear animal model. In the first step, survival analysis is performed to compute the so-called pseudorecords and their associated weights. Pseudo-records can be regarded as the result in the data of a linearization of the model. When analysed with a simple linear animal model, pseudo-records weighted appropriately lead to the same estimated genetic values as the initial survival model used to compute them. In the second step, genetic parameters on pseudo-records with their associated weights are estimated using a linear animal model.

In this paper, we apply a similar two-step approach to estimate genetic parameters for direct and associative effects on survival time in laying hens. In the second step, we will use the linear animal model including associative effects to estimate genetic parameters [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF][START_REF] Bijma | Multilevel Selection 2: Estimating the genetic parameters determining inheritance and response to selection[END_REF][START_REF] Muir | Incorporation of competitive effects in forest tree or animal breeding programs[END_REF]. For the remaining part of the paper, we will refer to the linear animal model including associative effects as LAM and to the two-step approach as 2STEP. Cross validation will be used to compare 2STEP with LAM [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF][START_REF] Bijma | Multilevel Selection 2: Estimating the genetic parameters determining inheritance and response to selection[END_REF]. LAM was applied directly to estimate genetic parameters for social effects on observed survival days. For the cross validation, the predicted hazard rate will be estimated using 2STEP and the predicted phenotype will be estimated using LAM. To judge the performance of both methods, predicted phenotypes or hazard rates will be compared with the observed phenotype.

Methods

For this study, the same data were used as described in Ellen et al. [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF]. The main characteristics are summarized below and further details are in [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF].

Population and housing

Data of three purebred White Leghorn layer lines from Institut de Sélection Animale B.V., a Hendrix Genetics company, were used in this study. The three lines were coded: W1, WB, and WF. For each line, observations on survival time of a single generation were used. Chickens of each line were hatched in two batches, each batch consisting of four age groups, differing by two weeks each. All chickens had intact beaks.

When the hens were on average 17 weeks old, they were transported to two laying houses with traditional four-bird-battery cages. Each batch was placed in another laying house. In both laying houses, the 17week-old hens were allocated to laying cages, with four birds of the same line and age in a cage. The individuals making up a cage were combined at random. In both laying houses, cages were grouped into eight double rows. Each row consisted of three levels (top, close to the light; middle; and bottom). A feeding trough was in front of the cages, and each pair of back-to-back cages shared two drinking nipples.

Pedigree

Sires used for both laying houses were largely the same while dams were different. For all three lines, sires and dams were mated at random. Each sire was mated to approximately eight dams, and each dam contributed on average 12.3 female offspring. Five generations of pedigree were included in the calculation of the relationship matrix (A). To avoid pedigree errors, hens with unknown identification or double identification were coded as having an unknown pedigree (n = 101). The observations on these hens were included in the analysis to better estimate fixed effects.

Data

All hens were observed daily. Dead hens were removed from the cages and not replaced, and wing band number and cage number were recorded. The study was ended when hens were on average 75 weeks old. For each hen, information was collected on survival and number of survival days. Survival was defined as alive or dead (0/1) at the end of the study. From these data, the survival rate was calculated as the percentage of laying hens still alive at the end of the study. Survival days were defined as the number of days from the start of the study (day of transport to laying houses) till either death or the end of the study. Hens that died before the end of the study were referred to as a failure (event = 1), whereas hens still alive at the end of the study were referred to as censored (event = 0). In total, 196 hens were removed from the study, due to reasons other than mortality. These hens were referred to as censored (event = 0). For the statistical analysis, 6,276 records were used for line W1; 6,916 for line WB; and 3,588 for line WF.

Data analysis

Data were analysed separately for each line. Two methods to estimate genetic parameters were compared: 1) LAM, a linear animal model including direct and associative effects applied directly to the observed survival days; this procedure is described in detail in [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF], and 2) 2STEP, a two-step approach [START_REF] Tarrés | Validation of an approximate approach to compute genetic correlations between longevity and linear traits[END_REF]. In the first step of 2STEP, data were analysed using survival analysis as implemented in the survival kit V5 [START_REF] Ducrocq | The survival kit a Fortran package for the analysis of survival data[END_REF], to produce pseudo-records as defined below. Survival analysis allows the combination of information from hens still alive at the end of the study (censored records) as well as hens that died (uncensored records). In the second step, genetic parameters for direct and associative effects on pseudo-records were estimated using a linear animal model [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF][START_REF] Bijma | Multilevel Selection 2: Estimating the genetic parameters determining inheritance and response to selection[END_REF], implemented in ASReml [START_REF] Gilmour | ASReml Users Guide Release 1[END_REF].

Step 1: Survival analysis Data were analysed using the Cox animal model [START_REF] Cox | Regression models and life tables[END_REF]. The Cox model can deal with non-linearity, censoring, and non-normal residuals. The model included a fixed effect for each combination of laying house, row, and level, and for average survival days in the back cage to account for a possible effect of the back neighbours [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF]. Age was fully confounded with laying house and row and, therefore, not included as a fixed effect. All the fixed effects were significant.

Using survival analysis results in a breeding value (a i ) and an associated weight (ω i ) for each hen i. It can be shown that ω i is the estimated cumulative risk of animal i from time 0 to censoring time or death, and is therefore a function of the (possibly censored) length of life of hen i, her censoring code (δ i = 0/1), and the fixed effects in the model [START_REF] Tarrés | Validation of an approximate approach to compute genetic correlations between longevity and linear traits[END_REF]. The pseudo-record for survival time of animal i was [START_REF] Ducrocq | A general approach for international genetic evaluations robust to inconsistencies of genetic trends in national evaluations[END_REF]:

y i i a i i * , = + -   1 (1)
where δ i is the censoring code of individual i (δ i = 1 if animal i is uncensored; δ i = 0 if animal i is censored); a i is the estimated direct breeding value of individual i; and ω i is the associated weight of individual i. Pseudorecords are functions of the data and of the effects estimated in the survival model, such that when a straightforward BLUP animal genetic evaluation is applied on these pseudo-records, the same estimated breeding values are obtained as in the initial survival model.

To verify 2STEP, pseudo-records with appropriate weights were analysed to estimate breeding values with a univariate BLUP animal model, with a heterogeneous residual variance 1 /  ∧ i for animal i. The correlation between the estimated breeding values of 2STEP and the estimated breeding values of the survival analysis was calculated [START_REF] Tarrés | Validation of an approximate approach to compute genetic correlations between longevity and linear traits[END_REF]. As expected, this correlation was one and the estimated breeding values were the same. Thus the computation of pseudo-records in 2STEP was correct.

Step 2: Associative effects model

To estimate variances and covariances for direct and associative effects, using the pseudo-records and associated weights from step 1, the model of Muir [START_REF] Muir | Incorporation of competitive effects in forest tree or animal breeding programs[END_REF] and Bijma et al. [START_REF] Bijma | Multilevel Selection 2: Estimating the genetic parameters determining inheritance and response to selection[END_REF] was used:

y Z a Z a e D D S S = + + , (2)
where y is a vector of the pseudo-records y i *; a D is a vector of direct breeding values, with incidence matrix Z D linking observations on individuals to their direct breeding value; a S is a vector of associative breeding values, with incidence matrix Z S linking observations on individuals to the associative breeding values of their group members (i.e., individuals in the same cage); and e is a vector of residuals, where Var e i e i ( ) = 1 2   . A weighted analysis was performed using the associated weight (ω i ) and the !WT statement in ASReml [START_REF] Gilmour | ASReml Users Guide Release 1[END_REF] and fixing  e 2 to one [START_REF] Ducrocq | Implementation of an approximate multitrait BLUP evaluation to combine production traits and functional traits into a total merit index[END_REF]. The covariance structure of genetic terms is

Var D s a a C A ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ⊗ , where C = ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥     A A A A D D S DS S 2 2 , in which  A D 2 is the direct genetic variance,  A S 2
is the associative genetic variance, and  A DS is the direct-associative genetic cov- ariance. Bijma et al. [START_REF] Bijma | Multilevel Selection 2: Estimating the genetic parameters determining inheritance and response to selection[END_REF] have shown that residuals of group members are correlated due to non-genetic associative effects. The covariance structure of the residual term, e, is given by Var( ) e R =  e 2 , where R ij = 1 when i = j, R ij = r when i and j are in the same group (i ≠ j), and R ij is zero otherwise. The value of r was estimated in the analysis, using a CORU statement in the residual variance structure in ASReml [START_REF] Gilmour | ASReml Users Guide Release 1[END_REF].

Heritable variation

When social interactions exist among individuals, each individual interacts with n -1 group members. In this study, n = 4. The total heritable impact of an individual on the population, referred to as its total breeding value (TBV), equals the sum of its direct breeding value and n -1 times its associative breeding value: [START_REF] Muir | Incorporation of competitive effects in forest tree or animal breeding programs[END_REF]. The total heritable variation equals the variance of the TBV among individuals,    

TBV i = A D , i + (n -1) A S,i
TBV A A D A DS S n n 2 2 2 2 2 1 1 = + - + - ( ) (
) [START_REF] Bijma | Multilevel Selection 2: Estimating the genetic parameters determining inheritance and response to selection[END_REF][START_REF] Bijma | Multilevel Selection 1: Quantitative genetics of inheritance and response to selection[END_REF]. With unrelated group members, the phenotypic variance equals    

P A A e D S n 2 2 2 2 1 = + - + ( )
. The total heritable variance expressed relative to the phenotypic variance equals

T TBV P 2 2 2 =  
. The T 2 expresses the total heritable variance relative to the phenotypic variance and is, therefore, a generalisation of the conventional h 2 to account for social interactions.

Cross validation

We compared 2STEP to LAM using cross validation [START_REF] Stone | Cross-validatory choice and assessment of statistical predictions[END_REF]. With cross validation, known phenotypes are set to missing and their value is predicted and compared with their observed phenotype. Validation was applied separately to each of the three lines. For this purpose, a random number was allocated to each cage within a fixed effect class. For each line, phenotypes of animals from 20% of the cages from each fixed effect class were set to missing, which resulted in five subsets, each containing 80% of the data. In this way, each cage was once removed from the total dataset, and each fixed effect class was present in all five subsets. The phenotypes set to missing were predicted using a combination of the direct breeding value of the individual itself and the associative breeding values of its group members

DSBV = + ( ) - A A D i n S j , , Σ 1
of either 2STEP or LAM. Comparing the predicted phenotypes of both methods is difficult for two reasons. First, a scale difference exists between estimated breeding values (EBV) of 2STEP and EBV of LAM. EBV of LAM are on the observed scale for survival days, whereas EBV of 2STEP are on the hazard rate scale. Transforming EBV of 2STEP into survival days is somewhat difficult, because the transformation is non-linear. Therefore, the predicted phenotypes using 2STEP are on the hazard rate scale, whereas the predicted phenotypes using LAM are on the observed scale for survival days. Second, in our dataset approximately 50-70% of the data were censored (animals that were still alive at the end of the testing period). These animals do not have an observed phenotype. In other words, a large proportion of the "observed" phenotypes is censored, and cannot be compared directly to their prediction. However, we know that their observed phenotypes are larger than those of animals that are not censored, which is highly relevant information.

To deal with these two difficulties, we used two approaches to evaluate both methods. The first approach is based on using groups of animals rather than single individuals. In this approach, for each subset and method, 25% of the animals with the best predicted phenotypes or hazard rates were selected as the best groups (best refers to animals with the highest predicted phenotypes using LAM or lowest predicted hazard rates using 2STEP), and 25% of the animals with the worst predicted phenotypes or hazard rates were selected as the worst groups. The Kaplan-Meier estimate of the survival curve was plotted for the best and worst groups based on the observed phenotypes. It was expected that the best groups would yield the best Kaplan-Meier estimate of the survival curve, whereas the worst groups would yield the worst one. Moreover, for both methods the mean observed survival days were calculated for the best and worst groups. From these, the difference in survival days between the best and worst group was calculated. For the best groups the percentage of overlapping animals, between 2STEP and LAM, was calculated.

To quantify the contribution of social effects to the predicted phenotype, phenotypes or hazard rates were predicted using different EBV: 1) classical BV (CBV); 2) direct BV of the individual itself (DBV = A D, i ); associative BV of the group members

( SBV = - Σ n S j A 1 ,
) and a combination of the direct BV of the individual itself and the associative BV of its group members

DSBV = + ( ) - A A D i n S j , , Σ 1
CBV were estimated using a classical linear animal model given in [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF] or using survival analysis (first step of 2STEP). DBV, SBV and DSBV were estimated using LAM or 2STEP.

The second approach is based on using ranks of individuals rather than their observed phenotypes. In this approach, the rank correlation between the observed phenotype and predicted phenotype or hazard rate and between predicted phenotype and predicted hazard rate was calculated. Due to the scale difference and censoring, Pearson correlations cannot be used. However, in both methods animals with the highest predicted phenotype (LAM) or lowest predicted hazard rate (2STEP) have the highest expected value for observed survival days. Therefore, the rank correlation between predicted and observed values can be used for both methods. Hence, the use of rank rather than Pearson correlation solves the scale issue. The remaining problem is animals with censored records, which have an unknown rank for the observed phenotype. However, the fact that animals were censored, represents an important information because those animals had the highest observed survival days. Animals with known phenotypes have a rank 1 through n, whereas censored animals have rank n+1 through N, but with an unknown order. For the censored animals, we assumed that their ranks are in random order between n+1 and N, so that the rank among the censored animals does not contribute to the estimated rank correlation. In this case, the rank correlation can be calculated by giving censored animals the average rank of all the censored animals [see Additional file 1]. In this way, we use the information that animals were censored, but make as little assumptions as possible about their order.

Before calculating rank correlations, observed phenotypes were corrected for the fixed effects [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF], ( ) P P i -. Next, for the 20% missing data, the correlation was calculated between the rank of the observed phenotypes corrected for the fixed effects, accounting for censoring as described above, and the rank of the predicted phenotypes or hazard rate: corr rank P P rank P . Furthermore, to quantify similarity of both methods, the rank correlation between the P i ∧ of 2STEP and the P i ∧ of LAM was calculated.

The rank correlation between predicted and observed phenotypes depends not only on the accuracy of the estimated breeding values underlying the predictions, but is also affected by non-genetic components of the observed phenotype. If breeding values underlying predicted phenotypes were estimated with full accuracy, the correlation between predicted and observed phenotypes would be equal to the square root of the proportion of phenotypic variance explained by breeding values,

r n A A S p D 2 2 2 2 1 = + - [ ( ) ]/    .
For any accuracy of predicted breeding values (r IH ), the expected correlation between predicted and observed phenotypes would be equal to corr r r lH = 2 (Figure 1), where r IH is the

accuracy of A A D i S j n ∧ ∧ + - ∑ , , 1
. Because animal breeders are interested in predicting breeding values rather than phenotypes, we calculated an approximate accuracy as r corr rank P P rank P r , were predicted. This accuracy is only approximate because it refers to the ranks rather than the phenotypes, and because the prediction from 2STEP refers to the scale of the hazard rate rather than the observed phenotype. For line W1 r 2 = 0.32, for line WB r 2 = 0.37, and for line WF r 2 = 0.17, when using the genetic parameters (see Table 1) given in Ellen et al. [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF].

lH i i i ∧ ∧ = - [ ( ) , (

Results

Survival

The Kaplan-Meier estimate of the survival function [START_REF] Kaplan | Nonparametric estimation from incomplete observations[END_REF] was plotted for the survival of the three layer lines (Figure 2). The survival function represents the proportion of laying hens that survived up to time t. The survival rate differed significantly between lines in both laying houses (p < 0.01). Line WF showed the highest survival rate i.e. 74.6%, whereas line WB showed the lowest survival rate i.e. 52.9%.

Genetic parameters

The estimated genetic parameters for direct and associative effects using 2STEP are given in Table 1. For all three lines, both the direct genetic variance ( )  A D 2 and the associative genetic variance ( )  A S 2 were significantly different from zero. The  A D 2 was lowest in line WF and highest in line W1, ranging from 0.12 through 0.31, whereas the  A S 2 was lowest in line WB and highest in line WF, ranging from 0.028 through 0.049. The total heritable variance  TBV 2 ranged from 0.44 (WB) through 0.81 (WF) and was significantly different from zero. Line WB showed the lowest total heritable variance in survival days expressed relative to the phenotypic variance (T 2 ), whereas line WF showed the highest T 2 , ranging from 32% through 64%. The estimated genetic correlation between direct breeding value and associative breeding value (r A ) was positive but not significantly different from zero in line W1 (0.13) and A line WF (0.55), and negative and not significantly different from zero in line WB (-0.20). Table 1 shows also the genetic parameters using LAM [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF].

Cross validation

Figure 3 shows the Kaplan-Meier estimate of the survival curves for the groups with the best and worst predicted phenotypes, using either 2STEP or LAM. As expected, for all three layer lines, groups with the best predicted phenotypes or hazard rates yielded the best observed survival curves, whereas groups with the worst predicted phenotypes or hazard rates yielded the worst observed survival curves. Both line W1 and WB showed a large difference in survival curves between the best and worst groups, whereas this difference was smaller in line WF. For all three lines, there was hardly any difference in survival curves between 2STEP and LAM. Meaning that both predicted phenotypes or hazard rates are good indicators for observed survival days. Table 2 shows the average survival days of the best and worst group for both methods and each line. Again, there was hardly any difference in average survival days between 2STEP and LAM. Furthermore, Table 2 shows the difference in survival days between the best and worst group. The difference was largest in line WB (67 days, for both methods) and smallest in line WF (16 days, for both methods). These results are in accordance with the difference in survival curves (Figure 3).

To quantify the contribution of social effects to the predicted phenotype or hazard rate, the phenotype or hazard rate is predicted using different breeding values, CBV, DBV, SBV and DSBV. Again, for each of the three lines, 25% of the animals with best predicted phenotypes or hazard rates were selected as the best group, and 25% of the animals with the worst predicted phenotypes or hazard rates were selected as the worst group. Table 3 shows the difference in survival days between the best and worst groups for both methods and each line. Besides, the Kaplan-Meier estimate of the survival curves for the groups with the best and worst predicted phenotypes based on CBV, using the two methods, is given in Figure 4. For both lines WB and W1, animals selected on predicted phenotypes or hazard rates using DSBV gave the largest difference in survival days between the best and worst group. Furthermore, for is the total heritable variance:

    TBV A A A D D S S n n 2 2 2 2 2 1 1 = + - + - ( ) ( )
;

 P 2 is the phenotypic variance:     P A A e D S n 2 2 2 2 1 = + - + ( )
, where  e 2 = 1; T 2 expresses the total heritable variance relative to the phenotypic variance:

T TBV P 2 2 2 =   /
; r A is the genetic correlation between direct breeding value and associative breeding value; r is the correlation between the residuals of group members both lines, using CBV to predict the phenotype or hazard rate gave similar difference in survival days as the DBV (approximately 45 days for line W1 and 58 days for line WB). For line WF, the difference in survival days depends on the method used. For 2STEP, the difference was largest when CBV was used, whereas for LAM the difference was largest when DSBV was used. Furthermore, for each layer line, the overlap of animals in the best group between 2STEP and LAM was calculated. The average overlap was 85% for both lines W1 and WB, and 74% for line WF. These results show that a large proportion of the animals selected for the best group, using either 2STEP or LAM, are the same.

The rank correlations between the observed phenotype, adjusted for fixed effects and censoring, and the predicted phenotype, corr rank P P rank

P i i ⎣ ⎦ ( ) , ( ) -
, are given in Table 4. The rank correlations were low and approximately the same for both methods. For line W1 (0.149 vs. 0.144) and WB (0.174 vs. 0.170), they were slightly, but not significantly, better for 2STEP, whereas for line WF (0.039 vs. 0.042) it was slightly, but not significantly, better for LAM. The rank correlations between the predicted phenotype using 2STEP and LAM were high and ranged from 0.879 (line WF) through 0.962 (line WB). These results are in line with the survival curves (Figure 3). Furthermore, approximate Best group = 25% of the animals with best predicted phenotypes or hazard rates; worst group = 25% of the animals with worst predicted phenotypes or hazard rates; phenotypes or hazard rates are predicted using CBV, DBV (A D, i ), For each subset and method, phenotypes or hazard rates were predicted based on DSBV. 25% of the animals with best predicted phenotypes or hazard rates were selected as the best groups (best refers to animals with the highest predicted phenotypes using LAM or lowest predicted hazard rates using 2STEP), and 25% of the animals with the worst predicted phenotypes or hazard rates were selected as the worst groups. Black solid line: best group using 2STEP; black dotted line: best group using LAM; gray solid line: worst group using 2STEP; gray dotted line: worst group using LAM. Results are averages of five subsets, each containing 20% of the data. Figures represent line W1 accuracies were calculated (Table 4) and were moderate, and approximately the same for both methods.

SBV A S j n ( ) , - ∑ 1 , or DSBV A A D i S j n ( ) , , + - ∑ 1 0

Discussion

We have estimated genetic parameters for direct and associative effects using 2STEP, combining survival analysis and a linear animal model including associative effects. Using 2STEP, the total heritable variance, including both direct and associative genetic effects, expressed as the proportion of phenotypic variance (T 2 ), ranged from 32% (line WB) through 64% (line WF). Using 2STEP, T 2 is substantially larger than using LAM. However, results of the cross validation do not show any difference between the two methods. Using cross validation, we showed that the difference in survival days between two extreme groups is largest when selecting on DSBV. Furthermore, we showed that social genetic effects contribute substantially to the difference in survival days (Table 3). These results indicate that there could be quite some gain in survival days when selecting on the combination of the direct breeding value and the associative breeding values of the group members (DSBV). Comparing genetic parameters of 2STEP and LAM is not straightforward. For 2STEP, genetic parameters are given on the hazard rate scale, the probability that an animal has a failure at a given time t, whereas genetic parameters of LAM are on the observed scale for survival days. The difference in genetic parameters between 2STEP and LAM originates from the fact that there is a scale difference, just like the difference in heritabilities for a 0/1trait between linear and threshold models [START_REF] Lynch | Genetics and analysis of quantitative traits Sunderland[END_REF]. Using 2STEP, the total heritable variance is 1.5 to 7-fold greater than the classical direct genetic variance using survival analysis. For both lines W1 and WB, this increase in total heritable variance is comparable with results found using LAM [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF][START_REF] Bijma | Multilevel Selection 2: Estimating the genetic parameters determining inheritance and response to selection[END_REF] (1.5 to 3-fold). For line WF, the increase is much larger using 2STEP (7-fold) than using LAM (3fold), which could be due to the fact that censoring is higher in line WF and 2STEP takes this better into account.

Theoretically, 2STEP would be a better method to analyse survival data, based on fewer assumptions known to be incorrect. We used two approaches to compare the two methods; selection of animals with the best and worst predicted phenotypes or hazard rates and the rank correlation between the predicted phenotypes or hazard rates and observed phenotypes. Both approaches show that there is hardly any difference between 2STEP and LAM. This applies to all three lines. At first glance, the difference in T 2 between both methods might suggest that using 2STEP would yield greater genetic improvement than using LAM [START_REF] Kadarmideen | Linear and threshold model genetic parameters for disease, fertility and milk production in dairy cattle[END_REF]. However, as explained above, this difference arises from For each subset and method, phenotypes or hazard rates were predicted based on CBV. 25% of the animals with best predicted phenotypes or hazard rates were selected as the best groups (best refers to animals with the highest predicted phenotypes using classical linear animal model or lowest predicted hazard rates using survival analysis), and 25% of the animals with the worst predicted phenotypes or hazard rates were selected as the worst groups. a difference in scale. The cross validation clearly demonstrates that both methods yield very similar rates of genetic improvement. Note that the rank correlation is low for all three lines, whereas the approximate accuracy is moderate for lines W1 and WB and low to moderate for line WF. Even though the approximate accuracy seems low, it is in accordance with the accuracy for methods that contain only half-or full-sib information (at least for lines W1 and WB) [START_REF] Ellen | Genetic improvement of traits affected by interactions among individuals: Sib selection schemes[END_REF][START_REF] Falconer | Introduction to Quantitative Genetics Harlow[END_REF]. Furthermore, a high rank correlation was found between the predicted hazard rates of 2STEP and the predicted phenotypes of LAM. Using selection of the best and worst predicted phenotypes or hazard rates, approximately 80% of the animals selected for the best predicted phenotypes were overlapping between 2STEP and LAM. Based on the high overlap of animals between 2STEP and LAM and the similar rank correlation, it implies that, for both methods, a similar genetic progress will be achieved.

We made a number of assumptions in the cross validation, when using the rank correlation, that may have affected the results. First, observed phenotypes were corrected for fixed effects using LAM, which may have favoured LAM compared to 2STEP. Second, when calculating the rank correlation, we assumed that ranks of censored records were in random order. This will probably not be true if censored animals were given the opportunity to actually produce a record. Alternatively, we could have used the ranks of the uncensored records only. However, in that case we would have ignored the information that the censored records are actually the "best records".

For all three layer lines, censoring occurred at the same time, at the end of the study. It could be that when censoring occurs at different times during the study period, differences may occur between the two methods. To investigate this, 50% of the censored records of line W1 were censored half way the study period (at 200 days). Again cross validation was used to compare the two methods. For both methods, the difference in survival days between the group with best predicted phenotypes or hazard rates and the group with the worst predicted phenotypes or hazard rates was calculated. For 2STEP the difference was 25.9 days, whereas for LAM the difference was 15.7 days. This indicates that, when the censoring times differ between individuals, 2STEP can better identify the genetically superior individuals than LAM. Thus both methods are practically equivalent when all animals are censored at the same survival time; with variation in censoring time, the 2STEP is superior.

Conclusion

This study shows that it is possible to use 2STEP, a combination of survival analysis and a linear animal model including associative effects, to estimate genetic parameters for the direct and associative effects on survival time in laying hens. We used cross validation to compare 2STEP with LAM. Based on the results in this paper, we can conclude that both 2STEP and LAM are practically equivalent when all animals are censored at the same survival time. Cross validation showed that selecting on a combination of the direct BV and the associative BV of the group members (DSBV) gave the largest difference in survival days between two extreme groups.

Furthermore, this study showed that social genetic effects contribute substantial to the difference in survival days between two extreme groups, which means that social genetic effects do exist.

Additional material

Additional file 1: Example and mathematical proof of rank correlation with censoring.

List of abbreviations 2STEP: two-step approach; LAM: linear animal model including direct and associative effects; CBV: classical breeding value; DBV: direct breeding value: SBV: associative breeding value; DSBV: combination of the direct breeding value of the individual itself and the associative breeding value of its group members. ; approximate accuracy = r IH = corr r / 2 , where corr is the rank correlation between observed phenotypes and predicted phenotypes or hazard rates, and r 2 = 0.32 for line W1, 0.37 for line WB, and 0.17 for line WF, when using the genetic parameters given in Ellen et al. [START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF]; results are averages of five subsets, each containing 20% of the data

∧

  denotes the predicted phenotype in case of LAM and the predicted hazard rate times -1 in case of 2STEP. Note, the P i ∧ of individual i is the sum of the estimated direct breeding value (or hazard rate) of heni AD i ( ) , ∧and the estimated associative breeding values (or hazard rates) of its group members j AS j n

Figure 1

 1 Figure 1 accuracy.

Figure 2

 2 Figure 2 Survival curve of the three layer lines. Survival curve is shown for the three lines W1, WB, and WF housed in laying house 1 (a) and laying house 2 (b).

Figure 3

 3 Figure3Survival curves using 2STEP or LAM. Kaplan-Meier nonparametric estimate of the observed survival curve of two extreme groups, based on the predicted phenotypes (LAM) or predicted hazard rates (2STEP). For each subset and method, phenotypes or hazard rates were predicted based on DSBV. 25% of the animals with best predicted phenotypes or hazard rates were selected as the best groups (best refers to animals with the highest predicted phenotypes using LAM or lowest predicted hazard rates using 2STEP), and 25% of the animals with the worst predicted phenotypes or hazard rates were selected as the worst groups. Black solid line: best group using 2STEP; black dotted line: best group using LAM; gray solid line: worst group using 2STEP; gray dotted line: worst group using LAM. Results are averages of five subsets, each containing 20% of the data. Figures represent line W1 (a), WB (b), and WF (c).

  Figure3Survival curves using 2STEP or LAM. Kaplan-Meier nonparametric estimate of the observed survival curve of two extreme groups, based on the predicted phenotypes (LAM) or predicted hazard rates (2STEP). For each subset and method, phenotypes or hazard rates were predicted based on DSBV. 25% of the animals with best predicted phenotypes or hazard rates were selected as the best groups (best refers to animals with the highest predicted phenotypes using LAM or lowest predicted hazard rates using 2STEP), and 25% of the animals with the worst predicted phenotypes or hazard rates were selected as the worst groups. Black solid line: best group using 2STEP; black dotted line: best group using LAM; gray solid line: worst group using 2STEP; gray dotted line: worst group using LAM. Results are averages of five subsets, each containing 20% of the data. Figures represent line W1 (a), WB (b), and WF (c).

Figure 4

 4 Figure 4 Survival curves based on CBV, using survival analysis or classical linear animal model. Kaplan-Meier non-parametric estimate of the observed survival curve of two extreme groups, based on the predicted phenotypes (classical linear animal model) or predicted hazard rates (survival analysis). For each subset and method, phenotypes or hazard rates were predicted based on CBV. 25% of the animals with best predicted phenotypes or hazard rates were selected as the best groups (best refers to animals with the highest predicted phenotypes using classical linear animal model or lowest predicted hazard rates using survival analysis), and 25% of the animals with the worst predicted phenotypes or hazard rates were selected as the worst groups. Black solid line: best group using survival analysis; black dotted line: best group using classical linear animal model; gray solid line: worst group using survival analysis; gray dotted line: worst group using classical linear animal model. Results are averages of five subsets, each containing 20% of the data. Figures represent line W1 (a), WB (b), and WF (c).

  Figure 4 Survival curves based on CBV, using survival analysis or classical linear animal model. Kaplan-Meier non-parametric estimate of the observed survival curve of two extreme groups, based on the predicted phenotypes (classical linear animal model) or predicted hazard rates (survival analysis). For each subset and method, phenotypes or hazard rates were predicted based on CBV. 25% of the animals with best predicted phenotypes or hazard rates were selected as the best groups (best refers to animals with the highest predicted phenotypes using classical linear animal model or lowest predicted hazard rates using survival analysis), and 25% of the animals with the worst predicted phenotypes or hazard rates were selected as the worst groups. Black solid line: best group using survival analysis; black dotted line: best group using classical linear animal model; gray solid line: worst group using survival analysis; gray dotted line: worst group using classical linear animal model. Results are averages of five subsets, each containing 20% of the data. Figures represent line W1 (a), WB (b), and WF (c).

Table 1

 1 Estimates of genetic parameters for direct and associative effects on survival time in three layer lines using 2STEP or LAM[START_REF] Ellen | Survival of laying hens: Genetic parameters for direct and associative effects in three purebred layer lines[END_REF]

			W1		WB		WF
			2STEP	LAM	2STEP	LAM	2STEP	LAM
	2  A D 2  A S 2  TBV 2  P T 2	0.31 ± 0.05 0.041 ± 0.01 0.77 ± 0.13 1.44 ± 0.06 0.53 ± 0.08	915 134 2,490 12,847 0.19	0.30 ± 0.05 0.028 ± 0.01 0.44 ± 0.09 1.38 ± 0.05 0.32 ± 0.06	1,917 273 3,007 20,111 0.15	0.12 ± 0.06 0.049 ± 0.02 0.81 ± 0.26 1.27 ± 0.08 0.64 ± 0.17	246 60 910 13,999 0.06
	r A		0.13 ± 0.15	0.18	-0.20 ± 0.14	-0.31	0.55 ± 0.28	0.11
	r		-0.003 ± 0.0003	0.08	-0.005 ± 0.0001	0.08	-0.004 ± 0.0003	0.10
	 A D 2	and  A S 2 are estimates of direct genetic variance and associative genetic variance;  TBV 2

Table 2

 2 Mean survival days of best and worst groups using 2STEP or LAM for three layer lines

		W1		WB		WF
		2STEP	LAM	2STEP	LAM	2STEP	LAM
	Mean	354 ± 2	326 ± 2	375 ± 2
	Best	377 ± 3 377 ± 3 359 ± 2 357 ± 3 384 ± 3 383 ± 5
	Worst	327 ± 2 327 ± 3 292 ± 6 290 ± 6 368 ± 7 367 ± 5
	Difference	50	50	67	67	16	16
	Best group = 25% of the animals with best predicted phenotypes (LAM) or
	hazard rates (2STEP); worst group = 25% of the animals with worst predicted
	phenotypes or hazard rates; difference = mean survival days of best group -
	mean survival days of worst group; results are averages of five subsets, each
	containing 20% of the data

Table 3

 3 Difference in survival days between best and worst groups using 2STEP or LAM for three layer lines

		W1		WB		WF
		2STEP	LAM	2STEP	LAM	2STEP	LAM
	CBV	46	43	58	57	19	14
	DBV	45	43	59	57	15	11
	SBV	25	26	32	33	12	9
	DSBV	50	50	67	67	16	16

Table 4

 4 Rank correlation and approximate accuracy based on 2STEP or LAM for three layer linesRank correlation is calculated between observed phenotypes and predicted hazard rates of 2STEP, between observed phenotypes and predicted phenotypes of LAM, and between predicted phenotypes and predicted hazard rates (2STEP; LAM); observed phenotype = phenotype corrected for fixed effects (

									Rank correlation	Approximate accuracy
								W1	WB	WF	W1	WB	WF
	2STEP					0.149 ± 0.011	0.174 ± 0.020	0.039 ± 0.017	0.47	0.47	0.22
	LAM					0.144 ± 0.010	0.170 ± 0.020	0.042 ± 0.012	0.45	0.46	0.24
	2STEP; LAM					0.954 ± 0.003	0.962 ± 0.004	0.879 ± 0.007	-	-	-
									) i -; P P
	predicted phenotype or hazard rate ( P i ∧ ) = sum of estimated direct breeding value of hen i AD i ( , ∧	)	and estimated associative breeding values of its group
	members j A	S j	P	i	A	D i	A	S j
	n					n

Acknowledgements

We would like to thank the employees of the laying houses for taking good care of the hens and for collecting the data. Johan van Arendonk is acknowledged for helpful comments on earlier versions of the manuscript. This research is part of a joint project of Institut de Sélection Animale B.V., a Hendrix Genetics Company, and Wageningen University on 'Genetics of robustness in laying hens', which is financially supported by SenterNovem. Both EDE and PB are financially supported by the Dutch science council (NWO) and part of this work was co-ordinated by the Netherlands Technology Foundation (STW).

Author details 1 Animal Breeding and Genomics Centre, Wageningen University, Marijkeweg 40, 6709PG Wageningen, The Netherlands. 2 UMR 1313 GABI, INRA, 78352 Jouy-en-Josas, France. 3 Animal Breeding and Genomics Centre, Wageningen UR Livestock Research, 8200AB Lelystad, The Netherlands.

Authors' contributions EDE performed the data analysis and the cross validation, wrote and prepared the manuscript for submission. VD helped with the data analysis and cross validation and reviewed the manuscript. BJD helped with the data analysis and reviewed the manuscript. RFV helped with the data analysis and reviewed the manuscript. PB was the principal supervisor of the study and assisted with data analysis, cross validation and preparation of the manuscript. All authors read and approved the manuscript.

Competing interests

The authors declare that they have no competing interests.