c
c glm.f
c
c This program can be the basis for a linear models
c package. Currently it creates no plots, checks no
c residuals, and deals only with fixed effects models.
c On the other hand, it can handle unbalanced data, nested
c effects, analysis of covariance, and arbitrary interactions.
c (The X matrix must be provided in a file by the user).
c It also performs power and sample size calculations.
c
c
c Steve Verrill
c
c
c 10/9/84
c 2/8/85
c 2/16/85
c 12/5/90 (and earlier) - modified to run on a Sun
c 5/91 - more modifications
c 7/99 - some cleanup
c
c The parameter and dimension statements in the main, ccheck,
c hyptest, estim, and power subroutines must be in accord.
parameter (nrdx=400,numpar=75,nrdc=numpar,numest=numpar)
implicit double precision (a-h,o-z)
dimension x(nrdx,numpar),y(nrdx)
dimension c(nrdc,numest)
dimension est(numest),cli(numest),cui(numest),
x clj(numest),cuj(numest)
dimension cb1(5),sigma(5),cb(10)
character*80 dataf,results
common /hec2/ resss,iresdof
common /ataic/ atai(numest,numest)
common /alphac/ alpha(3)
common /nrowc/ nrow
common /threshmc/ threshm
c
c Determine machine precision
c
xmult = .99d0
delta = 1.0d0
1 delta = xmult*delta
if ((1.0d0 + delta) .gt. 1.0d0) go to 1
delta = delta/xmult
c threshm = sqrt(delta)
threshm = delta*10.0d0
write(6,2)
2 format(/,1x,'What name do you want for your results file?',/)
read(5,3) results
3 format(a80)
open(unit=9,file=results,status='unknown')
write(6,4)
write(9,4)
4 format(/,1x,'What is the name of the data file?',/)
read(5,3) dataf
write(9,3) dataf
open(unit=3,file=dataf,status='old')
open(unit=4,file='runf',status='unknown')
open(unit=7,file='runfc',status='unknown')
open(unit=8,file='csandas',status='unknown')
iend = 0
c The input that is used to produce the current run
c is written to runfc.
c If, later, you want to reproduce the run with
c a different data file,
c copy runfc to runf and answer yes to the
c following question:
write(6,10)
write(9,10)
10 format(/,1x,'Do you want to read the run control cards',
x' from the file runf?',/,1x,'0 - no 1 - yes',/)
read(5,*) irunf
write(9,*) irunf
if (irunf .eq. 1) go to 5000
write(6,1030)
write(9,1030)
1030 format(/,1x,'How many rows are there in the X matrix',
x' (400 or fewer)?',/)
read(5,*) nrow
write(9,*) nrow
write(7,*) nrow
write(6,1040)
write(9,1040)
1040 format(/,1x,'How many columns are there in the X matrix',
x' (75 or fewer)?',/)
read(5,*) ncol
write(9,*) ncol
write(7,*) ncol
c
c Read the X matrix
c
do 1100 i = 1,nrow
read(3,*) (x(i,j),j = 1,ncol)
write(9,*) (x(i,j),j = 1,ncol)
1100 continue
write(6,1152)
write(9,1152)
1152 format(/,1x,'Do you want to evaluate the design of an experiment',
x', or',/,1x,'do you want to analyze the results from an ',
x'experiment?',//,1x,'Design - 0 Analysis - 1',/)
read(5,*) idfl
write(9,*) idfl
write(7,*) idfl
if (idfl .eq. 0) go to 3004
c
c Read the y vector
c
do 1160 i = 1,nrow
read(3,*) y(i)
write(9,*) y(i)
1160 continue
c
c Hypothesis testing and Estimation
c
1170 write(6,1200)
write(9,1200)
1200 format(/,1x,"Do you wish to test the hypothesis",
x" c1'b = c2'b = ... = cq'b = 0 ",/,1x,"for some q and some c's?",
x" 0 - no 1 - yes",/)
read(5,*) ihypf
write(9,*) ihypf
write(7,*) ihypf
if (ihypf .eq. 0) then
if (iend .eq. 1) stop
iend = 1
go to 2000
endif
write(6,1210)
write(9,1210)
1210 format(/,1x,'What is q?',/)
read(5,*) nq
write(9,*) nq
write(7,*) nq
write(6,1220)
write(9,1220)
1220 format(/,1x,'What are the c primes?',/)
write(8,1225)
1225 format(//,1x,"The c's are (as rows)",/)
do 1300 iii = 1,nq
read(5,*) (c(k,iii), k = 1,ncol)
write(9,*) (c(k,iii), k = 1,ncol)
write(7,*) (c(k,iii), k = 1,ncol)
write(8,*) (c(k,iii), k = 1,ncol)
1300 continue
call ccheck(c,ncol,nq,notind,nrdc,numest)
iend = 0
if (notind .eq. 1) go to 1170
c
c Perform the hypothesis tests
c
call hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
x notestf,0)
iend = 0
if (notestf .eq. 1) go to 1170
c
c Write the results of the hypothesis tests
c
write(6,1402) hypss
write(9,1402) hypss
1402 format(///,1x,'The hypothesis sum of squares is ',2x,e12.4,/)
write(6,1410) ihypdof
write(9,1410) ihypdof
1410 format(/,1x,'The hypothesis dof is',2x,i3,/)
write(6,1420) resss
write(9,1420) resss
1420 format(//,1x,'The residual sum of squares is ',2x,e12.4,/)
write(6,1430) iresdof
write(9,1430) iresdof
1430 format(/,1x,'The residual dof is',2x,i3,/)
xnum = ihypdof
xden = iresdof
f = (hypss/xnum)/(resss/xden)
write(6,1440) f
write(9,1440) f
1440 format(//,1x,'The f statistic is',2x,e12.4,/)
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.html
c
call cdff(1,p,q,f,xnum,xden,istatus,bound)
if (istatus .ne. 0) then
write(6,1442) istatus
1442 format(/,1x,'Error 1442: The istatus value from an',
x ' F routine was ',i2,'.',/)
stop
endif
pvalue = q
write(6,1445) pvalue
write(9,1445) pvalue
1445 format(/,1x,'The p value is',2x,e12.4,/)
c
c Estimation
c
2000 write(6,2200)
write(9,2200)
2200 format(/,1x,"Do you wish to estimate linear combinations",
x" of the parameters c1'b, c2'b, ... , cq'b",/,1x,
x"for some q and for some c's? 0 - no 1 - yes",/)
read(5,*) iestf
write(9,*) iestf
write(7,*) iestf
if (iestf .eq. 0) then
if (iend .eq. 1) stop
iend = 1
go to 1170
endif
isame = 0
if (ihypf .eq. 0) go to 2209
write(6,2205)
write(9,2205)
2205 format(/,1x,"Are the c's the same as those used in the",
x" hypothesis testing? 0 - no 1 - yes",/)
read(5,*) isame
write(9,*) isame
write(7,*) isame
if (isame .eq. 1) go to 2500
2209 write(6,2210)
write(9,2210)
2210 format(/,1x,'What is q?',/)
read(5,*) nq
write(9,*) nq
write(7,*) nq
write(6,2220)
write(9,2220)
2220 format(/,1x,'What are the c primes?',/)
write(8,2225)
2225 format(//,1x,"The c's are (as rows)",/)
do 2300 iii = 1,nq
read(5,*) (c(k,iii), k = 1,ncol)
write(9,*) (c(k,iii), k = 1,ncol)
write(7,*) (c(k,iii), k = 1,ncol)
write(8,*) (c(k,iii), k = 1,ncol)
2300 continue
2500 write(6,2510)
write(9,2510)
2510 format(/,1x,'What are the individual and the joint alphas?',/,
x1x,'(confidence level = (1 - alpha)*100 %)',/)
read(5,*) alphai,alphaj
write(9,*) alphai,alphaj
write(7,*) alphai,alphaj
c
c Do the estimations
c
call estim(x,y,nrdx,nrow,ncol,c,nrdc,nq,alphai,alphaj,
x est,cli,cui,clj,cuj,notestf)
iend = 0
if (notestf .eq. 1) go to 2000
c
c Write the results of the estimations
c
conli = (1.0 - alphai)*100.0
conlj = (1.0 - alphaj)*100.0
write(6,2590) conli
write(9,2590) conli
2590 format(///,1x,'The confidence level for the individual',
x' confidence intervals is ',f5.1,' percent',/)
write(6,2595) conlj
write(9,2595) conlj
2595 format(/,1x,'The confidence level for the joint',
x' confidence intervals is ',f5.1,' percent',/)
write(6,2600)
write(9,2600)
2600 format(/,1x,'The estimates, individual confidence intervals,',
x' and joint confidence intervals are',/)
write(6,2602)
write(9,2602)
2602 format(/,12x,'estimate',14x,'individual',21x,'joint',/)
do 2620 i = 1,nq
write(6,2610) i,est(i),cli(i),cui(i),clj(i),cuj(i)
write(9,2610) i,est(i),cli(i),cui(i),clj(i),cuj(i)
2610 format(/,1x,"c",i2,"'b",4x,e12.4,4x,e12.4,1x,e12.4,
x 4x,e12.4,1x,e12.4)
2620 continue
go to 1170
c
c Experimental design material --- estimability and power
c
3004 alpha(1) = .10
alpha(2) = .05
alpha(3) = .01
3005 write(6,3010)
write(9,3010)
3010 format(/,1x,"Do you wish to check whether some c'bs are estimable"
x,"? no - 0 yes - 1",/)
read(5,*) iestable
write(9,*) iestable
write(7,*) iestable
if (iestable .eq. 0) then
if (iend .eq. 1) stop
iend = 1
go to 4000
endif
write(6,3020)
write(9,3020)
3020 format(/,1x,'How many c primes are there?',/)
read(5,*) nq
write(9,*) nq
write(7,*) nq
write(6,3030)
write(9,3030)
3030 format(/,1x,'What are the c primes?',/)
do 3050 iii = 1,nq
read(5,*) (c(k,iii), k = 1,ncol)
write(9,*) (c(k,iii), k = 1,ncol)
write(7,*) (c(k,iii), k = 1,ncol)
3050 continue
call ccheck(c,ncol,nq,notind,nrdc,numest)
iend = 0
if (notind .eq. 1) go to 3005
c
c Check the estimability of the c's
c
call hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
x notestf,1)
iend = 0
if (notestf .eq. 1) go to 3005
c
c Perform a power calculation
c
4000 write(6,4010)
write(9,4010)
4010 format(/,1x,'Do you want to perform a power calculation?',
x' no - 0 yes - 1',/)
read(5,*) ipowfl
write(9,*) ipowfl
write(7,*) ipowfl
if (ipowfl .eq. 0) then
if (iend .eq. 1) stop
iend = 1
go to 3005
endif
ialready = 0
if (iestable .eq. 0) go to 4040
write(6,4030)
write(9,4030)
4030 format(/,1x,'Are the c primes the same as those for which',
x' you just performed estimability calculations?',/,1x,
x'no - 0 yes - 1',/)
read(5,*) ialready
write(9,*) ialready
write(7,*) ialready
if (ialready .eq. 1) go to 4100
4040 write(6,4050)
write(9,4050)
4050 format(/,1x,"What is the q in the hypothesis c1'b=...=cq'b",
x"=0 ?",/)
read(5,*) nq
write(9,*) nq
write(7,*) nq
write(6,4060)
write(9,4060)
4060 format(/,1x,'What are the c primes?',/)
do 4070 iii = 1,nq
read(5,*) (c(k,iii), k = 1,ncol)
write(9,*) (c(k,iii), k = 1,ncol)
write(7,*) (c(k,iii), k = 1,ncol)
4070 continue
call ccheck(c,ncol,nq,notind,nrdc,numest)
iend = 0
if (notind .eq. 1) go to 4000
4100 if (nq .gt. 1) go to 4109
write(6,4101)
write(9,4101)
4101 format(/,1x,"How many c'b values (5 or fewer) do you want ",
x"in the power table?",/)
read(5,*) numcb
write(9,*) numcb
write(7,*) numcb
write(6,4102)
write(9,4102)
4102 format(/,1x,"What are the c'b values?",/)
read(5,*) (cb1(i),i = 1,numcb)
write(9,*) (cb1(i),i = 1,numcb)
write(7,*) (cb1(i),i = 1,numcb)
go to 4125
4109 write(6,4110)
write(9,4110)
4110 format(/,1x,"What are the values of the c'bs under the ",
x"alternative hypothesis?",/)
read(5,*) (cb(i),i = 1,nq)
write(9,*) (cb(i),i = 1,nq)
write(7,*) (cb(i),i = 1,nq)
4125 write(6,4130)
write(9,4130)
4130 format(/,1x,'How many values (5 or fewer) of sigma ',
x'do you want in the power table?',/)
read(5,*) numsig
write(9,*) numsig
write(7,*) numsig
write(6,4140)
write(9,4140)
4140 format(/,1x,'What are the sigma values?',/)
read(5,*) (sigma(i),i = 1,numsig)
write(9,*) (sigma(i),i = 1,numsig)
write(7,*) (sigma(i),i = 1,numsig)
write(6,4160)
write(9,4160)
4160 format(/,1x,'What power do you wish to achieve?',/)
read(5,*) powerl
write(9,*) powerl
write(7,*) powerl
if (ialready .eq. 1) go to 4170
call hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
x notestf,2)
iend = 0
if (notestf .eq. 1) go to 4000
4170 call power(nq,cb1,numcb,sigma,numsig,cb,powerl,iresdof)
iend = 0
go to 3005
c
c End of experimental design material
c
c
c Reading the run control cards from a file
c
5000 write(6,5030)
write(9,5030)
5030 format(/,1x,'How many rows are there in the X matrix',
x' (400 or fewer)?',/)
read(4,*) nrow
write(6,*) nrow
write(9,*) nrow
write(6,5040)
write(9,5040)
5040 format(/,1x,'How many columns are there in the X matrix',
x' (75 or fewer)?',/)
read(4,*) ncol
write(6,*) ncol
write(9,*) ncol
c
c Read the X matrix
c
do 5100 i = 1,nrow
read(3,*) (x(i,j),j = 1,ncol)
write(9,*) (x(i,j),j = 1,ncol)
5100 continue
iend = 0
write(6,5152)
write(9,5152)
5152 format(/,1x,'Do you want to evaluate the design of an experiment',
x', or',/,1x,'do you want to analyze the results from an ',
x'experiment?',//,1x,'Design - 0 Analysis - 1',/)
read(4,*) idfl
write(6,*) idfl
write(9,*) idfl
if (idfl .eq. 0) go to 7004
c
c Read the y vector
c
do 5160 i = 1,nrow
read(3,*) y(i)
write(9,*) y(i)
5160 continue
c
c Hypothesis testing and Estimation
c
5170 write(6,5200)
write(9,5200)
5200 format(/,1x,"Do you wish to test the hypothesis",
x" c1'b = c2'b = ... = cq'b = 0 ",/,1x,"for some q and some c's?",
x" 0 - no 1 - yes",/)
read(4,*) ihypf
write(6,*) ihypf
write(9,*) ihypf
if (ihypf .eq. 0) then
if (iend .eq. 1) stop
iend = 1
go to 6000
endif
write(6,5210)
write(9,5210)
5210 format(/,1x,'What is q?',/)
read(4,*) nq
write(6,*) nq
write(9,*) nq
write(6,5220)
write(9,5220)
5220 format(/,1x,'What are the c primes?',/)
write(8,5225)
5225 format(//,1x,"The c's are (as rows)",/)
do 5300 iii = 1,nq
read(4,*) (c(k,iii), k = 1,ncol)
write(6,*) (c(k,iii), k = 1,ncol)
write(9,*) (c(k,iii), k = 1,ncol)
write(8,*) (c(k,iii), k = 1,ncol)
5300 continue
call ccheck(c,ncol,nq,notind,nrdc,numest)
iend = 0
if (notind .eq. 1) go to 5170
c
c Perform the hypothesis tests
c
call hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
x notestf,0)
iend = 0
if (notestf .eq. 1) go to 5170
c
c Write the results of the hypothesis tests
c
write(6,1402) hypss
write(9,1402) hypss
write(6,1410) ihypdof
write(9,1410) ihypdof
write(6,1420) resss
write(9,1420) resss
write(6,1430) iresdof
write(9,1430) iresdof
xnum = ihypdof
xden = iresdof
f = (hypss/xnum)/(resss/xden)
write(6,1440) f
write(9,1440) f
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.html
c
call cdff(1,p,q,f,xnum,xden,istatus,bound)
if (istatus .ne. 0) then
write(6,5442) istatus
5442 format(/,1x,'Error 5442: The istatus value from an',
x ' F routine was ',i2,'.',/)
stop
endif
pvalue = q
write(6,1445) pvalue
write(9,1445) pvalue
c
c Estimation
c
6000 write(6,6200)
write(9,6200)
6200 format(/,1x,"Do you wish to estimate linear combinations",
x" of the parameters c1'b, c2'b, ... , cq'b",/,1x,
x"for some q and for some c's? 0 - no 1 - yes",/)
read(4,*) iestf
write(6,*) iestf
write(9,*) iestf
if (iestf .eq. 0) then
if (iend .eq. 1) stop
iend = 1
go to 5170
endif
isame = 0
if (ihypf .eq. 0) go to 6209
write(6,6205)
write(9,6205)
6205 format(/,1x,"Are the c's the same as those used in the",
x" hypothesis testing? 0 - no 1 - yes",/)
read(4,*) isame
write(6,*) isame
write(9,*) isame
if (isame .eq. 1) go to 6500
6209 write(6,6210)
write(9,6210)
6210 format(/,1x,'What is q?',/)
read(4,*) nq
write(6,*) nq
write(9,*) nq
write(6,6220)
write(9,6220)
6220 format(/,1x,'What are the c primes?',/)
write(8,6225)
6225 format(//,1x,"The c's are (as rows)",/)
do 6300 iii = 1,nq
read(4,*) (c(k,iii), k = 1,ncol)
write(6,*) (c(k,iii), k = 1,ncol)
write(9,*) (c(k,iii), k = 1,ncol)
write(8,*) (c(k,iii), k = 1,ncol)
6300 continue
6500 write(6,6510)
write(9,6510)
6510 format(/,1x,'What are the individual and the joint alphas?',/,
x1x,'(confidence level = (1 - alpha)*100 %)',/)
read(4,*) alphai,alphaj
write(6,*) alphai,alphaj
write(9,*) alphai,alphaj
c
c Do the estimations
c
call estim(x,y,nrdx,nrow,ncol,c,nrdc,nq,alphai,alphaj,
x est,cli,cui,clj,cuj,notestf)
iend = 0
if (notestf .eq. 1) go to 6000
c
c Write the results of the estimations
c
conli = (1.0 - alphai)*100.0
conlj = (1.0 - alphaj)*100.0
write(6,6590) conli
write(9,6590) conli
6590 format(///,1x,'The confidence level for the individual',
x' confidence intervals is ',f5.1,' percent',/)
write(6,6595) conlj
write(9,6595) conlj
6595 format(/,1x,'The confidence level for the joint',
x' confidence intervals is ',f5.1,' percent',/)
write(6,6600)
write(9,6600)
6600 format(/,1x,'The estimates, individual confidence intervals,',
x' and joint confidence intervals are',/)
write(6,6602)
write(9,6602)
6602 format(/,12x,'estimate',14x,'individual',21x,'joint',/)
do 6620 i = 1,nq
write(6,6610) i,est(i),cli(i),cui(i),clj(i),cuj(i)
write(9,6610) i,est(i),cli(i),cui(i),clj(i),cuj(i)
6610 format(/,1x,"c",i2,"'b",4x,e12.4,4x,e12.4,1x,e12.4,
x 4x,e12.4,1x,e12.4)
6620 continue
go to 5170
c
c Experimental design material --- estimability and power
c
7004 alpha(1) = .10
alpha(2) = .05
alpha(3) = .01
7005 write(6,7010)
write(9,7010)
7010 format(/,1x,"Do you wish to check whether some c'bs are estimable"
x,"? no - 0 yes - 1",/)
read(4,*) iestable
write(6,*) iestable
write(9,*) iestable
if (iestable .eq. 0) then
if (iend .eq. 1) stop
iend = 1
go to 8000
endif
write(6,7020)
write(9,7020)
7020 format(/,1x,'How many c primes are there?',/)
read(4,*) nq
write(6,*) nq
write(9,*) nq
write(6,7030)
write(9,7030)
7030 format(/,1x,'What are the c primes?',/)
do 7050 iii = 1,nq
read(4,*) (c(k,iii), k = 1,ncol)
write(6,*) (c(k,iii), k = 1,ncol)
write(9,*) (c(k,iii), k = 1,ncol)
7050 continue
call ccheck(c,ncol,nq,notind,nrdc,numest)
iend = 0
if (notind .eq. 1) go to 7005
c
c Check the estimability of the c's
c
call hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
x notestf,1)
iend = 0
if (notestf .eq. 1) go to 7005
c
c Perform a power calculation
c
8000 write(6,8010)
write(9,8010)
8010 format(/,1x,'Do you want to perform a power calculation?',
x' no - 0 yes - 1',/)
read(4,*) ipowfl
write(6,*) ipowfl
write(9,*) ipowfl
if (ipowfl .eq. 0) then
if (iend .eq. 1) stop
iend = 1
go to 7005
endif
ialready = 0
if (iestable .eq. 0) go to 8040
write(6,8030)
write(9,8030)
8030 format(/,1x,'Are the c primes the same as those for which',
x' you just performed estimability calculations?',/,1x,
x'no - 0 yes - 1',/)
read(4,*) ialready
write(6,*) ialready
write(9,*) ialready
if (ialready .eq. 1) go to 8100
8040 write(6,8050)
write(9,8050)
8050 format(/,1x,"What is the q in the hypothesis c1'b=...=cq'b",
x"=0 ?",/)
read(4,*) nq
write(6,*) nq
write(9,*) nq
write(6,8060)
write(9,8060)
8060 format(/,1x,'What are the c primes?',/)
do 8070 iii = 1,nq
read(4,*) (c(k,iii), k = 1,ncol)
write(6,*) (c(k,iii), k = 1,ncol)
write(9,*) (c(k,iii), k = 1,ncol)
8070 continue
call ccheck(c,ncol,nq,notind,nrdc,numest)
iend = 0
if (notind .eq. 1) go to 8000
8100 if (nq .gt. 1) go to 8109
write(6,8101)
write(9,8101)
8101 format(/,1x,"How many c'b values (5 or fewer) do you want ",
x"in the power table?",/)
read(4,*) numcb
write(6,*) numcb
write(9,*) numcb
write(6,8102)
write(9,8102)
8102 format(/,1x,"What are the c'b values?",/)
read(4,*) (cb1(i),i = 1,numcb)
write(6,*) (cb1(i),i = 1,numcb)
write(9,*) (cb1(i),i = 1,numcb)
go to 8125
8109 write(6,8110)
write(9,8110)
8110 format(/,1x,"What are the values of the c'bs under the ",
x"alternative hypothesis?",/)
read(4,*) (cb(i),i = 1,nq)
write(6,*) (cb(i),i = 1,nq)
write(9,*) (cb(i),i = 1,nq)
8125 write(6,8130)
write(9,8130)
8130 format(/,1x,'How many values (5 or fewer) of sigma ',
x'do you want in the power table?',/)
read(4,*) numsig
write(6,*) numsig
write(9,*) numsig
write(6,8140)
write(9,8140)
8140 format(/,1x,'What are the sigma values?',/)
read(4,*) (sigma(i),i = 1,numsig)
write(6,*) (sigma(i),i = 1,numsig)
write(9,*) (sigma(i),i = 1,numsig)
write(6,8160)
write(9,8160)
8160 format(/,1x,'What power do you wish to achieve?',/)
read(4,*) powerl
write(6,*) powerl
write(9,*) powerl
if (ialready .eq. 1) go to 8170
call hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
x notestf,2)
iend = 0
if (notestf .eq. 1) go to 8000
8170 call power(nq,cb1,numcb,sigma,numsig,cb,powerl,iresdof)
iend = 0
go to 7005
c
c End of experimental design material
c
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)
write(9,10)
10 format(1x,'Matrix multiplication error:',/,
x 1x,'The a column and b row dimensions do not match.')
stop
endif
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
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)
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 ccheck(c,ncol,nq,notind,nrdc,numest)
c This subroutine checks whether the columns of c are
c linearly independent
c The parameter and dimension statements in the main, ccheck,
c hyptest, estim, and power subroutines must be in accord.
parameter (numpar=75,numestd=numpar)
implicit double precision (a-h,o-z)
dimension c(nrdc,numestd),uc(numpar,numestd),
xvc(numestd,numestd),
xxlamc(numestd),e(numestd)
dimension ctemp(numpar,numestd)
dimension work(numpar)
common /threshmc/ threshm
notind = 0
do 5 i = 1,ncol
do 4 j = 1,nq
ctemp(i,j) = c(i,j)
4 continue
5 continue
c
c ssvdc is the LINPACK singular value decomposition routine
c
call dsvdc(ctemp,numpar,ncol,nq,xlamc,e,uc,numpar,
xvc,numestd,work,21,info)
xnorm = 0.0
do 10 i = 1,nq
xnorm = xnorm + xlamc(i)**2
10 continue
xnorm = sqrt(xnorm)
thresh = xnorm*threshm
izerof = 0
do 20 i = 1,nq
if (xlamc(i) .gt. thresh) go to 20
ifirst = i
izerof = 1
go to 30
20 continue
30 if (izerof .eq. 0) return
irank = ifirst - 1
write(6,40)
write(9,40)
40 format(/,1x,'The columns of c = (c1 ... cq) appear',
x' to be linearly dependent.',/)
write(6,50) irank
write(9,50) irank
50 format(/,1x,'The rank of c appears to be ',i2,/)
write(6,60)
write(9,60)
60 format(/,1x,'The set of all d such that (c1 ... cq)d = 0',
x' (approximately)',/,1x,'is the space spanned by the ',
x'following vectors',//)
do 80 i = 1,nq
write(6,70) (vc(i,j),j = ifirst,nq)
write(9,70) (vc(i,j),j = ifirst,nq)
70 format(1x,10(e10.2,2x))
80 continue
notind = 1
return
end
subroutine hyptest(x,y,nrdx,nrow,ncol,c,nrdc,nq,hypss,ihypdof,
x notestf,idesign)
c
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
c The parameter and dimension statements in the main, ccheck,
c hyptest, estim, and power subroutines must be in accord.
parameter (nrowd=400,numpar=75,numest=numpar)
implicit double precision (a-h,o-z)
dimension x(nrdx,ncol),y(nrdx),c(nrdc,ncol),e(numpar),work(nrowd)
dimension f(numpar,numpar),ft(numpar,numpar),fft(numpar,numpar)
dimension vtemp(nrowd),temp1(nrowd),temp2(nrowd)
dimension g(nrowd,numpar),h(nrowd,numpar)
dimension ua(nrowd,numest),xlama(numest),va(numest,numest)
dimension notest(numest)
dimension at(numest,nrowd),ata(numest,numest)
dimension det(2)
common /hec/ ux(nrowd,numpar),vx(numpar,numpar),xlamx(numpar),
x a(nrowd,numest)
common /hec2/ resss,iresdof
common /ataic/ atai(numest,numest)
common /ixrankc/ ixrank
common /threshmc/ threshm
external vnorm
c
c notestf is set to 1 if some of the c's are not estimable
c
notestf = 0
do 5 i = 1,nrow
do 4 j = 1,ncol
g(i,j) = x(i,j)
4 continue
5 continue
c
c dsvdc is the LINPACK singular value decomposition routine
c
call dsvdc(g,nrdx,nrow,ncol,xlamx,e,ux,nrdx,vx,numpar,
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
7 format(1x,'ERROR hyptest 7: The info value from a singular',
x ' value decomposition was ',i3,'.',/)
stop
endif
c
c Obtain the threshold that is used to determine whether singular values
c are "essentially zero"
c
xnorm = 0.0
do 10 i = 1,ncol
xnorm = xnorm + xlamx(i)**2
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
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
do 52 j = 1,ncol
do 42 i = 1,ncol
fft(i,j) = 0.0
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
call vector(c,ncol,i,vtemp,nrdc)
thresh = vnorm(vtemp,ncol)*threshm
call matmult(fft,ncol,ncol,vtemp,ncol,1,temp1,numpar,nrowd)
call matsub(vtemp,temp1,temp2,ncol,1,nrowd)
if (vnorm(temp2,ncol) .lt. thresh) go to 200
mm = mm + 1
notest(mm) = i
200 continue
if (mm .eq. 0) go to 300
write(6,210) (notest(i),i = 1,mm)
write(9,210) (notest(i),i = 1,mm)
210 format(///,1x,"The following c's do not appear to be estimable",
1//,(1x,10(i3,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(9,302)
302 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,nrowd,numpar)
call matmult(h,nrow,ncol,c,ncol,nq,a,nrowd,nrdc)
if (idesign .eq. 0) go to 323
call mattran(a,nrow,nq,at,nrowd,numest)
call matmult(at,nq,nrow,a,nrow,nq,ata,numest,nrowd)
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.
c See the LINPACK manual
c for details. (Dongarra, Moler, Bunch, and Stewart,
c LINPACK Users' Guide, SIAM, 1979)
c
call dpofa(ata,numest,nq,info)
if (info .ne. 0) then
write(6,9999) info
9999 format(/,1x,'ERROR 9999: The info value from a matrix',
x ' inversion was',1x,i3,'.',/)
stop
endif
c
c Calculate the inverse of ata
c
job = 1
call dpodi(ata,numest,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
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
write(8,326)
326 format(//,1x,"The a's are (as columns)",/)
do 328 i = 1,nrow
write(8,327) (a(i,j),j = 1,nq)
327 format(10(1x,e11.3))
328 continue
c
c Perform a singular value decomposition of A
c in order to obtain an orthonormal basis of the
c linear span of A.
c dsvdc is a LINPACK singular value decomposition routine.
c
call dsvdc(g,nrowd,nrow,nq,xlama,e,ua,nrowd,va,numest,
xwork,20,infoa)
if (infoa .ne. 0) then
write(6,329) infoa
329 format(/,1x,'ERROR 329: The info value from a singular',
x ' value decomposition was ',i3,'.',/)
stop
endif
c
c obtain the hypothesis sum of squares
c
hypss = 0.0
do 340 j = 1,nq
temp0 = 0.0
do 330 i = 1,nrow
temp0 = temp0 + ua(i,j)*y(i)
330 continue
hypss = hypss + temp0**2
340 continue
ihypdof = nq
c
c obtain the residual sum of squares
c
xss = 0.0
do 440 j = 1,ixrank
temp0 = 0.0
do 430 i = 1,nrow
temp0 = temp0 + ux(i,j)*y(i)
430 continue
xss = xss + temp0**2
440 continue
totss = 0.0
do 450 i = 1,nrow
totss = totss + y(i)**2
450 continue
resss = totss - xss
return
end
subroutine estim(x,y,nrdx,nrow,ncol,c,nrdc,nc,alphai,alphaj,
x est,cli,cui,clj,cuj,notestf)
c
c This subroutine checks whether the c's are estimable and
c if they are it obtains confidence intervals.
c
c The parameter and dimension statements in the main, ccheck,
c hyptest, estim, and power subroutines must be in accord.
parameter (nrowd=400,numpar=75,numest=numpar)
implicit double precision (a-h,o-z)
dimension x(nrdx,ncol),y(nrdx),c(nrdc,numest),e(numpar),
xwork(nrowd)
dimension f(numpar,numpar),ft(numpar,numpar),fft(numpar,numpar)
dimension vtemp(nrowd),temp1(nrowd),temp2(nrowd)
dimension g(nrowd,numpar),h(nrowd,numpar)
dimension est(numest),cli(numest),cui(numest),clj(numest),
xcuj(numest)
dimension notest(numest)
c
c Notice the restrictions on the number of c's
c
dimension xlamc(numpar),uc(numpar,numpar),vc(numpar,numpar)
common /hec/ ux(nrowd,numpar),vx(numpar,numpar),xlamx(numpar),
x a(nrowd,numest)
common /hec2/ resss,iresdof
common /ixrankc/ ixrank
common /threshmc/ threshm
external vnorm
notestf = 0
do 5 i = 1,nrow
do 4 j = 1,ncol
g(i,j) = x(i,j)
4 continue
5 continue
c
c ssvdc is the LINPACK singular value decomposition routine
c
call dsvdc(g,nrowd,nrow,ncol,xlamx,e,ux,nrowd,vx,numpar,work,
x21,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
7 format(1x,'ERROR estim 7: The info value from a singular',
x ' value decomposition was ',i3,'.',/)
stop
endif
c
c Obtain the threshold that is used to determine whether singular values
c are "essentially zero"
c
xnorm = 0.0
do 10 i = 1,ncol
xnorm = xnorm + xlamx(i)**2
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 1000
120 continue
ixrank = ncol
c
c Determine the number of linearly independent c's
c
1000 if (ncol .gt. nc) then
do 1003 i = 1,nrowd
do 1002 j = 1,numpar
g(i,j) = 0.0
1002 continue
1003 continue
c
c C is made up of columns of c's
c
do 1005 i = 1,ncol
do 1004 j = 1,nc
g(i,j) = c(i,j)
1004 continue
1005 continue
c
c ssvdc is the LINPACK singular value decomposition routine
c
call dsvdc(g,nrowd,ncol,nc,xlamc,e,uc,numpar,vc,numpar,work,
x 21,info)
c
c Obtain the threshold that is used to determine whether singular values
c are "essentially zero"
c
xnorm = 0.0
do 1010 i = 1,nc
xnorm = xnorm + xlamc(i)**2
1010 continue
xnorm = sqrt(xnorm)
thresh = xnorm*threshm
c
c Determine the rank of C
c
do 1120 i = 1,nc
if (xlamc(i) .gt. thresh) go to 1120
icrank = i - 1
go to 130
1120 continue
icrank = nc
else
c
c nc >= ncol
c
do 2003 i = 1,nrowd
do 2002 j = 1,numpar
g(i,j) = 0.0
2002 continue
2003 continue
do 2005 i = 1,ncol
do 2004 j = 1,nc
g(j,i) = c(i,j)
2004 continue
2005 continue
c
c ssvdc is the LINPACK singular value decomposition routine
c
call dsvdc(g,nrowd,nc,ncol,xlamc,e,uc,numpar,vc,numpar,work,
x 21,info)
c
c Obtain the threshold that is used to determine whether singular values
c are "essentially zero"
c
xnorm = 0.0
do 2010 i = 1,ncol
xnorm = xnorm + xlamc(i)**2
2010 continue
xnorm = sqrt(xnorm)
thresh = xnorm*threshm
c
c Determine the rank of C
c
do 2120 i = 1,ncol
if (xlamc(i) .gt. thresh) go to 2120
icrank = i - 1
go to 130
2120 continue
icrank = ncol
endif
c
c Determine whether the c's are estimable (there exists an a
c such that X'a = c)
c
130 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
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
do 52 j = 1,ncol
do 42 i = 1,ncol
fft(i,j) = 0.0
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,nc
call vector(c,ncol,i,vtemp,nrdc)
thresh = vnorm(vtemp,ncol)*threshm
call matmult(fft,ncol,ncol,vtemp,ncol,1,temp1,numpar,nrowd)
call matsub(vtemp,temp1,temp2,ncol,1,nrowd)
if (vnorm(temp2,ncol) .lt. thresh) go to 200
mm = mm + 1
notest(mm) = i
200 continue
if (mm .eq. 0) go to 300
write(6,210) (notest(i),i = 1,mm)
write(9,210) (notest(i),i = 1,mm)
210 format(///,1x,"The following c's do not appear to be estimable",
x//,(1x,10(i2,2x)),/)
notestf = 1
return
c
c Calculate the ai's that satisfy X'ai = ci
c
300 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,nrowd,numpar)
call matmult(h,nrow,ncol,c,ncol,nc,a,nrowd,nrdc)
c
c Print the a vectors
c
write(8,326)
326 format(//,1x,"The a's are (as columns)",/)
do 328 i = 1,nrow
write(8,327) (a(i,j),j = 1,nc)
327 format(10(1x,e11.3))
328 continue
c
c Obtain the residual vector
c
xss = 0.0
do 440 j = 1,ixrank
temp0 = 0.0
do 430 i = 1,nrow
temp0 = temp0 + ux(i,j)*y(i)
430 continue
xss = xss + temp0**2
440 continue
totss = 0.0
do 450 i = 1,nrow
totss = totss + y(i)**2
450 continue
resss = totss - xss
iresdof = nrow - ixrank
c
c Calculate the individual and the joint confidence intervals
c
500 s = sqrt(resss/iresdof)
xxx = icrank
yyy = 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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.html
c
p = 1.0 - alphai/2.0
q = 1.0d0 - p
call cdft(2,p,q,t,yyy,istatus,bound)
if (istatus .ne. 0) then
write(6,502) istatus
502 format(/,1x,'ERROR 502: The istatus value from a T inverse',
x ' routine was ',i2,'.',/)
stop
endif
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.html
c
p = 1.0 - alphaj
q = alphaj
call cdff(2,p,q,fvalue,xxx,yyy,istatus,bound)
if (istatus .ne. 0) then
write(6,504) istatus
504 format(/,1x,'Error 504: The istatus value from an',
x ' F inverse',
x ' routine was ',i2,'.',/)
stop
endif
delta1 = s*t
delta2 = s*sqrt(icrank*fvalue)
do 600 i = 1,nc
aa = 0.0
xxx = 0.0
do 550 j = 1,nrow
aa = aa + a(j,i)*a(j,i)
xxx = xxx + a(j,i)*y(j)
550 continue
aa = sqrt(aa)
est(i) = xxx
delta = aa*delta1
cli(i) = xxx - delta
cui(i) = xxx + delta
delta = aa*delta2
clj(i) = xxx - delta
cuj(i) = xxx + delta
600 continue
return
end
subroutine power(nq,cb1,numcb,sigma,numsig,cb,powerl,iresdof)
c
c This subroutine performs power calculations.
c
c The parameter and dimension statements in the main, ccheck,
c hyptest, estim, and power subroutines must be in accord.
parameter (numest=75)
implicit double precision (a-h,o-z)
dimension cb1(5),sigma(5),cb(numest)
dimension powera(5,5,3)
dimension numrep(5,5,3)
dimension halfwidt(3,5)
common /alphac/ alpha(3)
common /ataic/ atai(numest,numest)
common /nrowc/ nrow
common /ixrankc/ ixrank
if (nq .gt. 1) go to 1000
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.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
10 format(/,1x,'ERROR 10: The istatus value from a T inverse',
x ' routine was ',i2,'.',/)
stop
endif
do 100 i = 1,numsig
do 90 j = 1,numcb
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.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
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.0)
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.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
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.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
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.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
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.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
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(9,205)
205 format(///,1x,"Power tables",//,1x,"(Here q = 1, and we consider",
x" a range of c'b values and a range of sigma values.)",/)
do 300 ii = 1,3
write(6,310) alpha(ii)
write(9,310) alpha(ii)
310 format(///,1x,'The significance level is ',f4.2,//)
write(6,320)
write(9,320)
320 format(/,1x,"The c'b values are on the top. The sigma values",
x " are along the left side.",//)
write(6,330) (cb1(i),i = 1,numcb)
write(9,330) (cb1(i),i = 1,numcb)
330 format(/,18x,5(e12.4,3x))
write(6,332)
write(9,332)
332 format(/)
do 340 i = 1,numsig
write(6,335) sigma(i),(powera(i,j,ii),j = 1,numcb)
write(9,335) sigma(i),(powera(i,j,ii),j = 1,numcb)
335 format(1x,e12.4,5x,5(e12.4,3x))
340 continue
300 continue
write(6,355)
write(9,355)
355 format(///,1x,"Halfwidth tables",//,1x,"(Here q = 1,",
x" and we consider",
x" a range of sigma values.)",/)
do 390 ii = 1,3
write(6,360) alpha(ii)
write(9,360) alpha(ii)
360 format(///,1x,'The significance level is ',f4.2,//)
write(6,365)
write(9,365)
365 format(/,1x,"The sigma values",
x " are along the left side.",//)
do 375 i = 1,numsig
write(6,370) sigma(i),halfwidt(ii,i)
write(9,370) sigma(i),halfwidt(ii,i)
370 format(1x,e12.4,5x,5(e12.4,3x))
375 continue
390 continue
write(6,405)
write(9,405)
405 format(///,1x,'Replication tables',//,
x1x,'These tables give us the number of replications of the ',
x'current design',/,1x,'that are needed to achieve the desired',
x' power.',/)
write(6,410) powerl
write(9,410) powerl
410 format(/,1x,'The desired power is ',f4.2,//)
do 500 ii = 1,3
write(6,420) alpha(ii)
write(9,420) alpha(ii)
420 format(///,1x,'The significance level is ',f4.2,//)
write(6,430)
write(9,430)
430 format(/,1x,"The c'b values are on the top. the sigma values",
x " are along the left side.",//)
write(6,440) (cb1(i),i = 1,numcb)
write(9,440) (cb1(i),i = 1,numcb)
440 format(/,18x,5(e12.4,3x))
write(6,442)
write(9,442)
442 format(/)
do 450 i = 1,numsig
write(6,445) sigma(i),(numrep(i,j,ii),j = 1,numcb)
write(9,445) sigma(i),(numrep(i,j,ii),j = 1,numcb)
445 format(1x,e12.4,5x,5(5x,i2,8x))
450 continue
500 continue
return
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.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(/,1x,'Error 1005: The istatus value from an',
x ' F inverse',
x ' routine was ',i2,'.',/)
stop
endif
do 1100 i = 1,numsig
xncp = 0.0
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)**2)
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.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(/,1x,'Error 1025: The istatus value from',
x ' a non-central F cdf',
x ' routine was ',i2,'.',/)
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.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(/,1x,'Error 1035: The istatus value from an',
x ' F inverse',
x ' routine was ',i2,'.',/)
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.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(/,1x,'Error 1040: The istatus value from',
x ' a non-central F cdf',
x ' routine was ',i2,'.',/)
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.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(/,1x,'Error 1062: The istatus value from an',
x ' F inverse',
x ' routine was ',i2,'.',/)
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 The library can be found at
c
c http://odin.mdacc.tmc.edu/anonftp/source.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(/,1x,'Error 1064: The istatus value from',
x ' a non-central F cdf',
x ' routine was ',i2,'.',/)
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,1205)
write(9,1205)
1205 format(///,1x,'Power tables',//)
write(6,1210) (cb(i),i = 1,nq)
write(9,1210) (cb(i),i = 1,nq)
1210 format(/,1x,"The c1'b, ... , cq'b values are",//,1x,
x5(e12.4,2x),/,1x,5(e12.4,2x),/)
do 1300 ii = 1,3
write(6,1310) alpha(ii)
write(9,1310) alpha(ii)
1310 format(///,1x,'The significance level is ',f4.2,//)
write(6,1320)
write(9,1320)
1320 format(/,1x,'The sigma and associated power values are',//)
do 1340 i = 1,numsig
write(6,1335) sigma(i),(powera(i,1,ii))
write(9,1335) sigma(i),(powera(i,1,ii))
1335 format(1x,e12.4,5x,e12.4)
1340 continue
1300 continue
write(6,1405)
write(9,1405)
1405 format(///,1x,'Replication tables',//,
x1x,'These tables give us the number of replications of the ',
x'current design',/,1x,'that are needed to achieve the desired',
x' power.',/)
write(6,1410) powerl
write(9,1410) powerl
1410 format(/,1x,'The desired power is ',f4.2,//)
do 1500 ii = 1,3
write(6,1420) alpha(ii)
write(9,1420) alpha(ii)
1420 format(///,1x,'The significance level is ',f4.2,//)
write(6,1430)
write(9,1430)
1430 format(/,1x,'The sigma and associated',
x ' replication values are',//)
do 1450 i = 1,numsig
write(6,1445) sigma(i),numrep(i,1,ii)
write(9,1445) sigma(i),numrep(i,1,ii)
1445 format(1x,e12.4,10x,i2)
1450 continue
1500 continue
return
end