About Rick_SAS

Rick_SAS · ‎08-16-2011

I would be careful of using some of the older papers. Although the code/macros should still work, SAS has added a lot of Bayesian features in the last few releases. These will likely make the analysis simpler to program. An overview of survival analysis procs is at http://support.sas.com/documentation/cdl/en/statug/63962/HTML/default/viewer.htm#introsurv_toc.htm The link includes a section titled "Bayesian Survival Analysis with SAS/STAT Procedures."

Rick_SAS · ‎08-15-2011

By default, weights are used to fit the model, but not to compute the ROC. See if the ROCOPTIONS(WEIGHTED) option on the PROC LOGISTIC statement addresses your needs. Rick

Rick_SAS · ‎08-11-2011

Do you have this in a macro loop. Otherwise, I don't see how you could be filling up the log. Maybe the Forum software "ate" some of your code? I am also seeing an empty VAR statement, which doesn't look right. Without the code to guide me, try this: %macro ODSOff(); ods graphics off; ods <dest> exclude all; ods noresults; %mend; %macro ODSOn(); ods results; ods <dest> exclude none; ods graphics on; %mend; Replace <dest> with whatever destingation you are using (sounds like LISTING.) Before the UNIVARIATE code, use %ODSOff, and after the proc call use %ODSOn. (Of course, you also want ODS TRACE OFF.)

Rick_SAS · ‎08-11-2011

Here's how to do it, but I hope you realize that this is not actually predicting anything. What you are doing is simulating Y, given the predicted probabilities. A prediction would use a statement such as PredY = (p>0.5); Anyway, here's what you asked for. /* Simulate random values for p. Pretend the PRED data is the output from LOGISTIC */ data pred; drop j; do j = 1 to 8e6; p = rand("Uniform"); output; end; run; /* use predicted values to simulate SimY ~ Bernoulli(p) */ data Sim; call streaminit(123); set pred; SimY = rand("Bernoulli", p); run;

Rick_SAS · ‎08-10-2011

No offense, but I wish you'd stop saying "SAS can't do this" and "SAS takes days to run." The following code runs in less than two seconds: data b; drop j; call streaminit(123); do j = 1 to 8e6; p = rand("Uniform"); x = rand("Bernoulli", p); output; end; run; proc print data=b(obs=10);run; proc means data=b mean std q1 median q3; var x; run; I need to do some other work now. Good luck with your project.

Rick_SAS · ‎08-10-2011

By the way, you can also sovle this problem in the DATA step. It takes about 2.5 seconds to generate the random values and to write all 8 million obs to a data set: data b; drop p1-p5 i j; call streaminit(123); array p [5] (0.05 0.1 0.5 0.75 0.99); array x [5] x1-x5; do i = 1 to 8e6; do j = 1 to 5; x{j} = rand("Bernoulli", p{j}); end; output; end; run; proc print data=b(obs=10);run; proc means data=b; run;

Rick_SAS · ‎08-10-2011

No loop over observtions is required. Like R, SAS/IML has the ability to fill many random values with a single call. Go to http://blogs.sas.com/iml and browse some of the articles under the "Sampling and Simulation" topic in the right-hand sidebar. For a particularly easy example, see http://blogs.sas.com/iml/index.php?/archives/9-Efficient-Sampling.html. Also see p. 29 and Chapter 13 of my book, Statistical Programming with SAS/IML Software.(Ch. 2 is available as a free download.) For your example, you say you want to sample from Bernoulli (output will be 0 or 1), but in R you are sampling from a binomial distribution (with "number of trials"=size=1) to get this. In SAS, you can sample from Bernoulli directly. Here's how you can do this in SAS/IML. It takes about 0.8 seconds: proc iml; call randseed(123); p = {0.05,0.1,0.5,0.75,0.99}; /* each column is for a separate probability */ x = J(8e6, nrow(p)); /* 8 million obs */ /* Use y to fill with random Bernoulli draws with prob=p */ y = x[,1]; t0 = time(); /* time the computation */ do i = 1 to nrow(p); call randgen(y, "Bernoulli", p); x[,i] = y; /* assign into the i_th column of x */ end; t = time()-t0; print t; /* print elapsed time = 0.8 seconds */

Rick_SAS · ‎08-08-2011

An interesting idea, but I think there is an easier way. For simple models like REG and GLM, you can use the SCORE procedure to score new data. You use the OUTEST= option (reg) or the OUTSTAT= option (GLM and other procs) to create an output data set that holds the parameters estimates that PROC SCORE will use. For MIXED, the model is more complicated. The way to score a mixed model is to use the STORE statement to output the model, and them use the PLM procedure to score the model on new data. The PLM procedure was introduced in SAS 9.2M3 for exactly this purpose. See http://support.sas.com/resources/papers/proceedings10/258-2010.pdf

