[Gsdesign-commits] r321 - in pkg/gsDesign: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Aug 20 16:40:30 CEST 2011
Author: keaven
Date: 2011-08-20 16:40:29 +0200 (Sat, 20 Aug 2011)
New Revision: 321
Modified:
pkg/gsDesign/DESCRIPTION
pkg/gsDesign/NAMESPACE
pkg/gsDesign/R/gsDesign.R
pkg/gsDesign/R/gsMethods.R
pkg/gsDesign/R/gsSpending.R
pkg/gsDesign/R/gsqplot.R
pkg/gsDesign/src/gsbound.c
pkg/gsDesign/src/gsbound1.c
Log:
Allow not testing one bound at an interim (see sfTruncated)
Modified: pkg/gsDesign/DESCRIPTION
===================================================================
--- pkg/gsDesign/DESCRIPTION 2011-08-19 05:40:20 UTC (rev 320)
+++ pkg/gsDesign/DESCRIPTION 2011-08-20 14:40:29 UTC (rev 321)
@@ -1,11 +1,10 @@
Package: gsDesign
-Version: 2.5-02
+Version: 2.5-03
Title: Group Sequential Design
Author: Keaven Anderson
Maintainer: Keaven Anderson <keaven_anderson at merck.com>
Depends: ggplot2, xtable
-Description: gsDesign is a package that derives group sequential designs and describes their properties. 2.5 additions focus on predictive power update and
-prediction intervals.
+Description: gsDesign is a package that derives group sequential designs and describes their properties.
License: GPL (>= 2)
Copyright: Copyright 2010, Merck Research Laboratories
Packaged: Thu Aug 6 17:28:16 2009; anderkea
Modified: pkg/gsDesign/NAMESPACE
===================================================================
--- pkg/gsDesign/NAMESPACE 2011-08-19 05:40:20 UTC (rev 320)
+++ pkg/gsDesign/NAMESPACE 2011-08-20 14:40:29 UTC (rev 321)
@@ -7,6 +7,6 @@
export(normalGrid)
export(plot.gsDesign, plot.gsProbability, print.gsProbability, print.gsDesign)
export(print.nSurvival, xtable.gsDesign, gsBoundSummary)
-export(sfPower, sfHSD, sfExponential, sfBetaDist, sfLDOF, sfLDPocock, sfPoints, sfLogistic, sfExtremeValue, sfExtremeValue2, sfLinear)
+export(sfPower, sfHSD, sfExponential, sfBetaDist, sfLDOF, sfLDPocock, sfPoints, sfLogistic, sfExtremeValue, sfExtremeValue2, sfLinear, sfTruncated)
export(sfCauchy, sfNormal, sfTDist, spendingFunction)
export(checkScalar, checkVector, checkRange, checkLengths, isInteger)
Modified: pkg/gsDesign/R/gsDesign.R
===================================================================
--- pkg/gsDesign/R/gsDesign.R 2011-08-19 05:40:20 UTC (rev 320)
+++ pkg/gsDesign/R/gsDesign.R 2011-08-20 14:40:29 UTC (rev 321)
@@ -50,13 +50,15 @@
# check input arguments
checkVector(I, "numeric", c(0, Inf), c(FALSE, TRUE))
- checkVector(trueneg, "numeric", c(0,1), c(FALSE, FALSE))
- checkVector(falsepos, "numeric", c(0,1), c(FALSE, FALSE))
+ checkVector(trueneg, "numeric", c(0,1), c(TRUE, FALSE))
+ checkVector(falsepos, "numeric", c(0,1), c(TRUE, FALSE))
checkScalar(tol, "numeric", c(0, Inf), c(FALSE, TRUE))
checkScalar(r, "integer", c(1, 80))
checkLengths(trueneg, falsepos, I)
k <- as.integer(length(I))
+ if (trueneg[k]<=0.) stop("Final futility spend must be > 0")
+ if (falsepos[k]<=0.) stop("Final efficacy spend must be > 0")
r <- as.integer(r)
storage.mode(I) <- "double"
storage.mode(trueneg) <- "double"
@@ -81,7 +83,7 @@
checkScalar(theta, "numeric")
checkVector(I, "numeric", c(0, Inf), c(FALSE, TRUE))
checkVector(a, "numeric")
- checkVector(probhi, "numeric", c(0,1), c(FALSE, FALSE))
+ checkVector(probhi, "numeric", c(0,1), c(TRUE, FALSE))
checkScalar(tol, "numeric", c(0, Inf), c(FALSE, TRUE))
checkScalar(r, "integer", c(1, 80))
checkScalar(printerr, "integer")
@@ -89,6 +91,7 @@
# coerce type
k <- as.integer(length(I))
+ if (probhi[k]<=0.) stop("Final spend must be > 0")
r <- as.integer(r)
printerr <- as.integer(printerr)
Modified: pkg/gsDesign/R/gsMethods.R
===================================================================
--- pkg/gsDesign/R/gsMethods.R 2011-08-19 05:40:20 UTC (rev 320)
+++ pkg/gsDesign/R/gsMethods.R 2011-08-20 14:40:29 UTC (rev 321)
@@ -94,11 +94,11 @@
"print.gsDesign" <- function(x, ...)
{
- if (x$nFixSurv > 0)
- { cat("Group sequential design sample size for time-to-event outcome\n",
- "with sample size ", x$nSurv, ". The analysis plan below shows events\n",
- "at each analysis.\n", sep="")
- }
+ if (x$nFixSurv > 0)
+ { cat("Group sequential design sample size for time-to-event outcome\n",
+ "with sample size ", x$nSurv, ". The analysis plan below shows events\n",
+ "at each analysis.\n", sep="")
+ }
if (x$test.type == 1)
{
@@ -294,7 +294,7 @@
}
"sfprint" <- function(x)
-{
+{
# print spending function information
if (x$name == "OF")
{
@@ -308,6 +308,13 @@
{
cat("Wang-Tsiatis boundary with Delta =", x$param)
}
+ else if (x$name == "Truncated")
+ { cat(x$param$name,"spending function truncated at",x$param$trange)
+ if (!is.null(x$parname) && !is.null(x$param))
+ {
+ cat("\n with", x$parname, "=", x$param$param)
+ }
+ }
else
{
cat(x$name, "spending function")
Modified: pkg/gsDesign/R/gsSpending.R
===================================================================
--- pkg/gsDesign/R/gsSpending.R 2011-08-19 05:40:20 UTC (rev 320)
+++ pkg/gsDesign/R/gsSpending.R 2011-08-20 14:40:29 UTC (rev 321)
@@ -14,9 +14,9 @@
# sfLogistic
# sfNormal
# sfPoints
-# sfLinear
# sfPower
# sfTDist
+# sfTruncated
# spendingFunction
#
# Hidden Functions:
@@ -440,10 +440,15 @@
if (k > 1)
{ inctime <- x$param[1:k] - c(0, x$param[1:(k-1)])
incspend <- x$param[(k+1):j]-c(0, x$param[(k+1):(j-1)])
- if ((j > 2) && (min(inctime) <= 0 || min(incspend)<= 0))
+ if ((j > 2) && (min(inctime) <= 0))
{
- stop("Timepoints specified and cumulative spending must be strictly increasing in sfLinear")
+ stop("Timepoints must be strictly increasing in sfLinear")
}
+ if ((j > 2) && (min(incspend) < 0))
+ {
+ stop("Spending must be non-decreasing in sfLinear")
+ }
+
}
s <- t
s[t<=0]<-0
@@ -614,6 +619,31 @@
x
}
+"sfTruncated" <- function(alpha, t, param){
+ checkScalar(alpha, "numeric", c(0, Inf), c(FALSE, FALSE))
+ checkVector(t, "numeric", c(0, Inf), c(TRUE, FALSE))
+ if (!is.list(param)) stop("param must be a list. See help(sfTruncated)")
+ if (!max(names(param)=="trange")) stop("param must include trange, sf, param. See help(sfTruncated)")
+ if (!max(names(param)=="sf")) stop("param must include trange, sf, param. See help(sfTruncated)")
+ if (!max(names(param)=="param")) stop("param must include trange, sf, param. See help(sfTruncated)")
+ if (!is.vector(param$trange)) stop("param$trange must be a vector of length 2 with 0 <= param$trange[1] <param$trange[2]<=1. See help(sfTruncated)")
+ if (length(param$trange)!=2) stop("param$trange parameter must be a vector of length 2 with 0 <= param$trange[1] <param$trange[2]<=1. See help(sfTruncated)")
+ if (param$trange[1]>=1. | param$trange[2]<=param$trange[1] | param$trange[2]<=0)
+ stop("param$trange must be a vector of length 2 with 0 <= param$trange[1] < param$trange[2]<=1. See help(sfTruncated)")
+ if (class(param$sf) != "function") stop("param$sf must be a spending function")
+ if (!is.numeric(param$param)) stop("param$param must be numeric")
+ spend<-array(0,length(t))
+ spend[t>=param$trange[2]]<-alpha
+ indx <- param$trange[1]<t & t<param$trange[2]
+ s <- param$sf(alpha=alpha,t=(t[indx]-param$trange[1])/(param$trange[2]-param$trange[1]),param$param)
+ spend[indx] <- s$spend
+ param$name <- s$name
+ x<-list(name="Truncated", param=param, parname=s$parname,
+ sf=sfTruncated, spend=spend, bound=NULL, prob=NULL)
+ class(x) <- "spendfn"
+ x
+}
+
"spendingFunction" <- function(alpha, t, param)
{
checkScalar(alpha, "numeric", c(0, Inf), c(FALSE, FALSE))
Modified: pkg/gsDesign/R/gsqplot.R
===================================================================
--- pkg/gsDesign/R/gsqplot.R 2011-08-19 05:40:20 UTC (rev 320)
+++ pkg/gsDesign/R/gsqplot.R 2011-08-20 14:40:29 UTC (rev 321)
@@ -149,20 +149,25 @@
ntx <- "n="
if (is.null(xlab)) xlab <- "Sample size"
}
- if (x$test.type > 1)
- { z <- fn(z=c(x$upper$bound,x$lower$bound), i=c(1:x$k, 1:x$k), x=x,
- ratio=ratio, delta0=delta0, delta=delta)
- Ztxt <- as.character(c(round(z[1:(x$k-1)],dgt[2]), round(z[x$k],max(dgt)),
- round(z[(x$k+1):(2*x$k-1)], dgt[1]), round(z[2*x$k],max(dgt))))
+ if (x$test.type > 1)
+ { z <- fn(z=c(x$upper$bound, x$lower$bound), i=c(1:x$k,1:x$k), x=x,
+ ratio=ratio, delta0=delta0, delta=delta)
+ Ztxt <- as.character(c(round(z[1:x$k],dgt[2]), round(z[x$k+(1:x$k)], dgt[1])))
+ # use maximum digits for equal final bounds
+ if (x$upper$bound[x$k]==x$lower$bound[x$k])
+ { Ztxt[c(x$k,2*x$k)] <- round(z[x$k],max(dgt))
+ }
+ indxu <- (1:x$k)[x$upper$bound < 20]
+ indxl <- (1:x$k)[x$lower$bound > -20]
y <- data.frame(
- N=as.numeric(c(x$n.I,x$n.I)),
- Z=as.numeric(z),
- Bound=c(array("Upper", x$k), array("Lower", x$k)),
- Ztxt=Ztxt)
+ N=as.numeric(c(x$n.I[indxu],x$n.I[indxl])),
+ Z=as.numeric(z[c(indxu,x$k+indxl)]),
+ Bound=c(array("Upper", length(indxu)), array("Lower", length(indxl))),
+ Ztxt=Ztxt[c(indxu,x$k+indxl)])
}else{
z <- fn(z=x$upper$bound, i=1:x$k, x=x,
ratio=ratio, delta0=delta0, delta=delta)
- Ztxt <- as.character(round(z,dgt[1]))
+ Ztxt <- as.character(round(z),dgt[1])
y <- data.frame(
N=as.numeric(x$n.I),
Z=as.numeric(z),
Modified: pkg/gsDesign/src/gsbound.c
===================================================================
--- pkg/gsDesign/src/gsbound.c 2011-08-19 05:40:20 UTC (rev 320)
+++ pkg/gsDesign/src/gsbound.c 2011-08-20 14:40:29 UTC (rev 321)
@@ -1,113 +1,121 @@
-#define DEBUG 0
-#define EXTREMEZ 20
-#include "R.h"
-#include "Rmath.h"
-/* Group sequential probability computation per Jennison & Turnbull
- xnanal- # of possible analyses in the group-sequential designs
- (interims + final)
- I - statistical information available at each analysis
- a - lower cutoff points for z statistic at each analysis (output)
- b - upper cutoff points for z statistic at each analysis (output)
- probhi- input vector with probability of rejecting (Z>bj) at
- jth interim analysis, j=1...nanal
- problo- input vector with probability of rejecting (Z<aj) at
- jth interim analysis, j=1...nanal
- xtol - relative change between iterations required to stop for 'convergence'
- xr - determinant of # of grid points for numerical integration
- r=17 will give a max of 201 points which is what they recommend
- retval- error flag returned; 0 if convergence; 1 indicates error
- printerr- 1 if error messages to be printed - other values suppress printing
-*/
-void gsbound(int *xnanal,double *I,double *a,double *b,double *problo,double *probhi,
- double *xtol,int *xr,int *retval,int *printerr)
-{ int i,ii,j,m1,m2,r,nanal;
- double plo,phi,dplo,dphi,btem=0.,atem=0.,atem2,btem2,rtdeltak,rtIk,rtIkm1,xlo,xhi;
- double adelta,bdelta,tol;
-/* note: should allocat zwk & wwk dynamically...*/
- double zwk[1000],wwk[1000],hwk[1000],zwk2[1000],wwk2[1000],hwk2[1000],
- *z1,*z2,*w1,*w2,*h,*h2,*tem,rt2pi;
- void h1(double,int,double *,double,double *, double *);
- void hupdate(double,double *,int,double,double *, double *,
- int,double,double *, double *);
- int gridpts(int, double,double,double,double *, double *);
- r=xr[0]; nanal= xnanal[0]; tol=xtol[0]; rt2pi=2.506628274631;
-/* compute bounds at 1st interim analysis using inverse normal */
- if (nanal < 1 || r<1 || r>83)
- { retval[0]=1;
- if (*printerr)
- { Rprintf("gsbound1 error: illegal argument");
- if (nanal<1) Rprintf("; nanal=%d--must be > 0",nanal);
- if (r<1 || r> 83) Rprintf("; r=%d--must be >0 and <84",r);
- Rprintf("\n");
- }
- return;
- }
- a[0]=qnorm(problo[0],0.,1.,1,0);
- b[0]=qnorm(probhi[0],0.,1.,0,0);
-/* set up work vectors */
- z1=zwk; w1=wwk; h=hwk;
- z2=zwk2; w2=wwk2; h2=hwk2;
- m1=gridpts(r,0.,a[0],b[0],z1,w1);
- h1(0.,m1,w1,I[0],z1,h);
- rtIk=sqrt(I[0]);
- /* use Newton-Raphson to find subsequent interim analysis cutpoints */
- for(i=1;i<nanal;i++)
- { /* set up constants */
- rtIkm1=rtIk; rtIk=sqrt(I[i]); rtdeltak=sqrt(I[i]-I[i-1]);
- atem2=qnorm(problo[i],0.,1.,1,0); adelta=1.;
- btem2=qnorm(probhi[i],0.,1.,0,0); bdelta=1.;
- j=0;
- while((adelta>tol || bdelta>tol) && j++ < 20)
- { plo=0.; phi=0.; dplo=0.; dphi=0.;
- atem=atem2; btem=btem2;
- /* compute probability of crossing boundaries & their derivatives */
- for(ii=0;ii<=m1;ii++)
- { xlo=(z1[ii]*rtIkm1-atem*rtIk)/rtdeltak;
- xhi=(z1[ii]*rtIkm1-btem*rtIk)/rtdeltak;
- plo += h[ii]*pnorm(xlo,0.,1.,0,0);
- phi += h[ii]*pnorm(xhi,0.,1.,1,0);
- dplo+=h[ii]*exp(-xlo*xlo/2)/rt2pi*rtIk/rtdeltak;
- dphi-=h[ii]*exp(-xhi*xhi/2)/rt2pi*rtIk/rtdeltak;
- }
- /* use 1st order Taylor's series to update boundaries */
- /* maximum allowed change is 1 */
- /* maximum value allowed is z1[m1]*rtIk to keep within grid points */
- adelta=problo[i]-plo;
- if (adelta>dplo) atem2=atem+1.;
- else if (adelta < -dplo) atem2=atem-1.;
- else atem2=atem+(problo[i]-plo)/dplo;
- if (atem2>EXTREMEZ) atem2=EXTREMEZ;
- else if (atem2 < -EXTREMEZ) atem2= -EXTREMEZ;
- bdelta=probhi[i]-phi;
- if (bdelta<dphi) btem2=btem+1.;
- else if (bdelta > -dphi) btem2=btem-1.;
- else btem2=btem+(probhi[i]-phi)/dphi;
- if (btem2>EXTREMEZ) btem2=EXTREMEZ;
- else if (btem2< -EXTREMEZ) btem2= -EXTREMEZ;
- if (atem2>btem2) atem2=btem2;
- adelta=atem2-atem; if (adelta<0) adelta= -adelta;
- bdelta=btem2-btem; if (bdelta<0) bdelta= -bdelta;
- }
- a[i]=atem; b[i]=btem;
- /* if convergence did not occur, set flag for return value */
- if (adelta>tol ||bdelta > tol)
- { if (*printerr)
- { printf("gsbound1 error: No convergence for boundary for interim %d; I=%7.0lf",i+1,I[i]);
- if (bdelta>tol) printf("\n last 2 upper boundary values: %lf %lf\n",btem,btem2);
- if (adelta>tol) printf("\n last 2 lower boundary values: %lf %lf\n",atem,atem);
- }
- retval[0]=1;
- return;
- }
- if (i<nanal-1)
- { m2=gridpts(r,0.,a[i],b[i],z2,w2);
- hupdate(0.,w2,m1,I[i-1],z1,h,m2,I[i],z2,h2);
- m1=m2;
- tem=z1; z1=z2; z2=tem;
- tem=w1; w1=w2; w2=tem;
- tem=h; h=h2; h2=tem;
- } }
- retval[0]=0;
- return;
-}
-
+#define DEBUG 0
+/* note: EXTREMEZ > 3 + log(r) + Z(1-alpha) + Z(1-beta)
+ per bottom of p 349 in Jennison and Turnbull */
+#define EXTREMEZ 20
+#define MAXR 83
+#include "R.h"
+#include "Rmath.h"
+/* Group sequential probability computation per Jennison & Turnbull
+ xnanal- # of possible analyses in the group-sequential designs
+ (interims + final)
+ I - statistical information available at each analysis
+ a - lower cutoff points for z statistic at each analysis (output)
+ b - upper cutoff points for z statistic at each analysis (output)
+ probhi- input vector with probability of rejecting (Z>bj) at
+ jth interim analysis, j=1...nanal
+ problo- input vector with probability of rejecting (Z<aj) at
+ jth interim analysis, j=1...nanal
+ xtol - relative change between iterations required to stop for 'convergence'
+ xr - determinant of # of grid points for numerical integration
+ r=17 will give a max of 201 points which is what they recommend
+ retval- error flag returned; 0 if convergence; 1 indicates error
+ printerr- 1 if error messages to be printed - other values suppress printing
+*/
+void gsbound(int *xnanal,double *I,double *a,double *b,double *problo,double *probhi,
+ double *xtol,int *xr,int *retval,int *printerr)
+{ int i,ii,j,m1,m2,r,nanal;
+ double plo,phi,dplo,dphi,btem=0.,atem=0.,atem2,btem2,rtdeltak,rtIk,rtIkm1,xlo,xhi;
+ double adelta,bdelta,tol;
+/* note: should allocat zwk & wwk dynamically...*/
+ double zwk[1000],wwk[1000],hwk[1000],zwk2[1000],wwk2[1000],hwk2[1000],
+ *z1,*z2,*w1,*w2,*h,*h2,*tem,rt2pi;
+ void h1(double,int,double *,double,double *, double *);
+ void hupdate(double,double *,int,double,double *, double *,
+ int,double,double *, double *);
+ int gridpts(int, double,double,double,double *, double *);
+ r=xr[0]; nanal= xnanal[0]; tol=xtol[0]; rt2pi=2.506628274631;
+/* compute bounds at 1st interim analysis using inverse normal */
+ if (nanal < 1 || r<1 || r>MAXR)
+ { retval[0]=1;
+ if (*printerr)
+ { Rprintf("gsbound1 error: illegal argument");
+ if (nanal<1) Rprintf("; nanal=%d--must be > 0",nanal);
+ if (r<1 || r> MAXR) Rprintf("; r=%d--must be >0 and <84",r);
+ Rprintf("\n");
+ }
+ return;
+ }
+ if (problo[0] <= 0) a[0] = -EXTREMEZ;
+ else a[0]=qnorm(problo[0],0.,1.,1,0);
+ if (probhi[0] <= 0) b[0] = EXTREMEZ;
+ else b[0]=qnorm(probhi[0],0.,1.,0,0);
+/* set up work vectors */
+ z1=zwk; w1=wwk; h=hwk;
+ z2=zwk2; w2=wwk2; h2=hwk2;
+ m1=gridpts(r,0.,a[0],b[0],z1,w1);
+ h1(0.,m1,w1,I[0],z1,h);
+ rtIk=sqrt(I[0]);
+ /* use Newton-Raphson to find subsequent interim analysis cutpoints */
+ for(i=1;i<nanal;i++)
+ { /* set up constants */
+ rtIkm1=rtIk; rtIk=sqrt(I[i]); rtdeltak=sqrt(I[i]-I[i-1]);
+ if (problo[i]<=0.) atem2= -EXTREMEZ;
+ else atem2=qnorm(problo[i],0.,1.,1,0);
+ if (probhi[i]<=0.) btem2= EXTREMEZ;
+ else btem2=qnorm(probhi[i],0.,1.,0,0);
+ adelta=1.; bdelta=1.;
+ j=0;
+ while((adelta>tol || bdelta>tol) && j++ < EXTREMEZ)
+ { plo=0.; phi=0.; dplo=0.; dphi=0.;
+ atem=atem2; btem=btem2;
+ /* compute probability of crossing boundaries & their derivatives */
+ for(ii=0;ii<=m1;ii++)
+ { xlo=(z1[ii]*rtIkm1-atem*rtIk)/rtdeltak;
+ xhi=(z1[ii]*rtIkm1-btem*rtIk)/rtdeltak;
+ plo += h[ii]*pnorm(xlo,0.,1.,0,0);
+ phi += h[ii]*pnorm(xhi,0.,1.,1,0);
+ dplo+=h[ii]*exp(-xlo*xlo/2)/rt2pi*rtIk/rtdeltak;
+ dphi-=h[ii]*exp(-xhi*xhi/2)/rt2pi*rtIk/rtdeltak;
+ }
+ /* use 1st order Taylor's series to update boundaries */
+ /* maximum allowed change is 1 */
+ /* maximum value allowed is z1[m1]*rtIk to keep within grid points */
+ adelta=problo[i]-plo;
+ if (adelta>dplo) atem2=atem+1.;
+ else if (adelta < -dplo) atem2=atem-1.;
+ else atem2=atem+(problo[i]-plo)/dplo;
+ if (atem2>EXTREMEZ) atem2=EXTREMEZ;
+ else if (atem2 < -EXTREMEZ) atem2= -EXTREMEZ;
+ bdelta=probhi[i]-phi;
+ if (bdelta<dphi) btem2=btem+1.;
+ else if (bdelta > -dphi) btem2=btem-1.;
+ else btem2=btem+(probhi[i]-phi)/dphi;
+ if (btem2>EXTREMEZ) btem2=EXTREMEZ;
+ else if (btem2< -EXTREMEZ) btem2= -EXTREMEZ;
+ if (atem2>btem2) atem2=btem2;
+ adelta=atem2-atem; if (adelta<0) adelta= -adelta;
+ bdelta=btem2-btem; if (bdelta<0) bdelta= -bdelta;
+ }
+ a[i]=atem; b[i]=btem;
+/* if convergence did not occur, set flag for return value */
+ if (adelta>tol ||bdelta > tol)
+ { if (*printerr)
+ { printf("gsbound1 error: No convergence for boundary for interim %d; I=%7.0lf",i+1,I[i]);
+ if (bdelta>tol) printf("\n last 2 upper boundary values: %lf %lf\n",btem,btem2);
+ if (adelta>tol) printf("\n last 2 lower boundary values: %lf %lf\n",atem,atem);
+ }
+ retval[0]=1;
+ return;
+ }
+ if (i<nanal-1)
+ { m2=gridpts(r,0.,a[i],b[i],z2,w2);
+ hupdate(0.,w2,m1,I[i-1],z1,h,m2,I[i],z2,h2);
+ m1=m2;
+ tem=z1; z1=z2; z2=tem;
+ tem=w1; w1=w2; w2=tem;
+ tem=h; h=h2; h2=tem;
+ } }
+ retval[0]=0;
+ return;
+}
+
Modified: pkg/gsDesign/src/gsbound1.c
===================================================================
--- pkg/gsDesign/src/gsbound1.c 2011-08-19 05:40:20 UTC (rev 320)
+++ pkg/gsDesign/src/gsbound1.c 2011-08-20 14:40:29 UTC (rev 321)
@@ -1,105 +1,108 @@
-#define DEBUG 0
-#define EXTREMEZ 2000
-#include "R.h"
-#include "Rmath.h"
-/* Group sequential probability computation per Jennison & Turnbull
- Computes upper bound to have input crossing probabilities given fixed input lower bound.
- xnanal- # of possible analyses in the group-sequential designs
- (interims + final)
- xtheta- drift parameter
- I - statistical information available at each analysis
- a - lower cutoff points for z statistic at each analysis (input)
- b - upper cutoff points for z statistic at each analysis (output)
- problo- output vector with probability of rejecting (Z<aj) at
- jth interim analysis, j=1...nanal
- probhi- input vector with probability of rejecting (Z>bj) at
- jth interim analysis, j=1...nanal
- tol - relative change between iterations required to stop for 'convergence'
- xr - controls # of grid points for numerical integration per Jennison & Turnbull
- they recommend r=17 (r=18 is default - a little more accuracy)
- retval- error flag returned; 0 if convergence; 1 indicates error
- printerr- 1 if error messages to be printed - other values suppress printing
-*/
-void gsbound1(int *xnanal,double *xtheta,double *I,double *a,double *b,double *problo,
- double *probhi,double *xtol,int *xr,int *retval,int *printerr)
-{ int i,ii,j,m1,m2,r,nanal;
- double plo=0.,phi,dphi,btem=0.,btem2,rtdeltak,rtIk,rtIkm1,xlo,xhi,theta,mu,tol,bdelta;
-/* note: should allocat zwk & wwk dynamically...*/
- double zwk[1000],wwk[1000],hwk[1000],zwk2[1000],wwk2[1000],hwk2[1000],
- *z1,*z2,*w1,*w2,*h,*h2,*tem;
- void h1(double,int,double *,double,double *, double *);
- void hupdate(double,double *,int,double,double *, double *,
- int,double,double *, double *);
- int gridpts(int,double,double,double,double *, double *);
- r=xr[0]; nanal= xnanal[0]; theta= xtheta[0]; tol=xtol[0];
- if (nanal < 1 || r<1 || r>83)
- { retval[0]=1;
- if (*printerr)
- { Rprintf("gsbound1 error: illegal argument");
- if (nanal<1) Rprintf("; nanal=%d--must be > 0",nanal);
- if (r<1 || r> 83) Rprintf("; r=%d--must be >0 and <84",r);
- Rprintf("\n");
- }
- return;
- }
- rtIk=sqrt(I[0]);
- mu=rtIk*theta; /* mean of normalized statistic at 1st interim */
- problo[0]=pnorm(mu-a[0],0.,1.,0,0); /* probability of crossing lower bound at 1st interim */
- b[0]=qnorm(probhi[0],mu,1,0,0); /* upper bound at 1st interim */
- if (nanal==1) {retval[0]=0; return;}
-/* set up work vectors */
- z1=zwk; w1=wwk; h=hwk;
- z2=zwk2; w2=wwk2; h2=hwk2;
- if (DEBUG) Rprintf("r=%d mu=%lf a[0]=%lf b[0]=%lf\n",r,mu,a[0],b[0]);
- m1=gridpts(r,mu,a[0],b[0],z1,w1);
- h1(theta,m1,w1,I[0],z1,h);
- /* use Newton-Raphson to find subsequent interim analysis cutpoints */
- retval[0]=0;
- for(i=1;i<nanal;i++)
- { rtIkm1=rtIk; rtIk=sqrt(I[i]); mu=rtIk*theta; rtdeltak=sqrt(I[i]-I[i-1]);
- btem2=qnorm(probhi[i],mu,1.,0,0); bdelta=1.; j=0;
- while((bdelta>tol) && j++ < 20)
- { phi=0.; dphi=0.; plo=0.;
- btem=btem2;
- if (DEBUG) Rprintf("i=%d m1=%d\n",i,m1);
- /* compute probability of crossing boundaries & their derivatives */
- for(ii=0;ii<=m1;ii++)
- { xhi=(z1[ii]*rtIkm1-btem*rtIk+theta*(I[i]-I[i-1]))/rtdeltak;
- phi += pnorm(xhi,0.,1.,1,0)*h[ii];
- xlo=(z1[ii]*rtIkm1-a[i]*rtIk+theta*(I[i]-I[i-1]))/rtdeltak;
- plo += pnorm(xlo,0.,1.,0,0)*h[ii];
- dphi-=h[ii]*exp(-xhi*xhi/2)/2.506628275*rtIk/rtdeltak;
- if (DEBUG) Rprintf("m1=%d ii=%d xhi=%lf phi=%lf xlo=%lf plo=%lf dphi=%lf\n",m1,ii,xhi,phi,xlo,plo,dphi);
- }
- /* use 1st order Taylor's series to update boundaries */
- /* maximum allowed change is 1 */
- /* maximum value allowed is EXTREMEZ */
- if (DEBUG)
- Rprintf("i=%2d j=%2d plo=%lf btem=%lf phi=%lf dphi=%lf\n",i,j,plo,btem,phi,dphi);
- bdelta=probhi[i]-phi;
- if (bdelta<dphi) btem2=btem+1.;
- else if (bdelta > -dphi) btem2=btem-1.;
- else btem2=btem+(probhi[i]-phi)/dphi;
- if (btem2>EXTREMEZ) btem2=EXTREMEZ;
- else if (btem2< -EXTREMEZ) btem2= -EXTREMEZ;
- bdelta=btem2-btem; if (bdelta<0) bdelta= -bdelta;
- }
- b[i]=btem;
- problo[i]=plo;
- /* if convergence did not occur, set flag for return value */
- if (bdelta > tol)
- { if (*printerr) printf("gsbound1 error: No convergence for boundary for interim %d; I=%7.0lf; last 2 upper boundary values: %lf %lf\n",
- i+1,I[i],btem,btem2);
- retval[0]=1;
- }
- if (i<nanal-1)
- { m2=gridpts(r,mu,a[i],b[i],z2,w2);
- hupdate(theta,w2,m1,I[i-1],z1,h,m2,I[i],z2,h2);
- m1=m2;
- tem=z1; z1=z2; z2=tem;
- tem=w1; w1=w2; w2=tem;
- tem=h; h=h2; h2=tem;
- } }
- return;
-}
-
+#define DEBUG 0
+#define EXTREMEZ 20
+#define MAXR 83
+#include "R.h"
+#include "Rmath.h"
+/* Group sequential probability computation per Jennison & Turnbull
+ Computes upper bound to have input crossing probabilities given fixed input lower bound.
+ xnanal- # of possible analyses in the group-sequential designs
+ (interims + final)
+ xtheta- drift parameter
+ I - statistical information available at each analysis
+ a - lower cutoff points for z statistic at each analysis (input)
+ b - upper cutoff points for z statistic at each analysis (output)
+ problo- output vector with probability of rejecting (Z<aj) at
+ jth interim analysis, j=1...nanal
+ probhi- input vector with probability of rejecting (Z>bj) at
+ jth interim analysis, j=1...nanal
+ tol - relative change between iterations required to stop for 'convergence'
+ xr - controls # of grid points for numerical integration per Jennison & Turnbull
+ they recommend r=17 (r=18 is default - a little more accuracy)
+ retval- error flag returned; 0 if convergence; 1 indicates error
+ printerr- 1 if error messages to be printed - other values suppress printing
+*/
+void gsbound1(int *xnanal,double *xtheta,double *I,double *a,double *b,double *problo,
+ double *probhi,double *xtol,int *xr,int *retval,int *printerr)
+{ int i,ii,j,m1,m2,r,nanal;
+ double plo=0.,phi,dphi,btem=0.,btem2,rtdeltak,rtIk,rtIkm1,xlo,xhi,theta,mu,tol,bdelta;
+/* note: should allocat zwk & wwk dynamically...*/
+ double zwk[1000],wwk[1000],hwk[1000],zwk2[1000],wwk2[1000],hwk2[1000],
+ *z1,*z2,*w1,*w2,*h,*h2,*tem;
+ void h1(double,int,double *,double,double *, double *);
+ void hupdate(double,double *,int,double,double *, double *,
+ int,double,double *, double *);
+ int gridpts(int,double,double,double,double *, double *);
+ r=xr[0]; nanal= xnanal[0]; theta= xtheta[0]; tol=xtol[0];
+ if (nanal < 1 || r<1 || r>MAXR)
+ { retval[0]=1;
+ if (*printerr)
+ { Rprintf("gsbound1 error: illegal argument");
+ if (nanal<1) Rprintf("; nanal=%d--must be > 0",nanal);
+ if (r<1 || r> MAXR) Rprintf("; r=%d--must be >0 and <84",r);
+ Rprintf("\n");
+ }
+ return;
+ }
+ rtIk=sqrt(I[0]);
+ mu=rtIk*theta; /* mean of normalized statistic at 1st interim */
+ problo[0]=pnorm(mu-a[0],0.,1.,0,0); /* probability of crossing lower bound at 1st interim */
+ if (probhi[0] <= 0.) b[0]=EXTREMEZ;
+ else b[0]=qnorm(probhi[0],mu,1,0,0); /* upper bound at 1st interim */
+ if (nanal==1) {retval[0]=0; return;}
+/* set up work vectors */
+ z1=zwk; w1=wwk; h=hwk;
+ z2=zwk2; w2=wwk2; h2=hwk2;
+ if (DEBUG) Rprintf("r=%d mu=%lf a[0]=%lf b[0]=%lf\n",r,mu,a[0],b[0]);
+ m1=gridpts(r,mu,a[0],b[0],z1,w1);
+ h1(theta,m1,w1,I[0],z1,h);
+ /* use Newton-Raphson to find subsequent interim analysis cutpoints */
+ retval[0]=0;
+ for(i=1;i<nanal;i++)
+ { rtIkm1=rtIk; rtIk=sqrt(I[i]); mu=rtIk*theta; rtdeltak=sqrt(I[i]-I[i-1]);
+ if (probhi[i] <= 0.) btem2=EXTREMEZ;
+ else btem2=qnorm(probhi[i],mu,1.,0,0); bdelta=1.; j=0;
+ while((bdelta>tol) && j++ < 20)
+ { phi=0.; dphi=0.; plo=0.;
+ btem=btem2;
+ if (DEBUG) Rprintf("i=%d m1=%d\n",i,m1);
+ /* compute probability of crossing boundaries & their derivatives */
+ for(ii=0;ii<=m1;ii++)
+ { xhi=(z1[ii]*rtIkm1-btem*rtIk+theta*(I[i]-I[i-1]))/rtdeltak;
+ phi += pnorm(xhi,0.,1.,1,0)*h[ii];
+ xlo=(z1[ii]*rtIkm1-a[i]*rtIk+theta*(I[i]-I[i-1]))/rtdeltak;
+ plo += pnorm(xlo,0.,1.,0,0)*h[ii];
+ dphi-=h[ii]*exp(-xhi*xhi/2)/2.506628275*rtIk/rtdeltak;
+ if (DEBUG) Rprintf("m1=%d ii=%d xhi=%lf phi=%lf xlo=%lf plo=%lf dphi=%lf\n",m1,ii,xhi,phi,xlo,plo,dphi);
+ }
+ /* use 1st order Taylor's series to update boundaries */
+ /* maximum allowed change is 1 */
+ /* maximum value allowed is EXTREMEZ */
+ if (DEBUG)
+ Rprintf("i=%2d j=%2d plo=%lf btem=%lf phi=%lf dphi=%lf\n",i,j,plo,btem,phi,dphi);
+ bdelta=probhi[i]-phi;
+ if (bdelta<dphi) btem2=btem+1.;
+ else if (bdelta > -dphi) btem2=btem-1.;
+ else btem2=btem+(probhi[i]-phi)/dphi;
+ if (btem2>EXTREMEZ) btem2=EXTREMEZ;
+ else if (btem2< -EXTREMEZ) btem2= -EXTREMEZ;
+ bdelta=btem2-btem; if (bdelta<0) bdelta= -bdelta;
+ }
+ b[i]=btem;
+ problo[i]=plo;
+ /* if convergence did not occur, set flag for return value */
+ if (bdelta > tol)
+ { if (*printerr) printf("gsbound1 error: No convergence for boundary for interim %d; I=%7.0lf; last 2 upper boundary values: %lf %lf\n",
+ i+1,I[i],btem,btem2);
+ retval[0]=1;
+ }
+ if (i<nanal-1)
+ { m2=gridpts(r,mu,a[i],b[i],z2,w2);
+ hupdate(theta,w2,m1,I[i-1],z1,h,m2,I[i],z2,h2);
+ m1=m2;
+ tem=z1; z1=z2; z2=tem;
+ tem=w1; w1=w2; w2=tem;
+ tem=h; h=h2; h2=tem;
+ } }
+ return;
+}
+
More information about the Gsdesign-commits
mailing list