Automated generator of optimal models for the statistical analysis of data

[0001] The U.S. Government has a paid-up license in this invention as provided for by the terms of contract numbers N00014-99-M-0031 and N00014-99-C-0424 awarded by the Office of Naval Research.

[0002] 1. Technical Field of the Invention

[0003] This invention relates generally to statistical modeling and analysis. This invention relates specifically to the automated process of building, analyzing and creating a hierarchy of statistical models fit to a test data set.

[0004] 2. Description of the Related Art

[0005] There has been a longstanding need for a method which would provide automated statistical analysis employable in a wide variety of applications, including detailed analysis of initial data mining results. The method would produce answers to several key questions, such as

• • [0006] a. How different factors affect the probability of a given event;
  • [0007] b. Which factors are most influential in determining the probability of a given event or the number of occurrences of a given event; and
  • [0008] c. How different factors affect the average number of occurrences of a specific event.

[0009] Examples of these questions are:

[0010] How do race, age, and gender affect the probability that a person will suffer a heart attack before age 40? What is the probability for a particular profile?

[0011] How do age range, gender, education level, and work experience affect the probability that a person's income will fall into a particular salary range? What is the most determining factor for women in particular?

[0012] How many traffic accidents (per 5 year period) should we expect from a specific category of insured drivers? Which factor most heavily affects the expected number of accidents?

[0013] There are countless situations in which businesses and researches need answers to these or similar questions. In general, they have available only limited sample data to estimate the desired probabilities or average number of events. Unless generated using a large number of test cases, these sample statistics tend to reflect anomalies in the data rather than true probability trends.

[0014] It is accordingly a primary object of the present invention to provide what is not available in the art.

[0015] Provided is an automated process for producing accurate statistical models to fit sample data tables. The process constructs and analyzes all reasonable statistical models that fit the sample data table provided, solves for the optimal parameters of the statistical models, and computes test statistics to check the fit of each model with the sample data table. The process uses these test statistics, combined with information about the degrees of freedom in the model, to compute a model score that allows the process to generate a complete ordering of the statistical models analyzed and to present a model that is deemed to be most suitable to the sample data table. In cases where the number of possible statistical models is prohibitively high, the process begins by constructing and solving more general models and then constructs and solves those more detailed models that are similar to those general models which are highest in the generated ordering. The process then presents the model with the highest ordering, as well as a set of average table values corresponding to the statistical model, including average values for table entries that were lacking initial information.

[0016] Another embodiment of the present invention additionally includes an automation of modeling using general linear models, where the only required inputs are the sample data table and, in the case of logistic modeling, an indication of which variable is to be considered as the response variable.

[0017] Another embodiment of the present invention additionally includes specifying which covariates must interact in the model, specifying which covariates are considered to be ordinal and specifying at what stage(s) of the process ordinal models should be considered.

[0018] Another embodiment of the present invention additionally includes automatically constructing and solving for the optimal parameters of all possible hierarchical nominal or ordinal models.

[0019] Another embodiment of the present invention additionally includes automatically determining model test statistics and residual degrees of freedom for each statistical model and using these test statistics and degrees of freedom to establish a complete ordering of statistical models.

[0020] Another embodiment of the present invention additionally includes specifying that a heuristic, rather than exhaustive, exploration of statistical models be performed on the nominal models, or on the nominal and ordinal models, or on the ordinal models only

[0021] This invention finds application in the area of credit scoring, where covariates such as age, profession, gender, and credit history are used to determine the likelihood that an individual will default on a loan. It also finds application in analyzing the effectiveness of many types of tools as they are used in various environments (e.g., the effectiveness of radar when used in different weather conditions). It also finds application in the area of insurance, where one wishes to estimate the future number of claims against a specific insurance policy based on a database of past insurance claims.

[0027] The present invention employs discrete category regression methods similar to linear regression methods. Both categorical and linear regression methods fall into a general category known as General Linear Models. The models dealing with categorical data are more advanced than those dealing with linear data, using more sophisticated modeling and requiring more computationally intensive approaches.

[0028] The statistical and numerical tools required for analysis of General Linear Models are not at all straightforward. Even advanced statistical software, while taking away some of the burden of computation, still requires the user to form and interpret individual models. In addition, the number of models to choose from grows exponentially in the number of factors that are considered.

[0029] The present invention completely automates the process of constructing, analyzing and comparing statistical models fit to a sample data set. Once the user has entered the sample data, the invention does all of the work, forming and solving tens, hundreds, or thousands of possible statistical models and determining the quality of each model. The user can view a ranking of which factors are scored as being most important and can view a variety of graphs showing how different factors affect the overall outcome. The user can also view the results in a variety of display formats.

[0030] Illustration

[0031] To illustrate how our method can answer two example questions, consider the following scenario:

