Suggested code [Two-Stage / GS Designs]

posted by ElMaestro  – Belgium?, 2017-08-21 16:13 (1078 d 10:25 ago) – Posting: # 17729
Views: 24,363

(edited by ElMaestro on 2017-08-21 16:37)

Hi Hötzi,


I imagine I was not able to explain what I meant??
I looked at the code available here.


Do this:

1. In the R code for the Potvin function, include this code:

Power.calc=function(Nps, GMR, CV, Is.St2=0)
{
  s=sqrt(log(1+CV*CV))
  nom1=log(1.25/GMR)
  nom2=-log(1.25*GMR)
  den=s*sqrt(2/(2*Nps))
  df=Nps*2-2
  if (Is.St2==T) df=df-1
  ##cat("1-a=", 1-alpha, "df=",df,"\n")
  if (Is.St2==0) q=vecQT1[df]
    else         q=vecQT2[df]
  pw=pt( (nom1/den) -q, df) - pt( (nom2/den) +q, df)
  if (pw<0) pw=0.01
  return(pw)
}

FindCritCV=function(Pwr.T, GMR, Nps, Is.St2=1, toler=0.000001)
{
  CV1=0.0123
  CV2=8.7654
  while (abs(
         Power.calc(Nps, GMR, CV1, Is.St2) -
         Power.calc(Nps, GMR, CV2, Is.St2))>toler)
  {
   CV3=0.5*(CV1+CV2)
   tmp=Pwr.T-Power.calc(Nps, GMR, CV3, Is.St2)

   if (tmp>0) CV2=CV3
      else    CV1=CV3
  }
  return(0.5*(CV2+CV1) )
}



2. immediately before the line that says "if (setseed) set.seed(1234567)", include this code:
MAX=1000 ##Hötzi, specify the number of values you need,
         ##I don't know it from the top of my head.
vecQT1<<-c(1:MAX)
vecQT2<<-c(1:MAX)
for (i in 1:MAX)
   {
    vecQT1[i]<<-qt(p=1-alpha[1], df=i)
    vecQT2[i]<<-qt(p=1-alpha[2], df=i)
   }


vecCritCVs<<-rep(0,MAX) 
for (Nps in 6:length(vecCritCVs))
  vecCritCVs[Nps]<<-FindCritCV(0.8, 0.95, Nps, 1, 1e-10)
  ##you can tweak this line with user setting for
  ##assumed GMR and target power as you please,
  ##using the values from the function call


3. Replace the code for SampleN.Stg2 with the code from my previous post:

» sampleN.Stg2.EM <- function(CV)
» {
»   #all other variable info is already accumulated thin the critical values, so:
»   Nps=1
»   while (vecCritCVs[Nps]<=CV) Nps=Nps+1
»   return(Nps*2) 
» }

(and make sure you call this function, possibly rename it as you deem appropriate)

4. Tweak and benchmark it. The relevant benchmark is the runtime for the entire Potvin call with 1.000.000 iterations at e.g. CV=0.3 and N1=36 total or something which takes tangible computation time. You will also need to modify the Power.calc function to accomodate other acceptance ranges.

5. Eliminate all other qt calls by replacing them with lookups into the vecQT1 or vecQT2 vectors.



Note: I did not test it, I only looked at the code online, so not everything here can be assumed to be working out of the box.

I could be wrong, but...

Best regards,
ElMaestro

"Pass or fail" (D. Potvin et al., 2008)

Complete thread:

Activity
 Admin contact
20,964 posts in 4,373 threads, 1,459 registered users;
online 19 (0 registered, 19 guests [including 16 identified bots]).
Forum time: Tuesday 02:38 UTC (Europe/Vienna)

My doctor gave me six months to live,
but when I couldn’t pay the bill
he gave me six months more.    Walter Matthau

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