Introduction to the Tetrachoric and Polychoric Correlation Coefficients

• Introduction
• Summary
• Pros and Cons
• Intuitive Explanation
• Detailed Description
• Using the Polychoric Correlation to Measure Agreement
• Extensions and generalizations
  • Skewed distributions
  • Nonparametric distributions
• Examples
  • Tetrachoric correlation
  • Polychoric correlation
  • Polychoric correlation with relaxed assumptions
• Factor analysis and SEM
• Software Updated!
  • Tetrachoric correlation
  • Polychoric correlation
  • Tetrachoric/polychoric correlation coefficient in SAS
  • Tetrachoric/polychoric correlation coefficient in R
  • Generalized latent correlation
• Generalized Latent Correlation
• MCMC Estimation
• References

Go to Agreement Statistics main page
Go to Latent Structure Analysis main page

Introduction

A clear, concise description of the tetrachoric and polychoric correlation coefficients, including issues relating to their estimation, is found in Drasgow (1988). Olsson (1979) is also helpful.

What distinguishes the present discussion is the view that the tetrachoric and polychoric correlation models are special cases of latent trait modeling. (This is not a new observation, but it is sometimes overlooked). Recognizing this opens up important new possibilities. In particular, it allows one to relax the distributional assumptions which are the most limiting feature of the "classical" tetrachoric and polychoric correlation models.

Note in any case that the terms tetrachoric correlation and polychoric correlation are obsolete and arguably inaccurate. They refer to the tetrachoric series and polychoric series, numerical methods previously (before modern computers) used to facilitate calculations. Now these correlations are estimated by maximum likelihood or other means. Hence a term like latent correlation (or latent continuous correlation) is more appropriate.

Summary

These statistics make certain assumptions, however. With the polychoric correlation, the assumptions can be tested. The assumptions cannot be tested with the tetrachoric correlation if there are only two raters; in some applications, though, theoretical considerations may justify the use of the tetrachoric correlation without a test of model fit.

Pros and Cons: Tetrachoric and Polychoric Correlation Coefficients

Pros:

• These statistics express rater association in a familiar form--a correlation coefficient.
• They provide a way to separately quantify association and similarity of category definitions.
• They do not depend on number of rating levels; results can be compared for studies where the number of rating levels is different.
• They can be used even if different raters have different numbers of rating levels.
• The assumptions can be easily tested for the polychoric correlation.
• Estimation software is routinely available (e.g., SAS PROC FREQ, and PRELIS).

Cons:

• Model assumptions not always appropriate--for example, if the latent trait is truly discrete.
• For only two raters, there is no way to test the assumptions of the tetrachoric correlation.

Intuitive Explanation

Consider the example of two psychiatrists (Raters 1 and 2) making a diagnosis for presence/absence of Major Depression. Though the diagnosis is dichotomous, we allow that depression as a trait is continuously distributed in the population.

+---------------------------------------------------------------+
|                                                               |
|                                                               |
|     |                         *                               |
|     |                     *       *                           |
|     |                  *             *                        |
|     |                *               | *                      |
|     |               *                |  *                     |
|     |             **                 |   **                   |
|     |          ***                   |     ***                |
|     |       ***                      |        ***             |
|     |  *****                         |           *****        |
|     +--------------------------------+----------------> Y     |
|                 not depressed        t    depressed           |
|                                                               |
+---------------------------------------------------------------+

In diagnosing a given case, a rater considers the case's level of depression, Y, relative to some threshold, t: if the judged level is above the threshold, a positive diagnosis is made; otherwise the diagnosis is negative.

Figure 2 portrays the situation for two raters. It shows the distribution of cases in terms of depression level as judged by Rater 1 and Rater 2.

Figure 2. Joint distribution (ellipse) of depression
severity as judged by two raters (Y1 and Y2); and
discretizing thresholds (t1 an t2)

a, b, c and d denote the proportion of cases that fall in each region defined by the two raters' thresholds. For example, a is the proportion below both raters' thresholds and therefore diagnosed negative by both.

These proportions correspond to a summary of data as a 2 x 2 cross-classification of the raters' ratings.

+------------------------------------------------+
|                                                |
|                      Rater 1                   |
|                     -       +                  |
|                 +-------+-------+              |
|               - |   a   |   b   | a + b        |
|      Rater 2    +-------+-------+              |
|               + |   c   |   d   | c + d        |
|                 +-------+-------+              |
|                   a + c   b + d     1          |
|                                                |
+------------------------------------------------+

