Helmut ★★★ Vienna, Austria, 20061228 22:33 (6414 d 12:16 ago) Posting: # 422 Views: 53,237 

Dear all, following this post and just out of curiosity I implemented SAScode by Jones and Kenward (2000) ^{[1]} in R (v2.4.1).To quote the beginning of Jones'/Kenward's SASCode (which is copyrighted material and therefore not given here; but you can afford buying a book for sure): /*** WARNING : PROGRAM OFFERED FOR USE WITHOUT ANY GUARANTEES ***/ Their example (pp 336337) yields from a withinsubject standard deviation of 0.355 in 58 subjects and a point estimate of 1.00 power=90.4%. Following Rcode gives 90.3988%. # significance level The following code reproduces the powercurves for n=8/12/18/24/30/36/40/60 and CV=20% given in Fig.1c by Diletti et al. (1991) ^{[2]}
a < 0.05 # alpha Values calculated by the Rcode for n=24 match Diletti et al.'s Fig.1c: +++ Please note, this code is not validated, but only a toy you can play around with! With some ambition you can modify the code (essentially varying the point estimate instead of the sample size) in order to match also Diletti et al.'s Fig.2c. Now it should be only a small step to recalculate their table of sample sizes...
— Diftor heh smusma 🖖🏼 Довге життя Україна! _{} Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
Helmut ★★★ Vienna, Austria, 20070102 00:59 (6410 d 09:51 ago) (edited by HS on 20070604 12:41) @ Helmut Posting: # 428 Views: 43,297 

Dear all, I performed some comparisons of the code based on Jones/Kenward implemented in R, Diletti's table and results obtained from StudySize (v2.0.1). Quite interesting... All comparisons were done for CV=20% with ratios of 0.85–1.20. Dilleti reported sample sizes to obtain ≥80% power, calculated odd sample sizes were reported rounded up to the next even number ( underlined ):
+======+=====+===========+===========+ Power with sample sizes given by Diletti et al. at a GMR of 1.05 were below 80%, calculated both with R and StudySize… A Monte Carlo Simulation (1000000 runs) for GMR=1.05 and 18 subjects in StudySize resulted in: Power 0.8006 (95% coverage probability: 0.79980.8013 ).Differences may be due to the implementation of the algorithm to obtain the value of the noncentral tdistribution by numeric integration... Maybe somebody of you has access to SAS or software specialized in power analysis (e.g., PASS or nQuery Advisor) and would like to check these results? — Diftor heh smusma 🖖🏼 Довге життя Україна! _{} Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
Helmut ★★★ Vienna, Austria, 20070106 16:17 (6405 d 18:33 ago) @ Helmut Posting: # 441 Views: 43,142 

Dear all, Dieter Hauschke informed me that Jones’/Kenward’s code uses Satterthwaite’s approximation for degrees of freedom (for heterogenous variances), which was corrected in a recent book. ^{[1]} In power calculations degrees of freedom are fixed by the number of subjects (and sequences). If you want to play around again, please replace code line #22
Power (CV=20%) by the (old) new codes for n=24 are:
+======+=======+=======+ Power >80% (CV=20%) by the (old) new codes are:
+======+=====+===========+===========+ Power is >80% for all combinations tested. ❝ Maybe somebody of you has access to SAS or software specialized in power analysis (e.g., PASS or NQuery) and would like to check these results? According to Dieter Hauschke Owen’s exact method is implemented in nQuery Advisor. Since there is still some confusion about different methods for sample size estimation in bioequivalence (software, tables, approximations) an entire chapter of a new book ^{[2]} will be devoted to this issue.
— Diftor heh smusma 🖖🏼 Довге життя Україна! _{} Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
mathews ● 20080403 15:19 (5952 d 20:31 ago) @ Helmut Posting: # 1757 Views: 41,904 

Dear all, Please tell me how to calculate power for a replicate crossover design (2 treatment 2 sequence 4 period). Regards Matz 
Helmut ★★★ Vienna, Austria, 20080403 15:49 (5952 d 20:01 ago) @ mathews Posting: # 1758 Views: 42,918 

Dear Matz, you will get the same power in ~½N subjects (TRRTRTTR) as compared to a 2×2×2 design in N subjects (half the sample size  but same number of treatments!). If you perform a study in a 2treatment 2sequence 3period design (e.g., TRTRTR) the same power is obtained in ~¾N subjects. ^{1)}
— Diftor heh smusma 🖖🏼 Довге життя Україна! _{} Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
mathews ● 20080404 15:21 (5951 d 20:28 ago) @ Helmut Posting: # 1763 Views: 41,391 

