As you probably know, there is abundant need for big data analysis and machine learning modeling these days.  One concept that is important in learning big data analysis is the bias-variance trade-off.  In this blog, I’ll describe and illustrate the bias-variance trade-off and explain how understanding this concept can simplify your learning experience by giving you a framework to view hyperparameter tuning in your machine learning models.   Most users of machine learning models are focused on predictive modeling.  Predictive models, also called supervised learning models, involve a target variable (the focus of your business or research) and inputs (the predictors that hopefully inform you about the target).  Supervised machine learning models identify input-target associations in order to predict future values of the target.   The focus of predictive modeling is to make the most accurate predictions possible.  This is a different goal than statistical inference, which usually involves explaining the input-target relationships.  For predictive modelers, explaining the relationships among variables is of secondary importance at best.  To assess the accuracy of the predictions, two common measures are the mean square error (MSE) and average square error (ASE). Typically, we use MSE for measuring error in the training data and ASE when applying a model to a data set such as validation data that was not involved in creating the model. Whichever metric we use, we can view this as the prediction error that we want to minimize.   So, all else being equal, a predictive model with lowest error can be considered a better model.  How do we reduce the prediction error in our machine learning models?  By optimizing the bias-variance trade-off.  Both bias and variance contribute to the prediction error.   What is bias?  In the context of predictions, a model makes biased predictions when the predicted values systematically underestimate or overestimate the actual values.  It is the difference between the expected value of an estimate and its true value.  Bias cannot be measured in real data because it requires knowing the true values which we are trying to learn in the first place.  Models that are too simple to capture the true patterns in the data have high bias.  Typically, we can reduce bias by increasing the complexity of the model (more parameters, higher order terms, etc.).   What is variance? In the context of predictions, a model has high variance when the predicted values are highly variable when the model is applied to different samples from the same population.  You can also think of variance as the sensitivity of the predictions to the particular data set used to train the model. Variance cannot be measured from fitting the model to a single data set as is typically done, although it could be measured by fitting it to several bootstrap samples of the same data.  Note that high variance here does not imply anything about making correct predictions.  It is about the variability in the predictions, not how accurate they are.   Models that are overly complex for the data have high variance.  “Overly complex” here means that the model describes aspects of the observed data that are not representative of the population and are unique to the observed data.  That is, overly complex models model the “noise” in addition to the “signal”.  A model with higher complexity is more flexible and can fit the data with lower bias. But because of the increased flexibility, the model’s predictions can differ dramatically between data sets.   So high bias is associated with models that are too simple to describe the true patterns in the population and high variance is associated with models that are unnecessarily complex. You can imagine that as you increase model complexity from too simple to too complex, the model bias decreases and variance increases.  There is a trade-off, where lowering one increases the other.   Prediction error increases with both the bias and variance:   MSE = bias 2 + variance + irreducible error   , the irreducible error being a component of the MSE unrelated to model complexity that won’t be considered here.  Finding the model that has the lowest error involves finding the right amount of complexity, that is, the bias and variance that minimizes the error. This amount is not something you will know before some trial and error.  The trial and error should be guided by measuring model performance on a holdout validation data set.  It turns out that adding a small amount of bias often results in a large reduction in variance and an overall decrease in prediction error.   Here are illustrations of models with high bias or high variance.  We have a population with the relationship Y= 7+ 11X +3X 2 .  I’ll draw 3 samples from this population and fit models that are too simple (high bias) and models that are too complex model (high variance). Here’s how I generated the data for the first sample:   data sample; call streaminit (123); sample="sample 1"; do i=1 to 50; X=i; y=7+0.3*X+40*X*X +5000*RAND("normal"); true_y=7+0.3*X+40*X*X; output; end; run;   The CALL STREAMINIT routine specifies a seed that starts the random number generator used by the RAND function.  RAND(“normal”) generates a random value from a normal distribution with a mean of 0 and SD=1. Here’s a plot of the true relationship:   Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   When we approximate this true curvilinear relationship with a linear regression model, we are using a model that is too simple and the predictions are biased.  At X-values near zero and near fifty, the predicted values (shown by the red line) are too low—they underestimate the true values. At X-values near 25, the linear regression overestimates Y.     What happens when we use a model that is too complex? I’ll approximate this curvilinear relationship with a spline function.  Splines are flexible, piecewise polynomial functions joined together smoothly at join points called knots.  I’ll use a spline model based on 10 th degree polynomials.  The predictions are no longer biased too low at the extremes and too high near the center.  But this spline model is overly complex, and it captures the noise in the data.  Each little “bump” in the spline is a pattern found in this specific sample and not likely to be found in the next sample drawn from the same population.  The spline model predictions have low bias, but they have high variance.     Let’s see a few more high bias and high variance models fit to other samples from the same population. To generate additional samples, I changed the seed of the CALL STREAMINIT statement in the SAS code above and concatenated the 3 samples into a data set called “combined”.    I used the following PROC SGPANEL code to visualize models fit to these samples:   title "high bias"; proc sgpanel data=bias_var.combined; panelby sample / novarname columns=3; scatter x=X y=y /; reg x=X y=y / nomarkers lineattrs=(color=red); run; title "high variance"; proc sgpanel data=bias_var.combined; panelby sample / novarname columns=3; scatter x=X y=y /; pbspline x=X y=y / nomarkers degree=10 lineattrs=(color=green); run;     We can see the linear regression is biased.  In all three samples, the average predictions are different than the true values.  So, this error is not a fluke.  It is systematic due to using a model that is too simple for these data.  Using a model that is too simple for the data is referred to as underfitting.     The 3 spline models demonstrate high variance.  At first, the curves might seem similar but look at each peak and valley across samples.  They don’t match up.  That is, the predictions vary considerably across different samples.  Models like these can be described as overfitting the data. These models are so specific to the data they were trained on that they will perform poorly at predicting new data. On average, the bias is low, but in practice, we will only have one sample for modeling. We may end up getting the sample and the model with predictions that are far from the average—there is no way of telling. Variance here is the sensitivity of the predicted values to the particular data set used to train the model.   So why does this matter? From a practical point of view, can’t we just focus on picking a model with lowest prediction error on a validation data set and ignore the bias-variance trade-off concept? Well, yes.  But understanding bias and variance gives a conceptual framework that will help with using machine learning models.  One situation this comes up in is tuning hyperparameters.  Hyperparameters are the settings for a machine learning model that are established before the model is fit to the data.  They set the guidelines for how the machine learning model goes about producing the predictions.   Here are some examples of hyperparameters that affect the bias and variance of the resultant model:     Many of the options in supervised machine learning SAS Viya procedures such as PROC BN, PROC FOREST, PROC GRADBOOST, PROC NNET, and PROC SVMACHINE are used to adjust the complexity of the model to optimize the amount of bias and variance.  The same hyperparameters can appear in several models (e.g. L1, also called LASSO).  Adjusting these will help find a model that avoids overfitting or underfitting. For example, a higher L1/LASSO parameter will make a model less sensitive to the training data to reduce the amount of noise being modeled and lower the chance of overfitting. Understanding bias and variance allows you to view a dozen different hyperparameters as methods to make models err on the side of too simple or too complex.   To learn more about the important concepts in machine learning modeling, try the SAS class “Statistics You Need to Know for Machine Learning”.  This course covers the bias-variance trade off as well as many other important concepts that will help you on your journey as a machine learning modeler.  It also is the best preparation for students interested in earning the SAS Credential “SAS Certified Associate: Applied Statistics for Machine Learning”.   See you at the next SAS class!     Links:   Course: Statistics You Need to Know for Machine Learning (sas.com)   SAS Certified Associate: Applied Statistics for Machine Learning     Find more articles from SAS Global Enablement and Learning here. ... View more

