ElMaestro
★★★

Denmark,
2009-10-03 21:05
(5684 d 04:14 ago)

Posting: # 4291
Views: 35,827
 

 non-central t, power, R: [🇷 for BE/BA]

Post reply
Dear all,

In a previous thread I asked what the power in a 2,2,2-BE trial is at N=6, CV=65% and T/R=95%.
The things is, Fartssie comes up with a negative reply. HS's implementation in R of the equations given by Julious (Statistics in Medicine, 2004 vol. 23, 1921-1986) also evaluates to a negative figure.
A central role here is played by the non-central t distribution. My guess was that we were torturing the algorithm, and I think this was HS's point as well (quote "I think we are pushing software to the limits, i.e., running into troubles of getting a reasonable value of the noncentral t-distribution (...)").

This I wanted to look a little into using R as my reference. Wikipedia tells us what the distribution's density function looks like. I made an independent implementation of this expression, plus an implementation of the integrated density function from t and downwards, which is basically the same as the function pt in R which is based on ASA243 according to the documentation.

The values from pt can be easily reproduced, it turns out. But I cannot "improve" them, although I have all the time in the world to integrate with ridiculously small stepsizes and tweak and fiddle. As it turns out in the concrete example, we are asking for probabilities from the noncentral t dist in a region of t-values plusminus 2.13 and noncentrality parameters 0.50 and -0.80. These are not at all extreme values. I think it sounds reasonable that "prob1" and "prob2" (using HS's terminology) are around 0.89 and 0.17 [try and plot the densities as function of t, eyeball the relative amount of curve left of the indicated values].

For this reason I am becoming a sceptic now. Either I am not able to make a better implementation than ASA243, or ...please forgive me... the equations given by Julious (and others) are not exact in certain situations (for example: I only see the issue when t2 < nc2), or something entirely else is in play.

Can anyone shed light on this?

Let me add: If the problem lies in ASA243 then I'd love to write one that gets it right, of course.)

Many thanks for any input and best regards,
EM.
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2009-10-03 23:45
(5684 d 01:34 ago)

@ ElMaestro
Posting: # 4292
Views: 32,561
 

 non-central t, power, R:

Post reply
Hi ElMaestro,

very interesting question! I'm too short in time to read Owen's paper again (and likely too stupid to understand it as well). nQuery Advisor 7 (ASA184, I guess) comes up with
The computed power was too low for the approximation to be considered reliable. Use a larger sample size.
The lowest sample size without an error message was 40 (power 5.871678%, more than 10times the value R comes up with)...

Have you considered contacting Russell Lenth (author of ASA243 and the VBA-routines in FARTSSIE)? He writes:

Accuracy

  • Formula accuracy
    Most computations are "exact" in the sense that they are based on exact formulas for sample size, power, etc.
  • Machine accuracy
    Even with exact formulas, computed values are inexact, as are all double-precision floating-point computations. Many computations (especially noncentral distributions) require summing one or more series, and there is a serious tradeoff between speed and accuracy. The error bound set for cdfs is 1E-8 or smaller, and for quantiles the bound is 1E-6. Actual errors can be much larger due to accumulated errors or other reasons. Quantiles, for example, are computed by numerically solving an equation involving the cdf; thus, in extreme cases, a small error in the cdf can create a large error in the quantile.
    […] in the case of quantile computations, no warning message is generated for extreme quantiles. If you want a power of .9999 at alpha=.0001, you can expect the computed sample size to not be accurate to the nearest integer! If you specify reasonable criteria, the answers will be pretty reliable.

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
ElMaestro
★★★

Denmark,
2009-10-04 23:51
(5683 d 01:29 ago)

@ Helmut
Posting: # 4293
Views: 32,370
 

 non-central t, power, R:

Post reply
Dear HS,

❝ Have you considered contacting Russell Lenth (author of ASA243 and the VBA-routines in FARTSSIE)?


Might be the way to go. One could hope that the author of Fartssie reads these discussions and looks into it. It is of course true, depending on how we look at this there is a lot of numerical integration and/or summation going on, and a lot of 'small' errors would indeed sooner or later amount to a 'large' error. To be honest I am just a little sceptical here, because we are not really involving extreme values of t or the noncentality parameters, so I am not at this point fully convinced that this is no more than just a numerical precision issue.
Must read some more papers, I guess. Or implement some C library which works with 117 decimals or something obscene and try anew.

If anyone has any insight into how Al Gore Rhythms can be validated when true results at the test conditions are not known then I'd be happy to learn it. Or even better, how to determine under which conditions ASA243 is accurate to 3 decimals or something like that.

EM.
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2009-10-05 01:58
(5682 d 23:22 ago)

@ ElMaestro
Posting: # 4294
Views: 32,380
 

 non-central t, power, R:

Post reply
Hi ElMaestro!

❝ Might be the way to go. One could hope that the author of Fartssie reads these discussions and looks into it.


Hhm, David logged in the last time in mid-March. But he has activated his e-mail in the forum and I had some nice conversations in the past.

❝ […] To be honest I am just a little sceptical here, because we are not really involving extreme values of t or the noncentality parameters, so I am not at this point fully convinced that this is no more than just a numerical precision issue.

❝ Must read some more papers, I guess. Or implement some C library which works with 117 decimals or something obscene and try anew.


Sorry I'm of no big help now (still a lot of slides to prepare for a two-day workshop next week), but I think the algo uses a Taylor-expansion which is truncated after the first (?) term. One thing coming into my mind is the possibility that when we run into numerical issues, they may be resolved by going 'deeper' into the Taylor series?

❝ If anyone has any insight [...]



Sancta simplicitas, non!
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
d_labes
★★★

Berlin, Germany,
2009-10-05 14:58
(5682 d 10:22 ago)

@ ElMaestro
Posting: # 4296
Views: 32,755
 

 non-central t is not the question

Post reply
Dear ElMaestro,

don't waste your time in attempting to improve Al Gore's rhythms.
The inventors of them for the non-central t distri had already made every effort for its accuracy and reliability, I think.

Your time is better invested for NLYW (blonde downstairs) or for the Seven seas :-D .

The formula given in Julious is definitely approximate, albeit a good one for the usual power range ≥70% as the sample size calculations with that formula show.

The exact algo for the power of the equivalence test (two one sided tests aka TOST) in a 2x2 cross-over is ([1],[2]):

Power=Prob(t1≥t1-alpha,(n-2) and t2≤-t1-alpha,(n-2)|bioequivalence)


The t1 and t2 values are the t-variates of the two one-sided t-tests of bioequivalence.

(t1, t2) have according to Owen a bivariate non-central t-distribution and the power can be calculated as the difference of two definite integrals (Owen's Q functions dependent on df=degrees of freedom and with four arguments)

Power=Qdf(-t1-alpha,(n-2),d2,0,R)-Qdf(t1-alpha,(n-2),d1,0,R)


See [1] and [2] for that and for the formulas for d1, d2 and R.

So you need not a better algo for the non-central t but rather an algo for Owen's Q-functions :-P .
My SAS beast has these functions implemented (but not documented!).

❝ In a previous thread I asked what the power in a 2,2,2-BE trial is at N=6, CV=65% and T/R=95%.


And my sasophylistic answer to the your ElMaestrophylistic question "what the power in a 2,2,2-BE trial is at N=38, CV=65% and T/R=95%." given here, power=4.5%, is based on these Q-functions.
Implementing the Julious formulas (according to Helmut's R-code) in "The power to know" results in power=-0.027799235, same as Fartssie!

[1] Diletti, E., Hauschke, D., and Steinijans, V.W. (1991),
"Sample Size Determination for Bioequivalence Assessment by Means of Confidence Intervals"
International Journal of Clinical Pharmacology, Therapy and Toxicology, Vol. 29, 1 -8.
[2] Hauschke, Steinijans, Pigeot
"Bioequivalence studies in drug development"
Wiley 2007, pages 120-121

—
Regards,

Detlew
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2009-10-05 16:53
(5682 d 08:26 ago)

@ d_labes
Posting: # 4297
Views: 32,508
 

 SAS -> R -> SAS

Post reply
Dear D Labes,

❝ Implementing the Julious formulas (according to Helmut's R-code) in "The power to know" results in power=-0.027799235, same as Fartssie!


Funny!
My R-code was just ported from Patterson's/Jones' SAS-code. Tweaking the code to calculate power (haha) at n=38:

options(digits=8)
alpha  <- 0.05; CV <- 0.65; ratio <- 0.95; n <- 38
Theta1 <- 0.80; Theta2 <- 1/Theta1
s      <- sqrt(2)*sqrt(log(1+CV^2))
t1     <- qt(1-alpha,n-2)
nc1    <- sqrt(n)*(log(ratio)-log(Theta1))/s
nc2    <- sqrt(n)*(log(ratio)-log(Theta2))/s
prob1  <- pt(+t1,n-2,nc1)
prob2  <- pt(-t1,n-2,nc2)
power  <- prob2-prob1
power

I get -0.027799235, which is either an Ouroboros or another strong hint that SAS and R follow the same algorithm.
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
d_labes
★★★

Berlin, Germany,
2009-10-06 16:31
(5681 d 08:49 ago)

@ Helmut
Posting: # 4301
Views: 39,487
 

 Lost in translation

Post reply
Dear Helmut,

❝ [...] which is either an Ouroboros or [...]



It is always a pleasure for me to learn from you. Thanks for that nice dragon I'm not aware of until now. :cool:

Seems to me as an allegory of some/many things I/we do in my/our daily work.
(delete as applicable)

Regarding the translation SAS->R->SAS I must notice that we are in great luck to have the same results after that cycle.
Others are not if I think about f.i. english->german->english translations.
Although such back and forth translations are seen as QA measure in doing translations of study documents, f.i. volunteers informed consent which must be in the language of the country a study is performed in.

To try out what happens with automatically translations visit Lost in translation. :-D
—
Regards,

Detlew
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2009-10-06 16:42
(5681 d 08:37 ago)

@ d_labes
Posting: # 4302
Views: 32,370
 

 Lost in translation

Post reply
Dear D Labes!

❝ ❝ […] an Ouroboros

❝

❝ Thanks for that nice dragon I'm not aware of until now. :cool:


As a trained chemist I heard of Kekulé's daydream which lead him to the discovery of the unsaturated ring structure of benzene. That's why I know this beast.

❝ Seems to me as an allegory of some/many things I/we do in my/our daily work.

❝ (delete as applicable)

Nothing to delete. :sleeping:
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
d_labes
★★★

Berlin, Germany,
2009-10-06 17:14
(5681 d 08:05 ago)

@ Helmut
Posting: # 4303
Views: 32,255
 

 Macacus cynocephalus

Post reply
Dear Helmut!

❝ As a trained chemist I heard of Kekulé's daydream


Here comes the true Kekulé's daydream:
[image]
Source:
Berichte der durstigen Chemischen Gesellschaft
Unerhörter Jahrgang, Nr. 20 zum 20.9.1886

Translated by babelfish:
Reports of the thirsty chemical society
outrageous class, No. 20 to 20.9.1886

Back to German
Reports der unverschämten Kategorie der durstigen chemischen Gesellschaft,
Nr. 20 bis 20.9.1886

To English again:
Report of the impudent category of the thirsty chemical society,
No. 20 to 20.9.1886 ....
—
Regards,

Detlew
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2009-10-06 23:46
(5681 d 01:34 ago)

@ d_labes
Posting: # 4306
Views: 32,289
 

 Lost in translation

Post reply
Dear D Labes!

[image] translate:
  1. DE → EN
    Reports of the thirsty Chemical Society
    Unheard-Year, No. 20 at 20.9.1886
  2. EN → DE
    Berichte der durstigen Chemical Society
    Unerhört Jahrgang, Nr. 20 am 20.9.1886
  3. DE → EN
    Reports of the thirsty Chemical Society
    Unheard-Volume, No. 20 on 20.9.1886
  4. EN → DE
    Berichte der durstigen Chemical Society
    Unerhört-Band, Nr. 20 am 20.9.1886
  5. DE → EN
    Reports of the thirsty Chemical Society
    Unheard-Volume, No 20 on 20.9.1886
What looked like a strange attractor first, stabilised after a few iterations and reached convergence in step 5. :-D
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
ElMaestro
★★★

Denmark,
2009-10-05 19:32
(5682 d 05:47 ago)

@ d_labes
Posting: # 4298
Views: 32,696
 

 non-central t is not the question

Post reply
Dear dlabes,

❝ Your time is better invested for NLYW (blonde downstairs) or for the Seven seas :-D .


I agree. One of the options you mention would get higher priority then the other if I had a choice.

❝ The formula given in Julious is definitely approximate (...)


Cool. It becomes more and more difficult to find out who to trust.

❝ So you need not a better algo for the non-central t but rather an algo for Owen's Q-functions :-P .


Excellent input. The non-central t dist is hereby dead and buried.
Seems this is a good place to start, and this power document gives the Qv function. Must see if I can implement that.

Thanks for the good post!

EM.
 
d_labes
★★★

Berlin, Germany,
2009-10-06 18:45
(5681 d 06:35 ago)

@ ElMaestro
Posting: # 4304
Views: 32,514
 

 Power of TOST in R

Post reply
Dear ElMaestro!

❝ Excellent input. The non-central t dist is hereby dead and buried. Seems this is a good place to start, and this power document gives the Qv function. Must see if I can implement that.


Once again: Your time is better invested for whatever you prefer :-D .
Also I confess that implementing it itself can be very satisfactory!
('Oh what a great guy I am which can solve such complicated definite integral!')

Have a look at the R package MBESS.
If I had seen it right there is a function power.equivalence.md that will give you some spare time.
Could not test it by myself. No spare time at moment.
—
Regards,

Detlew
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2009-10-07 04:18
(5680 d 21:02 ago)

@ d_labes
Posting: # 4309
Views: 32,342
 

 Power of TOST (R package MBESS)

Post reply
Dear D Labes!

❝ Have a look at the R package MBESS. If I had seen it right there is a function power.equivalence.md that will give you some spare time.


Quick tests

require(MBESS)
n  <- c(8,12,18,24,30,40,60) # total sample size
df <- n-2
power.equivalence.md.plot(.05, logscale=TRUE, .8, 1.25, .20,
                          n, df, 'Diletti, Figure 1c')

[image]
In the text output power is the theoretical value 0.05 (=α) at the borders of the acceptance range.

Now for ElMaestro's candidates:
n <- 6
power.equivalence.md(.05, logscale=TRUE, 0.8, 1.25, 0.95, 0.65, n, n-2)
[1] 0.001048352
n <- 38
power.equivalence.md(.05, logscale=TRUE, 0.8, 1.25, 0.95, 0.65, n, n-2)
[1] 0.01392505
Yet some other values for our collection. ;-)
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
d_labes
★★★

Berlin, Germany,
2009-10-07 10:46
(5680 d 14:34 ago)

@ Helmut
Posting: # 4310
Views: 32,356
 

 Power of TOST (R package MBESS)

Post reply
Dear Helmut!

❝ Quick tests


Oh actually quick! But dirty :-P .
I assume that your 5th parameter is meant as CV.
But the function requires sigma, I think.

With that modification

library(MBESS)
CV <- 0.65
sigma<-sqrt(log(1+(CV)**2))

n <- 6
power.equivalence.md(0.05,TRUE, 0.8, 1.25, 0.95, sigma, n, n-2)
[1] 0.001620383

n <- 38
power.equivalence.md(0.05,TRUE, 0.8, 1.25, 0.95, sigma, n, n-2)
[1] 0.04509401

The answer for n=38 seems to me very sasophylistic home-brewed magic and Q-functionalistic :-D .
Not so astonishing if you look at the author of that pice of R-software.

PS: I have just seen that the documentation suggest that CV can be used in case of log. But this seems erroneous.
—
Regards,

Detlew
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2009-10-07 14:38
(5680 d 10:41 ago)

@ d_labes
Posting: # 4314
Views: 32,213
 

 Power of TOST (R package MBESS)

Post reply
Dear D Labes!

❝ Oh actually quick! But dirty :-P .


I always thought of it as a nice combination. ;-)

❝ The answer for n=38 seems to me very sasophylistic home-brewed magic and Q-functionalistic :-D .

❝ Not so astonishing if you look at the author of that pice of R-software.


Ah, Kem Phillips!

❝ PS: I have just seen that the documentation suggest that CV can be used in case of log. But this seems erroneous.


Exactly. The example in the documentation drove me to my quickshot.
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
d_labes
★★★

Berlin, Germany,
2009-10-07 13:32
(5680 d 11:48 ago)

@ d_labes
Posting: # 4312
Views: 32,426
 

 The Power at limits

Post reply
Dear All!

Back from the date with the Blonde downstairs an having some spare time ;-),

Seems ElMaestro's question is able to put the power to the limits.
Here results of power.equivalence.md around N=6:

R code:
library(MBESS)

alpha    <- 0.05
logscale <- TRUE
ltheta1  <- 0.8
ltheta2  <- 1/ltheta1
ldiff    <- 0.95
CV       <- 0.65                # added [HS]
sigma    <- sqrt(log(1+(CV)^2)) # added [HS]
n        <- c(4,6,8,10,12,14,16)
df       <- n-2
pow      <- mapply(power.equivalence.md,n=n,nu=df, MoreArgs =list(alpha, logscale, ltheta1, ltheta2, ldiff, sigma))
res      <- data.frame(n=n, power=pow)
res

Result:
   n        power
1  4 0.0045417784
2  6 0.0016203829
3  8 0.0009763342
4 10 0.0007993335
5 12 0.0007966979
6 14 0.0009067329
7 16 0.0011306148

I would expect power is increasing with increasing n.
Or miss I here something?

BTW: My SAS code shows the same effect.
All this is hairsplitting of course.

Edit: Added two lines above to allow copy-pasting of code to R-console. BTW, change the sample size vector to n<-c(4,6,8,10,12,14,16,38,40,152,154,344,346) and watch (the gamma-function needs the faculty, which runs out of steam - at least on 32-bit operating systems). [Helmut]
—
Regards,

Detlew
 
d_labes
★★★

Berlin, Germany,
2009-10-07 17:15
(5680 d 08:05 ago)

@ d_labes
Posting: # 4316
Views: 32,096
 

 Power at limits 2

Post reply
Dear All!

Must post to myself because Helmut has posted an important finding in the edit line of the post of mine.

❝

❝ [...] BTW, change the sample size vector to n<-c(4,6,8,10,12,14,16,38,40,152,154,344,346) and watch (the gamma-function needs the faculty, which runs out of steam - at least on 32-bit operating systems). [Helmut]


Seems there is an implementation flaw in the code. And to add insult to injury it runs without error. But power is zero for n=346! :surprised:

After inspecting the code at the R-console
> power.density.equivalence.md
function (power_sigma, alpha = alpha, theta1 = theta1, theta2 = theta2,
    diff = diff, sigma = sigma, n = n, nu = nu)
{
    pdim <- length(power_sigma)
    power_density <- c(1:pdim)
    power_t <- qt(1 - alpha, nu, lower.tail = TRUE, log.p = FALSE)
    a <- sigma/sqrt(nu)
    power_const <- -(nu/2 - 1) * log(2) - log(gamma(nu/2)) - log(a)
    for (isigma in 1:pdim) {
        power_d1 <- ((power_sigma[isigma] * power_t * sqrt(2/n)) +
            (theta1 - diff))/(sigma * sqrt(2/n))
        power_d2 <- ((-power_sigma[isigma] * power_t * sqrt(2/n)) +
            (theta2 - diff))/(sigma * sqrt(2/n))
        power_phi1 <- pnorm(power_d1, mean = 0, sd = 1, lower.tail = TRUE,  log.p = FALSE)
        power_phi2 <- pnorm(power_d2, mean = 0, sd = 1, lower.tail = TRUE,  log.p = FALSE)
        power_phi <- power_phi2 - power_phi1
        power_density[isigma] <- power_phi * exp(power_const +
            (nu - 1) * log(power_sigma[isigma]/(a)) - 0.5 * (power_sigma[isigma]/(a))^2)
    }
    return(power_density)
}
you can see where it is.
So roll your own power.density.equivalence.md and change the line with red marker to
    power_const <- -(nu/2 - 1) * log(2) - lgamma(nu/2) - log(a)
Now we have:
     n        power
[...]
11 154 0.8033132631
12 344 0.9839843189
13 346 0.9844217782
14 500 0.9982834559

This is one example of why open source is so beautiful :love:.
—
Regards,

Detlew
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2009-10-08 00:36
(5680 d 00:44 ago)

@ d_labes
Posting: # 4319
Views: 32,149
 

 package MBESS

Post reply
Dear D Labes!

❝ This is one example of why open source is so beautiful.


Have you considered contacting Kem Phillips ([image] kemphillips[image]comcast.net) – also about the ambiguities in the documentation?
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
d_labes
★★★

Berlin, Germany,
2009-10-15 16:57
(5672 d 08:23 ago)

@ Helmut
Posting: # 4366
Views: 32,110
 

 Package MBESS and Power curiosity

Post reply
Dear Helmut! Dear all R adepts!

Helmut, as you suggested I have contacted Kem Phillips.
He has confirmed that sigma in case of logscale=TRUE is sqrt(log(1+CV^2)).
He has also confirmed the error if n>344 and will change the code in MBESS accordingly.

With respect to the power curiosity in my post "The Power at limits" above we are still in an ongoing discussion. Kem's first supposition was numerical problems. But considering the results I will present below he is thoughtful.

Because PASS2008 gives the same results I have confronted Jerry Hinze with that curiosity. Here his answer:
"This is surprising to me, but since it occurs in a range that is of little interest, I don't see the need to investigate it further. It could be due to a number of numerical problems that occur when dealing with such a small sample size. However, since as you say, it matches SAS Analyst, I would conjecture that it is probably due to the definition of the equivalence test--that is, the definition probably leads to this result."

For me I had seen the need to investigate it further.:-P
I performed an "in-silico" experiment to obtain an empirical power estimate.
I have simulated data sets for 2x2 cross-over with equal subjects in the two sequences using my beasty number cruncher SAS, performed the TOST and recorded the incidences of stating bioequivalence.
Here my results:

ratio0=0.95, CV=0.65, BE limits 0.8 ... 1.25

                Power       SAS      EatMyShorts,    empirical
     Power  (R-code with  magic DL   Apfelstrudel,    power
 n  (MBESS)   Owen's Q)  (Owen's Q)  trapez meffud (50.000 simul.)
----------------------------------------------------------------
 4  0.004542  0.004542    0.004542     0.004542      0.00452.
 6  0.001620  0.001620    0.001620     0.001620      0.00146.
 8  0.000976  0.000976    0.000976     0.000976      0.00086.
10  0.000799  0.000799    0.000799     0.000799      0.00080.
12  0.000797  0.000797    0.000797     0.000797      0.00074.
14  0.000907  0.000907    0.000907     0.000907      0.00084.
16  0.001131  0.001131    0.001131     0.001131      0.00110.
18  0.001501  0.001501    0.001501     0.001501      0.00132.
20  0.002077  0.002077    0.002078     0.002077      0.00230.
22  0.002952  0.002952    0.002952     0.002952      0.00352.
24  0.004253  0.004253    0.004253     0.004253      0.00468.

30  0.012708  0.012708    0.012708     0.012708      0.01216.
40  0.058717  0.058717    0.058717     0.058715      0.05666.
50  0.157587  0.157587    0.157587     0.157587      0.15790.

BTW: ElMaestro, can you tell me how Apfelstrudel has found his way to Gent/Belgium :confused:?

As you can see, all results from the different implementations have the same behavior. So Jerry Hinze's presumption seems true that this behavior is inherent to the test definition.

BTW: The empirical power in that region badly converges to a stationary value. One has to use huge number of simulated datasets to obtain a reliable value to some degree. 50.000 is obviously yet to low.

Greetings to all other pettifoggers (Hair splitting needs sharp tools and haven't got anything better to do :-D ).
—
Regards,

Detlew
 
ElMaestro
★★★

Denmark,
2009-10-17 14:26
(5670 d 10:54 ago)

@ d_labes
Posting: # 4370
Views: 31,865
 

 Package MBESS and Power curiosity

Post reply
Ahoy d_labes,

❝ BTW: ElMaestro, can you tell me how Apfelstrudel has found his way to Gent/Belgium :confused:?


When I made port in Bremerhafen in 1953 I decided to give Apfelstrudel a try; it turned out to be a fully acceptable alternative to the famous Vlaamse Friet and with a pleasantly matching (lack of) nutritional value.

Regarding only 50000 simulations:
You can get the famous software EFG 2.01 from HS. It will do millions of sims in a few seconds. HS - put it up for download?

Here's an update to give you more decimals, if you like:

Apfelstrudel<-function(t, df, ta, tb, ncp, ACCURACY)
{
  L1 = sqrt(2*pi);
  P=0;
  step = (tb-ta)/ACCURACY;
  foo=ta;
  for (i in 1:ACCURACY)
    {
      y1 = EatMyShorts(foo, df, t, ncp);
      y2 = EatMyShorts(foo+0.5*step, df,t, ncp);
      y3 = EatMyShorts(foo+step, df,t, ncp);
      P=P+step*(y1+y2*4+y3)/6; ## homer simpson's meffud!
      foo=foo+step;
    }
  L1*P
}


Best regards
ElMaestro - sailor and simulant.

Edit: BBCoded. [Helmut]
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2009-10-17 15:44
(5670 d 09:36 ago)

@ ElMaestro
Posting: # 4371
Views: 32,465
 

 Apfelstrudel and EFG 2.01

Post reply
¡Hola ElMaestro!

❝ When I made port in Bremerhafen in 1953 I decided to give Apfelstrudel a try; it turned out to be […] a pleasantly matching (lack of) nutritional value.


'Apfelstrudel' in Bremerhafen where waitresses can't even pronounce the word correctly (I would except to hear something similar to Appelstrudel - whereas the correct pronounciation is Apflschdrudl)? :vomit:
Should give the real stuff (of course together with whipped cream) at try in Vienna. 400 kcal/serving are nourishing sailors well.

❝ You can get the famous software EFG 2.01 from HS […] – put it up for download?


OK. I bundeled the executable EFG 2.01 together with the manual of EFG 1.02. After download (MD5: 892f7df65887edc7e1e4a1817ff3d90c) unzip to a folder of your choice. Should run on Windows (≥95); no installation required. Couldn’t test Win95, but WfW 3.11 with Win32-extensions did not work. :cool:

❝ It will do millions of sims in a few seconds.


Well – a couple of seconds. ;-)
Don’t get nervous – it never crashed my systems (Windows 2000 Pro SP4, XP/Pro SP3, and Vista/Pro SP2)*, but doesn’t give a clue about it’s progress. Don’t expect the infamous turning hourclass or a beep when the sim is finished. As software engineers use to say:

‘Quick and dirty – or clean and never.’


BTW, no other software I know can give you a sample size estimation for Denmark’s requirement (CI includes 100%) out of the box!

ratio0=0.95, CV=0.65, BE limits 0.8 ... 1.25
    EFG 2.01   EFG 2.01
    Normal     Brute force
 n             (2×106)
--------------------------
 4  0.0045418  0.008385
 6  0.0016204  0.004700
 8  0.0009763  0.000478
10  0.0007993  0.0004685
12  0.0007967  0.0005085
14  0.0009067  0.000647
16  0.0011306  0.0008795
18  0.0015009  0.0012605
20  0.0020775  0.001691
22  0.0029517  0.002559
24  0.0042529  0.003802
30  0.0127083  0.011871
40  0.0587169  0.058030
50  0.1575871  0.1570075

[image]

  • Runs on Win7 SP1 64bit! [Helmut]
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
d_labes
★★★

Berlin, Germany,
2009-10-20 13:35
(5667 d 11:45 ago)

@ ElMaestro
Posting: # 4372
Views: 32,350
 

 EFG and another power curiosity

Post reply
Ahoy old simulant!

❝ Regarding only 50000 simulations:

❝ You can get the famous software EFG 2.01 from HS. It will do millions of sims in a few seconds.


I have given EFG2.1 a try. Great :clap:.
Thanks for that stress test for my CPU. But my SAS can better :-D .

What I'm a little bit curious of is that regardless how big the number of sims is, the power for low sample size does not converge to your 'normal' values (at least for CV=30%) .
(This is the 'exact' power I think, using your 'Abbelssstrudel' ?).
Here an example (power in %):
CV=30%, ratio=0.95, BE limits= 0.8 - 1.25 (fixed in EFG?)
      ------ EM's EFG2.1----------  poor mans
       normal     grand brute         SAS®
  N  ('exact')  100.000  1.000.000  10.000
  4    3.42      1.31     1.28       3.33
  6    3.95      2.74     2.64       3.93
  8    5.96      4.87     4.90       5.56
 10    9.54      8.72     8.67       9.32
 12   14.85     14.21    14.20      15.08 / 14.47 (two runs)
My own simulations (but of course with lower number of sims because of CPU overheating) do not show this effect.

I would expect a convergence to the 'exact' values because we simulate data with distributional characteristics that match those for which the exact power formula was derived. Do I miss somefink here?
—
Regards,

Detlew
 
ElMaestro
★★★

Denmark,
2009-11-23 20:03
(5633 d 04:16 ago)

@ d_labes
Posting: # 4384
Views: 31,800
 

 EFG and another power curiosity

Post reply

❝ Ahoy old simulant!


Ahoy, d_labes, perhaps I should get a contract in Hollywood?

❝ What I'm a little bit curious of is that regardless how big the number of sims is, the power for low sample size does not converge to your 'normal' values (at least for CV=30%) .

❝ (...) I would expect a convergence to the 'exact' values because we simulate data with distributional characteristics that match those for which the exact power formula was derived.


I guess you are right, thanks for pointing it out. And yes absolutely, simulations of this type MUST converge towards the true value with increasing iterations. I will head down in the programming department and punish the perpetrator right away! Nuts. At any rate, I will post a little later a proposal on how to make the Brute Force parts of EFG 2.01 obsolete.

Sorry I did not see your post before. Or should I say, I am sorry the forum admin screwed up, making it impossible to access the forum, thereby putting a halt to rapid progression of the science that works in the interest of the sick. :-D

EM.
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2009-11-23 20:17
(5633 d 04:02 ago)

@ ElMaestro
Posting: # 4386
Views: 31,795
 

 EFG and another power curiosity

Post reply
Hi ElMaestro,

nice to hear from you again!

❝ […] Or should I say, I am sorry the forum admin screwed up, […]



Don't dare! :angry:
It was a hell of work to get the f**ing stuff up again.

❝ […] thereby putting a halt to rapid progression of the science that works in the interest of the sick. :-D


Yeah, that's what it is all about in Big Pharma, isn't it?
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2009-10-07 02:08
(5680 d 23:12 ago)

@ ElMaestro
Posting: # 4307
Views: 32,252
 

 non-central t, power, R:

Post reply
Dear all

I got power = 0.0016 with PASS as the attached picture. Funny.

[image]

❝ In a previous thread I asked what the power in a 2,2,2-BE trial is

❝ at N=6, CV=65% and T/R=95%.

❝ The things is, Fartssie comes up with a negative reply. HS's

❝ implementation in R of the equations given by Julious...

—
All the best,
-- Yung-jin Lee
bear v2.9.2:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
 
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2009-10-07 03:16
(5680 d 22:04 ago)

@ yjlee168
Posting: # 4308
Views: 33,021
 

 PASS

Post reply
Dear Yung-jin!

❝ I got power = 0.0016 with PASS as the attached picture. Funny.


Crazy. From the manual I guessed that PASS also uses AS 243 and would give results comparable with Patterson's/Jones' SAS code (also mine ported to R) and FARTSSIE. But from this post we know that PASS' result agrees with D Labes' magic SAS-homebrew...
That's another example supporting Jaime's advice: Never trust in any piece of software you haven't written yourself (and even then you should be cautious...) Was this the reason to write bear?
—
Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
 
yjlee168
★★★
avatar
Homepage
Kaohsiung, Taiwan,
2009-10-13 01:17
(5675 d 00:03 ago)

@ Helmut
Posting: # 4344
Views: 32,482
 

 PASS

Post reply
Dear Helmut,

❝ [..] you should be cautious...) Was this the reason to write bear?


