Suppose you own a small grocery store. You have recently heard a few rumors that people don’t want to shop at your store because they have to wait too long to check out. But hiring more cashiers is an expensive option. And how many extra cashiers should you hire? If you hire too few, you might lose customers. If you hire too many, you will lose money. You need to find a way to determine the best number of cashiers to hire. One way to handle this situation is to simulate having more cashiers. But a monte carlo simulation is not appropriate since the system is not stable. This process changes at a countable number of points in time called events like every time a customer enters the store. For example, the distribution of waiting times for a cashier for the second customer is different than the distribution of waiting times for the first customer. Discrete event simulation is designed for just these types of situations. SAS® Simulation Studio is a discrete event simulation tool that can be used to model a system as it changes over time. This post will discuss, at a high level, how to create a discrete event simulation model of a simple grocery store.   To start creating a discrete event simulation model, creating a process flow diagram is a good first step. This diagram illustrates all of the steps of the process that are relevant to the purpose of the study. In this case, we are interested in how long people wait for a cashier as well as how many cashiers are needed for the store. The process, in words at a high level is as follows:   The store opens at 8:00 am and closes at 6:00 pm. A customer enters the store. The customer shops for some time. The customer enters a checkout queue. The single cashier scans the purchases of the customer and collects the money. The customer leaves the store.   In order to create the process flow diagram, some more information is needed. For example, how often do customers come to the store? How long do they shop? How long does it take for the cashier to scan the purchases and collect the money? The answers to these questions are typically answered by using a statistical distribution. Typically, the distributions are determined by examining past data. For our example we will use the following: Customers enter the store according to an exponential distribution with a mean of 4 minutes. Customer shopping follows a uniform distribution from 5 minutes to 30 minutes. The checkout time follows an exponential distribution with mean of 3 minutes.     With the diagram in hand, we need to discuss a little bit of terminology. For discrete event simulations the item (or items) that moves through the process is called an entity. Entities can have properties or attributes such as the age of the customer. The entities need to be generated according to a distribution of interarrival times rather than the number of entities. An exponential distribution, which is commonly used for interarrival times, will be used with a mean of 4 minutes. While the customer is shopping, that appears as a delay in our simulation since the customer is not using any of our services. The customer then waits in a queue for the cashier. In fact, the cashier is a resource that is needed by the customer to complete their shopping. Resources control the flow of the entity through the process. In our case, the customer is delayed by using the resource (cashier). That delay or checkout time is called Service in the diagram also uses an exponential distribution with a mean of 5 minutes. Finally, we are interested in two outputs from this model: the customer’s average time waiting for the cashier and the utilization of the cashier. In other words, how busy is that cashier?   Building the model Simulation Studio involves dragging blocks to build the model and using links to connect the ports that are attached to each block. A block is the most fundamental unit of the discrete event simulation model. Each block has a well-defined and specialized function. There are ports on each block that allow communication between the blocks. Red ports are on the left and right side of the blocks that allow entities to pass through the block. Blue ports are on the top and bottom of the blocks and allow information or data to pass through the block. The links define the path of the flow of entities and/or information. You end up with a graphical picture for a model. We won’t be able to go through all of the details here, but we will build our simple grocery store model. Simulation Studio is based around projects. When you create a new project, you are presented with this interface.     You drag blocks from the templates into the model window to build the model. For this model we will only need six blocks.     We will start by dragging an entity generator block into the model to generate the customer entities. We attach a numeric source block to specify the exponential distribution with a mean of 4. After the entity generator block, the customer moves to a delay block to indicate their time spent shopping. Another numeric source block is attached to the delay block to specify the uniform distribution between 5 and 30 minutes. Next is a queue block which shows the customer entering the line for the cashier. Note that the queue block is only used if the cashier is already busy. Otherwise the customer passes directly to the server block which indicates the check out time with the cashier. Yet another numeric source block is used to indicate the distribution of checkout times. Finally, the customer leaves the store as indicated by the disposer block. The randomness of this model is all handled by the numeric source blocks which are pulling random numbers from the specified distributions. The final model looks as follows.     The details of this model, such as the different distributions that are used, and their parameters are all specified in the numeric source blocks. Similarly, the server block specifies the number of cashiers working in the store. A few additional modifications are made to this model so that we can see some of the statistics that we are interested in observing. I have added two number holder blocks which are blocks that will store the maximum queue length and the average amount of time spent waiting in the queue. I ran the simulation, closing the store doors at 6:00 pm and allowing the customers in the store to finish their shopping.     The results of running this simulation show that the maximum number of people waiting for the cashier was 13. We also see that on average people waited 19.9 minutes for the cashier. Finally, once we closed the doors (which occurred after 600 minutes), the system clock in the upper-right of the screen capture displays 652.02004. This means that it took the cashier just over 52 minutes to checkout the people in line after the store closed. Now these results are alarming, but this was a simulation. Maybe the next day would look better? So, let’s run the simulation for 10 days. This is accomplished by entering 10 replicates in the Experiment window. You will also see that some other responses have been added like cashier utilization.     Running the simulation for 10 days shows that the 19.9 minute average waiting time was indeed a bit unusual. On average across all of the days (the top line of the table), customers are waiting 9.6577 minutes. Maybe things are not as bad as they first seemed. But what if we did hire a second cashier? Running the simulation for 10 days with two cashiers results in lower waiting times, as expected.     The average wait time is now under a minute. However, the cashier utilization is now 0.38 meaning that cashiers are only busy just over 38% of the time. Two cashiers are not really needed. Perhaps better customer service can be achieved by having one cashier and maybe a manager to come and help checkout people during rush hours. Or possibly have a part-time cashier to come in during the busiest times. We could modify our model to handle either of these situations, if we want.   Hopefully this simple example illustrates the value of discrete event simulation. Discrete event simulation is a tool that can be used to model real-life scenarios that have a countable number of changes to the system over time. Simulation Studio is a flexible tool that can be used to allow analysts to make decisions on many hypothetical situations by creating a model of real-world scenarios.     Find more articles from SAS Global Enablement and Learning here. ... View more

  I had an analyst contact me the other day. He was conducting what seemed to be a straight-forward one way analysis of variance (ANOVA) and ran into a problem. He requested the LINESTABLE option and saw this result:   Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   These results concerned him, and he was confused on how to best explain this. I will discuss some ways to understand and explain these results in this post.   The analyst in this situation had a small dataset. Although I do not have the data, we can make a small example that will yield similar results. This program is one simple way to illustrate the issue:   data example; input group $ response; datalines; Level_1 46 Level_1 50 Level_1 54 Level_2 56 Level_2 60 Level_2 64 Level_3 66 Level_3 70 Level_3 74 ; proc print data=example; proc anova data=example; class group; model response=group; means group/tukey lines linestable; run;   Notice that I have added the LINES and LINESTABLE options to the means statement. The LINES options will provide a graphical representation of the results while the LINESTABLE option will provide the results in a tabular form. Having an example similar to this, the analyst correctly noted that there was a statistically significant difference between the means of the three groups of data, as indicated by the ANOVA table.     To determine where the differences were, a Tukey multiple comparison test was conducted by using the MEANS statement. This led to the connecting letters report from the LINESTABLE option shown earlier and repeated here.     I also added the LINES statement that shows this information in a graphical format since some people prefer a picture to a table with letters.     The analyst said that he could clearly see that Level_3 was different than Level_1. But Level_2 was causing some concern. As he stated: "If I had only studied Level_1 and Level_2 (or Level_3 and Level_2), I would conclude that no differences were found and that Level_2 observations are indistinguishable from Level_1 (or Level_3) observations. The only differences are due to noise. Without testing all three levels, I would have concluded that the mean of Level_2 is the same as the mean of Level_1 or Level_3. But what if I did not have enough resources to test all three levels? Would I have been misled? Does this imply that you should always test every possible level of a grouping variable in order to achieve the correct results?"   These are very good questions, but there is a flaw in the logic that is buried in his questions. The Tukey multiple comparison test, like most statistical tests, has a null hypothesis of equality. For this situation, there are three tests performed by the Tukey procedure. The null hypotheses are   Mean of Level_1 = Mean of Level_2, Mean of Level_1 = Mean of Level_3, and Mean of Level_2 = Mean of Level_3.   Remember that for statistical hypothesis tests, the null hypothesis is assumed to be true until the data tells us otherwise. If the data does not tell us otherwise does NOT make the null hypothesis true. In other words, we do not want to conclude that the null hypothesis is true. This test can only declare groups to have different means. It does not declare groups to have the same mean.   Going back to the questions that he raised, this point was made: “…I would have concluded that the mean of Level_2 is the same as the mean of Level_1 or Level_3.” That statement is the flaw. That statement is believing the null hypothesis to be true. Our statistical test does not allow you to make that conclusion.   This test allows you to say that there is a statistically significant difference between Level_1 and Level_3 at the 95% confidence level. However, we are not able to see a statistically significant difference between Level_1 and Level_2 or between Level_2 and Level_3. There may be a difference, but we could not see it using 95% confidence. Perhaps a larger sample size is needed to better assess not only the means of each group, but also the variances of the groups. Finally, maybe one of the best ways to see that we cannot conclude that the means are equal is to create confidence intervals for the group differences. This would be the CLDIFF option on the MEANS statement. So, the MEANS statement would look like this:   means group/tukey cldiff lines linestable;   Adding this to the code will result in a table with every pairwise mean difference.     This table shows that even for Level_1 versus Level_2, the difference is likely contained in the interval -20.021 and 0.021. Although the interval contains 0, the difference could still be a large negative value. Remember that sometimes when interpreting results, it can be very easy to fall into the trap of accepting the null hypothesis. So be on the lookout for that when forming your own conclusions from your analyses.     Find more articles from SAS Global Enablement and Learning here. ... View more

