Bioequivalence and Bioavailability Forum • Power at limits 2

Power at limits 2 [🇷 for BE/BA]

posted by d_labes – Berlin, Germany, 2009-10-07 17:15 (6089 d 15:13 ago) – Posting: # 4316
Views: 45,372

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

Complete thread:

non-central t, power, R: ElMaestro 2009-10-03 19:05 [🇷 for BE/BA]
- non-central t, power, R: Helmut 2009-10-03 21:45
  - non-central t, power, R: ElMaestro 2009-10-04 21:51
    - non-central t, power, R: Helmut 2009-10-04 23:58
- non-central t is not the question d_labes 2009-10-05 12:58
  - SAS -> R -> SAS Helmut 2009-10-05 14:53
    - Lost in translation d_labes 2009-10-06 14:31
      - Lost in translation Helmut 2009-10-06 14:42
        
        Macacus cynocephalus d_labes 2009-10-06 15:14
        
        Lost in translation Helmut 2009-10-06 21:46
  - non-central t is not the question ElMaestro 2009-10-05 17:32
    - Power of TOST in R d_labes 2009-10-06 16:45
      - Power of TOST (R package MBESS) Helmut 2009-10-07 02:18
        
        Power of TOST (R package MBESS) d_labes 2009-10-07 08:46
        
        Power of TOST (R package MBESS) Helmut 2009-10-07 12:38
      - The Power at limits d_labes 2009-10-07 11:32
        
        Power at limits 2d_labes 2009-10-07 15:15
        
        package MBESS Helmut 2009-10-07 22:36
        
        Package MBESS and Power curiosity d_labes 2009-10-15 14:57
        
        Package MBESS and Power curiosity ElMaestro 2009-10-17 12:26
        
        Apfelstrudel and EFG 2.01 Helmut 2009-10-17 13:44
        EFG and another power curiosity d_labes 2009-10-20 11:35
        
        EFG and another power curiosity ElMaestro 2009-11-23 19:03
        
        EFG and another power curiosity Helmut 2009-11-23 19:17
- non-central t, power, R: yjlee168 2009-10-07 00:08
  - PASS Helmut 2009-10-07 01:16
    - PASS yjlee168 2009-10-12 23:17
- R code for power ElMaestro 2009-10-09 11:32
  - Belgian beer ? Ohlbe 2009-10-09 11:50
  - Another one - R code for power of TROST d_labes 2009-10-09 13:50