Principal component analysis (PCA) is an analytical technique for summarizing the information in many quantitative variables. Often it is used for dimension reduction, by discarding all but the first few principal component axes.  These first few usually explain most of the variability in the original variables.  How can we determine the number of principle component axes to keep for further analysis?  In a previous post (How many principal components should I keep? Part 1: common approaches), I described some of the common approaches to determining the number of PC axes to retain.   These approaches are all subjective and do not distinguish components that summarize meaningful correlations among variables from components that capture only sampling error.  In this post, I will demonstrate how to construct a significance test for the number of PC axes to retain for analysis that avoids these limitations.     How many components to retain from PCA?   In a previous post (How many principal components should I keep? Part 1: common approaches), I described some of the common approaches to determining the number of PC axes to retain. These included keeping components that account for a set percentage of variability, retaining components with eigenvalues greater than 1, and components up to the "elbow" of a Scree plot. These approaches are easy to implement and may be adequate for many research goals.  But each of these approaches is somewhat subjective, and more importantly, they are unable to distinguish principal components that represent correlations among variables from sampling error.  For example, if we create data that is made of pure noise with no true correlation among variables, PCA will still produce eigenvalues and eigenvectors.   To illustrate this, I used PCA on randomly generated data with no structure.  The code below creates a data set randnumbers (n=100) with four variables X1 through X4, each taking on integer values between 1 and 100.  These variables are randomly generated with no underlying structure or associations among them.  Any correlation among them is purely due to sampling error (random chance).  I used PROC PRINCOMP on these data to carry out PCA on the randnumbers data.   data randnumbers; call streaminit(99); do i=1 to 100; x1=rand("integer", 1, 100); x2=rand("integer", 1, 100); x3=rand("integer", 1, 100); x4=rand("integer", 1, 100); output; end; run; proc princomp data=randnumbers; var x1 x2 x3 x4; run;   The PCA on the random data produced the following eigenvalues and Scree plot:   Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.     Saving PCs that explain at least 75% of the variation would result in retaining PC1 through PC3, which collectively explain almost 80% of the variability in the original predictors.  Retaining components with eigenvalues >1 results in treating only PC1 as meaningful. Using the Scree plot bend method suggests retaining the first 2 PCs.  But we know that none of these components represents true structure in the underlying data.     Randomization tests to determine the significant number of principal components   An alternative approach is to use a randomization test to determine the number of principal components to keep.  The idea here is that we want to distinguish meaningful PCs from noise.  That is, we want to know if the PCs explain more variability than similar data with no real structure that summarize no real correlations. To do this, the values of each variable can be randomly shuffled, and the PCA can be recalculated for the new structureless data set.  All the original data is retained in these permuted data sets, but the values of each of the variables are randomized relative to one another.  This randomization followed by PCA is repeated many times (say 100 to 1000 times) and the eigenvalues from each randomly permuted data set are retained.   This produces a null distribution of eigenvalues that can be used to test the hypothesis that the original data has no structure, that is, that the correlations represented by components is entirely due to sampling error.  To carry out the test using alpha=0.95, the first eigenvalue of the real data PCA would be compared to the 95 th percentile of the null distribution of the first eigenvalue.  If it is smaller than the 95 th percentile of the distribution of null distribution of PC1 eigenvalues, then the first PC represents sampling error and not true structure in the original data and none of the PC axes should be retained.  If the first PC eigenvalue is greater than 95% of the null distribution, it represents significantly more correlation among variables than you would expect from random chance.  If the first PC is significant, then proceed to testing PC2 the same way (and so on).  The test could be constructed with other alphas, such as alpha=0.01 by comparing to the 99 th percentile of the null distribution.   Here is how I conducted a randomization test with the PCA on Fisher’s iris data. These data include 4 variables measuring floral morphology for 3 species of iris.  The variables are sepal length, sepal width, petal length, and petal width.  Fifty flowers of each of 3 Iris species were measured for a total of n=150. These data are included with SAS in sashelp.iris.   First, I added a variable called sort to the data which has the consecutive values 1 through 150, then used PROC PRINCOMP to carry out PCA. The results were shown in my previous post (How many principal components should I keep? Part 1: common approaches).   In the macro program irisRand, I used PROC PLAN to randomly shuffle the integers values of sort, while keeping the original order of the rest of the variables.  The output data were sorted by sort and only sort and one floral variable were retained.  This was repeated for each of the 4 floral variables and the randomly shuffled values were merged by the macro program null_dist_PCA.  PROC PRINCOMP was run on these randomized data and the eigenvalues were saved, transposed, and appended to a data set null_EVs.  This was repeated 1000 times and PROC UNIVARIATE was used to find the 99 th and 95 th percentiles of the null distribution of each of the 4 principal components.   Here is the SAS code:   /* create variable sort for randomization */ data iris; set sashelp.iris; sort=_n_; run; /*find eigenvalues for original iris data */ proc princomp data=iris; var SepalLength SepalWidth PetalLength PetalWidth; run; /*randomize the order of sort variable */ %macro irisRand (trait=, traitno=); proc plan; factors sort=150/noprint; output data=iris out=randomized; run; /* sort by sort number to randomize the rows. Drop everything but one trait*/ proc sort data=randomized out=random&traitno (keep=&trait sort); by sort; run; %mend irisRand; /* randomize values of 4 traits and merge the data sets, carry out PCA and save eigenvalues, transpose and append eigenvalues, repeat &iter=1000 times */ %macro null_dist_PCA (iter=); %Do j=1 %to &iter; %irisRand (trait=sepallength, traitno=1); %irisRand (trait=sepalwidth, traitno=2); %irisRand (trait=petallength, traitno=3); %irisRand (trait=petalwidth, traitno=4); data randomized_iris; merge random1-random4; drop sort; run; ods select none; proc princomp data=randomized_iris; var SepalLength SepalWidth PetalLength PetalWidth; ods output Eigenvalues=ev(keep=eigenvalue); run; ods select all; proc transpose data=ev out=ev_transposed (drop=_name_) prefix=eigen; run; proc append data=ev_transposed base=null_EVs; run; %end; %mend null_dist_PCA; /*prevents overflow of log*/ proc printto log='randomization_log.txt' new; run; %null_dist_PCA (iter=1000); /*resets log to original destination*/ proc printto log=log; run; proc univariate data=null_EVs noprint; var eigen1 - eigen4; output out=crit pctlgroup=byvar pctlpts=99 95 pctlpre=eigenvalue1_ eigenvalue2_ eigenvalue3_ eigenvalue4_ pctlname=P99 P95; run; proc print data=crit; run;   Here are the eigenvalues from the original iris data:      And here are the 99th and 95th percentile of the null distribution for each of the four eigenvalues:     So, using the cutoffs for statistical significance of the 99 th percentile (or 95 th ), only the first eigenvalue is significant.  That is, we cannot reject the null hypothesis that the correlations summarized in PC2 through PC4 are due to random noise.  This is a different result than what would be achieved from retaining PCs that explained at least 75% of the variability in the original predictors.  And importantly, it is less subjective than several other approaches (despite the choice of alpha itself being subjective).  Of course, the decision on which method to use will depend on the specifics of your research question.  Regardless of your favorite method, it can be useful to have this tool in your statistical toolkit.   For more information on PCA and other multivariate techniques, try the SAS course Multivariate Statistics for Understanding Complex Data.  You can access this course as part of a SAS learning subscription (linked below).   See you at the next SAS class!     Links:   Course: Multivariate Statistics for Understanding Complex Data (sas.com)   SAS Learning Subscription | SAS     Find more articles from SAS Global Enablement and Learning here. ... View more

