krippendorffsalpha: An R Package for Measuring Agreement Using Krippendorff's Alpha Coefficient

R package krippendorffsalpha provides tools for measuring agreement using Krippendorff's Alpha coefficient, a well-known nonparametric measure of agreement (also called inter-rater reliability and various other names). This article first develops Krippendorff's Alpha in a natural way, and situates Alpha among statistical procedures. Then the usage of package krippendorffsalpha is illustrated via analyses of two datasets, the latter of which was collected during an imaging study of hip cartilage. The package permits users to apply the Alpha methodology using built-in distance functions for the nominal, ordinal, interval, or ratio levels of measurement. User-defined distance functions are also supported. The fitting function can accommodate any number of units, any number of coders, and missingness. Bootstrap inference is supported, and the bootstrap computation can be carried out in parallel.


Introduction
Krippendorff's α (Hayes and Krippendorff, 2007) is a well-known nonparametric measure of agreement (i.e., consistency of scoring among two or more raters for the same units of analysis (Gwet, 2014)). In R (Ihaka and Gentleman, 1996), Krippendorff's α can be applied using function kripp.alpha of package irr (Gamer et al., 2012), function kripp.boot of package kripp.boot (Proutskova and Gruszczynski, 2020), function krippalpha of package icr (Gamer et al., 2012), and functions krippen.alpha.raw and krippen.alpha.dist of package irrCAC (Gwet, 2019). But these packages fail to provide a number of useful features. In this article we present package krippendorffsalpha, which improves upon the above mentioned packages in (at least) the following ways. Package krippendorffsalpha • offers commonly used built-in distance functions for the nominal, ordinal, interval, and ratio levels of measurement, and also supports user-defined distance functions; • conforms to the R idiom by providing S3 methods confint, influence, plot, and summary; • supports embarrassingly parallel bootstrap computation; and • supports verbose communication with the user, including the display of a progress bar during the production of the bootstrap sample. The remainder of this article is organized as follows. In Section 2 we locate Krippendorff's α among statistical procedures. In Section 2.1 we first develop a special case of Krippendorff's α (call it α sed ) in a well-known parametric setting, and then we present α in its most general (i.e., nonparametric) form. In Section 2.2 we show that α is a type of multiresponse permutation procedure. In Section 2.3 we generalize α sed in a fully parametric fashion, arriving at Sklar's ω. In Section 3 we describe our package's bootstrap inference for α and compare the performance of our procedure to that of two alternative approaches. In Section 4 we briefly discuss robustness and influence. In Section 5 we provide a thorough demonstration of krippendorffsalpha's usage before concluding in Section 6.

