{'fraction_non_alphanumeric': 0.03192511804772053, 'fraction_numerical': 0.009715431866616522, 'mean_word_length': 4.820625076006324, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Data-Driven Analytics for Benchmarking and Optimizing Retail Store Performance

Growing competitiveness and increasing availability of data is generating tremendous interest in data-driven analytics across industries. In the retail sector, stores need targeted guidance to improve both the efficiency and effectiveness of individual stores based on their specific locations, demographics, and environment. We propose an effective data-driven framework for internal benchmarking that can lead to targeted guidance for individual stores. In particular, we propose an objective method for segmenting stores using a model-based clustering technique that accounts for similarity in store performance dynamics. The proposed method relies on an effective Finite Mixture of Regressions technique based on competitive learning for carrying out the model-based clustering with `must-link' constraints and modeling store performance. We also propose an optimization framework to derive tailored recommendations for individual stores within store clusters that jointly improves profitability for the store while also improving sales to satisfy franchiser requirements. We validate the methods using synthetic experiments as well as a real-world automotive dealership network study for a leading global automotive manufacturer.


Introduction
Increasing global competition combined with product proliferation, dropping customer loyalties, and shrinking product life-cycles is changing the environment facing most companies today. None of the industries seem to be immune and retail and franchise sectors seem to be particularly hurting. The number of retailers filing for bankruptcy protection in the U.S. is headed toward its highest annual tally since the Great Recession in the 1920s (Gustafson, 2017). While the reasons for bankruptcies and difficulties are several, many companies within these sectors lack comprehensive and effective performance management systems (Yu and Ramanathan, 2009). These sectors have no doubt seen tremendous efficiencies from employing in-store technologies (e.g., scanning systems, enhanced point-of-sale systems, self-service lines) and information technology to drive upstream operations (e.g., warehousing, logistics, and manufacturing) (King et al., 2004). While these technologies and the associated Retail Information Systems (RIS) are effective in helping to manage core store activities (like inventory control and logistics), these chains also need analytics platforms for performance management (Rigby, 2011;Nash et al., 2013).
With the increasing availability of data and technologies to store and process the data, business analytics has been experiencing growing attention among the researchers and is transforming the way businesses operate in many sectors (Sun et al., 2017;Henke et al., 2016). In the context of retail and franchise networks, stores need customized analytical guidance to improve profitability and sales of individual stores based on their specific location, demographics, and environment. Executives are unanimous in voicing their concerns over the lack of methods to assess store specific issues and derive equally specific insights (Bucklin and Gupta, 1999). Such guidance is key to important strategic management decisions, including evaluation, promotion and development. Strategic resource-allocation decisions, such as advertising budgets, store expansions/closings, are also based on management's understanding of store performance drivers (Thomas et al., 1998).
Parsons makes the distinction between store efficiency and effectiveness. Efficiency refers to the relationship between inputs and outputs, while effectiveness focuses on outputs relative to a particular objective (Parsons, 1994) . Productivity is the combination of efficiency and effectiveness and is the focus of this study. Thomas et al. note that productivity studies demand several careful considerations (Thomas et al., 1998) . First, relevant individual store differences must be considered within the platform to take into account advantages and disadvantages of particular stores, e.g., location, competitive intensity (Kamakura and Ratchfrord, 1996). Second, development is much more effective when specific practices can be observed and transferred to other stores. Effective practices should be identified, described, and used as benchmarks for less efficient stores. This is a key focus of this study and the proposed method attempts to attribute (through modeling) variation in performance across a set of stores to the way the stores are managed. In particular, we will demonstrate how the key performance indicators (KPIs) are managed using resources. Third, a distinction must be made between resources under the control of store vs. those they have little or no influence over (e.g., local land/rent costs). In our proposed approach, we emphasis deriving recommendations around "actionable" KPIs (e.g., increasing advertising budget) vs. not so actionable KPIs (e.g., changing store location or reducing overhead). Fourth, more than one outcome usually needs to be considered because stores are responsible for multiple and sometimes conflicting performance measures (e.g., dealership sales might be more important to an automotive manufacturer vs. profits to the owner's of the dealership). This is also addressed in this study through the use of multi-objective optimization methods for deriving recommendations.
Numerous methods have been proposed for evaluating retail efficiency of individual stores (Balakrishnan et al., 1994;Kamakura and Ratchfrord, 1996). Store performance could be influenced by trade area demographic factors (Ingene and Lusch, 1980), level of competition (Craig et al., 1984), retail atmospherics (Jain and Mahajan, 1979), and promotions (Walters and MacKenzie, 1988. The effect of internal retail environment including level of service and extended store hours as well as overlapping trading areas have been studied (Pauler et al., 2009;Kumar and Karande, 2000). As for methods, Data Envelopment Analysis (DEA) has been used extensively for benchmarking performance of retail stores (Kamakura and Ratchfrord, 1996;Donthu and Yoo, 1998;Vyt, 2008). DEA is a nonparametric method in operations research and economics to empirically measure productive efficiency of decision making units. The main advantage of DEA is its ability to accommodate a multiplicity of inputs and outputs and the absence of a need to explicitly specify a mathematical form for the production function. DEA develops a function whose form is determined by the most efficient producers and differs from the Ordinary Least Squares (OLS) statistical regression technique that bases comparisons relative to an average producer. A major drawback of DEA however is that the model specification and inclusion/exclusion of variables can affect the results. It allows individual units to employ nearly arbitrary weights/importance to inputs and performance metrics to showcase their relative efficiency. For these reasons, it is more applicable for activities such as outlet termination (Vyt and Cliquet, 2017). Instead, we take the statistical regression approach for deriving improved and actionable recommendations for retail performance management.
In summary, the focus of this manuscript is to facilitate improvement in the performance of individual stores by relying on a data-driven approach to internal benchmarking. In particular, the goal is to identify factors driving automotive dealership performance in comparison with "similar" dealerships and relying on optimization to derive tailored recommendations. The problem was brought to our attention by a global leader in providing automotive dealership location and network analysis to many automotive original equipment manufacturers (OEMs). Our contributions are as follows: 1) We propose an objective method for segmenting stores using a model-based clustering technique that accounts for similarity in store performance dynamics -this is being done to better account for the considerations promoted by (Thomas et al., 1998), 2) We propose an effective Finite Mixture of Regressions (FMR) technique based on competitive learning for carrying out the model-based clustering and modeling store performance called Mixture Models with Competitive Learning (MMCL), and 3) We propose an optimization framework to derive tailored recommendations for individual stores within store clusters that jointly improves profitability for the store while also improving sales to satisfy OEM/franchiser requirements. We illustrate the methods using synthetic experiments and a real-world dataset from a leading global OEM.
The rest of this manuscript is organized as follows: Section 2 provides an overview of automotive dealership performance management and the proposed overall approach for performance management. Section 3 provides background information regarding finite mixture models. Section 4 describes the proposed mixture model with competitive learning (MMCL) for the problem of finite mixture of regressions (FMR) under group structure constraints. Section 5 presents results from synthetic experiments to validate the effectiveness of MMCL. Section 6 illustrates our method for deriving tailored recommendations using cluster specific component models from MMCL using multi-objective optimization. Section 7 describes results from a dealership case study. Finally, section 8 offers some concluding remarks and directions for future research.

Automotive Dealership Performance Management
In the automotive industry, dealer efficiency and effectiveness are key factors for obtaining and maintaining competitiveness for OEMs. This is just as critical for the well-being and durability of dealers, for most dealerships tend to be franchises (in the U.S. and much of the world) that have a contract with an automotive OEM that allows them to sell its products. It is critically important to establish analytics platforms for assessing the productivity of the dealer network that not only is useful for the OEM but also provides customized guidance to individual dealerships. Automotive OEMs usually assess dealership performance according to market share and plan incentive systems by assigning annual sales targets to each dealership. However, performance assessment based on a simplistic comparison between the dealership and national or state average market share can not only lead to ineffective sales targets but could also compromise the productivity of the dealer and the competitiveness of the OEM (Biondi et al., 2013). As noted by Biondi and co-authors, this kind of assessment does not take into account either the availability or the utilization of resources. Dealership A may be more efficient than dealership B according to the market share method, although A can obtain a higher output than B (i.e., sell more vehicles) merely because A has a more consolidated presence in the territory (i.e., has had a sales mandate for longer) and/or is located in a more favorable geographical market (e.g., where the brand enjoys more loyalty). Therefore, a more objective modeling and analysis methodology is necessary for evaluating and improving the performance of dealerships.
The company wanted to develop a dealership performance management analytics platform to analyze monthly operations and financial data (including information on sales staffing levels/tenure, product assortment/mix, dealer services (e.g., financing, trade-ins, collision repair), advertising budgets/mix, service bays/technicians etc.) from thousands of dealers in the U.S. to understand factors that can jointly improve profitability for the dealership while also improving vehicle sales to satisfy OEM requirements. In the absence of objective data-driven analytics platforms, dealerships mostly rely on experienced consultants and ad hoc guidance from field personnel. We propose a model-based "competitive learning" method for clustering the stores into similar groups for benchmarking. The proposed method is generic and can be utilized in combination with a variety of both internal (e.g. inventory, advertising) and external data sources (e.g., demographics and local competition). Figure 1 illustrates our overall approach to retail performance management analytics. The process starts with acquiring appropriate data from both internal and external sources. Next, the stores are clustered by our proposed MMCL algorithm and the important KPIs that have the highest influence on the key process objectives (KPOs) (e.g. dealership profitability and sales effectiveness) are identified. The resulting relationships are employed by the proposed multi-objective optimization (MOO) method to determine the optimal settings for each KPI to achieving the desired objectives. The outcome of the optimization and MMCL phases provides tailored guidance to stores for improving their performance and productivity. Details of each phase are discussed in more detail within later sections of the manuscript.

Finite Mixture Models (FMM)
This section provides an overview of finite mixture models for the purpose of segmenting/clustering stores using a model-based clustering technique that accounts for similarity in store performance dynamics. The objective is to cluster the stores into a number of homogeneous groups for benchmarking and deriving more effective recommendations. For example, it might be inappropriate to benchmark a rural dealer with an urban dealer in a large metropolitan city. This is above and beyond the normal practice of examining stores based on regional location. For example, it is a common practice in the U.S. for the automotive OEMs to look at the continental U.S. as several major regions (e.g., North-East, Midwest, South-East, South-West, and West) due to significant differences in weather and other purchasing patterns. While we too recommend regional analysis, there is still room for improving the performance modeling by further clustering the regional stores into a number of smaller homogeneous store groups.
In mixture models, "components" are introduced into the mixture model to allow for greater flexibility in modeling a heterogeneous population that is unable to be modeled by a single model. The hope is that this form of clustering would allow for more effective modeling and comparison of stores within individual clusters and for recognition/extraction of effective practices to be used as benchmarks for less efficient stores.

Mixture Models for Clustering
Clustering is the process of finding subsets of a dataset based on "similarity, where members of the subsets are similar and members across subsets are dissimilar (Guha and Mishra, 2016). There are several algorithms that have been proposed for the clustering problem (Xu and Wunsch, 2005). Traditional clustering methods for the most part are heuristic techniques derived from empirical methods and have difficulty taking into account the characteristics of clusters (shapes, proportions etc.). Finite mixture models have attracted much attention in recent years for clustering. McLachlan and Basford were the first to highlight the usefulness of mixture models as a way of providing an effective clustering of various datasets under a variety of experimental designs (McLachlan and Basford, 1988). FMM offer considerable flexibility and permit certain classical criteria for vigorous analysis and have been widely used for market segmentation and similar studies (Andrews et al., 2011;Wedel and DeSarbo, 2002;Sarstedt, 2008;Tuma and Decker, 2013).
In mixture models, it is assumed that observations are generated according to several probability distributions (a.k.a. components) with certain parameters. Data points in each distribution are assumed to form a cluster. The general framework of an FMM is of the following form: where k is the number of mixture components, α k is the mixing weights i.e. the proportion of component k (α k > 0 and k α k = 1), θ k is the set of parameters for the kth component, and f (x k ; θ k ) is the distribution of the kth component. Each componentm, k, is assumed to come from a unique probability distribution f (x k ; θ k ), with probability α k that an observation comes from component k. In the case of mixture of Gaussian distributions, f (x k ; θ k ) ∼ N (µ, Σ), where µ and Σ denote the mean and covariance matrix for each component distribution, respectively.

Finite Mixture of Regressions (FMR)
Given our desire to model the performance of dealerships as a function of the KPIs, we need to explicitly handle the presence of dependent variables (e.g., standardized dealership sales and profits calculated relative to averages). Quadnt and Ramsey proposed a method of moment algorithm to estimate the parameters of FMM with the presence of a dependent variable (Quandt and Ramsey, 1978). In this setup, it is assumed that the data points (X ∈ IR n ) have an associated dependent variable (Y ∈ IR), and the relation between X and Y is linear: X = βY . The finite mixture of regression (FMR) models have been studied and used in various fields and applications (Bar-Shalom, 1978;Andrews and Currim, 2003;Bierbrauer et al., 2004), particular in the form of a mixture of "linear" regression models mainly because of advantages such as simplicity and interpretability. It is assumed that the response variable, y, can be predicted using a finite (K) number of linear regression models: There are a number of methods to estimate the parameters of FMR models. Maximum likelihood, method of moments, Gibbs sampling, and Bayesian approach are some examples (Quandt and Ramsey, 1978;Viele and Tong, 2002;Faria and Soromenho, 2010;Hurn et al., 2003), each with their own advantages and disadvantages (Melnykov and Maitra, 2010).

Mixture Models with Constraints
In some applications, however, instead of individual observations, "groups" of observations need to be clustered together. In other words, it may be desirable to force a group of observations to stay in the same cluster. Given the current context of clustering dealerships, we need all the data from a particular dealer to be clustered within the same model and not distributed across different models. Let D be the set of data points, and M be the number of all distinct groups within D. Define S as the set that holds all the group members: Each group s i , with n i observations, has to retain its members when assigned to a cluster. The goal is to assign each s i to one of K clusters. This problem is similar to what is known as "clustering with must-link constraint", which is introduced by Wagstaff and Cardie (2001) in the literature (Wagstaff et al., 2001). We refer to (Davidson and Basu, 2007) and (Basu, 2009) for further details and applications of constrained clustering.
To the best of our knowledge, all the existing algorithms have solved the problem of clustering with mustlink constraint in unsupervised/semi-supervised settings, meaning that the observations lack (or partially lack) the dependent variable. We could not find any work that addresses FMR with must-link constraints. Hence, we propose a non-parametric "competitive learning" approach for assigning predefined groups/blocks of observations to clusters during model-based regression. This method is outlined in the next section.

Mixture Model with Competitive Learning (MMCL)
In this section, we propose an algorithm that relies on "competitive learning" for FMR with group structure constraints, labeled Mixture Model with Competitive Learning (MMCL). Later sections illustrate how this algorithm can be employed to address the problem of automotive dealership clustering and performance management. We also offer guidance on parameter initialization for MMCL.

Competitive Learning
Competitive learning is a form of unsupervised learning originating in the domain of artificial neural networks, in which nodes of the network compete for the right to respond to a subset of the input data (Rumelhart et al., 1988). There are three basic elements to the standard competitive learning rule (Haykin et al., 2009): 1) A set of neurons that are all the same except for some randomly distributed synaptic weights, and therefore respond differently to a given set of input patterns; 2) A limit imposed on the 'strength' of each neuron, and 3) A mechanism that permits neurons to compete for the right to respond to a given subset of inputs, such that only one output neuron is active at a time.
In what follows, we adapt the standard competitive learning algorithm to the more general model setting of FMR (with grouping structure).

MMCL Algorithm
As noted in section 3.3, let D be the complete set of data points and M the number of distinct groups within D (e.g., data from each dealer forms an observation group). Define S as a set that holds all the groups: Each group s i (that has n i observations in it), has to retain all its members when assigned to a cluster, forming the group structure constraints. If the desired number of clusters is K, the goal is to assign each s i to one of K clusters. Ideally, the dataset is partitioned into training and testing datasets, where a subset of the data from each group is stored in the testing dataset for testing and improving model robustness.
The Psuedo code for MMCL is provided in Algorithm #1.

Algorithm 1 MMCL for FMR with Group Structure Constraints
procedure -Competitive Learning Input the dataset D Initialize AIC to a large value (e.g., 10 10 ) and to a small value (e.g., 0.001) Randomly select K groups (s i ) for initializing each of the cluster component models repeat AIC Old = AIC Learn K component models using observations assigned to K clusters for each group s i ∈ S do Make predictions for selected group using each of K models and record AIC Assign s i to cluster with the least AIC end for Assuming that the number of components, K, is known, the proposed MMCL algorithm starts by randomly selecting K groups (out of M ) for initializing the clusters, and fitting one model using observations in each group to get a function f i (x; θ) for each i ∈ 1, ..., K. In the event the individual groups are too small to learn the initial model for each cluster component, one can randomly assign multiple groups to each component for initialization. This approach is generic in a sense that the component model can be of any type such as linear regression, LASSO, multi-layer perceptron, support vector machine, etc. Next, all the groups s i ∈ S are selected one at a time, and are predicted using each of the K models. The group is assigned to one of the K cluster components that produces the best performance; this is "competition" portion of the MMCL.
For performance assessment, one can employ several criteria based on the component model, such as SSE of prediction, Bayesian Information Criteria (BIC) (Schwarz et al., 1978), or Akaike Information Criteria (AIC) (Akaike, 1998), on a held out testing dataset. Performance is evaluated on the testing dataset and not the training dataset to reduce the chance of over fitting. In our experiments with linear regression component models, AIC provided robust performance. Founded on information theory, given a collection of models for the data, AIC estimates the quality of each model relative to each of the other models (Akaike, 1998). In doing so, it deals with the trade-off between the goodness of fit of the model and the complexity of the model and is one of the most common model selection procedures available in most statistical software packages (Chaurasia and Harel, 2012). Akaike stated that modeling is not only about finding a model which describes the behavior of the observed data, but its main aim is predicated as a possible good, and the future of the process is under investigation (M., 2014). The AIC is calculated as: ( 3) where l f (θ) is the maximum value of the likelihood function of the model with parametersθ. Given that models with minimum AIC are preferred, AIC employs the term 2k to penalize complex models with more parameters.
Without loss of generality, in the rest of the manuscript, we assume that AIC is the criterion for cluster component competition within MMCL.
After assigning each group to the cluster with the best (i.e., minimum) AIC, the "Overall AIC" is calculated as the sum of the Cluster AICs: Overall AIC = K j=1 Cluster.AIC j . In the next iteration, the process is repeated and this time the K component models are learnt using all the aggregated observations (that consist of several groups) within each of K clusters; this is the "learning" portion of MMCL. The updated models will compete in the same fashion as described above, to select the cluster members that produce the lowest AIC. This process is repeated until algorithm convergence (e.g., no changes in cluster group memberships between iteration t and t+1 and/or the resulting Overall AIC) or if the maximum number of iterations has been reached. This is a common approach to stop searching in most of the meta-heuristic optimization methods (Safe et al., 2004).
Let us denote M j as the j th model (that is built using the observations in j th cluster) and A Mj ,si as the AIC value resulting from predicting the response values in group s i using model M j . The ultimate goal is to predict a label c (i) for each group s i . We can write the algorithm in the form of an optimization problem as follows: It can be seen that MMCL is a special version of the famous k-means clustering (Steinhaus, 1956), with the distance defined as the AIC of predictive models.
In the case of a linear regression component model, Eqn. 4 can be written as: where p is the number of parameters in the model and RSS Mj ,si = ni l=1 (ŷ Mj ,s il − y s il ) 2 is the residual sum of squares resulting from predicting response variables in group s i using model M j . It is easy to observe that the with a minor modification, any stopping criteria and model selection technique can be used.

Initializing MMCL Parameters
There are three main parameters that need to be selected prior to applying MMCL: 1) Number of clusters K: An existing method such as BIC, Calinsky criterion, Gap Statistics, etc. can be used. With the existence of response variables (supervised learning), it is recommended to divide the data to train and test sets and select K that yields the best result. 2) Initializing K clusters: This is a critical step as it effects both the convergence and effectiveness of the algorithm. One of the existing methods such as k-means++ (Arthur and Vassilvitskii, 2007) can be adopted to develop an algorithm that wisely selects the initial cluster groups. Inspired by k-means++, the following initialization algorithm (labeled MMCL++) is developed to select the initial groups. It tries to smartly choose the groups so that the selected groups have maximum dissimilarity. This is achieved by selecting a group s j ∈ S at random and predicting all the remaining groups using the model M j that is developed by the observations in that group. The quality of prediction for all the remaining groups using M j is evaluated, and the group s i , i = j that M j has the least power predicting it is identified as the candidate that has the maximum distance with s i . Again, different criteria such as correlation between s i and s j , RSS of M j , s i , etc. can be used for this purpose.

MMCL Validation: Synthetic Experiments
To evaluate the effectiveness of the proposed MMCL algorithm for FMR with group structure constraints, we employ Monte Carlo simulation experiments.

Experiment Setup
For the synthetic experiments, for ease of exposition, we use linear regression as the modeling technique and investigate the different parameters on the result is investigated. The experiment is conducted for the case K = 2 (number of clusters). Covariates (X) for each cluster are generated by drawing samples from a bivariate Gaussian distribution: X ∼ N (µ, Σ), with zero mean and a diagonal covariance matrix with unit variance. Referring to the Monte Carlo simulation parameters outlined in Table 1, N = 300 is the total number of observations in each cluster and S is the number of groups (blocks) in each cluster. Essentially, there will be 60 observations per group in cluster 1 and 20 observations per group in cluster 2. The response variable for each observation is generated by: y i = X i β + N oise level( ). parameter is used to control the amount of noise (uncertainty) added to the response variable y. It can be also seen as a parameter that controls the signal to noise ratio (SNR).
To study the effect of the degree of similarity between βs, the Euclidean distance (d 2 ) between β 1 and β 2 is calculated to control the level of separation for the two clusters (in the response domain).

Simulation Results
The Monte Carlo simulations are repeated 1000 times for each pair of d 2 and . Normalized Mutual Information (NMI) is used for assessing the clustering accuracy. NMI is a widely used technique in evaluating the clustering result when the true labels are available. The advantage of using NMI is that it is independent of permutation, meaning that the label switching does not affect the NMI score. It is bounded between zero and one. The closer the value to zero, the higher the indication that the cluster assignments are largely independent, while NMI close to one shows substantial agreement between the clusters. An NMI value of zero simply means that the label assignment is purely random. Figure 2a shows the average NMI for different levels of noise and d 2 , while Figure 2c shows the distribution of NMI among 1000 runs.
It is evident from Figures 2a and 2c that MMCL is able to achieve a high NMI value (about 0.8) when the is low. However, as expected, as the level of noise added to the response increases, it significantly affects the accuracy of clustering. With = 4, the clustering is almost done in a random fashion. We can also observe the effect of d 2 . As mentioned earlier, smaller d 2 means that βs are closer and more similar to each other. Figures 2a and 2c confirm that as d 2 gets smaller, it is harder to correctly cluster the observations. Figure 3a shows the average number of iterations it took for MMCL to converge under different scenarios. It can be seen that on average the algorithm converges within few iterations. Note that the number of iterations is highly affected by the stopping criteria. As mentioned earlier, the stopping criteria is based on the relative change in the overall AIC in two consecutive iterations, i.e. | AIC−AIC Old AIC Old | < . In this experiment, is set to 0.001 and the maximum number of iterations is set to 10.
To observe the impact of careful initial group selection (MMCL++), we can compare Figure 2a with 2b and 2c with 2d and note that there is a slight improvement in the result when MMCL++ is utilized, especially when the noise level is low (0.5). Figures 3a and 3b compare the average number of iterations until convergence between random initial group selection (a) and MMCL++ (b). The graphs clearly show the effectiveness of MMCL++ in reducing the number of iterations, especially under low noise level. This reduction may not seem to be notable in this simulation experiment, but it could significantly reduce the run time when dealing with large datasets where every iteration may take considerable time.

Deriving Recommendations under MMCL: Multi-objective Optimization
As noted earlier, the focus of this manuscript is to develop methods to improve the performance of individual stores by relying on a data-driven approach to internal benchmarking. In particular, the goal is to identify factors driving automotive dealership performance in comparison with "similar" dealerships and relying on optimization to derive tailored recommendations. However, as noted by (Thomas et al., 1998), more than one performance outcome usually needs to be considered because stores are responsible for multiple and sometimes conflicting performance measures (e.g., sales and profits). In addition, it is often the case that KPIs are competing for resources and cannot be adjusted independently at will (e.g., cash flow constraints might force a dealer to choose between adding more new vehicle sales staff or more service technicians but not both).
In the context of MMCL, each cluster yields specific component models depending on certain KPIs (predictors) for each of dependent variables (e.g., one model for sales and another model for profit). Performance management entails deriving recommendations (i.e., suggested levels for the predictor variables) consistent with the dynamics identified by the component models to attain balanced performance. Conflicts in performance measures will be revealed in the form of common predictor variables for multiple dependent variables and opposing coefficients (e.g., raising inventory at a store can increase sales but also increases inventory costs and might compromise profitability). The question becomes, "how to adjust the KPIs in order to jointly improve the performance of the dependent variables"? Without loss of generality, we discuss methods for addressing this problem under the assumption that there are two performance measures of interest y 1 and y 2 (e.g., vehicle sales and profit). Essentially, we need an optimization formulation to optimize both y 1 and y 2 . This is in fact a Multi-objective Optimization (MOO) problem with the possibility of multicollinearity (i.e., relations exist between KPIs or predictor variables) as well as potential relationships among dependent variables.
Let us suppose that we have a dataset with n observations and two dependent variables: y 1 and y 2 , with a set of KPI predictors {x 1 , x 2 , . . . , x p }, x i ∈ IR n . Let X denote the matrix of covariates (the design matrix), X ∈ IR n.p , where there are relations between x i s and (y 1 , y 2 ). Moreover, we have to check to see if there are relations between the predictors (x i s) themselves (i.e., multicollinearity), and consider (i.e., satisfy) these relations when formulating the MOO problem. This is the key requirement in formulating this problem. The reason is that without considering the multicollinearity among the independent variables, the model will not take any interdependencies into consideration and therefore it might produce ineffective results/recommendations.

Formulating MOO
The goal here is to find the optimal values for the actionable predictor KPI variables in order to jointly improve dependent variables y 1 and y 2 . Let x i ∈ IR n , i ∈ {1, 2, ..., p} denote a predictor KPI variable. We can setup the multi-objective optimization problem as follows: 1) Regression: Model relationships between independent (i.e., KPIs) and dependent variables Regress y 1 and y 2 on x 1 , x 2 , ..., x p to obtain f y1 (x 1 , x 2 , . . . , x p ) and f y2 (x 1 , x 2 , ..., x p ). If linear regression is employed for modeling, the result is: y 2 = f y2 (x 1 , x 2 , . . . , x p ) = γ 0 + γ 1 x 1 + γ 2 x 2 + · · · + γ p x p + y2 2) Multicollinearity: Regress each independent variable as a function of remaining independent variables Again, in the case of linear regression, we will have: x i = f (x 1 , x 2 , . . . , x r , y 1 , y 2 ) = α i0 + α i1 x 1 + α i2 x 2 + · · · + α ip x p + xi α ii = 0, i ∈ {1, 2, . . . , p} 3) Optimization: Formulation The MOO problem can be tackled using several different approaches. The classical means of solving such problems were primarily focused on scalarizing multiple objectives into a single objective (Deb and Jain, 2014). In many cases, there does not exist a single solution that simultaneously optimizes each objective. In that case, the objective functions are said to be conflicting, and there exists a (possibly infinite) number of Pareto optimal solutions. A solution is called nondominated, Pareto optimal, Pareto efficient or noninferior, if none of the objective functions can be improved in value without degrading some of the other objective values. There exist different solution philosophies and goals when setting and solving multi-objective optimization problems. Both exact and heuristic methods (e.g., genetic algorithms) have been extensively investigated in the literature. See (Deb, 2001) for a good overview on the topic.
Given that our independent and dependent variables are (mostly) continuous and we have two dependent variables (y 1 and y 2 ), there exists a Pareto curve of optimal solutions. Let us assume that the primary interest is to maximize y 1 , while satisfying an acceptable value for y 2 .
If we define t = β 0 + β 1 x 1 + β 2 x 2 + · · · + β p x p , then we arrive at the following formulation: The accuracy and effectiveness of the formulation results rely on the accuracy of the regression models.
Since the regression models cannot be assumed to be perfect, constraints (9b-9d) are designed to take the imperfection of the regression models into account. For this purpose, we allow slack for the regression models derived constraints proportional to k * σ, where σ denotes regression model standard error. Small values of k lead to strict constraints, i.e., strong agreement with regression models at the risk of recommendations that limit performance. Constraints (9e) are optional and limit the decision variables to practical bounds.

Case Study: Deriving Recommendations for Automotive Dealerships
We focus here on applying the proposed MMCL method for modeling the productivity of automotive dealerships across the U.S. for a particular OEM. The raw data mostly come from monthly financial statements of dealers to OEMs, along with dealership demographic information. The financial statements contain hundreds of input and output metrics regarding the various departments in the dealership (new vehicle sales, finance, parts/service, used department, body-shop etc), resources (staffing levels/tenure, service bays etc), inventory (new and used inventory, inventory mix, age of inventory etc), marketing (type/volume), and performance (sales, variable/fixed costs, and profits for each of the departments). To maintain confidentiality of information, we are not able to reveal the OMEM. There are 3074 dealerships with monthly data for five years (2010-1015) with 281 key performance indicators (KPIs) deemed important by the domain experts. The data is aggregated to construct a design matrix X ∈ IR (3074×60)×281 . After cleaning the data and calculating the KPIs using the raw data, we found that about 6% of the data is missing. We used a matrix completion via soft thresholding SVD technique to impute the missing data, using the "softImpute" package in R (Hastie and Mazumder, 2015). The variables are then standardized to carry a mean of zero and standard deviation of one for improved performance modeling.

Applying MMCL to the Dealership Performance Problem
Because of the large size of the dataset and specially the large number of predictors, Least Absolute Shrinkage and Selection Operator (LASSO) (Tibshirani, 1996) is used for regression modeling of both the sales as well as the profitability of each dealership, for each month. As discussed in Section 4, the main parameters should be selected before running the algorithm. In this case, the parameters are: number of clusters k, selection of k dealers for initializing the clusters, and LASSO regularizer (λ). To find the best values for these parameters, the data is split into training (first four years for each dealer), and testing (financial data from last year). The parameters are then selected using cross-validation, by evaluating the quality of the models on predicting the testing data. The parameters that produce the best result (highest R 2 value) on the testing data are selected. Once the parameters are selected, the MMCL algorithm is applied to the dataset. It is observed that in most cases, the algorithm converges in less than 15 iterations.

Results
To evaluate the effectiveness of the algorithm, the data from each region is further divided into three groups based on the expected number of vehicle registrations in each dealer's territory. We present here results for one region of the U.S. Note that there are two dependent variables (y 1 denoting profitability or return on sales and y 2 denoting new vehicle sales effectiveness). Tables 2 reports the effectiveness of the proposed MMCL method versus employing a single LASSO regression model for each registration group. As the results suggest, especially in the case of SE, the MMCL has improved the accuracy of the models in predictin g both profitability and sales effectiveness. For example, for medium registrations group, using MMCL++, we are able to achieve 26% R 2 on the testing dataset for sales effectiveness in comparison with just 11% for the case of a single model. This result also suggests that there is heterogeneity among dealers and by clustering them, one can improve the analysis and generate better recommendations to dealers for improving their performance.

Applying MOO for Dealership Performance Improvement
As noted earlier, there are two dealer performance characteristics of interest: Profitability (P : y 1 ) and Sales Effectiveness (SE: y 2 ). It is desired that both y 1 and y 2 be maximized for every dealer to improve the profitability of the dealership and satisfy the needs of the OEM in selling more new vehicles.
To apply the MOO, we first regress P (y 1 ) against X to obtain: P = f P (x 1 , . . . , x p ) = β 0 + β 1 x 1 + β 2 x 2 + · · · + β p x p + P where f P (x 1 , . . . , x p ) is the objective function that we want to maximize (ignoring the error term P ). We then regress y 2 against X to obtain the constraint that explains the relation between y 2 and x i s: Also, to guarantee that SE stays in the accepted range enforced by OEM, we add the following constraint:

SE ≥SE
Lastly, we regress each x i against other x i s, to account for the multicollinearity constraints. As explained in Section 6, to consider the fact that the regression models are imperfect, we allow a slack of k * σ for each regression model. The client also provided valid bounds for each of the predictor variables but are not reported here for confidentiality. Assuming t = f P (x 1 , . . . , x p ), the final formulation takes the following form: Figure 4: Pareto optimal frontier for a specific dealership group: (a) Blue dots report average SE and P for dealers; benchmark dealer is shown in red. (b) Zoomed-in region of Pareto optimal frontier

Generating Pareto Optimal Frontier
To construct the Pareto optimal front for each dealership group, the formulation above is solved repeatedly by changing the value ofSE to obtain the Pareto optimal points for P and SE, as reported in Figure 4a. This result identifies what is potentially possible in terms of performance for the dealers within the reference group. As expected, the frontier also reveals the trade-off between how much profit a dealer can generate (using all the potential resources such as new and used vehicle sales, service, body shop, parts, etc.) versus new vehicle sales (not as profitable these days with respect to other dealership operations such as service and used vehicle sales).

Assessing the Quality of Recommendations Derived through MOO
To further assess the effectives of the proposed methods, we compared the operational signature of a high performance dealership (shown in Figure 4a) with the recommendations derived through MOO. Radar charts are generated for some of the important KPIs (x i s) for different values ofSE (see Figure 5). The reasonable agreement between the two sets of KPI values between the high performance dealer and the derived recommendations further validate the effectiveness of our proposed algorithm. For example, it can be seen that if a dealer in the reference group wants to achieve a better SE (on average), they have to on average lower KPI3, KPI26, KPI31, and KPI15, and increase KPI7 and KPI32. The recommended KPI levels can be compared with each dealer within the group in order to show the potential strengths and deficiencies of that dealer. It provides tailored guidance for the dealership management on how to manage and operate their business in order to simultaneously please the OEM by selling new vehicles (keep SE in an acceptable range) and also increase their profits.

Conclusion
Increasing availability of data combined with improvements in computational platforms and technology is enabling more comprehensive and in-depth data analysis in the world of business. In the retail sector, individual stores need to utilize the available data to improve both their efficiency and effectiveness for survival and dominance.
Our objective was to develop a complete data-driven methodology for analyzing, processing, and modeling data from retail industry in order to understand the behavior of network of stores and provide scientific managerial guidance on how to improve and operate an individual as well as groups of stores. To achieve this goal, we addressed the problem of mixture models with group structure, and noted that this has not been addressed in the literature with the existence of a dependent variable (mixture of regressions with group structure). We proposed a solution to this problem called MMCL for segmenting stores using model based clustering. MMCL is an iterative, heuristic algorithm based on Competitive Learning to cluster groups of observations and provide a model for each group. It is a non parametric approach that can be combined with any underlying regression modeling technique. An extension to MMCL, labeled MMCL++ is introduced to smartly select the initial groups when performing the clustering. MMCL++ increases both the accuracy and speed of the algorithm. The other contribution is the development of a multi-objective optimization (MOO) formulation to improve the profitability of the retail unit while controlling other performance metrics to meet the expectations of different stakeholders. Our MOO formulation is based on statistical relations between the decision variables and is designed to handle multicollinearity among the variables. The outcome of MOO provides recommended values for the decision variables (e.g. key store performance drivers) and can be used as managerial operational guidance for the stores.
The proposed methods are validated using synthetic experiments as well as data from a real-world automotive dealership network case study for a leading global automotive manufacturer. The results of both synthetic and real-world experiments proved the accuracy and effectiveness of our methodologies. We believe that this research is a good starting point for developing an intensive and complete process for benchmarking and managing the performance of retail stores.
There are several avenues for future research. We can extend the proposed methods to embrace a maximum likelihood approach to FMR with must-link constraints and also allow joint determination of the optimal number of clusters. Also, it can be combined with recommendation process, meaning to judge the models based on the quality of recommendations. This can be an online reinforcement learning framework that assesses the result of recommendation and utilizes the result and feedback to adjust and improve the models. As for MOO, the formulation can be extended to form a pure multi-objective optimization formulation that jointly optimizes both objective functions (rather than maximizing a single objective while imposing constraints on the performance of the remaining objectives). The proposed methods should also be further validated using data from other retail settings/applications.