A GENERAL BIVARIATE WEIBULL MODEL APPLIED TO SHELF …

A GENERAL BIVARIATE WEIBULL MODEL APPLIED TO SHELF LIFE DETERMINATION OF FOOD PRODUCTS

Marta Afonso FreitasDepto. de Engenharia de Produção – Escola de Engenharia – UFMG

LADEC – Laboratório de Apoio à Decisão e Confiabilidade.Av. Antonio Carlos 6627- Cep 30161-970 – Belo Horizonte - MG

[email protected]

Fernanda Carla Álvares Dias Depto. de Estatística- ICEx - UFMG

Av. Antonio Carlos 6627- Cep 30161-970 – Belo Horizonte - [email protected]

RESUMO

Quando se fala da determinação da vida de prateleira de produtos alimentícios, análises química, físicas, microbiológicas e nutricionais são fundamentais, mas, igualmente importantes, são as características sensoriais do produto. Por esta razão, o papel de avaliações sensoriais para determinação de vida de prateleira está se tornando mais importante. Nestes experimentos, uma amostra de unidades do produto é armazenada sob certas condições e periodicamente, em tempos pré-especificados uma amostra de unidades é retirada e submetida a avaliações sensoriais por painelistas treinados. Devido a natureza destrutiva destas avaliações, nunca se sabe o exato “tempo de falha” de uma dada unidade, levando a observações censuradas à direita ou à esquerda. Neste artigo apresentamos um modelo para determinação da vida de prateleira. Os dados são dicotomizados e a distribuição subjacente dos “tempos de falha” de dois atributos é modelada por uma distribuição de Weibull bivariada. A abordagem é aplicada a dados reais.

PALAVRAS CHAVE. Distribuição Weibull bivariada. Vida de Prateleira. Confiabilidade. Área de classificação principal : Aplicações à Indústria.

ABSTRACTWhen one talks about determining the shelf life of food products chemical, physical microbiological and nutritional analysis are fundamental but equally important are the sensory characteristics of the product. For this reason, the role of sensory evaluations for shelf life determination is becoming more important. In such experiments, a sample of product units is stored under certain conditions and periodically, at pre-specified evaluation times a sample of units is taken from the ones stored and subjected to sensory evaluations by trained panelists. Due to the destructive nature of these evaluations, one never knows the exact “failure time” of a given unit, leading to either right or left censored observations. In this article we present a model for shelf life determination. The data is dichotomized and the underlying “failure times” of two attributes are jointly modeled by a bivariate Weibull distribution. The approach is applied to a real data set.

KEYWORDS. Bivariate Weibull distribution.Reliability.Shelf life. Main area: Applications to the Industry.

XLI SBPO 2009 - Pesquisa Operacional na Gestão do Conhecimento Pág. 648

1. Introduction Problems related to storage stability are common to the food industry and that is why storage studies are an essential part of every product development, improvement, or maintenance program. Some studies center on the rate of degradation, and others on estimating the shelf life: the length of time required for the product to be unfit for human consumption. By unfit for human consumption it is meant that the product exhibits either physical, chemical, microbiological, or sensory characteristics that are unaccepted for regular consumption.

The manufacturer attempts to develop a product with the longest shelf life practical, consistent with costs and the pattern of handling and use by distributors, retailers, and consumers. Inadequate shelf life determination will lead to consumer dissatisfaction or complains. At best, such dissatisfaction will eventually affect the acceptance and sales of brand name products. At worst, it can lead to malnutrition an even illness. For these reasons, food processors pay great attention to adequate storage stability or shelf life. When one talks about determining the shelf life, chemical, physical microbiological and nutritional analysis are fundamental but equally important are the sensory characteristics of the product. It is a fact that many consumers purchase a product on the basis of the sensory experience which it delivers for instance, sweetness, softness, chocolateness, aspect, odor, flavor, aftertaste, etc. For this reason, the role of sensory evaluations for shelf life determination is becoming more important as the value of these tests becomes recognized and as the consumer industry increases its use of these methods. In such experiments, a sample of product units is stored under certain conditions and periodically, at pre-specified evaluation times a sample of units is collected from the ones stored in a given condition and subjected to sensory evaluations by trained panelists. A frequent test procedure is to ask each panelist to judge each attribute separately by reference to a rating scale, for instance, a seven point rating scale varying from 0 (“no difference”) to 6 (“total difference”). Each “test unit” is compared to a “control unit” and the scores are the result of this comparison. Because of the destructive nature of evaluations, units that have already been evaluated at a given time cannot be restored to be evaluated later on. Consequently data coming from these sensory evaluations are either left or right censored (Meeker and Escobar, 1998).

An usual statistical approach to the analysis of these kind of data is to fit a parametric lifetime model such as Weibull, lognormal, to the “failure” data. Examples of this approach can be found in Gacula and Kubala (1975) and Gacula and Singh (1984). Evaluation times (which are in fact fixed) are used as approximations to the real failure times for those units scored less than the cut-off point. A parametric model is then fitted to the “approximated” and to the right censored failure times. Two problems with this procedure can be pointed out here: (1) the parametric failure time models are not being used appropriately since the evaluation times are fixed and, as a consequence, (2) there is a great risk of overestimating the shelf life. In an attempt to minimize the effect of these two mentioned problems, Gacula (1975) suggested the use of a “staggered design” to implement weekly sensory evaluations. In that design the intervals between evaluations are larger at the beginning of the experiment and smaller towards the end. The number of units sampled increases towards the end of the experiment. The idea behind that is to concentrate the majority of units in the experimental periods where the maximum information is desired. This time is the period where the units are likely to “fail” or are approaching the limit of acceptable quality. With this procedure the author intended to reduce the gap between the fixed evaluation times and the real failure times and also get a cut-back in experimentation cost.

