c This program performs sample size calculations c for balanced ANOVAs c in association with the power Web program. c Currently, it will handle at most 5 factors, at most 5 c levels per factor, and at most 5000 observations. c c It will also handle at most 5 cb1's, 5 cv's, and 5 sigma's. c Steve Verrill, 3/14/94 c Steve Verrill, 11/16/95 c Steve Verrill, 1/11/97 c The parameter and dimension statements in the anocpow, hyptest, c power, and kofn subroutines must be in accord. parameter (maxfac=5,maxlev=25,mxlevp1=maxlev+1, x mxlevm1=maxlev-1,nrdx=5000, x ncdx=mxlevp1,ncdc=mxlevm1) real temp(ncdx) real pred(nrdx),predsort(nrdx) c Note that a copy of xmat is created in g in hyptest. real*8 xmat(nrdx,ncdx) real*8 y(nrdx) real*8 c(ncdx,ncdc) real*8 hypss c??? real*8 resss ??? never assigned a value or used ??? real fracs(maxlev) real*8 cb1(5),sigma(5),cb(mxlevm1),powerl real*8 alpha c c If I were more careful, c I could probably use fewer id vectors. c integer iorigord(nrdx) character*20 id(nrdx),idorig(nrdx) integer idoink(nrdx) integer id2(nrdx) character*20 idbase(nrdx) integer idlevel(nrdx,maxfac) integer numlev(maxfac) character*80 datafile,results character*10 level(maxfac,maxlev) character*80 argv,intwrite common /alphac/ alpha(3) common /nrowc/ nrow common /cvc/ cv(5) call getarg(10,results) open(unit=7,file=results,access="append",status="unknown") call getarg(11,results) write(6,9900) results 9900 format("The results from this run have been appended to", x " the file`",a80,"`

To see how to download", x " this file,", x ' read this description of anonymous ftp', x '.','Alternatively, you should be able to do a', x ' "Save As" or "Print" in your browser.

') write(7,9700) 9700 format(////,1x,"**************************************",//,1x, x "Output from power.2.f",//,1x, x "**************************************",/) c c number of factors c call getarg(1,argv) write(intwrite,7) argv 7 format(a80) read(intwrite,*) numfac c read(3,*) numfac write(6,9501) numfac 9501 format("

The number of factors was ",i4,".

") write(7,9601) numfac 9601 format(/,1x,'The number of factors was ',i2,'.',/) call getarg(2,argv) write(intwrite,7) argv read(intwrite,*) (numlev(i),i = 1,numfac) c read(3,*) (numlev(i),i = 1,numfac) write(6,9502) 9502 format("The numbers of levels were:

") write(7,9602) (numlev(i),i = 1,numfac) 9602 format(/,1x,'The numbers of levels were:',//,1x,5(i2,2x),/) do 9504 i = 1,numfac write(6,9503) numlev(i) 9503 format("

",i4) 9504 continue call getarg(3,argv) write(intwrite,7) argv read(intwrite,*) numrep c read(3,*) numrep write(6,9505) numrep 9505 format("

The number of replicates was ",i4, x".

") write(7,9603) numrep 9603 format(/,1x,'The number of replicates was ',i4,'.',/) c c generate the 0,1 portion of the design matrix c c need dynamic memory allocation??? Yes --- probably don't need c to handle more than 5,000 specimens, but should allow up to 10 c levels per factor (say) and at least five factors and 1000 reps c which would require ... if (numfac .eq. 1) then nrow = 0 nlevtot = numlev(1) do 150 i = 1,numlev(1) do 142 j = 1,numlev(1) temp(j) = 0.0 142 continue temp(i) = 1.0 do 146 ii = 1,numrep nrow = nrow + 1 do 144 j = 1,numlev(1) xmat(nrow,j) = temp(j) 144 continue 146 continue 150 continue elseif (numfac .eq. 2) then nrow = 0 nlevtot = numlev(1) + numlev(2) do 170 i = 1,numlev(1) do 168 j = 1,numlev(2) do 152 k = 1,nlevtot temp(k) = 0.0 152 continue temp(i) = 1.0 temp(numlev(1)+j) = 1.0 do 156 ii = 1,numrep nrow = nrow + 1 do 154 k = 1,nlevtot xmat(nrow,k) = temp(k) 154 continue 156 continue 168 continue 170 continue elseif (numfac .eq. 3) then c write(6,9506) c 9506 format("

Got to 9506.