Dear HS, Thanks for your reply. How we can modify the above R code for a replicate crossover design. The above code used within subject s.d for a 2x2 design, but in a replicate crossover design we will get within subject s.d for both test and reference. So which value should be used?... Regards, Matz 
Helmut ★★★ Vienna, Austria, 20090925 21:20 (5412 d 14:29 ago) @ Helmut Posting: # 4255 Views: 40,469 

Dear all, I made an investment and was able to compare results from R with the current version 7 of nQuery Advisor:*
+======+=====+===========+========+==========+ nQuery requires sqrt(MSE) rather than CV as input = sqrt(ln(CV²+1)). The maximum precision is 6 decimal places. Therefore 0.1980422 rather than the value in full precision (0.1980422004353650284627675024887) like in R is used. But you see this value only if clicking the respective cell; both displayed and printed is 0.198042. nQuery gives power in percent with a maximum number of two decimal digits. That's a pity and values in the first column were the best one could get (display/print). If clicking in each cell power is given in percent to 5 decimal digits (7 significant digits); second column. It's clear that results are not rounded, but truncated. This looked weird first, but at least for power it makes sense. Values agree nicely with the exception of a small discrepancy at GMR=1 (83.31980% vs. 83.32001%). If I recall it correctly Dieter Hauschke wrote something about numerical problems in power estimation at T/R=1 in one of his papers. I love validating software.
— Diftor heh smusma 🖖🏼 Довге життя Україна! _{} Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
ElMaestro ★★★ Denmark, 20090928 01:27 (5410 d 10:23 ago) @ Helmut Posting: # 4256 Views: 40,037 

Hi HS, ❝ [Blahdeeblah] Therefore 0.1980422 rather than the value in full precision (0.1980422004353650284627675024887) [more blah] Try to put yourself in the shoes of the software writer. There are two known options to getting the calculation of power. One is the bruteforce approach, the other goes via probabilities of the noncentral t distribution. The former can be cumbersome and is generally not preferred (takes time but converges asymptotically to the true power with the number of resamples if the random number generator works (and trust me, it is not that hard to make one that does!)), but there is still preference for the latter. Now, actually programming a stable Al Gore Rhythm which get the probs from the nct dist is a nightmare, so the developer of course googles it and finds some really neat source code. This is what the developers of R did, and a dollar gets ten that this is what others do as well. And here's the deal: If you look at that source you can easily see that it is based on some (necessary) helper constants and shortcuts and convergence criteria and what not. In other words, whatever the algo does it is not exact. It is approximate, but how good that approximation actually is is unknown (unknown error will be a function of the conditions, by the way). We know R's power is an approximation since it uses an adapted ASA243. So when two or more applications get the same result it could be either because they are both right or because they both make the same errors (as would be the case if they are based on mods of the same algo etc), or it could be pure chance. And we have absolutely no way of telling which one is the right one until some clever guy works out the integrals exactly. I'll buy you a Mozartkugel if you do it. ❝ I love validating software. Decompile your recent acquisition and check if they one way or another use the ASA243. If they do, why validate? EM. 
Helmut ★★★ Vienna, Austria, 20090928 17:21 (5409 d 18:28 ago) @ ElMaestro Posting: # 4260 Views: 41,009 