Again, a, b, c and d in Figure 3 represent proportions (not frequencies).

Once we know the observed cross-classification proportions a, b, c and d for a study, it is a simple matter to estimate the model represented by Figure 2. Specifically, we estimate the location of the discretizing thresholds, t1 and t2, and a third parameter, rho, which determines the "fatness" of the ellipse. Rho is the tetrachoric correlation, or r*. It can be interpreted here as the correlation between judged disease severity (before application of thresholds) as viewed by Rater 1 and Rater 2.

The principle of estimation is simple: basically, a computer program tries various combinations for t1, t2 and r* until values are found for which the expected proportions for a, b, c and d in Figure 2 are as close as possible to the observed proportions in Figure 3. The parameter values that do so are regarded as (estimates of) the true, population values.

The polychoric correlation, used when there are more than two ordered rating levels is a straightforward extension of the model above. The difference is that there are more thresholds, more regions in Figure 2, and more cells in Figure 3. But again the idea is to find the values for thresholds and r* that maximize similarity between model-expected and observed cross-classification proportions.

Detailed Description

Introduction

In many situations, even though a trait may be continuous, it may be convenient to divide it into ordered levels. For example, for research purposes, one may classify levels of headache pain into the categories none, mild, moderate and severe. Even for trait usually viewed as discrete, one might still consider continuous gradations--for example, people infected with the flu virus exhibit varying levels of symptom intensity.

The tetrachoric correlation and polychoric correlation coefficients are appropriate when the latent trait that forms the basis of ratings can be viewed as continuous. We will outline here the measurement model and assumptions for the tetrachoric correlation. The model and assumptions for the polychoric correlation are the same--the only difference is that there are more threshold parameters for the polychoric correlations, corresponding to the greater number ordered rating levels.

Measurement Model

We begin with some notation and definitions. Let:

X1 and X2 be the manifest (observed) ratings by Raters (or procedures, diagnostic tests, etc.) 1 and 2; these are discrete-valued variables;

Y1, Y2 be latent continuous variables associated with X1 and X2; these are the pre-discretized, continuous "impressions" of the trait level, as judged by Raters 1 and 2;

T be the true, latent trait level of a case.

A rating or diagnosis of a case begins with the case's true trait level, T. This information, along with "noise" (random error) and perhaps other information unrelated to the true trait which a given rater may consider (unique variation), leads to each rater's impression of the case's trait level (Y1 and Y2). Each rater applies discretizing thresholds to this judged trait level to yield a dichotomous or ordered-category rating (X1 and X2).

Stated more formally, we have:

Y1 = bT + u1 + e1,
Y2 = bT + u2 + e2,

where b is a regression coefficient, u1 and u2 are the unique components of the raters' impressions, and e1 and e2 represent random error or noise. It turns out that unique variation and error variation behave more or less the same in the model, and the former can be subsumed under the latter. Thus we may consider the simpler model:

Y1 = b1T + e1,
Y2 = b2T + e2.

The tetrachoric correlation assumes that the latent trait T is normally distributed. As scaling is arbitrary, we specify that T ~ N(0, 1). Error is similarly assumed to be normally distributed (and independent both between raters and across cases). For reasons we need not pursue here, the model loses no generality by assuming that var(e1) = var(e2). We therefore stipulate that e1, e2 ~ N(0, sigma_e). A consequence of these assumptions is that Y1 and Y2 must also be normally distributed. To fix the scale, we specify that var(Y1) = var(Y2) = 1. It follows that b1 = b2 = b = the correlation of both Y1 and Y2 with the latent trait.

We define the tetrachoric correlation, r*, as

r* = b²

A simple "path diagram" may clarify this:

+-------------------------------------+
|                                     |
|                                     |
|               b    b                |
|          Y1 <--- T ---> Y2          |
|                                     |
|                                     |
+-------------------------------------+

One may wish to test the null hypothesis of no correlation between raters. There are at least two ways to do this. The first makes use of the estimated standard error of the polychoric correlation under the null hypothesis of r* = 0. At least for the tetrachoric correlation, there is a simple closed-form expression for this standard error (Brown, 1977). Knowing this value, one may calculate a z value as:

                 r*
        z = -----------
             sigmar*(0)

          --------------------------
                      Rater 2
                     ---------
          Rater 1     -     +   Total
          ---------------------------
            -        40    10     50

            +        20    30     50
          --------------------------
          Total      60    40    100
          --------------------------
          Table 1 (draft)

        rho       0.6071  (0.1152)

