Bioequivalence and Bioavailability Forum • R code for power

R code for power [🇷 for BE/BA]

posted by ElMaestro – Denmark, 2009-10-09 13:32 (5731 d 04:31 ago) – Posting: # 4332
Views: 33,995

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.

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 2 d_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 powerElMaestro 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