Principal component analysis (PCA) is an analytical technique for summarizing the information in many quantitative variables. Often it is used for dimension reduction, by discarding all but the first few principal component axes.  The first few components usually explain most of the variability in the original variables.  How can we determine the number of principal component axes to keep?  There are several commonly used approaches which I will describe in this post.  The number of axes used often depends on a percentage of variability explained, component axes with eigenvalues greater than 1, or visual inspection of a Scree plot.   In a follow up post, I will describe how to construct a significance test for the number of principal component axes to retain for further analysis.     What is PCA?   PCA is a widely used summarization and dimension reduction technique used with multivariate data. It starts with a set of original variables to be summarized.  New derived variables are constructed that are weighted linear combinations of these original variables. These derived variables are called principal components [PCs], and the number of PCs created are the same as the number of original variables.   For example, let’s look at a scatter plot of the heights and weights of 19 students from the sashelp.class data set:   Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   Height and weight are positively correlated.  A principal component analysis of these data would result in 2 principal component axes, the first going through the widest spread of the scatterplot and the second axis perpendicular to the first:     These new variables, the PC axes, are functions of the original variables, in this case PC1= 0.71*height+0.71*weight and PC2=0.71*height−0.71*weight.  When the PCs are functions of only a few variables as in this case, they can be interpreted.  The first PC which is based on positive weights or loadings of the height and weight variables could be interpreted as a measure of overall size or stature. In one sense, PC1 could be considered the “real” variable, and height and weight could be considered two ways of measuring overall size.  The second PC could potentially describe how slim vs stocky the students are for their given size, but we can’t directly see that from the data. When the PC are functions of more than a few variables they generally become uninterpretable.   With a third variable such as age included in the data, a third PC would be produced. This third PC axis would be perpendicular to the first two.  When PCs are used for further analysis, usually only the first few PCs are retained, as they explain most of the variability in the original predictors. Running PCA on height, weight, and age shows that the first two components explain over 96% of the variation in the three original variables.  This is the basis for using PCA for dimension reduction—there is minimal loss of information when PC3 is not used in future analyses.     Dimension reduction through PCA can be an early step in developing machine learning models.  For data scientists with many predictors, PCA can be used for replacing many of the original inputs with relatively few PCs.  Additionally, using PCs in place of the original inputs has the added benefit of removing any potential collinearity problems.  Collinearity describes strong correlations among sets of predictors, and it can lead to unstable parameter estimates with large standard errors. This is not a problem with PCs because they are constructed in a way that ensures all PCs are uncorrelated.     How principal components are constructed   The first principal component, PC1, is constructed such that it has the greatest variance of any possible linear transformation of the original variables.  Further PCs are constructed such that PC2 accounts for the second greatest proportion of the variability among the original variables, PC3 accounts for the third greatest proportion and so on.  The variances of the PCs are called eigenvalues.  The PCs have the property that the correlation among any principal components is exactly zero.   PCs can be constructed based on correlations among variables or based on covariances.  Using correlations is a much more common approach.  PCA based on covariances is highly sensitive to the scaling of variables.  The original variable with the greatest variance will typically be associated with the first PC axis to the exclusion of all other variables.  To remove the effect of scale, correlations are preferred unless the original variables are already on a similar scale.  Principal component analysis can be carried out programmatically using the SAS 9 procedure PROC PRINCOMP or the SAS Viya procedure PROC PCA.   Because the first few PCs tend to explain the majority of variation, often the remaining PCs can be discarded with minimal loss of information.  The assumption is that early components summarize meaningful information, while variation due to noise, experimental error, or measurement inaccuracies is contained in later components.  When PCs are discarded in this manner, the number of dimensions (variables) can be greatly reduced. And unlike the original variables, it is impossible to have collinearity problems with the PCs since they are perfectly uncorrelated.  But how many of the PCs should be kept for further analysis? If too few components are retained, future analyses using the components will suffer from loss of the relevant information.  If too many components are retained, the additional noise can distort the underlying pattern of correlations among variables that the components are supposed to summarize.     How many PCs to retain?   One method to decide how many PC to retain for analysis is to pick a proportion of variance to explain that sounds reasonable.  Would you be content to keep PCs until 75% of the variability in the original variables is accounted for?  How about at least 90%? You may not even need to decide advance, and simply look at your PCA results and pick the number of PC axes that save most of the original variability. While being arbitrary, it may be reasonable to do this for some research goals.   Another approach is to retain all PCs that have eigenvalues greater than 1.  When PCA is based on correlations, the sum of the variances of PCs (the eigenvalues) will equal the number of predictors being summarized.  This means the average of the eigenvalues equals 1. So, using eigenvalue>1 as a threshold for considering a PC to be meaningful amounts to keeping all PCs that are account for more variability than the average.   A third approach involves examining a plot of the variance of each PC vs the PC number, called a Scree plot. “Scree” means an accumulation of rocks at the base of a mountain or cliff.  If you picture the profile of a cliff, bending to meet the ground, this is the typical shape of a Scree plot and where it gets its name (thanks, Wikipedia).   Typically, the curve connecting the components will have a bend, making an “elbow” shape.  This point at which the curve flattens out indicates the maximum number of PCs to retain, and those below the elbow are discarded.  What’s the justification for this? Recall that the eigenvalues decrease sequentially from the first component to the last.   So, the Scree plot typically shows a steep decline that asymptotically approaches zero. Where the graph straightens out, the components are explaining a relatively small proportion of variance.   To illustrate these approaches, I used the SAS 9 procedure PROC PRINCOMP to analyze Fisher’s famous iris data set.  These data contain 4 variables measuring floral morphology for 3 species of iris.  The variables are sepal length, sepal width, petal length, and petal width.  Fifty flowers of each of 3 iris species were measured for a total of n=150. These data are included with SAS in sashelp.iris.   proc princomp data=sashelp.iris plots=all; var _numeric_; run;   The eigenvalues table shows the variances (eigenvalues) of the four principal components, along with their differences, their proportion of the total variance, and the cumulative variance explained:     If we want to save components that represent at least 75% of the variability, the first two PCs would be retained. These account for over 95% of the variability in the original 4 predictors.  If instead, eigenvalues >1 were retained, only the PC1 with an eigenvalue of 2.9 would be kept.  All remining PCs have eigenvalues below 1.   The Scree plot shows a steep decline in eigenvalue from PC1 to PC2 and a flattening out of the curve between PC3 and PC4:     Depending on where one sees the elbow, examining this Scree plot suggests keeping 3 or possibly only 2 PC for further analysis. Note that some authors suggest discarding the elbow point as well and only retaining the PCs before the bend.   There are several other graphs that can be produced to help interpret the principal components.  To learn more generally about the graphical output of PROC PRINCOMP, see Rick Wicklin’s blog How to interpret graphs in a principal component analysis - The DO Loop (sas.com).   An advantage of these approaches is that they are all easy to implement.  If the goal is dimension reduction, any of these approaches may be reasonable.  But when focus of the PCA is summary of the data with meaningful axes of variation, the choices presented here are might not be satisfactory.  They rely on subjective, somewhat arbitrary choices.  An additional downside of these methods is that none of them consider that the correlation captured by the principal components could be entirely due to sampling error.   In my next post, I’ll demonstrate an approach for constructing a test for significant principal components that can overcome these limitations.   For more information on PCA and other multivariate techniques, try the SAS course Multivariate Statistics for Understanding Complex Data.  You can access this course as part of a SAS learning subscription (linked below).   See you at the next SAS class!     Links:   How to interpret graphs in a principal component analysis - The DO Loop (sas.com).   Course: Multivariate Statistics for Understanding Complex Data (sas.com)   SAS Learning Subscription | SAS     Find more articles from SAS Global Enablement and Learning here. ... View more

