A Regression-constrained Optimization Approach to Estimating Suppressed Information using Time-series Data

Abstract

Most regional economic databases (e.g., US Economic Census and County Business Patterns [CBP]) have some employment records suppressed and then represented as ranges, in order to guarantee the confidentiality of the data. This article incorporates the implicit temporal relationships between annual employment data over several years into an optimization model designed to estimate suppressed records. This model minimizes (1) the sum of the deviations between the estimates and target values within the corresponding ranges and (2) the sum of the deviations between the estimates and an employment trend curve endogenously determined through absolute-value regression. The 1999–2006 CBP data for Arizona are used to test the model. Two decision-theoretic criteria (Pareto frontier and concordance–discordance analysis) are used to analyze the results, pointing to a specific set of parameters yielding the best estimates.

Keywords

data suppression county business patterns employment trend curve absolute-value regression optimization

Introduction

Regional economic databases, such as the Economic Census (EC) and County Business Patterns (CBP) in the United States and the Labour Force Survey in the United Kingdom, provide essential inputs to economic analyses. These databases vary in terms of their geographical detail, down to the state (EC-mining), county (EC and CBP), and place with 2,500+ people (EC-manufacturing), and in terms of their publication frequency (annual for CBP, every five years for EC). They also vary in terms of content and economic sector disaggregation, with employment as the most common data. An additional feature of these databases is that a significant number of data records are suppressed and represented by interval flags.

The issue of economic data suppression has generated a limited literature on methods to estimate suppressed records (e.g., Gardocki and Baj 1985; Sechrist 1986; Kreahling, Smith, and Frumento 1996; Isserman and Westervelt 2006; Zhang and Guldmann 2009). These methods focus on annual cross sections of data. However, all economic databases are longitudinal, and a given record, defined by sector and geography, may be suppressed in some years, but not in others. Accounting for the available time series may improve the estimation. This article develops and implements an optimization methodology for estimating suppressed employment data in a longitudinal database, extending the approach developed by Zhang and Guldmann (2009). The 1999–2006 CBP data for Arizona are used as a case study.

The remainder of the article is organized as follows. The second section presents an overview of CBP data and some of their applications, clarifies and illustrates the data suppression problem, and reviews existing estimation methods. The third section presents the estimation methodology, including the narrowing of the flag intervals and the optimization model. The fourth section presents an application of the methodology to Arizona CBP data, including a search for combinations of model parameters that provide better estimates. The fifth section consists of a discussion of the results and caveats. The sixth section concludes the article and outline areas for further research.

Statement of the Problem