") nrow = 0 nlevtot = numlev(1) + numlev(2) + numlev(3) do 200 i = 1,numlev(1) do 198 j = 1,numlev(2) do 196 k = 1,numlev(3) do 172 l = 1,nlevtot temp(l) = 0.0 172 continue temp(i) = 1.0 temp(numlev(1)+j) = 1.0 temp(numlev(1)+numlev(2)+k) = 1.0 do 176 ii = 1,numrep nrow = nrow + 1 do 174 l = 1,nlevtot xmat(nrow,l) = temp(l) 174 continue 176 continue 196 continue 198 continue 200 continue elseif (numfac .eq. 4) then nrow = 0 nlevtot = numlev(1) + numlev(2) + numlev(3) + numlev(4) do 240 i = 1,numlev(1) do 239 j = 1,numlev(2) do 238 k = 1,numlev(3) do 237 l = 1,numlev(4) do 202 m = 1,nlevtot temp(m) = 0.0 202 continue temp(i) = 1.0 temp(j+numlev(1)) = 1.0 temp(k+numlev(1)+numlev(2)) = 1.0 temp(l+numlev(1)+numlev(2)+numlev(3)) = 1.0 do 206 ii = 1,numrep nrow = nrow + 1 do 204 m = 1,nlevtot xmat(nrow,m) = temp(m) 204 continue 206 continue 237 continue 238 continue 239 continue 240 continue else c c numfac = 5 c nrow = 0 nlevtot = numlev(1)+numlev(2)+numlev(3)+numlev(4)+numlev(5) do 290 i = 1,numlev(1) do 289 j = 1,numlev(2) do 288 k = 1,numlev(3) do 287 l = 1,numlev(4) do 286 m = 1,numlev(5) do 242 mm = 1,nlevtot temp(mm) = 0.0 242 continue temp(i) = 1.0 temp(j+numlev(1)) = 1.0 temp(k+numlev(1)+numlev(2)) = 1.0 temp(l+numlev(1)+numlev(2)+numlev(3)) = 1.0 temp(m+numlev(1)+numlev(2)+numlev(3)+numlev(4)) x = 1.0 do 246 ii = 1,numrep nrow = nrow + 1 do 244 mm = 1,nlevtot xmat(nrow,mm) = temp(mm) 244 continue 246 continue 286 continue 287 continue 288 continue 289 continue 290 continue endif c c do the power calculations c ncol = nlevtot 3004 alpha(1) = .10 alpha(2) = .05 alpha(3) = .01 call getarg(4,argv) write(intwrite,7) argv read(intwrite,*) idfac write(6,9506) idfac 9506 format("The factor tested was ",i1,".

", x"This was a test for a difference in two", x" pre-specified levels.