Last year, a student in the SAS Visual Statistics: Interactive Modeling Building class I was teaching asked about modeling a continuous, skewed response variable with a large number of zeros.  He asked for some possible approaches using SAS Visual Statistics, and we ended up talking about two-stage modeling.  In this post, I’ll walk through how we used SAS Visual Statistics to fit a two-stage model to the course data set (VS_Bank).   VS_bank data   VS_Bank is a simulated data set from a large financial services firm’s accounts.  It has 1.1 M rows, with 12 interval-scaled predictors (Logi_RFM1–Logi_RFM12) and 2 categorical predictors (Cat_input1, Cat_input2).  The data set has 3 target variables: a binary target indicating if a purchase was made or not (B_tgt), an interval target indicating the total value of the purchases (Int_tgt) and a count target indicating the number of purchases (Cnt_tgt, not used in this post).  The data set has roughly 20% purchasers (B_tgt=1) and 80% non-purchasers (B_tgt=0).  The large number of zeros for the total value of purchases were coded as missing for the interval target (Int_tgt had 80% missing values).  All the missing predictor values were imputed (Int_tgt was not imputed) and interval predictors were log transformed to reduce skewness prior to modeling (hence the “logi_” prefix).   Two-stage modeling   The goal we chose for modeling with the VS_bank data was to predict the most valuable potential customers for the bank’s marketing department to target. The two-stage modeling approach involves modeling separately whether a customer is a likely purchaser (the stage-1 model using the binary target, B_tgt), and to use the predicted probability of making a purchase as an input for modeling the total value of purchases (the stage-2 model using the interval target, Int_tgt).  Incorporating the predictions from the stage-1 model into the stage-2 model can likely improve its predictive power, particularly when the $0 total value customers (i.e., the non-purchasers) were coded as missing values for Int_tgt.   Note that the same or different predictors can be used for the two models.   The predicted values from both models can be multiplied.  Then customers can be ranked from highest to lowest value for marketing.  Why multiply the probability of purchasing by the predicted total value of purchases? The product can be thought of as the expected future profit, and multiplying can give an improved assessment of customer profitability, compared to looking at the predictions for total value of purchases (the stage-2 model target) alone.  This is because the customers who are predicted to make the largest dollar amount of purchases might not be very likely to make any purchase at all.  If the predicted probability of B_tgt=1 and Int_tgt are negatively correlated, using only the predicted value of Int_tgt would paint a misleading picture of the customers most valuable to target for marketing.   This is essentially what we found after producing a stage-1 model and plotting the relationship between the stage-1 model predictions and the stage-2 target:     Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   Stratified partitioning for honest assessment   We decided to partition the data into training and validation sets and asses each potential stage-1 and stage-2 model on its validation performance.  The simplest approach to partitioning is to use a random sample of say 50% for training and 50% for validation.  Instead, I recommend using a stratified random sample for each partition.  For the stage-1 model (predicting a binary target), the partitions should be stratified based on the target.   Stratification means that when making the training and validation data sets, the exact same frequency of purchasers will be in each partition.  This makes the performance on the partitions more comparable than if there were say 22% purchasers in the training data and 18% purchasers in the validation data.  A stratified partition can be created using the New Data Item option in Visual Analytics, then choosing Partition:     For the interval target models, care should be taken with creating the partitions.  Monetary values are often highly skewed, and this is the case with the interval target, Int_tgt.  If, by chance, all the rare high value purchasers ended up in the validation data with none in training, the model performance will suffer.  To avoid this, a partition can be created that is stratified on a binned version of the interval target.  We did this using the New Data Item option to create a measure variable “binned_purchases”: missing values were put in bin 0, total purchases <= $50K went into bin 1, total purchases <= $100K went into bin 2, and purchases greater than $100K were put into bin 3 (the max was $500K):   We then used binned_purchases to create a categorical variable (bins) by choosing New Data Item, then Custom Category.  This categorical variable was used for stratification of the second partition we created.  SAS Visual statistics will hold one partition variable at a time, so the stage-1 model partition was hidden when the stage-2 partition was created.   Stage-1 models   We tried three stage-1 models for predicting the binary target B_tgt: logistic regression with backwards elimination based on SBC, logistic regression with stepwise selection based on AIC and a decision tree using default settings. The models were compared on their misclassification rate when applied to the validation data:   Model Validation misclassification Logistic backwards (SBC) 0.1615 Logistic stepwise (AIC) 0.1614 Decision Tree 0.1645   The AIC-based stepwise logistic model had the best performance and the Derive predicted option was used to add the predicted probabilities to the Data pane.   Stage-2 models   All of the stage-2 models used the predicted probability from the final stage-1 model as an input. How should we model the interval target? A linear regression model would be inappropriate for these data because of the high skewness of the target, Int_tgt.  We decided to try two generalized linear models: lognormal and gamma regressions.  These models both use continuous distributions that can account for skewness in the response.  Below is a picture of the skewed lognormal and gamma distributions, compared to a (truncated) normal distribution.  Note that the lognormal distribution has heavier tails than gamma and are more skewed than a truncated normal distribution: Image modified from SAS Statistics 2: ANOVA and regression course notes   Another advantage to using lognormal and gamma regression models is that these distributions are strictly positive.  So, unlike a linear regression model, these models will never produce a nonsensical prediction of a negative valued purchase.  The lognormal model and the gamma model are available using a Generalized Linear Model object and choosing the normal distribution with a log link and the gamma distribution with the log link, respectively from the options pane:     How do we get predictions for the 80% of the customers that did not make a purchase?  The $0 purchasers were coded as missing values for the target.  When a model is built in SAS on data that include missing values of the target, the model is fit on the only the rows with complete data (200K rows here), but predictions are made on the complete data set (1.1 M rows). Sometimes a researcher will use this to get predictions from a model by concatenating a data set to be scored (missing the target) with the data used for fitting a model, then calculating the predictions.  This is sometimes called the missing-Y trick.   It’s likely the predicted probability from the stage-1 model is non-linearly associated with the interval target.  How should we account for this non-linear relationship? One way to model a non-linear relationship is to use a spline function. A spline is a piecewise polynomial function where the pieces are smoothly joined at knots. Splines are very flexible and can be used to model complex non-linear relationships that are difficult to model with polynomials or interactions.  The spline effects are typically not interpretable, but the non-spline terms in a regression model still retain their interpretability.  To add a spline to the lognormal and gamma models, we used a Generalized Additive Model (GAM) object and created a spline of the predicted probability of B_tgt using the New data item menu and choosing Spline effect:     We let SAS Visual Statistics do the hard work of iteratively finding a spline function for our generalized additive model using the default settings.  For more information on creating generalized additive models in SAS programmatically or using SAS Model Studio, see Beth Ebersol’s post here: https://communities.sas.com/t5/SAS-Communities-Library/How-to-Use-Generalized-Additive-Models-in-SAS-Viya/ta-p/895086.   The lognormal and gamma models used backwards elimination with variables removed based on the validation average squared error.  The GAM models used the same predictors as the corresponding generalized linear models except the stage-1 predictions were removed from the Continuous effects role and a spline function of these stage-1 predictions was added. We also tried a decision tree as a stage-2 model.  The default decision tree settings were used, and the validation Average Squared Error was compared for each model:   Model Validation ASE Generalized linear model: normal distribution, log link 55,646,643 GAM: normal distribution, log link 55,539,930 Generalized linear model: gamma distribution, log link 56,074,933 GAM: gamma distribution, log link 60,839,972 Decision Tree (default settings) 51,307,883   The decision tree performed the best of the 4 models tested. The predictions from the decision tree model were saved using the Derive predicted option.  We then created a list table with the customer account IDs, the predictions from the stage-1 logistic regression, the predictions from the stage-2 decision tree, and the expected future profit (the product of the predictions).  Finally, we sorted the table by expected future profit:     To use the two-stage model with new data, the Export model option could be used for both models to produce SAS data step scoring code which could be implemented in the SAS Studio interface.   Summary and more information   This was a brief example of two-stage modeling using SAS Visual Statistics. Would you like to know more? Generalized linear models are a great tool to expand your modeling skills beyond linear regression.  These models are discussed in the course Statistics 2: ANOVA and Regression and the SAS Viya course SAS Visual Statistics: Interactive Modeling Building.  Generalized Additive Models are addressed in the course Regression Methods Using SAS Viya.  Hope this post has given you some new ideas about how to model your data.   See you at the next SAS class!     Find more articles from SAS Global Enablement and Learning here. ... View more

