Attribute VB_Name = "tetra_corr" 'Option Explicit 'tetrachoric correlation を求めるマクロ(tetra_corr)と関数(tetracorr) '関数は =tetracorr(a1:b2)のように範囲を指定する。 'マクロは相関係数以外の出力もする。 '範囲は2×2のクロス表(必ず頻度) 'http://ourworld.compuserve.com/homepages/jsuebersax/tetra.htm 'にあるJohn S. Uebersax氏のプログラムをvba 化した。2000/11/24日採取 'Keizo Hori 2000/12/4 Public Sub tetra_corr() Dim x As Range Dim n As Integer Set x = Application.InputBox(prompt:="2x2行列指定", Title:="tetrachoric correlation", Type:=8) n = x.Rows.Count If n <= 0 Then MsgBox "範囲を指定してください" Exit Sub ElseIf n <> 2 Or x.Columns.Count <> 2 Then MsgBox "2x2 の範囲を指定してください" Exit Sub End If Dim A As Double, B As Double, C As Double, D As Double, r As Double, SDR As Double, SDZERO As Double, ITYPE As Integer, IFAULT As Integer A = x(1, 1) B = x(1, 2) C = x(2, 1) D = x(2, 2) Call tetra(A, B, C, D, r, SDR, SDZERO, ITYPE, IFAULT) If (IFAULT = 1) Then MsgBox "<> Routine did not converge. Try increasing the maximum number of iterations (NITER) in Subroutine TETRA (see Brown, 1977)." Exit Sub End If If (IFAULT = 2) Then MsgBox "<> A negative frequency entered, or a row or column marginal is 0 (see Brown, 1977)." Exit Sub End If 'MsgBox = variant(r) ActiveCell.Value = "tetrachoric r = " ActiveCell.Offset(0, 1).Activate ActiveCell.Value = r ActiveCell.Offset(1, -1).Activate ActiveCell.Value = "std error =" ActiveCell.Offset(0, 1).Activate ActiveCell.Value = SDR Dim total As Double total = A + B + C + D ActiveCell.Offset(1, -1).Activate ActiveCell.Value = "threshold for item 1 = " ActiveCell.Offset(0, 1).Activate ActiveCell.Value = PPND((A + B) / total, 1) ActiveCell.Offset(1, -1).Activate ActiveCell.Value = "threshold for item 2 = " ActiveCell.Offset(0, 1).Activate ActiveCell.Value = PPND((A + C) / total, 1) ActiveCell.Offset(1, -1).Activate ActiveCell.Value = "phi coefficient= " ActiveCell.Offset(0, 1).Activate ActiveCell.Value = (A * D - B * C) / Sqr((A + B) * (C + D) * (A + C) * (B + D)) End Sub Function tetracorr(x As Range) As Variant Dim n As Integer n = x.Rows.Count If n <= 0 Then tetracorr = "範囲を指定してください" Exit Function ElseIf n <> 2 Or x.Columns.Count <> 2 Then tetracorr = "2x2 の範囲を指定してください" Exit Function End If Dim A As Double, B As Double, C As Double, D As Double, r As Double, SDR As Double, SDZERO As Double, ITYPE As Integer, IFAULT As Integer A = x(1, 1) B = x(1, 2) C = x(2, 1) D = x(2, 2) Call tetra(A, B, C, D, r, SDR, SDZERO, ITYPE, IFAULT) tetracorr = r If (IFAULT = 1) Then tetracorr = "<> Routine did not converge. Try increasing the maximum number of iterations (NITER) in Subroutine TETRA (see Brown, 1977)." Exit Function End If If (IFAULT = 2) Then tetracorr = "<> A negative frequency entered, or a row or column marginal is 0 (see Brown, 1977)." End If End Function Sub tetra(A As Double, B As Double, C As Double, D As Double, r As Double, SDR As Double, SDZERO As Double, ITYPE As Integer, IFAULT As Integer) 'C---------------------------------------------------------------------- 'C Applied Statistics algorithm AS 116 'C 'C Reference: Brown MB. Algorithm AS 116: the tetrachoric 'C correlation and its standard error. Applied Statistics, 1977, 26, 'C 343-351. 'C 'C All lowercase code in this subroutine added by me (but I have 'C commented these lines out.) 'C 'C John Uebersax 'C April 16, 2000 ' ' Fortran to VBA ' Keizo Hori ' Dec 4, 2000 'C---------------------------------------------------------------------- ' implicit double precision (a-h, o-z) 'C 'C ALGORITHM AS 116 APPL. STATIST. (1977) VOL.26, NO.3 'C 'C TO COMPUTE THE TETRACHORIC CORRELATION (R) AND ITS STANDARD 'C ERRORS (SDR AND SDZERO) FROM THE FREQUENCIES OF A 2*2 TABLE 'C (A, B, C AND D) 'C X AND W ARE CONSTANTS USED IN GAUSSIAN QUADRATURE Dim zero As Double, half As Double Dim one As Double, two As Double, four As Double, six As Double Dim twopi As Double, sqt2pi As Double, rlimit As Double, rcut As Double, Uplim As Double Dim cons0 As Double, CHALF As Double, CONV As Double, CITER As Double, NITER As Integer 'C Dim kdelta As Integer, aa As Double, bb As Double, cc As Double, dd As Double Dim delta As Double, tot As Double, ee As Double Dim probac As Double, probab As Double, PROBAA As Double, ie As Integer Dim tempa As Double, tempb As Double, KSIGN As Integer, ITER As Integer Dim zac As Double, zab As Double, ss As Double, rr As Double Dim wa As Double, wb As Double, va As Double, vb As Double, waa As Double, vaa As Double Dim term As Double, iterm As Integer, RRPREV As Double Dim Sum As Double, DERIV As Double, sr As Double, dr As Double, cof As Double Dim SUMPRV As Double, Prob As Double, RRSQ As Double, AMID As Double Dim XLEN As Double, IQUAD As Integer, xla As Double, xlb As Double, RREST As Double Dim PDF As Double, PAC As Double, pab As Double zero = 0# one = 1# two = 2# four = 4# six = 6# half = 0.5 twopi = 6.28318531 sqt2pi = 2.50662827 rlimit = 0.9999 rcut = 0.95 Uplim = 5# cons0 = 1E-36 CHALF = 1E-18 CONV = 0.00000001 CITER = 0.000001 NITER = 25 'C; 'C data conv /1E-8/, citer /1E-6/, niter /999/; 'C; Dim x(1 To 16) As Double, w(1 To 16) As Double '// DIMENSION X(16), W(16); x(1) = 0.9972638618 x(2) = 0.9856115115 x(3) = 0.9647622556 x(4) = 0.9349060759 x(5) = 0.8963211558 x(6) = 0.8493676137 x(7) = 0.794483796 x(8) = 0.7321821187 x(9) = 0.6630442669 x(10) = 0.5877157572 x(11) = 0.5068999089 x(12) = 0.4213512761 x(13) = 0.3318686023 x(14) = 0.2392873623 x(15) = 0.1444719616 x(16) = 0.0483076657: w(1) = 0.00701861 w(2) = 0.0162743947 w(3) = 0.0253920653 w(4) = 0.0342738629 w(5) = 0.042835898 w(6) = 0.0509980593 w(7) = 0.0586840935 w(8) = 0.0658222228 w(9) = 0.0723457941 w(10) = 0.0781938958 w(11) = 0.0833119242 w(12) = 0.087652093 w(13) = 0.0911738787 w(14) = 0.0938443991 w(15) = 0.0956387201 w(16) = 0.0965400885 'C 'C data conv /1E-8/, citer /1E-6/, niter /999/ 'C 'C 'C INITIALIZATION 'C r = zero SDZERO = zero SDR = zero ITYPE = 0 IFAULT = 0 'C 'C CHECK IF ANY CELL FREQUENCY IS NEGATIVE 'C If (A < zero Or B < zero Or C < zero Or D < zero) Then GoTo 92 'C 'C CHECK IF ANY FREQUENCY IS ZERO AND SET KDELTA 'C kdelta = 1 delta = zero If (A = zero Or D = zero) Then kdelta = 2 If (B = zero Or C = zero) Then kdelta = kdelta + 2 'C 'C KDELTA=4 MEANS TABLE HAS ZERO ROW OR COLUMN, RUN IS TERMINATED 'C Select Case kdelta Case 1: GoTo 4 Case 2: GoTo 1 Case 3: GoTo 2 Case 4: GoTo 92 End Select 'C 'C DELTA IS 0.0, 0.5 OR -0.5 ACCORDING TO WHICH CELL IS ZERO 'C 1 delta = half If (A = zero And D = zero) Then r = -one GoTo 4 2 delta = -half If (B = zero And C = zero) Then r = one 4 If (r <> zero) Then ITYPE = 3 'C 'C STORE FREQUENCIES IN AA, BB, CC AND DD 'C aa = A + delta bb = B - delta cc = C - delta dd = D + delta tot = aa + bb + cc + dd 'C 'C CHECK IF CORRELATION IS NEGATIVE, ZERO, POSITIVE 'C ' IF (AA * DD - BB * CC) 7, 5, 6 ee = aa * dd - bb * cc Select Case ee Case Is = 0: GoTo 5 Case Is < 0: GoTo 7 Case Is > 0: GoTo 6 End Select 5 ITYPE = 4 'C 'C COMPUTE PROBABILITIES OF QUADRANT AND OF MARGINALS 'C PROBAA AND PROBAC CHOSEN SO THAT CORRELATION IS POSITIVE. 'C KSIGN INDICATES WHETHER QUADRANTS HAVE BEEN SWITCHED 'C 6 PROBAA = aa / tot probac = (aa + cc) / tot KSIGN = 1 GoTo 8 7 PROBAA = bb / tot probac = (bb + dd) / tot KSIGN = 2 8 probab = (aa + bb) / tot 'C 'C COMPUTE NORMAL DEVIATES FOR THE MARGINAL FREQUENCIES 'C SINCE NO MARGINAL CAN BE ZERO, IE IS NOT CHECKED 'C zac = PPND(probac, ie) zab = PPND(probab, ie) ss = Exp(-half * (zac ^ 2 + zab ^ 2)) / twopi 'C 'C WHEN R IS 0.0, 1.0 OR -1.0, TRANSFER TO COMPUTE SDZERO 'C If (r <> zero Or ITYPE > 0) Then GoTo 85 'C 'C WHEN MARGINALS ARE EQUAL, COSINE EVALUATION IS USED 'C If (A = D And B = C) Then GoTo 60 'C 'C INITIAL ESTIMATE OF CORRELATION IS YULES Y 'C rr = (Sqr(aa * dd) - Sqr(bb * cc)) ^ 2 / Abs(aa * dd - bb * cc) ITER = 0 'C 'C IF RR EXCEEDS RCUT, GAUSSIAN QUADRATURE IS USED 'C 10 If (rr > rcut) Then GoTo 40 'C 'C TETRACHORIC SERIES IS COMPUTED 'C 'C INITIALIZATION 'C va = one vb = zac wa = one wb = zab term = one iterm = 0 Sum = probab * probac DERIV = zero sr = ss 15 If (Abs(sr) > cons0) Then GoTo 20 'C 'C RESCALE TERMS TO AVOID OVERFLOWS AND UNDERFLOWS 'C sr = sr / cons0 va = va * CHALF vb = vb * CHALF wa = wa * CHALF wb = wb * CHALF 'C 'C FORM SUM AND DERIVATIVE OF SERIES 'C 20 dr = sr * va * wa sr = sr * rr / term cof = sr * va * wa 'C 'C ITERM COUNTS NO. OF CONSECUTIVE TERMS < CONV 'C iterm = iterm + 1 If (Abs(cof) > CONV) Then iterm = 0 Sum = Sum + cof DERIV = DERIV + dr vaa = va waa = wa va = vb wa = wb vb = zac * va - term * vaa wb = zab * wa - term * waa term = term + one If (iterm < 2 Or term < six) Then GoTo 15 'C 'C CHECK IF ITERATION CONVERGED 'C If (Abs(Sum - PROBAA) > CITER) Then GoTo 25 'C 'C ITERATION HAS CONVERGED, SET ITYPE 'C ITYPE = term GoTo 70 'C 'C CALCULATE NEXT ESTIMATE OF CORRELATION 'C 25 ITER = ITER + 1 'C 'C write (*,*) iter, rr 'C 'C 'C IF TOO MANY ITERATlONS, RUN IS TERMINATED 'C If (ITER >= NITER) Then GoTo 93 delta = (Sum - PROBAA) / DERIV RRPREV = rr rr = rr - delta If (ITER = 1) Then rr = rr + half * delta If (rr > rlimit) Then rr = rlimit If (rr < zero) Then rr = zero GoTo 10 'C 'C GAUSSIAN QUADRATURE 'C 40 If (ITER > 0) Then GoTo 41 'C 'C INITIALIZATION, IF THIS IS FIRST ITERATION 'C Sum = probab * probac RRPREV = zero 'C 'C INITIALIZATION 'C 41 SUMPRV = probab - Sum Prob = bb / tot If (KSIGN = 2) Then Prob = aa / tot ITYPE = 1 'C 'C LOOP TO FIND ESTIMATE OF CORRELATION 'C COMPUTATION OF INTEGRAL (SUM) BY QUADRATURE 'C 42 RRSQ = Sqr(one - rr ^ 2) AMID = half * (Uplim + zac) XLEN = Uplim - AMID Sum = zero For IQUAD = 1 To 16 xla = AMID + x(IQUAD) * XLEN xlb = AMID - x(IQUAD) * XLEN 'C 'C TO AVOID UNDERFLOWS, TEMPA AND TEMPB ARE USED 'C tempa = (zab - rr * xla) / RRSQ If (tempa >= -six) Then Sum = Sum + w(IQUAD) * Exp(-half * xla ^ 2) * alnorm(tempa, False) tempb = (zab - rr * xlb) / RRSQ If (tempb >= -six) Then Sum = Sum + w(IQUAD) * Exp(-half * xlb ^ 2) * alnorm(tempb, False) Next IQUAD '44 CONTINUE Sum = Sum * XLEN / sqt2pi 'C 'C CHECK IF ITERATION HAS CONVERGED 'C If (Abs(Prob - Sum) <= CITER) Then GoTo 70 ITER = ITER + 1 'C 'C write (*,*) iter, rr 'C 'C 'C IF TOO MANY ITERATIONS, RUN IS TERMINATED 'C If (ITER >= NITER) Then GoTo 93 'C 'C ESTIMATE CORRELATION FOR NEXT ITERATION BY LINEAR INTERPOLATION 'C RREST = ((Prob - Sum) * RRPREV - (Prob - SUMPRV) * rr) / (SUMPRV - Sum) 'C 'C IS ESTIMATE POSITIVE AND LESS THAN UPPER LIMIT 'C If (RREST > rlimit) Then RREST = rlimit If (RREST < zero) Then RREST = zero RRPREV = rr rr = RREST SUMPRV = Sum 'C 'C IF ESTIMATE HAS SAME VALUE ON TWO ITERATIONS, STOP ITERATION 'C If (rr = RRPREV) Then GoTo 70 GoTo 42 'C 'C WHEN ALL MARGINALS ARE EQUAL THE COSINE FUNCTION IS USED 'C 60 rr = -Cos(twopi * PROBAA) ITYPE = 2 'C 'C COMPUTE SDR 'C 70 r = rr RRSQ = Sqr(one - r ^ 2) If (kdelta > 1) Then ITYPE = -ITYPE If (KSIGN = 1) Then GoTo 71 r = -r zac = -zac 71 PDF = Exp(-half * (zac ^ 2 - two * r * zac * zab + zab ^ 2) / RRSQ ^ 2) / (twopi * RRSQ) PAC = alnorm((zac - r * zab) / RRSQ, False) - half pab = alnorm((zab - r * zac) / RRSQ, False) - half SDR = (aa + dd) * (bb + cc) / four + pab ^ 2 * (aa + cc) * (bb + dd) + PAC ^ 2 * (aa + bb) * (cc + dd) + two * pab * PAC * (aa * dd - bb * cc) - pab * (aa * bb - cc * dd) - PAC * (aa * cc - bb * dd) If (SDR < zero) Then SDR = zero SDR = Sqr(SDR) / (tot * PDF * Sqr(tot)) 'C 'C COMPUTE SDZERO 'C 85 SDZERO = Sqr((aa + bb) * (aa + cc) * (bb + dd) * (cc + dd) / tot) / (tot ^ 2 * ss) If (r = zero) Then SDR = SDZERO GoTo 99 'C 'C ERROR TERMINATIONS 'C 92 IFAULT = 1 93 IFAULT = IFAULT + 1 'C 99 Exit Sub End Sub Function PPND(P As Double, IER As Integer) As Double '=NORMSINV(p)と同じ 'C 'C ALGORITHM AS 111, APPL.STATIST., VOL.26, 118-121, 1977. 'C 'C PRODUCES NORMAL DEVIATE CORRESPONDING TO LOWER TAIL AREA = P. 'C 'C See also AS 241 which contains alternative routines accurate to 'C about 7 and 16 decimal digits. 'C ' IMPLICIT DOUBLE PRECISION (A-H,O-Z) Dim zero As Double, one As Double, half As Double Dim a0 As Double, a1 As Double, a2 As Double, a3 As Double Dim b1 As Double, b2 As Double, b3 As Double, b4 As Double Dim c1 As Double, c2 As Double, c3 As Double, c0 As Double Dim d1 As Double, d2 As Double, split As Double Dim up As Boolean Dim q As Double, r As Double split = 0.42 a0 = 2.50662823884 a1 = -18.61500062529 a2 = 41.39119773534 a3 = -25.44106049637 b1 = -8.4735109309 b2 = 23.08336743743 b3 = -21.06224101826 b4 = 3.13082909833 c0 = -2.78718931138 c1 = -2.29796479134 c2 = 4.85014127135 c3 = 2.32121276858 d1 = 3.54388924762 d2 = 1.63706781897 zero = 0# one = 1# half = 0.5 'C IER = 0 q = P - half If (Abs(q) > split) Then GoTo 10 'C 'C 0.08 < P < 0.92 'C r = q * q PPND = q * (((a3 * r + a2) * r + a1) * r + a0) / ((((b4 * r + b3) * r + b2) * r + b1) * r + one) Exit Function 'C 'C P < 0.08 OR P > 0.92, SET R = MIN(P,1-P) 'C 10 r = P If (q > zero) Then r = one - P If (r <= zero) Then GoTo 20 r = Sqr(-Log(r)) PPND = (((c3 * r + c2) * r + c1) * r + c0) / ((d2 * r + d1) * r + one) If (q < zero) Then PPND = -PPND Exit Function 20 IER = 1 PPND = zero End Function Function alnorm(x As Double, upper As Boolean) As Double 'false の場合, =NORMSDIST(x) と同じ 'C 'C Algorithm AS66 Applied Statistics (1973) vol22 no.3 'C 'C Evaluates the tail area of the standardised normal curve 'C from x to infinity if upper is .true. or 'C from minus infinity to x if upper is .false. 'C Dim zero As Double, one As Double, half As Double Dim con As Double, z As Double, y As Double Dim P As Double, q As Double, r As Double Dim a1 As Double, a2 As Double, a3 As Double Dim b1 As Double, b2 As Double Dim c1 As Double, c2 As Double, c3 As Double, c4 As Double, c5 As Double, c6 As Double Dim d1 As Double, d2 As Double, d3 As Double, d4 As Double, d5 As Double Dim up As Boolean 'C*** machine dependent constants Dim ltone As Double, utzero As Double zero = 0# one = 1# half = 0.5 ltone = 7# utzero = 18.66 con = 1.28 P = 0.398942280444 q = 0.39990348504 r = 0.398942280385 a1 = 5.75885480458 a2 = 2.62433121679 a3 = 5.92885724438 b1 = -29.8213557807 b2 = 48.6959930692 c1 = -0.000000038052 c2 = 0.000398064794 c3 = -0.151679116635 c4 = 4.8385912808 c5 = 0.742380924027 c6 = 3.99019417011 d1 = 1.00000615302 d2 = 1.98615381364 d3 = 5.29330324926 d4 = -15.1508972451 d5 = 30.789933034 'C up = upper z = x If (z >= zero) Then GoTo 10 up = Not up z = -z 10 If (z <= ltone Or up And z <= utzero) Then GoTo 20 alnorm = zero GoTo 40 20 y = half * z * z If (z > con) Then GoTo 30 'C alnorm = (P - q * y / (y + a1 + b1 / (y + a2 + b2 / (y + a3)))) alnorm = z * alnorm alnorm = half - alnorm Debug.Print alnorm GoTo 40 30 alnorm = (z + c3 + d3 / (z + c4 + d4 / (z + c5 + d5 / (z + c6)))) alnorm = r * Exp(-y) / (z + c1 + d1 / (z + c2 + d2 / alnorm)) 40 If (Not up) Then alnorm = one - alnorm End Function