Suppose you are responsible for establishing the warranty for a product. You want less than 5% of the products to fall into a warranty claim situation. You know the target variable is related to several possible predictors that you can use to estimate how often a claim may occur. How do you handle such a situation? You could build a regression model for the target variable, but regression models predict the mean. If the regression assumptions are correct though, you could still estimate the 5th quantile to determine the right conditions for the warranty. But what if those assumptions are not correct? How can you estimate that 5th quantile to establish the warranty conditions? One possible approach is to use quantile regression. This post will talk about how to perform quantile regression, highlight some of its advantages as well as some disadvantages.   To understand a situation where quantile regression might help, let’s create some data. For this situation we will create a situation with a linear trend, but non-constant variance. This situation violates one of the assumptions of linear regression. The following code will generate 50 observations with a random input variable and a target or response variable that is generated from a linear model with a variance that is proportional to the input variable level.   data qtrreg; do i=1 to 50; predictor_x = 3 + rannor(1); target_y = 10 + 8 * predictor_x + rannor(1) * predictor_x * 2.2; output; end; run; proc sgplot data=qtrreg; reg x=predictor_x y=target_y; run; proc reg data=qtrreg; model target_y=predictor_x; run;   The graph shows the relationship between predictor_x and target_y. You can also see that the variance appears to be increasing as the value of predictor_x increases.   Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   We still fit the linear regression model to this data, but the residual plot shows the increasing variance.       If we wanted to predict the 10th quantile for when predictor_x is 4, the regression line would predict the mean. We would then need to use the assumption of a normal distribution to estimate the 10th quantile. Unfortunately, since the variance is changing, the regression assumption of a constant variance is violated.   There are several different ways to handle this nonconstant variance. One of the most common approaches is to transform the target variable. But since we are wanting to predict a specific quantile, we would need to feel comfortable with the normal distribution assumption on the transformed data, obtain the prediction for the mean, determine the quantile estimate, and transform the predicted quantile back into the original units. A different approach is to use quantile regression.   To illustrate how to perform quantile regression, we will estimate three quantiles: the 10th, 50th, and 90th. We will use PROC QUANTREG. If you wish to use a SAS Viya procedure, it would be PROC QTRSELECT. The syntax for our example is this:   proc quantreg data=qtrreg; model target_y=predictor_x / quantile=0.1 0.5 0.9; run;   By specifying three quantiles, we actually are fitting three different models. As a result, we have three sets of parameter estimates.         With three models there are three different intercepts and three different slope estimates. An advantage of quantile regression is that we do not need to make any distributional assumptions. A model is fit for each quantile, regardless of the distribution. Another advantage of quantile regression is that each model can have a different slope. You can see this on the graph.     As you can see, quantile regression does have advantages. But there are more advantages than just handling the heteroscedasticity or nonconstant variance. Quantile regression makes no distributional assumptions, is robust to outliers in the Y-space, and can provide a more complete picture of the predictor effects since a separate model is fit for each quantile.   But like any statistical procedure, there are some drawbacks. Since you are estimating quantiles, you need to have sufficient data to estimate those quantiles. In our example, we had 50 observations. What if we wanted the first quantile? There is not sufficient data to estimate to that level of precision. You can get a sense of this from the graph of our example. The extreme quantile (10 and 90) models both predict the left-most data point perfectly. Also notice that the median line is NOT between the 10th and 90th percentiles. There is not sufficient data at low values of X to adequately model that area.   Another drawback is that quantile regression, while robust to outliers in the Y-space, it is not robust to outliers in the x-space. Those points have high leverage. And finally, as you might expect, quantile regression is computationally intensive. Hopefully you can see that quantile regression is quite easy to perform, and has several advantages. As long as you can avoid the disadvantages, quantile regression can be a very useful tool to have in your arsenal of analysis techniques.     Find more articles from SAS Global Enablement and Learning here. ... View more

