Hello @A_Kh,   The documentation of the PANELBY statement says that the COLHEADERPOS= option "has no effect if the panel uses the PANEL layout" (which is the default of the LAYOUT= option). So use a different layout: COLUMNLATTICE looks very similar to you current layout. panelby origin/novarname columns=3 rows=1 colheaderpos=bottom layout=columnlattice; ... View more

Hi @Balli,   @Balli wrote: (...) I think this can be achieved using first and last variables ... Here is an approach using the FIRST.STATE and LAST.STATE variables: The observation numbers of the first and last observation of each state in dataset HAVE_A (assuming they are grouped) are stored in a temporary lookup table (hash object). While reading dataset HAVE_B they are retrieved and a (uniformly distributed) random number between them is selected using the RAND function. This random number, in turn, is used in the POINT= option of the SET statement retrieving a phone number from dataset HAVE_A. /* Create sample data for demonstration */ %let states='NSW','QLD','VIC','WA','ACT','NT','SA','TAS'; data have_A; do State=&states; do _n_=1 to 23-mod(rank(md5(State)),7); Phone=ranuni(2718)*(2**31-1); output; end; end; run; data have_B; call streaminit(27182818); do _n_=1 to 40000; AccountID=ranuni(3142)*(2**31-1); length State $3; State=choosec(rand('integer',8),&states); output; end; run; proc sql; insert into have_B values(123456789,'XYZ'); quit; /* Assign a randomly selected phone number from have_A to each account in have_B, based on their state */ data want(drop=n1 n2); call streaminit(31415927); if _n_=1 then do; if 0 then set have_B have_A; dcl hash h(); /* lookup table for the first and last obs. number per state in have_A */ h.definekey('State'); h.definedata('n1','n2'); h.definedone(); do _n_=1 by 1 until(last); set have_A end=last; by state notsorted; /* Dataset have_A must be grouped by State. */ if first.state then n1=_n_; if last.state then do; n2=_n_; h.add(); end; end; end; set have_b; if h.find()=0 then do; _n_=rand('integer',n1,n2); set have_A(keep=phone) point=_n_; end; else call missing(phone); run; Alternatively, you could store the phone numbers in the hash object and retrieve them from there (e.g., using a sequence number as a second key item).   But I think a solution using PROC SURVEYSELECT (@Ksharp's idea) is easier to understand.   Edit: Inserted the ELSE statement at the end of the program to avoid retaining a phone number from the previous observation if an unexpected value of State occurs. ... View more

@Ksharp wrote: proc surveyselect data=tablea out=phone seed=123 sampsize=SampleSize selectall; strata state; run; Great idea, @Ksharp, to use PROC SURVEYSELECT, but I'm pretty sure you actually wanted to use the options method=urs outhits outrandom instead of selectall. ... View more

Hello @k26_23,   Try this: /* Create temporary variable describing block structure */ data _tmp / view=_tmp; set have; by subject site; _block=ifn(first.site,1,sum(_block, (id=' ') ne (lag(id)=' '))); retain _block; run; /* Use the temporary variable to obtain the desired sort order */ proc sql; create table want(drop=_block) as select * from _tmp order by subject, site, _block, id, sort_order; drop view _tmp; quit; ... View more

Hello @sasgorilla,   I would probably create a temporary variable the sort order of which allows for using a HAVING clause selecting the minimum (or maximum) value of that variable (as a "tie-breaker"). So the new variable would reflect your criterion "the observation with the lowest b, c and d."   Here is an example using _x=b+c+d for that purpose: proc sql; create table merged(drop=_x) as select s2.id, eventstart, physcat, a, b, c, d, e, b+c+d as _x from set2 as s2 left join set1 as s1 on s1.id = s2.id and eventstart between d_start and d_end group by s2.id, eventstart having _x=min(_x); quit; If your interpretation of "the lowest b, c and d" is rather hierarchical, i.e., lowest b is most important, followed by lowest c (if the values of b are equal) and so on, you could define _x as something like 10000*b+100*c+d (assuming non-missing integer values b, c, d between 0 and 99). And many more definitions are possible (e.g., adding criteria for missing values). ... View more

I noticed a similar change when I received the most recent setinit file for my personal SAS license in February: The "grace period" was cancelled (i.e., set to zero). It had been constantly 62 days for the past ten years. The "warning period" (which follows the grace period) was cut down to 15 days (from 31 - 34 in previous years).   So I would have entered the shortened warning period two months earlier if I hadn't renewed the license in time. ... View more

Hi @twedwards and @FoxHope,   Thank you for providing evidence that this phenomenon still occurs with recent SAS versions. I observed it in May 2012 when I worked on a freelance project at a pharma company. My documentation indicates that it occurred with SAS 9.2 TS2M3 on their Windows 2008 server, but not on their Unix server. Interestingly, that Windows Server 2008 was a 32-bit operating system (the "W32_SRV08 platform"), which might be crucial for this issue.   Here is an example that I documented back then: data test; a=0.1; b=0.7; c=a+b; d=0.8; e=d-c; f=d-c; proc print; run; Result (from May 2012): Obs A B C D E F 1 0.1 0.7 0.8 0.8 8.3267E-17 1.1102E-16 As in your examples, it was always only the first of two or more identical calculations which had an unexpected result. The result in variable F was the expected "classic numeric precision" issue, a rounding error in the least significant bit, which can be explained in detail by performing the calculation in the binary system according to the IEEE 754 standard. Only the result in E was surprising. (One of the striking consequences of the issue was that, after defining a=0.1; b=0.7; c=a+b; d=0.8; the seemingly exhaustive IF condition d-c>1e-16 | d-c<1e-16 | d-c=1e-16 was not met, presumably because the result of d-c changed from 8.3267E-17 in the first part to 1.1102E-16 in the second and third part of the Boolean expression!)   Today, my SAS 9.4M5 on a 64-bit Windows workstation yields 1.1102E-16 in both variables, E and F, as expected. Similarly, it yields correctly 0 in your examples.   I agree with you and @Quentin that this looks very much like what is described in the Usage Note 6214 (which I wasn't aware of back in 2012) linked by @Kathryn_SAS.   A strong indication that it has to do with those "12 extra bits for the mantissa" (in FPU registers or somewhere else in the processor) mentioned in that Usage Note is this: The least significant bit of the internal 64-bit floating-point representation of 1/3 in a numeric variable (under Windows or Unix, including 32-bit platforms) has a place value of 2**-54=5.5511...E-17. Therefore, any "classic" numeric precision issue involving numbers close to 1/3 (namely >=0.25 and <0.5) can only create differences which are integer multiples of 2**-54. The absolute value of your example difference r=-1.8513E-17, however, is obviously less than that. In fact, it is close to one third of 2**-54 (1.85037...E-17) and an internal binary representation of it would require at least 11 additional mantissa bits: -683*2**-65 = -1.85127520951899882...E-17.   I would be curious if your difference r in BEST32. format looks like -1.851275209519E-17 and, even more precisely, like BC75580000000000 in HEX16. format. ... View more

@SunnySong wrote: How come I can only take one post as a solution? You have truly helped me as well. I'm going to give Tom's post an extra like to compensate for it not being the accepted solution. :-) ... View more

Thanks, @Tom, for the reminder.   Yes, I would say using %let MESE_TEXT = %sysfunc(putn("&DATA_PERIO"d, itadfmn));  is even better than my suggested nldatemn, as it doesn't rely on the LOCALE setting.   It just seems that SAS discourages us from using those language-specific format names. First, they are a bit hard to find in the documentation: They are not contained in the "Dictionary of Formats" and not mentioned in the "Dictionary of Formats for NLS" either. A few occurrences (including itadfmn) can be found in the "Example" sections of the documentation of the EUR... formats such as EURDFMNw. in the appendix "Additional NLS Language Elements". But the introductory section of that appendix says:   The following EUR language elements have been replaced with NL language elements. The EUR elements are supported in SAS 9.3 [sic!], but SAS recommends that you use the NL elements.   ... View more

Hello @VioletaB,   If you have two independent samples with a location shift between their (otherwise ideally equal) distributions, I would suggest the TWOSAMPLEWILCOXON statement of PROC POWER. Note, however, that the expected median difference and estimated standard deviation (together with power and significance level) do not provide sufficient information for the sample size calculation. Sample size for the Wilcoxon rank-sum test also depends on the shape of the distribution. This is why there is the VARDIST option where you need to specify the distributions of the two samples.   Example: Normal distributions with SD=1.2, shifted by |m 1 −m 2 |=1 proc power; twosamplewilcoxon vardist('Group 1')=normal(0, 1.2) vardist('Group 2')=normal(1, 1.2) variables = 'Group 1' | 'Group 2' sides = 2 alpha = 0.05 power = 0.85 npergroup = .; run; Result: n=28 per group.   With the more heavy-tailed Laplace distribution vardist('Group 1')=laplace(0, 0.848528) vardist('Group 2')=laplace(1, 0.848528) the result is only n=21 per group, although the median difference and standard deviation (sqrt(2)*0.848528...=1.2) are the same as for the normal distributions above. ... View more

[continuation of the previous post]   Now let's turn to your PROC POWER step where you obtained a sample size of 69 subjects per group. The corresponding power estimate is even more suitable for comparisons: ods output output=ppow; proc power; twosamplemeans test=equiv_diff dist=normal lower = -0.223 upper = 0.223 meandiff = 0.05 stddev = 0.4 npergroup = 69 power = .; run; proc print data=ppow; format power best16.; run; /* 0.80179614325271 */ /* The same result (up to 14 decimals!) can be obtained by using a formula from section 3.2.3 of the book (p. 62): */ data _null_; a=0.05; /* significance level alpha */ dU=0.223; /* bioequivalence limit (log(1.25)) */ d=0.05; /* mean difference delta */ s=0.4; /* common standard deviation in both groups */ n=69; /* sample size per group */ power=1-cdf('T',quantile('T',1-a,2*n-2),2*n-2,(dU-d)/(s*sqrt(2/n))) -cdf('T',quantile('T',1-a,2*n-2),2*n-2,(dU+d)/(s*sqrt(2/n))); put power= best16.; run; The power computed above is that of an equivalence test based on a much simpler model than model (10.2.1) in the book: Just two independent normal distributions (ideally) with equal variances (which are assumed to be unknown). No cross-over design, just overall mean plus fixed treatment effect plus normal random error ("two-sample parallel design").   /* Simulate data of 500,000 experiments in a two-sample parallel design */ %let n=69; /* sample size per group */ %let mu=1; /* overall mean, arbitrarily selected */ %let d=0.05; /* fixed treatment effect */ %let s=0.4; /* common standard deviation in both groups */ data sim2(drop=F:); call streaminit('MT64',31415927); retain sim k i y e; array F[*] FT FR (%sysevalf(&d/2) -%sysevalf(&d/2)); do sim=1 to 500000; do k=1 to 2; do i=1 to &n; e=rand('normal',0,&s); y=&mu+F[k]+e; output; end; end; end; run; /* Perform equivalence test */ options nonotes; ods select none; ods noresults; ods output EquivLimits=eqlim2(where=(method='Pooled')); proc ttest data=sim2 tost(-0.223, 0.223); by sim; class k; var y; run; options notes; ods select all; ods results; proc freq data=eqlim2; tables assessment / binomial; run; /* proportion assessment="Equivalent", i.e., estimated power: 0.8011 (95%-CI: [0.8000, 0.8022]) */ Again, the result of the simulation matches the computed power 0.8017... obtained before. ... View more

As promised, here is the SAS code I developed during my investigation.   /* Simulation of 500,000 clinical trials with 2x2 crossover design to establish average bioequivalence between two treatments using log-transformed data (cf. the example on p. 262 f. of Chow, S.C., Shao, J., and Wang, H. (2008): Sample Size Calculations in Clinical Research, 2nd ed., Chapman & Hall/CRC Press/Taylor & Francis Group, New York and the derivation of the sample size estimate in this book) */ %let n=18; /* sample size per group */ %let mu=1; /* overall mean in the statistical model (10.2.1) used in the book, arbitrarily selected */ %let d=0.05; /* fixed treatment effect ("delta" in the book) */ %let rho=0.5; /* Between-subject std. deviations (sBT, sBR) and within-subject std. dev. (sWT, sWR), */ %let sBT=0.35; /* denoted "sigma_BT" etc. in the book, as well as the correlation rho of each subject's random */ %let sBR=0.25; /* effect between test (T) and reference (R) treatment have been arbitrarily selected, but */ %let sWT=0.2; /* aiming at the value 0.4 for the derived standard deviation parameter s11 (see below) used in */ %let sWR=0.15; /* the book example. */ %let P1=0.07; /* fixed effect of the first treatment period */ %let Q1=0.03; /* fixed effect of the first sequence of treatments (T, R) */ /* Alternative set of parameters: %let n=18; %let mu=1; %let d=0.05; %let rho=0.15; %let sBT=0.10; %let sBR=0.05; %let sWT=0.37; %let sWR=0.11; %let P1=0.12; %let Q1=0.13; */ %let s11=%sysfunc(sqrt(&sBT**2+&sBR**2-2*&rho*&sBT*&sBR+&sWT**2+&sWR**2)); /* in the book: "sigma_1,1", described */ %put &=s11; /* =0.4 */ /* as "standard deviation for intra-subject comparison" */ data sim(drop=FT--sWR); call streaminit('MT64',27182818); retain sim k i j l S e y; array F[*] FT FR (%sysevalf(&d/2) -%sysevalf(&d/2)); array P[2] (&P1 -&P1); array Q[2] (&Q1 -&Q1); array S_[*] ST SR; array sW[*] sWT sWR (&sWT &sWR); do sim=1 to 500000; /* Reduce this number for shorter run times. */ do k=1 to 2; /* treatment sequence k=1: T, R; sequence k=2: R, T */ do i=1 to &n; ST=rand('normal',0,&sBT); /* i-th subject's random effect; its bivariate normal distribution is simulated ... */ SR=rand('normal',%sysevalf(&rho*&sBR/&sBT)*ST,%sysevalf(&sBR*%sysfunc(sqrt(%sysevalf(1-&rho**2))))); /* ... using the conditional distribution of SR given ST; alternatively, PROC SIMNORMAL could be used. */ do j=1 to 2; /* first and second treatment period */ l=2-(j=k); /* Values l=1, 2 correspond to treatments T and R, respectively. */ S=S_[l]; e=rand('normal',0,sW[l]); /* random errors */ y=&mu+F[l]+P[j]+Q[k]+S+e; /* see equation (10.2.1) in the book */ output; end; end; end; end; run; /* Compute mean values of the (log-transformed) measurements per period-sequence combination */ proc summary data=sim nway; by sim k; class j; var y; output out=stats(drop=_:) mean=; run; /* Estimate the treatment effect ("delta hat" in the book) */ data delta(keep=sim d_hat); set stats; by sim; if first.sim then d_hat=0; d_hat+((-1)**(j~=k)*y); if last.sim; d_hat=d_hat/2; run; /* Optional: Check the (normal) distribution of "delta hat" (formula (10.2.3) in the book) proc univariate data=delta; var d_hat; histogram / normal; run; */ /* Prepare estimation of "sigma_1,1 hat" squared */ proc transpose data=stats out=stats2(drop=_:) prefix=y; by sim k; id j; var y; run; /* Estimate "sigma_1,1 hat" squared using formula (10.2.5) of the book */ data s112(keep=sim s112); merge sim(keep=sim k i j y) stats2; by sim k; if first.sim then s112=0; lag_y=lag(y); if j=2 then s112+(lag_y-y-y1+y2)**2; if last.sim; s112=s112/(2*&n-2); run; /* Optional: "sigma_1,1 hat" squared and "delta hat" are independent; check uncorrelatedness and distribution of "sigma_1,1 hat" (cf. p. 261 of the book) data chk; merge delta s112; by sim; run; proc corr data=chk; var d_hat s112; run; data s112n; set s112; s112n=s112*(2*&n-2)/&s11**2; * This quantity should have a chi-square distribution with 2*n-2 degrees of freedom.; run; proc univariate data=s112n; var s112n; histogram / gamma(alpha=%eval(&n-1) sigma=2); * chi-square(df) distribution equals gamma(df/2, 2) distribution; run; */ /* Compute 90% confidence interval for delta and perform the equivalence test ("two one-sided tests procedure") using delta_U=log(1.25)=0.223... and delta_L=log(0.8)=-delta_U (bioequivalence limits as per 2001 FDA guidance "Statistical Approaches to Establishing Bioequivalence" [https://www.fda.gov/media/70958/download]) */ data confd; merge delta s112; by sim; hw=quantile('T',0.95,2*&n-2)*sqrt(s112)/2*sqrt(2/&n); /* half-with of the CI according to book p. 261 */ c=(abs(&d-d_hat)<=hw); /* =1 in case of coverage of the true parameter, 0 otherwise */ rej=(d_hat-hw>log(0.8) & d_hat+hw<log(1.25)); /* =1 if null hypotheses H_01 and H_02 are rejected */ run; /* at the 5% significance level */ /* Check coverage of the 90% CI for delta and estimate the power of the equivalence test */ proc freq data=confd; tables c rej / binomial(level='1'); run; /* proportion c=1: 0.8994 (95%-CI: [0.8985, 0.9002]) prop. rej=1: 0.8074 (95%-CI: [0.8063, 0.8085]) */ /* with &n=17: */ /* proportion c=1: 0.8995 (95%-CI: [0.8987, 0.9003]) prop. rej=1: 0.7837 (95%-CI: [0.7825, 0.7848]) */ The above estimates of the test's power indicate that n = 18 is the minimum sample size required for 80% power, given the assumed values of parameters d and s 1,1 . Very similar results are obtained with the "alternative set of parameters" (see that comment in the code).   /* Suggested better approximation of the power than the book's 2Phi(...)-1 formula (p. 261) */ data _null_; do n=17, 18; pow= cdf('normal', log(1.25)-quantile('T',0.95,2*n-2)*sqrt(&s11**2/(2*n)),&d,sqrt(&s11**2/(2*n))) -cdf('normal',-log(1.25)+quantile('T',0.95,2*n-2)*sqrt(&s11**2/(2*n)),&d,sqrt(&s11**2/(2*n))); put n= "--> approx. power=" pow 6.4; end; run; The above approximation yields 0.8095 for n = 18 and 0.7857 for n = 17, suggesting the same minimum sample size as the simulation.   As mentioned in my previous post, the equivalence test can also be implemented with PROC TTEST rather than in several steps as above: /* Restructure simulated data for application of PROC TTEST */ proc transpose data=sim out=simt(drop=_:) prefix=y; by sim k i; id j; var y; run; data simtt; set simt; trt1=k; trt2=3-k; run; /* Perform the equivalence test with PROC TTEST */ options nonotes; ods select none; ods noresults; ods output EquivLimits=eqlim(where=(period=' ' & method='Pooled')); proc ttest data=simtt tost(%sysfunc(log(0.8)), %sysfunc(log(1.25))); by sim; var y1 y2 / crossover=(trt1 trt2); run; options notes; ods select all; ods results; proc freq data=eqlim; tables assessment / binomial; run; /* Compare previous results to those from PROC TTEST */ data cmptt; merge confd(keep=sim rej) eqlim(keep=sim assessment); by sim; run; proc freq data=cmptt; tables rej*assessment; run; /* Results are identical: rej=1 <==> Assessment="Equivalent" */   My suggestion for an exact calculation of the power is to compute the integral of the joint probability density of "d hat" and "s 1,1  hat" over the rejection region, i.e., the set of ("d hat" and "s 1,1  hat") pairs satisfying the inequality |"d hat"| < d U  − t a,2n−2 "s 1,1 hat"/√(2n) where the two null hypotheses H 01 and H 02 are rejected. As those two random variables are independent (see book, p. 261), their joint pdf is just the product fg of the pdfs of a normal distribution (f ) and a scaled chi distribution (g ). I found the pdf of a chi distribution in section 8.3 of Evans, M., Hastings, N. and Peacock, B. (2000), Statistical Distributions, 3rd ed., John Wiley & Sons, New York. %let n=18; %let d=0.05; %let s11=0.4; %let s112=%sysevalf(&s11**2); %let alpha=0.05; /* Select a subset of points where the joint pdf does not almost vanish */ data intset; do i=-5000 to 6000; x=i/10000; do j=1000 to 8000; y=j/10000; fg= sqrt(&n/constant('PI'))/&s11*exp(-&n/&s112*(x-&d)**2) *2*((&n-1)/&s112)**(&n-1)/fact(&n-2)*y**(2*&n-3)*exp(-(&n-1)/&s112*y**2); if fg>1e-10 then output; end; end; run; /* 52891353 obs. */ proc means data=intset min max; var i j; run; /* Min and max are well inside the ranges, so it's safe to limit the integration to this set. */ /* Compute the integral numerically */ proc means data=intset sum; where abs(x)<log(1.25)-quantile('T',1-&alpha,2*&n-2)*y/sqrt(2*&n); var fg; run; /* --> Power = 80703401.67/1e8 = 0.80703... . */ This value matches the result of the simulation, 0.8074 (95%-CI: [0.8063, 0.8085]), quite well.   I will continue in a subsequent post to avoid that this one is getting too long. ... View more

Hi, Agnieszka (@A_Tomczyk), and welcome to the SAS Support Communities!   Thank you very much for posting this interesting question. I have investigated the issue using the second edition (2008) of your book (in which I found basically the same texts and formulas as shown in your photographs of the third edition, so I assume the preceding paragraphs of that section "10.2 Average Bioequivalence" are very similar as well).   Everything seemed relatively straightforward to me (note the typo −d L , which should read d L ), until I stumbled across the sentence "Under the alternative hypothesis that |e| < d, ..." ending with an approximative power formula of the form 2F(...)−1. As it turned out, the authors switched the notation there: It would have been more consistent to write that alternative hypothesis like |d| < d U  .   Also, the approximation 2F(...)−1 seems unnecessarily coarse to me: I say "unnecessarily" because their equations "do not have an explicit solution" (quote from p. 261) anyway, so they could also have suggested a better approximation such as F(d U  −t a,2n-2  s 1,1 /√(2n)) − F(−d U  +t a,2n-2  s 1,1 /√(2n)) with the cumulative distribution function F of "d hat", i.e., the cdf of a normal distribution (see formula 10.2.3 in the book if the numbering is the same in your edition). I use this in my code (which I will post separately). It is still an approximation because s 1,1 is used instead of its estimate "s 1,1 hat", but unlike theirs it does not require |d| (i.e., their |e|) to be very small. Using their approximation I sometimes got nonsensical negative values for the estimated power.   For an exact computation you can evaluate an integral of the joint probability density of "d hat" and "s 1,1 hat", which is not very difficult to do in a DATA step (also shown in my code) thanks to the independence of these two random variables.   In the next step of their derivation, the authors of the book use the upper b/2 quantile of the t distribution with 2n-2 degrees of freedom, although the corresponding normal quantile z b/2 would better fit the previous "2F(...)−1" expression. But, as we know, those two quantiles don't differ much if n is sufficiently large. (And they do use z b/2 in their "further simplified" sample size formula.)   Surprisingly, they made two mistakes in their example calculation (which you show in one of your photographs): z 0.05 =1.64..., not 1.96, and z 0.10 =1.28..., not 0.84. Using the correct quantiles, the result for the (approximate) estimated minimum sample size is 23, not 21. Interestingly, using the "better approximation" I've suggested further above, I obtain n = 18 -- and this is even confirmed by both the exact calculation (using the integral mentioned above) and the simulation of 500,000 cross-over studies (each for n = 17 and n = 18) that I did in my code.   The relevant hypothesis test can be implemented with PROC TTEST using the TOST option of the PROC TTEST statement in conjunction with the CROSSOVER= option of the VAR statement (also shown in my code). However, it seems to me that the corresponding sample size calculation is beyond the capabilities of PROC POWER (at least in my SAS/STAT 14.3 version), which does support calculations corresponding to using either TOST or CROSSOVER= in PROC TTEST. Your PROC POWER code uses the "TOST only" variant. In my code I show the corresponding PROC TTEST step as well as a simulation.   The GLMPOWER procedure doesn't help either (I think) because in the documentation it says: "Tests and contrasts that involve random effects are not supported." But the model used in Chow/Shao/Wang does involve a random effect.   So, unless someone else has other ideas, you can use one of the approximations from the Chow et al. book or the "better" approximation mentioned above or the exact integral calculation or a simulation  to obtain power and sample size for the model in question with SAS.   I am so sorry that I have to leave the office now, but I will post my code tomorrow (CEST time zone). ... View more

@PanShuyao wrote: I think RANNOR and RANGAM use the older RNG. This is true, but I'm not surprised that the results still don't match. The (pseudo-)random numbers depend on the order of operations and at least I don't know how PROC MI works internally, so it is not obvious how to replicate its computations manually.   Please see the addendum I added to my first reply for a more promising approach. ... View more