The reason to develop bear is to simply try to prove that R can do BE/BA too. We have a lot of funs from developing work. Lots to read and learn. That's been good enough for us. At least we have an open sourced package now for BE/BA community.
—
All the best,
-- Yung-jin Lee
bear v2.9.2:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here
 
ElMaestro
★★★

Denmark,
2009-10-09 13:32
(5678 d 11:48 ago)

@ ElMaestro
Posting: # 4332
Views: 32,317
 

 R code for power

Post reply
What comes out of lying sleepless at night:


EatMyShorts <- function (x, V, T, W)
{
 
  l1 = pnorm((T*x/sqrt(V)) - W);
 
  l2 = (V-1.0)*log(x) -lgamma(0.5*V) - (0.5*(V-2.0))* log (2.0);
  l3 = exp(l2);
  l4 =dnorm(x);
  l1*l3*l4;
 
}

Apfelstrudel<-function(t, df, ta, tb, ncp)
{
  L1 = sqrt(2*pi);
  P=0;
  step = (tb-ta)/500;
  foo=ta;
  for (i in 1:500)
    {
      y1 = EatMyShorts(foo, df, t, ncp);
      y2 = EatMyShorts(foo+step, df,t, ncp);
      P=P+step*0.5*(y1+y2); ## the trapez meffud!
      foo=foo+step;
    }
  L1*P
}