Any datapoint with an unusual or unexpected value can be considered an outlier. Outliers typically cause problems with any statistical analysis. But with the advent of machine learning techniques and larger datasets, do we still need to worry about outliers? The answer to this question will depend on the machine learning algorithm that is being used. In a previous post, I discussed how outliers might affect a neural network model. In this post we will explore the affect outliers have on a tree model. Specifically, we will be looking at a regression tree meaning that we will use a tree model to predict a continuous target variable.   To see how outliers might affect a regression tree, we will generate some data that follows an s-shaped curve. The data will have some error associated with it so that we would not expect any model to provide a perfect prediction, but there will only be 19 observations. These data will be group 0. Now to this shape, one of three outliers will be added to the data. The first outlier type is an outlier that follows the same trend as the rest of the data, but has a large x-value. This point will be the diamond in the plot of the data and is group 1. The second outlier type will be an observation that does not follow the trend of the data and also has a large x-value. This second type is the square in the plot and is identified as group 2. Finally, the third outlier occurs in the middle of the x-space and is far off the curve of the rest of the data. This point is identified by the triangle in the plot and is group 3.   Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   The code used to generate this data is below.   /* Creating the Original 19 Observations */ data casuser.curve; call streaminit(8675309); do i=1 to 19; x=rand('uniform')*2; error=rand('uniform')*0.025; group=0; y=0.2 + 0.6*x - 0.86*x**2 + 0.28*x**3 + error; output; end; run; /* Creating the Three Different Outlier Datasets that will be merged with the original data */ data casuser.outlier1; /* Outlier that follows the same trend as the data */ x=2.1; error=rand('uniform')*0.025; y=0.2 + 0.6*x - 0.86*x**2 + 0.28*x**3 + error; group=1; run; data casuser.outlier2; /* Outlier that does NOT follow the trend of the data */ x=2.1; error=rand('uniform')*0.025; y=0.1 + error; group=2; run; data casuser.outlier3; /* Outlier in the middle of the x-space and does not follow data trend */ x=1; error=rand('uniform')*0.025; y=0.45; group=3; run;   Now that we have some data we can start fitting some models.   I will use PROC HPSPLIT with all its defaults to build the models, but with a limited depth of 2 layers. Without a validation set, a tree will continue to split until the number of leaves is close to or equal to the number of observations. I want to limit the depth and see how the observations are distributed among the limited number of leaves. This code will limit the depth of the tree and because pruning is turned off, every tree will have a depth of 2 with a total of four leaves.   * Fitting the Model to the Original Data; proc hpsplit data=casuser.curve minleafsize=1 maxdepth=2; id x; model y=x; prune off; output out=hpsplout; run; proc print data=hpsplout; run; *Plotting the Raw Data and Overlaying Predicted Values; proc sgplot data=hpsplout; scatter X=x Y=y; scatter X=x Y=P_y / markerattrs=(symbol=asterisk color=red); run;   Although I have specified the minimum leaf size of 1, that is the default for PROC HPSPLIT. We use the ID statement to ensure that the group number and the x values are stored along with the predictions in the hpsplout file. The pruning is turned off so that every tree model will have two layers with four leaves. Finally, we plot the data along with the predictions as red asterisks using PROC SGPLOT. The PROC HPSPLIT output indicates that there are 4 leaves of the final tree and 19 observations were used.       The tree structure shows the four final leaves. For our purposes here, notice that number of observations in each leaf. They are 5, 5, 1, and 8. So other than the one leaf, which is node 5, the 19 observations are pretty close to evenly distributed. This can also be seen in the plot of the data overlaid with the predictions.     Now that we see how a regression tree fits our pattern we will add an outlier. Let’s start by adding an outlier that has a larger X value, but still follows the trend in the data, outlier 1 that was mentioned earlier.   /* Creating and Analyzing the Data with an Outlier that Follows the Trend */ title1 'Original Data Plus an Outlier that Follows the Trend'; data casuser.followtrend; set casuser.curve casuser.outlier1; proc print data=casuser.followtrend; run; ods graphics / attrpriority=none; * Plotting the Original Data and the Outlier; proc sgplot data=casuser.followtrend; styleattrs datasymbols=(CircleFilled DiamondFilled); scatter X=x Y=y / group=group; run; ods graphics / attrpriority=color; * Fitting the Model with Outlier; proc hpsplit data=casuser.followtrend minleafsize=1; id group x; model y=x; prune off; output out=hpsplout; run; proc print data=hpsplout; run; * Plotting Data with Overlaid Predictions; proc sgplot data=hpsplout; scatter X=x Y=y / group=group; scatter X=x Y=P_y / markerattrs=(symbol=asterisk color=red);; run;         Only a portion of the output is shown here, but you can see that the model has a node with only one observation. Looking at the plot of the data overlaid with the predictions, we can see that the outlier is the point that was split off on its own. By splitting that outlier into its own leaf node, the outlier does not affect the model much at all.   Perhaps if we had outlier number 2, the one with a large X value that does NOT follow the trend of the rest of the data the results will be different. The code for fitting a model to that data and some of the selected output is below.   /* Creating/Analyzing Data with an Outlier that does NOT follow the trend */ title1 'Original Data Plus an Outlier that Does NOT Follow Trend'; data casuser.offtrend; set casuser.curve casuser.outlier2; ods graphics / aatrpriority=none; proc sgplot data=casuser.offtrend; styleattrs datasymbols=(CircleFilled SquareFilled); scatter X=x Y=y / group=group; run; ods graphics / attrpriority=color; proc hpsplit data=casuser.offtrend minleafsize=1 maxdepth=2; id group x; model y=x; prune off; output out=hpsplout; run; proc print data=hpsplout; run; proc sgplot data=hpsplout; scatter X=x Y=y / group=group; scatter X=x Y=P_y / markerattrs=(symbol=asterisk color=red); run;         The model also has a leaf node with only one observation. That one leaf node contains the observation in the middle of the X space. So this model seems very similar to the original data with no outlier. In fact, the presence of the outlier in this situation does not affect the model much at all either.   One last situation to consider is the outlier 3 case, where the outlier is in the middle of the X space. As before, here is the code for this situation and some selected output.   /* Creating and Analyzing Data with an Outlier that is in the Middle of the X Space */ title1 'Original Data Plus an Outlier in Middle of X Space'; data casuser.midx; set casuser.curve casuser.outlier3; ods graphics / attrpriority=none; proc sgplot data=casuser.midx; styleattrs datasymbols=(CircleFilled TriangleFilled); scatter X=x Y=y / group=group; run; ods graphics / attrpriority=color; proc hpsplit data=casuser.midx minleafsize=1 maxdepth=2; id group x; model y=x; prune off; output out=hpsplout; run; proc print data=hpsplout; run; proc sgplot data=hpsplout; scatter X=x Y=y / group=group; scatter X=x Y=P_y / markerattrs=(symbol=asterisk color=red); run;         For this situation we have a model with two different leaf nodes that each contain only a single observation. Looking at the plot, we see that the observation in the middle of the X space that follows the trend is on its own in a leaf node. The other observation in its own leaf node is the outlier. The end result of this model is also quite similar to the model with all of the original data. The outlier is “quarantined” in its own node which will limit its influence on the model.   So ultimately, a tree model is quite robust to outliers. The outliers will often be included in a node of their own which limits the influence on the remainder of the model. However, be aware that if one changes the minimum leaf size away from one, outliers will be more likely to have an influence on the model as the outlier will be forced to be combined with non-outlying observations.     Find more articles from SAS Global Enablement and Learning here. ... View more

  Mixture experiments are certainly unique. The constraint that the sum of the ingredients must add to 1 causes multicollinearity that will always be present. That multicollinearity will necessarily lower the power of any experimental design. My previous post showed how to create a mixture design but did not discuss the power. This post will explain how to determine the power of a mixture experiment as well as why the power values may be extremely low.   Let’s start by looking at an example that we can use to understand power for a mixture design. For this example we have three components (X1, X2, and X3) with their associated ranges: X1: 0 – 0.10, X2: 0 – 1, X3: 0 – 1. Using the approach outlined in a previous post, an experiment can be created that might result in a 12-run design (in sorted order):   Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   The twelve points cover the design space quite well and contains some replicated points.     Although the design looks good, there is no indication of power. To calculate power, an estimate of the error is needed. One way to obtain this estimate is to use subject matter knowledge of what the expected responses might be for each design point. With the expected responses, a preliminary analysis can be completed to obtain the error estimate. Since this is a hypothetical example, I will use a formula to generate the expected responses.   YSim = 6*X1 + X2 + X3 + 8*X1*X2 + 8*X1*X3 + Error   This model is a quadratic model as indicated by the cross-product terms. This fact will be important a bit later. The error term is a bit of random error from a normal distribution. Applying the formula yields these realized values:     Using these values, we can fit a model to estimate the error. Recall in an earlier post that we can use a noint option to fit the Scheffe mixture model.   data dsgn; input run x1 x2 x3 ysim; datalines; 1 0.0 1.0 0.0 1.16 2 0.0 0.0 1.0 1.08 3 0.1 0.9 0.0 2.30 4 0.0 0.5 0.5 1.06 5 0.1 0.0 0.9 2.13 6 0.05 0.95 0.0 1.75 7 0.1 0.9 0.0 2.17 8 0.05 0.0 0.95 1.69 9 0.0 0.0 1.0 0.95 10 0.05 0.475 0.475 1.63 11 0.0 1.0 0.0 1.08 12 0.1 0.0 0.9 2.21 ; /* Fitting Mixture Model to Assess Std. Dev. */ title "Regression Analysis to Assess Error"; proc glm data=dsgn; model ysim = x1 x2 x3 x1*x2 x1*x3 x2*x3/noint; run;   The results of the analysis show an error standard deviation (labeled as Root MSE) of 0.066187.     Now we can turn our attention to the power of this design. The PROC GLMPOWER procedure will provide the power of the design.   /* Full Mixture Model */ title "Full Mixture Model"; proc glmpower data=dsgn; model ysim=x1 x2 x3 x1*x2 x1*x3 x2*x3; power alpha = 0.05 ntotal = 12 stddev=0.066 power = .; run;   PROC GLMPOWER requires the model statement where we get to specify our full model. We need to add the POWER statement with options (for this example) for alpha, number of observations, standard deviation, and power. I specified an alpha of 0.05 for the typical 95% confidence level. The design had 12 runs, so ntotal is 12. The standard deviation is 0.066 as we noted earlier, so that is the value that is entered on the STDDEV option. Finally, the power is set equal to a missing value which tells the procedure to solve for power. Unfortunately PROC GLMPOWER does not allow a noint option to truly fit the Scheffe model. This will become apparent on the results.     Because of the inherent singularity in the design matrix, the power for the main effects, X1, X2, and X3, cannot be calculated. The power for all the quadratic/cross-product terms is fairly low. We will explore why this is a bit later. Rather, I would like to get a power estimate for the main effects. This is where the prior post on analyzing a mixture experiment comes into play again. We can estimate a Scheffe mixture model by removing one of the main effects. I will start by removing the X3 term. All other aspects of the PROC GLMPOWER procedure will remain the same.   /* Mixture Model with X3 removed */ title "Mixture Model with X3 Effect Removed"; proc glmpower data=dsgn; model ysim = x1 x2 x1*x2 x1*x3 x2*x3; power alpha = 0.05 ntotal=12 stddev=0.066 power=.; run;   The results will now display the power for the X1 and X2 main effects along with the power estimates of the cross-product terms as before.     We can see that the power for X1 is very low, only 0.095. Why is this so low? The power for X2 is 0.27, which is quite a bit better, but still very low. And how can we determine the power for the X3 term? The X3 power can be determined by repeating the same approach, but removing the X2 term (or X1) and putting X3 back into the model.   /* Mixture Model with X2 removed */ title "Mixture Model with X2 Effect Removed"; proc glmpower data=dsgn; model ysim = x1 x3 x1*x2 x1*x3 x2*x3; power alpha = 0.05 ntotal=12 stddev=0.066 power=.; run;   The results are the same as before, but with the X3 power displayed, which is 0.27, instead of the power for X2.     Now we can finally turn our attention to why the power is so low.Since we know what this model is, we can explore the low power. Certainly some of the low power is due to the multicollinearity found in the data, but there is more. The narrow range for the X1 component also causes some difficulties. Look at the plot of the YSim response and X1 with an overlaid linear regression line and the actual quadratic relationship that was used to generate the response.     The curves are certainly very different. However, X1 only ranges from 0 to 0.1. If we zoom in on just that portion of the X1 range we see a very different picture.     Over the range where the data will be collected the difference between the linear fit and the quadratic fit is almost indistinguishable. The model is not really able to distinguish between the line and a quadratic fit due to the very narrow range. For this reason it is important to have the ranges for your components be as wide as possible in order to magnify the differences between the different possible models.     Find more articles from SAS Global Enablement and Learning here. ... View more

  Outliers have always been an issue for data analysis. Much effort has been spent on identifying outliers and remedying their effects on data analysis. But with the advent of machine learning techniques, do we still need to worry about outliers? The purpose of this post is to discuss a few issues regarding outliers and their effects on neural network models.   Let’s start with a small dataset and a known relationship. Doing this will ensure that we can create a true outlier. The following code will create a 19 run dataset with a random x variable for an input that will range uniformly from 0 to 2, an error term from a normal distribution with mean 0 and a standard deviation of 0.02, and a y variable that is based on the function y=0.2 + 0.6*x – 0.86*x**2 + 0.28*x**3 + error.   data curve; call streaminit(8675309); do i=1 to 19; x=rand('uniform')*2; error=rand('normal')*0.02; y=0.2 + 0.6*x - 0.86*x**2 + 0.28*x**3 + error; output; end; run; title1 'Plot of Original Data'; proc sgplot data=curve; scatter X=x Y=y; run;   Because the error term is small, a plot of the data shows the 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.   This data will be modified with a single outlier. But not all outliers are the same. So we will add that single outlier in three different ways. The first outlier will be at an x value of 2.1, which is higher than any other x value, but the y value will follow the same formula as the original data. Outlier number 2 will have an x value of 2.1, but the y value will NOT follow the rest of the data trend. The y value will be much lower than the trend. And finally, outlier number 3 will have an x value of 1.0, which is right in the middle of the x range, but the y value will be much larger than the formula would provide. This code will generate the three outliers which will each be joined to the original 19 points.   data outlier1; /* Outlier that follows the same trend as the data */ x=2.1; error=rand('uniform')*0.025; y=0.2 + 0.6*x - 0.86*x**2 + 0.28*x**3 + error; output; run; data outlier2; /* Outlier that does NOT follow the trend of the data */ x=2.1; error=rand('uniform')*0.025; y=0.1 + error; run; data outlier3; /* Outlier in the middle of the x-range and does not follow data trend */ x=1; y=0.45; run;   Here are the plots showing these outliers with the original data.      How would a neural network handle building a model with each of these scenarios? And would the fact that these are small datasets cause a problem? Let’s try fitting a neural network model to each of these situations to see if the outlier causes the algorithm any problems. This is not a large problem, so I will use PROC NEURAL with 5 hidden units in the hidden layer and many of the default options.   title1 'Neural Network when Outlier Follows the Trend'; proc dmdb data=followtrend dmdbcat=catlog; var x y; run; proc neural data=followtrend dmdbcat=catlog; input x / level=interval id=inputs; hidden 5 / id=HU1; target y / level=interval id=target; connect inputs HU1; connect HU1 target; train maxiter=100 estiter=1 outest=weights outfit=stats out=Predictions; quit; run;   I won’t bother looking at the fitting results. Instead, I will just focus on a graph that plots x versus y overlaid with the neural network predictions (in red with filled dots).     As expected, this neural network seems to fit the data well, and can predict the outlier quite well. So in this situation, the outlier does not seem to have any impact on the fitting results of this neural network.   Let’s move on to the next situation where the x value is 2.1, but the y value does not follow the rest of the data trend. As before, we will not worry about the actual fitting results and just focus on the graph of the results.     Here we see that the neural network model is modified by the outlier. Neural networks are very flexible, and therefore, may be more sensitive to outliers. They may be reacting to noise. For this situation, remember that 5% of the data is made up by this outlier. Later we will see what happens when we have a larger dataset for the neural network.   Let’s visit our last scenario now with the outlier in the middle of the x range.     Even with the outlier being “surrounded” by similar x value datapoints, the neural network model is significantly affected by the outlier. From these scenarios, it seems as though a neural network model is sensitive to outliers. But as mentioned earlier, perhaps this is really due to such a small dataset. Machine learning models typically work better with more data. Let’s try the same situation, but the original data will have 1999 observations rather than 19. This means that the outlier will only be 0.05% of the total dataset. Following the same approach, here are the three neural network models for these three scenarios.     As expected, when the outlier follows the trend, the neural network follows the signal in the data quite well.     When the outlier does not follow the trend, the neural network model still looks pretty good. It looks like the outlier does not really impact the results at all, but if we zoom in on the area where this outlier is, the story is slightly different.     When you zoom in, you can see that the neural network model is actually still biased because of this one outlier (which is in the lower right on the graph). There are more data points above the neural network prediction line than below it as the x value increases. The results may not be too bad, but even with only one outlier out of the 2000 observations, the model is affected.   You may not be surprised by this result. After all, the x value of 2.1 will have pretty high leverage since it is a larger x value than any other observation. Let’s look at the third scenario to see what happens.     Finally, this looks pretty good. The outlier is at the very top of the graph, but none of the neural network predictions seem to be affected by that point. Let’s zoom in around 0.95 to 1.05 to just make sure there is not a hidden problem.     Even when zoomed in, there does not seem to be a problem with the neural network model.   To summarize, outliers can be a problem with a neural network model. This is especially true if the dataset is small. Even with a larger dataset, if the outlier is an outlier in the x space AND the point does not follow the trend of the rest of the data, the neural network model will still likely be affected. If there are other observations in the input space near that outlier, the impact on the results will likely be minimal. Ultimately, even with a neural network model, care should be taken to identify and possibly mitigate outliers. Some possible actions are to examine the data or using a more robust error function for the neural network. In a future post I will look into the impact outliers have on a decision tree.     Find more articles from SAS Global Enablement and Learning here. ... View more

  In my previous post I discussed how to properly analyze a mixture experiment. But how do you choose which mixture combinations to run to properly estimate the model you want to fit? This is where the design of experiments is helpful. In this post I will discuss an approach to use when creating a designed experiment for a mixture situation.   So what is a mixture situation? A mixture situation is when the “independent variables” or factors have a constraint that they add to 100%. For example, suppose I want to make a mixed drink, a Harvey wallbanger. This drink is made up of three ingredients: orange juice, vodka, and a liqueur. When making this drink, I could make the drink in a small glass or in a large barrel. In other words, the amount of each ingredient does not matter. The only thing that matters is the proportion of each ingredient. So as one proportion goes up, another one (or more than one of the other ingredients) must go down. This constraint of adding to 100% makes the situation unique and causes some challenges in creating the design and in the analysis.   There are several different types of mixture designs: simplex-centroid, simplex-lattice, and extreme vertices are considered the “classic” mixture designs. However, these designs typically require an unconstrained design space or require too many runs to be of practical use. Therefore, optimal designs are typically used which is what will be discussed here.   Optimal designs find the best combination of runs according to a given criterion. There are many choices of criteria, such as A-optimality, D-optimality, U-optimality, and S-optimality. Regardless of the criteria chosen, you must supply the model form that you wish to use for the criterion to be evaluated and the number of runs that you desire for the experiment. For simplicity, I will only focus on D-optimal designs in this post, but you could use whichever optimality criterion you wish.   To get started, PROC OPTEX requires a list of candidate runs. These are possible runs to use in the final design. The algorithm will search through these candidate runs to find the best subset combination of runs to estimate the desired model. The approach is called the row-exchange algorithm since it will exchange one row of the candidate set for a row that is currently in the design. The algorithm continues in this fashion until the best set of runs is found.   As an example, suppose I wish to design an experiment for three mixture variables, A, B, and C. Component A can only be in the range 0.05 to 0.75, Component B has a range of 0.1 to 0.8, and Component C is in the range 0.15 to 0.6. To get started on creating the design, we need to form the candidate set. Typically, a fine grid of points is created that matches the criteria of the design space. This will provide a rich candidate set for PROC OPTEX to find a good design. You can make the grid as fine as you would like, but realize that the finer the grid is, the larger the candidate set, and the longer it may take for PROC OPTEX to run.   data candidates; do a=0.05 to 0.75 by 0.05; do b=0.1 to 0.8 by 0.05; c=1-a-b; if (0.15 <= c <= 0.6) then output; end; end; run;   This code will create a good candidate set. Notice that by forcing C=1-A-B, we are ensured of only getting valid mixtures. However, we also need to add the range of component C before outputting the run to ensure that we meet those constraints as well. If our constraints are well-defined, this step may not be necessary, but this is just a quick and easy way to ensure our constraints are consistent. For this situation there are 97 runs in our candidate set.   With the candidate set created we can now turn our attention to PROC OPTEX. In order to create a design, we need to specify the model that we wish to estimate. This is where knowledge of how to fit a mixture model is needed. For this situation, suppose that we wish to fit the special cubic model, meaning that our model should be Y=b1*A+b2*B+b3*C+b12*A*B+b13*A*+b23*B*C+123*A*B*C (let's call this the NOINT model). As discussed in my last post, this model is not really a no-intercept model, and another way to fit the model is to remove one of the main effects. So this model: Y=b0 + b1*A+b2*B+b12*A*B+b13*A*C+b23*B*C+b123*A*B*C (let's call this one the NOME model for no main effect) which has an intercept will be equivalent to the Scheffe mixture model. The intercept will be the estimate for component C.   proc optex data=candidates seed=27513 coding=none; model a b a*b a*c b*c a*b*c; generate n=12 criterion=D; output out=dsgn; run;   Two options were used on the PROC OPTEX line. The seed makes the design selection reproducible. Using this seed and running the code again will provide the same design. The coding=none option tells SAS to create the design without using any coding on the model effects. In many situations you would want to code the effects to make them orthogonal and on the same scale. However, for a mixture situation this coding is not typically used. You could still perform the coding, but the results are not going to be the same as what you would obtain when fitting the model.   The PROC OPTEX statements include a MODEL line which is fitting the intercept model that was discussed earlier. The GENERATE statement specifies that we are looking for a design with only 12 runs that uses the D-optimality criterion. Finally, the output line specifies the dataset that will contain the final design.   Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   PROC OPTEX shows the ranges for each of the components as well as the ten designs that were considered by the procedure. Considering 10 designs is the default for PROC OPTEX, but you could specify more (or fewer) by using the ITER=n options on the GENERATE statement. The output from PROC OPTEX can be seen with a PROC PRINT statement.     Or graphically:     Looking at this design you can see that it seems to cover the design space well. Although 9 points appear on the graph, by looking at the printed design you can see that observations 3 and 4 are the same conditions as are observations 8 and 9 and observations 11 and 12.   This approach to creating the design makes sense, but there are other ways to consider a mixture model. What if we want to use the NOINT model mentioned above? All of the code would look the same except for the MODEL statement in PROC OPTEX.   proc optex data=candidates seed=27513 coding=none; model a b c a*b a*c b*c a*b*c / noint; generate n=12 criterion=D; output out=dsgn2; run;   Alternatively you could also specify the model as   MODEL A|B|C / NOINT;   Running this code, printing out the design, and displaying it graphically will show that this model will give you the exact same design as the NOME model.   Regardless of how you wish to specify your mixture model, PROC OPTEX will find an appropriate mixture design given your optimization criteria and the number of runs.     Find more articles from SAS Global Enablement and Learning here. ... View more