") write(7,9604) idfac 9604 format(/,1x,'The factor tested was ',i1,'.',///, x1x,'This was a test for a difference in two', x' pre-specified levels.',/) c read(3,*) idfac c c comparing 2 pre-specified levels --- idtest = 0 c idtest = 0 c 3044 write(6,3045) c write(7,3045) c 3045 format(/,1x,'What are the 2 levels that you wish to compare?', c x/) c read(5,*) idlev1,idlev2 c write(7,*) idlev1,idlev2 c **************************************************** c c Since the design is balanced. The power c should not depend on which 2 levels of the factor c are chosen. c c***************************************************** idlev1 = 1 idlev2 = 2 c if ((idlev1 .lt. 1) .or. (idlev1 .gt. numlev(idfac))) then c write(6,3046) numlev(idfac) c write(7,3046) numlev(idfac) c 3046 format(/,1x,'The levels must lie between 1 and ',i2,/) c go to 3044 c elseif ((idlev2 .lt. 1) .or. (idlev2 .gt. numlev(idfac))) then c write(6,3046) numlev(idfac) c write(7,3046) numlev(idfac) c go to 3044 c elseif (idlev1 .eq. idlev2) then c write(6,3047) c write(7,3047) c 3047 format(/,1x,'You must specify 2 distinct levels.',/) c go to 3044 c endif do 3048 i = 1,ncol c(i,1) = 0.0 3048 continue numpre = 0 do 3049 i = 1,idfac-1 numpre = numpre + numlev(i) 3049 continue i1 = numpre + idlev1 i2 = numpre + idlev2 c(i1,1) = 1 c(i2,1) = -1 c 3052 write(6,3054) c write(7,3054) c 3054 format(/,1x,'How many differences do you want to consider?', c x' (5 or fewer)',/) c read(5,*) numdelta c write(7,*) numdelta call getarg(5,argv) write(intwrite,7) argv read(intwrite,*) numdelta numcb1 = numdelta c if ((numdelta .lt. 1) .or. (numdelta .gt. 5)) then c write(6,3056) c write(7,3056) c 3056 format(/,1x,'The number must lie between 1 and 5', c x ' (inclusive).',/) c go to 3052 c endif call getarg(6,argv) write(intwrite,7) argv read(intwrite,*) (cb1(i),i = 1,numdelta) c write(6,3060) c write(7,3060) c 3060 format(/,1x,'What are they? (e.g., .20 for a 20% difference', c x ' in means)',/) c read(5,*) (cb1(i),i = 1,numdelta) c write(7,*) (cb1(i),i = 1,numdelta) nq = 1 call getarg(7,argv) write(intwrite,7) argv read(intwrite,*) numcv c read(3,*) numcv numsig = numcv call getarg(8,argv) write(intwrite,7) argv read(intwrite,*) (cv(i),i = 1,numcv) c read(3,*) (cv(i),i = 1,numcv) call getarg(9,argv) write(intwrite,7) argv read(intwrite,*) powerl write(7,9605) powerl 9605 format(/,1x,'The power desired was ',f5.3,'.',/) c read(3,*) powerl c c prepare for the power calculations c call hyptest(xmat,y,nrdx,nrow,ncol,c,ncdx,nq,hypss,ihypdof, 1 notestf,2) c the notestf (not estimable flag) should never be anything other than 0 c since we are only feeding it standard designs and considering standard c tests, but could still test for notestf > 0 do 5810 jjj = 1,numcv sigma(jjj) = cv(jjj) 5810 continue call power(nq,cb1,numcb1,sigma,numsig,cb,powerl,idtest) stop end c c matrix addition c subroutine matadd(a,b,c,nrow,ncol,nrdim) implicit double precision (a-h,o-z) dimension a(nrdim,1),b(nrdim,1),c(nrdim,1) do 20 i = 1,nrow do 10 j = 1,ncol c(i,j) = a(i,j) + b(i,j) 10 continue 20 continue return end c c matrix subtraction c subroutine matsub(a,b,c,nrow,ncol,nrdim) implicit double precision (a-h,o-z) dimension a(nrdim,1),b(nrdim,1),c(nrdim,1) do 20 i = 1,nrow do 10 j = 1,ncol c(i,j) = a(i,j) - b(i,j) 10 continue 20 continue return end c c matrix multiplication c subroutine matmult(a,nra,nca,b,nrb,ncb,c,nrda,nrdb) implicit double precision (a-h,o-z) dimension a(nrda,1),b(nrdb,1),c(nrda,1) if (nca .ne. nrb) then write(6,10) c write(7,10) 10 format(1x,'Matrix multiplication error:',/, 1 1x,'The a column and b row dimensions do not match.') stop else do 50 j = 1,ncb do 40 i = 1,nra c(i,j) = 0.0 do 30 k = 1,nca c(i,j) = c(i,j) + a(i,k)*b(k,j) 30 continue 40 continue 50 continue endif return end c c matrix transposition c subroutine mattran(a,nra,nca,at,nrda,ncda) implicit double precision (a-h,o-z) dimension a(nrda,1),at(ncda,1) do 20 i = 1,nra do 10 j = 1,nca at(j,i) = a(i,j) 10 continue 20 continue return end function dot(x,y,n) implicit double precision (a-h,o-z) dimension x(n),y(n) dot = 0.0 do 10 i = 1,n dot = dot + x(i)*y(i) 10 continue return end function vnorm(x,n) implicit double precision (a-h,o-z) dimension x(n) external dot vnorm = sqrt(dot(x,x,n)) return end subroutine vector(c,nrow,j,v,nrdc) c This places the jth column of the matrix, c, c into vector v. implicit double precision (a-h,o-z) dimension c(nrdc,1),v(nrow) do 10 i = 1,nrow v(i) = c(i,j) 10 continue return end subroutine hyptest(x,y,nrdx,nrow,ncol,c,ncdx,nq,hypss,ihypdof, 1 notestf,idesign) c If idesign equals 0, c this subroutine checks whether the c's are estimable and c if they are it tests the hypothesis c1'b = ... = cq'b = 0. c If idesign equals 1 (estimability calculations requested) or c idesign equals 2 (power calculations requested), c this subroutine checks the estimability of the c's and provides for c power calculations. c The parameter and dimension statements in the anocpow, hyptest, c power, and kofn subroutines must be in accord. parameter (maxfac=5,maxlev=25, x nrowd=5000, x nrdc=maxlev+1) implicit double precision (a-h,o-z) dimension x(nrdx,ncol),y(nrdx),c(ncdx,ncol),e(nrdc),work(nrowd) dimension f(nrdc,nrdc),ft(nrdc,nrdc),fft(nrdc,nrdc) dimension vtemp(nrdc),temp1(nrdc),temp2(nrdc) dimension g(nrowd,nrdc),h(nrowd,nrdc) dimension ua(nrowd,nrdc),xlama(nrdc),va(nrdc,nrdc) c What is notest? It stands for "not estimable". It is c a listing of contrasts that do not appear to be estimable. c notestf lets us know whether there any non-estimable contrasts. dimension notest(nrdc) dimension at(nrdc,nrowd),ata(nrdc,nrdc) dimension det(2) c****** this common statement is unneccessary? need to dimension c****** but not in a common statement???? common /hec/ ux(nrowd,nrdc),vx(nrdc,nrdc),xlamx(nrdc), x a(nrowd,nrdc) common /ataic/ atai(nrdc,nrdc) common /ixrankc/ ixrank external vnorm c*** I'm almost certainly be wasteful of memory in this and c*** associated routines, but I don't currently have the c*** time or the incentive to clean it up. It works fine c*** on the machines that I use. 5000x26 --- lots of these c*** so why aren't I running out of space????? Would c*** on a PC? c determining machine double precision xmult = .99d0 delta = 1.0d0 2 delta = xmult*delta if ((1.0d0 + delta) .gt. 1.0d0) go to 2 delta = delta/xmult threshm = sqrt(delta) notestf = 0 do 5 i = 1,nrow do 4 j = 1,ncol g(i,j) = x(i,j) 4 continue c write(7,*) (g(i,j),j = 1,ncol) 5 continue c c ssvdc is the linpack singular value decomposition routine c call dsvdc(g,nrdx,nrow,ncol,xlamx,e,ux,nrdx,vx,ncdx, xwork,21,info) c c If info is not equal to 0 print an error message and stop. c if (info .ne. 0) then write(6,7) info c write(7,7) info 7 format('

ERROR 7: The info value from a singular', x ' value decomposition was ',i3,'.

Call Steve', x ' Verrill at 608-231-9375 for assistance.

