Time dependence in joint replacement to multi-products grouped. The case of hospital food service

This paper proposes a stochastic inventory policy of periodic review with joint replacement for multi-products applicable to food service. The products to supply are grouped by multivariate cluster analysis by features of interest in this service. The associated inventory model that describes the random demand for each product considers forecast temporal dependence through a generalized linear autoregressive moving average (GLARMA) model. Optimization of the cost function related to the inventory model is obtained considering the expected value of the forecast of demand per unit time and stochastic programming. The proposed policy is exemplified with real-world data from a Chilean food service hospital. Subjects: Business; Management and Accounting; Operational Research/Management Science; Operations Management; Production, Operations & Information Management; Research Methods in Management; Supply Chain Management


Introduction
Supply systems have been shown to reduce the vulnerability of managing the supply chain companies. This reduction is achieved by optimizing inventory levels to meet the demand for components and finished products in enterprises (Hillier & Lieberman, 2005).
In inventory models for random demands with periodic review (PR), this random variable (RV) can be described by a continuous or discrete distribution. In order to achieve savings in total inventory costs (TC) when the cost of placing an order to the supplier for a number of different products has two components: (a) a significant cost ordering regardless of the number of different products in the order, and (b) a minor order cost of depending on the number of different products in order, several products can be replenished grouped together with the same cycle of PR for all assortment of inventory (Khouja & Goyal, 2008). Van Eijs, Heuts, and Kleijnen (1992) has proposed algorithms to find solutions to the problem of joint replacement (JR). to solve the problem of JR can be classified into two types: a direct strategy group (DGS) and an indirect strategy group (IGS). Under DGS, the products are divided into a predetermined number of sets and products within each set are reset together with the same cycle time. Under IGS a reset is performed at regular time intervals and each product has a sufficient amount of replacement to last exactly an integer multiple of regular time interval. IGS groups indirectly formed by products with the same integer multipliers.
Multivariate analysis Cluster proposes a way to achieve groups of items, taking into account the minimum distance from each multiple variables that can be used to characterize the groups formed (Johnson & Wichern, 2007). In this case, we propose multivariate cluster analysis to use DGS in JR problem.
Once defined indicators and assumptions distribution of demand per unit time (DPUT) and demand during the delivery time (LT), named LTD, have been established, the expected value of the objective function based on TC must be optimized (Namit & Chen, 1999). Often RVs DPUT are not independent and identically distributed (IID) in the time (Cysneiros, Paula, & Galea, 2005). When you want to optimize inventory consists of several products, it should not be spared the possible time dependence of the demand for products that make up the assortment of inventory (Calfa, 2015). Such optimization can be improved if it is considered the random demand forecast by conditional probability distributions to the latest information, such as generalized linear model of autoregressive moving average (GLARMA) (Dunsmuir & Scott, 2015), this is achieve by the density forecast (DF) technique, with which is possible to obtain the parameters that will require stochastic programming of inventory models of PR in JR problem.
The aim of this paper is to propose a methodology based on stochastic inventory models considering forecast of the demand dependent on time for JR for several products grouped by various variables supplied from a single distributor. The presentation of this research is organized as follows: Section 2 makes a literature review about: (i) modeling DPUT of products conditionally to past information based on the GLARMA model, (ii) inventory models of PR with JP based in cluster analysis, (iii) stochastic programming and (iv) financial indicators of the inventory policy. Section 3 exposes the proposed methodology, whereas Section 4 illustrates it taking foods supply of a Chilean public food service hospital with real-world data; and Section 5 provides a discussion of the results obtained in this research, as well as its limitations and future research.