Hi ElMaestro! ❝ Try to put yourself in the shoes of the software writer. Oh, I took these boots on many years ago. And they got muddy almost immediatelly. I tried to get a numerical approximation of the cdf of the tdistribution which lead me to the gammafunction and headaches. Finally I ended up with numerical approximations.^{[1,2]} ❝ […] whatever the algo does it is not exact. It is approximate, but how good that approximation actually is is unknown. Sure. ❝ […] we have absolutely no way of telling which one is the right one until some clever guy works out the integrals exactly. Definitely not me. I’m getting more stupid every day. ❝ I'll buy you a Mozartkugel if you do it. THX, no needs. ❝ Decompile your recent acquisition and check if they one way or another use the ASA243. I guess the SW was built using the Borland C++ compiler. You know that reverse engineering of software is a breach of the license? BTW, according to the manual the ‘Illinois method’ (Algorithm AS 184, FORTRAN source) is used. I assume AS 243^{[3]} is ‘better’ than AS 184.^{[4]} But to which extent? Is it my job to go through pros and cons of algorithms? ❝ If they do, why validate? First, different algos are used. Second, I have to validate software  or not? At the product's site I found a nice statement: These sets of solutions were reviewed by Janet Elashoff, who checked for consistency, face validity, and for computational accuracy against other sources. Face validity?What I actually wanted to do was to check my R code (which is also part of bear) with a piece of SW of high reputation. Checking the sample size I wasn’t satisfied (only integers), therefore I opted to do it the other way 'round (power). It made me angry that there is no (easy) way to obtain more than 4 significant digits. By this external validation is effectively prevented. Another example from the dirt track of SW validation: Diletti et al.^{[5]} published exact samples size tables for BE ranges of 0.911 and 0.7–1.43. When trying to validate my code against these tables I found small discrepancies (only at low CVs and close to the acceptance range). By trialanderror I got the solution: I could reproduce tables only if setting the upper acceptance limit not to the reciprocal of the lower AL (0.9^{1}, 0.7^{1}) but to exactly 1.1111 or 1.4286… BTW, nitpicking Diletti’s Table 1^{[6]} (power 70%, T/R=1, CV=7.5%) and nQuery 7 give n=4, whilst my R code, StudySize 2.0.1, and FARTSSIE 1.6 come up with n=6. Now what? Go with a democratic vote of 3:2 for 6?
— Diftor heh smusma 🖖🏼 Довге життя Україна! _{} Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
ElMaestro ★★★ Denmark, 20090928 21:12 (5409 d 14:38 ago) @ Helmut Posting: # 4261 Views: 40,132 

Dear HS, ❝ Oh, I took these boots on many years ago. And they got muddy almost immediatelly. I tried to get a numerical approximation of the cdf of the tdistribution which lead me to the gammafunction and headaches. ❝ Finally I ended up with numerical approximations.^{[1,2]} Let me highfive you for this. I did the same, in C, and it took me quite some time, only to find out later that C actually comes with a beautiful function called lgamma which does it and returns the logarithm. I wasted a good part of my youth reinventing the wheel it seems ❝ BTW, nitpicking Diletti’s Table 1^{[6]} (power 70%, T/R=1, CV=7.5%) and nQuery 7 give n=4, whilst my R code, StudySize 2.0.1, and FARTSSIE 1.6 come up with n=6. Now what? Go with a democratic vote of 3:2 for 6? Two solutions:

Helmut ★★★ Vienna, Austria, 20090928 21:52 (5409 d 13:58 ago) @ ElMaestro Posting: # 4263 Views: 39,839 

Hi ElMaestro! ❝ Let me highfive you for this. BTW, great subject line! 42? (or look what _{} has to say). ❝ I wasted a good part of my youth reinventing the wheel Welcome to the club of wheel reinventors. ❝ ❝ BTW, nitpicking Diletti's Table 1 […] ❝ Two solutions: ❝ 1. Ask dlabes Ha! Have you read his post? With that code you can reproduce the block for power=80% […] with some minor deviations for CV=5% and 7.5%. ❝ 2. Get your power from a program that bases its evaluation on brute force rather than on numerical integration. Do you have one handy? — Diftor heh smusma 🖖🏼 Довге життя Україна! _{} Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
ElMaestro ★★★ Denmark, 20090928 22:18 (5409 d 13:31 ago) @ Helmut Posting: # 4264 Views: 40,090 

Hi HS, ❝ ❝ 2. Get your power from a program that bases its evaluation on brute force rather than on numerical integration. ❝ ❝ Do you have one handy? I think I can make you an elmaestrolophystic piece of software for this; will take "some time"^{*} to write it, I guess. EM. *: Remember Alan Turing? He proved that there is no way of knowing if an algorithm converges until it converges. If it converges. Unfortunately it works the same way when I fire up my C compiler. In spite of good intentions and fancy plans I may get personally stuck in an infinite loop. By the way, it is pretty much the same with women and telephones, but that's a completely different story. And another forum, I guess. 
yjlee168 ★★★ Kaohsiung, Taiwan, 20090928 23:30 (5409 d 12:19 ago) @ Helmut Posting: # 4265 Views: 39,980 

Dear Helmut, I just tried PASS and got the following (for power 70%, T/R=1, CV=7.5%). Power Analysis of 2x2 CrossOver Design for Testing Equivalence Using Ratios N = 4. Sorry for my previous message (I set power= 0.8, not 0.7!). ❝ [...] ❝ Do you have one handy? — All the best,  Yungjin Lee bear v2.9.1: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan https://www.pkpd168.com/bear Download link (updated) > here 
Helmut ★★★ Vienna, Austria, 20090929 01:18 (5409 d 10:32 ago) @ yjlee168 Posting: # 4266 Views: 40,219 