') stop endif c c obtain the threshold that is used to determine whether singular values c are "essentially zero" c xnorm = 0.0d0 do 10 i = 1,ncol xnorm = xnorm + xlamx(i)*xlamx(i) 10 continue xnorm = sqrt(xnorm) thresh = xnorm*threshm c c determine the rank of x c do 120 i = 1,ncol if (xlamx(i) .gt. thresh) go to 120 ixrank = i - 1 go to 130 120 continue ixrank = ncol c c determine whether the c's are estimable (there exists an a c such that x'a = c) c 130 iresdof = nrow - ixrank do 150 i = 1,ncol do 140 j = 1,ixrank f(i,j) = vx(i,j) 140 continue 150 continue c c obtain the transpose of f c c******* could just call mattran do 22 i = 1,ncol do 12 j = 1,ixrank ft(j,i) = f(i,j) 12 continue 22 continue c c multiply f by ft to obtain fft c c******* could just call matmult do 52 j = 1,ncol do 42 i = 1,ncol fft(i,j) = 0.0d0 do 32 k = 1,ixrank fft(i,j) = fft(i,j) + f(i,k)*ft(k,j) 32 continue 42 continue 52 continue c c fft is the projection operator onto the space spanned by c the columns of v (or the columns of x') c mm = 0 do 200 i = 1,nq c******* could replace vnorm with BLAS dnrm2 call vector(c,ncol,i,vtemp,ncdx) thresh = vnorm(vtemp,ncol)*threshm call matmult(fft,ncol,ncol,vtemp,ncol,1,temp1,ncdx,ncdx) call matsub(vtemp,temp1,temp2,ncol,1,ncdx) if (vnorm(temp2,ncol) .lt. thresh) go to 200 mm = mm + 1 notest(mm) = i 200 continue if (mm .eq. 0) go to 300 c****** note that the stuff below only allows 99 c's (i2 write format) write(6,210) (notest(i),i = 1,mm) write(7,9210) (notest(i),i = 1,mm) 210 format("

The following c's do not appear to be", x" estimable", 1"

",(1x,10(i2,2x)),"

") 9210 format(///,1x,"The following c's do not appear to be estimable", 1//,(1x,10(i2,2x)),/) notestf = 1 return c c calculate the ai's that satisfy x'ai = ci c 300 if (idesign .ne. 1) go to 305 write(6,302) write(7,9302) 302 format("

The c's all appear to be estimable.

