Helmut
★★★
avatar
Homepage
Vienna, Austria,
2015-02-20 02:47
(3324 d 10:29 ago)

Posting: # 14466
Views: 12,463
 

 Gimmick & survey [Two-Stage / GS Designs]

Dear all,

I’m still exploring R’s graphic capabilities…

[image]

Can you guess which plot shows the type I errors of which design?
Bonus question: What can we learn from these plots?


For initiates: CV 10–80% (step 2%), n1 12–72 (2) ⇒ 1.116·109 simulations/plot. Power2Stage is amazing.

Edit: Unfair, bloody unfair! Poor TOST.
[image]

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
nobody
nothing

2015-02-20 09:25
(3324 d 03:50 ago)

@ Helmut
Posting: # 14467
Views: 11,459
 

 Gimmick & survey

Hi Vienna!

Totally off topic:

http://link.springer.com/article/10.1007/s00228-015-1806-2?no-access=true

...you should publish open-source... ;-)

(Springer handled the additional material .doc quite loveless, btw)

OT:
First graph = Potvin B
Second graph = Potvin C
Third graph = conventional 2x2 ?

Kindest regards, nobody
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2015-02-20 14:01
(3323 d 23:15 ago)

@ nobody
Posting: # 14471
Views: 11,460
 

 Gimmick & survey

Hi Sweden (or Germany?)

❝ ...you should publish open-source... ;-)


Everybody – who could afford it – should do so. It was an invited review and assessing the endless list of references and running the simulations took me close to two months. My one-man show didn’t want to spend € 2,200 for OA.

❝ (Springer handled the additional material .doc quite loveless, btw)


Yes. :-(
I had some discussions with the publisher about the format of the plots in the paper only to find them downscaled making some legends almost illegible…

❝ First graph = Potvin B

❝ Second graph = Potvin C

❝ Third graph = conventional 2x2 ?


Almost – almost – yes.
Potvin’s adjusted α of 0.0294 for both methods was unfortunate. “Type 2” always needs more ad­just­ment than “type 1”. I used the optimized values of 0.0303 and 0.0282. No significant inflation of the TIE (maxima 0.05015 and 0.05010).

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
nobody
nothing

2015-02-20 15:18
(3323 d 21:58 ago)

@ Helmut
Posting: # 14472
Views: 11,446
 

 Gimmick & survey

Hi from nowhere!

❝ Potvin’s adjusted α of 0.0294 for both methods was unfortunate. “Type 2” always needs more ad­just­ment than “type 1”. I used the optimized values of 0.0303 and 0.0282. No significant inflation of the TIE (maxima 0.05015 and 0.05010).


Clever! :-D

❝ ...and running the simulations took me close to two months.


Sorry for OT again, but I heard that all the "for..." constructions are not so effective in R and better to use vectors for all "for ...to" and calculate them. Have you any experience on this?

(Read the R developers are thinking about an API for distributed computing. so we can share all our machines for large simulations :-D)

Kindest regards, nobody
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2015-02-20 16:38
(3323 d 20:38 ago)

@ nobody
Posting: # 14473
Views: 11,454
 

 Speeed!

Hi to nowhere!

❝ ❝ I used the optimized values of 0.0303 and 0.0282. No significant inflation of the TIE (maxima 0.05015 and 0.05010).


❝ Clever! :-D


THX. Not new (see this thread).

❝ ❝ ...and running the simulations took me close to two months.


Last year some options were not implemented in Power2Stage and I had to write my own – clumsy – code…

❝ I heard that all the "for..." constructions are not so effective in R and better to use vectors for all "for ...to" and calculate them.


Yes. Sometimes the performance boost can by 100×.

runs  <- 1e5
run   <- 1:runs
ptm   <- proc.time()
v1    <- NULL
for(j in seq_along(run)) v1 <- c(v1, rnorm(1, mean=0, sd=1))
t1    <- proc.time()-ptm
ptm   <- proc.time()
v2    <- vector("numeric", length=runs)
for(j in seq_along(run)) v2[j] <- rnorm(1, mean=0, sd=1)
t2    <- proc.time()-ptm
ptm   <- proc.time()
v3    <- rnorm(number, mean=0, sd=1)
t3    <- proc.time()-ptm
cat("\n",
"increasing vector:", t1[3], "seconds\n",
"indexed vector   :", t2[3], "seconds\n",
"internal function:", t3[3], "seconds\n")


 increasing vector: 11.74 seconds
 indexed vector   : 0.29 seconds
 internal function: 0.02 seconds


❝ Have you any experience on this?


Functions in PowerTOST and Power2Stage are completely vectorized and pre-compiled. With my latest code finding a suitable adjusted α and validating for TIE and power generally takes less than 30 minutes on my machine. Example for my “type 2” at the location of the maximum TIE. Four seconds for 106 sim’s!

library(Power2Stage)
power.2stage(method="C", alpha=rep(0.0282, 2), CV=0.1, n1=16,
             theta0=1.25, nsims=1e6)

Total time consumed (secs):
   user  system elapsed
    3.9     0.1     4.0


Method C: alpha0 = 0.05, alpha (s1/s2) = 0.0282 0.0282
Target power in power monitoring and sample size est. = 0.8
BE margins = 0.8 ... 1.25
CV = 0.1; n(stage 1)= 16; GMR = 0.95
GMR = 0.95 and mse of stage 1 in sample size est. used
Futility criterion Nmax = Inf

1e+06 sims at theta0 = 1.25 (p(BE)='alpha').
p(BE)    = 0.0501
p(BE) s1 = 0.0501
Studies in stage 2 = 0%

Distribution of n(total)
- mean (range) = 16 (16 ... 22)
- percentiles
 5% 50% 95%
 16  16  16


❝ (Read the R developers are thinking about an API for distributed computing. so we can share all our machines for large simulations :-D)


Yes, that would be nice! It’s already possible to run R in a local cloud-computing environment. We tried that at the Uni Lancaster and the speed boost was 20–50. With a SETI-like approach we could even think about TSDs for RSABE…

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
nobody
nothing

2015-02-20 17:37
(3323 d 19:39 ago)

@ Helmut
Posting: # 14475
Views: 11,351
 

 Speeed!

❝ ❝ Clever! :-D


❝ THX. Not new (see this thread).


Yeeaaaaahhhh, now I remember that... even ran the code in R sometimes... ;-)