Good morning, bears! Thanks for the PASSrun. Now the game is tied at 3:3 for sample sizes of 4 or 6 (target power 70%, T/R=1, CV=7.5%). power % n Underlined sample sizes are rounded up to the next even integer (calculated results in parentheses). From the manual I would say that PASS uses also AS 243. — Diftor heh smusma 🖖🏼 Довге життя Україна! _{} Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
d_labes ★★★ Berlin, Germany, 20090929 15:37 (5408 d 20:13 ago) @ Helmut Posting: # 4267 Views: 40,083 

Dear guys, let me take part in your hairsplitting party . Hope you have already waited for me. Here are the results of the SAS Proc Power. Code can be found here, of course appropriately adapted. The POWER Procedure The power is given with 3 decimals only by default. No simple way to change. Seems the power is calculated in this case totally wrong! This is astonishing enough, because the same Proc Power code can reproduce most of the sample sizes given in Diletti et.al., at least for power 80%. Ok, Proc Power is called experimental by SAS, usually means the user shall do the testing and report bugs. Others call such software beta. My own rolled code (relying on Owen's Qfunctions, undocumented but available in SAS) gives:
Coeff. of — Regards, Detlew 
Helmut ★★★ Vienna, Austria, 20090929 16:52 (5408 d 18:58 ago) @ d_labes Posting: # 4269 Views: 40,402 

Dear D Labes! ❝ Hope you have already waited for me. Yesss! ❝ Here are the results of the SAS Proc Power. ❝ Seems the power is calculated in this case totally wrong! Fascinating. ❝ My own rolled code (relying on Owens Qfunctions, undocumented but available in SAS) gives: That’s the more interesting part (homebrew)! n % power n % power n=4 takes the lead by 4:3! Observations: #5 (based on Algorithm AS 243) and #6 give identical results. This is not surprising, because the VBAroutine within FARTSSIE was written by Russell V. Lenth, the author of AS 243. #2/3 agree; for n=6 power is comparable with #5/6. #4 calculates lower power than #13, but still would suggest n=4. The difference between #13 and #5/6 might be due to different algorithms for estimating the noncentral tdistribution (#1/3: ?, #2: AS 243?, #4: AS 184, #5/6: AS 243). #7 uses approximations (Ref. 2 in this post); however, yields conservative results. Maybe I should contact Dieter Hauschke to join the party… Edit: Many, many years later.
— Diftor heh smusma 🖖🏼 Довге життя Україна! _{} Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
ElMaestro ★★★ Denmark, 20090929 20:14 (5408 d 15:36 ago) @ Helmut Posting: # 4270 Views: 39,866 

Dear HS, ❝ n=4 takes the lead by 4:3! Eighth observation: I made a little application that goes under the codename EFG (which is short for ElF#ckupGrande) and which implements brute force in conjunction with Mersenne Twister and some very ugly elmaestrolophystic code clutter. At N=4, CV=0.075 and T/R=1.0 I get power = 72.7191% with 2000000 iterations. At N=6 power is 98.97435%. EM. 
yjlee168 ★★★ Kaohsiung, Taiwan, 20090930 01:55 (5408 d 09:55 ago) @ Helmut Posting: # 4272 Views: 39,910 

Dear Helmut, Please read my message of this thread. Sorry for my mistake. I find that PASS can reach the results that are pretty close to D Labes's SAS run. ❝ That's the more interesting part (homebrew)! n % power n % power — All the best,  Yungjin Lee bear v2.9.1: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan https://www.pkpd168.com/bear Download link (updated) > here 
yjlee168 ★★★ Kaohsiung, Taiwan, 20090930 01:49 (5408 d 10:01 ago) @ yjlee168 Posting: # 4271 Views: 39,745 

Dear Helmut, I made some correction with PASSrun for the power when N = 6. Please help to update your collection data. Thanks. Power Analysis of 2x2 CrossOver Design for Testing Equivalence — All the best,  Yungjin Lee bear v2.9.1: created by Hsinya Lee & Yungjin Lee Kaohsiung, Taiwan https://www.pkpd168.com/bear Download link (updated) > here 
KDA ☆ India, 20190123 05:24 (2006 d 05:26 ago) @ Helmut Posting: # 19800 Views: 16,388 

Dear all, I am a amateur in BA/BE area and trying to learn from this forum. I would appreciate if anyone kindly answer few of my queries for the rcode used to reproduce Fig.1c by Diletti et al. (1991). Q1) If it is 2x2 crossover study, test and ref drug would be given only once. Can we get true intrasubject coefficient of variance (CV)from 2x2 crossover study? Q2) If it is possible to calculate the intrasubject coefficient of variance (CV)from 2x2 crossover study, is this value coming from ref drug only? Q3) What is the equation for calculating this intrasubject coefficient of variance (CV)? Q4) In the equation "s = sqrt(2)*sigmaW", what 's' stands for? why '2' is squarerooted? Q5) I have used the following rcode for calculating point estimate and CI limit of AUC. A1 < data$AUC[which(data$TRT=="R1" & data$PER==1)]
Can I calculate CV by the following equation? CV=SD/(exp(mD)) Based on my understanding, this my calculated CV is not the same as the 'CV' that represent intrasubject coefficient of variance used to reproduce Fig.1c by Diletti et al. (1991) what code/formula should I use if I want to calculate intrasubject coefficient of variance? Regards KDA 
Ohlbe ★★★ France, 20190127 19:24 (2001 d 15:26 ago) @ KDA Posting: # 19814 Views: 15,997 