Literature review
Usually RVs DPUT are not independently distributed on the time. This dependence can be modeled by time series as autoregressive and moving average (ARMA) model (Cysneiros et al., 2005). GLARMA models consider ARMA components but transforming the mean of the data with a link function in the line of generalized linear models (GLM). GLARMA models are widely flexible, easy to estimate and interpret, and besides its prediction is straightforward (Dunsmuir & Scott, 2015). Prediction based on GLARMA models may be carried out by the density forecast (DF) technique; for more details, see Diebold, Gunther, andTay (1998), Diebold, Hahn, andTay (1999), Bauwens, Giot, Joachim, and David (2004) and Calfa (2015).
A multi-product problem that commonly occurs in inventory policies is to decide which optimal amounts of the products must be ordered simultaneously from a single source. The cost of placing an order for multiple products has two parts: (i) one that is independent of the number of units of products to be ordered (high cost of order), and (ii) another depending on items to order (order of minor cost). This problem is known as JR and considers the savings obtained with an inventory policy which higher cost of products ordered multiple to a supplier is less than the cost of ordering items.
Many algorithms have been proposed to find solutions to the problem of JR (Van Eijset al., 1992). However, when supply is a single distributor the major ordering cost for a group of products is equal to a minor ordering cost of one single product.
Multiple products can be jointly replenished with the same cycle (Khouja & Goyal, 2008) by a grouping them into a fixed amount of product sets. To achieve grouping jointly products supplied, you can use the cluster analysis, which includes individuals considering the distances between them, with respect to different variables, forming homogeneous groups but are heterogeneous between them. Because the technique assumes that the grouping variables should not be correlated with each other, it is possible to first apply principal component analysis or factor analysis (Rojas, 2016). For details about the cluster analysis, see Johnson and Wichern (2007).
Periodic review policy provides solutions to inventory management problems in many real-life situations. This policy indicates that the action has been completed at periodic intervals to reach a predefined level of inventory (Lau & Lau, 2003). PR for multiple products can include restrictions on storage volumes, maximum amount of orders per period and maximum availability of capital stocked inventory, among others (Lee & Lee, 2013).
The stochastic programming problems is mathematical programming whose in formulation appears some stochastic element whose value is unknown but can estimate its probability distribution (Shapiro, Dentcheva, & Ruszczynski, 2014). The expected value of the objective function based on the inventory total cost (TC) must be optimized in PR policy for JR of a assortment of products. To reach this expected value can be considered various future scenarios for this RV. This can be accomplished by knowing the conditional density over time of this RV as to project future values. Then, future scenarios are constructed by performing simulations in parallel variable interest based on the prediction of conditional density over time, see Gülpnar, Rustem, and Settergren (2004). Stochastic programming can be used to solve this optimization problem by the differential evolution (DE) algorithm, which belongs to the family of genetic algorithms, mimicking the process of natural selection in an evolutionary manner; for details about the DE algorithm, see Price, Storn, and Lampinen (2006), Thangaraj, Pant, Bouvry, and Abraham (2010) and Wanke and Leiva (2015).
The benefits of inventory models are evaluated by efficiency indicators. Some of them are financial, such as TCs type inventory or sales relative rotation of stocks available; while others are operational, such as the ability to meet demand with inventory planned by the inventory model (called "fill-rate"), or the expected shortage per cycle. In this paper will deal TCs as financial indicator of an assortment of products (Wanke, Ebwank, Leiva, & Rojas, 2016).
In this work we propose a methodology based on inventory models for JR for multiple products grouped by several variables, for supply from a single distributor with DGS. In application it is carried out by cluster analysis to identify homogeneous products considering multi variables that can serve to characterize the groups to form for supply a assortment of products. We consider the DPUT of each product is an RV following a probability distribution, which has dependence over time being described by GLARMA models. Optimization of the inventory TC is carried out by stochastic programming. We apply the proposed methodology to a food service study, where foods are segmented forming homogeneous clusters, but heterogeneous between them. Segmentation is used to determine unique replacement periods in each group of foods, according to the characteristics and appropriate periods of the members of each group supply from a single distributor with DGS.

Demand during lead time plus review period
Let Y t be an RV corresponding to the DPUT at time t. In addition, let the RV W be the time of review (R) plus the LT (L) between the ordering of a product and its delivery, which has mean is obtain as a weighted sum of future scenarios of this RV, for which it is necessary to know the DF at the time of this variable of interest. On the other hand, variance Var(Y t ) = 2 Y t . Moreover, let D be the DLT plus review period (DLT + RP) for the product, which is the random sum given by with forecast probability density function (PDF) f D defined on [0, ∞) (non-negative support), cumulative distribution function (CDF) and quantile function (QF) d(q) = F −1 D (q), for 0 < q < 1. The expectation and variance of D are, respectively, expressed as and Note that, in general, W and Y t can be modeled by any discrete or continuous distribution.

GLARMA model
The assumption of IID RVs for DPUTs Y t depending over time t, for t = 1, … , m, and for the DLT D can be violated. GLARMA models are derived from GLM, with the response variable belonging to the -parametric exponential family of distributions with PDF given by where a t and c t are sequences of constants. Information on past time Ω t−1 is summarized in the state variable t . The systemic component of model GARMA(p, q) is described by a link function g of the mean where, t = g( Y t | t ) corresponds to the mean of the variable of interest to be predicted conditional to past information Ω t−1 given by where and corresponds to the ARMA components of a model of orders p and q, respectively. Note that ⊤ t = ( 0 , 1 , … , n ) is a vector of coefficients associated with n covariates with dependence over time denoted by The link function of the GLARMA model generally is the identity or logarithmic function and the corresponding model variance is assumed to be constant over time (Benjamin, Rigby, & Stasinopoulos, 2003).

