Rick_SAS Tracker

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;

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

Re: predicting outcome values with lsmeans in quadratic regression

Rick_SAS — Fri, 13 Feb 2026 10:49:10 GMT

If hedonic is a "categorical predictor variable," you should include a CLASS statement? For a classification variable, you can only evaluate the model at levels that are in the model. If you model hedonic as a continuous variable, you can evaluate the model at values that are not in the input data.

From your description, it sounds like you want the predicted value of the model at specific values of the independent variable. As KSharp says, the simplest way to get that information is to "score" (or evaluate) the model at specific values of hedonic by using the STORE statement to save the model and the PLM procedure to evaluate the model at the specified values. Here is an example that uses the Sashelp.Class data:

proc glm data=sashelp.class plots=none;
/* class age; */ /* ??? Treat as categorical?  */
model Weight= age age*age;
store out=GLMModel;
quit;

data AtLevels;
input Age @@;
datalines;
12 12.5 14 16
;

proc plm restore=GLMModel;
  score data=AtLevels out=GLMOut;
run;

proc print; run;

The input data does not contain any subjects with Age=12.5. If you include the CLASS statement, the model cannot form a prediction at Age=12.5. However, if you treat Age as a continuous variable, the model can calculate the predicted value.

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

Rick_SAS — Thu, 12 Feb 2026 14:28:14 GMT

SAS regression procedures provide many diagnostic plots that you can use to assess the fit of the model and to identify potential outliers and influential observations. See An overview of regression diagnostic plots in SAS - The DO Loop

In the article, I discuss two diagnostic plot that use studentized residuals. The most important feature of a studentized residual is that it is standardized. Therefore, you can use rules such as "If the magnitude of the studentized residual exceeds 2, you should examine the observation as a possible outlier."

> are Studentized residual plots best used to find outliers, but should not be relied on to determine whether the specified model has a better "fit" to normal expectations?

My opinion is that the studentized residual plots are best used to identify outliers. I usually use the raw residuals for Q-Q plots and histogram of residuals, but that is primarily because I like to see the residuals on the scale of the response variable, not on some standardized scale. Because the studentized residual for x[i] excludes x[i] from the computation, I have a hard time visualizing studentized residuals, whereas the raw residuals (y[i] - predicted_y[i]) are very easy to visualize and understand.

Re: Burr distribution CDF

Rick_SAS — Thu, 29 Jan 2026 10:36:29 GMT

I am updating this old thread so that others who search for this information can find it. The following articles discuss using the Burr XII distribution in SAS:

Re: I want to write a SAS/IML module that invokes ODS statements at execution time, not compilation

Rick_SAS — Wed, 21 Jan 2026 14:45:49 GMT

It looks like you intend this functionality as part of library that others will use.

If so, I encourage you to avoid using ODS _ALL_ CLOSE, which will cause havoc with your users' open ODS destinations. You should use ODS PDF CLOSE only. If you want to exclude the output from other open ODS destinations, then use ODS EXCLUDE ALL prior to writing the PDF file, followed by ODS EXCLUDE NONE when you are ready to resume normal output. Or, if you know the name of the output table(s), you can use ODS EXCLUDE <nameoftable>.

For more details about why ODS CLOSE ALL is a bad programming practice in a public library, see Five reasons to use ODS EXCLUDE to suppress SAS output - The DO Loop.

Re: I want to write a SAS/IML module that invokes ODS statements at execution time, not compilation

Rick_SAS — Mon, 19 Jan 2026 11:25:19 GMT

I think your confusion is that ODS statements are not IML statements.

The ODS statements are GLOBAL statements in SAS. They are similar to TITLE, FOOTNOTE, LIBNAME, FILENAME, etc.

Global statements are never seen by PROC IML! When SAS parses a program, global statements are sent to other parts of SAS to execute. This is true in every procedure. So, for example, in PROC REG, the TITLE and ODS statements are never seen by PROC REG. Those statements are passed to something called the "supervisor" which handles them appropriately.

You might want to read the article "Calling a global statement inside a loop," which discusses the issue. The article discusses using a global statement inside a loop, but the concept of using a global statement inside a function is similar. IML only compiles and store IML statements. Global statements are stripped out and never seen by IML.

The article shows an example that uses CALL EXECUTE, which you are already using.