The thresholds (std. errors) for the two raters are estimated as:

        Rater 1   0.0000  (0.1253)
        Rater 2   0.2533  (0.1268)

Tallis suggested that the number of lambs born is a manifestation of the ewe's fertility--a continuous and potentially normally distributed variable. Clearly the situation is more complex than the simple "continuous normal variable plus discretizing thresholds" assumptions allow for. We consider the data simply for the sake of a computational example.

          -----------------------------------
           Lambs     Lambs born in 1952
          born in    ------------------
            1953      None    1    2    Total
          -----------------------------------
            None       58    52    1     111

             1         26    58    3      87

             2          8    12    9      29
          -----------------------------------
           Total       92   122   13     227
          -----------------------------------
          Table 2 (draft)

The data in Table 2 are analyzed with the POLYCORR program (Uebersax, 2000). Application of the joint ML approach produces the following estimates (standard errors):

          rho                  0.4192  (0.0761)
          threshold 2, 1952   -0.2421  (0.0836)
          threshold 3, 1952    1.5938  (0.1372)
          threshold 2, 1953   -0.0297  (0.0830)
          threshold 3, 1953    1.1331  (0.1063)

          rho                  0.4199  (0.0747)
          threshold 2, 1952   -0.2397
          threshold 3, 1952    1.5781
          threshold 2, 1953   -0.0276
          threshold 3, 1953    1.1371

          ---------------------------------------------
           Rating         Rating of Rater 1
          of Rater  ---------------------------
             2       1    2    3    4    5    6   Total
          ---------------------------------------------
             1      30    1    0    0    0    0     31
             2       0   10    2    0    0    0     12
             3       0    4    8    3    1    0     16
             4       0    3    3   37    9    0     52
             5       0    0    1   25   71   49    146
             6       0    0    0    2   20  181    203
          ---------------------------------------------
           Total    30   18   14   67  101  230    460
          ---------------------------------------------
          Table 3 (draft).
          Ratings of plant health by two judges

The distribution of the latent trait estimated by the model is follows:

         .5 +                                       *
    D       |                                       *
    e    .4 +                                       *
    n       |                                       *
    s    .3 +                                       *
    i       |                                *      *
    t    .2 +                                *      *
    y       |                                *      *
         .1 +                         *      *      *
            |    *             *      *      *      *
            +----*------*------*------*------*------*----
               -2.5   -1.5   -0.5    0.5    1.5    2.5

                          Latent Trait Level
    Figure 5 (draft).
    Estimated latent trait distribution

POLYCORR is a program I've written to estimate the polychoric correlation and its standard error using either joint ML or two-step estimation. Goodness-of-fit and a lot of other information are also provided.  Note: this program is just for a single pair of variables, or a few considered two at a time. It does not estimate a matrix of tetra- or polychoric correlations.

• Basic version. This handles square tables only (i.e., models where both items have the same number of levels).

• Advanced version. This allows non-square tables and has other advanced technical features, such as the ability to combine cells during estimation.

Note. At present the basic version doesn't run on 64-bit Windows operating systems. This will be addressed later in September 2015.

R
(See next section.)

SAS
A single polychoric or tetrachoric correlation coefficient can be calculated with the PLCORR option of SAS PROC FREQ. Example:

proc freq;
   tables var1*var2 / plcorr maxiter=100;
run;

Joint estimation is used. The standard error is supplied, but not thresholds. No goodness-of-fit test is performed.

As of SAS® Base version 9.3 and SAS/STAT version 13.1 there are two new ways to compute polychoric/tetrachoric correlation coefficients: with a new POLYCHORIC option in PROC CORR, and with PROC IRT. Given more than two variables, both procedures can supply a matrix of tetrachoric/polychoric correlation coefficients. PROC CORR will (if requested) write results listwise to an output dataset that an be easily manipualted to produce a matrix. PROC IRT will supply matrix output directly, but, unlike PROC CORR, does not supply information on model fit for individual correlations.

The following code uses PROC CORR to estimate the polychoric correlation coefficient for the lambs and ewes data in Table 2 above.