Situating Krippendorff's Alpha among statistical procedures
Since Krippendorff's α is defined in terms of discrepancies (Krippendorff, 2013), at first glance one might conclude, erroneously, that α is a measure of dis-agreement and so answers the wrong question. In Sections 2.1-2.3 we will show, by examining Krippendorff's α's place among statistical procedures, that α is, in fact, a sensible measure of agreement. Also, establishing a context for α may help practitioners make educated decisions regarding α's use.
The UML class diagram (Fowler et al., 2004) shown below in Figure 1 provides a conceptual roadmap for our development. Briefly, a special case of α (which we denote as Alpha(SED) or α sed ) arises naturally in the context of the one-way mixed-effects ANOVA model. Alpha(SED) can then be generalized in a nonparametric fashion to arrive at Krippendorff's α as it has been presented by Hayes and Krippendorff (see Gwet (2015) for a development of nonparametric α in terms of agreement rather than discrepancies); this nonparametric form of α is a (slightly modified) multiresponse permutation procedure. Alternatively, α sed can be generalized in a parametric fashion to arrive at Sklar's ω, a Gaussian copula-based methodology for measuring agreement.
2.1. Parametric genesis of Krippendorff 's Alpha coefficient. In this section we develop Krippendorff's α (Hayes and Krippendorff, 2007) in an intuitive and bottom-up fashion. Our starting point is a fully parametric model, namely, the classic one-way mixed-effects ANOVA model (Ravishanker and Dey, 2001). To ease exposition we will consider only a balanced version of the model. We have, for n u units and n c coders, scores where • µ (the population mean) is a fixed real number, • the τ i are independent N (0, σ 2 τ ) random variables, • the ε ij are independent N (0, σ 2 ε ) random variables, and • the τ i are independent of the ε ij . In this setup we have n c Gaussian codes Y i1 , . . . , Y inc for unit i ∈ {1, . . . , n u }. Conditional on τ i , said codes are N (µ+τ i , σ 2 ε ) random variables. Since the variables share the "unit effect" τ i , the variables are correlated. The correlation, which is usually called the intraclass correlation, is given by We use α to denote this quantity precisely because Krippendorff's α is the intraclass correlation for codes that conform to this model. That is, for the one-way mixedeffects ANOVA model, Krippendorff's α is the intraclass correlation. The reader may recall that the estimator of α for this model iŝ whereȲ i• andȲ•• denote the arithmetic means for the ith unit and for the entire sample, respectively. The form of this estimator is not surprising, of course, since it is well known that assuming Gaussianity leads to variance estimators involving sums of weighted squared deviations from sample arithmetic means.
We can eliminate the arithmetic means in (1) by employing the identity This givesα (2) Now, let d 2 (x, y) = (x − y) 2 and rewrite (2) aŝ where D o and D e denote observed and expected disagreement, respectively. This is Krippendorff's α for the squared Euclidean distance (which is not a metric but a Bregman divergence (Bregman, 1967)). (We will henceforth refer to this version of α as Alpha(SED) or α sed .) As we mentioned above, this form of Krippendorff's α arises quite naturally when the data at hand conform to the one-way mixedeffects ANOVA model, for which agreement corresponds to positive correlation. More generally, Krippendorff recommends this form of α for the interval level of measurement. For other levels of measurement, Krippendorff presents other distance functions d 2 (several possibilities are shown in Table 1). Note that package krippendorffsalpha supports user-defined distance functions as well as the interval, nominal, and ratio distance functions shown in the table. 2.2. Alpha as a multiresponse permutation procedure. In any case, (3) is nonparametric for arbitrary d 2 since then the estimatorα does not usually correspond to a well-defined population parameter α. This more general form of Krippendorff's α is in fact a special case of the so called multiresponse permutation procedure (MRPP). The MRPPs form a class of permutation methods for discerning differences among two or more groups in one or more dimensions (Mielke and Berry, 2007). Note, however, that although α can be viewed as an MRPP (as we are about to show), α has been modified for the purpose of measuring agreement rather than discerning differences. To show that Krippendorff's α is an MRPP, we first present the general form of the MRPP. Following Mielke and Berry, let Ω = {ω 1 , . . . , ω n } be a finite sample that is representative of some population of interest, let S 1 , . . . , S a+1 denote a partition of Ω into a+1 disjoint groups, and let ρ be a metric that makes sense for the objects of Ω. (Strictly speaking, ρ need not be a metric; a symmetric distance function will suffice.) To ease notation a bit, let ρ jk ≡ ρ(ω j , ω k ). Then the MRPP statistic can be written as is the average distance between distinct pairs of objects in group S i ; n i ≥ 2 is the number of objects in group i; l = a i=1 n i ; n a+1 = n − l ≥ 0 is the number of remaining unclassified objects in group S a+1 ; and 1 i is the indicator function for membership in group i.
Note that this formulation is quite general since the objects in Ω can be scalars, vectors, or more exotic objects; and we are free to choose the metric and the weights. In the case of Krippendorff's α we can produce δ = D o by letting ρ = d 2 for some appropriately chosen distance function d 2 , and choosing weights C i = 1/n u .

A parametric generalization of Alpha(SED).
In the preceding sections we generalized α sed in a nonparametric fashion by substituting other notions of distance for the squared Euclidean distance. Now we will present a fully parametric generalization of α sed , namely, Sklar's ω (Hughes, 2018).
The statistical model underpinning Sklar's ω is a Gaussian copula model (Xue-Kun Song, 2000). The most general form of the model is given by where Ω is a correlation matrix, Φ is the standard Gaussian cdf, and F i is the cdf for the ith outcome Y i . Note that U = (U 1 , . . . , U n ) is a realization of the Gaussian copula, which is to say that the U i are marginally standard uniform and exhibit the Gaussian correlation structure defined by Ω. Since U i is standard uniform, applying the inverse probability integral transform to U i produces outcome Y i having the desired marginal distribution F i .
To see that the one-way mixed-effects ANOVA model (and hence α sed ) is a special case of Sklar's ω, let the copula correlation matrix Ω be block diagonal, where the ith block corresponds to the ith unit (i = 1, . . . , n u ) and has a compound symmetry structure. That is, where Complete the specification by letting F ij (i = 1, . . . , n u ; j = 1, . . . , n c ) be the cdf for the Gaussian distribution with mean µ and variance σ 2 . Then ω = α, the intraclass correlation coefficient.