") 9302 format(///,1x,"The c's all appear to be estimable.",/) 305 do 320 i = 1,nrow do 310 j = 1,ixrank g(i,j) = ux(i,j)/xlamx(j) 310 continue 320 continue c c h is the matrix that transforms c's into a's c (h = u(nxr)*diaginv(rxr)*vt(rxp)) c call matmult(g,nrow,ixrank,ft,ixrank,ncol,h,nrdx,ncdx) call matmult(h,nrow,ncol,c,ncol,nq,a,nrdx,ncdx) if (idesign .eq. 0) go to 323 call mattran(a,nrow,nq,at,nrdx,ncdx) call matmult(at,nq,nrow,a,nrow,nq,ata,ncdx,nrdx) c c Taken together, dpofa and dpodi yield an inverse of the c ata matrix. dpofa makes use of the Cholesky algorithm c to factor the positive definite matrix ata. c dpofa and dpodi are linpack routines. Their source is c included in this distribution. See the linpack manual c for details. (Dongarra, Moler, Bunch, and Stewart, c Linpack Users' Guide, SIAM, 1979) c call dpofa(ata,ncdx,nq,info) if (info .ne. 0) then write(6,9999) info c write(7,9999) info 9999 format('

ERROR 9999: The info value from a matrix', x ' inversion was',1x,i3,'

Call Steve Verrill at', x ' 608-231-9375 for assistance.

') stop endif c c calculate the inverse of ata c job = 1 call dpodi(ata,ncdx,nq,det,job) do 322 i = 1,nq do 321 j = i,nq atai(i,j) = ata(i,j) atai(j,i) = ata(i,j) 321 continue 322 continue return c c idesign = 0 c 323 do 325 i = 1,nrow do 324 j = 1,nq g(i,j) = a(i,j) 324 continue 325 continue c c print the a vectors c c write(8,326) c326 format(//,1x,"The a's are (as columns)",/) c do 328 i = 1,nrow c write(8,327) (a(i,j),j = 1,nq) c327 format(10(1x,e11.3)) c328 continue c c perform a singular value decomposition of a c in order to obtain an orthonormal basis of L(a) c (the linear span of the columns of a) c call dsvdc(g,nrdx,nrow,nq,xlama,e,ua,nrdx,va,ncdx, xwork,20,infoa) c It appears to me that can get the non-centrality parameter without c dpofa and dpodi if have ua --- no if only know c'b values c but yes if we know Xb values. if (infoa .ne. 0) then write(6,327) infoa c write(7,327) infoa 327 format('

ERROR 327: The info value from a singular', x ' value decomposition was ',i3,'.

Call Steve', x ' Verrill at 608-231-9375 for assistance.

') stop endif c c obtain the hypothesis sum of squares c hypss = 0.0d0 do 340 j = 1,nq temp0 = 0.0d0 do 330 i = 1,nrow temp0 = temp0 + ua(i,j)*y(i) 330 continue hypss = hypss + temp0*temp0 340 continue ihypdof = nq c c obtain the residual sum of squares c xss = 0.0d0 do 440 j = 1,ixrank temp0 = 0.0d0 do 430 i = 1,nrow temp0 = temp0 + ux(i,j)*y(i) 430 continue xss = xss + temp0*temp0 440 continue totss = 0.0d0 do 450 i = 1,nrow totss = totss + y(i)*y(i) 450 continue c 11/18/95 resss never returned ? resss = totss - xss return end c c c subroutine power(nq,cb1,numcb1,sigma,numsig,cb,powerl,idtest) c c This subroutine performs power calculations. c c The parameter and dimension statements in the anocpow, hyptest, c power, and kofn subroutines must be in accord. parameter (maxlev=25,mxlevm1=maxlev-1,ncdx=maxlev+1) implicit double precision (a-h,o-z) dimension cb1(5),sigma(5),cb(mxlevm1) c What is the 5,5,3? c first five -- sigmas second five -- cb1's three -- alphas dimension powera(5,5,3) dimension numrep(5,5,3) c What is the 3,5? 3 -- alphas 5 -- sigmas dimension halfwidt(3,5) real cv common /alphac/ alpha(3) common /ataic/ atai(ncdx,ncdx) common /nrowc/ nrow common /ixrankc/ ixrank common /cvc/ cv(5) iresdof = nrow - ixrank if (idtest .eq. 1) go to 1000 c c comparison of two particular levels of a factor c d = sqrt(atai(1,1)) do 200 ii = 1,3 f = iresdof c c cdft is a routine from DCDFLIB. It performs calculations c associated with the t distribution. DCDFLIB is a c public domain library of "routines for cumulative distribution c functions, their inverses, and their parameters." It was produced c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and c Kathy Russell of the Department of Biomathematics, c M.D. Anderson Cancer Center, The University of Texas. c One way to obtain it is through statlib. c To get started on statlib, send the e-mail message c c send index c c to statlib@lib.stat.cmu.edu. c The statlib URL is c http://lib.stat.cmu.edu/ c or, if that doesn't work c http://lib.stat.cmu.edu/.index.html c p = 1.0 - alpha(ii)/2.0 call cdft(2,p,1.0d0-p,xbase,f,istatus,bound) if (istatus .ne. 0) then write(6,10) istatus c write(7,10) istatus 10 format('

The istatus value from a T inverse', x ' routine was ',i2,'.

Please call Steve Verrill', x ' at 608-231-9375 for assistance.

') stop endif do 100 i = 1,numsig do 90 j = 1,numcb1 xncp = (cb1(j)*d)/sigma(i) c c cdffnc is a routine from DCDFLIB. It performs calculations c associated with the non-central F distribution. DCDFLIB is a c public domain library of "routines for cumulative distribution c functions, their inverses, and their parameters." It was produced c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and c Kathy Russell of the Department of Biomathematics, c M.D. Anderson Cancer Center, The University of Texas. c One way to obtain it is through statlib. c To get started on statlib, send the e-mail message c c send index c c to statlib@lib.stat.cmu.edu. c The statlib URL is c http://lib.stat.cmu.edu/ c or, if that doesn't work c http://lib.stat.cmu.edu/.index.html c f = iresdof call cdffnc(1,p,q,xbase*xbase,1.0d0,f,xncp*xncp, x istatus,bound) if (istatus .ne. 0) then write(6,1025) istatus c write(7,1025) istatus stop endif powera(i,j,ii) = 1.0d0 - p c c find the number of replications needed to achieve the c desired power ("1 replication" means the original design, c "k replications" means k copies of the original design) c numrep(i,j,ii) = 1 if (powerl .le. powera(i,j,ii)) go to 90 numl = 1 numh = 50 xncpt = xncp*sqrt(50.0d0) idof = 50*nrow - ixrank f = idof c c cdft is a routine from DCDFLIB. It performs calculations c associated with the t distribution. DCDFLIB is a c public domain library of "routines for cumulative distribution c functions, their inverses, and their parameters." It was produced c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and c Kathy Russell of the Department of Biomathematics, c M.D. Anderson Cancer Center, The University of Texas. c One way to obtain it is through statlib. c To get started on statlib, send the e-mail message c c send index c c to statlib@lib.stat.cmu.edu. c The statlib URL is c http://lib.stat.cmu.edu/ c or, if that doesn't work c http://lib.stat.cmu.edu/.index.html c p = 1.0 - alpha(ii)/2.0 call cdft(2,p,1.0d0-p,x,f,istatus,bound) if (istatus .ne. 0) then write(6,10) istatus c write(7,10) istatus stop endif c c cdffnc is a routine from DCDFLIB. It performs calculations c associated with the non-central F distribution. DCDFLIB is a c public domain library of "routines for cumulative distribution c functions, their inverses, and their parameters." It was produced c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and c Kathy Russell of the Department of Biomathematics, c M.D. Anderson Cancer Center, The University of Texas. c One way to obtain it is through statlib. c To get started on statlib, send the e-mail message c c send index c c to statlib@lib.stat.cmu.edu. c The statlib URL is c http://lib.stat.cmu.edu/ c or, if that doesn't work c http://lib.stat.cmu.edu/.index.html c call cdffnc(1,p,q,x*x,1.0d0,f,xncpt*xncpt, x istatus,bound) if (istatus .ne. 0) then write(6,1025) istatus c write(7,1025) istatus stop endif pt = 1.0d0 - p if (pt .ge. powerl) go to 50 c c numrep is set to 100 so that ** is printed in the power table. c ** indicates that the necessary number of replications is c greater than 50. c numrep(i,j,ii) = 100 go to 90 50 idel = numh - numl if (idel .lt. 2) then numrep(i,j,ii) = numh go to 90 endif numt = numl + idel/2 xnumt = numt xncpt = xncp*sqrt(xnumt) idof = numt*nrow - ixrank f = idof c c cdft is a routine from DCDFLIB. It performs calculations c associated with the t distribution. DCDFLIB is a c public domain library of "routines for cumulative distribution c functions, their inverses, and their parameters." It was produced c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and c Kathy Russell of the Department of Biomathematics, c M.D. Anderson Cancer Center, The University of Texas. c One way to obtain it is through statlib. c To get started on statlib, send the e-mail message c c send index c c to statlib@lib.stat.cmu.edu. c The statlib URL is c http://lib.stat.cmu.edu/ c or, if that doesn't work c http://lib.stat.cmu.edu/.index.html c p = 1.0 - alpha(ii)/2.0 call cdft(2,p,1.0d0-p,x,f,istatus,bound) if (istatus .ne. 0) then write(6,10) istatus c write(7,10) istatus stop endif c c cdffnc is a routine from DCDFLIB. It performs calculations c associated with the non-central F distribution. DCDFLIB is a c public domain library of "routines for cumulative distribution c functions, their inverses, and their parameters." It was produced c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and c Kathy Russell of the Department of Biomathematics, c M.D. Anderson Cancer Center, The University of Texas. c One way to obtain it is through statlib. c To get started on statlib, send the e-mail message c c send index c c to statlib@lib.stat.cmu.edu. c The statlib URL is c http://lib.stat.cmu.edu/ c or, if that doesn't work c http://lib.stat.cmu.edu/.index.html c call cdffnc(1,p,q,x*x,1.0d0,f,xncpt*xncpt, x istatus,bound) if (istatus .ne. 0) then write(6,1025) istatus c write(7,1025) istatus stop endif pt = 1.0d0 - p if (pt .ge. powerl) go to 70 numl = numt go to 50 70 numh = numt go to 50 90 continue c c calculate approximate half-widths of confidence intervals c halfwidt(ii,i) = xbase*sigma(i)/d 100 continue 200 continue write(6,205) write(7,9205) 205 format('Power tables

') 9205 format(///,1x,'Power tables',//) do 300 ii = 1,3 write(6,310) alpha(ii) write(7,9310) alpha(ii) 310 format('

The significance level is ',f4.2, x '

') 9310 format(///,1x,'The significance level is ',f4.2,//) write(6,320) write(7,9320) 320 format('

The differences are on top.

', x 'The coefficients of variation are along the left side.', x '

') 9320 format(/,1x,'The differences are on top.',/,1x, x 'The coefficients of variation are along the left side.',//) write(6,8010) 8010 format("") write(6,8012) 8012 format("

") do 330 k = 1,numcb1 write(6,8014) cb1(k) 8014 format(" | ",f5.3) 330 continue write(7,9330) (cb1(i),i = 1,numcb1) 9330 format(/,18x,5(e12.4,3x)) write(7,9332) 9332 format(/) do 340 i = 1,numsig write(6,8020) cv(i) 8020 format(" |
---|---|

",f5.3) do 335 j = 1,numcb1 write(6,8030) powera(i,j,ii) 8030 format(" | ",f5.3) 335 continue 340 continue do 9340 i = 1,numsig write(7,9335) cv(i),(powera(i,j,ii),j = 1,numcb1) 9335 format(1x,e12.4,5x,5(e12.4,3x)) 9340 continue write(6,8040) 8040 format(" |

", x"The values", x" in the tables are the half widths that one would

", x"expect for confidence intervals on the" x" differences.

") 9355 format(///,1x,"Half width tables",//,1x,"The values", x" in the tables are the half widths that one would", x/,1x, x"expect for confidence intervals on the" x" differences.",/) do 390 ii = 1,3 write(6,360) alpha(ii) write(7,9360) alpha(ii) 360 format('

The significance level is ',f4.2, x '

') 9360 format(///,1x,'The significance level is ',f4.2,//) write(6,365) write(7,9365) 365 format("

The coefficients of variation", 1 " are along the left side.

") 9365 format(/,1x,"The coefficients of variation", 1 " are along the left side.",//) write(6,8010) do 375 i = 1,numsig write(6,8020) cv(i) write(6,8050) halfwidt(ii,i) 8050 format("",e12.4) write(7,9370) cv(i),halfwidt(ii,i) 9370 format(1x,e12.4,5x,5(e12.4,3x)) 375 continue write(6,8040) 390 continue write(6,405) write(7,9405) 405 format('

', 1'These tables give us the number of replications of the ', 1'current design

that are needed to achieve the desired', 1' power.

') 9405 format(///,1x,'Replication tables',//, 11x,'These tables give us the number of replications of the ', 1'current design',/,1x,'that are needed to achieve the desired', 1' power.',/) write(6,410) powerl write(7,9410) powerl 410 format('

The desired power is ',f4.2,'

') 9410 format(/,1x,'The desired power is ',f4.2,//) do 500 ii = 1,3 write(6,420) alpha(ii) write(7,9420) alpha(ii) 420 format('

The significance level is ',f4.2, x '

') 9420 format(///,1x,'The significance level is ',f4.2,//) write(6,430) write(7,9430) 430 format('

The differences are on top.

', x 'The coefficients of variation', x ' are along the left side.

') 9430 format(/,1x,'The differences are on top.',/,1x, x 'The coefficients of variation', x ' are along the left side.',//) write(6,8010) write(6,8012) do 440 k = 1,numcb1 write(6,8014) cb1(k) 440 continue write(7,9440) (cb1(i),i = 1,numcb1) 9440 format(/,18x,5(e12.4,3x)) write(7,9442) 9442 format(/) do 450 i = 1,numsig write(6,8020) cv(i) do 445 j = 1,numcb1 write(6,8035) numrep(i,j,ii) 8035 format("",i2) 445 continue write(7,9445) cv(i),(numrep(i,j,ii),j = 1,numcb1) 9445 format(1x,e12.4,5x,5(5x,i2,8x)) 450 continue write(6,8040) 500 continue return c c c 1000 do 1200 ii = 1,3 p = 1.0 - alpha(ii) d1 = nq d2 = iresdof c c cdff is a routine from DCDFLIB. It performs calculations c associated with the F distribution. DCDFLIB is a c public domain library of "routines for cumulative distribution c functions, their inverses, and their parameters." It was produced c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and c Kathy Russell of the Department of Biomathematics, c M.D. Anderson Cancer Center, The University of Texas. c One way to obtain it is through statlib. c To get started on statlib, send the e-mail message c c send index c c to statlib@lib.stat.cmu.edu. c The statlib URL is c http://lib.stat.cmu.edu/ c or, if that doesn't work c http://lib.stat.cmu.edu/.index.html c call cdff(2,p,1.0d0-p,xbase,d1,d2,istatus,bound) if (istatus .ne. 0) then write(6,1005) istatus 1005 format('

Error 1005: The istatus value from an', x ' F inverse', x ' routine was ',i2,'.

Please call Steve Verrill', x ' at 608-231-9375 for assistance.

') stop endif do 1100 i = 1,numsig xncp = 0.0d0 do 1020 j = 1,nq do 1010 k = 1,nq xncp = xncp + atai(j,k)*cb(j)*cb(k) 1010 continue 1020 continue xncp = xncp/(sigma(i)*sigma(i)) c c call the non-central f routine and obtain powera(i,1,ii) c c c cdffnc is a routine from DCDFLIB. It performs calculations c associated with the non-central F distribution. DCDFLIB is a c public domain library of "routines for cumulative distribution c functions, their inverses, and their parameters." It was produced c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and c Kathy Russell of the Department of Biomathematics, c M.D. Anderson Cancer Center, The University of Texas. c One way to obtain it is through statlib. c To get started on statlib, send the e-mail message c c send index c c to statlib@lib.stat.cmu.edu. c The statlib URL is c http://lib.stat.cmu.edu/ c or, if that doesn't work c http://lib.stat.cmu.edu/.index.html c d2 = iresdof call cdffnc(1,p,q,xbase,d1,d2,xncp,istatus,bound) powera(i,1,ii) = 1.0d0 - p if (istatus .ne. 0) then write(6,1025) istatus 1025 format('

Error 1025: The istatus value from', x ' a non-central F cdf', x ' routine was ',i2,'.

Please call Steve Verrill', x ' at 608-231-9375 for assistance.

') stop endif c c find the number of replications needed to achieve the c desired power ("1 replication" means the original design, c "k replications" means k copies of the original design) c numrep(i,1,ii) = 1 if (powerl .le. powera(i,1,ii)) go to 1100 numl = 1 numh = 50 xncpt = xncp*50.0d0 idof = 50*nrow - ixrank d2 = idof c c cdff is a routine from DCDFLIB. It performs calculations c associated with the F distribution. DCDFLIB is a c public domain library of "routines for cumulative distribution c functions, their inverses, and their parameters." It was produced c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and c Kathy Russell of the Department of Biomathematics, c M.D. Anderson Cancer Center, The University of Texas. c One way to obtain it is through statlib. c To get started on statlib, send the e-mail message c c send index c c to statlib@lib.stat.cmu.edu. c The statlib URL is c http://lib.stat.cmu.edu/ c or, if that doesn't work c http://lib.stat.cmu.edu/.index.html c p = 1.0 - alpha(ii) call cdff(2,p,1.0d0-p,x,d1,d2,istatus,bound) if (istatus .ne. 0) then write(6,1035) istatus 1035 format('

Error 1035: The istatus value from an', x ' F inverse', x ' routine was ',i2,'.

Please call Steve Verrill', x ' at 608-231-9375 for assistance.

') stop endif c c cdffnc is a routine from DCDFLIB. It performs calculations c associated with the non-central F distribution. DCDFLIB is a c public domain library of "routines for cumulative distribution c functions, their inverses, and their parameters." It was produced c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and c Kathy Russell of the Department of Biomathematics, c M.D. Anderson Cancer Center, The University of Texas. c One way to obtain it is through statlib. c To get started on statlib, send the e-mail message c c send index c c to statlib@lib.stat.cmu.edu. c The statlib URL is c http://lib.stat.cmu.edu/ c or, if that doesn't work c http://lib.stat.cmu.edu/.index.html c call cdffnc(1,p,q,x,d1,d2,xncpt,istatus,bound) pt = 1.0d0 - p if (istatus .ne. 0) then write(6,1040) istatus 1040 format('

Error 1040: The istatus value from', x ' a non-central F cdf', x ' routine was ',i2,'.

Please call Steve Verrill', x ' at 608-231-9375 for assistance.

') stop endif c if (pt .ge. powerl) go to 1050 c c numrep is set to 100 so that ** is printed in the power table. c ** indicates that the necessary number of replications is c greater than 50. c numrep(i,1,ii) = 100 go to 1100 1050 idel = numh - numl if (idel .lt. 2) then numrep(i,1,ii) = numh go to 1100 endif 1060 numt = numl + idel/2 xnumt = numt xncpt = xncp*xnumt idof = numt*nrow - ixrank d2 = idof c c cdff is a routine from DCDFLIB. It performs calculations c associated with the F distribution. DCDFLIB is a c public domain library of "routines for cumulative distribution c functions, their inverses, and their parameters." It was produced c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and c Kathy Russell of the Department of Biomathematics, c M.D. Anderson Cancer Center, The University of Texas. c One way to obtain it is through statlib. c To get started on statlib, send the e-mail message c c send index c c to statlib@lib.stat.cmu.edu. c The statlib URL is c http://lib.stat.cmu.edu/ c or, if that doesn't work c http://lib.stat.cmu.edu/.index.html c p = 1.0 - alpha(ii) call cdff(2,p,1.0d0-p,x,d1,d2,istatus,bound) if (istatus .ne. 0) then write(6,1062) istatus 1062 format('

Error 1062: The istatus value from an', x ' F inverse', x ' routine was ',i2,'.

Please call Steve Verrill', x ' at 608-231-9375 for assistance.

') stop endif c c cdffnc is a routine from DCDFLIB. It performs calculations c associated with the non-central F distribution. DCDFLIB is a c public domain library of "routines for cumulative distribution c functions, their inverses, and their parameters." It was produced c by Barry Brown (bwb@odin.mda.uth.tmc.edu), James Lovato, and c Kathy Russell of the Department of Biomathematics, c M.D. Anderson Cancer Center, The University of Texas. c One way to obtain it is through statlib. c To get started on statlib, send the e-mail message c c send index c c to statlib@lib.stat.cmu.edu. c The statlib URL is c http://lib.stat.cmu.edu/ c or, if that doesn't work c http://lib.stat.cmu.edu/.index.html c call cdffnc(1,p,q,x,d1,d2,xncpt,istatus,bound) pt = 1.0d0 - p if (istatus .ne. 0) then write(6,1064) istatus 1064 format('

Error 1064: The istatus value from', x ' a non-central F cdf', x ' routine was ',i2,'.

Please call Steve Verrill', x ' at 608-231-9375 for assistance.

') stop endif if (pt .ge. powerl) go to 1070 numl = numt go to 1050 1070 numh = numt go to 1050 1100 continue 1200 continue write(6,205) write(7,9205) do 1300 ii = 1,3 write(6,310) alpha(ii) write(7,9310) alpha(ii) write(6,1320) write(7,9321) 1320 format('The coefficients of variation', x ' and associated power values are

') 9321 format(/,1x,'The coefficients of variation', x ' and associated power values are',//) write(6,8010) do 1340 i = 1,numsig write(6,8020) cv(i) write(6,8030) powera(i,1,ii) write(7,9336) cv(i),(powera(i,1,ii)) 9336 format(1x,e12.4,5x,e12.4) 1340 continue write(6,8040) 1300 continue write(6,405) write(7,9405) write(6,410) powerl write(7,9410) powerl do 1500 ii = 1,3 write(6,420) alpha(ii) write(7,9420) alpha(ii) write(6,1430) write(7,9431) 1430 format('The coefficients of variation', x ' and associated replication values are

') 9431 format(/,1x,'The coefficients of variation', x ' and associated replication values are',//) write(6,8010) do 1450 i = 1,numsig write(6,8020) cv(i) write(6,8035) numrep(i,1,ii) write(7,9446) cv(i),numrep(i,1,ii) 9446 format(1x,e12.4,10x,i2) 1450 continue write(6,8040) 1500 continue return end c c c subroutine kofn(k,n,id) c c This program draws (without replacement) a random sample of size c k from a set of size n. c c IMPORTANT: It assumes that uni has already been initialized c via a dum = uni(istart) call. c c c Steve Verrill, 5/11/89 c c c c The parameter and dimension statements in the anocpow, hyptest, c power, and kofn subroutines must be in accord. parameter (idim=5000) integer id(idim),iorigord(idim) c do 30 i = 1,n iorigord(i) = i 30 continue c do 100 i = 0,k - 1 numleft = n - i u = uni(0) l = numleft*u + 1 id(i + 1) = iorigord(l) do 90 j = l,numleft - 1 iorigord(j) = iorigord(j + 1) 90 continue 100 continue end