Parameter estimation
Given n observations of Y t , for t = 1, … , n, the likelihood function is constructed as the product of conditional PDFs of Y t given Ω t−1 . The state vector t is a function of the GLARMA parameters ⊤ , ⊤ and ⊤ , where each time embodies these conditioning variables. Therefore, the corresponding loglikelihood function is given by By differentiating (1) for each of the parameters, it is possible to obtain their maximum likelihood (ML) estimates.

Density forecast
When formulating a statistical model, sample data are used to estimate its parameters. Then, we can employ the estimated model for diverse purposes. However, one has no idea about how good and/or reliable its predictions are out of this sample. Thus, one could evaluate the out-of-sample performance of the model by splitting this sample into two parts, one of them being a training set that could correspond to the first two thirds of the data. Hence, one can estimate the model parameters with the training set and then evaluate out-of-sample performance with the third remaining of the data. This is known as out-of-sample testing. If one estimates the model parameters and verifies the assumptions on which this model relies using the data sample with all of the thirds, then this is an in-sample testing. Thus, out-of-sample testing should be considered in the model validation process. Out-of-sample predictive performance can be evaluated by the DF technique. Specifically, let denote for a sequence of PDFs defining the data generating process. In the DF technique, it must be checked whether or not. However, to test the truth of (2) is not an easy task, because {p Y t (y t | t )} is not observed. We can use the probability integral transform given by to carry out the mentioned task. Diebold et al. (1998) demonstrated that, under a null hypothesis based on (2), the sequence of probability (3), is generated from IID U(0, 1) RVs. Thus, first, one can consider a graphical analysis by using the histogram and autocorrelation function (ACF) and/or partial ACF (PACF) plots of the sequence {z t } to detect departures from the IID U(0, 1) assumption. Then, the Anderson-Darling (AD) and Kolmogorov-Smirnov (KS) tests and the Ljung-Box (LB) test can be used to corroborate goodness of fit and independence, respectively, detected in the graphical analysis of {z t }. This indicates to us that the empirical sequence of probability integral transforms implied by the model given in (3) can be used for testing the hypothesis based on (2) and, consequently, the out-of-sample performance of this model. We propose Algorithm 1 for out-of-sample testing, which summarizes the steps of the DF technique.
Algorithm 1 DF technique Forecasts using GLARMA models may be carried out in an analogous way to that employed for GARCH models; see Tsay (2009). Thus, based on GLARMA models and supposing that the forecast origin is j = m and its horizon is h, we have the h-step ahead forecast is obtained from y m+h , with initial prediction value y m at the origin m and forecast error e m (h) = y m+h − y m (h), for h ≥ 1.

JR with inventory model PR to group of product using multivariate cluster analysis
We will compare cost savings obtained for an individual and JR for a set of products, using a inventory model of RP. This inventory system checks inventory position periodically and reorder the DLT plus review period. The optimal values for the cycle time (R) and the maximum inventory level (which is a function of safety stock factor, k p ) should be determined to minimize joint cost function of order, storage and shortage. However, some functions and demand distributions make difficult to solve the inventory problem analytically. As a result, approximations and heuristic are often used in PR inventory system. We assume normality for DLT plus the review period.
The objective is to develop an effective heuristic to solve a multi-item PR inventory in JR problem with a review period constraint under stochastic demand.
To form clusters of products with homogeneous behavior intra groups and heterogeneous behavior inter groups for several grouping variables, it is necessary to evaluate the correlation between these variables. If correlations exist, it is necessary to reduce dimensionality by factor analysis and form new uncorrelated grouping variables. The extraction method of variance with varimax rotation ensures that each extracted factor possesses maximum shared variance with the original variables that composes it. The criteria to carry out the reduction of dimensionality are: themselves > 1 for each factor, the sum must collect inertia values above 70%, confirming the statistical significance of the factor loadings with magnitude greater than 0.7, and discarding loads < 0.5. With the factor scores for each product, we proceed to conduct the cluster analysis. First, a hierarchical clustering of standardized factor scores of each product can be performed by linking the Ward algorithm, which minimizes the intra-group variance and finds a dendrogram finds with many preliminary groups shape. Knowing that, at a smaller distance, clusters are more homogeneous, the process is stopped when the horizontal lines are very long. Members and characteristics of each hierarchical clustering preliminary factorial scores are found occupying no hierarchical cluster analysis in k-means, which must be confirmed by an ANOVA (Rojas, 2016).