Inference for Krippendorff's Alpha
Mielke and Berry describe hypothesis testing for MRPPs. Specifically, they discuss three approaches: permutation, Monte Carlo resampling, and Pearson type III moment approximation; the latter of which has significant advantages. For Krippendorff's α, though, we are interested not in hypothesis testing but in interval estimation. This can be done straightforwardly and efficiently using Monte Carlo resampling. Since D e is invariant to permutation of the scores, our resampling procedure focuses on D o only. The algorithm proceeds as follows.
(1) Collect the scores in an n u × n c matrix, A, where each row corresponds to a unit.
(2) For i ∈ {1, . . . , n b }, form matrix A i by sampling, with replacement, n u rows from A. We carried out a number of realistic simulation experiments and found that this approach to interval estimation performs well in a wide variety of circumstances. When the true distribution ofα is (at least approximately) symmetric, Gwet's closed-form expression forV(α), which is implemented (for categorical data only) in package irrCAC, also performs well. By contrast, we found that the bootstrapping procedure recommended by Krippendorff (2016), which is implemented in packages kripp.boot and icr, generally performs rather poorly, producing intervals that are far too narrow (e.g., 95% intervals achieve 74% coverage).

Robustness and interpretation
For some levels of measurement, one may, in the interest of robustness, be tempted to replace squares with absolute values (in the distance function d 2 ). This would be advantageous if one aimed to do hypothesis testing. But for Krippendorff's α, using absolute values instead of squares proves disastrous, for the resulting estimatorα is substantially negatively biased and tends to lead to erroneous inference regarding agreement. All is not lost, however, since package krippendorffsalpha provides a means of investigating the influence onα of any unit or coder (see the next section for examples).