[0032] A pharmaceutical company is studying the distribution of a certain drug. The company has enlisted the help of a sample group of physicians to determine the average number of prescriptions written based on the following list of factors: patient gender, patient age bracket, physician specialization, and geographic region within the country.

[0033] The different factors were broken down as shown in the following Table 1:

[0034]

TABLE 1
Factors affecting number of prescriptions written
			Patient
Geographic Region	Patient Gender	Specialization	Age Bracket
Northeast	Male	Pediatrics	0–11
Mid-Atlantic	Female	General	12–18
South		Internal Medicine	19–35
Southwest		Geriatrics	36–65
Midwest			65+
Northwest
Pacific

[0035] Each physician reported the number of prescriptions written to patients in each category during the test period and the results were combined into a sample table. Our process started with this sample table and recommended a best statistical model, giving information about which factors are most important in determining how many prescriptions will be written, which groups of factors seem to operate together and which seem to operate independently, and what the true average number of prescriptions for each category seems to be.

[0036] In constructing this illustration, we simulated our test data by taking a random sample from a predetermined distribution. Thus, we have the “true average” number of responses for each combination of factors. We use this simulated data set as the measured data. We can compare both the sample and modeled data sets to this “true average,” giving us an idea of the benefits of using the modeled data as opposed to the sample data.

[0037] In this particular example, we must choose from among 166 possible nominal and 66,558 possible ordinal models. Our method can automatically generate and analyze all of these models, or we can ask for a quicker, heuristic search, in which case the method efficiently searches the space of all possible models to arrive at a “best” model. We can then display the groups of factors that interact in the resulting model. See FIG. 1.

[0038] From FIG. 1, we see that we can consider region independently, but that gender and patient age should be viewed together, as should physician specialization and patient age. In other words, we could well make a statement such as “More prescriptions tend to be given in the Southwest than in the Midwest.” Such a statement could be made without specifying a gender, age bracket or physician type. However, we could not say, “More prescriptions tend to be given to males than to females” unless we also include more detailed information about age, such as, “For patients over the age of 65, more prescriptions tend to be given to males than to females.”

[0039] We can get an idea of which of the four factors is most influential in the model output by our method by ranking the factor effects. The result is the menu displayed in FIG. 2.

[0040] We can proceed to view the total number of prescriptions as a function of specialization. See FIG. 3.

[0041] Note that the measured (sample) data and the modeled data agree. The other curve is the actual mean number of responses which were used to form the distribution of the sample data. In the aggregate data, summed over all ages and regions and genders, there is very little difference between the sample, modeled, and true curves. This is to be expected, since anomalies in the random data get smoothed out over large data samples. This changes when we look at a specific response category, and the benefits of using well-modeled data become more evident.

[0042] We proceed to analyze the data in more detail. The graphs in FIG. 4 show the number of prescriptions written for one specific group: males between the ages of 19 and 35, living in the Southwest.

[0043] The jagged measured data curve in FIG. 4 bears the marks of the variance of the distribution, but the modeled curve matches the true mean number of prescriptions almost exactly. It is at this level of detail that we most clearly see the advantages of using the generated “best” statistical model to analyze sample data, rather than simply analyzing the sample data directly.

[0044] Those familiar with Linear Regression will recall how a statistical model can be used to smooth out raw data and determine trends. Linear Regression, however, requires both continuity and ordering in the variables studied. Categorical Data Analysis is used with data which falls into distinct categories. In this analysis, we work with two types of statistical models, each answering a different question about the data studied.

[0045] Logistic Regression Models

[0046] In Logistic Regression Models, the dependent variable follows a multinomial distribution (J possible responses) with probabilities determined by the independent variables. Define π_j(x_i) as the probability of an event in the i'th scenario falling into the j'th response category (each scenario consists of a different combination of categories of the independent variables). The logistic regression model is

[0047] ${\pi }_j\left(x_i\right)=\frac{e^{g_j\left(x_i\right)}}{1+\sum _{k=1}^J{e}^{g_k\left(x_i\right)}},$
where g_j(x_i) is a function which depends on the details of the model (more on this later) and g₀(x_i)=0. Note that we can recover g(x) using the logit transformation, defined as

[0048] $g_j\left(x_i\right)=\ln \left(\frac{{\pi }_j\left(x_i\right)}{{\pi }_0\left(x_i\right)}\right).$

[0049] Once we have determined the interactions in the model, we need to determine a set of model parameters. These are constants that determine how heavily each factor or combination of factors is weighed. These parameters contribute to the calculation of g_j(x_i) and thus to the value of π_j(x_i). There is a set of parameters corresponding to each g_j(x_i) (j>0).

[0050] For fixed j, several scenarios will share some of the model parameters (such as the constant term). The Newton-Raphson method is used to determine which parameters β maximize the likelihood function