Stochastic programming
The problem of stochastic programming to be formulated should consist of: (i) a single decision variables of first stage of stochastic programming involved in each model adminstration inventory associated with the product to optimize; and (ii) maximize the objective function corresponding to TCs, that depend on multiple future scenarios can be modeled by simulation in parallel once known the forecast demand PDF (Gülpnar et al., 2004). As the decision variables depend on the forecast of the random demand for future stochastic programming stages, and ensure that these variables will have a single value, the conditions of non-anticipivity of the problem it is relax and then the results of the decision variables are restricted to be equal in all scenarios predicted. The TC of the product i is expressed as the weighted sum of the costs (in monetary units) in different scenarios: (i) generation of purchase orders per order (O i ) per year; (ii) storage per unit of product (H i ) per year; and (iii) of shortage or rupture per unit of the component (S i ) a year. It is assumed that demand is continuing every day of the year. The optimization problem of the expected inventory TC can be visualized as a model of stochastic programming formulated as min{Z j = C(R j , k p i )}, subject to a R j > 0, where the group j is a subset of products i = 1, … , n. Here, R j is a review period for the group j, and k p i is a secure factor of the product i, in a PR system for all products i of the group j. The review period R j may require restrictions. For example, in supply of foods in public hospitals, replacement time cannot exceed the period of expiration of the products (E i ), then R j ≤ E i . As supply is a single distributor, only the major ordering cost for a group of products (O j ) is considered, and is equal to a minor ordering cost of one single product (O i ) that not considered in supply of group of products. Other constraints are: S i > 0, H i > 0, O i > 0, i > 0, where i is the rate of DPUT of a product i conditional on past information and i the standard deviation which it remains constant at all times.

Inventory model PR for multiple products
The inventory model of PR for all products i = 1, … , n in the group j in JR minimizes the expected annual costs min E(TC i ) = min E(O i + H i + S i ) . This sum of terms is described as http://dx.doi.org/10.1080/23311916.2016.1251029 subject to R j , H i , O j , D i , i > 0, and R j ≤ E i , for all group j, where k p i is the optimal standardized QF of the product i associated with a service level of p × 100 (or amount of standard deviations-SDs-of the DLT plus review period). Both R j as k p i correspond to decision variables that minimize the inventory TC.
For the members of the product i of each group j, each element has a cycle unrestricted supply (R j ). D i corresponds to the DLT + RP for product i, while the function of shortage S(max i ), can be calculated as where mD i corresponds the maximum DLT + RP, f D i is the PDF of the DLT + RP and the maximum supply of product i corresponds to Algorithm 2 Main methodological steps 1: Collect demand data for the product i (i = 1, … , n) in each month during 2 years 2: For the statistical analysis 2.1 Carry out an autocorrelation study for data collected in Step 1 examining plot of the ACF and partial autocorrelation (PACF) to detect possible dependence to past information 2.2 Propose GLARMA models considering negative binomial distribution for the response variable (which is a monthly count is sometimes zero), for each product, selecting those through the Akaike criteria was lower Note that for modeling the shortages caused by the DLT plus review period, we assume that the variable D i is normally distributed according a best god fit distribution, because although data of demand for products have a discrete nature, usually, these are modeled by continuous distributions.

Summary of the methodology
Algorithm 2 summarizes our methodology in four main steps divided in 15 sub-steps based on the aspects detailed in Subsections 3.1 to 3.7, from the collection of data until the establishment of the TCs to evaluate the optimized group JR system in relation to the individual optimized system. We recall this algorithm considers the demand dependent to past information for products, but once all the products are considered, the total contribution of the products used in the service are maximized.

Application
The foods supply in food service units of Chilean primary health centers is channeled through their central warehouse, which acts as an intermediary between suppliers and output units (OU). The OUs receive the demand for foods, including its own food service, which performs dispensing of prescriptions foods to patients. This warehouse needs the storage, conservation and distribution of such foods. Supply of warehouse is carried out by a supplier. The warehouse delivers the products on a monthly basis to all OUs by using aggregated demand requirements for each of them in the same period.

The data set
To validate the proposed methodology, we use real-world monthly demand data for an assortment of 223 foods products, which it has extracted a set of nine products autocorrelated monthly demands used in this example. They are shipped from the warehouse and delivered to a family health center, located at the city of Concon, Chile, during 48 months of the years 2012-2015 (from January 1 to December 31). The data set collected for each variable considers the following grouping of products by variables: (a) generic product (only identification), (b) expected value of the quantities demanded depending on the time, (c) SD of the monthly demands, (d) unit costs of purchase, (e) food form (ordered according to the storage volume), and (f) diet-therapeutic use (ordered according to the urgency of use).