In my previous post I discussed no-intercept regression models, when to use them, and how to interpret the results.  The purpose of this post is to build on those concepts but present them in the unique situation of formulation data that changes the interpretation of the models.   A formulation is where the predictor variables, called ingredients or components, sum to a constant value, usually 1 or 100%. An example would be a situation where three different fruit juices need to be blended to create a new better tasting juice. You could create a glass of the blend or a tank truck full of the blend and it will not matter. Because of this, we measure the proportion of the ingredients in the blend rather than their amounts. These situations are common in chemistry and pharmaceuticals to name just two industries.   Analyzing formulation data has some inherent challenges. Because of the constraint that the ingredients add to one, a typical regression model cannot be fit. The constraint or multicollinearity makes it impossible to find the parameter estimates. As an example, this program creates a small example of some formulation trials and fits a regression model.   data mix; input x1 x2 x3 y; datalines; 1 0 0 38.17 1 0 0 37.00 0 1 0 41.10 0 1 0 41.65 0 0 1 71.27 0 0 1 69.50 ; run; title1 'Regression with Intercept'; proc reg data=mix plots=(fit(nolimits)); model y=x1 x2 x3; run;   The results from this regression show the issues with formulation data. A note in the output states that the model is not full rank. To provide some analysis, SAS will set a parameter to 0. The parameter that is set to zero will always be the last main effect in the model. This also leads to other parameter estimates being biased. The output even includes a line indicating the multicollinearity that is found.    Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   So how should the formulation data be analyzed? The models are adjusted using the constraint. The Scheffe canonical form is the most common adjustment. The original regression model for this example is Y=B0+B1*X1+B2*X2+B3*X3. But since X1+X2+X3=1, then we can rewrite this as Y=B0*(X1+X2+X3)+B1*X1+B2*X2+B3*X3. Simplifying and regrouping the terms leads to Y=(B0+B1)*X1+(B0+B2)*X2+(B0+B3)*X3. We can then rewrite the terms in parentheses as B1', B2', and B3' respectively which yields Y=B1'*X1+B2'*X2+B3'*X3. This reparameterization of the model allows the model to be fit despite the constraint. So although the model looks like a no-intercept model, there is an intercept. The intercept is included with the main effect parameter estimates.   With this unusual model one might think that it would be difficult to fit this Scheffe canonical mixture model. But, in fact, it is not difficult. You have3 three choices. The first approach would be to use the NOINT option of PROC REG.   title1 'Regression with NOINT Option'; proc reg data-mix plots=(fit(nolimits)); model y=x1 x2 x3/noint; run;   This code would yield these results:     With this output, you will notice that the model is considered a no-intercept model, as expected. This means that the Analysis of Variance table will not be correct for this revised model. You can see this by looking at the degrees of freedom for the model. There are three parameters in the model, so the correct degrees of freedom should be 3-1=2 since there is still a mean to estimate. However, the parameter estimates are correct. This model would be Y=37.585*X1+41.375*X2+70.385*X3. The tests on those parameter estimates are using the null hypothesis of the parameter estimate=0. Although the test is performed correctly, the test is of limited usefulness since the parameter estimate contains the intercept. We would not expect the parameter estimate to be equal to 0 in this situation.   The second possibility is to use the RESTRICT INTERCEPT=0 option. The code and the results are here.   title1 'Regression with RESTRICT Option'; proc reg data=mix plots=(fit(nolimits)); model y=x1 x2 x3; restrict intercept=0; run;     These results are essentially the same as what was seen when using the NOINT option.   Finally, the last option may not make much sense at first, but it does work. Instead of fitting a no-intercept model, fit a model with an intercept, but leave off one of the main effects. The intercept line of the parameter estimates table will be the parameter estimate for the removed term. This approach will have the added benefit of providing the proper ANOVA table.   title1 ' Regression with Intercept, but One Term Removed'; proc reg data=mix plots=(fit(nolimits)); model y=x1 x2; run;     Starting with the Parameter Estimates table, you will see that the intercept has the same estimate as X3 from the other analyses and the X1 and X2 estimates are the same also. We also have the added benefit of the ANOVA table now being correct since there are only 2 degrees of freedom for the model.   Formulation models are different from other regression models. They appear to be a no-intercept model, but in reality, there is an intercept that is not expressed directly due to the constraint the components must add to 1. Although NOINT and RESTRICT options could be used, the most informative results can be obtained by simply removing one of the effects from the model and leaving the intercept as the placeholder for the missing effect.     Find more articles from SAS Global Enablement and Learning here. ... View more

  Sometimes when using a linear regression model, you might wish to remove the intercept from the model. This situation commonly occurs when you theoretically know that the model should go through the origin. Another way of stating that the model goes through the origin is that the intercept of the model should be zero. Forcing the intercept to be zero has some consequences, which is what this blog is about.   Removing the intercept from a model using SAS code is straight-forward and can be done in two common ways with PROC REG: using either the RESTRICT statement or the NOINT option.   To see this in action, here is some sample code that creates a small dataset. We are measuring the height and weight of items of various sizes. We will start by fitting the typical regression model with an intercept and examining those results.   data demo; input height weight; datalines; 24 170 31 197 29 194 42 242 58 305 53 289 ; title1 ‘Regression WITH the Intercept’; proc reg data=demo plots=(fit(nolimits)); model weight=height; run;   Looking at the results of the model with an intercept, notice that the Analysis of Variance table has one degree of freedom for the model, 4 degrees of freedom for the error, and 5 degrees of freedom for the Corrected Total. This makes sense because we have two parameters estimated in the model (an intercept and a slope). These parameters can be seen in the parameter estimates table.   Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   The statistical test that is performed by the Analysis of Variance table has a null hypothesis of the model equal to the mean of the response variable. Because the Pr>F is small, we would decide that the model is significant. The model is different than just using the overall mean of the response.     The plot of the model shows a pretty good fit to the data.   Now let’s fit a no intercept model. This model will be fit in two ways. First will be using the NOINT option of PROC REG. The second method will be to use the RESTRICT statement in PROC REG will be shown later. The code using NOINT would look like this.   title1 'Regression Using the NOINT Option'; proc reg data=demo plots=(fit(nolimits)); model weight=height / noint; run;   Let’s explore the no intercept model created with the NOINT option first.     There is a note in the output that the model R-Square is redefined because this is a no-intercept model. R-Square is usually calculated as SS(Model) divided by SS(Total) where the SS(Total) is the corrected total sum of squares. But the ANOVA table shows an uncorrected total sum of squares is being used. So the R-Square is redefined to be SS(Model) divided by SS(Uncorrected Total). As a result, the interpretation of R-Square is no longer the proportion of variance explained and does not have a clear interpretation.   One must also be careful in interpreting the residuals from a no-intercept model. Because the residuals may not sum to zero, the residuals may not be evenly distributed around the zero line. Residual plots may look different than expected.   Finally, the plot of the model shows that this model does not fit the data very well. Why not? Clearly, the model should go through the origin as you cannot have weight without a height.     There may be a few good reasons that this model does not fit the data very well. The most likely cause is that the model may not actually be linear. The lowest height value that we have is 24. We have no data to tell us what the relationship between height and weight is from 0 to 24. We probably should not assume the linear relationship seen with this data holds all the way to the origin. That is an extrapolation. Another possible cause is that the variance of the error terms may not be constant. Perhaps the error term grows as you get closer to the origin. That could lead to a large intercept that would still be within the bounds of the error term.   Now that we have fully explored the no-intercept regression model using the NOINT option, what do we see if we use the RESTRICT statement? Let’s submit this code.   title1 "Regression Using the RESTRICT Option"; proc reg data=demo plots=(fit(nolimits)); model weight=height; restrict intercept=0; run;   The results are almost the same as the NOINT option. However, the RESTRICT statement adds one additional line, the RESTRICT line, to the Parameter Estimates table. This line is testing the sensitivity of the sum of squared errors to the restriction. If the restriction estimate is small, then the restricted estimates are the same as the unrestricted estimates. The p-value that is provided indicates that the slope for height has a statistically significant difference from the unrestricted estimate. This might lead one to question if using the no-intercept model is really appropriate in this situation.   With all the cautions, why would you even fit a no-intercept model? There are still situations where you may need this capability. Perhaps you want to build your own custom design matrix that does not include the intercept column. Or maybe you are trying to compare two different measurement systems. You consider measurement system one as the independent variable and measurement system two as the response. Fitting a regression line with a slope of 1 and an intercept of 0 would provide an appropriate comparison of the two systems. One last possibility is if you need to fit a Scheffe mixture model which is not exactly a no-intercept model, but we will see more about this situation in my next blog post. These are just a few situations where you might want to fit a no-intercept regression model.   Ultimately, one should be very careful when using a no-intercept model as the typical regression interpretations no longer apply. In general, the safest course of action is to not use the no-intercept model. Interpretations with an intercept model will be more understandable. Even if the model does go through the origin, an intercept model should only have a small intercept due to sampling error and should be of no consequence. There is no drawback to using the intercept model. ... View more