[0051] $l\left(\stackrel{->}{\beta }\right)=\prod _{i=1}^n\prod _{j=0}^J{{\pi }_j\left(x_i\right)}^{y_{i,j}}.$

[0052] Here y_i,j is the number of occurrences of the j'th event in the i'th scenario in the raw data.

[0053] Loglinear Models

[0054] Loglinear Models estimate the number of events that occur given a combination of factors. Let n_i denote the number of recorded events in the i'th scenario. The number of events in each scenario is modeled as a Poisson distribution with mean m_i, where log m_i=g(x_i).

[0055] Again, Newton-Raphson is used to maximize the likelihood function. Ignoring constants and taking the logarithm, it is equivalent to maximize the function

[0056] $L\left(\stackrel{->}{m}\right)=\sum _i\left[n_i\log \left(m_i\right)-m_i\right].$

[0057] Model Building

[0058] We now go into more detail about the function g(x) (or g_j(x)). For both logistic regression and loglinear modeling, there is a corresonding design matrix X, which gives a linear transformation on the parameters used in the model. Each row of X contains the multipliers for the parameters to indicate how they are used in each particular scenario. These multipliers are 0 or 1 for nominal models and 0 or a scaling constant for ordinal models. The height of X is the number of scenarios and the width is the number of parameters in the model.

[0059] We now define g(x_i):

[0060] $g\left(x_i\right)=\sum _{j=0}^k{\beta }_j{x}_{\mathrm{ij}},$
where k is the number of parameters in the model. The term g_j(x) is defined in the same way, except that the sum is only over those parameters which correspond to the j'th response. If the width of X were split evenly into J-1 rows, each group of columns would be identical to the others, and g_j(x) would be defined using the product of the j'th group with the corresponding β terms.

[0061] We have a wide variety of models to choose from, depending on which factors and which combination of factors we incorporate into the model.

[0062] The linear function g(x) has the general form:

[0063] $\begin{array}{c}g\left(x\right)=\mathrm{constant}\mathrm{term}+\\ \mathrm{term}\mathrm{depending}\mathrm{on}\mathrm{which}\mathrm{category}\mathrm{of}\mathrm{the}\mathrm{factor}\mathrm{we}\mathrm{are}\\ \mathrm{considering}\left(\mathrm{for}\mathrm{each}\mathrm{factor}\mathrm{included}\mathrm{in}\mathrm{the}\mathrm{model}\right)+\\ \mathrm{term}\mathrm{depending}\mathrm{on}\mathrm{which}\mathrm{combination}\mathrm{of}\mathrm{the}2\mathrm{factors}\mathrm{we}\mathrm{are}\\ \mathrm{considering}\left(\mathrm{for}\mathrm{each}2\text{-}\mathrm{factor}\mathrm{interaction}\mathrm{term}\mathrm{in}\mathrm{the}\mathrm{model}\right)+\\ \mathrm{term}\mathrm{depending}\mathrm{on}\mathrm{which}\mathrm{combination}\mathrm{of}\mathrm{the}3\mathrm{factors}\mathrm{we}\mathrm{are}\\ \mathrm{considering}\left(\mathrm{for}\mathrm{each}3\text{-}\mathrm{factor}\mathrm{interaction}\mathrm{term}\mathrm{in}\mathrm{the}\mathrm{model}\right)\dots \\ \mathrm{etc}.\end{array}$

[0064] The parameters β correspond to the terms added above. Suppose, for example, that factor X has three categories, factor Y has four categories, and that XY interaction is included in the model. For a particular scenario, we would compute g(x) as follows:

[0065] $\begin{array}{c}g\left(x\right)=\mathrm{constant}+\\ \mathrm{term}\mathrm{corresponding}\mathrm{to}X\mathrm{category}\left(3\mathrm{possibilities}\right)+\\ \mathrm{term}\mathrm{corresponding}\mathrm{to}Y\mathrm{category}\left(4\mathrm{possibilities}\right)+\\ \mathrm{term}\mathrm{corresponding}\mathrm{to}\mathrm{XY}\mathrm{combination}\left(12\mathrm{possibilities}\right)+\\ \mathrm{any}\mathrm{other}\mathrm{terms}\mathrm{in}\mathrm{the}\mathrm{model}\end{array}$

[0066] We only consider hierarchical models, requiring that any model that includes an n-factor interaction term also include the (n−1)-factor terms from which it is built. So, for example, a model with XYZ interaction must also have XY, XZ and YZ interactions. Parameters satisfy zero-sum constraints at each interaction. In the example above, once we have determined two of the possible X category parameters, then the third must be the negative of the sum of the other two. Likewise the first eleven of the XY interaction parameters determine the last parameter.

[0067] A variable may be considered to be ordinal if its categories have a natural ordering. It is considered to be nominal if the model does not regard its categories as ordered. Taking ordering into account can help smooth the model without the addition of unnecessary parameters. We may use models where variables are considered to be ordinal or strictly nominal. We can also use models where variables behave as ordinals in some interactions and as nominals in others.