library(Power2Stage)

power.2stage(method="C", alpha=rep(0.0282, 2), CV=0.1, n1=16, theta0=1.25, nsims=1e6)


Total time consumed (secs):

   user  system elapsed

    3.9     0.1     4.0


OK, 4.9s here...almost 25% down...but not a problem at all

❝ ...With a SETI-like approach we could even think about TSDs for RSABE…


Let's start with distributed computing in, let's say 10 to 20 instances on different computers..

Kindest regards, nobody
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2015-02-20 17:58
(3323 d 19:17 ago)

@ nobody
Posting: # 14476
Views: 11,436
 

 Speeed!

Hi nobody,

❝ Yeeaaaaahhhh, now I remember that... even ran the code in R sometimes... ;-)


Wouldn’t recommend that any more. If you want I can send you my current version (in­cluding futility criteria). Parallel designs not implemented yet.
  • In the primary step instead of assessing the maxima for each n1, I first look at the global maximum and run the second step for each TIE exceeding f.i. 95% of the glo­bal one. That’s better especially for “type 2” designs which have a flat TIE surface.
  • Anders’s linear model is not always the best. Now I fit both a linear and quadratic and select the better one based on the minimum AIC.

❝ ❝ ...With a SETI-like approach we could even think about TSDs for RSABE…


❝ Let's start with distributed computing in, let's say 10 to 20 instances on different computers..


Once the API is ready… On my machine I often use three instances of R without problems.

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
nobody
nothing

2015-02-20 18:01
(3323 d 19:14 ago)

@ Helmut
Posting: # 14477
Views: 11,396
 

 Speeed!

❝ If you want I can send you my current version (in­cluding futility criteria). Parallel designs not implemented yet.


Would be nice, something to play with over the weekend! :-D

Many thanx in advance...

Kindest regards, nobody
ElMaestro
★★★

Denmark,
2015-02-21 00:48
(3323 d 12:27 ago)

@ Helmut
Posting: # 14480
Views: 11,326
 

 Speeed!

Hi Hötzi,

❝ for(j in seq_along(run)) v1 <- c(v1, rnorm(1, mean=0, sd=1))


This is a funny way to do it. There will be constant memory re-allocation of the vector this way. I guess that's what takes the time.

Try and check if a technique called loop unrolling also helps speed things up in R. Sometimes it works miracles but is very implementation-specific. It has a little effect in my C implementations, but nothing really substantial so far.

