Hi @Rick_SAS !  Thank you for following up.  Yes, it looks like the reference that mentioned standard deviation uses the geometric mean in their BC formula, and your example didn't.  Interestingly enough, in your blog post from last year, you added the geometricmean option on your PROC TRANSREG statement.  I think that only affects the final transformed values of Y, but not the way PROC TRANSREG computes the optimal lambda.  I'll have to see if I can find the Draper and Smith paper the PROC TRANSREG code is based on.  They don't seem to be using the log likelihood for either the geometric mean version or the non-geometric mean version of the BC formula.  When I used the geometric mean BC formula log likelihood on your example data, I got the same lambda result as PROC TRANSREG (-0.68), and the maximum log likelihood occurred for the same lambda (again, -0.68) as the minimum standard deviation (computed using N rather than N-1).   By the way, in my data step code I got the lambda = 0 term wrong, dividing by the geometric mean when I should have multiplied.  I got a discontinuity in the output.  But it wasn't because of that error -- my code never got to the lambda = 0 subsection, because the loop of -2 to +2 by 0.01 is never exactly zero in SAS, so it must have used the power law for a very small absolute value of lambda.  Live and learn :-). ... View more

Hi @WarrenKuhfeld !  Thank you for the great tip.  I ran @Rick_SAS 's PROC TRANSREG example with and without the geometricmean option, and printed out the output.  The lambda is the same and all the log likelihood values are the same, so internally it seems to be doing the same thing to arrive at the optimal lambda.  The only thing different is the set of final transformed Y values, which differ by the geometric mean constant (raised to the appropriate power). ... View more

Hi @WarrenKuhfeld !  Thank you for the reference.  Although my log likelihood numbers are different from the ones in PROC TRANSREG, I'm at least reassured that on the example I ran, I got the same answer for the BC lambda as the one from PROC TRANSREG. ... View more