[0068] Model Assessment

[0069] Each statistical model has a test statistic which is asymptotically Chi-squared if the model truly reflects the distribution of the data. The degree of freedom of the Chi-squared test statistic is equal to the residual degrees of freedom (rdf) of the model. The rdf of the model is defined as the difference between the number of scenarios and the number of free parameters in the model (not counting those determined by the zero-sum constraints). In other words, the rdf is the ‘free play’ left in the model. If we had one parameter per scenario, we would exactly fit the model to the empirical data. Checking the statistic value against the corresponding Chi-squared table gives us an idea of whether the model is a reasonable fit for the data.

[0070] Logistic Regression

[0071] The test statistic we use for Logistic Regression is defined as X²=Σ_ie_i², where

[0072] $e_i=\frac{y_i-n_i{\hat{\pi }}_i}{{\left(n_i{\hat{\pi }}_i\left(1-{\hat{\pi }}_i\right)\right)}^{\frac{1}{2}}}.$

[0073] If there are more than two possible outcomes in the Logistic Regression model, we must modify the test statistic to assure that the e_i terms are still normally distributed. Let n_i,j denote the number of outcomes j in the i'th scenario. If we have only one random event, then Var(n_i,j)=p_i,j(1−p_i,j) and Cov(n_i,j, n_j,k)=−P_i,jp_j,k. If n_i is the vector (n_i,1, . . . n_i,J) (note that the first response is not included), p_i is the vector (p_i,1, . . . p_i,J) and N_i is the total number of responses for scenario i, then (n_i−N_ip_i)/sqrt(N_i) converges to a multidimensional normal distribution, N(0,Σ), where Σ=diag(p_i)−(p_ip_j)_ij is the covariance matrix of the multinomial distribution as described above.

[0074] Note that if Y˜N(0,Σ) then MY˜N(0,MΣM′). We diagonalize the matrix Σ=VDV′, and then we get V′Y˜N(0,D). The term corresponding to each scenario is then (V′Y)_i/sqrt(D_ii)˜N(0,1). We add the square of each of these terms to get a Chi Squared distribution. A formal proof of the actual number of degrees of freedom has yet to be completed, but we estimate this number as (Number of scenarios)*(Number of response categories−1)−(Number of model parameters). Recall that the probability assigned to the first response category of any scenario is determined by the remaining probabilities in that scenario, and thus we disregard it for all computational purposes.

[0075] Loglinear Modeling

[0076] There are a variety of statistics available for Loglinear Modeling, but we use a general statistic, developed by Cressie and Read.

[0077] $\begin{array}{cc}{X}^2=\frac{2}{\lambda \left(\lambda +1\right)}\sum n_i\left({\left(\frac{n_i}{m_i}\right)}^{\lambda }-1\right)& \left(-\infty <\lambda <\infty \right).\end{array}$

[0078] Recall that the n_i are the empirical counts and the m_j are the modeled counts. This statistic generalizes other statistics in that it is equal to the Pearson statistic when λ=1 and approaches the likelihood-ratio statistic as λ approaches zero. We have chosen to use λ=⅔, a value which Cressie and Read recommend to reduce the effects of sparseness in the empirical data.

[0079] Model Selection and Heuristics

[0080] Although methods exist for making comparisons between different models, it is not clear how to select the ‘best’ of all possible models. In fact, it is only the surge in computational power that makes it feasible to compare as many models as we do. Our method chooses the model which maximizes the score p*rdf, where p refers to the right probability tail of the Chi-squared statistic for the model. In this way, we prefer models that seem to fit the data well without over-fitting.

[0081] Because the number of possible models grows exponentially in the number of factors, we use several heuristics to arrive at near-optimal solutions.

[0082] The heuristic for selecting a near-optimal model with strictly nominal interactions is as follows. If there are n factors, we start with the model with all (n−1)-way interactions. We record the score and compare it to the score with all (n−2)-way interactions. This shrinking process continues until the models stop improving after the m-way interaction model. We return to the model with all m-way interactions and add the (m+1)-way variable which most improves the model's score. We keep adding (m+1)-way variables so long as the model keeps improving. We then return to the model with all m-way interactions and delete the m-way variable whose deletion results in the greatest improvement in score. We continue to delete m-way variables so long as the models improve. We conclude by choosing the best of all models analyzed in this process.

[0083] Once we have chosen a nominal model, we may begin to add ordinal interactions. The heuristic used is to give ordinal interaction to one factor in one variable at the highest level of interaction used in the model. The factor and variable are chosen so as to maximize improvement. We repeat this process so long as the models improve. When improvement stops, we perform the same algorithm one level down, but only allow those ordinal interactions which are preserved at the higher level variables which are built from the variable in question.