data one;
  input x1 x2 f;
  datalines;
0 0 58
0 1 52
0 2 1
1 0 26
1 1 58
1 2 3
2 0 8
2 1 12
2 2 9
;

data two; set one;
  do i = 1 to f; output; end;
  drop i f;
run;

proc corr data = two polychoric;
  var x1 x2;
* optional: save results in dataset named pc;
  ods output PolychoricCorr = pc;
run;

PROC CORR supplies tetrachoric/polychoric correlation coefficients, standard errors, and two tests of statistical significance (i.e., that the correlation is significantly different from 0). To test goodness of fit (for a single polychoric correlation coefficient only) one can use PROC IRT:

proc irt data = two link = probit polychoric;
  model x1-x2 / resfunc = graded;
  equality x1-x2 / parm = [slope];
run;

Likelihood ratio and Pearson chi-squared statistics and df are reported in the Model Fit Statistics section of the output. Again, this only works for a single polychoric correlation coefficient (model fit for tetrachoric correlation coefficients is untestable). One can use PROC IRT to estimate a matrix of tetrachoric/polychoric correlation coefficients, but in this case the model fit statistics have a somewhat different interpretation.

For those with older SAS version, a SAS macro, %POLYCHOR, can construct a matrix of polychoric or tetrachoric correlations. The macro is relatively slow.

For tetrachoric correlations, if there is a single 0 frequency in the 2×2 crossclassification table for a pair of variables (see Figure 3 above), the PLCORR option of PROC FREQ and %POLYCHOR might unnecessarily supply a missing value result, at least if maxiter is left at the default value of 20. So far I have found this problem is avoided by setting maxiter higher, e.g., to 40, 50 or 100.

SPSS
SPSS has no intrinsic procedure to estimate polychoric correlations. As noted above, Dirk Enzmann has written an SPSS macro to estimate a matrix of tetrachoric correlations.

PRELIS. A useful program for estimating a matrix of polychoric or tetrachoric correlations is PRELIS. It includes a goodness-of-fit test for each pair of variables. Standard errors can be requested. PRELIS uses two-step estimation. Because it is supplied with LISREL, PRELIS is widely available. Most university computation centers probably already have copies and/or site licenses.

Mplus can estimate a matrix of polychoric and tetrachoric correlations and estimate their standard errors. Two-step estimation is used. Features similar to PRELIS/LISREL.

Stata
Stata's internal function for tetrachoric correlations is a very rough approximation (e.g., actual tetrachoric correlation = .5172, Stata reports .6169!) based on Edwards and Edwards (1984) and is unsuitable for many or most applications. A more accurate external module has been written by Stas Kolenikov to estimate a matrix of polychoric or tetrachoric correlations and their standard errors.

MicroFACT will estimate polychoric and tetrachoric correlations and standard errors. Provisions for smoothing an improper correlation matrix are supplied. No goodness-of-fit tests. A free student version that handles up to thirty variables can be downloaded. Also does factor analysis.

WinBUGS
As Albert (1992) and others have shown, it is possible to estimate the polychoric correlation using Bayesian (MCMC) methods. The WinBUGS program can probably be used for this. (If anyone has code for this, I'll place it here.)

Calculating the Tetrachoric/Polychoric Correlation Coefficient in R

John Fox has contributed an R library, polycor, with functions to estimate the polychoric (and tetrachoric) correlation.

How to Calculate the Polychoric Correlation Coefficient using R

1. Download R from the CRAN website (free) and install on your computer: http://cran.r-project.org/

2. Download and install RStudio (free): https//www.rstudio.com

3. Open RStudio.

4. At command prompt (">") type: install.packages("polycor")

5. When this action completes, type: library("polycor")

6. Supply frequencies for a two-way crossclassification table as a vector, e.g.,

x <- c(10, 20, 30, 40, 50, 60, 70, 80, 90)

7. Arrange column as a matrix, e.g.,

y <- matrix(x, 3, 3)

8. Estimate polychoric correlation using the 'polychor' (note spelling) function.

polychor(y, ML=T, std.err=T)

Results:

Polychoric Correlation, ML est. = -0.1183 (0.06098)
Test of bivariate normality: Chisquare = 1.216, df = 3, p = 0.7491

Row Thresholds
Threshold Std.Err.
1 -0.6227 0.06345
2 0.2535 0.05976

Column Thresholds
Threshold Std.Err.
1 -1.11100 0.07450
2 -0.08275 0.05914

The estimated value of the polychoric correlation coefficient is -0.1183 with an estimated standard error of 0.06098.

If x has only two levels (rows and columns), the tetrachoric correlation is estimated.

It is also possible to supply raw data in the form of two vectors, x and y. For example, let x = c(1, 2, 1, 1, 2, 1, 1, 1, 2, 1), y = c(2, 1, 2, 2, 1, 1, 1, 2, 1, 2). Note that these are two variables measured on the same subjects/objects, so their lengths and orderings must match.

x <- c(1, 2, 1, 1, 2, 1, 1, 1, 2, 1, 2, 2, 1, 1, 2, 2)
y <- c(1, 1, 2, 1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 2, 1)
polychor(x, y, ML=T, std.err=T)

Results:

Polychoric Correlation, ML est. = 0.3672 (0.3574)

Row Threshold
Threshold Std.Err.
0.1573 0.3147

Column Threshold
Threshold Std.Err.
0.1573 0.3147

(Note however that n=10) is far too few observations from which to estimate the tetrachoric correlation coefficient; an n of at least 50 is recommended.)