Pwr222 <- function (N, alpha, AccLimitLo, AccLimitHi, TRRatio, CV)
{
  df = N-2;
  tc = qt( 1-alpha, df);
  k  = sqrt(2/N);
  d1 = (log(TRRatio)-log(AccLimitLo))/(k * sqrt(log(1+CV*CV)));
  d2 = (log(TRRatio)-log(AccLimitHi))/(k * sqrt(log(1+CV*CV)));
  R  = sqrt(df) * (d1-d2) / (2*tc);
 
  P2 = Apfelstrudel(-tc, df, 0,R,d2);
  P1 = Apfelstrudel(tc, df, 0,R,d1);
 
  P2-P1
}


Increasing the integration steps (to some number higher than 500) may ímprove accuracy, as will e.g. a method like Simpson's formula. Can be played around with.

Best regards
EM.
 
Ohlbe
★★★

France,
2009-10-09 13:50
(5678 d 11:30 ago)

@ ElMaestro
Posting: # 4333
Views: 32,178
 

 Belgian beer ?

Post reply
Dear ElMAestro

❝ What comes out of lying sleepless at night:

❝ EatMyShorts <- function (x, V, T, W)

❝ Apfelstrudel<-function(t, df, ta, tb, ncp)


:-D

Regards
Ohlbe
—
Regards
Ohlbe
 