Rick_SAS · ‎08-08-2011

Maybe I'm confused. I think all you have to do to get the groups is to assign group=1 when the predicted value is greater than 0.5 and group=0 when the predicted value is less than 0.5. I don't see why RANDGEN comes into play or why there would be a loop. Here is some code that generates some fake data and calls logistic to output predicted probabilities. The PROC IML code just assigns 1 or 0 depending on the predicted probabilities. You can use the DATA step to do the same thing. data a(drop = i prob); call streaminit(321); do i = 1 to 1000; x = rand("normal"); prob = exp(x) / (1 + exp(x)); y = rand("Bernoulli", 1-prob); output; end; proc logistic data=a; model y(event='1') = x; output out=out pred=pred; run; proc iml; use out; read all var {pred y}; close out; class = (pred>= 0.5); print (sum(class=y));

Rick_SAS · ‎08-07-2011

One suggestion: solve the normal equations more efficiently by using the SOLVE function instead of using the INV function to form the general inverse. See the Programming Tip on p. 51 of Statistical Programming with SAS/IML Software. You can download Chapter 2 of the book for free at the books Web site. The module that you post is from a previous Discussion Forum post and my blog (although you've deleted my helpful comments and did not cite the source). For details, interested readers can see http://blogs.sas.com/iml/index.php?/archives/66-Creating-a-Matrix-with-All-Combinations-of-Zeros-and-Ones.html . The rows of this matrix (call it m) enumerate the 2##p subsets of the p variables in the model. The first row is the intercept-only model. For the other rows, set Xi = X[,loc(m[i,])] and solve the i_th regression problem with the data matrix Xi. This is the brute force method; there are more efficient ways (such as not storing all those zeros/ones, but rather computing them in the same loop that does the regression). Notice, that 2#20 is more than 1 million, so even if you spend only 0.01 seconds on each regression, we're talking about 10,000 seconds = 175 minutes = 3 hours. If you are getting preliminary results in 40 minutes, I'm impressed. (I guess it depends on the size of your data.)

Rick_SAS · ‎07-29-2011

Assuming that the function has K maxima, I think that a "heat map" is a way to go. For each initial condition, the NLP routine will converge to one of the K extrema. Label them 1, 2, 3,...K and assign them different colors. Then make a scatter plot of the inial conditions, and color each point by the extremum that it converged to. To be concrete, suppose that your function has three maximums: m1, m2, and m3, and that you've identified these location. You choose a bunch of initial conditions on a rectangular grid and record which max each one converges to. For example, on a 11x11 integer grid, do the following: x=.; y=.; m=.; create results var {x y m}; do x = -5 to 5; do y = -5 to 5; /** use NLP function to find optimum, z **/ /** if z is near m1, then m=1; if z is near m2, then m=2; etc; **/ append; end; end; close results; PROC SGPLOT data=results; scatter x=x y=x / group=m; run; I have not verified that the previous code compiles/runs, but I hope this gets you started.

Rick_SAS · ‎07-26-2011

Obviously, the categories need to be ordinal so that you can order them. Maybe you have levels such as "bad", "ok", "good", and "excellent"? The way that you convert these to numbers will affect your answer. You can recode those levels as 1,2,3,4, or you might decide that "bad" and "excellent" have more extreme values, and so you recode as 0, 2,3,5. There is some terse discussion of this in the doc for the SCORES statement of the FREQ procedure: http://support.sas.com/documentation/cdl/en/procstat/63963/HTML/default/viewer.htm#procstat_freq_a0000000560.htm#procstat.freq.freqscores For a technical consideration of whether you should do this, along with instructions how, see Analysis of Ordinal Categorical Data by Alan Agresti. You can also try an internet seach using terms such as Ordinal Categories Rank Ridit Scores.

Rick_SAS · ‎07-25-2011

You can do scoring and much more with the new PROC PLM procedure. See http://support.sas.com/resources/papers/proceedings10/258-2010.pdf An example of scoring starts on p.7.

Rick_SAS · ‎07-21-2011

You can format a variable at run-time by using the PUTN function: x = '14Jul2011'd; y = putn(x, "DATE9.");

Rick_SAS · ‎07-20-2011

Several SAS procedure support the "GLM notation" that enables you to specify a model that includes main effects and all two-way interactions. You use the "|" operator to specify interactions and the "@2" operator to limit the interactions to tw0-way. For example, a GLM analysis might be proc glm data=a; class s1-s7; model y = s1|s2|s3|s4|s5|s6|s7 @2; run; If the procedure you want to use does not support the GLM notation, you can use PROC GLMMOD with the OUTDESIGN= option to generate the design matrix. For more on specifying interactions, see http://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_introcom_a0000002993.htm

Online Status	Online
Date Last Visited	28m ago