Illustrations
Here we illustrate the use of krippendorffsalpha by applying Krippendorff's α to a couple of datasets. We will interpret the results according to the ranges given in Table 2, but we suggest-as do Krippendorff and others (Artstein and Poesio, 2008;Landis and Koch, 1977)-that an appropriate reliability threshold may be context dependent.
5.1. Nominal data analyzed previously by Krippendorff. Consider the following data, which appear in (Krippendorff, 2013  Note that the scores for all units except the sixth are constant or nearly so. This suggests near-perfect agreement, and so we should expectα to be greater than 0.8.
To apply Krippendorff's α to these data, first we load package krippendorffsalpha.
Now we create the dataset as a matrix such that each row corresponds to a unit and each column corresponds to a coder.
R> nominal = matrix(c (1,2,3,3,2,1,4,1,2 Next we apply Krippendorff's α for the nominal level of measurement. If argument level is set to "nominal", the discrete metric d 2 (x, y) = 1{x = y} is used by default. We do a bootstrap with sample size n b = 1,000 (argument confint defaults to TRUE, and control parameter bootit defaults to 1,000). We set control parameter parallel equal to FALSE because the dataset is too small to warrant parallelization of the bootstrap computation. Finally, we set argument verbose equal to TRUE so that a progress bar is shown during the bootstrap computation. The computation took less than one second. As is customary in R, one can view a summary by passing the fit object to summary.krippendorffsalpha, an S3 method. If krippendorffs.alpha was called with confint = TRUE, summary displays a 95% confidence interval by default. The confidence level can be specified using argument conf.level. In any case, the quantile method (Davison and Hinkley, 1997) is used to estimate the confidence limits. Any arguments passed to summary.krippendorffsalpha via . . . are passed on to R's quantile function; this allows the user to control, for example, how the sample quantiles are computed.  We see thatα = 0.74 and α ∈ (0.46, 1.00). This point estimate indicates only substantial agreement, which is not what we expected. At least the interval is consistent with near-perfect agreement, but we should not take this interval too seriously since the interval is rather wide (owing to the small size of the dataset).
Perhaps the substantial disagreement for the sixth unit was influential enough to yieldα ≤ 0.8. We can use influence.krippendorffsalpha, another S3 method, to investigate. This function, like other R versions of influence (e.g., influence.lm, influence.glm), computes DFBETA statistics (Young, 2017), as illustrated below.
R> (inf.6 = influence(fit.full, units = 6)) $dfbeta.units 6 -0.1141961 Leaving out the sixth unit yields a DFBETA statistic of -0.11, which implies thatα would have been 0.86. This is consistent with our initial hypothesis. We see that excluding the sixth unit leads to α ∈ (0.66, 1.00). The new 95% interval was returned by S3 method confint.krippendorffsalpha, whose level argument defaults to 0.95, in keeping with R's other confint methods. Note that confint.krippendorffsalpha, like summary.krippendorffsalpha, passes any . . . arguments on to the quantile function.
We conclude this example by producing a visual display of our results (Figure 3). (The figure was produced via a call to S3 method plot.krippendorffsalpha, which in turn calls hist and abline, and does not show a kernel density estimate. Function plot.krippendorffsalpha is capable of producing highly customized plots; see the package documentation for details.) Sinceα is close to 1 and the dataset is small, the bootstrap distribution is substantially skewed to the left. Thus these data provide a textbook example of the importance of bootstrapping.

Nominal Data
Bootstrap Estimates Since the dataset used in this example has missing values, we take this opportunity to explain how the package handles missingness. First, the scores for a given unit of analysis are included in the computation only if two or more scores are present for that unit. Otherwise the unit's row of the data matrix is simply ignored. Second, if two or more scores are present for a given unit, each NA for that unit is ignored in the computations for that row. This is handled both by the loop (adjusted denominator) and by the distance function, which should return 0 if either of its arguments is NA. In the next example we illustrate this by way of a user-defined distance function, and of course the package's built-in distance functions take the same approach.

5.2.
Interval data from an imaging study of hip cartilage. The data for this example, some of which appear in Figure 4, are 323 pairs of T2* relaxation times (a magnetic resonance quantity) for femoral cartilage (Nissi et al., 2015) in patients with femoroacetabular impingement (Figure 5), a hip condition that can lead to osteoarthritis. One measurement was taken when a contrast agent was present in the tissue, and the other measurement was taken in the absence of the agent. The aim of the study was to determine whether raw and contrast-enhanced T2* measurements agree closely enough to be interchangeable for the purpose of quantitatively assessing cartilage health.

12JOHN HUGHES DEPARTMENT OF STATISTICS THE PENNSYLVANIA STATE UNIVERSITY
First we load the cartilage data, which are included in the package. The cartilage data are stored in a data frame; we convert the data frame to a matrix, which is the format required by krippendorffs.alpha.
R> data(cartilage) R> cartilage = as.matrix(cartilage) Now we computeα for the interval level of measurement, i.e., squared Euclidean distance. We also produce a bootstrap sample of size 10,000. Since this dataset is much larger than the dataset analyzed in the preceding section, we parallelize the bootstrap computation. We use three CPU cores (of the four available on the author's computer). Setting argument verbose to TRUE causes the fitting function to display a progress bar once again. The computation took five seconds to complete.

|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=05s
A call of function summary.krippendorffsalpha produced the output shown below. We see thatα = 0.84 and α ∈ (0.81, 0.86). Thus these data suggest that raw T2* measurements agree almost perfectly with contrast-enhanced T2* measurements, perhaps rendering gadolinium-based contrast agents (GBCAs) unnecessary in T2*-based cartilage assessment. This finding could have clinical significance since the use of GBCAs is not free of risk to patients, especially pregnant women and patients with impaired kidney function. For much additional information regarding the potential risks associated with the use of GBCAs, we refer the interested reader to the University of California, San Francisco's policy on MRI with contrast: https://radiology.ucsf.edu/patient-care/patient-safety/ contrast/mri-with-contrast-gadolinium-policy. Figure 6 provides a visual display of the cartilage results. The histogram and kernel density estimate show the expected large-sample behavior ofα, i.e., the estimator is approximately Gaussian-distributed and has a small variance.

Cartilage Data (SED)
Bootstrap Estimates  Figure 6. A plot of the results from our analysis of the cartilage data. The histogram and kernel density estimate (dotted black curve) show the bootstrap sample, the solid orange line marks the value ofα, and the dashed blue lines mark the 95% confidence limits.
We mentioned above that attempting to robustify Krippendorff's α by using absolute values in place of squares may prove problematic. This is evident for the cartilage data, as we now demonstrate.
First, define a new distance function as follows. Note that any user-defined distance function must deal explicitly with NAs if the data at hand exhibit missingness. There are no missing values in the cartilage data, but we illustrate the handling of NA anyway. Control parameter 'type' must be "SOCK", "PVM", "MPI", or "NWS". Setting it to "SOCK".

Summary and discussion
In this article we described Krippendorff's α methodology for measuring agreement, and illustrated the use of R package krippendorffsalpha. We first established α's context among statistical procedures. Specifically, the one-way mixedeffects ANOVA model provides a natural, intuitive genesis for α as the intraclass correlation coefficient. This form of α can be generalized in a parametric fashion to arrive at Sklar's ω, or in a nonparametric fashion to arrive at the form of α presented by Krippendorff, which is a special case of the multiresponse permutation procedure.
We demonstrated the use of krippendorffsalpha version 1.1 by analyzing two datasets: a nominal dataset previously analyzed by Krippendorff, and a sample of raw and contrast-enhanced T2* values from an MRI study of hip cartilage. These analyses highlighted the benefits of the package, which include the use of S3 methods, parallel bootstrap computation, support for user-defined distance functions, and a means of identifying influential units and/or coders.

Computational details
The results in this paper were obtained using R 4.0.3 for macOS, and the pbapply 1.4-2 package. R itself and all packages used (save kripp.boot) are available from the Comprehensive R Archive Network (CRAN) at http://CRAN. R-project.org. Package krippendorffsalpha may be downloaded from CRAN or from the author's GitHub repository, which can be found at https://github. com/drjphughesjr/krippendorffsalpha. Information about the author's other R packages can be found at http://www.johnhughes.org/software.html.