Rick_SAS Tracker

Re: What is collinearity and why does it matter?

Rick_SAS — Mon, 27 Jul 2026 12:58:42 GMT

For more details and examples of collinearity in regression data, and how to use the COLLIN option and diagnostic statistics in SAS, see

Re: Model-Based Clustering (Part 1): Exploring Its Significance

Rick_SAS — Thu, 23 Jul 2026 13:20:51 GMT

To learn more about Model-Based Clustering and the algorithm behind it (the E-M algorithm), see Fit a multivariate Gaussian mixture model by using the expectation-maximization (EM) algorithm - The DO Loop The article uses SAS 9.4. I do not claim it is a replacement for PROC MBC, but it shows the main ideas.

Re: Mean and Standard Deviation for Large Amounts of Data

Rick_SAS — Sat, 20 Jun 2026 15:28:53 GMT

The Welford method is more useful for simulation studies than for data analysis.

The MEANS and SUMMARY procedures are the usual methods in SAS 9.4 for descriptive statistics. They excel at producing statistics across BY groups and CLASS variables (including n-way crosstabulations), which is a common situation in demographics. They also compute standard errors and confidence intervals for the population quantities.

The SUMMARY and HPSUMMARY procedures use multithreading to divide the computations across threads. You can get almost perfect parallelization since both procedures are essentially summations.

data Have;
array mu[100] _temporary_ (1:100);
array x[100];
do i=1 to 1E6;
  do j=1 to dim(x);
     x[j] = rand("Normal", mu[j], 20);
     if (rand("Bern", 0.01)) then x[j]=.;  /* add some missing values */
  end;
  output;
end;
run;

options fullstimer;
proc hpsummary data=Have mean std;
   var x:;
   output out=Summary;
run;
options nofullstimer;

proc print data=Summary;
   where _STAT_ in ("N" "MEAN" "STD");
   format _numeric_ best6.;
run;

Re: How to use ODS OUTPUT with PROC DATASETS

Rick_SAS — Thu, 04 Jun 2026 21:05:23 GMT

I agree with Quentin and Tom. An ODS SELECT/OUTPUT statement must go before the statement that generates the table (or graph). In this case, it is the PROC DATASETS statement itself that creates the output. Run this example WITHOUT any RUN or QUIT statement:

data class1 class2;
  set sashelp.class(obs=3);
run;
ods trace on;
proc datasets lib=work memtype=data;

Even though there is no RUN or QUIT statement, you can see from the ODS trace statement that the output is created.

I was not aware of this behavior. In my 2024 article about interactive procedures (Interactive procedures in SAS - The DO Loop), I stated the following information from the PROC DATASETS documentation:
The interactive nature of PROC DATASETS is a different from other SAS procedures. The documentation states that it supports RUN-group processing, but "the DATASETS procedure supports four types of RUN groups.... Some statements in PROC DATASETS act as implied RUN statements because they cause the RUN group preceding them to execute."

I have boldfaced the statement that applies in this case. When I wrote the article, I don't think I fully appreciated what that boldfaced statement was telling us.

Re: CV for ridge regression

Rick_SAS — Mon, 11 May 2026 09:45:28 GMT

I haven't done this myself, but there is an example in the GLMSELECT documentation that appears to use the ELASTICNET method to do what you want. Here are some links that I hope will be useful:

SAS Help Center: Elastic Net Selection (ELASTICNET)
SAS Help Center: MODEL Statement See the SELECTION=ELASTICNET method and the various "L2" options. The default search option for the optimal ridge parameter are L2LOW=0, L2HIGH=1, L2STEPS=50, and L2SEARCH=GRID. These look like good defaults to me.
SAS Help Center: Elastic Net and External Cross Validation

Re: MLEs of GE

Rick_SAS — Mon, 13 Apr 2026 09:44:37 GMT

Is there anything else we can help you with? If not, please close this thread.

Re: MLEs of GE

Rick_SAS — Wed, 08 Apr 2026 14:06:32 GMT

> I still get the same error:

ERROR: NEWRAP Optimization cannot be completed.

WARNING: Optimization routine cannot improve the function value.