Hi @Rick_SAS !  Thank you once again for taking the time to assist me.  Mucho appreciated, and I totally understand if you're done with this for now.   I would like to post some follow-up.   1.  You mentioned my "original question" and my "latest question."  Sorry if I was unclear, but my question didn't really change.  I'm still trying to figure out whether Box-Cox is equivalent to minimizing the standard deviation.  I'm now pretty sure that it IS equivalent, once we understand which standard deviation we are referring to:  the STD for the BC formula WITH THE GEOMETRIC MEAN.   2.  No, my log likelihood formula was not missing a factor of 1/2.  It's just that I wrote it a little differently than you usually see it.  I wrote:   -(N * ln(sigma) )  Usually you see:  -(N * 0.5 * ln(sigma**2)).  They're the same.  And my formula produced the same numbers as the SAS LOGPDF function did.   3.  I don't think your example proves what you think it does concerning standard deviation.  You didn't include the geometric mean term in your BC formula.  Not everyone does.  It's just a constant.  So, you're correct about the standard deviation without the geometric mean term.  But if you don't somehow include the geometric mean in the search for lambda I don't think you'll get the right result.  If you include the geometric mean term, then the lambda that maximizes the log likelihood is also the lambda that minimizes the maximum likelihood standard deviation estimate (which differs from the usual unbiased STD estimator by using N instead of N-1).  When I run your data by hand, I still get different values of the log likelihood than PROC TRANSREG reports, but the maximum occurs at -0.68, just like PROC TRANSREG found.  I don't think that's a coincidence.  If I run without the geometric mean, the maximum likelihood over the range from -2 to +2 occurs at -2, and so does the minimum standard deviation.  I don't think that's the correct BC lambda, and if we tried a wider range we'd get an even lower value.   This has been a really helpful discussion for me.  I still would like to find out why my log likelihood numbers are different from the ones in PROC TRANSREG.  Perhaps the references you provided will help there.  And I see that @WarrenKuhfeld provided the reference for where he got the method, so I'll try to locate that one.  My code to run log likelihoods from lambda -2 to +2 is below, both with and without the geometric mean term.   data have_llln_search_2 ; set have end = last ; keep lambda lambda_mean lambda_std_ml lambda_ll_n lambda_n_logsigml lambda_u_mean lambda_u_std_ml lambda_u_ll_n lambda_u_n_logsigml ; array ty{1:10} _temporary_ ; array tyul{1:10} _temporary_ ; array tyl{1:10} _temporary_ ; ty{_N_} = Y ; if last then do ; n_obs = _N_ ; n_fact = sqrt((n_obs-1)/n_obs) ; y_gmn = geomean(of ty{*}) ; lln = n_obs * log(2 * constant("pi")) / 2 ; min_yl_stdn = 1000 ; max_yl_lln = -1000 ; lambda_minimize = 999 ; lambda_maximize = 999 ; min_yul_stdn = 1000 ; max_yul_lln = -1000 ; lambdau_minimize = 999 ; lambdau_maximize = 999 ; do lambda = -2 to +2 by 0.01 ; if lambda ne 0 then do ; do ix = 1 to n_obs ; tyl{ix} = ((ty{ix})**(lambda) - 1) / (lambda * (y_gmn)**(lambda-1)) ; tyul{ix} = ((ty{ix})**(lambda) - 1) / lambda ; end ; end ; else do ; do ix = 1 to n_obs ; tyl{ix} = (log(ty{ix}) / y_gmn) ; tyul{ix} = log(ty{ix}) ; end ; end ; yl_mean = mean(of tyl{*}) ; yl_std = std(of tyl{*}) ; yl_stdn = yl_std * n_fact ; yl_lln = 0 ; yul_mean = mean(of tyul{*}) ; yul_std = std(of tyul{*}) ; yul_stdn = yul_std * n_fact ; yul_lln = 0 ; do ix = 1 to n_obs ; yl_lln + logpdf('NORMAL',tyl{ix},yl_mean,yl_stdn) ; yul_lln + logpdf('NORMAL',tyul{ix},yul_mean,yul_stdn) ; end ; lambda_mean = yl_mean ; lambda_std_ml = yl_stdn ; lambda_ll_n = yl_lln ; lambda_n_logsigml = -n_obs * log(yl_stdn) ; lambda_u_mean = yul_mean ; lambda_u_std_ml = yul_stdn ; lambda_u_ll_n = yul_lln ; lambda_u_n_logsigml = -n_obs * log(yul_stdn) ; output ; if (yl_stdn < min_yl_stdn) then do ; min_yl_stdn = yl_stdn ; lambda_minimize = lambda ; end ; if (yl_lln > max_yl_lln) then do ; max_yl_lln = yl_lln ; lambda_maximize = lambda ; end ; if (yul_stdn < min_yul_stdn) then do ; min_yul_stdn = yul_stdn ; lambdau_minimize = lambda ; end ; if (yul_lln > max_yul_lln) then do ; max_yul_lln = yul_lln ; lambdau_maximize = lambda ; end ; end ; put min_yl_stdn= lambda_minimize= ; put max_yl_lln= lambda_maximize= ; put min_yul_stdn= lambdau_minimize= ; put max_yul_lln= lambdau_maximize= ; end ; run ; title2 "proc print data = &syslast." ; proc print data = &syslast. ; run ; title2 ;   ... View more