Do you frequently use binary logistic regression? If you do, you may have encountered a problem called quasi-complete separation. This problem is often associated with categorical variables with many levels and it can cause the parameter estimates and p-values for your model to be untrustworthy. In this post, I’ll explain the problem and my favorite approach to dealing with it: converting the categorical predictor to a continuous predictor using smoothed weight of evidence coding.   So, what is quasi-complete separation? Before describing quasi-complete separation, it may be helpful to understand complete separation in logistic regression. For this example, I’m using the AmesHousing3 data set from the SAS course Statistics 1: Introduction to ANOVA, Regression, and Logistic Regression (these data are a modified subset of the data from De Cock 2011). In these data, each observation is a house that was sold in Ames, Iowa between the years 2006–2010. The predictor variables in these data are descriptors of the houses such as the basement area in square feet (Basement_Area) and the number of bedrooms (Bedroom_AbvGr). The response variable Bonus is a binary variable with the value 1 indicating that the realtors selling houses were eligible to receive a bonus for the sale and 0 indicating they were not.   Using PROC LOGISTIC to fit the model for Bonus with the predictors Basement Area and SalePrice of the houses results in the following warning in the log and convergence status message in the results:   proc logistic data=STAT1.AmesHousing3; model Bonus(event='1')= Basement_Area SalePrice; run;   WARNING: There is a complete separation of data points. The maximum likelihood estimate does not exist.   WARNING: The LOGISTIC procedure continues in spite of the above warning. Results shown are based on the last maximum likelihood iteration. Validity of the model fit is questionable.   Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   This is accompanied by unusual p-values for the overall model. Assuming alpha=0.05, the likelihood ratio p-value below says the model is highly significant (p<0.0001), but the Wald chi-squared p-value is at least 500X bigger than Likelihood Ratio p-value (p=0.0504).     Have you ever seen this complete separation problem come up in your own logistic regression analyses? Possibly not, because this is actually a rare problem. The related problem “quasi-complete separation” is much more common. With complete separation, SAS is telling us that the maximum likelihood estimate does not exist. What caused this?   Well, I caused it. I was able to create the complete separation problem by running a binary logistic regression problem with a “perfect” predictor. In the AmesHousing3 data set, Bonus =1 when the price of a house (SalePrice) is greater than $175,000 otherwise Bonus=0. So, SalePrice completely separates the observations into the two possible outcomes for the response: bonus eligible and bonus ineligible. Complete separation means that there is some linear combination of the explanatory variables that perfectly predicts the dependent variable. When this happens, the maximum likelihood algorithm does not converge and the results of the model cannot be trusted.   Quasi-complete separation   Quasi-complete separation is very similar to complete separation. When a continuous predictor nearly completely separates the binary response into different categories, it can produce quasi-complete separation. But the more common cause is having a categorical predictor with levels that have all events or non-events for each case. Let’s see an example of this.   Now I’m changing my logistic regression model to have Bedrooms_AbvGr (bedrooms above ground) as the sole predictor. This categorical variable has 5 levels, corresponding to 0–4 bedrooms.   proc logistic data=STAT1.ameshousing3; class Bedroom_AbvGr/param=ref ref=first; model Bonus(event='1')= Bedroom_AbvGr; run;   Here is the warning from the log and the model status message in the results:   WARNING: There is possibly a quasi-complete separation of data points. The maximum likelihood estimate may not exist.   WARNING: The LOGISTIC procedure continues in spite of the above warning. Results shown are based on the last maximum likelihood iteration. Validity of the model fit is questionable.     Partial output:     Bedrooms is not a significant predictor but with the quasi-complete separation warning, the parameter estimates, odds ratios, and p-values can’t be trusted. Even before we ran the model, we could tell that there will be a problem when we look at a cross-tabulation of Bonus and Bedrooms:     Of the houses with 0 or 4 bedrooms above ground, 0 were bonus eligible. This is problematic for logistic regression because the logit for each level is calculated as the natural logarithm of (# bonus eligible houses/# bonus ineligible houses) and the natural log of zero is undefined. This could make the parameter estimates invalid for the levels with zero events (i.e., 0 bonus eligible houses). Even worse, if one of these levels is the reference level for the parameter estimates (the last level 4 is the default but I used level 0) then the parameter estimates for all levels will be affected. The same problem would occur if a level had no bonus ineligible houses because calculating the logit would involve dividing by zero which of course is also not defined.   Are the levels with zero events really perfect predictors? No. Small sample size in a level can result in zero events or non-events by chance. Perhaps with a larger sample size for say 4-bedrooms, there might be some bonus eligible houses. This quasi-complete separation problem can come up with small sample sizes or when categorical predictors have many levels, since there are more chances to have zero events or non-events. How can we avoid this problem?   One way to remove the quasi-complete separation problem from your logistic regression is to convert the categorical predictor into a continuous predictor. The continuous predictor can be calculated using smoothed weight of evidence (SWOE) coding. SWOE is similar to averaging the log odds in a particular level with the log odds of the overall sample. Let’s see the formula for calculating smoothed weight of evidence.   If we replaced the levels of bedrooms with the natural log of the odds of being bonus eligible (i.e., the logit), we would be using “weight of evidence” (WOE) coding:     WOE coding replaces each level of the categorical predictor with the logit, a continuous variable. But the logit is undefined when there are zero events or non-events. So we need to modify WOE to get rid of zeros in the numerator and denominator. One approach is to add a small number to both the numerator and the denominator (say 0.5) to remove any zeros. Another approach is to calculate smoothed weight of evidence (SWOE):     where ρ 1 = the proportion of events in the entire sample and c is a smoothing parameter. Since ρ 1 /(1-ρ 1 ) is the odds of the event for the whole sample, the equation for SWOE is like a weighted average of the logit for each level and the logit for the overall sample. Higher values of c reduce the variability of the logits for each level and make them more similar to the logit for the overall sample.   What’s a good choice for the smoothing parameter c? This is an empirical question that depends on your specific data set. Several values of c can be tested to see which produces the SWOE that has the strongest relationship with the response variable. Ideally, c could be found using a hold-out validation data set to reduce overfitting.   I tried a few values of the smoothing parameter (1, 5, 10, 15, 25) and found c=1 produced the best SWOE. Bear in mind that AmesHousing3 is a small data set (n=300) and I did not use a hold-out validation data set to find c. Here’s the code I used to calculate SWOE and add it to the data. The PROC MEANS step creates a data set (Counts) that has the number of bonus eligible houses (events) and the total number of houses (_FREQ_) for each level of Bedrooms. In the Data step, I calculate SWOE (bedrooms_swoe) using ρ 1 =0.176 and the smoothing parameter c=1.   proc means data=STAT1.AmesHousing3 sum nway; class Bedroom_AbvGr; var Bonus; output out=counts sum=events; run; %let rho1= 0.176; *45 bonus eligible/255 bonus ineligible houses; %let c=1; data counts; set counts; nonevents=_FREQ_-events; bedrooms_swoe=log((events+&c*&rho1)/(nonevents+&c*(1-&rho1))); run; proc print data=counts noobs; var Bedroom_AbvGr events nonevents bedrooms_swoe; run;   The PROC PRINT output:     Then to add SWOE to the data and fit the model for Bonus with the continuous predictor, I used the following code:   data AmesHousing3; set Stat1.AmesHousing3; if Bedroom_AbvGr=0 then bedrooms_swoe=-2.772588722; if Bedroom_AbvGr=1 then bedrooms_swoe=-1.393311934; if Bedroom_AbvGr=2 then bedrooms_swoe=-1.676328118; if Bedroom_AbvGr=3 then bedrooms_swoe=-1.747823684; if Bedroom_AbvGr=4 then bedrooms_swoe=-3.912023005; run; proc logistic data=AmesHousing3; model Bonus(event='1')=bedrooms_swoe; run;   Partial output:       The model converged and bedrooms_swoe has more reasonable parameter estimates, odds ratios, and p-values.   Not only does SWOE coding remove any quasi-complete separation problem, it also reduces the number of parameters being estimated in the model. In these models, 4 parameters are being estimated for Bedrooms_AbvGr but only 1 is estimated for bedrooms_swoe. What if instead of a 5-level categorical predictor, we had a 1,000 level categorical predictor? SWOE coding would convert this to a continuous predictor and still require only 1 parameter estimate instead of nearly 1,000.   Are you interested in learning more about SWOE and other methods for dealing with quasi-complete separation? Then consider taking the SAS class Predictive Modeling Using Logistic Regression. Not only will you learn about other approaches to dealing with the quasi-complete separation problem, this class is great preparation for attaining the SAS Statistical Business Analyst credential.   References   De Cock 2011. “Ames, Iowa: Alternative to the Boston Housing Data as an End of Semester Regression Project.” Journal of Statistics Education Volume 19, Number 3 ... View more

A question about model hierarchy came up recently while I was teaching the SAS class Predictive Modeling Using Logistic Regression.  In class, we were using PROC LOGISTIC to investigate all 1.12 quadrillion possible statistical models given 50 effects by using the BEST=SCORE option.  This option produced a list of models of varying sizes, from 1 effect to 50 effects, showing the single best model of each size ("best" meaning having the highest score chi-squared). The 3-effect model with the highest score chi-squared statistic was y= X 1 + X 2 + X 3 *X 4 , but we ignored it because it is not a hierarchically well-formulated model.  Not everyone was familiar with model hierarchy and why it is important.  So, in this post, I’ll explain why it’s best practice to maintain model hierarchy in your regression models.   What does it mean to maintain model hierarchy?   Maintaining model hierarchy means regression models with interactions (e.g., X 1 *X 2 ) include the corresponding main effects (e.g., X 1 and X 2 ) even if they are non-significant.  Model hierarchy also applies to polynomial models in that a model with higher order polynomial terms (e.g., X 3 ) includes all lower order terms (e.g., X 2 and X) even if these are nonsignificant. Models that maintain model hierarchy are sometimes referred to as “hierarchically well-formulated models”.   So model hierarchy is maintained in these models:   Y= A + B + C + D + C*D   Y= A + B + C + B*C   Y= X 3 + X 2 + X   Y= X 2 + X   but not these models:   Y= C*D   Y= A + B + B*C   Y= X 4   Y= X 3 + X 2     Why is it best practice to maintain model hierarchy?   There are several reasons to use hierarchically well-formulated models in your business or research.  In brief, (1) they allow fine tuning of the fit of the model to your data, (2) omitting the lower order terms makes strong (often incorrect!) assumptions about relationships in your data, (3) they make the interaction effect inference invariant to coding, and (4) not maintaining hierarchy opens the door to absurdly complex models.   Maintaining hierarchy allows fine-tuning the fit of the model to your data   Often, we’re trying to fit a regression model to our data without having strong theory telling us what the model should be.  We see a Y vs X scatter plot showing a straight pattern, so we include the variable X. If we see a bend in the pattern, we include X 2 . But while X 2 captures the curvature, including X provides additional useful information about the shape of the response function.  For example, in this picture below, we can see how changing the coefficient for X in y= X 2 + bX changes the shape of the X-Y relationship:    Modified from Wikipedia https://commons.wikimedia.org/wiki/File:Quadratic_equation_coefficients.png)   Select any image to see a larger version. Mobile users: If you do not see this image, scroll to the bottom of the page and select the "Full" version of this post.   Keeping bX in the model provides greater flexibility to the model shape. So, including the lower order terms will likely allow better fit of your models than leaving it out.   Omitting the lower order terms makes strong assumptions about relationships in your data   Let’s start off with a regression model with 2 predictors and their interaction:   y=b 0 + b 1 X 1 +b 2 X 2 + b 3 X 1 X 2 + e   When dealing with an interaction model such as the one above, we are often taught to not interpret the coefficients b 1 and b 2 when the interaction is significant.  This is because a significant 2-way interaction means the effect of one predictor on the response depends on another predictor.  So it doesn’t make sense to interpret either of the interacting predictors alone. While this is true, the coefficients b 1 and b 2 do have meaning in this interaction model, and it is not the average effect of the predictors. The meanings of the coefficients are:   b 1 = the effect of a unit-increase in X 1 when X 2 =0   b 2 = the effect of a unit-increase in X 2 when X 1 =0.   [Note: the above interpretation only makes sense when X 1 and X 2 have true zero points (i.e., are on ratio scale) as in measures of size or distance. Interval scaled continuous variables have arbitrary zero points, so it would not be meaningful to interpret the effect of one variable while another variable was set to its arbitrary zero point.]   Leaving out X 1 and X 2 from your model is equivalent to insisting that b 1 and b 2 both equal zero.  This is a strong assumption that X 1 has zero effect on Y when the other X 2 equals zero (yet X 1 is still important to the interaction).  Even if the effect is not statistically different from zero, if it is not exactly zero, leaving it out will bias your other parameter estimates. Remember that a non-significant statistical test does not mean that the null hypothesis is true.  So, to sum it up, this assumption shouldn’t be made without a justification for it.   When is it justified? I found it difficult to think of an example when this would be reasonable but here goes: you know those epoxy adhesives that come in two tubes?  When the two liquids are combined, they form a strong epoxy adhesive that will glue practically anything to anything.  A model of adhesive strength (y) that included the amount of epoxy components (X 1 and X 2 ) could have zero coefficients of the main effects (since they are inert alone) but still have a strong interaction effect.  Without strong justification, you’re better off maintaining model hierarchy and keeping the main effects in your model.   Keeping lower order terms makes interaction effect inferences invariant to coding   Model hierarchy is desirable because models that are hierarchically well formulated have inferences that are invariant to the coding that you choose for your predictors.  One common way to recode predictors is to center them.  Centering means subtracting the average from each value so that the transformed variable has a mean of zero. There are a variety of reasons to center data including reducing collinearity when interactions are present and making the intercept term more meaningful.  Let’s see how centering affects inferences for models that do and don’t maintain model hierarchy.   For this example, I’ll use the data set sashelp.cars.  I’ll use PROC STDIZE to create centered versions of the predictors Cylinders and MPG_City and save them to an output data set called centered. Then I’ll create 4 regression models, each using the response variable MSRP.  The predictors for the first model are Cylinders, MPG_City, and their interaction.  The second model is a version of the first model that uses centered variables.  Model three is not a hierarchically well-formulated model—it only has the Cylinders * MPG_City interaction and is missing main effects. The fourth model is the same as model three, except that centered variables are used.   The code used is shown below.  What happens to the coefficient of determination (R 2 ) and the p-values for the interactions when the variables are centered?   proc stdize data=sashelp.cars out=centered (keep=msrp Cylinders MPG_City centered:) sprefix=centered_ oprefix; Var Cylinders MPG_City; run; *hierarchically well-formulated model; proc glm data=centered; model msrp=cylinders|mpg_city; run; *hierarchically well-formulated model (centered variables); proc glm data=centered; model msrp=centered_cylinders|centered_mpg_city; run; *no model hierarchy; proc glm data=centered; model msrp=cylinders*mpg_city; run; *no model hierarchy (centered variables); proc glm data=centered; model msrp=centered_cylinders*centered_mpg_city; run;   Here is the output from models 1 through 4.  Focus your attention on the R 2 and the interaction p-values from the Type 3 SS tables:     Model 1     Model 2     Model 3     Model 4     Centering the variables does not affect the interaction p-value or the model R 2 for the hierarchically well-formulated model.  But in the interaction models without main effects, centering the predictors changes the R 2 from 0.07 to 0.005 and changes the interaction p-value from p=0.01 to p=0.14. This is not a desirable outcome!   We can see why omitting constituent terms from an interaction is problematic if we can tolerate a little algebra.  Let’s start with the model y= b 0 +b 1 X 1 +b 2 X 2 +b 3 X 1 X 2 , then add constant c to X 1 and constant d to X 2 .   y= b 0 +b 1 (X 1 +c) +b 2 (X 2 +d) +b 3 (X 1 +c) (X 2 +d)     = b 0 +b 1 X 1 +b 1 c+b 2 X 2 +b 2 d+b 3 X 1 X 2 +b 3 dX 1 +b 3 cX 2 +b 3 cd   Next, collect all the constant terms (everything that does not involve an X), then all the functions of X 1 , functions of X 2 , then functions of X 1 X 2 :     = b 0 + b 1 c + b 3 cd + b 2 d + b 1 X 1 + b 3 dX 1 + b 2 X 2 + b 3 cX 2 + b 3 X 1 X 2   and simplify:     =(b 0 + b 1 c  + b 3 cd + b 2 d) + (b 1 + b 3 d)X 1   + (b 2 + b 3 c)X 2   + (b 3 )X 1 X 2   The sums in parentheses are an intercept, the coefficient for X 1 , the coefficient for X 2 and the coefficient for the X 1 X 2 interaction using the transformed variables. We can fit this model the same way we fit the original model y= b 0 +b 1 X 1 +b 2 X 2 +b 3 X 1 X 2 , since it has 2 main effects and their interaction. Compare this to the model y= b 3 X 1 X 2 (with no main effects) when we add the constants c and d to X 1 and X 2 respectively:   y= b 3 (X 1 + c)(X 2 + d)     = (b 3 cd) + (b 3 d)X 1 + (b 3 c)X 2 + (b 3 )X 1 X 2   The parentheses were added in the last step above to make it easier to see that this is a regression with an intercept, a main effect of X 1 , a main effect of X 2 and in X 1 X 2 interaction effect.  But we can’t estimate this model if we try to fit y= b 3 X 1 X 2 because it has no main effects.  Adding a constant in the interaction model without hierarchy creates linear effects of the predictors that aren’t captured by the model being fit. These main effects aren’t being estimated and will get lumped in with the interaction.  As we’ve already seen, adding a constant will change the R 2 and p-values for then non-hierarchical model.  We don’t want this!   The above explanation can be found in several papers such as Paul D. Allison’s (1977) “Testing for Interaction in Multiple Regression”. This is the same Allison who wrote my favorite books on logistic regression and survival analysis (Logistic Regression Using SAS: Theory and Application, Second Edition and Survival Analysis Using SAS: A Practical Guide, Second Edition).  Check out Allison’s paper for a more thorough explanation of model hierarchy and invariance to coding.   Not maintaining model hierarchy opens the door to absurdly complex models   Let’s look back at the p-value for the model 1 interaction from the type 3 sums of squares table.  The p- value is 0.0119 and if we’re using alpha=0.05 as our significance threshold, this tells us that the interaction is statistically significant.  But remember what type 3 (aka partial) sums of squares are showing.  It shows the additional sums of squares explained by an effect when added to a model containing all the other model effects.   In other words, the p-value tells us if the effect is significant after explaining as much as possible using all the remaining variables. So, the significant Cylinders*MPG_City interaction is telling us that after accounting for the main effects, the interaction is still important.   What if we leave out the main effects, as in model 3? The significance of the interaction is no longer assessed after accounting for main effects.  If model 3 is all we look at, we miss knowing whether a simpler, main effects only model could have explained the data.   Researchers normally follow the principle of parsimony or Occam’s Razor which basically says to use the simplest model that gets the job done. There’s no sense in using the model msrp=MPG_City 5 even though it is statistically significant (p=0.029) when the hierarchically well-formulated quadratic model msrp= MPG_City + MPG_City 2 (p<0.001) gets the job done and results in all higher order terms being nonsignificant.  If we use a 5 th degree polynomial without the lower order terms, why not use y= X 17 ? The model y= X 17 is needlessly, even absurdly complex when a much simpler model does the job.  The same logic applies to modeling an interaction without the main effects.   Now you know why to use hierarchically well-formulated models. If what you read here was intriguing, there are several places to go for more information. For learning effective model selection approaches with PROC LOGISTIC, check out the course Predictive Modeling Using Logistic Regression. For a deeper understanding of the do’s and don’ts of statistical inferences in regression modeling try the course Statistics 2: ANOVA and Regression. For a foundational course on many of the same topics, try Statistics 1: Introduction to ANOVA, Regression, and Logistic Regression.  If you are interested in a foundational statistics course geared towards predictive modeling and machine learning, try out Statistics You Need to Know for Machine Learning. Live Instructor Dates for these classes, and more, can be found by using the links provided.   See you in the next SAS class!   References   Allison, Paul D. 1977 “Testing for Interaction in Multiple Regression” American Journal of Sociology, Vol. 83: 1, pp. 144-153   Allison, Paul D. 2010. Survival Analysis Using SAS ® : A Practical Guide, Second Edition. Cary, NC: SAS Institute Inc.   Allison, Paul D. 2012. Logistic Regression Using SAS ® : Theory and Application, Second Edition. Cary, NC: SAS Institute Inc.     Find more articles from SAS Global Enablement and Learning here. ... View more

Most predictive modelers are familiar with data splitting and honest assessment.  When developing a predictive model, data is typically split into training, validation, and test sets for developing and assessing models.  Despite having knowledge of partitioning, many modelers make the mistake of imputing in a manner that allows data to leak from the validation and test sets into the training data. In this post, I’ll describe the problem of data leakage, explain how to impute to avoid leakage, and give advice for using SAS tools that make imputation easy to accomplish as well as some options to avoid when preparing data and modeling.   Data partitioning   First, let’s review data partitioning.   The goal of predictive modeling is generalization, that is, to make predictions or score new data.  So predictive modelers need to build models that function well at making predictions on new data sets.  Data splitting is what best allows one to evaluate a model’s ability to accomplish this goal.  When developing a predictive model, the researcher will typically split their data into a set for training the model to make predictions (training), a data set for choosing a model out of the candidates (validation), and a data set for final unbiased assessment of the champion model’s performance (test).  Without data splitting, we can only assess how well the model predicts the historic data used to train the model.  When we assess a model on the training data, we get an optimistically biased assessment of performance.  The model will have been chosen because it does well on the particular data set used for training, which may have features not shared by data to be scored later. That’s not what we want.  We want to know how accurately the model functions, how well we can expect the model to accurately make predictions when applied to previously unseen data.  Data splitting allows honest assessment, that is, an unbiased assessment of the chosen model’s performance. It allows the modeler to simulate applying a model built today to new data seen tomorrow. For this reason, the modeler must ensure that any information in the validation and test sets is excluded from any type of fitting and pre-modeling processing steps such as imputation and transformations.   For a review of data partitioning using SAS Visual Analytics/Visual Statistics, see Beth Ebersol’s blog [1].  For steps on partitioning data using DATA step code and PROC SURVEY SELECT see Rick Wicklin’s blog [2].   In some cases, researchers use only training and validation sets, omitting the test set.  When this occurs, it is with the understanding that the validation data model assessment measures should be considered the upper bounds of model performance.  For simplicity, I’ll consider data that has been partitioned into only training and validation for the rest of this post. I’ll also assume mean imputation is the chosen imputation approach.   Data leakage   Data leakage, also called information leakage, occurs when information from the validation data contaminates the training data.  This information that is introduced during training would not be accessible in a real-world application of the predictive model.  This is different than a related issue, sometimes called target-leakage, in which information about the target is used as a predictor. Data leakage can lead to optimistically biased assessments of model performance, possibly making unacceptable model performance appear satisfactory.  If any information from the validation data is used during training, the model gains an unfair advantage that won’t exist when it is deployed in a real-world scenario.  To avoid data leakage, all data preprocessing including imputation of missing values is done solely using the training data.  This simulates as closely as possible the conditions under which the model will be used in practice.   So how does leakage actually occur?  Often a researcher will impute missing values prior to partitioning their data. This seems reasonable at first because imputing after partitioning requires imputing twice instead of once.  But this results in leakage, with the model being trained on data that contains information that the model should not be able to see. It is effectively cheating, showing the model some of the answers it is trying to predict.  Kapoor and Narayanan [3] show that several kinds of data leakage are surprisingly common in published literature across several scientific disciplines.   The best practice for imputing to avoid leakage is as follows.  First, split the data into training and validation sets.  For mean imputation, find the training data means, then impute the training means for missing values in both the training and validation sets.   For example, if a data set of 1 million observations was partitioned into 50% for each training and validation, and the variable means were:   Data set Whole data (n = 1M) Training (n = 500K) Validation (n = 500K) Variable mean 75 80 70   we would use 80 to replace missing values in both the training data and the validation data.  When researchers impute missing values before partitioning, they are using 75 for the whole data set.  This means that the model has access to information in the training data that it would not have in practice.  The model is cheating by seeing information about data that should be unseen.   After partitioning then imputation, the training and validation sets may need to be recombined.  Many SAS Viya statistical and machine learning modeling procedures will compute fit statistics on all data partitions if contained within the same data set.   How to impute while avoiding leakage   There are several ways to impute in SAS, and the two methods that I use most are shown below.  PROC VARIMPUTE is a nice procedure for imputation using in-memory data in SAS Viya.  If you’re using SAS 9, PROC STDIZE will do the job.  Both allow you to impute the same values on the validation data that were used for imputing the training data.   PROC VARIMPUTE can be used easily using the point-and-click SAS Studio Imputation task. You can find it located under the SAS Viya Prepare and Explore Data tasks. After identifying variables to impute and the data sets involved, the task-generated code will need to be edited to include a CODE statement.  This will produce SAS code that can impute the training means onto other data such as the validation or test sets.  The SAS code produced then can be referenced in a data step using %include to impute the training means. PROC VARIMPUTE can also be used to impute medians, random numbers between the minimum and maximum values, or custom values.   proc varimpute data=mycas.training; input var1 var2 var3 / ctech=mean; output out=mycas.training_imputed copyvars=(_all_); code file='/home/student/imputed_vars.sas'; run; data mycas.validation_imputed; set mycas.validation; %include '/home/student/imputed_vars.sas'; run;   PROC STDIZE has an OUTSTAT option which produces a data set containing location and scale measures as well as other computed statistics. The REPONLY option in the code below replaces missing values with the variable means, without altering non-missing values.  A second PROC STDIZE step gets used with the METHOD=IN() option to replace missing validation data values with the training data means. PROC STDIZE has many additional options for transforming and imputing values.   proc stdize data=work.training method=mean reponly outstat=training_means; var var1 var2 var3; run; proc stdize data=work.validation reponly method=in(training_means) out=work.validation_imputed; var var1 var2 var3; run;   Further recommendations   Many of the SAS Viya regression procedures have an informative missingness option that can automatically apply mean imputation for missing values.  These procedures include GENSELECT, LOGSELECT, PHSELECT, QTRSELECT, and REGSELECT and this option can be applied by specifying the INFORMATIVE option in their MODEL statements.  Applying the informative missingness option does 3 things.  First, it imputes the mean for missing continuous variables. Second, it creates missing indicator variables which can potentially capture the relationship between missingness and the target.  Third, it treats missing values for categorical predictors as a legitimate level for analysis. These procedures also have PARTITION statements, which allow random partitioning by the procedure or specification of a partition indicator variable if one is already present in the data. Using PARTITION results in fit statistics being calculated on both training and validation sets as well as allowing a champion model to be chosen automatically based on validation data performance.   I recommend not using the PARTITION statement and INFORMATIVE option with the same PROC step for these procedures. The informative missingness option will use the whole data set means for imputation, then randomly partition the data into training and validation sets.  This leads to leakage as explained earlier.  Instead, use the informative missingness option when working only with training data, not a data set that contains all partitions.  When I use the PARTITION statement to identify partitions and get fit statistics on each partition, it is with data that has already been imputed, using the training means for the whole data set.   In summary, for predictive modeling, we don’t want the unconditional variable means imputed for the whole data set.  Instead, we want the training data means imputed for the whole data set. Imputation needs to be done after partitioning, not on the unpartitioned data.  Also, in several SAS Viya regression procedures, it is best to avoid using the PARTITON statement and the INFORMATIVE option together.  The INFORMATIVE option makes sense if modeling the training data only.   Learning more   Are you interested in learning more about data preparation for predictive modeling including partitioning and imputation? Then consider taking the SAS 9 class Predictive Modeling Using Logistic Regression, which covers all the topics mentioned here (and many more!) in more detail. Not only will you learn about the care needed in data preparation, this class is great preparation for attaining the SAS Statistical Business Analyst credential (https://www.sas.com/en_us/certification/credentials/advanced-analytics/statistical-business-analyst.html).  If you’re modeling using SAS Viya, Supervised Machine Learning Procedures Using SAS Viya in SAS Studio is a great option too.   See you in the next SAS class!   References   [1] Beth Ebersol 2019 “Training, Validation, and Testing for Supervised Machine Learning Models” Training, Validation, and Testing for Supervised Machine Learning Models (sas.com)   [2] Rick Wicklin 2019 “Create training, validation, and test data sets in SAS” Create training, validation, and test data sets in SAS - The DO Loop   [3] Kapoor and Narayanan  (2022) “Leakage and the Reproducibility Crisis in ML-based Science https://doi.org/10.48550/arXiv.2207.07048     Find more articles from SAS Global Enablement and Learning here. ... View more