A more accurate statement is, "For about 10% of the random samples of size N=25, the optimization does not converge."

This is typical behavior for three-parameter models, especially those with a threshold parameter. In some cases, the log-likelihood does not have a maximum on the feasible domain. The parameter estimates (especially the xi parameter for GEV) are very sensitive to samples whose shapes deviate from the distribution. For example, if the random sample generates too many or too few observations from the tail of the distribution, or few observations close to the threshold, then the MLE might not converge.

You can find many papers and textbooks that discuss these issues. For example, see Rosian, Na, and Gabda (2020), whose abstract begins with the sentence, "The standard method of the maximum likelihood has poor performance in GEV parameter
estimates for small sample data."

I guess the important question to ask is, "what are you trying to accomplish or learn by running this simulation study?"

Re: SAS Code Errors

Rick_SAS — Thu, 19 Mar 2026 09:56:08 GMT

First, I made a mistake. The statement that subsets the good indices should be
Est = Est[good_idx, ]; /* Note the comma! */

To answer your questions:
> 1. the mean of estimates of betas are far away from the true values. Is it acceptable?

I don't know what "acceptable" means. What is your goal? Your estimates will be far away from the true mean regardless of which optimization method you use.

If your goal is to show that the regression estimates for this model are highly variable and should not be trusted for samples of size n=25, then, yes, I think this simulation demonstrates that fact very well. If you plot your parameter estimates, there are many extreme values, which will, of course, strongly influence the Monte Carlo mean and standard deviation. A more robust Monte Carlo estimate would be the median of the columns. You should always graph the parameter estimates to visualize their variability.

colMedian=t(median(Est));
print colMedian;

title "Graph of parameter estimates from NLPNRA";
call scatter(Est[,1], Est[,2]);

The graph might cause you to ask whether the computations in the simulation are wrong. You can answer that question by performing the same computation for larger sample sizes. Many of your "problems" with nonconvergence and highly variable estimates will go away if you increase the sample size, to, say, n=100.

2. Yes, the Monte Carlo estimate is based on the 997 estimates that converged

3. Yes, you can increase the value of the 'replicate' variable to obtain additional Monte Carlo samples. But the math doesn't care if you use 997 or 1000. The math says that the mean of the estimates (however many you have) is the Monte Carlo estimate.

In summary, your "problems" are caused because the sample size in the simulation is so small. You can "fix" your problems by using a larger sample size.

Re: SAS Code Errors

Rick_SAS — Thu, 19 Mar 2026 09:56:46 GMT

> I got the following error 3 times in the log window when n=25:

ERROR: NEWRAP Optimization cannot be completed.

ERROR: NEWRAP needs more than 200 iterations or 500 function calls.

How can I fix this?

This error message appears for 3 out of 1000 simulation trials. Errors like this appear in simulation studies when a random sample contains extreme values or is otherwise not representative. For those samples, the MLE step does not converge. This does not indicate a problem in your code, but rather a hazard when you are simulating small data sets (in your example, n=25). For details and a discussion, see Monitor convergence during simulation studies in SAS - The DO Loop

There are two ways to address your concern:

Unless you are specifically interested in Newton-Raphson optimization, you can switch optimization schemes. For example, the quasi-Newton optimization in NLPQN might be faster and more robust. Since quasi-Newton methods do not require Hessian computations, you can replace your call to NLPNRA with
call nlpqn(rc, res, "SUMLOGLHF", initial_param,opt) grd="driv"; /*Use quasi-Newton algorithm that minimizes negative log likelihood function */
This is what I would do.
If, for some reason, Newton iteration is important, you can keep track of the return code (rc) and remove the results from samples for which the MLE does not converge. This is the suggestion in the blog post, as applied to your IML program. This requires allocating a vector for the return codes, then using the LOC function to keep only the results from MLE that converged:

rc_status = j(replicate,1);
do k=1 to replicate;
...
   call nlpnra(rc, res, "SUMLOGLHF", initial_param, opt) grd="driv" hes="HessLogLik";
   Est[k,]=res;
   rc_status[k] = rc;
end;