Hi @Rick_SAS !  Thank you so much for responding.  I ran your example, but I guess it would help me to see the equations, because when I tried to replicate some of the log-likelihood values by hand, I got very different results.  First of all, is this the normal log likelihood formula?   -(N * (0.5) * ln(2 * pi))  -(N * ln(sigma))  -((0.5 * sum((x[j] - mu)**2) / (sigma**2))   where x[j] are the data values, N is the number of values, and mu and sigma are the parameters to determine.   Then, if mu = (sum(x[j]) / N)  and  sigma_ml = sqrt(sum((x[j] - mu)**2) / N), when you substitute back into the log likelihood you get:   -(N * (0.5) * ln(2 * pi))  -(N * ln(sigma_ml))  -(0.5 * N)   In this case, for each lambda, the x[j] values are the lambda transformed values of the original set of x values.  I tried the logpdf function in SAS, and got numbers corresponding to the formula I gave above.  But that isn't what came out of PROC TRANSREG for your example.  For your original set of Y values, I get ll = -17.4666119.  But PROC TRANSREG says ll = -3.80403 for lambda = +1.  For your Ym2, I got ll = 51.408472367 and for your Y2, I got ll = -40.72112937.  As you remark, the BC formula without the geometric mean term puts the data on different scales for different lambdas.  But if I use the geometric mean term in the BC formula for your example data, for lambda = -2 I get ll = -17.38343536 and for lambda = +2 I get ll = -17.79049346, which are in the same ballpark.  For lambda = -2 PROC TRANSREG says ll = -3.72085 and for lambda = +2, PROC TRANSREG says ll = -4.12791.  My code is below.   When you look at the log likelihood formula with the maximum likelihood parameter value substitutions, it seems clear that the lambda that minimizes the value of sigma_ml will maximize the log likelihood.  So, what am I getting wrong?     data handy ; set have end = last ; keep n_obs n_fact y_gmn lln y_mean y_std y_stdn y_ll y_lln ym_mean ym_std ym_stdn ym_ll ym_lln yp_mean yp_std yp_stdn yp_ll yp_lln yl_mean yl_std yl_stdn yl_ll yl_lln yml_mean yml_std yml_stdn yml_ll yml_lln ypl_mean ypl_std ypl_stdn ypl_ll ypl_lln ; array ty{1:10} _temporary_ ; array tym{1:10} _temporary_ ; array typ{1:10} _temporary_ ; array tyl{1:10} _temporary_ ; array tyml{1:10} _temporary_ ; array typl{1:10} _temporary_ ; ty{_N_} = Y ; tym{_N_} = Ym2 ; typ{_N_} = Y2 ; if last then do ; n_obs = _N_ ; n_fact = sqrt((n_obs-1)/n_obs) ; y_gmn = geomean(of ty{*}) ; lln = n_obs * log(2 * constant("pi")) / 2 ; do ix = 1 to n_obs ; tyl{ix} = ty{ix} - 1 ; tyml{ix} = ((ty{ix})**(-2) - 1) / (-2 * (y_gmn)**(-3)) ; typl{ix} = ((ty{ix})**(2) - 1) / (2 * (y_gmn)**(+1)) ; end ; y_mean = mean(of ty{*}) ; y_std = std(of ty{*}) ; y_stdn = y_std * n_fact ; ym_mean = mean(of tym{*}) ; ym_std = std(of tym{*}) ; ym_stdn = ym_std * n_fact ; yp_mean = mean(of typ{*}) ; yp_std = std(of typ{*}) ; yp_stdn = yp_std * n_fact ; yl_mean = mean(of tyl{*}) ; yl_std = std(of tyl{*}) ; yl_stdn = yl_std * n_fact ; yml_mean = mean(of tyml{*}) ; yml_std = std(of tyml{*}) ; yml_stdn = yml_std * n_fact ; ypl_mean = mean(of typl{*}) ; ypl_std = std(of typl{*}) ; ypl_stdn = ypl_std * n_fact ; y_ll = 0 ; y_lln = 0 ; ym_ll = 0 ; ym_lln = 0 ; yp_ll = 0 ; yp_lln = 0 ; yl_ll = 0 ; yl_lln = 0 ; yml_ll = 0 ; yml_lln = 0 ; ypl_ll = 0 ; ypl_lln = 0 ; do ix = 1 to n_obs ; y_ll + logpdf('NORMAL',ty{ix},y_mean,y_std) ; y_lln + logpdf('NORMAL',ty{ix},y_mean,y_stdn) ; ym_ll + logpdf('NORMAL',tym{ix},ym_mean,ym_std) ; ym_lln + logpdf('NORMAL',tym{ix},ym_mean,ym_stdn) ; yp_ll + logpdf('NORMAL',typ{ix},yp_mean,yp_std) ; yp_lln + logpdf('NORMAL',typ{ix},yp_mean,yp_stdn) ; yl_ll + logpdf('NORMAL',tyl{ix},yl_mean,yl_std) ; yl_lln + logpdf('NORMAL',tyl{ix},yl_mean,yl_stdn) ; yml_ll + logpdf('NORMAL',tyml{ix},yml_mean,yml_std) ; yml_lln + logpdf('NORMAL',tyml{ix},yml_mean,yml_stdn) ; ypl_ll + logpdf('NORMAL',typl{ix},ypl_mean,ypl_std) ; ypl_lln + logpdf('NORMAL',typl{ix},ypl_mean,ypl_stdn) ; end ; yl_ll_h = -lln -(0.5*(n_obs-1)) - (n_obs*log(yl_std)) ; yl_lln_h = -lln -(0.5*n_obs) - (n_obs*log(yl_stdn)) ; yml_ll_h = -lln -(0.5*(n_obs-1)) - (n_obs*log(yml_std)) ; yml_lln_h = -lln -(0.5*n_obs) - (n_obs*log(yml_stdn)) ; ypl_ll_h = -lln -(0.5*(n_obs-1)) - (n_obs*log(ypl_std)) ; ypl_lln_h = -lln -(0.5*n_obs) - (n_obs*log(ypl_stdn)) ; put n_obs= n_fact= y_gmn= lln= ; put y_mean= y_std= y_stdn= y_ll= y_lln= ; put ym_mean= ym_std= ym_stdn= ym_ll= ym_lln= ; put yp_mean= yp_std= yp_stdn= yp_ll= yp_lln= ; put yl_mean= yl_std= yl_stdn= yl_ll= yl_lln= ; put yml_mean= yml_std= yml_stdn= yml_ll= yml_lln= ; put ypl_mean= ypl_std= ypl_stdn= ypl_ll= ypl_lln= ; put yl_ll_h= yl_lln_h= ; put yml_ll_h= yml_lln_h= ; put ypl_ll_h= ypl_lln_h= ; output ; end ; run ; title2 "proc print data = &syslast." ; proc print data = &syslast. ; run ; title2 ;   ... View more

Hi @Rick_SAS !  Thank you for the quick response.  And for pointing out where to look.  It seems pretty clear now, but I somehow didn't realize what it meant when I first saw it.  Much appreciated!   Okay, now, apologies if this is a dumb follow-up question, but doesn't finding the lambda that maximizes the normal log likelihood of an intercept only regression amount to the same thing as minimizing the sum of squares of the residuals around the average of the lambda-transformed data values?  If so, would that not be the lambda that produces the minimum standard deviation (of the lambda-transformed data)?   Thanks! ... View more

Hi @WarrenKuhfeld !  Thank you for the quick response.  It's helpful to be told where to look by someone who knows.  Much appreciated. ... View more

Hi @FreelanceReinh !  Thank you so much for your code!  Beautiful, just what I wanted, it ran like a charm and gave me the picture that you displayed.  Much appreciated! ... View more

Hi @Zard !  Thank you so much for your code.  From the picture you posted, it looks like it would do the trick.  There was a small typo ("output" misspelled as "ouput") that the compiler understood how to compensate for, but ultimately I was unable to run it, got the following error, which might be peculiar to my installation:   WARNING: Unsupported device 'SVG' for TAGSETS.SASREPORT13(EGSR) destination. Using default device 'PNG'. NOTE: 18810 bytes written to ... ERROR: Insufficient authorization to access /sasgrid/prd/config/app/Lev1/SASApp/sasgraph.svg. WARNING: Output 'histogram' was not created. Make sure that the output object name, label, or path is spelled correctly. Also, verify that the appropriate procedure options are used to produce the requested output object. For example, verify that the NOPRINT option is not used.   Do you know how I could get around this?   ... View more

Hi @Rick_SAS!   Thank you for responding.  I can see your point about the formula, but I still don't have a good intuitive feel for how the event frequency will affect the ASE and CI size, nor how reliable the ASE and CI are for very low event counts.    I ran some code which I have copied into this message below (I can't upload files, sorry).  It does some testing with 10,000 observations (actually 10,001), one set of "predictions" (ordervar) and one or two events scattered in different "dependent variables" (the s_1_* and s_2_* variables).    For each of the single events, placed at one end or the other end (Somers' D +/- 1) at Q1 or Q3 (SD +/- 0.5) or at the median (SD 0), the confidence intervals are very tight.    But if you drop in a second event the Somers' D can change drastically and the CI can blow wide open.  So the low event count results are not very stable and don't seem trustworthy to me.    I'm wondering whether there is any published guidance on how many events are needed to stabilize the results (like a jackknife test, so that adding or removing one event doesn't completely change the picture).   SELF-CONTAINED CODE:   %let dsn = 02 ; %let loval&dsn. = 0 ; %let hival&dsn. = 10000 ; %let med&dsn. = %sysfunc(floor(%sysevalf((&&hival&dsn.. - &&loval&dsn..) / 2))) ; %let medm1&dsn. = %sysevalf(&&med&dsn.. - 1) ; %let medp1&dsn. = %sysevalf(&&med&dsn.. + 1) ; %let q1&dsn. = %sysfunc(floor(%sysevalf(&&loval&dsn. + ((&&hival&dsn.. - &&loval&dsn..) / 4)))) ; %let q3&dsn. = %sysfunc(floor(%sysevalf(&&hival&dsn. - ((&&hival&dsn.. - &&loval&dsn..) / 4)))) ; %let him1&dsn. = %sysevalf(&&hival&dsn.. - 1) ; %let lop1&dsn. = %sysevalf(&&loval&dsn.. + 1) ; %put med&dsn. = &&med&dsn.. ; %put medm1&dsn. = &&medm1&dsn.. ; %put medp1&dsn. = &&medp1&dsn.. ; %put q1&dsn. = &&q1&dsn.. ; %put q3&dsn. = &&q3&dsn.. ; %put him1&dsn. = &&him1&dsn.. ; %put lop1&dsn. = &&lop1&dsn.. ; /**/ data test_smdcr_&dsn. ; keep ordervar s_0 s_1_l s_1_q1 s_1_m s_1_q3 s_1_h s_2_l_p1 s_2_l_q1 s_2_l_m s_2_l_q3 s_2_l_h s_2_q1_m s_2_q1_q3 s_2_q1_h s_2_m_q3 s_2_m_h s_2_q3_h s_2_m1_h s_2_m1_p1 ; do ordervar = &&loval&dsn.. to &&hival&dsn.. ; s_0 = 0 ; s_1_l = 0 ; s_1_q1 = 0 ; s_1_m = 0 ; s_1_q3 = 0 ; s_1_h = 0 ; s_2_l_p1 = 0 ; s_2_l_q1 = 0 ; s_2_l_m = 0 ; s_2_l_q3 = 0 ; s_2_l_h = 0 ; s_2_q1_m = 0 ; s_2_q1_q3 = 0 ; s_2_q1_h = 0 ; s_2_m_q3 = 0 ; s_2_m_h = 0 ; s_2_q3_h = 0 ; s_2_m1_h = 0 ; s_2_m1_p1 = 0 ; if (ordervar = &&loval&dsn..) then do ; s_1_l = 1 ; s_2_l_p1 = 1 ; s_2_l_q1 = 1 ; s_2_l_m = 1 ; s_2_l_q3 = 1 ; s_2_l_h = 1 ; end ; else if (ordervar = &&lop1&dsn..) then do ; s_2_l_p1 = 1 ; end ; else if (ordervar = &&q1&dsn..) then do ; s_1_q1 = 1 ; s_2_l_q1 = 1 ; s_2_q1_m = 1 ; s_2_q1_q3 = 1 ; s_2_q1_h = 1 ; end ; else if (ordervar = &&medm1&dsn..) then do ; s_2_m1_p1 = 1 ; end ; else if (ordervar = &&med&dsn..) then do ; s_1_m = 1 ; s_2_l_m = 1 ; s_2_q1_m = 1 ; s_2_m_q3 = 1 ; s_2_m_h = 1 ; end ; else if (ordervar = &&medp1&dsn..) then do ; s_2_m1_p1 = 1 ; end ; else if (ordervar = &&q3&dsn..) then do ; s_1_q3 = 1 ; s_2_l_q3 = 1 ; s_2_q1_q3 = 1 ; s_2_m_q3 = 1 ; s_2_q3_h = 1 ; end ; else if (ordervar = &&him1&dsn..) then do ; s_2_m1_h = 1 ; end ; else if (ordervar = &&hival&dsn..) then do ; s_1_h = 1 ; s_2_m1_h = 1 ; s_2_l_h = 1 ; s_2_q1_h = 1 ; s_2_m_h = 1 ; s_2_q3_h = 1 ; end ; output ; end ; run ; /**/ title2 "proc freq data = test_smdcr_&dsn. s_*ordervar cl" ; proc freq data = test_smdcr_&dsn. ; tables (s_0 s_1_l s_1_q1 s_1_m s_1_q3 s_1_h s_2_l_p1 s_2_l_q1 s_2_l_m s_2_l_q3 s_2_l_h s_2_q1_m s_2_q1_q3 s_2_q1_h s_2_m_q3 s_2_m_h s_2_q3_h s_2_m1_h s_2_m1_p1) * ordervar / measures cl noprint ; test smdcr ; output smdcr out = smdcr_test_&dsn. ; ; run ; title2 ;     ... View more

Hi @sbxkoenk !  Thank you for responding.  Computing a bootstrap CI could be good confirmative information, but it doesn't solve the sample size sufficiency issue, does it?  Doesn't a bootstrap still require a certain number of observations to be reliable?  I think I would still need to know how the number of events affects the reliability, if at all. ... View more

Thank you for responding @SASKiwi .  If what you say is true, and as of now I have no reason to doubt it, I do feel disappointed.  After all, Enterprise Guide can communicate between my local drive and the server.  Why can't SAS do it too?  I understand that EG sits on my computer and knows how to locate the grid server because I gave it the grid server address. But then EG should be able to tell SAS where my computer is located; of course, SAS would have to have a way to use that information to reach back from the server to my computer, and it doesn't appear to have that capability.  Is my understanding correct? ... View more