Dear KDA, ❝ Q1) If it is 2x2 crossover study, test and ref drug would be given only once. Can we get true intrasubject coefficient of variance (CV)from 2x2 crossover study? No. You need a replicate design study. I'll let others answer Q5. — Regards Ohlbe 
Brus ★ Spain, 20211202 16:00 (961 d 18:50 ago) @ Helmut Posting: # 22683 Views: 6,206 

Dear Helmut, I am not an expert in R, and I am trying to adapt your example of power calculation with plot replesentation, but for the calculation of the sample size. My objetive is to plot a graphicala representation with ratio in X axis and sample size in Y axis and setting the CV and power as a fixed value. Can you or someone help me? Best regards, 
ElMaestro ★★★ Denmark, 20211202 16:15 (961 d 18:34 ago) @ Brus Posting: # 22684 Views: 6,324 

Hello Brus, ❝ I am not an expert in R, and I am trying to adapt your example of power calculation with plot replesentation, but for the calculation of the sample size. My objetive is to plot a graphicala representation with ratio in X axis and sample size in Y axis and setting the CV and power as a fixed value. Here's something basic to get you started. I assume you already have the PowerTOST package on your system. Try.this=function(CV=0.2, tpow=0.8) Can be refined with designs and multipanel options and what not. PS: Correct, I am not entirely comfy with the tapply, sapply, apply family of functions.— Pass or fail! ElMaestro 
Helmut ★★★ Vienna, Austria, 20211202 18:30 (961 d 16:20 ago) @ ElMaestro Posting: # 22686 Views: 6,168 

Hi ElMaestro, ❝ A suggestion. ❝ PS: Correct, I am not entirely comfy with the You are not alone. To complicate things, in PowerTOST many functions are not or only partly vectorized. See the script in this article for a remedy (click on to show it).— Diftor heh smusma 🖖🏼 Довге життя Україна! _{} Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 
Helmut ★★★ Vienna, Austria, 20211202 17:57 (961 d 16:52 ago) @ Brus Posting: # 22685 Views: 6,201 

Hi Brus, ❝ I am trying to adapt your example of power calculation with plot replesentation, but for the calculation of the sample size. My objetive is to plot a graphicala representation with ratio in X axis and sample size in Y axis and setting the CV and power as a fixed value. I suggest this goody first. May sound picky but we can’t get the exact desired power, only one which is at least that high.* Try this one:
The sample size is a staircase function and due to the symmetry in \(\small{\log_{e}}\)scale we require the same sample size for any given \(\small{\theta_0}\) and \(1/\small{\theta_0}\) (see also there).
Diving deeper into the matter: Power is a hypersurface depending on \(\small{\theta_0}\), \(\small{n}\), and \(\small{CV}\). We can only show spatial projections at certain \(\small{CV}\text{}\)values (the panels) with \(\small{\log_{e}\theta_0}\) (xaxis), \(\small{n}\) (yaxis), and power (zaxis). \(\small{n=12292}\) and \(\small{\theta_0=0.801.25}\) to demonstrate that at the limits of the acceptance range power (which is there the Type I Error) is ≤0.05. The intersections with the horizontal plane (the target power) show the Ushaped dependency of \(\small{n}\) on \(\small{\theta_0}\) as in the plot above. The slices at certain \(\small{\theta_0}\) and \(\small{n}\)values give the plots like in the very old post above. Don’t try to imagine fourdimensional objects. Even a simple hypercube (tesseract) may fry you brain. — Diftor heh smusma 🖖🏼 Довге життя Україна! _{} Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮 Science Quotes 