See also functions for the polychoric and polyserial correlation coefficients in the psych R library of William Revelle at Northwestern.

Another R function for the polychoric correlation coefficient has been written by David Duffy.

Generalized Latent Correlation

An obvious potential limitation of the polychoric correlation coefficient is the assumption that the latent variables have a bivariate Gaussian distribution. Ekström (2011) reports that Pearson was little bothered by this, evidently believing the polychoric correlation could be used effectively even when distributional assumptions weren't strictly met.

In any case, several papers have examined extensions based on relaxation of this assumption. Non-Gaussian distributions, mixtures of distributions and skewed Gaussian distributions have all been considered (Ekström, 2011; Kottas, Müller & Quintana, 2005; Quiroga, 1992; Roscino & Pollice, 2006; Timofeeva & Khailenko, 2016).

The approach I've taken in several papers reparameterizes the problem as a one-dimensional latent trait model. If the latent trait has a Gaussian distribution, results produced are identical to those of the usual polychoric correlation coefficient. However other distributions of the latent trait can be assumed — such as skewed distributions, a mixture of Gaussian distributions (Uebersax & Grove, 1993) or a semi-parametric distribution (Uebersax, 1993). Because integration in only one dimension is required, this method of computationally attractive. Some computer programs to estimate these models are described below.

The glc program generalizes the tetrachoric correlation to estimate the latent correlation between binary variables assuming a skewed latent trait distribution.

The skewed distribution is modeled as a mixture of two Gaussian distributions, the parameters of which the user supplies; that is, one specifies in advance the shape of the latent trait distribution, based on prior beliefs/knowledge. This program is much simpler to use than those described below. Several sets of data (summarized as a series of 2×2 tables) can be analyzed in a single run.

The LTMA program can similarly be used to estimate a generalized polychoric correlation, based on a latent trait mixture model (Uebersax & Grove, 1993). This is basically a fancier version of the glc program: (1) it handles ordered-categorical as well as dichotomous variables, and (2) it will estimate the shape of the latent trait distribution from the data (again, modeling it as a mixture of two component Gaussians).

The LLCA program can be used to estimate a polychoric correlation with nonparametric distributional assumptions. The latent trait is represented as a sequence of latent classes on a single continuum (Uebersax 1993b). That is, the latent trait distribution is modeled as a "histogram," where the densities at each point are estimated, rather as a continuous parametric distribution.

MCMC Estimation

Following Albert (1992), there have been many papers applying Markov Chain Monte Carlo (MCMC) methods to estimate the polychoric correlation coefficient. That the MCMC approach has many benefits (such as that ability to incorporate informative prior information into estimation) cannot be denied, and research will doubtless continue in this promising area. (More references will be added in the future.)

However at the same time it's important to bear in mind distinct advantages of a traditional (i.e., direct maximum likelihood) approach. These include: (1) more exact results, (2) ability to test model fit using chi-squared statistics, (3) ability to compare nested models via likelihood-ratio and similar statistics, and (4) faster computation, which, among other things, facilitiates estimation of an entire matrix of correlations.

Helpful Links

• "Polychoric correlation" - Wikipedia
• Correlation page (Statnotes) - by G. David Garson

References