Chapter Contents |
Previous |
Next |
Language Reference |
proc iml; /* function to insert in m the geometry column a at level k*/ start table(m,a,k); if ncol(m) = 0 & nrow(m) = 0 then m = j(nrow(a),k,.); if nrow(m) < nrow(a) then m = m// j(nrow(a)-nrow(m),ncol(m),.); if ncol(m) < k then m = m || j(nrow(m),k-ncol(m),.); m[1:nrow(a),k] = a; finish; call table(m,{-6,2},1); call table(m,{-6,3},2); call table(m,{-6,4,5,6},3); call table(m,{-6,4},4); call seq(prob,m) eps = 1.e-8 den="density"; print m; print prob; print density;The following output displays the values returned for m, prob and den, respectively.
proc iml; reset nocenter; /* function to insert in m the geometry column a at level k*/ start table(m,a,k); if ncol(m) = 0 & nrow(m) = 0 then m = j(nrow(a),k,.); if nrow(m) < nrow(a) then m = m// j(nrow(a)-nrow(m),ncol(m),.); if ncol(m) < k then m = m || j(nrow(m),k-ncol(m),.); m[1:nrow(a),k] = a; finish; call table(m,{-20,2},1); call table(m,{-20,20},2); call table(m,{-3,3},3); /**************************************/ /* TSCALE has the default value of 1 */ /**************************************/ call seq(prob1,m) eps = 1.e-8 den="density"; print m[format=f5.] prob1[format=e12.5]; call table(mm,{-20,2},1); call table(mm,{-3,3},2); /* We can show a 2-step separation between the levels */ /* while dropping the intermediate level at 2 */ tscale = { 2 }; call seq(prob2,mm) eps = 1.e-8 den="density" TSCALE=tscale; print mm[format=f5.] prob2[format=e12.5];
proc iml; /*****************************************/ /* first check whether the numbers yield */ /* 0.95 for the alpha level */ /*****************************************/ bm ={-3.663 -2.884 -2.573 -2.375 -2.037, -2.988 -2.537 -2.407 -2.346 -2.156, -2.598 -2.390 -2.390 -2.390 -2.310, -2.446 -2.404 -2.404 -2.404 -2.396}; bplevel = { 0.5 0.25 0.1 0.05}; level = 0.95; /* this the required alpha value */ sigma = diag(sqrt(1:5)); /* global sigma matrix */ do i = 1 to 4; m = bm[i,]; plevel = bplevel[i]; geom = (m//(-m))*sigma; /***************************/ /* Try the null hypothesis */ /***************************/ call seq(prob,geom) eps = 1.e-10; palpha = (prob[2,]-prob[1,])[5]; /**********************************/ /* Try the alternative hypothesis */ /**********************************/ call seqshift(prob,mean,geom,plevel); beta = (prob[2,] -prob[1,])[5]; p = prob[3,]-prob[2,]+prob[1,]; /**********************************/ /* Number of patients per group */ /**********************************/ tn = 4*mean##2; maxn = 5*tn; /*************************************/ /* compute the average sample number */ /*************************************/ asn = tn *( 5 - (4:0) * p`); summary = summary // ( palpha || level || beta || plevel || tn || maxn ||asn); end; print summary[format=10.5];Note that the variables eps and tscale have been internally set to their default values. The following values are returned for the matrix SUMMARY:
proc iml; start func(delta,k) global(level); m = ((1:k))##delta; mm = (-m//m); /*******************************/ /* meet the significance level */ /* by scaling */ /*******************************/ call seqscale(prob,scale,mm,level); return(scale); finish; /*********************************/ /* alpha levels of 0.05 and 0.01 */ /*********************************/ blevel = {0.95 0.99}; do i = 1 to 2; level = blevel[i]; free summary; do delta = 0 to .7 by .1; free row; do k=2 to 5; x = func(delta,k); row = row || x; end; summary = summary //row; end; print summary[format=10.5]; end;The value of SUMMARY for the 0.95 level is as follows.
proc iml; %let nl = 5; start func(plevel) global(level,scale,mean,palpha,beta,tn,asn); m = sqrt((1: &nl)); mm = -m //m; /*******************************/ /* meet the significance level */ /* by scaling */ /*******************************/ call seqscale(prob,scale,mm,level); palpha = (prob[2,]-prob[1,])[&nl]; mm = mm *scale; /*******************************/ /* meet the power condition */ /*******************************/ call seqshift(prob,mean,mm,plevel); return(mean); finish; /****************/ /* alpha = 0.95 */ /****************/ level = 9.50000E-01; bplevel = { 0.5 .25 .1 0.05 0.01}; free summary; do i = 1 to 5; summary = summary || func(bplevel[i]); end; print summary[format=10.5];The value returned for SUMMARY are shown in the following table, and the entries agree with Table 2 of Pocock (1977).
proc iml; %let nl=5; start func(delta) global(level,plevel,mean, scale,alpha,beta,tn,asn); m = ((1: &nl))##delta; mm = (-m//m); /*******************************/ /* meet the significance level */ /*******************************/ call seqscale(prob,scale,mm,level); alpha = (prob[2,]-prob[1,])[&nl]; mm = mm *scale; /*******************************/ /* meet the power condition */ /*******************************/ call seqshift(prob,mean,mm,plevel); beta = (prob[2,]-prob[1,])[&nl]; /*************************************/ /* compute the average sample number */ /*************************************/ p = prob[3,]-prob[2,]+prob[1]; tn = 4*mean##2; /* number per group */ asn = tn *( &nl - p *(%eval(&nl-1):0)`); return(asn); finish;
/**********************************************/ /* set up the global variables needed by func */ /**********************************************/ level = 0.95; plevel = 0.01; /*****************************************/ /* set up the controlling options of the */ /* optimization routine */ /*****************************************/ opt = {0 2 0 1 6}; tc = repeat(.,1,12); tc[1] = 100; tc[7] = 1.e-4; par = { 1.e-13 . 1.e-10 . . .} || . || epsd; /*****************************/ /* provide the initial guess */ /* and let nlpdd do the work */ /*****************************/ delta = 0.5; call nlpdd(rc,rx,"func",delta) opt=opt tc=tc par=par;The following output displays the results.
Optimization Start Parameter Estimates Gradient Objective N Parameter Estimate Function 1 X1 -1.500000 -8.09752 Value of Objective Function = 35.232023082 Double Dogleg Optimization Dual Broyden - Fletcher - Goldfarb - Shanno Update (DBFGS) Without Parameter Scaling Gradient Computed by Finite Differences Number of Parameter Estimates 1 Parameter Estimates 2 Functions (Observations) 2 Optimization Start Active Constraints 0 Criterion = 35.232 Max Abs Gradient Element 8.098 Radius = 1.000 Function Active Objective Iter Restart Calls Constraints Function 1 0 3 0 34.8914 2* 0 4 0 34.8774 3* 0 5 0 34.8774 Iter difcrit maxgrad lambda slope 1 0.3406 1.644 49.273 -0.830 2* 0.0140 0.0440 0 -0.0144 3* 0.00001 0.00013 0 -1E-5 Optimization Results Iterations 3 Function Calls 6 Gradient Calls 5 Active Constraints 0 Criterion 34.877417 Max Grad Element 0.000126832 Slope -0.0000100034 Radius 1 NOTE: FCONV convergence criterion satisfied. Optimization Results Parameter Estimates N Parameter Estimate Gradient 1 X1 0.586554 -0.0001268 Value of Objective Function = 34.877416815
The optimal function value of 34.88 agrees with the entry in Table 2 of Wang and Tsiatis (1987) for five groups, , and .Note that the variables eps and tscale are internally set to their default values. For more information on the NLPDD subroutine, see the section "NLPDD Call". For details on the opt, tc, and par arguments in the NLPDD call, see the "Options Vector" section, the "Termination Criteria" section, and the section "Control Parameters Vector", respectively.
You can replicate other values in Table 2 of Wang and Tsiatis (1987) by changing the values of the variables NL and PLEVEL. You can obtain values from Table 3 by changing the value of the variable LEVEL to 0.99 and specifying NL and PLEVEL accordingly.
This example illustrates how to find the boundaries that minimize ASN given the required significance level and the required power. It replicates some of the results published in Table 3 of Pocock (1982). The IML program computes the domain that
proc iml ; %let nl=5; start func(m) global(level,plevel,sigma,epss, geometry,stgeom,gscale,mean,alpha,beta,tn,asn); m = abs(m); mm = ( -m // m)*sigma; /*******************************/ /* meet the significance level */ /*******************************/ call seqscale(prob,gscale,mm,level) iguess=gscale eps=epss; stgeom = gscale*m; geometry= mm*gscale; alpha = (prob[2,]-prob[1,])[&nl]; /*******************************/ /* meet the power condition */ /*******************************/ call seqshift(prob,mean,geometry,plevel) iguess=mean eps=epss; beta = (prob[2,]-prob[1,])[&nl]; p = prob[3,] - prob[2,]+prob[1,]; /*************************************/ /* compute the average sample number */ /*************************************/ tn = 4*mean##2; /* number per group */ asn = tn *( &nl - p *(%eval(&nl-1):0)`); return(asn); finish; /**********************************************/ /* set up the global variables needed by func */ /**********************************************/ epss = 1.e-8; epso = 1.e-5; level = 9.50000E-01; plevel = 0.05; sigma = diag(sqrt(1:5));
/*****************************************/ /* set up the controlling options of the */ /* optimization routine */ /*****************************************/ opt = {0 2 0 1 6}; tc = repeat(.,1,12); tc[1] = 100; tc[7] = 1.e-4; par = { 1.e-13 . 1.e-10 . . .} || . || epso; /************************************/ /* provide the constraint matrix */ /* we need monotonically increasing */ /* significance levels */ /************************************/ con = { . . . . . . . , . . . . . . . , 1 -1 . . . 1 0 , . 1 -1 . . 1 0 , . . 1 -1 . 1 0 , . . . 1 -1 1 0 }; /*****************************/ /* provide the initial guess */ /* and let nlp do the work */ /*****************************/ m = { 1 1 1 1 1 }; call nlpdd(rc,rx,"func",m) opt=opt blc = con tc=tc par=par; print stgeom;Note that while eps has been set to eps=10^{-8}, tscale has been internally set to its default value. You may choose to run the IML program with and without the specification of the keyword IGUESS to see the effect on the execution time.
Note the following about the optimization process:
Chapter Contents |
Previous |
Next |
Top |
Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.