Optimal Sampling Schedule of Diet Components: Model Robustness to Departure from Assumptions1
Article Outline
- Abstract
- Introduction
- Materials and Methods
- Results and Discussion
- Conclusions
- Acknowledgment
- Supplementary data
- References
- Copyright
Abstract
Monitoring the nutritional composition of forage can be conceptualized as a quality control process that can be accomplished using control charts such as the Shewhart X-bar chart. A sampling schedule for an X-bar chart is defined by 3 parameters: the number of samples n taken at each sampling time, the sampling interval h, and the control limits L. All 3 parameters affect the performance of the chart, and thus, the total quality cost (TQC). A TQC function consists of cost per cycle while the process is in-control, cost per cycle while the process is out-of-control, cost per cycle for sampling and analyses, and the expected duration of a cycle, with a cycle defined as the time between the start of successive in-control periods. A general TQC function was derived for a renewal reward process. Optimization of this TQC function allows for the determination of the optimal n, h, and L values that minimize the total daily quality costs. The model assumes an abrupt change in composition of the forage when the process goes out-of-control. It also assumes a normal distribution of the measurements when the process is in-control and an absence of outlier measurements. The objective of this research was to evaluate model robustness to departure from these 3 basic assumptions. A series of Monte Carlo simulations was performed while varying the average time that the process is in control from 5 to 90 d using 1) errors of measurements that follow a standard normal distribution (SN); 2) SN with ± 3.5 SD outliers inserted with a frequency of 1, 5, and 10%; 3) log normal error of measurements with SD = 1; and finally 4) SN with a gradual shift of the mean from 0 to 1.5 SD over 7, 14, and 28 d. The model was very robust to the presence of outliers; the average change in TQC was less than 1%, even with a frequency of outliers of 10%. The model is also very robust to asymmetry in the distribution of the measurements (i.e., probability distribution function with a long right tail): the log normal distribution, as opposed to the assumed normal, resulted in an increase in TQC of less than 1.4%. Finally, the gradual shift in mean composition did not result in an increased TQC but in a 17.3% decrease compared with an abrupt change. The model appears very robust to departure from normality, presence of outliers, and a skewed distribution of measurements. Gradual changes in the process are readily detected by the optimum X-bar chart with the conventional decision criterion, and monitoring performance is not markedly improved by augmenting the number of decision criteria in the X-bar chart or by the addition of a cumulative sum chart. Because of its robustness, the model can be applied to optimize forage sampling on dairy farms, with expected savings ranging between $80 and $100/cow per yr.
Key words: optimal sampling schedule, forage composition, robustness, control charts
Introduction
Variation in the nutritional content of animal diets is receiving increased attention because of its effects on animal performance and profitability (Fadel et al., 2006; Kohn, 2006; St-Pierre and Weiss, 2006). Control of this variation can occur during diet formulation (St-Pierre and Harvey, 1986; Tozer, 2000), diet preparation (St-Pierre and Weiss, 2006), or at the time feed ingredients are sourced. Sourcing can reduce nutritional variation either by increasing purchases of commercial feed-stuffs from suppliers of less variable feedstuffs or from the implementation of quality control programs on farm (St-Pierre and Weiss, 2006). Forages are among the most variable feed ingredients, whether they are produced on farm or purchased (Mertens, 2006). They generally account for over 50% of total diet DM fed to lactating cows and are often produced on the dairy farm. In these instances, forage inventories are inevitably commingled because of the enormous capital requirement that would be associated with an effective separation of inventories of uniform composition. From a quality control standpoint, the interest is then to detect the inevitable abrupt composition change that occurs with forages. St-Pierre and Cobanov (2007) studied this stochastic process and proposed a model for a total quality cost (TQC) function based on renewal reward process theory (Lorenzen and Vance, 1986; Parzen, 1999).
The model derived by St-Pierre and Cobanov (2007) requires 13 input variables and is optimized with respect to 2 integer and 1 continuous sampling design parameters used to define a Shewhart X-bar control chart (Shewhart, 1931; St-Pierre, and Cobanov, 2007). Four assumptions were made during model derivation. First, it was assumed that the process conforms to a Poisson renewal reward process, implying that the process is memoryless and cannot return to an in-control state without intervention. The first part of this assumption seems rather obvious. As for the second part, as long as the interval between major compositional changes is much greater than the time to sample, analyze, and chart the results, the renewal reward process approximation should yield results that are very close to the true process (Knappenberger and Grandage, 1969). Second, it was assumed that when the process goes out-of-control, it shifts abruptly to a different state, meaning that the change in forage composition is sudden. Third, it was assumed that measurement errors are normally distributed, implying 2 other assumptions: the distribution of measurements while in-control is symmetrical, and there are no outliers.
There are 2 strategies used to assess a model with respect to its assumptions. The first consists of studying the physical process (i.e., make measurements) to demonstrate that the process does conform to the assumptions. The second is to determine whether a given assumption really matters; that is, to estimate the magnitude of error in the model and its answers when an assumption is clearly violated. In instances in which the model is not very sensitive to an assumption, it is declared robust and no further studies of the physical process are required in regard to the assumption in question (Ishikawa, 1976; Yashin, 1993).
The objective of this research was to determine the degree of robustness of the TQC model of St-Pierre and Cobanov (2007) to departure from 3 basic assumptions used during its derivation.
Materials and Methods
Model
The TQC function of St-Pierre and Cobanov (2007) is as follows:
[1]Table 1 summarizes the input variables and reports the input values used for the robustness assessment.
Table 1. Input variables to the total quality cost model
| Input variable | Description | Input value |
|---|---|---|
| 1/λ | Mean time process is in control (d) | 5 to 90 |
| Δ | Number of standard deviations slip when process moves out-of-control (abrupt change in forage composition) | 1.5 |
| T0 | Expected search time for false alarm (d) | 2 |
| T2 | Expected time to repair the process (i.e., to rebalance and implement new rations) (d) | 2 |
| E | Time to sample and analyze all items (i.e., time awaiting for lab results) (d) | 4 |
| Nc | Number of cows in the herd | 100 |
| Δm0 | Milk production loss due to white noise change in forage composition while the process is in-control (kg/cow per d) | 0.45 |
| Δm1 | Milk production loss when process goes out-of-control (abrupt change in forage composition) (kg/cow per d) | 2.27 |
| Pm | Price of milk ($/kg) | 0.286 |
| Y | Cost per false alarm ($) | 200 |
| W | Cost to locate and repair the assignable cause (i.e., to rebalance and implement new rations) ($) | 200 |
| a | Fixed cost per sampling time (preparation and mailing of samples) ($) | 15 |
| b | Cost per unit sampled ($) | 25 |
Robustness Assessment
SimulationA simulation model was constructed to study the consequences of a sampling design with parameters n, h, and L determined from the optimization of [1] but applied under conditions in which an assumption is clearly being violated. In each simulation, 1/λ was varied from 5 to 90 d, because this input variable was previously identified as one of the most influential in the determination of the optimum sampling design (i.e., in identifying the optimal n, h, and L; St-Pierre and Cobanov, 2007). For each value of 1/λ, a total of 5,000 simulation runs was executed. Each run consisted of 10,000 simulation days, yielding 50
million simulated days per result point. An event scheduling approach (Fishman, 1978) was used for simulating the basic process illustrated in Figure 1. Pseudo-random variables were generated for the Gaussian, negative exponential, and log-normal distributions using the algorithms of Fishman (1978).
Figure 1.
Graphical representation of a quality cycle; 1/λ = mean time process is in-control, τ = expected time of occurrence of abrupt change after last sampling while in-control, Tc = expected time between the occurrence of abrupt change and the next sampling time, h = sampling interval, Te = expected time until an out-of-control signal occurs, T1 = expected time to investigate the cause of the change, T2 = expected time to reformulate and implement new diets, and Tf = expected time to investigate the cause of the change, reformulate and implement new diets.
To assess model robustness in the presence of outliers, sample measurements were generated from a standard normal distribution with outlier frequency of 1, 5, and 10%. Outliers were set at values of ± 3.5 SD at random.
Robustness to the Assumption of a Symmetric DistributionTo test robustness to the assumption of symmetric distribution of measurements when incontrol, sample measurements were generated from a log-normal probability distribution with SD equal to 1. Because the log-normal probability distribution is defined only for the positive values, the whole probability distribution was centered to obtain an expectation of 0 (
). The probability distribution function of the resulting log-normal is compared with the standard normal distribution in Figure 2, showing the heavy right tail of the log-normal.
Figure 2.
Probability distribution functions of the standard normal (solid) and log-normal (dashed) distributions with equal means (= 0) and standard deviations (= 1); x is a standardized random variable and P(x) is the value of the probability functions at x.
Robustness to this assumption was assessed by changing the mean composition of the forage from 0 to +1.5 SD over 7, 14 or 28 d, and generating sample measurements from a standard normal distribution with a mean equal to a linearly interpolated value. For example, if forage composition was to change gradually from 0 to 1.5 SD over a period of 28 d, then a sample taken at 14 d would be generated from a normal distribution with a mean of 0 + 
, and an SD of σ. The lowest value of 7 d is representative of the progressive but relatively rapid change in a tower silo with a diameter of 6.1m, a vertical unloading rate of 0.15m/d, and an angle of repose of 30°. The greatest value of 28 d is representative of a large horizontal silo, filled with a 30° wedge, and unloaded at the rate of 0.2 horizontal m/d.
Augmenting the Charting Decision Criteria
Shewhart control charts are generally effective at detecting large shifts in processes, but relatively ineffective at detecting small shifts (Reynolds and Stoumbos, 2005). Many authors have promoted the idea of either augmenting the standard X-bar chart with additional decision criteria (Cook et al., 1992) or to use combinations of charts (Westgard et al., 1977; Lucas, 1982; Blackwell et al., 1994; Reynolds and Stoumbos, 2005). Because the shift is relatively small in the first half of a gradual shift (i.e., less than 
), it is possible that the X-bar chart is very ineffective during this time and that augmenting the chart could reduce the mean daily quality cost. Thus, an augmented charting scheme was constructed combining the following 3 approaches:
Results and Discussion
Performance of Standard X-Bar Chart
The performances of the standard X-bar chart and the simulation model when assumptions are met are reported first to verify the accuracy of the simulation model and to establish a framework of comparison for when the assumptions are not met.
Figure 3 reports the optimal parameter values (n*, h*, and L*, where the symbol * signifies the optimal values of n, h, and L) for the sampling design at various values of 1/λ. Total quality costs were calculated either analytically using [1] (an exact value under the assumptions) or from the simulation model. Results showed a very tight agreement between the analytical solution and that estimated by the simulation model, thus verifying that the simulator was likely correctly coded.
Figure 3.
Optimal sampling schedule, total quality cost calculated analytically, and cost obtained by the simulation model as a function of the mean time process is in control (1/λ): —— = the sampling interval (h); —Δ— = number of standard deviations used as control limits of the X-bar chart (L); — ♢ — = number of samples taken (n); —■— = the total daily quality cost calculated using the analytical formula (C); and —○— = the total daily quality cost calculated through the simulation process.
Figure 4 shows the quality outcome from the optimal sampling design for various values of 1/λ when assumptions were met. Across the full range of 1/λ, the frequency of false alarms remained small, peaking at a frequency of 172 (±12) occurrences per 104 d (i.e., an average of 58.1 d between false alarms) at 1/λ = 18 d. The sharp drop in the frequency of false alarms occurring between 1/λ = 18 and 19 d was due to a change in the optimal number of samples n* changing from 1 to 2 between these 2 values of 1/λ (Figure 4). The average run length while the process was out-of control (ARL1) was short across all values of 1/λ. In fact, ARL1 was always shorter than the optimal sampling interval (h*) across all values of 1/λ. This indicates that the X-bar chart with the optimal sampling design is very efficient at detecting changes in forage composition for the values of input variables used. This is a consequence of a relatively low cost for false alarms compared with out-of-control costs.
Figure 4.
Frequency of false alarms, average run length when out-of-control (ARL1), and total quality cost per day as a function of the mean time process is in control (1/λ) when all assumptions are met: — ♢ — = average run length; —Δ— = frequency of false alarms; and — ○ — = total quality cost calculated through the simulation process.
Robustness to Presence of Outliers
Figure 5 shows the quality outcome using the sampling design optimized without the presence of outliers when 10% outliers; (±3.5 SD) were introduced in the simulated samples. The frequency of false alarms reached about 190 (±15) occurrences per 104 d (i.e., an average of 53.0 d between false alarms) at 1/λ = 18 d. This is a small, nonsignificant (P>0.1) increase over what was observed in the absence of outliers (172 occurrences per 104 d). In addition, the average run length when the process goes out-of-control was even less than what was observed in the absence of outliers. This indicates that the X-bar control chart operates very effectively, even in the presence of 10% outliers; it readily detects when the process goes out-of-control without significantly increasing the number of false alarms when the process is in-control. Also, the resulting TQC is identical to what was observed without outliers across all values of 1/λ. Thus, the TQC model is very robust to the presence of outliers in the measurements.
Figure 5.
Effect of the presence of 10% outliers (±3.5 SD) on the frequency of false alarms, average run length when out-of-control (ARL1), and total quality cost per day as a function of the mean time process is in control (1/λ): —♢— = average run length; —— = frequency of false alarms; —○— = total quality cost calculated through the simulation process; and —Δ— = total quality cost calculated using the analytical formula.
Robustness to the Assumption of a Symmetric Distribution
Figure 6 shows the quality outcome using the sampling design optimized with the assumption of normally distributed samples while the true samples were from a heavily skewed log-normal distribution. The frequency of false alarm when the process is in-control reached a maximum of 147 (±10) occurrences per 104 d (i.e., an average of 68.0 d between false alarms) at 1/λ = 18 d, a value which is not significantly less (P>0.1) than the 172 occurrences per 104 d observed when the assumption were implicitly met (Figure 4). Also, the average run length when the process goes out-of-control was less than what was observed under the normality assumption (Figure 4). Thus, the X-bar chart with the optimized sampling design is very efficient at detecting out-of-control states, even when the distribution of the samples is heavily skewed to the right. The realized TQC under the log-normal distribution were very close to the analytical TQC calculated with the normality assumption. Consequently, we conclude that the TQC model is very robust to departure from the symmetry distribution of sample measurements.
Figure 6.
Effect of asymmetry (log-normal distribution of samples) on the frequency of false alarms, average run length when out-of-control (ARL1), and total quality cost per day as a function of the mean time process is in control (1/λ): ♢ = average run length; —— = frequency of false alarms; —○— = total quality cost calculated through the simulation process; and —Δ— = total quality cost calculated using the analytical formula.
Robustness to the Assumption of Abrupt Change in Forage Composition
Figure 7 shows the quality outcome when the shift to an out-of-control state occurs over 28 d as opposed to the abrupt change assumed when model [1] was derived. The frequency of false alarms peaked at 160 (±10) occurrences per 104 d (i.e., an average of 62.5 d between false alarms) at 1/λ = 18 d, which is numerically, but not statistically, less (P>0.1) than the 172 (±12) occurrences per 104 d observed when the assumption of abrupt change was met (Figure 4). The ARL1 when the change in forage composition was gradual over 28 d (Figure 7) was nearly equal to the ARL1 when the change occurred abruptly (Figure 4). Again, this indicates that the X-bar chart with the optimum design parameters is very efficient at detecting the out-of-control state without markedly increasing the number of false alarms when in an in-control state. The TQC with the gradual shift were significantly less for 1/λ less than 28 d compared with the same costs with the abrupt shift (Figure 7). The reason for this is that the costs when the process goes out-of-control are less for the gradual shift than the abrupt shift up until 28 d after the shift initiation. Because the X-bar chart is still very efficient at detecting an out-of-control state, even when the shift is gradual (less than 10 d on average from the time the process starts the out-of control state until detection across all values of 1/λ), the magnitude of production losses is actually less under the gradual shift than under the abrupt shift. This indicates that the model is very robust to departure from the assumption of an abrupt shift in forage composition.
Figure 7.
Effect of a gradual shift in mean composition over a period of 28 d on the frequency of false alarms, average run length when out-of-control (ARL1), and total quality cost per day as a function of the mean time process is in control (1/λ): ♢ = average run length; —— = frequency of false alarms; —○— = total quality cost calculated through the simulation process; and —Δ— = total quality cost calculated using the analytical formula.
Augmenting the Charting Decision Criteria
Augmenting the charting decision criteria with additional rules to the X-bar chart in combination with a CUSUM chart had little effect on the TQC across all values of 1/λ (Figure 8). This is in contrast to what has been observed by authors in other applications of control charts (Westgard et al., 1977; Lucas, 1982; Blackwell et al., 1994; Reynolds and Stoumbos, 2005). These other applications, however, were related to the monitoring of processes in which the control limits of the X-bar charts are arbitrarily set at 3 SD, in line with the 6-sigma philosophy. In our case, however, L is optimized and its optimal value was close to 1 SD across the whole range of conditions studied during the sensitivity analysis of the model. Thus, we would expect an X-bar chart to be considerably more sensitive to an out-of-control state under our conditions. In fact, the augmented criteria or the CUSUM were seldom the cause of an out-of-control signal. In our case, the cost of the extra false alarms resulting from L* smaller than traditionally had been used (i.e., L = 2 or 3) is largely offset by the shorter time to signal when there is a composition change and, thus, a much smaller cost for going out-of-control. Our results suggest that in our model, the augmentation of the X-bar chart with additional rules in conjunction with a CUSUM chart yields trivial benefits and is unnecessary.
Figure 8.
Comparison of standard X-bar chart and augmented X-bar in combination with a cumulative sum (CUSUM) chart on the total quality cost per day over a range of mean time that process is in control (1/λ): —Δ— = total quality cost calculated through the simulation process using the standard X-bar chart; and —○— = total quality cost calculated through the simulation process using the X-bar chart with augmented X-bar chart rules in conjunction with a CUSUM chart.
Limitations and Future Work
St-Pierre and Cobanov (2007) showed that total quality cost was most sensitive to Nc, Δm1, 1/λ, Pm, and E, thus requiring precise input values for these variables. Under field conditions, getting precise estimates of Nc, Pm, and E should be relatively easy. A precise estimate of 1/λ, however, would require extensive prior records of forage analyses done at a frequency seldom achieved in practice. Likewise, a precise estimate of Δm1 would be very difficult in practice. Issues of cow grouping strategies and uncertainty of requirements are poorly understood (St-Pierre and Thraen, 1999) and affect animal response to nutrients. In its current structure, Δm1 is a direct, fixed input in the model. That is, it has no distribution and no uncertainty, although it is basically an educated guess, albeit possibly the outcome of complex diet formulation systems. Subjective prior distributions of parameters are easily handled by Bayesian decision theory (De Groot, 1970), but subjectivity would certainly be a concern from a scientific point of view. In addition, optimization of [1] under a random Δm1 would be considerably more complicated than under a fixed value as is presently done.
The extent of forage variation on commercial farms is not well understood. Current estimates of composition variance (NRC, 2001) do not partition the total variance into meaningful components. For example, although analytical variance is part of the total variance reported in Table 15-1 of NRC (2001), it should be factored out when expressing the composition variance experienced by the animals (i.e., the variance in the laboratory is of no consequence to the animals). Thus, there is a need for intensive on-farm measurements to determine the components of composition variance in forages, to measure the size of the variation while in-control, and to establish the size and frequency of abrupt composition changes.
In its current formulation, our model approaches separately the sampling schedule of each nutrient in a given forage. This approach does not factor the association (correlation) between nutrients, which can be large and important. Thus, there is a need to expand the method from a univariate approach in which each composition element is monitored separately to a multivariate approach in which relationships between composition elements are factored in the monitoring process.
Conclusions
The TQC model of St-Pierre and Cobanov (2007) used to optimize the sampling schedule of forages on dairy farms appears very robust to departure from its basic assumptions of absence of outliers, symmetrical distribution of the measurements, and abrupt change in nutritional composition of forage.
Acknowledgment
This project was a component of the NC-1119 Regional Research Project. The authors thank Joanne Knapp, Maurice Eastridge, and 3 anonymous reviewers for helpful comments and suggestions on an earlier version of this manuscript.
Supplementary data
Interpretive summary.
References
- . Use of combined Shewhart-CUSUM control charts in internal quality control of enzyme-linked immunosorbent assays for typing of foot and mouth disease virus-antigen. Rev. -Off. Int. Epizoot. 1994;13:687–699
- . A knowledge-based approach to statistical process control. Comput. Electron. Agric. 1992;7:13–22
- . Optimal Statistical Decision. New York, NY: McGraw-Hill; 1970;
- . Statistical Methods for Industrial Process Control. New York, NY: Chapman and Hall; 1997;
- . Impact of variation in diet nutrient inputs on model output prediction. J. Dairy Sci. 2006;89(Suppl. 1):383;(Abstr.)
- . Principles of Discrete Event Simulation. New York, NY: John Wiley & Sons; 1978;
- . Guide to Quality Control. Tokyo, Japan: Asian Productivity Organization; 1976;
- . Minimum cost quality control tests. AIIE Trans. 1969;1:24–32
- . How can dairy nutrition models deal with uncertainty?. J. Dairy Sci. 2006;89(Suppl. 1):382;(Abstr.)
- . The economic design of control charts: A unified approach. Technometrics. 1986;28:3–10
- . Combined Shewhart-CUSUM quality control schemes. J. Qual. Tech. 1982;14:51–59
- . Quantifying assay variation in nutrient analysis of feedstuffs. J. Dairy Sci. 2006;89(Suppl. 1):383;(Abstr.)
- . Nutrient Requirements of Dairy Cattle. 7th rev. ed.. Washington, DC: Natl. Acad. Sci; 2001;
- . Stochastic Processes. Philadelphia, PA: Society for Industrial and Applied Mathematics; 1999;
- . Should exponentially weighted moving average and cumulative sum charts be used with Shewhart limits?. Technometrics. 2005;47:409–424
- . Statistical Methods for Quality Improvement. 2nd ed.. New York, NY: Wiley; 2000;
- . Economic Control of Quality of Manufactured Products. New York, NY: Van Nostrand; 1931;
- . A model to determine the optimal sampling schedule of diet components. J. Dairy Sci. 2007;90:5383–5394
- . Uncertainty in composition of ingredients and optimal rate of success for a maximum profit total mixed ration. J. Dairy Sci. 1986;69:3075–3089
- . Animal grouping strategies, sources of variation, and economic factors affecting nutrient balance on dairy farms. J. Anim. Sci. 1999;77(Suppl. 2):72–83
- . Managing feedstuff variation in nutritional practice. J. Dairy Sci. 2006;89(Suppl 1):383;(Abstr.)
- . Least-cost ration formulations for Holstein dairy heifers by using linear and stochastic programming. J. Dairy Sci. 2000;83:443–451
- . Combined Shewhart-CUSUM control charts from improved quality control in clinical chemistry. Clin. Chem. 1977;23:1881–1887
- . Statistical control schemes: Methods, applications and generalizations. Int. Stat. Rev. 1993;61:41–66
PII: S0022-0302(07)72012-X
doi:10.3168/jds.2006-728
© 2007 American Dairy Science Association. Published by Elsevier Inc. All rights reserved.