/* find trials that converged */
good_idx = loc(rc_status > 0);
Est = Est[good_idx, ];
nSim = nrow(Est);
colMean=t(mean(Est));
colStanddev=t(std(Est));
print nSim colMean colStanddev;

Re: SAS runs slower with each macro call simulated power

Rick_SAS — Tue, 17 Mar 2026 19:09:47 GMT

Chapter 6 "Strategies for Efficient and Effective Simulation" of Wicklin (2013) has a lot of techniques for efficient simulation. Another technique that applies to you is the value of running a small preliminary study (pp. 94-95). In your study, you could initially run 1000 (instead of 10,000) simulations while you develop and debug the program, and as you search for the optimal parameter values. That will enable you to conjecture that the optimal sample size is approximately N=470. You can then increase the number of simulations to 10,000 and re-run the simulations for a limited range of sample sizes such as [465, 475].

To quote Wicklin, p. 95:
In summary, when you begin a large simulation study, resist the urge to immediately write and
run the entire simulation. Start small. Debug and optimize the simple cases. Run a small-scale
simulation. These techniques provide important information that you can incorporate into the final
study.

Re: SAS runs slower with each macro call simulated power

Rick_SAS — Tue, 17 Mar 2026 18:39:36 GMT

> I'm using the old Display Manager. Only a two-way proc freq and the four observation dataset generated by it is going to the output window.

If you expand the control widgets for the PROC FREQ in the Results window, you will see that you are generating 2 x 10,000 little control widgets every time you run the macro.

Put

ods noresults;

at the top of your macro and

ods results;

at the end.

Re: SAS runs slower with each macro call simulated power

Rick_SAS — Tue, 17 Mar 2026 18:20:09 GMT

What client or editor are you using to run this program? SAS EG? SAS Studio? The old Windows Display Manager?

From your description, it sounds like some resource is accumulating, which is why the process slows down. If you are seeing output to the log, an ODS output destination, or even the Results window in the old Display Manager. See Turn off ODS when running simulations in SAS - The DO Loop for a discussion and some tricks.

Re: SAS Code Errors

Rick_SAS — Mon, 16 Mar 2026 18:02:51 GMT

The statement

z = max(-700, min(700, z));/*Prevents exp() overflow/underflow*/

is wrong.

The MIN and the MAX function return a scalar value. So instead of z being a vector, it is a scalar.

The correct way to clip a vector of values into an interval is shown in the article
"Clip values into an interval"

The article contains this little IML helper function:

/* use the elementwise minimum (><) and elementwise maximum (<>) operators to 
   truncate elements of a vector into an interval. See
   https://blogs.sas.com/content/iml/2026/02/04/clip-values.html
*/
start Clip(x, ab);
   a = ab[1]; b = ab[2];
   return( a <> (x >< b) );
finish;

Use that function instead of your nested MAX-MIN statement. For example:

start SUMLOGLHF(beta) global(n,x,y);
z = x*t(beta);
z = Clip(z, {-700 700}); /*Prevents exp() overflow/underflow*/
...

There might be other issues, but fix this issue first.

Re: PCA : correltaion

Rick_SAS — Mon, 16 Mar 2026 09:26:18 GMT

If you need a refresher on the matrix operations behind creating scores plots, see "A classical principal component analysis in SAS/IML" in the article, "Robust PCA in SAS."

Re: proc princomp : labels

Rick_SAS — Tue, 24 Feb 2026 11:06:54 GMT

Ah! Thanks @PaigeMiller for pointing out that the OP wants to plot variable LABELS, not names. I agree that there is no built-in option to make the ODS graphics from PROC PRINTCOMP use labels. The OP will need to use the VALIDVARNAME=ANY option to rename variables as n-literals. Here is an example:

/* What are the labels for the variables? */
proc contents data=sashelp.iris varnum;
   ods select Position;
run;

/* create a copy of the data and use VALIDVARNAME=ANY to rename the 
   variables to their label values */
option validvarname=ANY;
data Iris;
   set sashelp.iris;
   rename 
   PetalLength = 'Petal Length (mm)'n
   PetalWidth  = 'Petal Width (mm)'n
   SepalLength = 'Sepal Length (mm)'n
   SepalWidth  = 'Sepal Width (mm)'n;