A different approach was developed by Freitas, Borges e Ho (2003). At each (fixed) evaluation time ),,1( kii =τ , the authors defined a new random variable ijY , Bernoulli distributed with probability of “failure” (“unit considered unfit”) given by

0()1( PYPp ijij === < )iijT τ≤ )(1 ijR τ−= where ijT is the failure time of the thj unit

evaluated at iτ , with reliability function )()( tTPtR ijj >= . The authors went one step further


and assumed that the underlying distribution of ijT was a Weibull with scale and form

parameters given respectively by { }0 1 1expj j jq qX Xα β β β= + + + and { }exp , 0δ γ γ= ≥ , with 1(1, , , )j j jqX X X= being a )1(1 +× q vector of variables (covariates or factors of a planned experiment).

The likelihood function was then written for the Bernoulli random variables, incorporating the equations for jα and δ given by the underlying Weibull distribution. Parameter estimates were obtained by maximum likelihood and the shelf life time distribution was characterized through percentiles estimates. The authors applied the model to real data coming from sensory evaluations performed on a food product. Further in Freitas and Costa (2006) , the basic model developed by Freitas, Borges and Ho (2003) was expanded allowing the Weibull form parameter δ to depend on variables (covariates or factors of a planned experiment).

All the approaches mentioned so far have been trying to deal with the fact that the results of the experiments are either left or right censored. But right and left censored data are in fact particular cases of interval censoring (Kalbfleisch and Prentice, 2002). So, without loss of generality, we can say that the data arising from sensory evaluations for shelf life determination are interval censored. Other examples of interval-censored data occur naturally when a lifetime response arises from a clinical trial or longitudinal studies that entail periodic follow-up (Kim, Gruttola and Lagakos, 1993; Sun, 1996). Some methods have been proposed for the analysis of interval-censored data, for instance Turnbull (1976), Peto (1973), Finkelstein (1986), and Rabinowitz (1995). A flexible family of regression models where the underlying hazard function is estimated using B splines was presented by Younes and Lachin (1997). It is worth noting though, that lifetime data with interval censoring are usually broadly grouped and therefore, methods for discrete models are somehow natural and have more appeal. Regression methods for grouped data were presented by Lawless (1982) and Collet (2003). The regression structure considers modeling the probability of a subject´s ou experimental unit´s “survival” past a visit or evaluation time conditional on his survival at the previous evaluation. Two approaches were presented: (1) assuming that lifetimes come from a continuous proportional hazards model (Prentice and Goeckler, 1978) and (2) assuming they come from a proportional odds model. Freitas and Gomes (2006) developed the expressions of a regression discrete model where the probability ( 1) (0ij ijp P Y P= = = < )ij iT τ≤ 1 ( )j iR τ= − was modeled adopting the proportional hazards model (Collet, 2003; Cox, 1972).

It is important to recall that in the case of the sensory evaluations of food products, more than one attribute is evaluated at the same time. The usual practice is to model “the failure time” of each one separately. The product´s shelf life is then determined based on chosen percentiles of the “failure time” distribution of each of them, usually by adopting the “worst case”. Note that all the approaches reviewed in this section can be used to do that. Nevertheless it would be interesting to be able to model more than one attribute jointly. This approach would allow one to first, account for possible correlations between the attributes (a situation that might occur in practice) and second, by doing so, obtain a better estimate of the product´s shelf life. In addition, this strategy would allow one to understand a little better the effect of the correlations in the final estimates of the shelf life.

In this paper, we present an approach for modeling the “failure times” of two attributes jointly. We use the bivariate Weibull model proposed by Hougaard (1986) combined with some of the previous ideas of Freitas, Borges and Ho (2003). A simulation study is also implemented in order to understand the effect of correlations on the estimate of the fraction of products unfit for consumption (“fraction of failures”). The model is applied to the data set presented in Freitas, Borges and Ho (2003).

The paper is organized as follows. Section 2, presents briefly the real situation described by Freitas, Borges and Ho (2003) . Section 3 introduces the bivariate modeling approach. Section 4 presents the results of the simulation study implemented in order to evaluate the proposed


model. Section 5 shows the results of the of application of the proposed model to the practical situation described in Section 2. Section 6 presents conclusions and suggests directions for future research.

2. Motivating situation revisited In this section we briefly present the situation described by Freitas, Borges and Ho (2003). Sensory evaluations were conducted by a food company at laboratory level in order to

determine the shelf life of a manufactured product, stored at three environmental conditions. Three attributes were evaluated by trained panelists: odor, flavor and general appearance. The main characteristics of the study are presented next.

• Experimental design : a lot of units was sampled from the product line and were randomly assigned to one of the following storage conditions:

1. refrigeration (approximately 4 o C): these units were used as reference (control) during the trials;

2. room temperature and humidity: temperature and humidity were monitored and registered continuously. Average weekly values were calculated and used in the study;

3. environmental chamber 1: temperature and humidity levels controlled at 30 o C and 80% respectively;

4. environmental chamber 2: temperature controlled at 37 o C. Humidity levels were not controlled but they were registered daily and average weekly values were computed

• Laboratory panel : 5 to 8 trained panelists were assigned to each week. • Follow-up time : room temperature, 51 weeks. chamber 1, 36 weeks and chamber 2, 18

weeks.• Test procedure : sensory evaluations were performed weekly. Scores were assigned to each

attribute separately by each one of the panelists designated to that particular week. Scores were based on a 7-point “off-flavor” scale. Units were compared with a reference (“ 0”, total difference; “6” no difference).

• Cut-off point : score “3”. If for a given unit, a particular attribute receives a score less than or equal 3 then that unit is considered inadequate for consumption regarding that particular attribute.

In the previous work, percentiles and fraction of defectives were estimated at different points in time for each one of the attributes separately. The purpose here is to model two attributes jointly.

3. The Bivariate Parametric Model for Interval Censored Data, based on a Bivariate Weibull Distribution.

In the previous methods discussed in Section 1, the study of the (shelf) life distribution was implemented according to a univariate approach. In other words, the results of the sensory evaluations (in the case of our practical situation) were modeled separately. The purpose now is to model two attributes jointly in order to take into account possible correlations that might exist between them. The model is presented next in the context of the practical situation.

Let ijkZ be the score assigned to the thk attribute ( 1, , )k r= of the thj product unit (

1,2, , )sj n= ⋯ evaluated at time , ( 1,2, , )i i sτ = . The thj product unit evaluated at iτ is

considered unfit for consumption (regarding the attribute being checked) if ijkZ c≤ and considered adequate for consumption otherwise. The value “c” is a cut-off point in the scale, which is usually defined by the food company.


Then, a new random variable ijkY is defined, through the dichotomization of the scores. In

other words, at each evaluation time ( 1, , )i i sτ = , 1ijkY = if the score of the thj( 1, , )sj n= unit evaluated at iτ is less than or equal “c” (the cut-off point) and 0ijkY = otherwise. Recall that, in our case, 2r = (two attributes will be modeled jointly) and c=3 (cut-off point).

Therefore, what we observe at each time iτ , is a sample of size sn from a random vector

1; 2( )ij ijY Y , bivariate-Bernoulli distributed with probabilities (00) ; (01) ; (10) ; (11)andij ij ij ijp p p p , in other

words, 1; 2 (00) ; (01) ; (10) ; (11)bivariate Bernoulli( ) ( )tij ij ij ij ij ijY Y p p p p:

where (00) (01) (10) (11) 1ij ij ij ijp p p p+ + + = and the bivariate probability function is given by1 2 1 2 1 2 1 2(1 )(1 ) (1 ) (1 )

1 1 2 2 (00) (01) (10) (11)( , ) ij ij ij ij ij ij ij ijy y y y y y y yij ij ij ij ij ij ij ijP Y y Y y p p p p− − − −= = = (1)

( ijijij ppp )11()10(1 += and ijijij ppp )11()01(2 += ).

It is also possible to write the probabilities in (1) as:

(11) 1 2 1 2 1 2( , ) ( 1, 1) (0 ,0 )ij ij ij ij ij ij i ij ip P Z c Z c P Y Y P T Tτ τ= ≤ ≤ = = = = ≤ ≤ ≤ ≤ (2)(10) 1 2 1 2 1 2( , ) ( 1, 0) (0 , )ij ij ij ij ij ij i ij ip P Z c Z c P Y Y P T Tτ τ= ≤ > = = = = ≤ ≤ > (3)

(01) 1 2 1 2 1 2( , ) ( 0, 1) ( , 0 )ij ij ij ij ij ij i ij ip P Z c Z c P Y Y P T Tτ τ= > ≤ = = = = > ≤ ≤ (4)

(00) 1 2 1 2 1 2( , ) ( 0, 0) ( , ) ( ; )ij ij ij ij ij ij i ij i i ip P Z c Z c P Y Y P T T Rτ τ τ τ= > > = = = = > > = (5)

where ijkT is the “failure time'” or “shelf life time” associated to the thk attribute ( 1,2)k = of

the thj unit ( 1, , )ij n= evaluated at iτ and ( ; )i iR τ τ in (5) is the jointly reliability function, evaluated at iτ ( 1, , )i s= . The likelihood function takes then the following form:

1 2 1 2 1 2 1 2(1 )(1 ) (1 ) (1 )(00) (01) (10) (11)

1 1

( )i

ij ij ij ij ij ij ij ij

nsy y y y y y y y

ij ij ij iji j

L p p p p− − − −

= =

= ∏ ∏θ (6)

where the probabilities are given by the expressions (2). to (5) and t t t t t00 01 10 11θ = (θ ,θ ,θ ,θ ) is

the parameter vector with dimension 41

1s

ii

n=

×

∑

Note that in (5), a relationship between the underlying distribution of the random vector

1 2( , )tij ijT T and the scores actually observed has been established. Therefore in order to be able to

do go further on the model development, it was necessary to search for families of bivariate distributions. It seemed reasonable to restrict our search to those distributions whose marginal distributions have proven satisfactory for describing a single univariate property (in our case, a single shelf life time for one attribute). One potential category of bivariate distributions is the family of bivariate Weibull distributions. That was our choice since the univariate Weibull distribution is widely used to model lifetimes and, in addition, the previous univariate modeling approaches for the sensory evaluation data used Weibull distributions. Gumbel (1960) was the first one to derive a Weibull bivariate model and the bivariate reliability function associated to this model. In this article, we consider the bivariate Weibull model proposed by Hougaard (1986). The marginal distributions are well understood and the reasonableness of this bivariate model has been confirmed in several practical situations. The marginal distributions are two-parameter Weibull distributions. We will be using basically the same idea developed in Freitas, Borges and Ho (2003). The authors dichotomized the results of the sensory evaluations and used univariate


Weibull distributions as the underlying distributions of the shelf life times of each attribute. Here, we will use the same idea but, a bivariate Weibull will be the underlying distribution.

Let us assume that the random vector of “failure times” or “shelf life times” 1 2( , )tij ijT T for

two attributes considered has a bivariate Weibull distribution with parameters 1 2 1 2, , , ,j j j jα α δ δ λ . The marginal distributions are each a Weibull ( jkα jkδ ) ( 1,2k = ) where

jkα and jkδ are the scale and shape parameters respectively. The probability density function is given by:

( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

2/ 1 / 1 / /1 2 1 2( , ), 1 2 1 1 11 2 2 2 2 11 2 21 2

/ / / /11 2 1 21 exp11 2 2 11 2 2

j j j jf t t t t t tT T j j j j j j j jij ij

j j j jt t t tj j j j

λδ λ δ λ δ λ δ λ

δ α α δ α α α α

λδ λ δ λ δ λ δ λ

α α α αλ

−− −

= +

+ + − − +

( ) 7

λ

Lu and Bhattacharyya (1990) give the expressions of the moments associated to the bivariate Weibull in (7). Note that the index “j” is associated with the thj ( 1, , )ij n= experimental unit evaluated at a

fixed time iτ ( 1, , )i s= . This notation allows one to incorporate possible differences between the experimental units. For example, the scale parameters might be modeled as functions of a vector

1( , , )tj j jpX X=X of covariates or factors of designed experiments associated to each specific unit j

(i.e., temperature, humidity, lot, supplier, etc.). In our case, we are assuming that all the experimental units under evaluation were are under the same conditions. Therefore, for the sake of simplicity, the index “j” will not be used in the remaining of this article when referring to the model parameters. The fifth parameter λ is related to the covariance between the two components of the random vector 1 2( , )t

ij ijT T (see Lu and

Bhattacharyya, 1990 expression for the expression of the covariance). If 1λ = then 1 2cov( , ) 0ij ijT T = .

In addition, if 1λ = in (7), then

( ) ( ){ } ( ) ( ){ } ( ) ( )1 111 1 2 2( , ) exp exp , 1 2 1 1 1 1 1 2 2 2 2 2 1 21 2 1 2(8)f t t t t t t f t f tT T T Tij ij ij ij

δδ δ δ δα δ α α δ α

− −= − − =

where ( )1 1ijTf t and ( )

2 2ijTf t are the probability density functions of random variables Weibull

distributed, with parameters 1 1( )α δ and 2 2( , )α δ . The correlation, 1 2co ( , )ij ijrr T T , is a

complicated function of 1,λ δ and 2δ . Lu and Bhattacharyya (1990) looked at the correlation when 1 2δ δ δ= = for 0.1,0.25,0.5,1,2,4δ = and 10 and showed that the correlation is almost linear in δ for equal shape parameters 1 2 1δ δ= = , which could enhance estimation and/or

simulation procedures. The vector of parameters 1 2 1 2( , , , , )tα α δ δ λ=η can be estimated by the Maximum Likelihood (ML) method (Casella and Berger, 2001), through maximization of the log-likelihood, where the likelihood function is given by (6). This maximization can be obtained by implementing suitable maximization algorithms such as the Newton-Raphson algorithm (Casella and Berger, 2001). The ML point estimates of the probabilities given by equations (2) to (5) can be obtained though the invariance property of the ML estimators. In other words it is possible to get the estimate ijp )00(ˆ of the joint reliability and of the “fractions of units unfit for

consumption” ijp )01(ˆ , ijp )10(ˆ , and ijp )11(ˆ through the expressions:


( ) ( )ˆ

ˆ ˆ ˆ ˆ/ /1 2ˆ ˆ ˆˆ ( , ) exp 1 2(00)p R i i i iij

λδ λ δ λ

τ τ α τ α τ

= = − +

(9)

1 2 ˆ(01) , 1 2 1 20

ˆ ( , ) |i

i

ij T Tp f t t dt dtτ

τ

∞

== ∫ ∫ η η 1 2 ˆ(10) , 1 2 1 20

ˆ ( , ) |i

i

ij T Tp f t t dt dtτ

τ

∞

= ∫ ∫ η=η 1 2 ˆ(11) , 1 2 1 2

0 0

ˆ ( , ) | i i

ij T Tp f t t dt dtτ τ

= ∫ ∫ η=η (10)

where 1 2 1 2

ˆ ˆ ˆˆ ˆ ˆ( , , , , )tα α δ δ λ=η is the ML point estimate of 1 2 1 2( , , , , )tα α δ δ λ=η . Confidence intervals for the quantities of interest can be obtained by the parametric resampling method Bootstrap (Efron and Tibshirani, 1986). However, the calculations necessary to obtain the ML estimator in this case are not straightforward. In order to get the analytical expressions of the reliability and “fraction of units unfit for consumption” as functions of η̂ (equations (9) and (10)) it is necessary, among other tasks, to integrate quite complicated functions. In addition, those integrals will end up in very extent equations whose first and second derivatives will have to be calculated in order to construct the Hessian Matrix, necessary in order to implement the Newton Raphson Method. These calculations would require additional computational efforts, especially during the simulation study (see Section 4). Consequently, a new parameterization was proposed. A new random variable ijW was defined as a result of the combination of the scores obtained by each one of the two attributes under consideration. The details are given in the next section.

3.1 Definition of the variable ijW

Let ijW be the random variable defined as:th

i

thi

1, if the j unit is adequate for consumption at time regarding the two attributes (11)

0, if the j unit is for consumption at regarding of the attributeijW

τ

τ=

inadequate at least one s

From the above definition we have:1 2 1 2 1 2 (00)( 1) ( , ) ( 0, 0) ( , )ij ij ij ij ij ij i ij i ijP W P Z c Z c P Y Y P T T pτ τ= = > > = = = = > > = .

Consequently, ijW has a Bernoulli distribution, with probability of “success” ijp )00( and the

following density (mass) function ( ) ( ) ( ) 1

( 00) ( 00 )1ij ijW W

ij ij ij ijP W w p p−

= = − . The log-

likelihood function takes now the form:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ){ }1 2 1 2/ / / /1 2 1 2

1 1

log ( ) 1 log 1 expins

ij i i ij i ij j

L W Wλ λδ λ δ λ δ λ δ λα τ α τ α τ α τ

= =

= − + + − − − + ∑ ∑η (12)

where ( )1 2 1 2, , , , tα α δ δ λ=η is the parameter vector to be estimated. As it was already

mentioned, the parameter vector ( )1 2 1 2, , , , tα α δ δ λ=η is estimated by ML using the log-

likelihood (12) and the point estimates of (00) 2 2( , ) ( , )ij ij i ij i i ip P T T Rτ τ τ τ= > > = and (00)1 ijp−

are obtained by substituting 1 2 1 2ˆ ˆ ˆˆ ˆ ˆ( , , , , )tα α δ δ λ=η (the ML point estimate) into equation (9).

Confidence intervals for the quantities of interest can be obtained by the parametric resampling method Bootstrap (Efron and Tibshirani, 1986). Note that, with this new parameterization it is possible to estimate the fraction of items unfit for consumption (according to the evaluation performed by the panelists) in a given (fixed) evaluation time, in other words, the fraction of products whose both attributes became unfit for consumption. In fact, according to the food company, this is the information necessary to establish a shelf life period for the product, as far as


sensory attributes are concerned. On the other hand, it is not possible to know which one of the attributes became “unfit” first. In the next section, the results of a simulation study are presented.

4. Simulation StudyThe purpose of the simulation study is two-fold. First, to study the “quality” of the estimates

of fraction of units unfit for consumption (or “fraction of defectives”) obtained with the model presented in Section 3. The basic idea is to generate observations under known conditions and use the proposed model to estimate the “fraction of defectives” at different points in time. By “known conditions” it is meant situations which emulate the behavior followed by the real data available. Those estimates are then compared to the real values. Second, to evaluate the effect of different levels of correlation among the attributes failure times on the estimated “fraction of defective units”. The results obtained with the bivariate model proposed in Section 3 are compared to the ones obtained by univariate models fitted to each one of the attributes independently. In the case of the univariate analysis, the fraction of defectives ( (00) 2 21 1 [ , ]ij ij i ij ip P T Tτ τ− = − > > ) is calculated as if they were uncorrelated. The details of the study are given next.

4.1 Parameter values used in the simulation.Table 1 presents the six scenarios considered. The values of jα and jδ ( 1,2)j = used in

the simulation are the same ones used by Freitas, Borges and Ho (2003) . The values were chosen in an attempt to make the degradation mechanism of the attribute 2 (j=2) more “aggressive” than the one for attribute 1 (j=1). The values of λ were chosen based on the work by Johnson, Evans and Green (1999). Recall that 1λ = implies zero covariance (zero correlation) among the failure times of the two attributes considered. On the other hand, 0λ = implies a perfect correlation (

1)ρ = between the two failure times. For each one of the cases listed in Table 1, 1000 samples of a bivariate Weibull distribution are generated. Details of the procedure used to generate bivariate Weibull distributions can be found in Lu and Bhattacharyya (1990). 4.2 Steps followed in the simulation study In the proposed model, the underlying distribution of the jointly failure time (shelf life time) is a bivariate Weibull but in a real situation, one does not observe the actual failure times. In fact only information available is the score assigned by a panelist to each one of the attributes of a given product at some (fixed) point in time. In the model proposed in Section 3, the results are dichotomized according to the cut-off point established by the company. In other words, the result is either one or zero depending on the failure times 1ijT and 2ijT being both located after

the evaluation time iτ or not. In the simulation study we assume that the evaluations are implemented weekly and the total follow-up period used is 36 (thirty six) weeks. It is also assumed that 7 (seven) panelists were requested weekly to compose the laboratory panel. This follow-up period was chosen since the ones observed in the real experimental data were 18, 32 and 51 weeks depending on the storage condition. We decided to fix a number located somewhere in the middle range.

Table 1 – Parameter values used in the simulation study.Scenarios Parameter values

α1 α2 δ1 δ2 λ1 0.035 0.035 1.2 1.6 0.22 0.035 0.035 1.2 1.6 0.53 0.035 0.035 1.2 1.6 0.94 0.025 0.035 1.2 2 0.25 0.025 0.035 1.2 2 0.56 0.025 0.035 1.2 2 0.9

4.3 Simulation Results Estimates of the “fraction of units unfit for consumption” were obtained for each one of the scenarios

listed in Table 1, for weeks 1 to 36, using the bivariate model and the independency assumption. Table 2 presents the results obtained for the 24th week. The table displays also the correlation values (real and estimated by the bivariate model) for each one of the six scenarios and the estimated fraction of units unfit,


calculated by univariate models fitted to each one of the attributes separately. Figures 1 (a) and (b) present the curves of the estimated “fraction of units unfit” vs. weeks for scenarios 1 and 2 .The figures show also the real values of those quantities. The estimated values of the “fraction of unfit units” obtained based on the assumption of zero correlation, were calculated using the estimated values produced by univariate models fitted to each one of the attributes separately, in other words, using the following expression:

1 1 ( ; ) 1 ( ) ( ) 1 [1 ( )][1 ( )]1 2 1 200( )fraction at p R R R F Fi i i i i i iiτ τ τ τ τ τ τ= − = − = − = − − − (13)

where, ( ) ( )1 1F P Ti iiτ τ= ≤ and ( ) ( )2 2F P Ti iiτ τ= ≤ are the “fraction of unfit units” at week iτ for attributes 1 and 2 respectively. Therefore, for scenario 1, the value 0.7906 for the “fraction of unfit units” at 24 24τ = wks was obtained as ( ) ( )0.7906 1 1 0.5588 1 0.5254= − − − Table 2 - Fraction of unfit at the 24th week (real and estimated values using the bivariate model and the independency assumption).

Quantities of interest(real and estimated values)

Scenarios

[1]

(λ=0.2)

[2]

(λ=0.5)

[3]

(λ=0.9)

[4]

(λ=0.2)

[5]

(λ=0.5)

[6]

(λ=0.9)

Correlation (real value) 0.9124 0.6154 0.1221 0.9047 0.6200 0.1261Correlation (estimated –bivariate

model )0.9198(0.0352) (1)

0.6360(0.0253)

0.2618(0.0321)

0.8732(0.0431)

0.6091(0.1030)

0.2075(0.0248)

Fraction of unfit (real ) 0.5945 0.6702 0.7684 0.5228 0.5892 0.6880

Fraction of unfit (estimated value - bivariate model)

0.5967(0.0346)

0.6711(0.0305)

0.7415(0.0694)

0.5423(0.0378)

0.5996(0.0429)

0.6768(0.0627)

Fraction of unfit (estimated value - independent fit (assuming zero covariance)

0.7906(0.1250)

0.7931(0.1020)

0.7958(0.1101)

0.7178(0.1315)

0.7191(0.1421)

0.7112(0.1022)

Fraction of unfit (estimated value for attribute 1 – univariate model)

0.5588(0.0346)

0.5570(0.0306)

0.5598(0.0285)

0.4294(0.0400)

0.4244(0.0412)

0.4156(0.0357)

Fraction of unfit (estimated value for attribute 2 – univariate model)

0.5254(0.0349)

0.5330(0.0372)

0.5362(0.0354)

0.5055(0.0420)

0.5120(0.0381)

0.5056(0.0369)

(1)= Standard deviation

(a) (b) Figure 1 – Fraction of unfit units (“fraction of defectives”) vs. evaluation weeks : (a) scenario 1 -

0, 035; 1, 2; 1, 6; 0, 21 2 1 2α α δ δ λ= = = = = ⇒ correlation =0.9124; (b) scenario 2 -

0, 035; 1, 2; 1, 6; 0, 51 2 1 2α α δ δ λ= = = = = ⇒ correlation =0.6154).From the analysis of the results presented in Table 2 and Figures 1 (a) and (b) some important aspects

can be highlighted: (1) for cases 1 and 2 where 0 0.5λ≤ ≤ (1 ≤ correlation ≤ 0.5) the curve of the estimated values of the “fraction of unfit units” obtained under the assumption of no covariance among the attributes (independent fit) is well above the other two curves. In other words, it seems that if the two attributes have shelf life times strongly correlated one always gets higher estimated values for the fraction of unfit units when using the independent fit (see values on Table 2 and Figures 1 (a) and (b). In addition, the curve of the estimated values calculated with the bivariate model is very close to the curve of the real values, indicating that the model is providing good point estimates in those cases. The standard deviation values do not show any specific pattern indicating that the precision


of the estimates is not associated with higher or lower covariance levels. The same observation can be made about the mean square error values (not shown here); (2) for very low values of correlation, the estimated values provided by the independent univariate fits are very close to the real ones. For the specific value used in the simulation ( 0,9λ = ⇒ correlation=0.1221) the curve of the estimated values is still above the other two but, as it was mentioned before, it is much closer to the other ones. This is result was already expected since the calculations based on the independent fits are a special case of the bivariate fit when 1λ = (correlation=0). This figure is not provided here but it is available upon request.

The simulation results for the other cases are similar and will not be presented here. These preliminary results indicate that if the level of correlation between the shelf life times of the attributes is low, one does not loose much information by performing the univariate analyses independently and then combining the results at the end in order to calculate the fraction of unfit units. Unfortunately, the practitioners still have to deal with the problem that in real situations like the one considered in this paper, one does not observe the actual failure times. Therefore it is not easy to get an estimate of the correlation beforehand. For cases like the one described in Section 2, one strategy might be to calculate the sample correlation between the scores values assigned to each attribute. That value might indicate the level of correlation between the non observable shelf life times and provide additional information to the researcher, in case he (she) intends to choose one modeling approach. In the next section, the results of the analysis of the shelf life data described in Section 2 will be presented. The data analyzed refers to the storage condition “environmental chamber 2” (temperature controlled at 37 o C) .

5. Application.

In this section, we return to the practical situation in Section 2. We present the results of the data analysis for the storage condition “environmental chamber 2”, using the scores assigned to the attributes “odor” and “flavor”. First, we give the results based on the dichotomized model with the bivariate Weibull as the underlying distribution. Then, we compare those results to the ones obtained with the independent univariate analyses.

5.1 Data analysis using a bivariate Weibull as the underlying distribution. Table 3 presents the parameter estimates of the underlying bivariate Weibull distribution for

the attributes “odor” and “flavor”. The table shows also the estimated value of the correlation between the shelf life times (or “failure times’) of the two attributes and p-values associated to some hypothesis testing.

Table 3 - Parameter estimates of the bivariate Weibull model (attributes “odor “ and “flavor” – storage condition: “environmental chamber 2”) Parameter Estimates (1) p values

1α 1δ 2α 2δ λ correlation0 : 1 2H α α= 0 1 2:H δ δ= 0 : 0H λ =

0.068(0.0086) (2)

1.380(0.0109)

0.067(0,0079)

2.550(0.0156)

0.453 0.684 0.04 0.0007 0.002

(1) subscript “1” = flavor; “2”= odor; (2) standard deviation

After having examined the results contained in Table 3 we can say that: (1) there is a (statistically) significant (median to high) correlation between the two shelf life times associated to the two attributes studied. The estimated value, 0.68, is high enough to make the bivariate analysis worth it (p value for 0 : 0H λ = is 0.002) ; (2) we can conclude at 5% level that the failure times of the two attributes do not have the same marginal Weibull distributions, although the evidence leading to the rejection of the hypothesis involving the α parameters is weak (p value or

0 : 1 2H α α= is 0.04 and for 0 1 2:H δ δ= p value = 0.0007). Figure 2 (a) shows the point


estimates 00( )ˆˆ1 1 ( , )i i ip R τ τ− = − (fraction of units unfit for consumption) and pointwise two sided

95% approximate confidence intervals for 1 ( , )i iR τ τ− .

5.2 Comparison with the data analysis for each attribute independently, using univariate Weibulls as underlying distributions. We assume here that 0λ = (no correlation). These results agree with the ones obtained by Freitas, Borges and Ho (2003). Figure 2 (b) presents point estimates of the “fraction of defective units” obtained with the bivariate Weibull and by combing the results of the univariate independent fits. In addition, Table 5 shows the specific estimated values for weeks 1, 10 and 35 for each one of the two approaches discussed in this paper. We used equation (13) to combine the results of the two univariate analysis. The estimated values provided by the bivariate Weibull model are smaller then the ones based on the univariate fittings. According to the results of the simulation study, the values provided by the bivariate fitting should be closer to the real values then the ones calculated with the univariate independent Weibull fittings. From Table 4 we can say that for example, that by the first week of storage (in conditions similar to those of the “environmental chamber 2”) about 2,4% of the units will be inappropriate for consumption regarding at least one of the attributes considered. Therefore, the food company should use the estimated values provided by the bivariate model as a reference to provide the shopper with “expiration date” information in the food package or open dating like information, for instance, “do not use after Aug.25”. Table 4 – Estimated Fraction of unfit units obtained by the two approaches (weeks 1, 10 and 35).

week Bivariate Weibull Univariate Weibulls1 0.024263 0.1711

10 0.430828 0.846835 0.865693 0.9658

(a) (b) Figure 2 – (a) Point estimates of the fraction of unfit units (bivariate Weibull model) with pointwise two-sided 95% confidence intervals ; (b) point estimates of the “fraction of defective units” obtained with the bivariate Weibull and by combing the results of the univariate independent fits.

6. Concluding remarks and future workIn this paper we presented the preliminary results obtained with a model to analyze data

generated by sensory evaluations of food products. The main purpose was to use the data analysis to provide the producers with information regarding the shelf life of sensory characteristics of the product. The shelf life distribution is the key information behind the “expiration date” provided in the food´s package. The model proposed, uses dichotomization of the results and a bivariate Weibull as the underlying distribution. With this model, we were able to study two attributes jointly, taking into account the correlation between the “failure times” associated to each one of them. The simulation results indicated that the bivariate model provides better estimates for the fraction of unfit units than the ones provided by separate univariate fits (as far as the measures bias and standard deviation are concerned), specially in the cases where there is a strong


correlation between the “failure times” of the attributes. There remain many open questions about how to collect and analyze this kind of data. Some of these include the following: (1) we used a bivariate Weibull as the underlying distribution. It would be interesting to investigate other families of bivariate distributions; (2) we realize that with the dichotomization proposed (variable W), one looses the information of which one of the attributes “fails” first. It would be interesting then to work with the distribution of the minimum of 1 2{ ; }T T in order to assess this information; (3) the simulation study needs to be expanded in order to accomplish other practical situations; (4) there are important questions about how to design an experiment to provide the most efficient use of one´s resources and to assess the precision that one can expect to achieve with a specific design. Important questions include how to choose the interval between inspections, the number of units that should be tested and (5) sensitivity to possible model departures if of concern to all who analyze data. Sensitivity can be assessed to some extent, by changing assumptions and reanalyzing. It can be studied via simulation.

AcknowledgementsThe authors are grateful to CNPq (no.471217/2007-7) and Fapemig for the financial support of this research.7. References.Casella, G., Berger, R.L. Statistical Inference. 2nd. Edition. Pacific Grove, Ca, Duxbury, 2001.Collet, D. Modelling Survival Data in Medical Research, 2nd edition Florida: Chapman & Hall/CRC Press LLC, 2003.Finkelstein, D.M.(1986). “A proportional hazards model for interval-censored failure time data”. Biometrics , no. 42, pp. 845-854.Efron, B., Tibshirani, R. (1986). Bootstrap methods for standard errors, confidence intervals and other measures of statistical accuracy. Statistical Science, 1, 54-77.Freitas, M.A., Borges, W, and Ho, L.L. (2003). A Statistical Model for shelf life estimation using sensory evaluation scores. Communications in Statistics – Theory and Methods, no. 32, 8, pp. 1559-1589.Freitas. M.A., Costa, J.C. (2006) . Shelf Life determination using sensory evaluation scores: a general Weibull modeling approach. Computers and Industrial Engineering, 51, pp. 652-670. Freitas, M.A., Gomes, R. C.D. A proportional hazards model for interval-censored data generated by sensory evaluations: application to shelf life determination of food products, Annals of the International Conference on Degradation , Damage, Fatigue and Accelerated Life Models in Reliability Testing (ALT 2006), pp. 232-240, Angers, France, May 22-24, 2006.Gacula, M.C. Jr and Kubala, J.J. (1975). Statistical Models for Shelf Life failures. Journal of Food Sciences, no. 40, pp. 404-409. Gacula M.C.Jr. The design of experiments for shelf life study. (1975). Journal of Food Science, no. 40, pp 399-403.Gumbel, E.J. (1960). Bivariate exponential distribution. Journal of the American Statistical Association, 55, 698-707.Hougaard, P. (1986). A class of multivariate failure time distributions. Biometrika, 73, 671-678.Johnson, R.A., Evans, W.J., Green, D.W. (1999). Some Bivariate Distributions for Modeling the Strength Properties of Lumber. United States Department of Agriculture, Forest Service, Forest Products laboratory, Research Paper FLP-RP-575. Kalbfleisch, J.K., Prentice, R.L. The Statistical Analysis of Failure Time Data, New York: Wiley, 2002. Kim, M.Y., De Gruttola, V.G., Lagakos, S.W. (1993). Analyzing doubly censored data with covariates, with application of AIDS, Biometrics, no. 49, pp. 13-22. Lawless, J.F. Statistical Models and Methods for Lifetime Data, New York:: Wiley, 1982.Lu, J.., Bhattacharyya, G.K. (1990). Some new constructions of bivariate Weibull models. Annals of the Institute of Statistical Mathematics, 42(3), 543-559.Meeker. W.Q., Escobar, L.A. Statistical Methods for Reliability Data, Wiley, New York, 1998.Peto, R. (1973). “Experimental survival curves for interval censored data”, App. Statistics, no.22, pp. 86-91.

Rabinowitz, D., Tsiatis, A., Aragon, J. (1995). Regression with interval censored data, Biometrika., no. 82, 3, pp. 501-513. Turnbull, B. (1976). The empirical distribution function with arbitrarily grouped censored and truncated data”. Journal of the Royal Statistical Society., no. 38, pp. 290-295. Younes, N., Lachin, J. (1997). Link-based models for survival data with interval and continuous time censoring, Biometrics, no. 53, pp. 1199-1211, 1997.


Documents

A GENERAL BIVARIATE WEIBULL MODEL APPLIED TO SHELF …