Thread View: bit.listserv.stat-l
1 messages
1 total messages
Started by Jerry Dallal
Wed, 28 Jul 1993 15:55
Intra-class correlation coefficient (LONG)
Author: Jerry Dallal
Date: Wed, 28 Jul 1993 15:55
Date: Wed, 28 Jul 1993 15:55
559 lines
20538 bytes
20538 bytes
ICC
A Program for Sample Size Determinations
for the Intra-Class Correlation Coefficient
Version 0.1
July 27, 1993
"One of many STATOOLS(tm)..."
by
Gerard E. Dallal
54 High Plain Road
Andover, MA 01810
NOTICE
Documentation and original code copyright 1993 by
Gerard E. Dallal. Reproduction of material for non-
commercial purposes is permitted, without charge, provided
suitable reference is made to ICC and its author.
Neither ICC nor its documentation should be modified in
any way without permission from the author, except for those
changes that are essential to move ICC to another
computer.
Please acknowledge ICC in any manuscript that uses its
calculations.
DISCLAIMER
STATOOLS are provided "as is" without warranty of any
kind. The entire risk as to the quality, performance, and
fitness for intended purpose is with you. You assume
responsibility for the selection of the program and for the
use of results obtained from that program.
DESCRIPTION
ICC calculates the sample size necessary so that, with
specified probability, the lower limit of a one-tailed
confidence interval for an intra-class correlation coefficient
based on a specified number or measurements per observational
unit will exceed a specified value for a specified underlying
value of the intra-class correlation and a specified number of
measurements per observational unit.
Example: ICC will calculate the sample size necessary to
have an 80% chance that the lower limit of a one-tailed 95-%
confidence interval for an intra-class correlation coefficient
(icc) will exceed 0.40 when the underlying icc is 0.80 and 3
measurements are made on each observational unit.
This formulation can be stated in terms of tests of
significance. This sample size is the sample size required
for a one-tailed fixed-level test of the null hypothesis that
the intra-class correlation coefficient is less than or equal
to some specified constant to have specified power for a
specified alternative.
Example: The earlier sample size will give 80-% power
that an 0.05 level test will reject the null hypothesis that
the icc is less than or equal to 0.40 when the underlying
value is 0.80 and 3 measurements are made on each
observational unit.
The theory behind the sample size estimation can be found
in Scheffe (1959, section 7.2). The critical result is that
the usual F-ratio (MS Between / MS Within) divided by (1+n r),
where n is the number of measurements per observational unit
and r is the ratio of the between units variance to the within
units variance, follows an F distribution with numerator
degrees of freedom equal to one less than the number of
observational units and denominator degrees of freedom equal
to the number of observational units multiplied by one less
than the number of measurements per unit.
The application of the result is straightforward, once an
expression involving the icc is recast in terms of r
(icc = r / (1+r).
Donner and Eliasziw (1987) give power contours for some
common special cases of the testing formulation (level = 0.05;
power = 0.80; lower limit = 0.0(0.2)0.8; axes: number of
measurements per observational unit, number of units; contours
for 3 to 5 values of the underlying icc).
OPERATION
ICC is prompt driven. Quantities appearing in brackets
(e.g., [0.800]) are defaults that can be obtained by merely
pressing the Enter key.
PROGRAMMER'S NOTES
The terminal $'s in FORMAT statements are used to prevent
the cursor from moving to the next line when input is
requested. This is standard VAX/VMS FORTRAN syntax.
Microsoft FORTRAN specifies the backslash for this purpose,
but it provides undocumented support for the $, too.
The program is set up for the IBM-PC. To move ICC to
another platform, it should be necessary only to reset the
input/output unit numbers in the first DATA statement and to
rewrite the OPEN statement for the output file. The
modifications required by VAX/VMS FORTRAN are contained in
comments in the source code.
ALGORITHMS
ICC makes use of the following published routines,
modified to run in double precision:
Cran GW, Martin KJ and Thomas GE (1977), "Remark AS R19 and
Algorithm AS 109. A Remark on Algorithms AS 63: The
Incomplete Beta Integral, and AS 64: Inverse of the
Incomplete Beta Function Ratio," Applied Statistics, 26,
111-114.
Majumder KL and Bhattacharjee GP (1973), "Algorithm AS 63.
The Incomplete Beta Integral," Applied Statistics, 22,
409-411.
and my FORTRAN translation of
Pike MC and Hill ID (1966), "Algorithm 291. Logarithm of the
Gamma Function," Communications of the Association for
Computing Machinery, 9, 684.
REFERENCES
Donner A and Eliasziw M (1987), "Sample Size Requirements for
Reliability Studies," Statistics in Medicine, 6, 441-448.
Scheffe H (1959), The Analysis of Variance. New York: John
Wiley and Sons, Inc.
C***********************************************************************
C ICC
C A Program for Sample Size Determinations
C for the Intra-Class Correlation Coefficient
C Version 0.1
C July 27, 1993
C
C "One of many STATOOLS(tm)..."
C by
C
C Gerard E. Dallal
C 54 High Plain Road
C Andover, MA 01810
C
C
C NOTICE
C
C Documentation and original code copyright 1993 by
C Gerard E. Dallal. Reproduction of material for non-
C commercial purposes is permitted, without charge, provided
C suitable reference is made to ICC and its author.
C
C Neither ICC nor its documentation should be modified in
C any way without permission from the author, except for those
C changes that are essential to move ICC to another
C computer.
C
C Please acknowledge ICC in any manuscript that uses its
C calculations.
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
LOGICAL QFILE
CHARACTER*1 QUERY
CHARACTER FNAME*40, XCORR*18
C IIN -- input unit number (screen)
C IOUT -- output unit number (screen)
C IWOUT -- saved output unit number
C Open statement located just before statement 50
DATA IIN /0/, IOUT /0/, IWOUT /11/, QFILE /.FALSE./
C For VAX/VMS FORTRAN
C DATA IIN /5/, IOUT /6/, IWOUT /11/, QFILE /.FALSE./
C Initialize choices
ALPHA = 0.05D0
RPOWER = 0.80D0
ECORR = 0.80D0
CORR = 0.60D0
NREP = 2
C
WRITE(IOUT,10)
10 FORMAT (//38X,'ICC'/
* 20X,'A Program for Sample Size Determinations'/
* 18X,'for the Intra-Class Correlation Coefficient'/
* 32X,'Version 0.10'/35X,'(c) 1993'//
* 26X,'"One of many STATOOLS(tm)..."'/
* 38X,'by'//
* 31X,'Gerard E. Dallal'/
* 30X,'54 High Plain Road'/
* 30X,'Andover, MA 01810'///
* 19X,'Please acknowledge ICC in any publication'/
* 22X,'that makes use of its calcluations.')
C
WRITE(IOUT,20)
20 FORMAT(//' Do you wish to save the results of this session? ',
* '[N]: '$)
READ (IIN,30) QUERY
30 FORMAT (A1)
IF (QUERY.EQ.'Y' .OR. QUERY.EQ.'y') QFILE = .TRUE.
IF (.NOT. QFILE) GOTO 50
WRITE (IOUT,40)
40 FORMAT (' Enter output filename: '$)
READ (IIN,'(A)') FNAME
OPEN (UNIT = IWOUT, FILE = FNAME, STATUS = 'UNKNOWN')
C For VAX/VMS FORTRAN
C OPEN (UNIT = IWOUT, FILE = FNAME, STATUS = 'UNKNOWN',
C * CARRIAGECONTROL = 'LIST')
C
50 WRITE (IOUT,60) 1.0D0 - ALPHA
60 FORMAT (/33X,'Enter confidence level [',F5.3,']: '$)
READ (IIN,70,ERR=50) DALPHA
70 FORMAT (BN,F18.0)
IF (DALPHA.LT.0.0D0 .OR. DALPHA.GE.1.0D0) GOTO 50
IF (DALPHA.NE.0.0D0) ALPHA = 1.0D0 - DALPHA
C Intra-class correlation cofficient
80 WRITE (IOUT,90) ECORR
90 FORMAT (
* ' Enter expected intra-class correlation coefficient [',
* F5.3,']: '$)
READ (IIN,70,ERR=80) DECORR
IF (DECORR.LT.0.0D0 .OR. DECORR.GE.1.D0) GOTO 80
IF (DECORR.NE.0.0D0) ECORR = DECORR
C
100 WRITE (IOUT,110) CORR
110 FORMAT (
* 14X,'Enter lower limit for confidence interval [',F5.3,']: '$)
READ (IIN,'(A)') XCORR
IF (XCORR.NE.' ') THEN
READ (XCORR,70,ERR=100) DCORR
IF (DCORR.LT.0.0D0 .OR. DCORR.GE.1.D0) GOTO 100
CORR = DCORR
ENDIF
IF (CORR.GE.ECORR) THEN
WRITE (IOUT,120)
120 FORMAT (' Lower limit must be less than expected icc.')
GOTO 80
ENDIF
C Convert to theta's
ETH = ECORR / (1.0D0 - ECORR)
TH = CORR / (1.0D0 - CORR)
C
130 WRITE (IOUT,140) RPOWER
140 FORMAT (' Enter probability'/
* 13X,'that lower limit is no less than specified [',F5.3,']: '$)
READ (IIN,70,ERR=130) DPOWER
IF (DPOWER.LT.0.0D0 .OR. DPOWER.GE.1.D0) GOTO 130
IF (DPOWER.NE.0.0D0) RPOWER = DPOWER
C
150 WRITE (IOUT,160) NREP
160 FORMAT (
* ' Enter number of measurements per observational unit [',
* I5,']: '$)
READ (IIN,70,ERR=150) XNREP
IF (XNREP.LT.0.0D0 .OR. XNREP.EQ.1.0D0 .OR.
* XNREP.NE.DINT(XNREP)) GOTO 150
IF (XNREP.NE.0.0D0) NREP = XNREP
FNREP = NREP
C
C = (1.0D0 + FNREP*TH) / (1.0D0 + FNREP*ETH)
C NS = 2
NS = 2
F2 = 2 * (NREP - 1)
XALPHA = XINBTA(0.5D0*F2, 0.5D0, ALPHA)
FALPHA = F2 * ((1.0D0 - XALPHA) / XALPHA)
POWER = FDFC(FALPHA*C, 1.0D0, F2)
IF (POWER.GE.RPOWER) GOTO 200
C Find sample sizes that bracket required power
180 NS = 2 * NS
F1 = NS - 1
F2 = NS * (NREP - 1)
AA = F2 / 2.0D0
BB = F1 / 2.0D0
XALPHA = XINBTA(AA, BB, ALPHA)
FALPHA = (F2 / F1) *((1.0D0 - XALPHA) / XALPHA)
POWER = FDFC(FALPHA*C, F1, F2)
IF (POWER.LT.RPOWER) GOTO 180
C Binary search
NHI = NS
HIPOWR = POWER
NLO = NS / 2
190 NS = (NLO + NHI) / 2
C
F1 = NS - 1
F2 = NS * (NREP - 1)
AA = F2 / 2.0D0
BB = F1 / 2.0D0
XALPHA = XINBTA(AA, BB, ALPHA)
FALPHA = (F2 / F1) *((1.0D0 - XALPHA) / XALPHA)
POWER = FDFC(FALPHA*C, F1, F2)
IF (POWER.LT.RPOWER) THEN
NLO = NS
ELSE
NHI = NS
HIPOWR = POWER
ENDIF
IF (NHI - NLO.GT.1) GOTO 190
C
200 WRITE (IOUT,210) NHI, 100.0D0-100.0D0*ALPHA,CORR, HIPOWR
210 FORMAT (/1X,I12,' units: ',
* 'Pr (',F4.1,'-% CI lies above',F5.3,') =',F7.5)
C
IF (QFILE) THEN
WRITE(IWOUT,220) ECORR
220 FORMAT (//
* 12X,'Expected intra-class correlation coefficient:',F6.3)
WRITE(IWOUT,230) NREP
230 FORMAT (12X,
* 'Number of measurements per observational unit:',I6)
WRITE (IWOUT,240) NHI, 100.0D0-100.0D0*ALPHA,CORR, HIPOWR
240 FORMAT (1X,I15,' units: ',
* 'Pr (',F4.1,'-% CI lies above',F5.3,') =',F7.5)
ENDIF
C
WRITE (IOUT,250)
250 FORMAT (/' Do you wish to continue? [Y]: '$)
READ (IIN,30) QUERY
IF (QUERY.EQ.'N' .OR. QUERY.EQ.'n') STOP ' '
GOTO 50
C
END
C***********************************************************************
DOUBLE PRECISION FUNCTION ALGAMA(S)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C CACM Algorithm 291
C Logarithm of gamma function
C by Pike, M.C. and Hill, I.D.
C
C Translated into FORTRAN by
C Gerard E. Dallal
C 54 High Plain Road
C Andover, MA 01810
C
X = S
ALGAMA = 0.0D0
IF (X .LE. 0.0D0) RETURN
F = 0.0D0
IF (X .GE. 7.0D0) GOTO 30
C
F = 1.0D0
Z = X - 1.0D0
C
10 Z = Z + 1.0D0
IF (Z .GE. 7.0D0) GOTO 20
X = Z
F = F * Z
GOTO 10
C
20 X = X + 1.0D0
F = -DLOG(F)
C
30 Z = 1.0D0 / X ** 2
ALGAMA = F + ( X - 0.5D0) * DLOG(X) - X + 0.918938533204673D0 +
1 (((-0.000595238095238D0 * Z + 0.000793650793651D0) * Z
2 - 0.002777777777778D0) * Z + 0.083333333333333D0) / X
RETURN
END
C***********************************************************************
DOUBLE PRECISION FUNCTION BETAIN(X,P,Q)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C Algorithm AS 63 Appl.Statist. (1973), VOL.22, NO.3
C Modified as per remark ASR 19 (1977), VOL.26, NO.1
C ***[Fault indicator removed by G.E.D.]***
C
C Computes incomplete beta function ratio for arguments
C X between 0 and 1, P and Q positive.
C Log of complete beta function, BETA, assumed to be known.
C ***[Calculation of BETA added by G.E.D]***
C
LOGICAL INDEX
C
C DEFINE ACCURACY AND INITIALIZE
C
DATA ACU /0.1D-7/
BETAIN = X
C
C TEST FOR ADMISIBILITY OF AGRUMENTS
C
IF (P .LE. 0.0D0 .OR. Q .LE. 0.0D0) STOP 40
IF (X .LT. 0.0D0 .OR .X .GT. 1.0D0) STOP 41
IF (X .EQ .0.0D0 .OR .X .EQ. 1.0D0) RETURN
C
C CHANGE TAIL IF NECESSARY AND DETERMINE S
C
PSQ = P + Q
CX = 1.0D0 - X
IF (P .GE .PSQ * X) GOTO 1
XX = CX
CX = X
PP = Q
QQ = P
INDEX = .TRUE.
GOTO 2
1 XX = X
PP = P
QQ = Q
INDEX = .FALSE.
2 TERM = 1.0
AI = 1.0
BETAIN = 1.0
NS = QQ + CX * PSQ
C
C USE SOPER'S REDUCTION FORMULAE
C
RX = XX / CX
3 TEMP = QQ - AI
IF (NS .EQ. 0) RX = XX
4 TERM = TERM * TEMP * RX / (PP + AI)
BETAIN = BETAIN + TERM
TEMP = ABS(TERM)
IF (TEMP .LE. ACU .AND. TEMP .LE. ACU * BETAIN) GOTO 5
AI = AI + 1.0D0
NS = NS - 1
IF (NS .GE. 0) GOTO 3
TEMP = PSQ
PSQ = PSQ + 1.0D0
GOTO 4
C
C CALCULATE RESULT
C
5 BETA = ALGAMA(P) + ALGAMA(Q) - ALGAMA(P+Q)
BETAIN = BETAIN *
* DEXP(PP * DLOG(XX) + (QQ - 1.0D0) * DLOG(CX) - BETA) / PP
IF (INDEX) BETAIN = 1.0D0 - BETAIN
RETURN
END
C***********************************************************************
DOUBLE PRECISION FUNCTION XINBTA(P, Q, ALPHA)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C Algorithm AS 109 Appl.Statist. (1977), Vol.26, NO.1
C (Replacing algorithm AS 64 Appl. Statist. (1973),
C Vol.22, NO.3)
C ***[Fault indicator removed by G.E.D.]***
C
C Computes inverse of incomplete beta function
C ratio for given positive values of the arguments
C P and Q, ALPHA between 0 and 1.
C Log of complete beta function, BETA, is assumed to be known.
C ***[Calculation of BETA added by G.E.D]***
C
LOGICAL INDEX
C
C DEFINE ACCURACY AND INITIALIZE
C
DATA ACU /1.0D-14/
XINBTA = ALPHA
C
C TEST FOR ADMISIBILITY OF PARAMETERS
C
IF (P.LE.0.0D0.OR.Q.LE.0.0D0) STOP 50
IF (ALPHA.LT.0.0D0.OR.ALPHA.GT.1.0D0) STOP 51
IF (ALPHA.EQ.0.0D0.OR.ALPHA.EQ.1.0D0) RETURN
C
BETA=ALGAMA(P)+ALGAMA(Q)-ALGAMA(P+Q)
C
C CHANGE TAIL IF NECESSARY
C
IF (ALPHA .LE. 0.5D0) GOTO 1
A = 1.0D0 - ALPHA
PP = Q
QQ = P
INDEX = .TRUE.
GOTO 2
1 A = ALPHA
PP = P
QQ = Q
INDEX = .FALSE.
C
C CALCULATE INITIAL APPROXIMATION
C
2 R = DSQRT(-DLOG(A * A))
Y = R - (2.30753D0 + 0.27061D0 * R) /
* (1.0D0 + (0.99229D0 + 0.04481D0 * R) * R)
IF (PP .GT. 1.0D0 .AND. QQ .GT. 1.0D0) GOTO 5
R = QQ + QQ
T = 1.0D0 / (9.0D0 * QQ)
T = R * (1.0D0 - T + Y * DSQRT(T)) ** 3
IF (T .LE. 0.0D0) GOTO 3
T = (4.0D0 * PP + R - 2.0D0) / T
IF (T .LE. 1.0D0) GOTO 4
XINBTA = 1.0D0 - 2.0D0 / (T + 1.0D0)
GOTO 6
3 XINBTA = 1.0D0 - DEXP((DLOG((1.0D0 - A) * QQ) + BETA) / QQ)
GOTO 6
4 XINBTA = DEXP((DLOG(A * PP) + BETA) / PP)
GOTO 6
5 R = (Y * Y - 3.0D0) / 6.0D0
S = 1.0D0 / (PP + PP - 1.0D0)
T = 1.0D0 / (QQ + QQ - 1.0D0)
H = 2.0D0 / (S + T)
W = Y * DSQRT(H + R) / H - (T - S) *
* (R + 5.0D0 / 6.0D0 - 2.0D0 / (3.0D0 * H))
XINBTA = PP / (PP + QQ * DEXP(W + W))
C
C SOLVE FOR X BY A MODIFIED NEWTON-RAPHSON METHOD,
C USING THE FUNCTION BETAIN.
C
6 R = 1.0D0 - PP
T = 1.0D0 - QQ
YPREV = 0.0
SQ = 1.0
PREV = 1.0
IF (XINBTA .LT. 0.0001D0) XINBTA = 0.0001D0
IF (XINBTA .GT. 0.9999D0) XINBTA = 0.9999D0
7 Y = BETAIN(XINBTA, PP, QQ)
Y = (Y - A) * DEXP(BETA +
* R * DLOG(XINBTA) + T * DLOG(1.0D0 - XINBTA))
IF (Y * YPREV .LE. 0.0D0) PREV = SQ
G = 1.0
9 ADJ = G * Y
SQ = ADJ * ADJ
IF (SQ .GE. PREV) GOTO 10
TX = XINBTA - ADJ
IF (TX .GE. 0.0D0 .AND. TX .LE. 1.0D0) GOTO 11
10 G = G / 3.0D0
GOTO 9
11 IF (PREV .LE. ACU) GOTO 12
IF (Y * Y .LE. ACU) GOTO 12
IF (TX .EQ. 0.0D0 .OR. TX .EQ. 1.0D0) GOTO 10
IF (TX .EQ. XINBTA) GOTO 12
XINBTA = TX
YPREV = Y
GOTO 7
12 IF (INDEX) XINBTA = 1.0D0 - XINBTA
RETURN
END
C***********************************************************************
DOUBLE PRECISION FUNCTION FDFC(FQUAN, DFN, DFD)
C The complement of the cdf of the F distribution
C with DFN numerator degrees of freedom and DFD denominator
C degrees of freedom evaluated at FQUAN.
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
ESQ = DFD / (DFD + DFN * FQUAN)
FDFC = BETAIN (ESQ, DFD/2.0D0, DFN / 2.0D0)
RETURN
END
C***********************************************************************
Thread Navigation
This is a paginated view of messages in the thread with full content displayed inline.
Messages are displayed in chronological order, with the original post highlighted in green.
Use pagination controls to navigate through all messages in large threads.
Back to All Threads