The idea is (and pardon me for going from a for loop to a while loop; this is just because I don't know of a way to make a for loop increase by more than one per iteration, but I am sure you can do that):

In stead of

i=0
while (i<10)
{
  Do.Something(i)
  i=i+1
}


we can try

i=0
while (i<10)
{
   Do.Something(i)
   Do.Something(i+1)
   i=i+2
}


Here in stead of looping 10 times, we just need to loop 5 times. You can of course unroll a lot more than just x2 as the above example. One obvious drawback is that code can appear visually bloated but this of course doesn't affect functionality.

Pass or fail!
ElMaestro
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2015-02-21 02:11
(3323 d 11:04 ago)

@ ElMaestro
Posting: # 14481
Views: 11,286
 

 Speeed!

Hi ElMaestro,

❝ ❝ for(j in seq_along(run)) v1 <- c(v1, rnorm(1, mean=0, sd=1))


❝ This is a funny way to do it. There will be constant memory re-allocation of the vector this way. I guess that's what takes the time.


Exactly. That’s really stupid (see the R Inferno Circle 2). Detlew cured me.

❝ The idea is (and pardon me for going from a for loop to a while loop; this is just because I don't know of a way to make a for loop increase by more than one per iteration, but I am sure you can do that):


I have learned that one never ever should try to manipulate the counter of a for-loop. ;-)

runs <- 1e7
ptm  <- proc.time()
i    <- 1
while (i <= runs) {
  # Do.Something(i)
  i <- i+1
}
t1   <- proc.time()-ptm
ptm  <- proc.time()
i    <- 1
while (i <= runs) {
  # Do.Something(i)
  # Do.Something(i+1)
  i <- i+2
}
t2   <- proc.time()-ptm
cat("\n",
"simple :", t1[3], "seconds\n",
"partial:", t2[3], "seconds\n")

 simple : 5.07 seconds
 partial: 2.55 seconds


But: Quite often the number of iterations is not known beforehand (here it must be even) and the time-consuming part is the Do.Something().
Reminds me on my first Pascal-lessons almost 40 years ago.

Instead of    better use
y := 2*x;     y := x+x;   (addition faster than multiplication)
y := x/2;     y := 0.5*x; (multiplication faster than division)

Not kiddin’.

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,
2015-02-21 02:32
(3323 d 10:44 ago)

@ Helmut
Posting: # 14482
Views: 11,292
 

 Unrolled example

Hi Hötzi,

just tried to implement something that more or less corresponded to the example you had.

I get about a 10% speedup difference if I unroll 4x:

runs  <- 1e5
run   <- 1:runs
ptm   <- proc.time()
i=1
v2    <- vector("numeric", length=runs)
while (i<=runs)
{
  v2[i] <- rnorm(1, mean=0, sd=1)
  i=i+1
}
 
t2    <- proc.time()-ptm
cat("Simple time: ", t2[3], "seconds.\n")
#
#



runs  <- 1e5
run   <- 1:runs
ptm   <- proc.time()
i=1
v2    <- vector("numeric", length=runs)
while (i<=runs)
{
  v2[i]    <- rnorm(1, mean=0, sd=1)
  v2[i+1]  <- rnorm(1, mean=0, sd=1)
  v2[i+2]  <- rnorm(1, mean=0, sd=1)
  v2[i+3]  <- rnorm(1, mean=0, sd=1)

  i=i+4
}
 
t2    <- proc.time()-ptm
cat("Unrolled x4 time: ", t2[3], "seconds.\n")
#
#



10% isn't bad, actually. One could probably even optimise the number of unrolls.

Pass or fail!
ElMaestro
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2015-02-21 02:50
(3323 d 10:26 ago)

@ ElMaestro
Posting: # 14483
Views: 11,344
 

 Unrolled example

Hi ElMaestro,

❝ I get about a 10% speedup difference if I unroll 4x:


Confirmed. ;-)

❝ 10% isn't bad, actually. One could probably even optimise the number of unrolls.


Yeah, but at the end of the day it’s not necessarily the loop-overhead that counts. Such a code is terrible to maintain.


PS: In your unrolled code set runs <- 5 and after execution type length(v2). Kinetica-style.

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,
2015-02-21 03:27
(3323 d 09:49 ago)

@ Helmut
Posting: # 14484
Views: 11,291
 

 Unrolled example

Hi Hötzi,


PS: In your unrolled code set runs <- 5 and after execution type length(v2). Kinetica-style.


Wow, I like it :-D
In C you can get punished real bad if you do stuff like that. Looks like R knows that users make occasional mistakes and will reallocate the array length when the user tries to set something that is out of bounds.
Stuff like this restores my faith in humanity.

Pass or fail!
ElMaestro
Helmut
★★★
avatar
Homepage
Vienna, Austria,
2015-02-21 14:13
(3322 d 23:03 ago)

@ ElMaestro
Posting: # 14485
Views: 11,209
 

 Increasing vector & another example

Hi ElMaestro,

❝ In C you can get punished real bad if you do stuff like that.


Makes sense.

❝ Looks like R knows that users make occasional mistakes and will reallocate the array length when the user tries to set something that is out of bounds.


Couldn’t find anything about in the man-pages or on R-help. IMHO, R should at least throw a warning.

Another surprising goodie:

x   <- runif(1e3)
ops <- 1e6
ptm <- proc.time()
for(i in 1:ops) { }
t0  <- proc.time()-ptm
ptm <- proc.time()
for(i in 1:ops) { y1 <- mean(x) }
t1  <- 1e6*(proc.time()-ptm-t0)/ops
ptm <- proc.time()
for(i in 1:ops) { y2 <- sum(x)/length(x) }
t2  <- 1e6*(proc.time()-ptm-t0)/ops
cat("\n",
  sprintf("%14s: %.2f %s", "mean()", t1[3], "\u03bcs/op\n"),
  sprintf("%14s: %.2f %s", "sum()/length()", t2[3], "\u03bcs/op\n\n"))


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
nobody
nothing

2015-02-21 22:02
(3322 d 15:14 ago)

@ Helmut
Posting: # 14486
Views: 11,134
 

 Speeed!

...funny!

Kindest regards, nobody
UA Flag
Activity
 Admin contact
22,957 posts in 4,819 threads, 1,639 registered users;
70 visitors (0 registered, 70 guests [including 5 identified bots]).
Forum time: 13:16 CET (Europe/Vienna)

Nothing shows a lack of mathematical education more
than an overly precise calculation.    Carl Friedrich Gauß

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