- Email: [email protected]
Bayesian survival analysis modeling applied to sensory shelf life of foods
Food Quality and Preference 17 (2006) 307–312 www.elsevier.com/locate/foodqual
Bayesian survival analysis modeling applied to sensory shelf life of foods M. Luz Calle a, Guillermo Hough
b,*,1
, Ana Curia
b,1
, Guadalupe Go´mez
c
a
Departament d’Informatica i Matema´tica, Escola Politecnica Superior, Universitat de Vic, 08570 Vic, Spain Instituto Superior Experimental de Tecnologı´a Alimentaria, Nueve de Julio, 6500 Buenos Aires, Argentina Departament d’Estadı´stica i Investigacio´ Operativa, Universitat Polite`cnica de Catalunya, Pau Gargallo 5, 08028 Barcelona, Spain b
c
Received 25 September 2004; received in revised form 7 March 2005; accepted 22 March 2005 Available online 13 May 2005
Abstract Data from sensory shelf-life studies can be analyzed using survival statistical methods. The objective of this research was to introduce Bayesian methodology to sensory shelf-life studies and discuss its advantages in relation to classical (frequentist) methods. A specific algorithm which incorporates the interval censored data from shelf-life studies is presented. Calculations were applied to whole-fat and fat-free yogurt, each tasted by 80 consumers who answered ‘‘yes’’ or ‘‘no’’ to whether they would consume samples with different storage times. The Bayesian approach was used to obtain the failure functions and the estimated shelf lives of the yogurts. Two approaches were considered: first non-informative priors were used, and then results from the fat-free yogurt were used to construct a more informative prior for the whole-fat product. For both approaches a Weibull distribution was assumed. For the non-informative priors the results were analogous to those obtained by the classical methods. When an informative prior was used to estimate the whole-fat shelf-life parameters, there was a tendency to closer credible limits in the percentile distributions. Advantages of the Bayesian approach are discussed: no need for asymptotical assumptions and the possibility of improving the precision of the results by introducing external information. Ó 2005 Elsevier Ltd. All rights reserved. Keywords: Shelf life; Survival analysis statistics; Bayesian model; Yogurt
1. Introduction Sensory evaluation is the key factor for determining the shelf life of many food products. Survival analysis (Go´mez, 2004; Klein & Moeschberger, 1997; Kleinbaum, 1996; Meeker & Escobar, 1998) is a branch of statistics used extensively in clinical studies, epidemiology, biology, sociology, and reliability studies. Hough, Langohr, Go´mez, and Curia (2003) applied survival analysis statistics to sensory shelf life of foods. Hough, *
Corresponding author. Fax: +54 2317 431309. E-mail address: [email protected] (G. Hough). 1 Authors Hough and Curia are research fellows of the Comisio´n de Investigaciones Cientı´ficas de la Provincia de Buenos Aires. 0950-3293/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.foodqual.2005.03.012
Garitta, and Sa´nchez (2004) used this methodology to calculate acceptance limits to sensory defects which can be applied to setting quality control specifications. In survival analysis, the failure function F(t) can be defined as the probability of an individual failing before time t (Klein & Moeschberger, 1997). Referring this definition to the sensory shelf life of the food, the ‘‘individual’’ would not be the food itself, but rather the consumer, that is the failure function would be defined as the probability of a consumer rejecting a product stored for a time lower than t. The hazard would not be focused on the product deteriorating, rather on the consumer rejecting the product. Methods of survival analysis have been developed to evaluate times until an event of interest, often called
308
M.L. Calle et al. / Food Quality and Preference 17 (2006) 307–312
survival times, taking into account the presence of censored data. These occur whenever a time of interest, for example, the time until death in epidemiological studies or the time until failure of machines in reliability studies, cannot be observed exactly. If we know this time to be superior to an observed time, we speak of rightcensored data, if it is known to be less than the observation time, the time is called left-censored. Interval-censoring is given when the time of interest falls into an observed interval (Go´mez, Calle, & Oller, 2004). This kind of censoring appears in shelf-life studies where samples with different storage times are presented to consumers. For example, suppose that the storage times are: 0, 5, 10 and 15 days. If a consumer accepts the sample with 5 days storage and rejects it with 10 days storage, the exact rejection time could be any value between 5 and 10. This is defined as interval censoring. A special type of interval censoring is when a consumer rejects the sample with 5 days storage, thus rejection time is 65 and this is called left censoring. If the consumer accepts all samples, rejection would occur for a storage time >15 and the data is right censored. Hough et al. (2003) defined parametric models to estimate the failure function in shelf-life studies. Concepts and calculations were applied to a data set obtained from 50 consumers who each tasted 7 yogurt samples with different storage times, answering yes or no to whether they would consume the samples. From this censored data set parametric models were obtained which allowed shelf-life estimations. Hough et al. (2003) used frequentist or classical methods, that is methods based on maximum likelihood estimations, performing calculations with procedures from the S-PLUS statistical software (Insightful Corp., Seattle, Washington, USA). An alternative approach to sensory shelf-life calculations are the Bayesian methods. These improve on the traditional, frequentist approach to statistical analysis in their ability to formally incorporate the prior opinion of one or more experimenters into the results via the prior distribution (Carlin & Louis, 1998, chap. 4). These methods have also proven to be very appropriate to combine information coming from different similar experiments. When no additional external information is available, Bayesian and classical approaches usually result in similar conclusions. Bayesian methods have been used in a number of survival analysis applications in other fields. For example in HIV epidemiological studies the time of HIV infection was estimated as a function of the duration of intravenous drug use (Go´mez, Calle, Egea, & Muga, 2000). These Bayesian methods have not been applied to sensory shelf-life calculations and the objective of this paper is to introduce this approach and discuss its advantages in relation to classical methods.
2. Bayesian modeling for survival data 2.1. Bayesian inference Bayesian inference consists of the combination of two elements: data and prior information. This combination is performed in probabilistic terms through the Bayesian theorem which is summarized as follows (Gamerman, 1997, chap. 2). Data x is extracted from a population that is distributed according to a density f(xjh), where h are unknown parameters (to be estimated) and f(xjh) is called the data model or observational distribution. Looking at f(xjh) as a function of h gives the likelihood function. Prior information on the parameter of interest h is expressed through a probability distribution p(h) which is called the prior density. The prior density expresses the uncertainty on h before the observation of the value of x. Inferences will be based on the probability distribution of h after observing the value of x. This distribution, denoted by p(hjx), is called posterior distribution and is obtained by Bayes theorem f ðxjhÞpðhÞ pðhjxÞ ¼ ; f ðxÞ where Z f ðxÞ ¼ f ðxjhÞpðhÞ dh. 2.2. Bayesian analysis of interval-censored data This general Bayesian approach cannot be directly applied when analyzing interval-censored survival data. The reason is that due to censoring, the exact survival times T1, . . . , Tn of each individual are not observed, rather only an interval [Li, Ri] is observed where the event has occurred. The way to deal with censored data is to treat the censored times T1, . . . , Tn as further unknowns. The advantage of this data augmentation method is that it can be easily implemented with the WinBugs program (Bayesian inference Using Gibbs Sampling) (Spiegelhalter et al., 1996). This program is specifically developed to deal with Bayesian models, and it provides many useful numerical and graphical aspects of the posterior distribution which facilitates the interpretation of the results. For the implementation it is necessary to specify both the data model and the prior distribution, from which the posterior distribution is calculated. 2.2.1. Data model The Weibull distribution is a very flexible right skewed distribution which is particularly appropriate for modeling survival data and it has been used previously in food shelf-life modeling (Cardelli & Labuza, 2001; Duyvesteyn, Shimoni, & Labuza, 2001; Hough, Puglieso, Sa´nchez, & Mendes da Silva, 1999). Thus a
M.L. Calle et al. / Food Quality and Preference 17 (2006) 307–312
Weibull distribution was chosen for the data model. The failure function of a Weibull distribution is given by lnðtÞ l F ðtÞ ¼ F sev r where Fsev(Æ) is the distribution function of the smallest extreme value distribution, Fsev(w) = 1 exp(ew ), and l (location parameter) and r (shape parameter) are the models parameters of interest. 2.2.2. Prior distributions Prior distributions for l and r are needed. In this work we will consider non-informative priors that do not favor any possible parameter value over any other. This is the standard way of expressing the lack of expert knowledge or that we dont want to use it in the analysis. Such an approach will allow a better comparison between the results obtained with the proposed Bayesian methodology and the results obtained with the classical approach (Hough et al., 2003). Usual non-informative priors for l were considered, specified in two stages in the following hierarchical form:
309
truncated because Ti is restricted to the Li, Ri interval. Step 2. Given the imputed times T1, . . . , Tn, sample l and r from their corresponding posterior distribution based on a complete observed sample: p(ljT1, . . . , Tn) and p(rjT1, . . . , Tn). The successive implementation of these two steps yields a sample of the parameters of interest l and r and any posterior estimate, such as the mean, median or percentile, can be obtained from this sample. This algorithm can be performed with the program WinBUGS (Bayesian inference Using Gibbs Sampling) (Spiegelhalter et al., 1996) which is a program specifically developed to deal with complex Bayesian models. The WinBUGS program can be obtained freely at the following site: http://www.mrc-bsu.cam.ac.uk/bugs/, and the code of the algorithm is available from the authors upon request.
3. Experimental data
Prior for l : pðlÞ N ða0 ; r0 Þ
3.1. Samples
Prior for a0 : pða0 Þ N ð0; 103 Þ
Commercial stirred, strawberry-flavored yogurt with strawberry pulp was used; both whole-fat (2.7% fat) and fat-free (0.05% fat) products were tested. Bottles (1000 ml) from different batches were stored at 10 °C in such a way as to have samples with different storage times ready on the same day. Storage times were chosen based on a preliminary experiment and were 0, 14, 28, 42, 56, 70 and 84 days. All batches were made with the same formulation and were checked to be similar to the previous batch by consensus among three expert assessors. Once all samples had reached their storage time at 10 °C, they were refrigerated at 2 °C, until they were tasted. All the measurements were made in a period not longer than a week, to guarantee no changes in the samples. Microbiological analysis (aerobic mesophiles, coliforms, yeasts, and molds) showed that the samples were fit for consumption. The Ethical Committee of our Institute decided that all samples were adequate for tests on humans in the quantities to be served.
Prior for r0 : pðr0 Þ IGð0.001; 0.001Þ where a0 and r0 are called hyperparameters and are given diffuse priors. For r the usual non-informative prior is the inverse gamma distribution Prior for r : pðrÞ IGð0.01; 0.01Þ More informative prior distributions can be constructed as explained below in the results section, where results from fat-free yogurts will be used to build a more informative prior for the whole-fat estimations. 2.2.3. Posterior distribution with censored data As T1, . . . , Tn are censored values treated as further unknowns, the vector h of parameters is h = (T1, . . . , Tn, l, r). The posterior distribution of h given the observations, [L1, R1] , . . . , [Ln, Rn], is in general not available analytically. We will use the iterative simulation method proposed by Calle (2003) to obtain a sample from the posterior distribution p(hj[L1, R1] , . . . , [Ln, Rn]). The algorithm consists in giving initial values for the unknown parameters and then iteratively perform steps 1 and 2: Step 0. Give initial values for l and r. Step 1. Given the current values of l and r, sample each censored time Ti from the corresponding posterior distribution, p(Tijl, r, [Li, Ri]), which in our setting is a truncated Weibull distribution. It is
3.2. Consumers Consumers who consumed either whole-fat or fatfree strawberry yogurt at least once a week were recruited from the population of the city of Nueve de Julio, Buenos Aires, Argentina. For each type of yogurt 80 consumers were recruited, that is a total of 160 consumers participated in the study. Each consumer received the 7 yogurt samples (corresponding to each storage time at 10 °C) monadically in random order. Fifty grams of each sample were served in 70-ml plastic
310
M.L. Calle et al. / Food Quality and Preference 17 (2006) 307–312
cups. Water was available for rinsing. For each sample they had to taste it and answer the question: ‘‘Would you normally consume this product? Yes or No?’’. It was explained that this meant that if they bought the product to drink it, or it was served to them at their homes, whether they would consume it or not. The tests were conducted in a sensory laboratory with individual booths with artificial daylight type illumination, temperature control (between 22 and 24 °C) and air circulation.
4. Results and discussion 4.1. Raw data and censoring considerations Table 1 presents data for 5 of the 80 whole-fat yogurt consumers to illustrate the interpretation given to each consumers data. For a discussion on this interpretation see Hough et al. (2003). For the 80 whole-fat consumers, 42 were interval censored, 14 right censored, 4 left censored and 20 rejected the fresh sample. For the 80 fat-free consumers, 53 were interval censored, 2 right censored, 18 left censored and 7 rejected the fresh sample. 4.2. Failure probability calculations Weibull model parameters and percentiles for wholefat and fat-free yogurts are in Table 2. Fig. 1 shows the failure function expressed as percent of consumers rejecting the products versus storage time. The parameters of interest in sensory shelf-life studies are usually the percentiles of the failure function distribution, thus these are presented in Table 2. For example, it is of interest to know how many days the yogurts can be stored so that less than 50% of the consumers reject the
product; or what is the difference in storage time between whole-fat and fat-free yogurts for 25% of the consumers rejecting the product at the end of its shelf life. In principle it seems out of the question to have 25–50% of consumers rejecting the product, but in actual fact it is 25–50% of the small proportion of consumers who taste the product at the end of its shelf life. Fig. 1 shows that fat-free yogurt had a shorter shelf life than whole-fat yogurt, with differences increasing as the percent consumers rejecting the product increased. From Table 2 shelf lives for whole-fat and fat-free yogurts for 25% consumer rejection, were 41 days and 36 days, respectively. It could be of interest to know whether this difference was significant. The P25 difference distribution is easily estimated by Bayesian methodology by sampling the difference as a new parameter together with l and r in the corresponding algorithm. These differences with their 95% credible limits are in Table 3 for P25 and P50. For the P25 difference the limits included zero, thus we accepted the null hypothesis of no difference between the shelf lives of whole-fat and fat-free yogurts corresponding to 25% consumer rejection. For the P50 difference the limits did not include zero, thus we rejected the null hypothesis and concluded that shelf lives between whole-fat and fat-free yogurts corresponding to 50% consumer rejection were different. Classical calculations as outlined by Hough et al. (2003) were also applied to the data. Results were practically identical. For example, l for whole-fat yogurt was 4.25 with lower and upper 95% confidence intervals of 4.13 and 4.38, respectively. Due to the non-informative priors used in the Bayesian approach, this coincidence in results was to be expected and it is a way of validating the Bayesian methodology.
Table 1 Acceptance/rejection data for 5 subjects who tasted whole-fat strawberry yogurt samples with different storage times at 10 °C Subject
1 2 3 4 5
Storage time (days)
Censoring
0
14
28
42
56
70
84
Yes Yes Yes No Yes
Yes No Yes Yes Yes
Yes Yes No Yes Yes
No Yes Yes Yes Yes
No No No Yes Yes
No No No No Yes
No No No Yes Yes
Interval: 28–42 Left: 656 Interval: 14–56 Not consider Right: >84
Table 2 Estimates and 95% credible intervals of the Weibull model parameters (l and r) and percentiles (P) for whole-fat and fat-free yogurts Parameter
l r P25 P50
Whole-fat strawberry
Fat-free strawberry
95% Lower limit
Estimated parameter
95% Upper limit
95% Lower limit
Estimated parameter
95% Upper limit
4.13 0.33 32.2 52.1
4.26 0.44 41.1 60.3
4.4 0.59 49.4 69.1
3.89 0.25 30.4 42.7
3.98 0.32 36.2 47.7
4.07 0.41 41.4 52.5
M.L. Calle et al. / Food Quality and Preference 17 (2006) 307–312
100
% of consumers rejecting
Fat-free Whole-fat
75
50
25
0 0
14
28
42 56 Storage time (days)
70
84
311
Since the posterior distribution of r2 corresponding to the fat-free yogurt was approximately an inverse gamma IG(17.12, 1.67), an inverse gamma distribution for the whole-fat yogurt was also assumed, that is IG(a1, b1), with b1 = 1.67 and a1 varying around 17.12, p(a1) G(13, 0.76). Results from using non-informative and informative priors for whole-fat yogurt are in Table 4. Results are similar with a tendency to obtain closer credible limits in the percentile distributions. Having used only one set of data to improve the priori of the whole-fat yogurt led to small improvements in the posterior distributions. However, it is to be expected that the use of this methodology based on large data sets, as can be constructed in industry over the years, will allow specifying more restrictive priors which would lead to more precise shelf-life estimations.
Fig. 1. Percent of consumers rejecting yogurts as a function of storage time for the Weibull model.
5. Conclusions Table 3 Percentile (P) differences between fat-free and whole-fat yogurts for 25% and 50% consumer rejection, with their corresponding 95% credible limits
P25 P50
95% Lower limit
Fat-free–whole-fat
95% Upper limit
5 3
5 13
15 23
As mentioned above, an advantage of Bayesian methods is the possibility of introducing prior information to improve the precision of the estimates of the shelf-life parameters. To illustrate this possibility the results obtained from the fat-free yogurt were used to construct the hyperpriors as prior information in calculating the whole-fat yogurt posterior distributions. The posterior distribution of l for the fat-free yogurt was approximately a normal distribution N(3.98, 0.044). Thus, the prior distribution for l for the whole-fat yogurt was also considered normal, that is N(a0, r0), with r0 = 0.044 which is the standard deviation of the posterior fat-free yogurt l distribution. The parameter a0 was allowed to vary around the corresponding estimated parameter of the fat-free yogurt, 3.98, by setting the hyperprior: p(a0) N(3.98, 0.25).
Lately Bayesian methods have been used in many research studies, especially in the field of medicine, as an alternative to classical or frequentist statistical methods. One of the reasons is that classical methods base their maximum likelihood estimations on asymptotic considerations that are usually only valid for a considerable data size. In practice the methods are applied to relatively small data sets where the validity of the asymptotic assumptions is doubtful. For example, in calculating the confidence intervals for the 50% percentile, a normal distribution of the percentile is assumed. This assumption may hold but may not necessarily be so. In the Bayesian approach no assumption is made as to the shape of the percentile distribution, rather the data themselves specify the distribution. Another advantage of the Bayesian approach is the possibility of improving the precision of the results by introducing external information in terms of the priori distribution. This can be useful for R&D where most developments are modifications of existing products, and thus previous data bases can used to better adjust calculations of new products. A Bayesian approach has been applied to data from sensory shelf life of foods, with results comparable to
Table 4 Estimates and 95% credible intervals of the Weibull model parameters (l and r) and percentiles (P) for whole-fat yogurt using non-informative and informative priors Parameter
l r P25 P50
Whole-fat strawberry non-informative
Whole-fat strawberry informative
95% Lower limit
Estimated parameter
95% Upper limit
95% Lower limit
Estimated parameter
95% Upper limit
4.13 0.33 32.2 52.1
4.26 0.44 41.1 60.3
4.40 0.59 49.4 69.1
4.14 0.30 35.2 53.4
4.25 0.39 43.5 61.1
4.36 0.51 51.1 68.1
312
M.L. Calle et al. / Food Quality and Preference 17 (2006) 307–312
the classical frequentist method. Future research in the application of Bayesian methods points to the introduction of external information to improve calculations for new products, and also to use these methods to combine information from different populations, for example data generated in different countries.
Acknowledgement Agencia Nacional de Promocio´n Cientı´fica y Tecnolo´gica PICT 98-09-04827 and the Universitat de Vic for a travel grant. References Calle, M. L. (2003). Parametric Bayesian analysis of interval-censored and doubly-censored data. Journal of Probability and Statistical Sciences, 1, 103–118. Cardelli, C., & Labuza, T. P. (2001). Application of Weibull hazard analysis to the determination of the shelf life of roasted and ground coffee. Lebensmittel-Wissenschaft und-Technologie, 34, 273–278. Carlin, B. P., & Louis, T. A. (1998). Bayes and empirical Bayes methods for data analysis. Boca Raton: CRC Press. Duyvesteyn, W. S., Shimoni, E., & Labuza, T. P. (2001). Determination of the end of shelf-life for milk using Weibull
hazard method. Lebensmittel-Wissenschaft und-Technologie, 34, 143–148. Gamerman, D. (1997). Markov Chain Monte Carlo. Stochastic simulation for Bayesian inference. Boca Raton: CRC Press. Go´mez, G. (2004). Ana´lisis de Supervivencia. 84-688-5607-X. Universitat Polite`cnica de Catalunya. Go´mez, G., Calle, M. L., Egea, J. M., & Muga, R. (2000). Estimation of the risk of HIV infection as a function of the length of intravenous drug use. A nonparametric Bayesian approach. Statistics in Medicine, 19, 2641–2656. Go´mez, G., Calle, M. L., & Oller, R. (2004). Frequentist and Bayesian approaches for interval-censored data. Statistical Papers, 45, 139–173. Hough, G., Garitta, L., & Sa´nchez, R. (2004). Determination of consumer acceptance limits to sensory defects using survival analysis. Food Quality and Preference, 15, 729–734. Hough, G., Langohr, K., Go´mez, G., & Curia, A. (2003). Survival analysis applied to sensory shelf life of foods. Journal of Food Science, 68, 359–362. Hough, G., Puglieso, M. L., Sa´nchez, R., & Mendes da Silva, O. (1999). Sensory and microbiological shelf-life of a commercial Ricotta cheese. Journal of Dairy Science, 82, 454–459. Klein, J. P., & Moeschberger, M. L. (1997). Survival analysis, techniques for censored and truncated data (502p). New York: Springer-Verlag. Kleinbaum, D. G. (1996). Survival analysis, a self-learning text (324p). New York: Springer-Verlag. Meeker, W. Q., & Escobar, L. A. (1998). Statistical methods for reliability data (680p). New York: John Wiley & Sons. Spiegelhalter, D. et al. (1996). Bayesian inference using Gibbs sampling, Version 0.5, (version ii). Cambridge: MRC Biostatistics Unit.
Bayesian survival analysis modeling applied to sensory shelf life of foods M. Luz Calle a, Guillermo Hough
b,*,1
, Ana Curia
b,1
, Guadalupe Go´mez
c
a
Departament d’Informatica i Matema´tica, Escola Politecnica Superior, Universitat de Vic, 08570 Vic, Spain Instituto Superior Experimental de Tecnologı´a Alimentaria, Nueve de Julio, 6500 Buenos Aires, Argentina Departament d’Estadı´stica i Investigacio´ Operativa, Universitat Polite`cnica de Catalunya, Pau Gargallo 5, 08028 Barcelona, Spain b
c
Received 25 September 2004; received in revised form 7 March 2005; accepted 22 March 2005 Available online 13 May 2005
Abstract Data from sensory shelf-life studies can be analyzed using survival statistical methods. The objective of this research was to introduce Bayesian methodology to sensory shelf-life studies and discuss its advantages in relation to classical (frequentist) methods. A specific algorithm which incorporates the interval censored data from shelf-life studies is presented. Calculations were applied to whole-fat and fat-free yogurt, each tasted by 80 consumers who answered ‘‘yes’’ or ‘‘no’’ to whether they would consume samples with different storage times. The Bayesian approach was used to obtain the failure functions and the estimated shelf lives of the yogurts. Two approaches were considered: first non-informative priors were used, and then results from the fat-free yogurt were used to construct a more informative prior for the whole-fat product. For both approaches a Weibull distribution was assumed. For the non-informative priors the results were analogous to those obtained by the classical methods. When an informative prior was used to estimate the whole-fat shelf-life parameters, there was a tendency to closer credible limits in the percentile distributions. Advantages of the Bayesian approach are discussed: no need for asymptotical assumptions and the possibility of improving the precision of the results by introducing external information. Ó 2005 Elsevier Ltd. All rights reserved. Keywords: Shelf life; Survival analysis statistics; Bayesian model; Yogurt
1. Introduction Sensory evaluation is the key factor for determining the shelf life of many food products. Survival analysis (Go´mez, 2004; Klein & Moeschberger, 1997; Kleinbaum, 1996; Meeker & Escobar, 1998) is a branch of statistics used extensively in clinical studies, epidemiology, biology, sociology, and reliability studies. Hough, Langohr, Go´mez, and Curia (2003) applied survival analysis statistics to sensory shelf life of foods. Hough, *
Corresponding author. Fax: +54 2317 431309. E-mail address: [email protected] (G. Hough). 1 Authors Hough and Curia are research fellows of the Comisio´n de Investigaciones Cientı´ficas de la Provincia de Buenos Aires. 0950-3293/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.foodqual.2005.03.012
Garitta, and Sa´nchez (2004) used this methodology to calculate acceptance limits to sensory defects which can be applied to setting quality control specifications. In survival analysis, the failure function F(t) can be defined as the probability of an individual failing before time t (Klein & Moeschberger, 1997). Referring this definition to the sensory shelf life of the food, the ‘‘individual’’ would not be the food itself, but rather the consumer, that is the failure function would be defined as the probability of a consumer rejecting a product stored for a time lower than t. The hazard would not be focused on the product deteriorating, rather on the consumer rejecting the product. Methods of survival analysis have been developed to evaluate times until an event of interest, often called
308
M.L. Calle et al. / Food Quality and Preference 17 (2006) 307–312
survival times, taking into account the presence of censored data. These occur whenever a time of interest, for example, the time until death in epidemiological studies or the time until failure of machines in reliability studies, cannot be observed exactly. If we know this time to be superior to an observed time, we speak of rightcensored data, if it is known to be less than the observation time, the time is called left-censored. Interval-censoring is given when the time of interest falls into an observed interval (Go´mez, Calle, & Oller, 2004). This kind of censoring appears in shelf-life studies where samples with different storage times are presented to consumers. For example, suppose that the storage times are: 0, 5, 10 and 15 days. If a consumer accepts the sample with 5 days storage and rejects it with 10 days storage, the exact rejection time could be any value between 5 and 10. This is defined as interval censoring. A special type of interval censoring is when a consumer rejects the sample with 5 days storage, thus rejection time is 65 and this is called left censoring. If the consumer accepts all samples, rejection would occur for a storage time >15 and the data is right censored. Hough et al. (2003) defined parametric models to estimate the failure function in shelf-life studies. Concepts and calculations were applied to a data set obtained from 50 consumers who each tasted 7 yogurt samples with different storage times, answering yes or no to whether they would consume the samples. From this censored data set parametric models were obtained which allowed shelf-life estimations. Hough et al. (2003) used frequentist or classical methods, that is methods based on maximum likelihood estimations, performing calculations with procedures from the S-PLUS statistical software (Insightful Corp., Seattle, Washington, USA). An alternative approach to sensory shelf-life calculations are the Bayesian methods. These improve on the traditional, frequentist approach to statistical analysis in their ability to formally incorporate the prior opinion of one or more experimenters into the results via the prior distribution (Carlin & Louis, 1998, chap. 4). These methods have also proven to be very appropriate to combine information coming from different similar experiments. When no additional external information is available, Bayesian and classical approaches usually result in similar conclusions. Bayesian methods have been used in a number of survival analysis applications in other fields. For example in HIV epidemiological studies the time of HIV infection was estimated as a function of the duration of intravenous drug use (Go´mez, Calle, Egea, & Muga, 2000). These Bayesian methods have not been applied to sensory shelf-life calculations and the objective of this paper is to introduce this approach and discuss its advantages in relation to classical methods.
2. Bayesian modeling for survival data 2.1. Bayesian inference Bayesian inference consists of the combination of two elements: data and prior information. This combination is performed in probabilistic terms through the Bayesian theorem which is summarized as follows (Gamerman, 1997, chap. 2). Data x is extracted from a population that is distributed according to a density f(xjh), where h are unknown parameters (to be estimated) and f(xjh) is called the data model or observational distribution. Looking at f(xjh) as a function of h gives the likelihood function. Prior information on the parameter of interest h is expressed through a probability distribution p(h) which is called the prior density. The prior density expresses the uncertainty on h before the observation of the value of x. Inferences will be based on the probability distribution of h after observing the value of x. This distribution, denoted by p(hjx), is called posterior distribution and is obtained by Bayes theorem f ðxjhÞpðhÞ pðhjxÞ ¼ ; f ðxÞ where Z f ðxÞ ¼ f ðxjhÞpðhÞ dh. 2.2. Bayesian analysis of interval-censored data This general Bayesian approach cannot be directly applied when analyzing interval-censored survival data. The reason is that due to censoring, the exact survival times T1, . . . , Tn of each individual are not observed, rather only an interval [Li, Ri] is observed where the event has occurred. The way to deal with censored data is to treat the censored times T1, . . . , Tn as further unknowns. The advantage of this data augmentation method is that it can be easily implemented with the WinBugs program (Bayesian inference Using Gibbs Sampling) (Spiegelhalter et al., 1996). This program is specifically developed to deal with Bayesian models, and it provides many useful numerical and graphical aspects of the posterior distribution which facilitates the interpretation of the results. For the implementation it is necessary to specify both the data model and the prior distribution, from which the posterior distribution is calculated. 2.2.1. Data model The Weibull distribution is a very flexible right skewed distribution which is particularly appropriate for modeling survival data and it has been used previously in food shelf-life modeling (Cardelli & Labuza, 2001; Duyvesteyn, Shimoni, & Labuza, 2001; Hough, Puglieso, Sa´nchez, & Mendes da Silva, 1999). Thus a
M.L. Calle et al. / Food Quality and Preference 17 (2006) 307–312
Weibull distribution was chosen for the data model. The failure function of a Weibull distribution is given by lnðtÞ l F ðtÞ ¼ F sev r where Fsev(Æ) is the distribution function of the smallest extreme value distribution, Fsev(w) = 1 exp(ew ), and l (location parameter) and r (shape parameter) are the models parameters of interest. 2.2.2. Prior distributions Prior distributions for l and r are needed. In this work we will consider non-informative priors that do not favor any possible parameter value over any other. This is the standard way of expressing the lack of expert knowledge or that we dont want to use it in the analysis. Such an approach will allow a better comparison between the results obtained with the proposed Bayesian methodology and the results obtained with the classical approach (Hough et al., 2003). Usual non-informative priors for l were considered, specified in two stages in the following hierarchical form:
309
truncated because Ti is restricted to the Li, Ri interval. Step 2. Given the imputed times T1, . . . , Tn, sample l and r from their corresponding posterior distribution based on a complete observed sample: p(ljT1, . . . , Tn) and p(rjT1, . . . , Tn). The successive implementation of these two steps yields a sample of the parameters of interest l and r and any posterior estimate, such as the mean, median or percentile, can be obtained from this sample. This algorithm can be performed with the program WinBUGS (Bayesian inference Using Gibbs Sampling) (Spiegelhalter et al., 1996) which is a program specifically developed to deal with complex Bayesian models. The WinBUGS program can be obtained freely at the following site: http://www.mrc-bsu.cam.ac.uk/bugs/, and the code of the algorithm is available from the authors upon request.
3. Experimental data
Prior for l : pðlÞ N ða0 ; r0 Þ
3.1. Samples
Prior for a0 : pða0 Þ N ð0; 103 Þ
Commercial stirred, strawberry-flavored yogurt with strawberry pulp was used; both whole-fat (2.7% fat) and fat-free (0.05% fat) products were tested. Bottles (1000 ml) from different batches were stored at 10 °C in such a way as to have samples with different storage times ready on the same day. Storage times were chosen based on a preliminary experiment and were 0, 14, 28, 42, 56, 70 and 84 days. All batches were made with the same formulation and were checked to be similar to the previous batch by consensus among three expert assessors. Once all samples had reached their storage time at 10 °C, they were refrigerated at 2 °C, until they were tasted. All the measurements were made in a period not longer than a week, to guarantee no changes in the samples. Microbiological analysis (aerobic mesophiles, coliforms, yeasts, and molds) showed that the samples were fit for consumption. The Ethical Committee of our Institute decided that all samples were adequate for tests on humans in the quantities to be served.
Prior for r0 : pðr0 Þ IGð0.001; 0.001Þ where a0 and r0 are called hyperparameters and are given diffuse priors. For r the usual non-informative prior is the inverse gamma distribution Prior for r : pðrÞ IGð0.01; 0.01Þ More informative prior distributions can be constructed as explained below in the results section, where results from fat-free yogurts will be used to build a more informative prior for the whole-fat estimations. 2.2.3. Posterior distribution with censored data As T1, . . . , Tn are censored values treated as further unknowns, the vector h of parameters is h = (T1, . . . , Tn, l, r). The posterior distribution of h given the observations, [L1, R1] , . . . , [Ln, Rn], is in general not available analytically. We will use the iterative simulation method proposed by Calle (2003) to obtain a sample from the posterior distribution p(hj[L1, R1] , . . . , [Ln, Rn]). The algorithm consists in giving initial values for the unknown parameters and then iteratively perform steps 1 and 2: Step 0. Give initial values for l and r. Step 1. Given the current values of l and r, sample each censored time Ti from the corresponding posterior distribution, p(Tijl, r, [Li, Ri]), which in our setting is a truncated Weibull distribution. It is
3.2. Consumers Consumers who consumed either whole-fat or fatfree strawberry yogurt at least once a week were recruited from the population of the city of Nueve de Julio, Buenos Aires, Argentina. For each type of yogurt 80 consumers were recruited, that is a total of 160 consumers participated in the study. Each consumer received the 7 yogurt samples (corresponding to each storage time at 10 °C) monadically in random order. Fifty grams of each sample were served in 70-ml plastic
310
M.L. Calle et al. / Food Quality and Preference 17 (2006) 307–312
cups. Water was available for rinsing. For each sample they had to taste it and answer the question: ‘‘Would you normally consume this product? Yes or No?’’. It was explained that this meant that if they bought the product to drink it, or it was served to them at their homes, whether they would consume it or not. The tests were conducted in a sensory laboratory with individual booths with artificial daylight type illumination, temperature control (between 22 and 24 °C) and air circulation.
4. Results and discussion 4.1. Raw data and censoring considerations Table 1 presents data for 5 of the 80 whole-fat yogurt consumers to illustrate the interpretation given to each consumers data. For a discussion on this interpretation see Hough et al. (2003). For the 80 whole-fat consumers, 42 were interval censored, 14 right censored, 4 left censored and 20 rejected the fresh sample. For the 80 fat-free consumers, 53 were interval censored, 2 right censored, 18 left censored and 7 rejected the fresh sample. 4.2. Failure probability calculations Weibull model parameters and percentiles for wholefat and fat-free yogurts are in Table 2. Fig. 1 shows the failure function expressed as percent of consumers rejecting the products versus storage time. The parameters of interest in sensory shelf-life studies are usually the percentiles of the failure function distribution, thus these are presented in Table 2. For example, it is of interest to know how many days the yogurts can be stored so that less than 50% of the consumers reject the
product; or what is the difference in storage time between whole-fat and fat-free yogurts for 25% of the consumers rejecting the product at the end of its shelf life. In principle it seems out of the question to have 25–50% of consumers rejecting the product, but in actual fact it is 25–50% of the small proportion of consumers who taste the product at the end of its shelf life. Fig. 1 shows that fat-free yogurt had a shorter shelf life than whole-fat yogurt, with differences increasing as the percent consumers rejecting the product increased. From Table 2 shelf lives for whole-fat and fat-free yogurts for 25% consumer rejection, were 41 days and 36 days, respectively. It could be of interest to know whether this difference was significant. The P25 difference distribution is easily estimated by Bayesian methodology by sampling the difference as a new parameter together with l and r in the corresponding algorithm. These differences with their 95% credible limits are in Table 3 for P25 and P50. For the P25 difference the limits included zero, thus we accepted the null hypothesis of no difference between the shelf lives of whole-fat and fat-free yogurts corresponding to 25% consumer rejection. For the P50 difference the limits did not include zero, thus we rejected the null hypothesis and concluded that shelf lives between whole-fat and fat-free yogurts corresponding to 50% consumer rejection were different. Classical calculations as outlined by Hough et al. (2003) were also applied to the data. Results were practically identical. For example, l for whole-fat yogurt was 4.25 with lower and upper 95% confidence intervals of 4.13 and 4.38, respectively. Due to the non-informative priors used in the Bayesian approach, this coincidence in results was to be expected and it is a way of validating the Bayesian methodology.
Table 1 Acceptance/rejection data for 5 subjects who tasted whole-fat strawberry yogurt samples with different storage times at 10 °C Subject
1 2 3 4 5
Storage time (days)
Censoring
0
14
28
42
56
70
84
Yes Yes Yes No Yes
Yes No Yes Yes Yes
Yes Yes No Yes Yes
No Yes Yes Yes Yes
No No No Yes Yes
No No No No Yes
No No No Yes Yes
Interval: 28–42 Left: 656 Interval: 14–56 Not consider Right: >84
Table 2 Estimates and 95% credible intervals of the Weibull model parameters (l and r) and percentiles (P) for whole-fat and fat-free yogurts Parameter
l r P25 P50
Whole-fat strawberry
Fat-free strawberry
95% Lower limit
Estimated parameter
95% Upper limit
95% Lower limit
Estimated parameter
95% Upper limit
4.13 0.33 32.2 52.1
4.26 0.44 41.1 60.3
4.4 0.59 49.4 69.1
3.89 0.25 30.4 42.7
3.98 0.32 36.2 47.7
4.07 0.41 41.4 52.5
M.L. Calle et al. / Food Quality and Preference 17 (2006) 307–312
100
% of consumers rejecting
Fat-free Whole-fat
75
50
25
0 0
14
28
42 56 Storage time (days)
70
84
311
Since the posterior distribution of r2 corresponding to the fat-free yogurt was approximately an inverse gamma IG(17.12, 1.67), an inverse gamma distribution for the whole-fat yogurt was also assumed, that is IG(a1, b1), with b1 = 1.67 and a1 varying around 17.12, p(a1) G(13, 0.76). Results from using non-informative and informative priors for whole-fat yogurt are in Table 4. Results are similar with a tendency to obtain closer credible limits in the percentile distributions. Having used only one set of data to improve the priori of the whole-fat yogurt led to small improvements in the posterior distributions. However, it is to be expected that the use of this methodology based on large data sets, as can be constructed in industry over the years, will allow specifying more restrictive priors which would lead to more precise shelf-life estimations.
Fig. 1. Percent of consumers rejecting yogurts as a function of storage time for the Weibull model.
5. Conclusions Table 3 Percentile (P) differences between fat-free and whole-fat yogurts for 25% and 50% consumer rejection, with their corresponding 95% credible limits
P25 P50
95% Lower limit
Fat-free–whole-fat
95% Upper limit
5 3
5 13
15 23
As mentioned above, an advantage of Bayesian methods is the possibility of introducing prior information to improve the precision of the estimates of the shelf-life parameters. To illustrate this possibility the results obtained from the fat-free yogurt were used to construct the hyperpriors as prior information in calculating the whole-fat yogurt posterior distributions. The posterior distribution of l for the fat-free yogurt was approximately a normal distribution N(3.98, 0.044). Thus, the prior distribution for l for the whole-fat yogurt was also considered normal, that is N(a0, r0), with r0 = 0.044 which is the standard deviation of the posterior fat-free yogurt l distribution. The parameter a0 was allowed to vary around the corresponding estimated parameter of the fat-free yogurt, 3.98, by setting the hyperprior: p(a0) N(3.98, 0.25).
Lately Bayesian methods have been used in many research studies, especially in the field of medicine, as an alternative to classical or frequentist statistical methods. One of the reasons is that classical methods base their maximum likelihood estimations on asymptotic considerations that are usually only valid for a considerable data size. In practice the methods are applied to relatively small data sets where the validity of the asymptotic assumptions is doubtful. For example, in calculating the confidence intervals for the 50% percentile, a normal distribution of the percentile is assumed. This assumption may hold but may not necessarily be so. In the Bayesian approach no assumption is made as to the shape of the percentile distribution, rather the data themselves specify the distribution. Another advantage of the Bayesian approach is the possibility of improving the precision of the results by introducing external information in terms of the priori distribution. This can be useful for R&D where most developments are modifications of existing products, and thus previous data bases can used to better adjust calculations of new products. A Bayesian approach has been applied to data from sensory shelf life of foods, with results comparable to
Table 4 Estimates and 95% credible intervals of the Weibull model parameters (l and r) and percentiles (P) for whole-fat yogurt using non-informative and informative priors Parameter
l r P25 P50
Whole-fat strawberry non-informative
Whole-fat strawberry informative
95% Lower limit
Estimated parameter
95% Upper limit
95% Lower limit
Estimated parameter
95% Upper limit
4.13 0.33 32.2 52.1
4.26 0.44 41.1 60.3
4.40 0.59 49.4 69.1
4.14 0.30 35.2 53.4
4.25 0.39 43.5 61.1
4.36 0.51 51.1 68.1
312
M.L. Calle et al. / Food Quality and Preference 17 (2006) 307–312
the classical frequentist method. Future research in the application of Bayesian methods points to the introduction of external information to improve calculations for new products, and also to use these methods to combine information from different populations, for example data generated in different countries.
Acknowledgement Agencia Nacional de Promocio´n Cientı´fica y Tecnolo´gica PICT 98-09-04827 and the Universitat de Vic for a travel grant. References Calle, M. L. (2003). Parametric Bayesian analysis of interval-censored and doubly-censored data. Journal of Probability and Statistical Sciences, 1, 103–118. Cardelli, C., & Labuza, T. P. (2001). Application of Weibull hazard analysis to the determination of the shelf life of roasted and ground coffee. Lebensmittel-Wissenschaft und-Technologie, 34, 273–278. Carlin, B. P., & Louis, T. A. (1998). Bayes and empirical Bayes methods for data analysis. Boca Raton: CRC Press. Duyvesteyn, W. S., Shimoni, E., & Labuza, T. P. (2001). Determination of the end of shelf-life for milk using Weibull
hazard method. Lebensmittel-Wissenschaft und-Technologie, 34, 143–148. Gamerman, D. (1997). Markov Chain Monte Carlo. Stochastic simulation for Bayesian inference. Boca Raton: CRC Press. Go´mez, G. (2004). Ana´lisis de Supervivencia. 84-688-5607-X. Universitat Polite`cnica de Catalunya. Go´mez, G., Calle, M. L., Egea, J. M., & Muga, R. (2000). Estimation of the risk of HIV infection as a function of the length of intravenous drug use. A nonparametric Bayesian approach. Statistics in Medicine, 19, 2641–2656. Go´mez, G., Calle, M. L., & Oller, R. (2004). Frequentist and Bayesian approaches for interval-censored data. Statistical Papers, 45, 139–173. Hough, G., Garitta, L., & Sa´nchez, R. (2004). Determination of consumer acceptance limits to sensory defects using survival analysis. Food Quality and Preference, 15, 729–734. Hough, G., Langohr, K., Go´mez, G., & Curia, A. (2003). Survival analysis applied to sensory shelf life of foods. Journal of Food Science, 68, 359–362. Hough, G., Puglieso, M. L., Sa´nchez, R., & Mendes da Silva, O. (1999). Sensory and microbiological shelf-life of a commercial Ricotta cheese. Journal of Dairy Science, 82, 454–459. Klein, J. P., & Moeschberger, M. L. (1997). Survival analysis, techniques for censored and truncated data (502p). New York: Springer-Verlag. Kleinbaum, D. G. (1996). Survival analysis, a self-learning text (324p). New York: Springer-Verlag. Meeker, W. Q., & Escobar, L. A. (1998). Statistical methods for reliability data (680p). New York: John Wiley & Sons. Spiegelhalter, D. et al. (1996). Bayesian inference using Gibbs sampling, Version 0.5, (version ii). Cambridge: MRC Biostatistics Unit.