Computational framework
R is a non-commercial and open source software for statistics and graphs, which can be obtained at no cost from http://www.r-project.org. The statistical software R is being currently very popular in the international scientific community. For a use of this software in inventory models, see Rojas, Leiva, Wanke, and Marchant (2015) and Rojas and Leiva (in press). Some R packages related to statistical distributions that may be useful in inventory models are available at http://CRAN.R-project. org (Barros, Paula, & Leiva, 2009;Dunsmuir & Scott, 2015;Leiva, Hernandez, & Sanhueza, 2008). The expected value for monthly demand products is calculated by using GLARMA models implemented in the R software by the packages glarma() for time series with different discrete probability distributions. To build groups of foods, we carry out factor and cluster analyses by variables, using an R package named cluster, while the stochastic programming of inventory models is performed by the packaged DEoptim() of the same software.  Table 1 shows an descriptive statistics and Table 2 shows an illustrative example for a set of products with the best fit to GLARMA models with parameters of the negative binomial distribution chosen by the Akaike information criteria (AIC), displaying: parameter of NBI distribution of DPUT ( Y t i ) non conditional over time and positive dispersion parameter ( i ), order ARMA, autoregressive (AR), moving average (MA) and intercept coefficients with standard error (SE) respective, and value of AIC. The models exclude explanatory variables.

(c)
To confirm the correct fit of the proposed GLARMA model are examined six plots for checking: Observed time series related to fixed effect of GLM estimation or GLARMA estimation, disposition of Pearson Residuals in the time, Histogram of Uniform PIT, Histogram of Randomized Residuals normalized (randomized for a discrete response distribution), Quantile-Quantile (QQ) Plots for randomized residuals of a fitted GLARMA object (plot(glarmamod)), referred to as residuals below: a plot of the ACF of the residuals, see example for product 1 in Figure 3. This results indicated that is possible apply a GLARMA model to original data, considering a covariate and logarithmic link function to the mean of response variable in each time of its realization. The selection of order the best GLARMA model is according Akaike criteria. Also check Box-Pierce test type Ljung to corroborate aleatory disposition to 2 lags with statistic X 2 = 1.21, degree of freedom (df) = 2, and p-value = 0.5456, and normal distributions of randomized residuals of model, with Shapiro-Wilk normality test obtain     Figure 4 between observations, their fitted mean at time t, and the monthly prediction for one year of the quantity demanded given its conditional probability distribution over time. The forecast PDF serves to simulate future scenarios of demand. Figure 5 show the dendogram of groupings of products formed from scoring of principal components non-correlated realized by method hierarchical of Ward. Note that the hierarchical clustering and subsequently corroborated by the k-means method results in the formation of four groups of  products. Tables 3 and 4 show a summary of policy for PR of products grouped and individual products, where order cost is 97.83 USD/order, holding cost is 0.34 USD/unit/year and penalty cost for shortage is 0.12 USD/unit/cycle. Note that the TC for grouped products are 6% lower than to individual policy and shortages units expected decrease significantly in grouped policy.

Discussion
In this investigation, we has explored the possibility to perform forecasts based on a conditional autoregressive model whose quality could be evaluated by the probability integral transformed, would update in n-steps ahead the myopic policy inventories. The policy inventories raised show an interesting way to achieve significant savings in TC and shortage expected for multi-product systems. The grouping of products that achieves multivariate cluster analysis considers multiple variables that influence the characterization of supply. The proposed supply cycle for groups of products is useful in institutions that have a supplier to meet your requirements, and that generally are characterized by high bureaucracy in your administrative systems as public hospitals, where the operation of a supply system as proposed would be facilitated (Bennett & Gilson, 2001;Birkett, Mitchell, & McManus, 2001). The parameterization of the mean of a PDF under a temporal structure of a conditional autoregressive model to information passed as GLARMA, achieves an efficient and dynamic characterization of a myopic policy involving the DPUT inventories and DLT plus PR. Although a normal distribution deal in this work, the proposal is valid for any PDF that can be parameterized with respect to the conditional mean of a time series model. This condition opens multiple areas of research in the area. In the future, it is possible to obtain improvements in the modeling raised occupying other than normal distributions, and even considering that in many cycles demand component is zero, requiring demand model using probability distributions for continuous and/or discrete data zero-inflated.