**********************************************************************; **********************************************************************; * GLUMIP20.IML GLUMIP V2.0 Modified 02/25/08 *; * *; *07/07/03: *; * Quantile Transform for Power and Type I error rate calculations *; * *; * Modified version of GLUMIP.IML that includes the new exact *; * calculation for the distribution of the internal pilot test *; * statistic *; * *; * This SAS code, used inside PROC IML via an INCLUDE statement, *; * computes power under an internal pilot design for a specified *; * method (combination of a sample size re-estimation rule and a *; * test) *; * *; * This code follows the notation introduced in Coffey and Muller, *; * Biometrics, 2001. *; * *; * NOTE: To Find type I error rate use the fact that it is the power *; * under the null and specify BETA_ALT as the matrix of primary *; * parameters under the null. *; * *; * WARNING: This code assumes full rank ESSENCEX and C ! *; * *; **********************************************************************; * REQUIRED INPUTS: *; * *; * 1) Design Parameters: *; * ESSENCEX = Essence design matrix *; * ALPHAT = Target type I error rate *; * POWERT = Target power *; * C = Between subject contrast matrix *; * BETA_PLN = Vector of primary parameters for power calculations *; * used for initial planning *; * SIGMA0 = Variance estimate used for planning (**NOT** st dev) *; * *; * 2) Sample Size Allocation: *; * N1 = Internal pilot sample size *; * NPLUSMIN = Minimum size of final sample *; * *; * 3) Fixed Parameters *; * GAMLIST = vector of ratios of true to pre-planned variance *; * BETA_ALT = vector of assumed true primary parameters *; * (specify under the null for type I error rate) *; * *; * 4) Specify Method to be used *; * RULE = Sample Size Re-estimation Rule (0,1,2) *; * TEST = Testing Procedure (0,1,2,3,4) *; * *; * OPTIONAL INPUTS: *; * WEIGHTS = Vector whose elements correspond to weights for each *; * row of ESSENCEX per replication. Default assumption is *; * a balanced design *; * (Sum of elements equals replication size) *; * NPLUSMAX = Maximum size of final sample. Default is no finite *; * maximum *; * ROUND = helps define ROUNDOFF=ALPHA#10##(-ROUND) *; * controls precision in computing Pr{sample size} (1-10).*; * Default = 3. *; * *; **********************************************************************; * OUTPUTS: *; * *; * _IPCALCS: Matrix with one row for each GAMMA value requested *; * and columns corresponding to: *; * (1) ALPHAT *; * (2) CRITVAL *; * (3) POWERT *; * (4) GAMMA *; * (5) N1 *; * (6) NPLUSMIN *; * (7) NPLUSMAX *; * (8) RULE *; * (9) TEST *; * (10) EXP_N: Expected sample size under the IP design *; * (11) POWER: Power under the IP design *; * *; * _IPNAMES: Matrix of column names corresponding to _IPCALCS *; **********************************************************************; **********************************************************************; ***************************************************************; ***************************************************************; * GLUMIP.IML *; * Define GLUMIP command. *; * User will only have to type RUN GLUMIP. The input matrices *; * will not have to be listed in the RUN statement. *; ***************************************************************; ***************************************************************; START GLUMIP; RUN _GLUMIP(_IPCALCS,_IPNAMES); FINISH GLUMIP; ***************************************************************; ***************************************************************; * MAXTS.IML *; * Define MAXTS command. *; * User will only have to type RUN MAXTS. The input matrices *; * will not have to be listed in the RUN statement. *; ***************************************************************; ***************************************************************; START MAXTS; RUN _MAXTSB(_MAXTS,_MAXTSNM); FINISH MAXTS; ****************************************************************; ****************************************************************; * FINDADJ.IML *; * Define FINDADJ command. *; * User will only have to type RUN FINDADJ. The input matrices *; * will not have to be listed in the RUN statement. *; ****************************************************************; ****************************************************************; START FINDADJ; RUN _FINDADJB(_FINDADJ,_FINDADJNM); FINISH FINDADJ; ***************************************************************; ***************************************************************; * N2CALC.IML *; * Define N2CALC command. *; * User will only have to type RUN N2CALC. The input matrices *; * will not have to be listed in the RUN statement. *; ***************************************************************; ***************************************************************; START N2CALC; RUN _N2CALC(_N2CALCS,_N2CALCNM); FINISH N2CALC; ***************************************************************; ***************************************************************; * _N2CALC.IML *; * Main Code *; ***************************************************************; ***************************************************************; START _N2CALC(_N2CALCS,_N2CALCNM) GLOBAL(SIGHAT1,ALPHAT,POWERT,BETA_PLN,_ZERO_,ESSENCEX, C,N1,NPLUSMIN,NPLUSMAX,WEIGHTS,RULE); ******************************************; *Define default values if not specified *; ******************************************; IF NROW(WEIGHTS)=0 THEN _WEIGHTS=J(NROW(ESSENCEX),1,1); ELSE IF NCOL(WEIGHTS)>1 THEN _WEIGHTS=WEIGHTS`; ELSE _WEIGHTS=WEIGHTS; IF NROW(NPLUSMAX)=0 THEN DO; _NPMAX=NPLUSMIN+1000#SUM(_WEIGHTS); _NPMAX_=.i; END; ELSE IF NPLUSMAX=. THEN DO; _NPMAX=NPLUSMIN+1000#SUM(_WEIGHTS); _NPMAX_=.i; END; ELSE DO; _NPMAX=NPLUSMAX; _NPMAX_=NPLUSMAX; END; *********************************; *Create initial conditions *; *********************************; R = NCOL(ESSENCEX); *Assume full rank design; A = NROW(C); *Numerator df; NU1 = N1-R; *Error df for internal pilot sample; REPSIZE = SUM(_WEIGHTS); *Replication size; IF RULE=1 THEN DO; FCRIT = FINV(1-ALPHAT,A,NU1); NONCENT = FNONCT(FCRIT,A,NU1,1-POWERT); END; THTA_POP = C*BETA_PLN; *True value of THETA in population; XPX_WESS = ESSENCEX`*(DIAG(_WEIGHTS)/SUM(_WEIGHTS))*ESSENCEX; INVXPXES = SOLVE(XPX_WESS,I(R)); *This will produce error if LTFR; SSH_WESS = (THTA_POP`)*INV(C*(INVXPXES)*C`)*THTA_POP; * CAUTION: FNONCT function may blow up at extremes; DONE=0; DO NPLUS=NPLUSMIN TO _NPMAX BY REPSIZE UNTIL (DONE=1); NUPLUS = NPLUS - R; N2CALC = NPLUS-N1; IF (NPLUS=_NPMAX) THEN DO; IF RULE=0 THEN FCRIT = FINV(1-ALPHAT,A,NUPLUS); IF RULE=2 THEN FCRIT = FINV(1-ALPHAT,A,N2CALC); DONE=1; END; ELSE DO; IF RULE=0 THEN DO; FCRIT = FINV(1-ALPHAT,A,NUPLUS); NONCENT = FNONCT(FCRIT,A,NUPLUS,1-POWERT); END; ELSE IF RULE=2 THEN DO; FCRIT = FINV(1-ALPHAT,A,N2CALC); NONCENT = FNONCT(FCRIT,A,N2CALC,1-POWERT); END; SIGMANP = (NPLUS#SSH_WESS)/NONCENT; IF SIGHAT1<=SIGMANP THEN DONE=1; END; END; NPLUSCALC=N1+N2CALC; ESTNONCENT=(NPLUSCALC#SSH_WESS)/SIGHAT1; IF RULE=0 THEN POWPROJ=1-PROBF(FCRIT,A,NUPLUS,ESTNONCENT); IF RULE=1 THEN POWPROJ=1-PROBF(FCRIT,A,NU1,ESTNONCENT); IF RULE=2 THEN POWPROJ=1-PROBF(FCRIT,A,N2CALC,ESTNONCENT); _N2CALCS=SIGHAT1||N1||RULE||N2CALC||NPLUSCALC||POWPROJ; _N2CALCNM= {"SIGMAHAT1" "N1" "RULE" "N2CALC" "NPLUSCALC" "Projected Power"}; FINISH _N2CALC; ***************************************************************; ***************************************************************; * _MAXTSB.IML *; * Main Code *; ***************************************************************; ***************************************************************; START _MAXTSB(_MAXTS,_MAXTSNM) GLOBAL(ALPHAT,POWERT,BETA_PLN,ROUND,_ZERO_,ESSENCEX, C,N1,NPLUSMIN,NPLUSMAX,SIGMA0,WEIGHTS,TEST); ******************************************; *Define default values if not specified *; ******************************************; IF NROW(TEST)>0 THEN IF ( (TEST=1)|(TEST=2)|(TEST=3)|(TEST=4) ) THEN PRINT "WARNING: MAXTS only computes the maximum type I error rate for TEST=0"; IF NROW(_ZERO_)=0 THEN _ZERO_=1E-9; *Taken as numeric zero; IF NROW(WEIGHTS)=0 THEN _WEIGHTS=J(NROW(ESSENCEX),1,1); ELSE IF NCOL(WEIGHTS)>1 THEN _WEIGHTS=WEIGHTS`; ELSE _WEIGHTS=WEIGHTS; IF NROW(NPLUSMAX)=0 THEN DO; _NPMAX=NPLUSMIN+1000#SUM(_WEIGHTS); _NPMAX_=.i; END; ELSE IF NPLUSMAX=. THEN DO; _NPMAX=NPLUSMIN+1000#SUM(_WEIGHTS); _NPMAX_=.i; END; ELSE DO; _NPMAX=NPLUSMAX; _NPMAX_=NPLUSMAX; END; IF NROW(ROUND)=0 THEN _ROUND=3; ELSE _ROUND=ROUND; ALPH_ADJ=ALPHAT; _MAXTS=_MAXTS(ALPHAT,ALPH_ADJ,POWERT,BETA_PLN,_ROUND,ESSENCEX, C,N1,NPLUSMIN,_NPMAX,SIGMA0,_Weights); _MAXTSNM = {"GAMMA(MAX)" " MAX TYPE I ERR. RATE" " RATIO: MAX/ALPHAT"}; FINISH _MAXTSB; ***************************************************************; ***************************************************************; * _FINDADJB.IML *; * Main Code *; ***************************************************************; ***************************************************************; START _FINDADJB(_FINDADJ,_FINDADJNM) GLOBAL(ALPHAT,POWERT,BETA_PLN,ROUND,_ZERO_,ESSENCEX, C,N1,NPLUSMIN,NPLUSMAX,SIGMA0,WEIGHTS,TEST); ******************************************; *Define default values if not specified *; ******************************************; IF NROW(TEST)>0 THEN IF ( (TEST=1)|(TEST=2)|(TEST=3)|(TEST=4) ) THEN PRINT "WARNING: MAXTS only computes the maximum type I error rate for TEST=0"; IF NROW(_ZERO_)=0 THEN _ZERO_=1E-9; *Taken as numeric zero; IF NROW(WEIGHTS)=0 THEN _WEIGHTS=J(NROW(ESSENCEX),1,1); ELSE IF NCOL(WEIGHTS)>1 THEN _WEIGHTS=WEIGHTS`; ELSE _WEIGHTS=WEIGHTS; IF NROW(NPLUSMAX)=0 THEN DO; _NPMAX=NPLUSMIN+1000#SUM(_WEIGHTS); _NPMAX_=.i; END; ELSE IF NPLUSMAX=. THEN DO; _NPMAX=NPLUSMIN+1000#SUM(_WEIGHTS); _NPMAX_=.i; END; ELSE DO; _NPMAX=NPLUSMAX; _NPMAX_=NPLUSMAX; END; IF NROW(ROUND)=0 THEN _ROUND=3; ELSE _ROUND=ROUND; _FINDADJ=_FINDADJ(ALPHAT,POWERT,BETA_PLN,_ROUND, ESSENCEX,C,N1,NPLUSMIN,_NPMAX, SIGMA0,_WEIGHTS); _FINDADJNM = {"Adjusted Alpha"}; FINISH _FINDADJB; ***************************************************************; ***************************************************************; * _GLUMIP.IML *; * Main Code *; ***************************************************************; ***************************************************************; START _GLUMIP(_IPCALCS,_IPNAMES) GLOBAL (GLOBAL1,ESSENCEX,ALPHAT,POWERT,C,BETA_PLN,SIGMA0, N1,NPLUSMIN,GAMLIST,BETA_ALT,RULE,TEST, WEIGHTS,NPLUSMAX,ROUND,_ZERO_); **********************************************; *1. checks for presence of required inputs *; **********************************************; *Could use CHECKER in POWERLIB to check existence, type, dimensions; NAMES={ESSENCEX ALPHAT POWERT C BETA_PLN SIGMA0 N1 NPLUSMIN GAMLIST BETA_ALT RULE TEST}`; ROWS=NROW(ESSENCEX)//NROW(ALPHAT)//NROW(POWERT)//NROW(C)// NROW(BETA_PLN)//NROW(SIGMA0)//NROW(N1)//NROW(NPLUSMIN)// NROW(GAMLIST)//NROW(BETA_ALT)//NROW(RULE)//NROW(TEST); IF MIN(ROWS)=0 THEN DO; PRINT "Undefined required input in GLUMIP-- # rows = " ROWS [ROWNAME=NAMES]; STOP; END; ******************************************; *2. Delete previous versions of _IPCALCS *; ******************************************; IF NROW(_IPCALCS)>0 THEN FREE _IPCALCS; *********************************************************; *3. Create ROW matrix (list) for GAMLIST if user did not*; *********************************************************; IF NROW(GAMLIST)>1 THEN _GAMLIST=GAMLIST`; ELSE _GAMLIST=GAMLIST; ********************************************; *4. Define default values if not specified *; ********************************************; IF NROW(_ZERO_)=0 THEN _ZERO_=1E-9; *Taken as numeric zero; IF NROW(WEIGHTS)=0 THEN _WEIGHTS=J(NROW(ESSENCEX),1,1); ELSE IF NCOL(WEIGHTS)>1 THEN _WEIGHTS=WEIGHTS`; ELSE _WEIGHTS=WEIGHTS; IF NROW(NPLUSMAX)=0 THEN DO; _NPMAX=NPLUSMIN+1000#SUM(_WEIGHTS); _NPMAX_=.i; END; ELSE IF NPLUSMAX=. THEN DO; _NPMAX=NPLUSMIN+1000#SUM(_WEIGHTS); _NPMAX_=.i; END; ELSE DO; _NPMAX=NPLUSMAX; _NPMAX_=NPLUSMAX; END; IF NROW(ROUND)=0 THEN _ROUND=3; ELSE _ROUND=ROUND; ************************************; *5. checks on validity of inputs *; ************************************; BAD="----------None----------"; IF NROW(ESSENCEX)=1) THEN BAD=BAD// {"(ALPHAT<=0)|(ALPHAT>=1)"}; IF (POWERT=1) THEN BAD=BAD// {"(POWERT=1)"}; IF NCOL(C)^=NROW(BETA_PLN) THEN BAD=BAD// {"C*BETA_PLN does not conform"}; IF NCOL(C)^=NROW(BETA_ALT) THEN BAD=BAD// {"C*BETA_ALT does not conform"}; IF NROW(BETA_PLN)^=NCOL(ESSENCEX) THEN BAD=BAD// {"ESSENCEX*BETA_PLN does not conform"}; IF (SIGMA0<=0) THEN BAD=BAD// {"SIGMA0<=0"}; IF ROUND(N1,1)^=N1 THEN BAD=BAD// {"N1 not integer"}; IF MOD(N1,SUM(_WEIGHTS))^=0 THEN BAD=BAD// {"N1 does not agree with replication size"}; IF ROUND(NPLUSMIN,1)^=NPLUSMIN THEN BAD=BAD// {"NPLUSMIN not integer"}; IF ROUND(_NPMAX,1)^=_NPMAX THEN BAD=BAD// {"NPLUSMAX not integer"}; IF (NPLUSMIN < N1) THEN BAD=BAD// {"NPLUSMIN < N1"}; IF (NPLUSMIN = N1)&(TEST=4) THEN BAD=BAD// {"NPLUSMIN must be > N1 for TEST=4.", "For your given inputs, TEST=1 could be appropriate."}; ELSE IF (TEST=4) THEN PRINT "WARNING: TEST=4 is only appropriate when NPLUSMIN = N_not, the" , "originally planned sample size. Otherwise interpret results with caution."; IF (NPLUSMIN = N1)&(TEST=2) THEN BAD=BAD// {"NPLUSMIN must be > N1 for TEST=2. " , "TEST=2 may still not run for NPLUSMIN close to N1"}; ELSE IF (TEST=2) & (NPLUSMIN <= (N1+SUM(_WEIGHTS)) ) THEN PRINT "WARNING: TEST=2 may not run for NPLUSMIN close to N1"; IF MOD(NPLUSMIN,SUM(_WEIGHTS))^=0 THEN BAD=BAD// {"NPLUSMIN does not agree with replication size"}; IF (_NPMAX < NPLUSMIN) THEN BAD=BAD// {"NPLUSMAX < NPLUSMIN"}; IF (_NPMAX < N1) THEN BAD=BAD// {"NPLUSMAX < N1"}; IF MOD(_NPMAX,SUM(_WEIGHTS))^=0 THEN BAD=BAD// {"NPLUSMAX does not agree with replication size"}; IF (MIN(_GAMLIST)<=0) THEN BAD=BAD// {"(GAMMA value(s) <=0)"}; IF NROW(BETA_ALT)^=NCOL(ESSENCEX) THEN BAD=BAD// {"ESSENCEX*BETA_ALT does not conform"}; IF (RULE^=0)&(RULE^=1)&(RULE^=2) THEN BAD=BAD// {"Undefined Rule"}; IF (TEST^=0)&(TEST^=1)&(TEST^=2)&(TEST^=3)&(TEST^=4) THEN BAD=BAD// {"Undefined Test"}; IF (_ROUND < 1 | _ROUND > 10) THEN BAD=BAD// {"User Specified ROUND must be between 1 and 10"}; IF NROW(BAD)>1 THEN DO; BAD=BAD[2:NROW(BAD),]; PRINT BAD; STOP; END; ************************************; *6. Create inputs to other modules *; ************************************; _ESSENCEX = ESSENCEX; _ALPHAT = ALPHAT; _POWERT = POWERT; _C = C; _BETAPLN = BETA_PLN; _SIGMA0 = SIGMA0; _N1 = N1; _NPMIN = NPLUSMIN; _BETAALT = BETA_ALT; _RULE = RULE; ***********************************************************; *7. Find Adjusted Critical Value if Using Bounding Method *; ***********************************************************; IF TEST=3 THEN _ALPHA = _FINDADJ(_ALPHAT,_POWERT,_BETAPLN,_ROUND, _ESSENCEX,_C,_N1,_NPMIN,_NPMAX, _SIGMA0,_WEIGHTS); ELSE _ALPHA=_ALPHAT; ***********************************************************; *8. Perform Calculations *; ***********************************************************; GAMMAMAX=MAX(_GAMLIST); _ISTOP=0; *Indicater for NDIST to use original stopping rule; RUN _NDIST(GAMMAMAX,_ALPHAT,_POWERT,_ESSENCEX,_BETAPLN,_C,_N1, _ROUND,_NPMIN,_RULE,_SIGMA0,_NPMAX,_WEIGHTS,_TABLE,_ISTOP); _TABLEMAX=_TABLE; NTABLE=NROW(_TABLEMAX); _NPMAX=_TABLEMAX[NTABLE,2]; _ISTOP=1; *Indicater for NDIST stopping rule new; DO K=1 TO NCOL(_GAMLIST); *Solves NDIST and GETPOW for each gamma; FREE _TABLE; _GAMMA = _GAMLIST[1,K]; IF _GAMMA ^= GAMMAMAX THEN RUN _NDIST(_GAMMA,_ALPHAT,_POWERT,_ESSENCEX,_BETAPLN,_C, _N1,_ROUND, _NPMIN,_RULE, _SIGMA0,_NPMAX,_WEIGHTS, _TABLE,_ISTOP); ELSE _TABLE=_TABLEMAX; EXP_N = ROUND(_TABLE[,2]`*_TABLE[,3],.1); IF TEST=3 THEN _TEST=0; ELSE _TEST=TEST; POWER=_GETPOW(_GAMMA,_BETAALT,_ALPHA,_TABLE,_ESSENCEX, _C,_N1,_NPMIN,_TEST,_SIGMA0,_WEIGHTS); _IPCALCS = _IPCALCS//(_ALPHAT||_ALPHA||_POWERT||_GAMMA||_N1||_NPMIN||_NPMAX_||_RULE||TEST||EXP_N||POWER); END; _IPNAMES = {"ALPHAT" "CRITVAL" "POWERT" "GAMMA" "N1" "NPLUSMIN" "NPLUSMAX" "RULE" "TEST" "E(N)" "POWER"}; FINISH _GLUMIP; **********************************************************************; **********************************************************************; * _FCTS1.IML *; * Determines the value of the first integrand for computing the *; * distribution of the internal pilot test statistic *; * *Reference: implment equation 2.9 in PAPER0921.WXP *; * *; **********************************************************************; * REQUIRED INPUTS: *; * P = Point at which integrand is to be evaluated *; * *; * (GLOBALS) *; * ALB1 = Constant multiplier of first chi-square R.V. *; * ANC1 = Noncentrality of first chi-square R.V. *; * A = Degrees of freedom for first chi-square R.V. *; * N2 = Degrees of freedom for second chi-square R.V. *; * NU1 = Degrees of freedom for internal pilot sample *; * TRUNC = Vector corresponding to the truncation region *; * *; **********************************************************************; **********************************************************************; START _FCTS1(P) GLOBAL (GLOBAL1,_ZERO_); *Here 0<=P<=1 is a probability from a quantile transform; ALB1 = GLOBAL1[1]; ANC1 = GLOBAL1[2]; A = GLOBAL1[3]; N2 = GLOBAL1[4]; NU1 = GLOBAL1[5]; TRUNC = GLOBAL1[6:7]; IF P <=_ZERO_ THEN RETURN({0}); IF ({1}-P)<=_ZERO_ THEN RETURN({0}); X=CINV(P,N2+NU1); PROBC=CDF("CHISQ",(X/ALB1) ,A ,ANC1); B=TRUNC[2]/X; IF B<=_ZERO_ THEN RETURN({0}); IF ({1}-B)<=_ZERO_ THEN PROBB={1}; ELSE PROBB = CDF("BETA", B ,NU1/2,N2/2); PRODUCT = PROBC#PROBB; RETURN(PRODUCT); FINISH _FCTS1; **********************************************************************; * _FCTS2.IML *; * Determines the value of the second integrand for computing the *; * distribution of the internal pilot test statistic *; * Reference: implment equation 2.9 in PAPER0921.WXP *; * *; **********************************************************************; * REQUIRED INPUTS: *; * P = Point at which integrand is to be evaluated *; * *; * (GLOBALS) *; * ALB1 = Constant multiplier of first chi-square R.V. *; * ANC1 = Noncentrality of first chi-square R.V. *; * A = Degrees of freedom for first chi-square R.V. *; * N2 = Degrees of freedom for second chi-square R.V. *; * NU1 = Degrees of freedom for internal pilot sample *; * TRUNC = Vector corresponding to the truncation region *; * *; **********************************************************************; **********************************************************************; START _FCTS2(P) GLOBAL (GLOBAL1,_ZERO_); *Here 0<=P<=1 is a probability from a quantile transform; ALB1 = GLOBAL1[1]; ANC1 = GLOBAL1[2]; A = GLOBAL1[3]; N2 = GLOBAL1[4]; NU1 = GLOBAL1[5]; TRUNC = GLOBAL1[6:7]; IF P <=_ZERO_ THEN RETURN({0}); IF ({1}-P)<=_ZERO_ THEN RETURN({0}); X=CINV(P,N2+NU1); PROBC=CDF("CHISQ",(X/ALB1) ,A ,ANC1); B=TRUNC[1]/X; *only line of code dif from _FCTS1; IF B<=_ZERO_ THEN RETURN({0}); IF ({1}-B)<=_ZERO_ THEN PROBB={1}; ELSE PROBB = CDF("BETA", B ,NU1/2,N2/2); PRODUCT = PROBC#PROBB; RETURN(PRODUCT); FINISH _FCTS2; **********************************************************************; **********************************************************************; * _FCTT.IML *; * Determines the value of the integrand for computing the *; * distribution of the internal pilot test statistic when the final *; * sample size equals the internal pilot sample size or *; * when TEST=1 *; * Reference: Equation (2.24) in Coffey and Muller, 1999 *; * *; **********************************************************************; * REQUIRED INPUTS: *; * S = Point at which integrand is to be evaluated *; * *; * (GLOBALS) *; * ALB1 = Constant multiplier of first chi-square R.V. *; * ANC1 = Noncentrality of first chi-square R.V. *; * A = Degrees of freedom for first chi-square R.V. *; * NU1 = Degrees of freedom for internal pilot sample *; * *; **********************************************************************; **********************************************************************; START _FCTT(S) GLOBAL (GLOBAL1,_ZERO_); ALB1 = GLOBAL1[1]; ANC1 = GLOBAL1[2]; A = GLOBAL1[3]; NU1 = GLOBAL1[5]; IF S<_ZERO_ THEN DO; PRODUCT=0; RETURN(PRODUCT); END; X1=CDF("CHISQUARE",S/ALB1,A,ANC1); X2=PDF("CHISQUARE",S,NU1,0); PRODUCT=X1#X2; RETURN(PRODUCT); FINISH _FCTT; **********************************************************************; **********************************************************************; * _QUAD2D01.IML *; * 2D numerical integration OF quantile transformed integrand for *; * univariate IP *; * Uses NDIMRI, p506, Davis and Rabinowitz (1984) Methods of *; * Numerical Integration *; * *; **********************************************************************; * REQUIRED INPUTS: *; * H = 2x1 fractions of form H[]=1/I for a positive integer I; *; * MAXITER = 1x1 maximum # of iterations allowed; *; * ERROR = 1x1 maximum error allowed; *; * *; * (GLOBALS) *; * D = FHI - FLOW = Pr{N+ = n+} *; * FLOW = CDF("CHISQ", Q1 , NU1 ) *; * NU1 = n1 -rank(X) *; * W = ( n(+min) - n1 ) / N2 *; * N2 *; * A hypoth df = numerator df *; * ANC1 omega+ , noncentrality, n+#omegaess *; * ALBMIN = ( n(+min) - rank(X) ) / (A# FCRIT) *; * Q1 *; * *; * OUTPUTS *; *AREA 1x1 returned value of integral *; *FLAG 1x1 -1=>failed to converge, 0=>invalid input, 1=>converged *; **********************************************************************; **********************************************************************; START _QUAD2D01(H,MAXITER,ERROR,AREA,FLAG) GLOBAL(_ZERO_,GLOBAL2); *Grab IP variables; D =GLOBAL2[1]; *= FHI - FLOW = Pr{N+ = n+}; FLOW =GLOBAL2[2]; *= CDF("CHISQ", Q1 , NU1 ); NU1 =GLOBAL2[3]; *= n1 -rank(X) ; W =GLOBAL2[4]; *= ( n(+min) - n1 ) / N2 ; N2 =GLOBAL2[5]; A =GLOBAL2[6]; *hypoth df = numerator df; ANC1 =GLOBAL2[7]; *omega+ , noncentrality, n+#omegaess; ALBMIN =GLOBAL2[8]; *= ( n(+min) - rank(X) ) / (A# FCRIT) ; Q1 =GLOBAL2[9]; UPPERMAX=15; *upperbound limit on maximum iterations; AL=2; *controls refinement of grid; P =J(2,1,.); FLAG=0; AREA=.; IF (MAXITER<{1}) | (MIN(H)<={0}) | (MAX(H)>{1}) | (ERROR<_ZERO_) THEN DO; PRINT "Invalid input causes STOP in _QUAD2D01" H ERROR MAXITER _ZERO_; STOP; RETURN; END; IF MAXITER>UPPERMAX THEN DO; PRINT "Program limits MAXITER to " UPPERMAX; STOP; END; V =J(MAXITER,1,.); AA=J(MAXITER,1,.); K =J(MAXITER,1,.); K[1:2]= {1}//AL ; *was {1,2}; C={1,1}; *ranges of integration; *bounds are unit square; L={1}; L21: U ={0}; G=H/K[L]; NN= FLOOR({1}/G+{.5}); *variables starting in I-N are integer in FORTRAN; DO I2=1 TO NN[2]; P[2]= G[2]#(I2-{.5}); DO I1=1 TO NN[1]; P[1]= G[1]#(I1-{.5}); *********************************************; *Function embedded to improve speed ; *********************************************; IF (Q1<=_ZERO_)&(MAX(P)<=_ZERO_) THEN DO; PROB12={0}; GO TO GOT_PROB; END; IF ({1}-P[1]) <=_ZERO_ THEN DO; PROB12={1}; GO TO GOT_PROB; END; X1=CINV(P[1],N2); PX2=P[2]#D+FLOW; IF ({1}-PX2) <=_ZERO_ THEN DO; PROB12={1}; GO TO GOT_PROB; END; X2=CINV(PX2,NU1); XVALUE= (X1#W + X2) /ALBMIN; PROB12=CDF("CHISQ", XVALUE , A, ANC1 ); GOT_PROB: ; ******************End of function evaluation***************; U=U+PROB12; END; END; U=U#G[#]; V[L]=U; IF L>={2} THEN GO TO L44; AA[1]=V[1]; L=L+1; GO TO L21; L44: EN=K[L]; DO LL=2 TO L; I=L+{1}-LL; V[I]= V[(I+{1})] + ( V[(I+{1})] - V[I] ) / ( ((EN/K[I])##{2}) - {1} ); END; AREA=V[1]; FLAG={1}; IF ABS(AREA-AA[L-{1}]) < ABS(AREA#ERROR) THEN RETURN; FLAG={-1}; IF L=MAXITER THEN RETURN; AA[L]=AREA; L= L + {1}; K[L]= AL#K[L-{1}] ; GO TO L21; FINISH _QUAD2D01; **********************************************************************; **********************************************************************; * _GETPOW.IML *; * Computes power under an internal pilot design *; * (Reference: Coffey and Muller, Biometrics 2001) *; * *; **********************************************************************; * REQUIRED INPUTS: *; * *; * ESSENCEX = Essence design matrix *; * WEIGHTS = Weights for each row of ESSENCEX per replication *; * ALPHAC = Type I error rate used for determining critical value *; * C = Between subject contrast matrix *; * SIGMA0 = Variance estimate used for planning (**NOT** std dev) *; * N1 = Internal pilot sample size *; * NPLUSMIN = Minimum size of final sample *; * TABLE = Q(N+)||N+||PDF||CDF *; * GAMMA = Ratio of true variance to pre-planned variance *; * BETA_ALT = True vector of primary parameters *; * TEST = Testing Procedure (0,1,2,3) *; * *; **********************************************************************; **********************************************************************; START _GETPOW(GAMMA,BETA_ALT,ALPHAC,TABLE, ESSENCEX,C,N1,NPLUSMIN,TEST, SIGMA0,WEIGHTS) GLOBAL (GLOBAL1,GLOBAL2,_ZERO_); *******************************************; * Create initial conditions *; *******************************************; A = NROW(C); R = NCOL(ESSENCEX); NU1 = N1-R; WARN = 0; TRUNCERR= 0; THETA = C*BETA_ALT; XPX_WESS = ESSENCEX`*(DIAG(WEIGHTS)/SUM(WEIGHTS))*ESSENCEX; INVXPXES = SOLVE(XPX_WESS,I(R)); OMEGAESS = ((THETA`)*INV(C*INVXPXES*C`)*(THETA))/(SIGMA0#GAMMA); ***********************************************************; * Define values which remain constant over the looping *; ***********************************************************; IF TEST=1 THEN DO; FCRIT = FINV(1-ALPHAC,A,NU1); ALB1 = NU1/(FCRIT*A); END; * Initialize Variables; Q2 = 0; NPTERMS = NROW(TABLE); *If TABLE has only one row then we effectively have a fixed sample size distribution; IF NPTERMS=1 THEN DO; NPLUS = TABLE[1,2]; ANC1 = NPLUS*OMEGAESS; IF TEST=0 THEN DO; NUPLUS = NPLUS - R; FCRIT = FINV(1-ALPHAC,A,NUPLUS); POWER = SDF("F",FCRIT,A,NUPLUS,ANC1); RETURN(POWER); END; ELSE IF TEST=1 THEN DO; POWER = SDF("F",FCRIT,A,NU1,ANC1); RETURN(POWER); END; ELSE IF TEST=2 THEN DO; N2 = NPLUS - N1; FCRIT = FINV(1-ALPHAC,A,N2); POWER = SDF("F",FCRIT,A,N2,ANC1); RETURN(POWER); END; ELSE IF TEST=4 THEN DO; IF NPLUS=NPLUSMIN THEN DO; NUPLUS = NPLUS - R; FCRIT = FINV(1-ALPHAC,A,NUPLUS); POWER = SDF("F",FCRIT,A,NUPLUS,ANC1); RETURN(POWER); END; END; END; **************************************; * Perform looping over NPLUS values *; **************************************; DO NPCASE = 1 TO NPTERMS; NPLUS = TABLE[NPCASE,2]; ANC1 = NPLUS*OMEGAESS; * Noncentrality Parameter; N2 = NPLUS-N1; IF TEST=0 THEN DO; NUPLUS = NPLUS-R; FCRIT = FINV(1-ALPHAC,A,NUPLUS); ALB1 = NUPLUS/(FCRIT*A); * Eigenvalue of AF; END; IF (TEST^=2) THEN DO; Q1 = Q2; * Lower Integration Value; Q2 = TABLE[NPCASE,1]; * Upper Integration Value; GLOBAL1 = ALB1||ANC1||A||N2||NU1||Q1||Q2; END; **************************************************; * Integrations for Test 0 *; **************************************************; IF TEST=0 THEN DO; IF NPLUS = N1 THEN DO; Q = GLOBAL1[6:7]; CALL QUAD(AREA,"_FCTT",Q,1E-7,,,"NO"); FIRST=AREA; SECOND=0; END; ELSE DO; QUANT1=PROBCHI(Q1,NU1+N2); QUANT2=PROBCHI(Q2,NU1+N2); QD1 = QUANT1||QUANT2||1; QD2 = QUANT1||1; CALL QUAD(FIRST,"_FCTS1",QD1,1E-7,,,"NO"); IF (NPCASE^=1) THEN CALL QUAD(SECOND,"_FCTS2",QD2,1E-7,,,"NO"); ELSE SECOND=0; AREA = SUM(FIRST)-SECOND; END; CUMTERM = AREA/TABLE[NPCASE,3]; OUT_AREA= OUT_AREA//(NPLUS||AREA||CUMTERM||TABLE[NPCASE,3]||SUM(FIRST)||SECOND); END; **************************************************; * Integrations for Test 1 *; **************************************************; IF TEST=1 THEN DO; Q = GLOBAL1[6:7]; CALL QUAD(AREA,"_FCTT",Q,1E-7,,,"NO"); OUT_AREA = OUT_AREA//(NPLUS||AREA); END; **************************************************; * Calculations for Test 2 *; **************************************************; IF TEST=2 THEN DO; N2 = NPLUS-N1; FCRIT = FINV(1-ALPHAC,A,N2); CUMTERM = SDF("F",FCRIT,A,N2,ANC1); AREA = TABLE[NPCASE,3]*CUMTERM; OUT_AREA= OUT_AREA//(NPLUS||AREA||CUMTERM||TABLE[NPCASE,3]); END; **************************************************; * Calculations for Test 4 *; **************************************************; IF TEST=4 THEN DO; N2 = NPLUS-N1; IF Q1=0 THEN F_LOW = 0; ELSE F_LOW = CDF("CHISQUARE",Q1,NU1,0); IF NPCASE=NPTERMS THEN F_HIGH=1; ELSE F_HIGH = CDF("CHISQUARE",Q2,NU1,0); F_DIFF = F_HIGH-F_LOW; W =(NPLUSMIN-N1)/N2; NUPLUS = NPLUS-R; FCRIT = FINV(1-ALPHAC,A,NPLUSMIN-R); *DF SHOULD BE (N_not-R), assuming N+,min=N_not; ALBMIN = (NPLUSMIN-R)/(A*FCRIT); GLOBAL2 = F_DIFF||F_LOW||NU1||W||N2||A||ANC1||ALBMIN||Q1; H ={0.1,0.1}; *MUST BE COLUMN VEC; ERR=ALPHAC#0.01; CALL _QUAD2D01(H,7,ERR,AREA,FLAG); AREA = AREA*F_DIFF; OUT_AREA = OUT_AREA//(NPLUS||AREA); END; END; IF WARN>0 THEN PRINT "WARNING: Bound on possible error due to truncation is:" GAMMA TRUNCERR; IF TEST=2 THEN POWER=OUT_AREA[+,2]; ELSE POWER=1-OUT_AREA[+,2]; RETURN(POWER); FINISH _GETPOW; **********************************************************************; **********************************************************************; * _NDIST.IML *; * Finds the distribution of the final sample size *; * (Reference: Coffey and Muller, Biometrics, 2001) *; * *; **********************************************************************; * REQUIRED INPUTS: *; * *; * ESSENCEX = Essence design matrix *; * WEIGHTS = Weights for each row of ESSENCEX per replication *; * ALPHAT = Target type I error rate *; * POWERT = Target power *; * C = Between subject contrast matrix *; * BETA_PLN = Vector of primary parameters for power calculation *; * SIGMA0 = Variance estimate used for planning (**NOT** std dev) *; * N1 = Internal pilot sample size *; * NPLUSMIN = Minimum size of final sample *; * NPLUSMAX = Maximum size of final sample *; * GAMMA = Ratio of true variance to pre-planned variance *; * RULE = re-estimation rule (0,1,2) *; * ROUND = helps define ROUNDOFF=ALPHA#10##(-ROUND) *; * controls precision in computing Pr{sample size} (1-10).*; * Default = 3. *; **********************************************************************; * OUTPUTS: *; * *; * _TABLE = Table of cutpoints, NPLUS values, pdf, cdf *; **********************************************************************; **********************************************************************; START _NDIST(GAMMA,ALPHAT,POWERT,ESSENCEX,BETA_PLN,C, N1,ROUND,NPLUSMIN,RULE, SIGMA0,NPLUSMAX,WEIGHTS, _TABLE,ISTOP) GLOBAL(_ZERO_); *********************************; *Create initial conditions *; *********************************; R = NCOL(ESSENCEX); *Assume full rank design; A = NROW(C); *Numerator df; NU1 = N1-R; *Error df for internal pilot sample; REPSIZE = SUM(WEIGHTS); *Replication size; IF RULE=1 THEN DO; FCRIT = FINV(1-ALPHAT,A,NU1); NONCENT = FNONCT(FCRIT,A,NU1,1-POWERT); END; ROUNDOFF0=ALPHAT#((0.1)##ROUND); ROUNDOFF1=ALPHAT#((0.1)##(ROUND+2)); IF (ISTOP=0) THEN ROUNDOFF = ROUNDOFF0; ELSE IF ISTOP=1 THEN ROUNDOFF = ROUNDOFF1; THTA_POP = C*BETA_PLN; *True value of THETA in population; XPX_WESS = ESSENCEX`*(DIAG(WEIGHTS)/SUM(WEIGHTS))*ESSENCEX; INVXPXES = SOLVE(XPX_WESS,I(R)); *This will produce error if LTFR; SSH_WESS = (THTA_POP`)*INV(C*(INVXPXES)*C`)*THTA_POP; PROB_HI=0; PDF_N2 =0; _TABLENM ={"Q_NPLUS" "NPLUS" "PDF" "CDF"}; **************************************************************; * Tabulate end point and probability for each REPN value *; **************************************************************; * CAUTION: FNONCT function may blow up at extremes; DO NPLUS=NPLUSMIN TO NPLUSMAX BY REPSIZE; IF RULE=0 THEN DO; NUPLUS = NPLUS - R; FCRIT = FINV(1-ALPHAT,A,NUPLUS); NONCENT = FNONCT(FCRIT,A,NUPLUS,1-POWERT); END; ELSE IF RULE=2 THEN DO; N2 = NPLUS-N1; FCRIT = FINV(1-ALPHAT,A,N2); NONCENT = FNONCT(FCRIT,A,N2,1-POWERT); END; IF (NPLUS=NPLUSMAX) THEN CDFNPLUS = 1; ELSE DO; SIGMANP = (NPLUS#SSH_WESS)/NONCENT; Q_NPLUS = (NU1#SIGMANP)/(SIGMA0#GAMMA); CDFNPLUS = PROBCHI(Q_NPLUS,NU1); *See (2.17) in Coffey and Muller, 1999; END; *Create table containing pdf of final sample size; IF ( ({1}-CDFNPLUS)<=ROUNDOFF ) THEN DO; CDFNPLUS = 1; * Collapse all probs. into maximum; Q_NPLUS = CINV({1}-ROUNDOFF1,NU1); * Let upper bound be chosen so as to leave an 'acceptable' error; PROB_LOW = PROB_HI; PR_NPLUS = 1-PROB_LOW; *See (2.18) in Coffey and Muller, 1999; TEMP = Q_NPLUS||NPLUS||PR_NPLUS||CDFNPLUS; _TABLE = _TABLE//TEMP; GO TO DONEGOOD; END; PROB_LOW = PROB_HI; PROB_HI = CDFNPLUS; PR_NPLUS = PROB_HI-PROB_LOW; *See (2.18) in Coffey and Muller, 1999; IF CDFNPLUS>ROUNDOFF1 THEN DO; TEMP = Q_NPLUS||NPLUS||PR_NPLUS||CDFNPLUS; _TABLE = _TABLE//TEMP; END; END; ********************************************************; * Let the user know if convergence was not obtained *; ********************************************************; NTABLE =NROW(_TABLE); START =MAX(1,NTABLE-20); LAST_FEW=_TABLE[START:NTABLE,*]; PRINT "ERROR: _NDIST Failed to converge with 1000 replications. " _ROUND , LAST_FEW [COLNAME=_TABLENM]; DONEGOOD: ; *Jump here for normal finish; FINISH _NDIST; **********************************************************************; **********************************************************************; * _FINDADJ.IML *; * Returns an adjusted level of alpha to ensure that the maximum test *; * size is no larger than the target type I error rate. This module *; * determines the modification to the critical value required to *; * implement Test 3 (Bounding Method) *; **********************************************************************; * REQUIRED INPUTS: *; * *; * ESSENCEX = Essence design matrix *; * WEIGHTS = Weights for each row of ESSENCEX per replication *; * ALPHAT = Target type I error rate *; * POWERT = Target power *; * C = Between subject contrast matrix *; * BETA_PLN = Vector of primary parameters for power calculations *; * SIGMA0 = Variance estimate used for planning (**NOT** std dev) *; * ROUND = helps define ROUNDOFF=ALPHA#10##(-ROUND) *; * controls precision in computing Pr{sample size} (1-10).*; * Default = 3. *; * N1 = Internal pilot sample size *; * NPLUSMIN = Minimum size of final sample *; * NPLUSMAX = Maximum size of final sample *; * *; **********************************************************************; **********************************************************************; START _FINDADJ(ALPHAT,POWERT,BETA_PLN,ROUND, ESSENCEX,C,N1,NPLUSMIN,NPLUSMAX, SIGMA0,WEIGHTS) GLOBAL(_ZERO_); ************************************************; * Determine maximum type I error rate at ALPHAT*; ************************************************; PRECISION=.001; REL_ALPHA=PRECISION#ALPHAT; ALPH_ADJ = ALPHAT; MAXT=_MAXTS(ALPHAT,ALPH_ADJ,POWERT,BETA_PLN,ROUND, ESSENCEX,C,N1,NPLUSMIN,NPLUSMAX, SIGMA0,WEIGHTS); ***********************************************; * Use a linear interpolation as a first guess *; * for the adjustment (include +/- 10%) *; ***********************************************; ADJFIRST = (ALPHAT)/MAXT[,3]; ADJ_LO = {0.90}#ADJFIRST; ADJ_HI = MIN({1.1}#ADJFIRST,ALPHAT); ***********************************************; * Based on interpolation begin search *; ***********************************************; ADJLIST = ADJ_LO||ADJFIRST||ADJ_HI; CNAMES = {"Adj. Alpha" "Max. Gamma" "Max. Type I Error Rate" "Ratio: Max. to Type I Error Rate"}; INAMES = {"L Adj" "L Gamma" "L Max." "L RATIO" "U Adj" "U Gamma" "U Max." "U RATIO"}; DO I = 1 TO NCOL(ADJLIST); ALPH_ADJ = ADJLIST[1,I]; MAXVAL=_MAXTS(ALPHAT,ALPH_ADJ,POWERT,BETA_PLN,ROUND, ESSENCEX,C,N1,NPLUSMIN,NPLUSMAX, SIGMA0,WEIGHTS); OUT = OUT//(ALPH_ADJ||MAXVAL); END; OUT = OUT//(ALPHAT||MAXT); *print out[format=9.7]; *Ensure have a true lower bound, i. e., < ALPHAT; DO I=1 TO 10 WHILE(OUT[1,4] > {1}); *Note that OUT[1,4]=OUT[1,3]/ALPHAT; ALPH_ADJ = (OUT[1,1]/OUT[1,4]); MAXVAL=_MAXTS(ALPHAT,ALPH_ADJ,POWERT,BETA_PLN,ROUND, ESSENCEX,C,N1,NPLUSMIN,NPLUSMAX, SIGMA0,WEIGHTS); OUT = (ALPH_ADJ||MAXVAL)//OUT; END; IF OUT[1,4]>{1} THEN DO; PRINT "Failed to find lower bound"; PRINT OUT; RETURN(.); END; *PRINT "Initial Values in Search for Adjusted Alpha: " OUT[COLNAME=CNAMES]; *Check to see if finished based on initial values; *ZERO= LOC( (ABS(OUT[*,3]-ALPHAT)0 THEN DO; *MODALPHA = OUT[MIN (ZERO),1]; MODALPHA = OUT[MAX (ZERO),1]; RETURN(MODALPHA); END; **********************; * Determine bounds *; **********************; ILOWER=MAX(LOC(OUT[*,4]<{1})); LOWER=OUT[ILOWER ,*]; UPPER=OUT[ILOWER+{1},*]; * Define Matrix to keep track of Iteration History for future purposes; ITERHIS = LOWER||UPPER; *********************************************; * Use Interval Halving to determine maximum *; * type I error rate *; *********************************************; DO I=1 TO 100; ALPH_ADJ=(LOWER[1,1]+UPPER[1,1])/{2}; MAXVAL=_MAXTS(ALPHAT,ALPH_ADJ,POWERT,BETA_PLN,ROUND, ESSENCEX,C,N1,NPLUSMIN,NPLUSMAX, SIGMA0,WEIGHTS); * Check to see if midpoint equals the desired value; IF ( MAXVAL[,2]<=(ALPHAT+_ZERO_) ) & (ABS(MAXVAL[,2]-ALPHAT) <= REL_ALPHA) THEN DO; *PRINT "Iteration History 1: " out[format=11.9] ITERHIS[COLNAME=INAMES] maxval; RETURN(ALPH_ADJ); END; * Define new boundary points based on maximum at midpoint; ELSE IF (MAXVAL[,3]<1) THEN LOWER=ALPH_ADJ||MAXVAL; ELSE IF (MAXVAL[,3]>1) THEN UPPER=ALPH_ADJ||MAXVAL; ITERHIS=ITERHIS//(LOWER||UPPER); * Check to see if adjusted ALPHA values are close enough together to stop iteration; IF ((ABS(LOWER[1,1] - UPPER[1,1]) ) /UPPER[1,1]) < PRECISION THEN DO; ALPH_ADJ = LOWER[1,1] ; RETURN(ALPH_ADJ); END; END; PRINT "Failed to find adjusted alpha value"; PRINT "Iteration History: " ITERHIS[COLNAME=INAMES]; RETURN(.); FINISH _FINDADJ; **************************************************************************; **************************************************************************; * _MAXTS.IML *; * Returns the maximum type I error rate as a function of GAMMA,the value *; * of GAMMA at which this maximum ocurs for a specified study design *; * using the unadjusted test (Test 0 in Coffey and Muller 2001), the *; * ratio of the maximum alpha to the ALPHAT-the desired type I error rate,*; * and the number of iterations required by the algorithm to converge *; **************************************************************************; * REQUIRED INPUTS: *; * ESSENCEX = Essence design matrix *; * WEIGHTS = Weights for each row of ESSENCEX per replication *; * ALPHA1 = Type I error rate for re-estimating sample size *; * ALPHA2 = Type I error rate for determining critical value *; * POWERT = Target power *; * C = between subject contrast matrix *; * BETA_PLN = Vector of primary parameters for power calculations *; * SIGMA0 = Variance estimate used for planning (**NOT** std dev) *; * ROUND = helps define ROUNDOFF=ALPHA#10##(-ROUND) *; * controls precision in computing Pr{sample size} (1-10). *; * Default = 3. *; * N1 = Internal pilot sample size *; * NPLUSMIN = Minimum size of final sample *; * NPLUSMAX = Maximum size of final sample *; * *; **************************************************************************; **************************************************************************; START _MAXTS(ALPHAT,ALPH_ADJ,POWERT,BETA_PLN,ROUND,ESSENCEX,C,N1,NPLUSMIN,NPLUSMAX,SIGMA0,WEIGHTS) GLOBAL(_ZERO_); NULLBETA = J(NCOL(ESSENCEX),1,0); NPMTEMP=NPLUSMAX; *Fibonacci Search; *Sets FIB_1-FIB_30; FIB={1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025 121393}; DO J=2 TO 26; *Creates needed ratio: FIB_(N-1)/FIB_N; G=G||(FIB[J-1]/FIB[J]); end; _ISTOP=0; EPSILON=1E-7; *N sets precision: btw 2-25, approx. (0.613)^(N-1) precision in betaspace; N=20; *at N=20 this is 1/10946=0.00009 precision in betaspace; LOWER={0}; *initial lower and upper for MAXBETA; UPPER={1}; INT=UPPER-LOWER; *interval in betaspace; X1=1-G[N]; X2=G[N]; BETALIST=X1||X2; *sets initial BETA values (BETA defined as GAMMA/(1+GAMMA)); GAMLIST=BETALIST/(1-BETALIST); *solves for corresponding GAMLIST; FREE ALPHALIST; DO I = 1 TO NCOL(GAMLIST); *Begin loop: solve for initial ALPHAs; GAMMA = GAMLIST[1,I]; FREE _TABLE; RUN _NDIST(GAMMA,ALPHAT,POWERT,ESSENCEX,BETA_PLN,C,N1,ROUND,NPLUSMIN,0, SIGMA0,NPMTEMP,WEIGHTS,_TABLE,_ISTOP); ALPHA = _GETPOW(GAMMA,NULLBETA,ALPH_ADJ,_TABLE,ESSENCEX,C,N1,NPLUSMIN,0,SIGMA0,WEIGHTS); ALPHALIST=ALPHALIST||ALPHA; END; *End loop: solve for initial ALPHAs; ITERATIONS={0}; *ineration counter; DO WHILE (N>1); *Begin loop: find max ALPHA and corresponding GAMMA; ITERATIONS=ITERATIONS+{1}; *determine case; IF (ALPHALIST[1] > ALPHALIST[2]) THEN CASE={1}; *LOWER ALPHALIST[1]) THEN CASE={2}; *B12 THEN DO; BETALIST[2]=BETALIST[1]; BETALIST[1]=LOWER+(1-G[N])*INT; END; IF N>1 THEN DO; GAMMA=BETALIST[1]/({1}-BETALIST[1]); FREE _TABLE; IF N=14 THEN _TABLE=_TABLEMAX14; ELSE IF N=6 THEN _TABLE=_TABLEMAX6; ELSE RUN _NDIST(GAMMA,ALPHAT,POWERT,ESSENCEX,BETA_PLN,C, N1,ROUND, NPLUSMIN,0,SIGMA0,NPMTEMP,WEIGHTS,_TABLE,_ISTOP); ALPHA1=_GETPOW(GAMMA,NULLBETA,ALPH_ADJ,_TABLE,ESSENCEX,C,N1,NPLUSMIN,0,SIGMA0,WEIGHTS); ALPHALIST=ALPHA1||ALPHALIST[1]; END; END; *End CASE=1 calcs; IF CASE={2} THEN DO; *Begin CASE=2 calcs; N=N-1; LOWER=BETALIST[1]; IF N=14 THEN DO; UPPER_GAM=UPPER/({1}-UPPER); FREE _TABLE; RUN _NDIST(UPPER_GAM,ALPHAT,POWERT,ESSENCEX,BETA_PLN,C,N1, ROUND,NPLUSMIN,0,SIGMA0,NPMTEMP,WEIGHTS,_TABLE,_ISTOP); _TABLEMAX14=_TABLE; NTABLE=NROW(_TABLE); NPMTEMP=_TABLE[NTABLE,2]; _ISTOP=1; END; IF N=6 THEN DO; _ISTOP=0; UPPER_GAM=UPPER/({1}-UPPER); FREE _TABLE; RUN _NDIST(UPPER_GAM,ALPHAT,POWERT,ESSENCEX,BETA_PLN,C,N1, ROUND,NPLUSMIN,0,SIGMA0,NPMTEMP,WEIGHTS,_TABLE,_ISTOP); _TABLEMAX6=_TABLE; NTABLE=NROW(_TABLE); NPMTEMP=_TABLE[NTABLE,2]; _ISTOP=1; END; INT=UPPER-LOWER; IF N=2 THEN DO; BETALIST[1]=BETALIST[2]; BETALIST[2]=LOWER+(G[N]*INT)+EPSILON; END; IF N>2 THEN DO; BETALIST[1]=BETALIST[2]; BETALIST[2]=LOWER+(G[N]*INT); END; IF N>1 THEN DO; GAMMA=BETALIST[2]/({1}-BETALIST[2]); FREE _TABLE; IF N=14 THEN _TABLE=_TABLEMAX14; ELSE IF N=6 THEN _TABLE=_TABLEMAX6; ELSE RUN _NDIST(GAMMA,ALPHAT,POWERT,ESSENCEX,BETA_PLN,C, N1,ROUND, NPLUSMIN,0,SIGMA0,NPMTEMP,WEIGHTS,_TABLE,_ISTOP); ALPHA2=_GETPOW(GAMMA,NULLBETA,ALPH_ADJ,_TABLE,ESSENCEX,C,N1,NPLUSMIN,0,SIGMA0,WEIGHTS); ALPHALIST=ALPHALIST[2]||ALPHA2; END; END; *End CASE=2 calcs; *PRINT BETALIST LOWER UPPER ALPHALIST INT N; END; *END LOOP: FIND MAX ALPHA AND CORRESPONDING GAMMA; MAXBETA=MEAN(LOWER,UPPER); *Calcs for 'Converged' case; MAXGAMMA=MAXBETA/({1}-MAXBETA); FREE _TABLE; RUN _NDIST(MAXGAMMA,ALPHAT,POWERT,ESSENCEX,BETA_PLN,C, N1,ROUND, NPLUSMIN,0,SIGMA0,NPMTEMP,WEIGHTS,_TABLE,_ISTOP); MAXALPHA=_GETPOW(MAXGAMMA,NULLBETA,ALPH_ADJ,_TABLE,ESSENCEX,C,N1,NPLUSMIN,0,SIGMA0,WEIGHTS); RATIO=MAXALPHA/ALPHAT; OUT=MAXGAMMA||MAXALPHA||RATIO; RETURN (OUT); FINISH _MAXTS;