run;

/* run PCA on the new data set */
ods graphics on;
proc princomp data=Iris
     plots=pattern(vector circles=25 50 75 100);
var _numeric_;
run;

PatternPlot30.png

I assume this is a one-time task, but if you need to incorporate it into a process, you can write a macro that uses the output of PROC CONTENTS tp automatically create the new variable (with long n-literal names) as copies of the original variables. See 35973 - Display variable labels instead of variable names in procedure results

Re: proc princomp : labels

Rick_SAS — Mon, 23 Feb 2026 21:12:46 GMT

It's hard to tell from the image you posted, but it looks like you want a "Pattern Plot" (somtimes called a "loadings Plot"). Try imitating this code:

ods graphics on;ods graphics on;
proc princomp data=sashelp.iris plots=pattern(vector circles=25 50 75 100);
var _numeric_;
run;

If you need a refresher on how to interpret a pattern plot, see https://blogs.sas.com/content/iml/2019/11/04/interpret-graphs-principal-components.html

Re: predicting outcome values with lsmeans in quadratic regression

Rick_SAS — Sat, 21 Feb 2026 11:23:40 GMT

In addition to what KSharp said in his reply, the observed DM values are not important when scoring the model. The predicted values depend only on the explanatory variable(s). The DM is only used during the fitting processes. It is also used to form residuals, which are useful for assessing the goodness of the fit.

Re: proc sgscatter display issues (for 18 variables)

Rick_SAS — Fri, 20 Feb 2026 16:55:40 GMT

Assuming that you want at least 200 pixels for each graph, a 1000x1000 plot can only hold 5 plots in each dimension. If you are trying to create a matrix that shows the pairwise scatter plot of 18 variables, you are only getting 55 pixels per plot, which is very small.

What are you hoping to see in this display, which is going to have 153 scatter plots? Depending upon your goals it might be better to plot subsets of your data instead of trying to see it all at once. For example, there are 10 numeric variables in the SASHELP.CARS data set. If I want to see how all variables are related to the MPG_CITY variable, I can create a big panel, but sometimes it is better to display multiple panels that each have 5-6 graphs, like this:

ods graphics / height=640px width=480px;
title "First group of vars";
proc sgscatter data=sashelp.cars;
  compare y=(MSRP Invoice EngineSize Cylinders Horsepower)
          x=(mpg_city);
run;

title "Second group of vars";
proc sgscatter data=sashelp.cars;
  compare y=(MPG_Highway Weight Wheelbase Length)
          x=(mpg_city);
run;

Re: predicting outcome values with lsmeans in quadratic regression

Rick_SAS — Thu, 19 Feb 2026 13:45:50 GMT

/* create fake data set that has the same names as the OP's data */
%let Orig_Data = sashelp.class;
data poly_data;
set &Orig_Data;
DM = Weight;
hedonic = Age;
hedonic2 = hedonic**2;
hedonic3 = hedonic**3;
keep DM hedonic hedonic2 hedonic3;
run;

proc glm data=poly_data;
    model DM = hedonic hedonic2 hedonic3/p solution alpha=.05 clparm;
    output out=GLMOut Predicted=Pred;
quit;

Of course, any obs with the same value of hedonic will have the same predicted value. If you want to see only the unique predictions, you can sort by hedonic and then output only the first observation for each level:

proc sort data=GLMOut;
by hedonic;
run;

data GLMOut_NoDups;
set GLMOut;
by hedonic;
if first.hedonic;
run;

proc print; run;

Re: Are Studentized Residuals a Better Measure for Assessing Model Normalcy Than Plain Vanilla Resid

Rick_SAS — Tue, 17 Feb 2026 16:33:14 GMT

> it conforms with our shared expectation that OLS basically requires residuals to “be iid normally distributed.”

Just to clarify, OLS does not depend on iid normal residual for prediction. That assumption is used for inferential statistics such as confidence limits for parameters. So "iid normal" is sometimes considered to be an optional assumption that is useful for some analyses.

For more on this topic, see On the assumptions (and misconceptions) of linear regression - The DO Loop