The US Census Bureau (CB) began to publish CBP data annually in 1964, and complete data sets are available for download from the Census website (http://www.census.gov/econ/cbp/download/index.htm) back to 1986. These data sets provide employment, payroll, and total number of establishment, by industry and by geography. Numbers of establishments by employment range are also available. The database includes a national file, a state file for each of the fifty states, and a county file for each of the 3,143 counties. The North American Industry Classification System (NAICS) has been used since 1998. It has a hierarchical structure. The grand total in any geography is disaggregated into twenty-one 2-digit sectors, each of which consists of a variable number of 3-digit sectors. This nested disaggregation continues down to the 6-digit level.

CBP data have been extensively used by academic researchers, governmental agencies, and private firms. For instance, CBP data have been used by Isserman (1977) and Guimaraes, Figueiredo, and Woodward (2009) to compute location quotients in order to assess regional economic impacts and to measure industrial agglomeration, and by Stevens and Moore (1980), Miller and Blair (1985), and Henderson (1997) for shift-share, regional input–output, and econometric analyses. They have also been used to evaluate government economic policies, such as enterprise zone programs (Dowall 1996; Moore 2003). Pagoulatos (2004) has combined CBP with other data to assess the environmental impacts of economic development.

To guarantee that individual establishment records cannot be retrieved, records have to be suppressed. Table 1 presents the distribution of suppressed records across the economic and geographic hierarchies of the 2006 CBP. More records are suppressed at higher disaggregation levels. Only 4.0 percent and 31.2 percent of two-digit records are suppressed in the state and county files, respectively, as compared to 41.2 percent and 73.9 percent at the six-digit level. More records are suppressed in smaller geographies. At the national level, twelve records (0.6 percent) are suppressed (of the 2,148 records). At the state level, 33.4 percent of the 99,082 records are suppressed. As the county level, 68.1 percent of the 2,203,501 records are suppressed. The CB provides a flag (A, B, C, E, … , or M) for each suppressed record. Each flag is associated with an employment interval. For instance, flag A corresponds to the interval [0, 19].

Table 1.

Numbers of Suppressed Records and Rates of Suppression in the 2006 County Business Patterns (CBP).

Sector Level	US CBP			State CBP			County CBP
Sector Level	# of Flagged Records	Total # of Records	Percentage Flagged	# of Flagged Records	Total # of Records	Percentage Flagged	# of Flagged Records	Total # of Records	Percentage Flagged
Grand total	—	1	—	—	51	—	16	3,188	0.50
2-digit	—	20	—	41	1,020	4.02	18,423	59,016	31.22
3-digit	—	87	—	640	4,400	14.55	99,757	187,023	53.34
4-digit	—	294	—	3,368	14,531	23.18	280,631	433,474	64.74
5-digit	2	654	0.31	9,418	31,413	29.98	495,893	701,533	70.69
6-digit	10	1,092	0.92	19,616	47,667	41.15	605,054	819,267	73.85
Total	12	2,148	0.56	33,083	99,082	33.39	1,499,774	2,203,501	68.06

The CB recommends that the midpoints of the flag intervals be used in place of the suppressed employment data. Table 2 illustrates the suppression issue when the midpoint approach is applied to the 2000 CBP data for Arizona (sixteen counties and twenty-one 2-digit sectors). The actual and estimated employments are presented for only four counties (Coconino, Graham, Greenlee, and Pinal) and for sector and county totals. A gray cell indicates a suppressed record. Note that county and sector totals are not suppressed. At least two records must be suppressed in every row and column to avoid possible disclosure of information. For instance, in Coconino county, if the data of only one sector (eleven or ninety-nine) is suppressed, this data can be retrieved by subtracting the sum of all the known data from the county total (38,917). Similarly, if the data for educational services is suppressed for only one county (Graham or Greenlee), it can be retrieved by subtracting the sum of the known data from the sector total (24,486). The midpoint approach creates inconsistencies in county and sector totals, as illustrated in the last row and column in Table 2, where the differences between the computed and actual totals are presented. This inconsistency is more acute in counties with large numbers of suppressed records, for instance, Greenlee (529 or 19.1 percent).

Table 2.

Illustration of Data Suppression and the Midpoint Approach—The Case of 2000 Arizona Two-digit Data.

North American Industry Classification System Sector	County				Sector Total
North American Industry Classification System Sector	Coconino	Graham	Greenlee	Pinal	Computed	Actual	Difference
Forestry, fishing, hunting, and agriculture support (11)	10	5		193	1,963	1,931	32
Mining (21)	117	60	1,750	750	9,319	9,956	−637
Utilities (22)	234	60	60	175	9,942	9,923	19
Construction (23)	2,557	349	60	1,395	167,671	167,401	270
Manufacturing (31–33)	2,345	256	10	3,747	200,754	200,859	−105
Wholesale trade (42)	1,029	203	48	656	86,222	86,222	0
Retail trade (44–45)	7,408	1,355	142	4,607	255,001	255,001	0
Transportation and warehousing (48)	1,605	60	60	348	70,397	70,337	60
Information (51)	595	90	60	311	47,452	47,411	41
Finance and insurance (52)	777	112	10	671	111,928	111,928	0
Real estate and rental and Leasing (53)	1,019	62	10	464	41,286	41,265	21
Professional, scientific, and technical services (54)	1,005	175	10	715	113,809	113,687	122
Management of companies and enterprises (55)	402	10	60	98	45,429	45,452	−23
Admin, support, waste management, remediation services (56)	1,376	126	60	2,354	200,430	200,421	9
Educational services (61)	581	60	10	134	24,502	24,486	16
Health care and social assistance (62)	4,975	970	175	2,886	199,974	199,894	80
Arts, entertainment, and recreation (71)	1,167	10	10	652	34,341	34,394	−53
Accommodation and food services (72)	9,454	681	128	3,774	202,756	202,756	0
Other services (except public administration) (81)	1,691	414	34	1,198	78,730	78,730	0
Auxiliaries (except corporate, subsidiary, and regional management) (95)	446	10		60	14,582	14,541	41
Unclassified establishments (99)	175	1	10	60	2,885	2,758	127
County total
Computed	38,968	5,069	2,707	25,248
Actual	38,917	4,939	3,236	25,131
Difference	51	130	−529	117

Several other methods have been proposed to estimate suppressed information. Gardocki and Baj (1985) take advantage of the lower suppression rates in large geographies, and use state or national data to estimate suppressed county employment. Weddleton and Olson (1988) test regression and ratio approaches to estimate undisclosed CBP earnings, using available disclosed data. The regressions estimate the relations between industry earnings and number of establishments, and the ratio approaches directly use the earning–establishment ratio, which may be available for larger geographies and is assumed to be stable within industries. The above methods do not consider internal consistency within the geographic and economic hierarchies. The method of indirect standardization (Sechrist 1986) guarantees internal consistency within the economic hierarchy (the sum of disclosed and estimated employment matches upper-level disclosed totals). Kreahling, Smith, and Frumento (1996) propose an approach similar to the indirect standardization but with improved accuracy. They consider both flag range and establishment size information (ESI) to create smaller intervals. Ellison and Glaeser (1997) further consider geographic consistency and use an ad hoc sequential procedure to estimate suppressed information. Isserman and Westervelt (2006) propose an ad hoc approach to provide employment ranges internally consistent with the available information. They claim that these are the smallest possible ranges. Zhang and Guldmann (2009) develop an optimization model that considers flag and establishment information and other constraints within the economic and geographic hierarchies, and generates employment estimates that minimize the sum of weighted squared deviations between the estimates and a set of target values. The approach was tested with 2000 Arizona CBP data.

Method

The optimization model presented in this section extends the static model proposed by Zhang and Guldmann (2009), by incorporating employment trend information and tighter bounds for the suppressed data. Prior to presenting this model, the approach to first tightening the bounds is discussed, following in part the approach of Isserman and Westervelt (2006). Second, actual employment trends are analyzed by industry and county to assess the potential value of incorporating temporal information into the model.

Narrowing the Flagged Employment Intervals

Several pieces of information are used to generate smaller employment ranges, including suppression flag, number of establishments by employment range, and hierarchical relationships across NAICS levels. This process is illustrated with Arizona data.

The narrowing process starts with the flag intervals. There are 2,832 two-digit records in the sixteen counties of Arizona over 1998–2006, 778 of which are suppressed, mostly flagged as A or B. The lower and upper bounds of a flag interval for sector i in county c in year t, as defined by the CB, are noted as FMIN_cit and FMAX_cit . The column “Interval 1” in Table 3 presents the interval size for each flag.

Table 3.

Descriptive Statistics for the Size of Flagged and Narrowed Intervals—The Case of Two-digit Arizona Data over 1998–2006.

	Suppressed Records		Interval 1^a	Interval 2^a				Interval 3^a				Improvements^b (%)
Flag	Number	Percentage	Interval 1^a	M	Min.	Max.	SD	M	Min.	Max.	SD	Improvements^b (%)
A (0–19)	282	36.25	19	14.0	4	19	5.9	13.9	4	19	5.9	26.84
B (20–99)	308	39.59	79	52.1	7	79	22.9	47	3	79	22.7	40.51
C (100–249)	113	14.52	149	113.3	3	149	42.6	86.8	3	149	49.2	41.74
E (250–499)	30	3.86	249	221.6	39	249	56.7	168.8	16	249	76.9	32.21
F (500–999)	24	3.08	499	447.8	208	499	88.7	280.8	69	479	133.8	43.73
G (1000–2499)	16	2.06	1,499	1,171.6	709	1,499	364.0	992.4	40	1,499	527.0	33.80
H (2500–4999)	4	0.51	2,499	2,302	1,711	2,499	394.0	2,279	1,619	2,499	440.0	8.80
L (50000–99999)	1	0.13	49,999	49,999	49,999	49,999		284	284	284		99.43
Total	778	100.00

Note: ^aInterval 1 is the flag interval itself. Interval 2 accounts for additional two-digit level establishment information. Interval 3 further considers information across hierarchical levels.

^bComputed as (mean interval 3 size/interval 1 size)/interval 1 size, that is, percentage reduction for the flagged interval.

ESI, which is never suppressed, allows for narrowing the flag interval. For instance, sector 21 (mining, quarrying, and oil and gas extraction) in county 021 (Pinal) is flagged as F [500, 999] in 2000. However, the ESI indicates a range of [347, 708].¹ The flag interval [500, 999] can then be narrowed down to [500,708]. For 23 of the 51,883 Arizona records flagged over 1998–2006 (across all NAICS levels), the lower bound generated by the ESI bound exceeds the upper bound of the flag. Census documentation indicates that the flag interval is based on mid-March employment, while establishment counts include any establishment with at least one employee during the year (http://www.census.gov/econ/cbp/definitions.htm). Therefore, employment inconsistencies may occur between flag and establishment intervals. The lower ESI should therefore be considered unreliable, and only the upper ESI bound is used to narrow down the flag interval. In Table 1, the “Interval 2” columns provide descriptive statistics for the sizes of these new intervals.

Finally, information at more disaggregated NAICS levels can be used to further narrow Interval 2. Appendix A illustrates a forward approach to find a narrower interval for sector 21 (mining) in county 21 (Pinal, Arizona) in 2000, considering flag intervals and ESI across all NAICS levels, starting at the six-digit level. In Table 3, “Interval 3” columns provide descriptive statistics for this new set of intervals. For every type of flag, the average interval size has been further reduced. Overall, reductions in interval size range from 27 percent to 44 percent, except for flag H.

Employment Temporal Patterns

Does sector/county employment display a stable trend? To test this hypothesis, several regression models, running from linear to quadratic to cubic, are estimated with actual (non-suppressed) Arizona employment data over 1998–2006. More details are presented in Appendix B. Overall, these regressions suggest that most sectors display a stable employment pattern over time. On average, the linear model accounts for 50.9 percent of employment variations, the quadratic model for 65.6 percent, and the cubic model for 77.3 percent. These patterns support the incorporation of employment time series in the optimization model.

Regression-constrained Optimization Model

Definitions

The indices are:

county;

two-digit NAICS sector;

time (=1 if 1998, =2 if 1999, … , and =9 if 2006).

The basic decision variable is:

X_cit :

employment estimate for two-digit sector i in county c at time t.

The model parameters are:

set of combinations (c, i, t) for which employment is flagged;

E_cit , E_ct , E_it :

actual employment, (c, i, t), (c, t), (i, t) ∉ M;

MIN_cit :

minimum employment (c, i, t) ∊ M;

MAX_cit :

maximum employment (c, i, t) ∊ M;

CMIN_ct :

minimum county total employment (c, t) ∊ M;

CMAX_ct :

maximum county total employment (c, t) ∊ M;

SMIN_it :

minimum sectoral total employment (i, t) ∊ M;

SMAX_it :

maximum sectoral total employment (i, t) ∊ M.

The narrow interval [MIN_cit , MAX_cit ] for a two-digit suppressed record is discussed in Narrowing the Flagged Employment Intervals section. As county totals and sector totals (across the whole state) may also be suppressed, [CMIN_ct , CMAX_ct ] and [SMIN_it , SMAX_it ] represent the narrow intervals for the corresponding records. When county or sector totals are not suppressed, the above intervals collapse onto the actual total employments, E_ct and E_it .

Basic Constraints

The following constraints enforce data consistency in terms of (1) suppressed data intervals, (2) county total employment, and (3) state total employment:

X_{c i t} = 0 (c, i, t) \notin M

M I N_{c i t} \leq X_{c i t} \leq M A X_{c i t} (c, i, t) \in M

\sum_{i} (E_{c i t} + X_{c i t}) = E_{c t} (c, t) \notin M

C M I N_{c t} \leq \sum_{i} (E_{c i t} + X_{c i t}) \leq C M A X_{c t} (c, t) \in M

\sum_{c} (E_{c i t} + X_{c i t}) = E_{i t} (i, t) \notin M

S M I N_{i t} \leq \sum_{c} (E_{c i t} + X_{c i t}) \leq S M A X_{i t} (i, t) \in M

Constraint (1) forces the estimation variable to be zero when the actual employment is available. Constraint (2) guarantees that any estimate must lie within its interval. Constraint (3) guarantees consistency between the sum of sectoral employments in a county and that county total employment. Constraint (4) replaces constraint (3) when the total county employment is flagged. Constraint (5) guarantees consistency between the sum of sectoral employments across all counties and that sector total employment at the state level. Constraint (6) replaces constraint (5) when the total state employment is flagged.

Estimation Objective Functions

The first objective function, F ₁, is based on Zhang and Guldmann (2009) and minimizes the sum of the absolute values of the deviations between an estimate and a target location:

F_{1} = \sum_{(c, i, t) \in M} |X_{c i t} - (M I N_{c i t} + l (M A X_{c i t} - M I N_{c i t}))|,

where l is a location parameter varying between 0 and 1. For instance, l = 0, 0.5, or 1 indicates that the target location is the lower bound, the midpoint, or the upper bound of [MIN_cit , MAX_cit ].

In order to illustrate the second objective function, F ₂, that involves employment trends, consider the hypothetical employment data pattern in Figure 1, corresponding to a specific county/sector (c, i). The eight time periods (e.g., years) include four periods ( $t \in T1$ ) with known employment, and four periods ( $t \in T2$ ) with suppressed employment, which is represented by its interval, as determined in Narrowing the Flagged Employment Intervals section. A first and natural approach would be to estimate the E_ci (t) curve using the actual employment data ( $t \in T1$ ) and ordinary least squares (OLS), and then compute employment for $t \in T2$ with this estimated function. In the case of Figure 1, a linear curve E(t) appears to closely match the actual employment data, and the intersections of E(t) with the suppressed data intervals are all located within these intervals, and are therefore acceptable as estimates. However, such feasible intersections are not at all guaranteed, and preliminary analyses of Arizona data show that there are many situations where these intersections would lie outside the intervals, and thus would not be feasible. Further, even if the intersection lies in the interval, it may be located at one or the other extreme of the interval (e.g., t = 6 in Figure 1), away from a more likely location within the interval (say, the midpoint).

Figure 1.

Hypothetical employment regression curve.

In order to avoid the above problem, the OLS estimation could incorporate both the actual data (E_cit ) and the estimation variables (X_cit ), which would be subject to constraints (1)–(6). However, this approach would then completely ignore the deviation function F ₁, thus creating a single-objective optimization model. It would take the form of a nonlinear program, due to the nonlinearity of the objective function, with the possibility of generating only locally optimal solutions. The alternative is to use an absolute-value deviation regression, which can be easily linearized. Consider the third order polynomial regression function:

{\hat{E}}_{c i t} = a_{c i} + b_{c i} t + d_{c i} t^{2} + g_{c i} t^{3},

where the regression coefficients ( $a_{c i}, b_{c i}, d_{c i}, g_{c i}$ ) are to be endogenously determined by the model. The aggregate deviation from the regression line is given by:

F_{2} = \sum_{c, i, t} |X_{c i t} + E_{c i t} - (a_{c i} + b_{c i} t + d_{c i} t^{2} + g_{c i} t^{3})| .

In equation (9), either X_cit = 0 or E_cit = 0, and the unknowns are X_cit , a_ci , b_c , d_ci , and g_ci . The form of the regression curve is endogenously determined. For instance, if a linear curve provides the best fit for both the estimated data (X_cit ) and the actual data (E_cit ), then the output will include d_ci = 0 and g_ci = 0.

Linearizing the Objective Function

In order to convert the objective functions (7) and (9) into linear forms, two sets of deviation variables are introduced: (1) PD_cit and ND_cit , to measure the deviation between an employment estimate and its location target and (2) PR_cit and NR_cit , to measure the deviation between an employment estimate or the actual employment and the regression curve. The two objective functions are reformulated as:

F_{1} = \sum_{(c, i, t) \in S} (P D_{c i t} + N D_{c i t}) .

F_{2} = \sum_{(c, i, t)} (P R_{c i t} + N R_{c i t}) .

The deviation constraints are as follows:

P D_{c i t} - N D_{c i t} = X_{c i t} - (M I N_{c i t} + l (M A X_{c i t} - M I N_{c i t})) (c, i, t) \in M .

P R_{c i t} - N R_{c i t} = X_{c i t} + E_{c i t} - (a_{c i} + b_{c i} t_{t} + d_{c i} t_{t}^{2} + g_{c i} t_{t}^{3}) All (c, i, t) .

P D_{c i t} \geq 0, N D_{c i t} \geq 0, P R_{c i t} \geq 0, N R_{c i t} \geq 0.

Minimizing either F ₁ or F ₂ is likely to lead to contrasted solutions, as F ₁ pulls the solution to a selected target within the intervals and F ₂ pulls it toward the regression curves. As it is impossible to determine, a priori, which criterion is more important, the two functions are combined into a single function representing a weighted convex combination of F ₁ and F ₂, using the weights (1 − w) and w, with $w \in [0, 1]$ :

F = (1 - w) F_{1} + w F_{2} .

w = 1 indicates that only temporal trends are considered, and w = 0 completely ignores the temporal trends. The latter implies independent optimization of each cross section, as in the Zhang-Guldmann (2009) model. The final model takes the form of a linear program, minimizing F (equation 15) subject to constraints (1)–(6) and (12)–(14).

Application to Arizona

Overview

The methodology presented in Method section is applied using CBP county data for Arizona. Because an inconsistency was uncovered in the source data between county total and two-digit level information ² in the 1998 CBP, all further analyses use only the 1999–2006 data. After a short graphical illustration of the estimation results (Estimation Illustration section), the approach to determining the “best” parameters l and w is outlined in Choosing the Parameters l and w section. The notations L and W are introduced as aliases for l and w, the location and weight parameters.

The 1999–2006 suppressed employment records are estimated for ranges of values for L and W, creating complete data sets. As the CB does not provide any information on its data suppression algorithm, it is impossible to model the exact data generation process that creates the suppressed records. The alternative is to use a Monte Carlo simulation (Kennedy 2003) to replicate the observed suppression patterns by generating several suppressed data sets by artificially suppressing records in the complete data sets. The optimization model then derives employment estimates for different location and weight combinations (l, w). The quality of these estimates is then assessed using different criteria, and dominant parameter combinations (l, w) are proposed, that would generate “best” estimates.

Estimation Illustration

To illustrate the applicability of the model, the parameter combinations (L = W = 0.5) and (L = 0.5, W = 0) are first used to estimate the originally suppressed two-digit Arizona records over 1999–2006. Figure 2 illustrates the improvement due to introducing employment trends, focusing on three sectors/counties. The upper-block diagrams (A) present results when L = 0.5 and W = 0, thus only considering the objective F ₁. The lower-block diagrams (B) present results when L = W = 0.5, thus weighting F ₁ and F ₂ equally. The continuous black curve is derived by OLS with actual employment data, while the absolute-deviation regression curve derived by the optimization appears only on the lower blocks in dashed red line. In case 1, the optimization-derived curve closely matches the OLS curve, and the unique estimated value (t = 5) shifts upward from A to B to be closer to the dashed curve, at the interval boundary. In case 2, with two suppressed records, the match between the two curves is less close, because the optimization-derived curve is adjusted to capture the midpoint of the interval in year t = 1. The OLS curve only accounts for the actual employment data and thus ignores the information provided by these intervals. Finally, in case 3, the OLS curve cannot be estimated because there is only one actual employment value. However, this is not an issue for the optimization-derived curve. As expected, the estimates are either at the midpoint or very close to the midpoint when W = 0 (A). However, when W = 0.5, some estimates are pulled away from the midpoints and toward the optimization-derived regression curve.

Figure 2.

Estimation illustration—Arizona data.

Choosing the Parameters l and w

Are there specific (l, w) combinations that will generate better estimates? Figure 3 presents the flowchart of the process used for delineating such combinations. Any (L, W) combination generates a set of estimates for the originally suppressed records, and therefore creates a complete data set, which is then randomly suppressed, generating ten suppressed data sets (D1–D10). For each of these data sets, any combination (l, w) generates a set of estimates, the quality of which is then assessed by comparing them with the “known” values in the complete data set, using various criteria. Several combinations (l, w) are tested, and the best (or dominant) ones are selected. However, these “best” (l, w) combinations may depend on the choice of (L, W). Hence, the above process is repeated for different (L, W) combinations. The best (l, w) combinations are compared across the range of (L, W) combinations, in order to identify a core of combinations that produce better estimation results, no matter how the complete data set is created. If this core exists, one could be confident that its combinations would achieve better estimation results.

Figure 3.

Assessing combinations (l, w).

Generating Suppressed Data Sets

Twenty-five combinations of $(L, W) \in (0. 1, 0. 3, 0. 5, 0. 7, 0. 9)$ are used. Let D_LW be the complete data set generated by a given (L, W) combination. For each D_LW , a stratified random sampling is implemented to create 10 suppressed data sets, thus a total of 250 data sets. Let D_LWv be such a data set (v = 1→10). These samples mimic the original suppression pattern in the source CBP file (Table 4). Note that county and state totals not suppressed in the source data were not considered for suppression in this process. State employment totals were never suppressed over 1999–2006, nor were county totals over 1999–2005. In 2006, two county totals (Greenlee and La Paz) were suppressed, and this suppression was replicated in all samples.

Table 4.

Distribution of Two-digit Suppressed Records in Arizona, 1999–2006.

	1999	2000	2001	2002	2003	2004	2005	2006
A	25	27	32	35	38	34	32	30
B	42	36	39	30	24	30	29	36
C	5	13	19	11	10	14	13	16
E	4	2	3	4	1	2	4	5
F	2	3	1	1	3	3	4	3
G		2	3	2	2	2	2	2
H	1							1
L								1
Total	79	83	97	83	78	85	84	94

The next step is to apply the optimization model, with different (l, w) combinations, to the 250 data sets D_LWv . Preliminary tests were carried out for data set $D_{L = W = 0.5}$ while considering eleven possible values (0, 0.1, 0.2, … , 1) for both l and w. The results showed that l values in the range [0.2–0.6] and w values in the range [0.4–0.8] produced better estimates. Therefore, additional l values [0.25, 0.35, 0.45, and 0.55] and w values [0.45, 0.55, 0.65, and 0.75] were considered. As a result, 225 combinations (fifteen l’s and fifteen w’s) were considered for each D_LWv data set. Thus, 225 optimization results are obtained for each D_LWv , and a total of 56,250 (225 × 250) sets of employment data are produced by the model.

Evaluation of the Estimated Data Sets

The estimation results are compared to the known values of data set D_LW to assess the effects of the different parameter combinations. Let D_LWvlw be the data set obtained by applying the optimization model with parameters (l, w) to estimate the suppressed data in data set D_LWv . This section presents the assessment of these estimates. A goodness-of-fit criterion is first defined, that helps identify preferred combination (l, w) across the (L, W) combinations. Next, Pareto frontier (PF) and concordance–discordance (CD) criteria are proposed and computed. A summary comparison across these criteria identifies core (l, w) values.

Goodness of Fit

The goodness-of-fit criterion (S) represents the sum of the squared differences between estimates and “known” values. In the case of location and weight parameters (l, w), sample v, and complete data set D_LW , the criterion is given by:

S_{L W v l w} = \sum_{c, i, t \in N_{L W v}} {(E s t_{L W v l w c i t} - K n o_{L W c i t})}^{2},

where $E s t_{L W v l w c i t}$ is the employment estimate for county c, sector i, time t, and data set D_LWv , $K n o_{L W c i t}$ is the corresponding “known” employment in data set D_LW , and N_LWv is the set of suppressed records in D_LWv . For a given data set, the lower S_LWvlw , the better the estimates.

The gray cells in Figure 4 represent the top ten best combinations (l, w) for any specific data set D_LW , using the average S criterion. For instance, for the data set D_LW estimated with L = W = 0.5, the combination (l = 0.3, w = 0.7) yields the lowest average S over the ten samples. Overall, the best combinations are generated by using $l \in [0. 2, 0. 55]$ and $w \in [0. 6, 0. 8]$ . The last column and last row of each D_LW block (AL columns and AW rows) indicate the average S value for each parameter. In case of parameter w, w = 0.8 yields the lowest average S for eighteen of the twenty-five data sets. For the other data sets, 0.7 and 0.75 yield the lowest values. In the case of parameter l, 0.2 is the best choice for fourteen data sets, 0.25 for six data sets, 0.3 for three data sets, and 0.35 for two data sets. When the weight w shifts from 0 (considering only interval deviation) to 0.1 (considering temporal information), the average S drops significantly in all the twenty-five data sets, by 41.7 percent to 48.2 percent, with an average of 43.9 percent. When w = 0, the model focuses on each year separately. The significantly large improvement achieved by only slightly considering employment time series provides a strong validation for the new model developed in this article.

Figure 4.

Best combinations of location and regression parameters (l and w) based on the average goodness of fit.

PF Analysis

All the 225 S values for each sample D_LWv are next normalized, using the lowest S value. Table 5 presents descriptive statistics for the original and normalized S (NS) in the case L = W = 0.5. For instance, in the case of sample 1, all the 225 S values are divided by the lowest one (144,489). More generally, a NS is formulated as:

N S_{L W v l w} = S_{L W v l w} / \underset{(l, w)}{\min (} S_{L W v l w}) .

Table 5.

Descriptive Statistics for Original and Normalized S for L = W = 0.5.

	Sample 1	Sample 2	Sample 3	Sample 4	Sample 5	Sample 6	Sample 7	Sample 8	Sample 9	Sample 10
Original S
M	267,681	208,749	243,722	301,490	195,040	274,079	297,156	254,213	213,450	258,193
Min.	144,489	114,576	154,234	205,910	110,990	145,353	187,894	153,832	136,149	136,869
Max.	770,551	590,930	697,590	623,705	554,373	888,488	856,134	605,995	686,904	804,361
SD	86,784	71,827	69,431	55,645	71,574	84,396	77,973	66,894	77,583	85,771
Normalized S (NS)
M	1.85	1.82	1.58	1.46	1.76	1.89	1.58	1.65	1.57	1.89
Min.	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
Max.	5.33	5.16	4.52	3.03	4.99	6.11	4.56	3.94	5.05	5.88
SD	0.60	0.63	0.45	0.27	0.64	0.58	0.41	0.43	0.57	0.63

The NS values enable direct comparisons of the performance of parameter combinations across samples. The closer NS is to 1 the better the choice, and the best choice corresponds to a value of 1. The best (l, w) combinations are expected to consistently yield close-to-one NS values for all the ten samples of a data set D_LW .

For each combination (L, W, l, w), there are ten NS values for the ten samples v. The mean value

M N S_{L W l w} = \frac{1}{10} \sum_{v = 1}^{10} N S_{L W v l w},

may be viewed as a good indicator of the performance of the combination (L, W, l, w) and the closer to one the better. However, the spread of these values, as measured by the standard deviation

S N S_{L W l w} = {(\sum_{v = 1}^{10} {[N S_{L W v l w} - M N S_{L W l w}]}^{2} / 10)}^{1/2}

should also be considered. A truly optimal combination (L, W, l, w) would minimize both the mean and the standard deviation.

A PF analysis is conducted to assess the trade-off between MNS and SNS across the ten samples for (L, W, l, w) combinations. Figure 5 presents the scatter plot of MNS and SNS for the 225 data sets created with L = 0.1 and W = 0.9. A PF emerges, pointing to three best combinations of parameters, [l = 0.25, w = 0.8], [l = 0.35, w = 0.8], and [l = 0.45, w = 0. 8]. They do not dominate each other, that is, when MNS decreases, SNS increases, and vice versa. Further, as shown in Figure 5, there is a clear cluster of points when w ≠ 0. The estimates produced with w = 0, that is, not considering employment time trends, are clearly inferior to the other estimates, confirming the results in Goodness of Fit section. Results produced with w = 0.5 and w = 1 are not optimal but acceptable. This pattern holds for the other twenty-four combinations of L and W.

Figure 5.

Pareto frontier analysis for L = 0.1, W = 0.9.

Table 6 lists the points (l, w) on the PFs of all the twenty-five data sets D_LW , together with the MNS and SNS values. There are five (L, W) combinations where the PF is reduced to one point, which corresponds to the minimizing of both MNS and SNS. The PF consists of three points for five (L, W) combinations, two points for fifteen combination, and one point (where both MNS and SNS are optimized simultaneously) for five combinations. The results point to the following ranges: $l \in [0. 25 - 0. 45]$ and $w \in [0.7 5 - 0.8]$ . These ranges are consistent with those delineated on Figure 4, but tighter. While there is a slight correlation between L and l (0.38), there is no correlation between W and w (0.19). Overall, the delineated parameter ranges are fairly narrow and insensitive to the way the data set D_LW is created. A model user selecting (l, w) values within these ranges to estimate the suppressed values for Arizona should then be reassured that the estimation results will be close to the (unknown) true values.

Table 6.

Summary of Pareto Frontier Analysis.

Data set		Point A				Point B				Point C
L	W	l	w	MNS	SNS	l	w	MNS	SNS	l	w	MNS	SNS
0.1	0.1	0.35	0.8	1.143	0.056	0.45	0.8	1.128	0.089
	0.3	0.35	0.8	1.156	0.060	0.45	0.8	1.120	0.077
	0.5	0.35	0.8	1.154	0.046	0.45	0.8	1.119	0.080
	0.7	0.4	0.8	1.149	0.077	0.45	0.8	1.123	0.099
	0.9	0.25	0.8	1.155	0.127	0.35	0.8	1.166	0.054	0.45	0.8	1.156	0.111
0.3	0.1	0.35	0.8	1.187	0.098
	0.3	0.3	0.75	1.176	0.131	0.35	0.8	1.178	0.094
	0.5	0.25	0.75	1.170	0.142	0.3	0.75	1.169	0.153	0.35	0.8	1.176	0.110
	0.7	0.3	0.8	1.151	0.080	0.35	0.8	1.138	0.081
	0.9	0.25	0.75	1.137	0.121	0.25	0.8	1.138	0.101	0.35	0.8	1.158	0.075
0.5	0.1	0.3	0.7	1.124	0.070	0.35	0.75	1.158	0.073
	0.3	0.3	0.7	1.124	0.070
	0.5	0.3	0.75	1.135	0.087	0.35	0.75	1.162	0.043
	0.7	0.35	0.75	1.131	0.052
	0.9	0.35	0.75	1.157	0.051	0.4	0.8	1.100	0.107	0.45	0.8	1.108	0.091
0.7	0.1	0.35	0.8	1.200	0.103	0.45	0.8	1.196	0.138
	0.3	0.35	0.8	1.192	0.088	0.4	0.8	1.177	0.127
	0.5	0.35	0.8	1.178	0.101	0.4	0.8	1.172	0.118
	0.7	0.35	0.8	1.182	0.129	0.45	0.8	1.191	0.113
	0.9	0.35	0.8	1.151	0.054	0.45	0.8	1.136	0.105
0.9	0.1	0.45	0.8	1.129	0.068
	0.3	0.45	0.8	1.138	0.081	0.5	0.8	1.142	0.081
	0.5	0.45	0.8	1.130	0.068
	0.7	0.45	0.8	1.147	0.073	0.4	0.8	1.173	0.071	0.5	0.8	1.143	0.085
	0.9	0.35	0.8	1.130	0.098	0.45	0.8	1.137	0.097

CD Analysis

Figure 4 points to the ranges [0.2–0.55] for l, and [0.6–0.8] for w, or a total of forty combinations for every data set (L, W). For each (l, w) combination, there are 250 S values. These forty combinations are first sorted by average S value from the smallest to the largest, denoted Y ₁, Y ₂, … , Y ₄₀. Given that these Y’s are not normally distributed, the Wilcoxon Signed-Rank test is used for the following hypothesis:

H_{0} : Y_{g} \geq Y_{j} H_{A} : Y_{g} < Y_{j},

where g = 1→40, j = 2→40, and g < j. The null hypotheses (H ₀) states that Y_g is larger than Y_j (i.e., Y_g is inferior to Y_j ). The alternative hypothesis (H_A ) states the opposite, that is, Y_g dominates Y_j . A CD analysis is conduced based on the Wilcoxon Signed-Rank test results.

When making decisions about the best parameter combinations, CD analysis (Tsoukias, Perny, and Vincke 2002) considers both positive and negative reasons. The positive reason is that a combination dominates a number of other combinations, which supports the superiority of this combination. Nevertheless, the same combination may also be dominated by other combinations. A concordance measure is defined as C_gj = 1 if Y_g dominates Y_j ; otherwise, C_gj = 0, where g = 1→42, j = 2→42, and g < j. Similarly, a discordance measure is defined as D_gj = −1 if Y_g is dominated by Y_j ; otherwise, D_gj = 0. Finally, the CD index for combination g is defined as $I_{g} = \sum_{j} C_{g j} + \sum_{j} D_{g j}$ . A positive value indicates that the positive measures overcome the negative ones. A negative value indicates the opposite.

CD indices are calculated based on the Wilcoxon Signed-Rank test results, at the 1 percent, 5 percent, and 10 percent significance levels, and for each (L, W) combination. Table 7 summarizes the results. At all significance levels, the parameter combinations with $l \in [0. 25 - 0. 5]$ and $w \in [0. 7 - 0. 8]$ are superior combinations.

Table 7.

Summary of Concordance/Discordance Analyses.

W	L	1%		5%		10%
W	L	w	l	w	l	w	l
0.1	0.1	0.75	0.35	0.80	0.35	0.80	0.45
	0.3	0.75	0.35	0.80	0.45	0.80	0.45
	0.5	0.75	0.35	0.80	0.45	0.80	0.45
	0.7	0.80	0.45	0.80	0.45	0.80	0.45
	0.9	0.80	0.25	0.80	0.45	0.80	0.45
0.3	0.1	0.70	0.35	0.75	0.30	0.75	0.30
				0.80	0.35	0.70	0.30
						0.80	0.35
	0.3	0.75	0.40	0.80	0.35	0.75	0.30
				0.75	0.35
	0.5	0.75	0.30	0.75	0.30	0.75	0.30
	0.7	0.75	0.35	0.80	0.35	0.80	0.35
	0.9	0.75	0.35	0.75	0.35	0.75	0.25
				0.75	0.30
0.5	0.1	0.70	0.30	0.70	0.30	0.70	0.30
		0.75	0.30
	0.3	0.70	0.30	0.70	0.30	0.70	0.30
	0.5	0.70	0.40	0.70	0.40	0.70	0.40
						0.70	0.35
	0.7	0.75	0.40	0.75	0.40	0.75	0.40
						0.75	0.35
	0.9	0.80	0.40	0.80	0.40	0.80	0.40
0.7	0.1	0.80	0.40	0.80	0.40	0.80	0.45
	0.3	0.80	0.40	0.80	0.40	0.80	0.40
						0.80	0.45
	0.5	0.80	0.35	0.80	0.45	0.80	0.40
	0.7	0.75	0.40	0.75	0.40	0.80	0.45
	0.9	0.80	0.35	0.80	0.45	0.80	0.45
0.9	0.1	0.80	0.45	0.80	0.45	0.80	0.45
	0.3	0.80	0.45	0.80	0.50	0.80	0.45
		0.80	0.30
	0.5	0.75	0.35	0.80	0.45	0.80	0.45
		0.75	0.40
	0.7	0.75	0.35	0.80	0.50	0.80	0.50
	0.9	0.80	0.45	0.80	0.45	0.80	0.45
Descriptive statistics
Min.		0.70	0.25	0.70	0.30	0.70	0.25
Mean		0.76	0.36	0.78	0.39	0.77	0.39
Median		0.80	0.40	0.80	0.43	0.80	0.40
Max.		0.80	0.45	0.80	0.50	0.80	0.50

Evaluation Summary

Table 8 recapitulates the “optimal” parameter combinations delineated in the previous sections. They are consistent with each other and provide a secure basis for parameter selection. Selecting l within [0.25–0.50] and w within [0.7–0.8] provides better estimates than other parameter combinations. However, this recommendation is made for Arizona only, based on Arizona data. The extension to other states is discussed in Discussion section.

Table 8.

Summary of Parameter Analyses.

	Top 10 S Values	Best MNSs Across w or l	Pareto Frontier	Concordance/Discordance Analysis Based on Wilcoxon Signed-Rank Test
	Top 10 S Values	Best MNSs Across w or l	Pareto Frontier	1%	5%	10%
w	0.6–0.8	0.75–0.8	0.75–0.8	0.7–0.8	0.7–0.8	0.7–0.8
l	0.2–0.55	0.2–0.3	0.25–0.45	0.25–0.45	0.3–0.5	0.25–0.5

Discussion

In this section, the focus is on several issues raised by the proposed methodology, namely: (1) what is the computational burden of implementing the optimization model? (2) can the results obtained for Arizona be extended to other states? (3) what is the applicability of the technique to the whole CBP database? (4) how could this method be applied to other regional economic databases?

The optimization model was solved using the General Algebraic Modeling System (GAMS) software installed on a personal computer with an installed random access memory of 4 GB and an Intel Core 2 CPU 6600@2.40 GHz processor. The basic model, for any given (l, w) combination, has been implemented for Arizona over eight years, sixteen counties, and twenty-one 2-digit level sector. It includes 10,910 single equations organized into fifteen blocks of equations and 13,130 single variables organized into fifteen blocks of variables. The basic run takes, on average, 5.3 seconds. The work space allocated to solve the model is around 35 Mb. The above suggests that a basic state model is solvable within a very short time. Of course, the computation time was much larger in this study, because the basic model was solved 56,250 (225 × 250) times. Preparing the data inputs to the GAMS model is a lengthier, though by no means difficult process. It is outlined in Appendix C.

When extending the research to other states, a similar suppression analysis should be conducted to find the best parameters. As a preliminary analysis for other states, the location parameter l could be obtained by analyzing unsuppressed information. For instance, in Arizona, the optimal range for the location parameter (0.25–0.5) is consistent with unsuppressed data (38.4 percent of Arizona 1998-006 records across all hierarchies). If these records were suppressed and replaced by the appropriate flag interval, the resulting average location parameter would be 0.31, which is consistent with the positively skewed distributions of firm sizes, as observed in the United States (Axtell 2001) and Portugal (Cabral and Mata 2003). With regard to parameter w, this research confirms that considering employment trends should significantly improve the accuracy of the estimates. But it is unclear whether [0.7, 0.8] is the optimal range for other states. Preliminary testing is recommended.

The application of the methodology to the whole CBP database (all states, all counties, all economic hierarchical levels, and multiple years) would, of course, be the holy grail of this research. There is no conceptual obstacle to this expansion, as the optimization model would be expanded to include blocks of equations and variables for each of the fifty US states. The new features of this expanded model would be constraints on (1) the summation of state sector employments, which should match the employment in the US CBP file and (2) the summation of employment at any economic hierarchical level, which should match the corresponding employment at the higher hierarchical level.

The methodology can be applied to the employment data of other US regional economic databases, such as the EC and the Bureau of Economic Analysis (BEA) data, with necessary revisions of the constraints. For instance, the EC uses the NAICS to classify economic activities and follows the same flag strategy as used in CBP. However, county-level employment data are not provided for several sectors,³ and therefore, the county total constraint (equation 3) would have to be modified. In some of these sectors, only state totals are available. Metropolitan and micropolitan area data are available for other sectors (e.g., utilities). In some sectors (e.g., manufacturing), data are available for places and counties. Therefore, the methodology proposed here would have to be tailored to each major sector of the EC, depending upon the geographical detail. In the BEA database, flag information is not available for a suppressed record, and the flag constraint (equation 2) would have to be removed from the model. BEA employment data are available for counties, metropolitan, and micropolitan areas, BEA economic areas, states, and the whole United States. However, the level of industry disaggregation varies across these geographies, and county-level employment data are available only at the one-digit NAICS level. Constraints could be introduced to enforce consistency between county and metropolitan/micropolitan employment data, in addition to consistency with county and state totals.

Conclusion

A regression-constrained optimization methodology for estimating suppressed records in a longitudinal economic database has been developed and implemented with 1999–2006 CBP data for Arizona. The actual suppressed records were first estimated to create a set of twenty-five complete data sets, using twenty-five combinations of the parameters used in the optimization model. Ten random data set samples were then created for each of these twenty-five complete data sets, mimicking the suppression pattern in the original data. The optimization model was then applied to these suppressed data sets while varying two parameters: the location l of the estimation target within its interval, and the weight w used to trade-off the two objective functions. A goodness-of-fit measure was defined and calculated to measure the quality of the estimates. The results show that accounting for employment temporal trends strongly improves the accuracy of the estimates. Two decision-theoretic approaches (PF and CD) were then implemented to uncover possibly dominant parameter combinations. These analyses consistently point to the ranges of $l \in [0. 25 - 0. 5]$ and $w \in [0. 7 - 0. 8]$ , and these ranges are independent of the way the complete data sets were created.

Several issues related to the applicability and possible extensions of the proposed methodology were discussed. Further research could focus on (1) the adaptation and application of the methodology to other databases (e.g., EC and BEA), accounting for their specific features and (2) assessing the benefits of using data in other databases to help reduce the uncertainty in the estimation of suppressed data in any given database. In addition, spatial econometric models accounting for the possible spatial autocorrelation of employment between any county and its surrounding counties could help obtain more precise estimates. This extension is supported by Desmet and Fafchamps (2005), who analyze the spatial concentration of employment across US counties and provide evidence of spatial spillover effects.

Footnotes

Appendix A

Appendix B

Appendix C

Authors’ Note

An earlier version of this article was presented at the Mid-Continent Regional Science Association 43rd Annual Conference, Bloomington, Minnesota, 2012. The comments and suggestions of five referees on earlier drafts of this article are very much appreciated.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Notes

References

Axtell

R. L.

2001. “Zipf Distribution of U.S. Firm Sizes.” Science 293(5536): 1818–20.

Cabral

L. M. B.

Mata

. 2003. “On the Evolution of the Firm Size Distribution: Facts and Theory.” The American Economic Review 93(4): 1075–90.

Desmet

Fafchamps

. 2005. “Changes in the Spatial Concentration of Employment across US Counties: A Sectoral Analysis 1972-2000.” Journal of Economic Geography 5 (3): 261–84.

Dowall

D. E.

1996. “An Evaluation of California’s Enterprise Zone Programs.” Economic Development Quarterly 10(4): 352–68.

Ellison

Glaeser

E. L.

. 1997. “Geographic Concentration in U.S. Manufacturing Industries: A Dartboard Approach.” The Journal of Political Economy 105(5): 889–927.

Gardocki

B. C.

Baj

. 1985. Methodology for Estimating Nondisclosure in County Business Patterns. DeKalb: Center for Governmental Studies, Northern Illinois University.

Guimaraes

Figueiredo

Woodward

. 2009. “Dartboard Tests for the Location Quotient.” Regional Science and Urban Economics 39(3): 360–64.

Henderson

1997. “Externalities and Industrial Development.” Journal of Urban Economics 42:449–70.

Isserman

A. M.

1977. “The Location Quotient Approach to Estimating Regional Economic Impacts.” Journal of the American Institute of Planners 43(1): 33–41.

10.

Isserman

A. M.

Westervelt

. 2006. “1.5 Million Missing Numbers: Overcoming Employment Suppression in County Business Patterns Data.” International Regional Science Review 29(3): 311–35.

11.

Kennedy

2003. A Guide to Econometrics. Cambridge, MA: MIT Press.

12.

Kreahling

K. S.

Smith

S. M.

Frumento

. 1996. “A Procedure for Estimating Nondisclosed Employment in the County Business Patterns Data Set.” Staff Paper Series No. 299, Department of Agricultural Economics and Rural Sociology, Pennsylvania State University, University Park, PA.

13.

Miller

R. E.

Blair

P. D.

. 1985. Input-output Analysis: Foundations and Extensions. Englewood Cliffs, NJ: Prentice Hall.

14.

Moore

W. S.

2003. “Enterprise Zones, Firm Attraction and Retention: A Study of the California Enterprise Zone Program.” Public Finance and Management 3(3): 376–92.

15.

Pagoulatos

2004. “Interactions between Economic Growth and Environmental Quality in U.S. Counties.” Growth and Changes 35(1): 90–108.

16.

Sechrist

R. P.

1986. “Procedures for Estimating Non-disclosed Data from County Business Patterns.” Research in Contemporary and Applied Geography 10(1): 16.

17.

Stevens

B. H.

Moore

C. L.

. 1980. “A Critical Review of the Literature on Shift-share as a Forecasting Technique.” Journal of Regional Science 20(4): 419–37.

18.

Tsoukias

Perny

Vincke

. 2002. “From Concordance/Discordance to the Modeling of Positive and Negative Reasons in Decision Aiding. Aiding Decisions with Multiple Criteria.” International Series in Operations Research and Management Science 44(3): 147–74.

19.

Weddleton

Olson

D. O.

. 1988. “Estimating Undisclosed Earnings Data.” The Annals of Regional Science 22(1): 84–90.

20.

Zhang

Guldmann

J.-M.

. 2009. “Estimating Suppressed Data in Regional Economic Databases: A Goal-programming Approach.” European Journal of Operational Research 192:521–37.