d_labes
★★★

Berlin, Germany,
2009-10-09 15:50
(5678 d 09:30 ago)

@ ElMaestro
Posting: # 4335
Views: 32,700
 

 Another one - R code for power of TROST

Post reply
Dear ElMaestro, dear All!

Also suffering from senile "Bettflucht" (insomnia) here is my quick shoot with Owen's Q:
(forgive me, I'm a very beginning adept of R)

# The functions of normal-, t-distributions and integrate()
require(stats)
# -----------------------------------------------------------------------------
# Owen's Q-function
OwensQ <- function (nu, t, delta, a, b)
{
   # .Machine$double.eps^.5 = 1.490116e-08 on my machine
   # What method integrate() uses is not clear to me. But seems sophisticated
   # to assure enough accuracy
   Qintegral <- integrate(Q.integrand, lower = a, upper = b,
                nu=nu, t=t, delta = delta,
                subdivisions = 1000,
                rel.tol = .Machine$double.eps^.5, abs.tol = .Machine$double.eps^.5,
                stop.on.error = TRUE, keep.xy = FALSE, aux = NULL)
   Q <- Qintegral[[1]]
   return(Q) 
}
# -----------------------------------------------------------------------------
# Integrand of the definit integral in Owen's Q. Used in the call of integrate()
# Not useful alone, I think ?     
Q.integrand <- function(x, nu, t, delta)
{
   # assure a vector answer, copied from Power.density.md MBESS
   # look for better implementation! is loop necessary?
   pdim <- length(x)
   dens <- c(1:pdim)
   # constants for the loop
   sqrtnu   <- sqrt(nu)
   lnQconst <- 0.5*log(2*pi)-lgamma(nu/2)-((nu-2)/2)*log(2)

   for (ix in 1:pdim) {
     xx <- x[ix]
     if (xx > 0) {
       d <- pnorm(t*xx/sqrtnu - delta, mean = 0, sd = 1, log = TRUE) +
            (nu-1)*log(xx) + dnorm(xx, mean = 0, sd = 1, log = TRUE)
       dens[ix] <- exp(lnQconst + d)
     } else {
       # do we need  < 0 ? or must xx>=0? for our power it is.
       dens[ix] <-0
     }
   }

   return(dens)
}

# -----------------------------------------------------------------------------
# Power of two-one-sided-t-tests using OwensQ
# In case of logscale=TRUE give ldiff, ltheata1 and ltheta2 as ratios
# f.i. ldiff=0.95, ltheta1=0.8, ltheta2=1.25
# in case of logscale=FALSE give ldiff, ltheata1 and ltheta2 as difference to 1
# f.i. ldiff=0.05 (5% difference) ltheata1=-0.2, ltheta2=0.2 (20% equivalence margins)
# CV is always the coefficient of variation but as ratio, not %

Power.TOST <- function(alpha, logscale, ltheta1, ltheta2, ldiff, CV, n, df)
{
    theta1 <- ltheta1
    theta2 <- ltheta2
    diff   <- ldiff
    s      <- CV
    if (logscale) {
        theta1 <- log(ltheta1)
        theta2 <- log(ltheta2)
        diff   <- log(ldiff)
        s <- sqrt(log(1+CV*CV))
    }
    tval <- qt(1 - alpha, df, lower.tail = TRUE, log.p = FALSE)
    delta1 <- (diff-theta1)/(s*sqrt(2/n))
    delta2 <- (diff-theta2)/(s*sqrt(2/n))
    R <- (delta1-delta2)*sqrt(df)/(2*tval)
    p1 <- OwensQ(df,  tval, delta1, 0, R)
    p2 <- OwensQ(df, -tval, delta2, 0, R)

    power=p2-p1

    return(power)
}

Some spare 'validation' results (Sample sizes from Diletti et. al. for target power 70%, 80%, 90%):
CV=0.3, BE margins 0.8-1.25, Null ratio=0.95, alpha=0.05
   n power(Q) Apfelstrudel
1 32 0.7177182  0.7177182
2 40 0.8158453  0.8158453
3 52 0.9019652  0.9019652

PS: I think it's enough R power in this thread to have it in the R category now.

Edit: Category changed. require(stats) not needed, because loaded by default. [Helmut]
—
Regards,

Detlew
 
UA Flag
Activity
 Admin contact
23,424 posts in 4,927 threads, 1,668 registered users;
19 visitors (0 registered, 19 guests [including 1 identified bots]).
Forum time: 01:20 CEST (Europe/Vienna)

A drug is that substance which, when injected into a rat,
will produce a scientific report.    Anonymous

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5