CTRHO
Estimates the bivariate normal correlation coefficient using a contingency table.
Required Arguments
TABLE — NROW by NCOL contingency table containing the observed counts. (Input)
RHO — Maximum likelihood estimate of the correlation coefficient. (Output)
VAR — Estimated asymptotic variance of RHO. (Output)
PLTMY — Vector of length NROW + NCOL ‑ 2 containing the points of polytomy of the marginal rows and columns of TABLE. (Output)
The first NROW ‑ 1 elements of PLTMY are the points of polytomy for the rows while the last NCOL ‑ elements are the points of polytomy for the columns.
PROB — NROW by NCOL matrix containing the bivariate normal probabilities corresponding to RHO and PLTMY. (Output)
DRIV — NROW by NCOL matrix containing the partial derivatives of the bivariate normal probability with respect to RHO. (Output)
Optional Arguments
NROW — Number of rows in the table. (Input)
Default: NROW = size (TABLE,1).
NCOL — Number of columns in the table. (Input)
Default: NCOL = size (TABLE,2).
LDTABL — Leading dimension of TABLE exactly as specified in the dimension statement of the calling program. (Input)
Default: LDTABL = size (TABLE,1).
EPS — Convergence criterion in the iterative estimation. (Input)
RHO will be within EPS of the maximum likelihood estimate unless roundoff errors prevent this precision. EPS must be less than 2. EPS less than or equal to zero defaults to 0.00001.
Default: EPS = .00001.
LDPROB — Leading dimension of PROB exactly as specified in the dimension statement in the calling program. (Input)
Default: LDPROB = size (PROB,1).
LDDRIV — Leading dimension of DRIV exactly as specified in the dimension statement in the calling program. (Input)
Default: LDDRIV = size (DRIV,1).
FORTRAN 90 Interface
Generic: CALL CTRHO (TABLE, RHO, VAR, PLTMY, PROB, DRIV [, …])
Specific: The specific interface names are S_CTRHO and D_CTRHO.
FORTRAN 77 Interface
Single: CALL CTRHO (NROW, NCOL, TABLE, LDTABL, EPS, RHO, VAR, PLTMY, PROB, LDPROB, DRIV, LDDRIV)
Double: The double precision name is DCTRHO.
Description
Routine CTRHO computes the maximum likelihood estimate and the asymptotic variance for the correlation coefficient of a bivariate normal population from a two-way contingency table. The maximum likelihood estimates are conditional upon the points of polytomy in the marginal distribution. The resulting estimate for the correlation coefficient should be very close to the unconditional estimate (see Martinson and Hamdan 1972).
The points of polytomy for the row and column marginal probabilities are first computed. If the i-th cumulative column marginal is denoted by pci, then the point of polytomy xi is given as
where Φ denotes the cumulative normal distribution. Let αi, i = 0, …, r denote these points for the row marginal cumulative probabilities where r = NROW, α0 = ‑∞, and αr = ∞. Similarly, let βj, j = 0, …, c denote the points of polytomy for the columns where c = NCOL. Then, the probability of the (i, j) cell in the table, pij, is defined as
where X and Y are the bivariate random variables. Maximum likelihood estimates for the correlation coefficient are computed based upon the bivariate normal density. The likelihood is specified by the multinomial distribution of the table using probabilities pij.
Routine CTRHO assumes that the row random variable decreases with increasing row number while the column variable increases with the column number. If this is not the case, the sign of the estimated correlation coefficient may need to be changed.
Example
The data are taken from Martinson and Hamdan (1972), who attribute it to Karl Pearson. The row variable is head breadth (in millimeters) for a human male while the column variable is the head breadth of his sister. Head breadth increases across the columns and decreases down the row. The row and column variables have been categorized into one of three intervals. The original table is as follows:
1.0 | 36.5 | 77.5 |
52.5 | 340.5 | 143.5 |
40.5 | 58.0 | 9.0 |
Note that routine CTRHO can accept other than integer counts. It is not clear from Martinson and Hamdan (1972) how the non-integral counts arise in the table here. The correlation is estimated to be 0.5502.
USE CTRHO_INT
USE UMACH_INT
USE WRRRN_INT
IMPLICIT NONE
INTEGER LDDRIV, LDPROB, LDTABL, NCOL, NROW
PARAMETER (LDDRIV=3, LDPROB=3, LDTABL=3, NCOL=3, NROW=3)
!
INTEGER NOUT
REAL DRIV(LDDRIV, NCOL), PLTMY(NROW+NCOL-2),PROB(LDPROB,NCOL), &
RHO, TABLE(LDTABL, NCOL), VAR
!
DATA TABLE/1.0, 52.5, 40.5, 36.5, 340.5, 58.0, 77.5, 143.5, 9.0/
!
!
CALL CTRHO (TABLE, RHO, VAR, PLTMY, PROB, DRIV)
!
CALL UMACH (2, NOUT)
WRITE (NOUT,*) 'RHO =', RHO, ' VAR =', VAR
CALL WRRRN ('PLTMY', PLTMY, 1, NROW+NCOL-2, 1, 0)
CALL WRRRN ('PROB', PROB)
CALL WRRRN ('DRIV', DRIV)
END
Output
RHO = 0.549125 VAR = 1.33199E-03
PLTMY
1 2 3 4
-1.073 1.030 -1.156 0.516
PROB
1 2 3
1 0.0015 0.0517 0.0983
2 0.0700 0.4398 0.1970
3 0.0523 0.0816 0.0077
DRIV
1 2 3
1 -0.0134 -0.0984 0.1118
2 -0.0717 0.1388 -0.0672
3 0.0851 -0.0404 -0.0447
