equate: An R Package for Observed-Score Linking and Equating

The R package equate contains functions for observed-score linking and equating under single-group, equivalent-groups, and nonequivalent-groups with anchor test(s) designs. This paper introduces these designs and provides an overview of observed-score equating with details about each of the supported methods. Examples demonstrate the basic functionality of the equate package.


Introduction
Equating is a statistical procedure commonly used in testing programs where administrations across more than one occasion and more than one examinee group can lead to overexposure of items, threatening the security of the test. In practice, item exposure can be limited by using alternate test forms; however, multiple forms lead to multiple score scales that measure the construct of interest at differing levels of difficulty. The goal of equating is to adjust for these differences in difficulty across alternate forms of a test, so as to produce comparable score scales.
Equating defines a functional statistical relationship between multiple test score distributions and thereby between multiple score scales. When the test forms have been created according to the same specifications and are similar in statistical characteristics, this functional relationship is referred to as an equating function, and it serves to translate scores from one scale directly to their equivalent values on another. The term linking refers to test forms which have not been created according to the same specifications, for example, forms which differ in length or content; in this case, the linked scales are considered similar but not interchangeable; they are related to one another via a linking function. Specific requirements for equating The equate package, available on the Comprehensive R Archive Network (CRAN) at https: //CRAN.R-project.org/package=equate, is designed for observed-score linking and equating. It differs from other packages primarily in its overall structure and usability, its plotting and bootstrapping capabilities, and its inclusion of more recently developed equating and linking functions such as the general-linear, synthetic, and circle-arc functions, and traditional methods such as Levine, Tucker, and Braun/Holland. Equating with multiple anchor tests and external covariates is also supported, as demonstrated below. Linking and equating are performed using a simple interface, and plotting and summary methods are provided to facilitate the comparison of results and the examination of bootstrap and analytic equating error. Sample data and detailed help files are also included. These features make the package useful in teaching, research, and operational testing contexts. This paper presents some basic linking and equating concepts and procedures. Equating designs are first discussed in Section 2. In Section 3, linear and nonlinear observed-score linking and equating functions are reviewed. In Section 4, methods are presented for linking and equating when examinee groups are not equivalent. Finally, in Section 5, the equate package is introduced and its basic functionality is demonstrated using three data sets.

Equating designs
Observed-score linking and equating procedures require data from multiple test administrations. An equating design specifies how the test forms and the individuals sampled to take them differ across administrations, if at all. For simplicity, in this paper and in the equate package, equating designs are categorized as either involving a single group, equivalent groups, or nonequivalent groups of examinees, and test forms are then constructed based on the type of group(s) sampled.
In the single-group design, one group, sampled from the target population T , takes two different test forms X and Y , optionally with counterbalancing of administration orders (one group takes X first, the other takes Y first). Any differences in the score distributions on X and Y are attributed entirely to the test forms themselves, as group ability is assumed to be constant; thus, if the distributions are not the same, it is because the test forms differ in difficulty. Related to the single-group design is the equivalent-groups design, where one random sample from T takes X and another takes Y . Because the samples are taken randomly, group ability is again assumed to be constant, and any differences in the score distributions are again identified as form difficulty differences.
Without equivalent examinee groups, two related problems arise: (1) the target population must be defined indirectly using samples from two different examinee populations, P and Q; and (2) the ability of these groups must then be accounted for, as ability differences will be a confounding factor in the estimation of form difficulty differences. In the nonequivalentgroups design these issues are both addressed through the use of what is referred to as an anchor test, V , a common measure of ability available for both groups. All non-equivalence in ability is assumed to be controlled or removed via this common measure. External covariates, such as scores from other tests, can also be used to control for group differences.
Equating procedures were initially developed using the single-group and equivalent-groups designs. In this simpler context, the traditional equating functions include mean, linear, and equipercentile equating; these and other equating functions are reviewed in Section 3. More complex procedures have been developed for use with the nonequivalent-groups design; these equating methods are presented in Section 4. Unless otherwise noted, additional details on the functions and methods described below can be found in Kolen and Brennan (2014).

Equating functions
Procedures for equating test forms to a common scale are referred to here and in the equate package as different types of equating functions. The equating function defines the equation for a line that expresses scores on one scale, or axis, in terms of the other. The available types of equating functions are categorized as straight-linear (i.e., linear), including identity, mean, and linear equating, and curvilinear (i.e., nonlinear), including equipercentile and circle-arc equating. The straight-line types differ from one another in intercept and slope, and the curvilinear lines differ in the number of coordinates on the line that are estimated, whether all of them or only one. Combinations of equating lines, referred to here as composite functions, are also discussed.
The goal of equating is to summarize the difficulty difference between X to Y . As shown below, each equating type makes certain assumptions regarding this difference and how it does or does not change across the X and Y score scales. These assumptions are always expressed in the form of a line within the coordinate system for the X and Y scales.

Identity functions
Linear functions are appropriate when test form difficulties change linearly across the score scale, by a constant b and rate of change a. Scores on X are related to Y as In the simplest application of Equation 1, the scales of X and Y define the line. Coordinates for scores of x and y are found based on their relative positions within each scale: Here, (x 1 , y 1 ) and (x 2 , y 2 ) are coordinates for any two points on the line defined by the scales of X and Y , for example, the minimum and maximum possible scale values. Solving Equation 2 for y results in the identity linking function: equate: Observed-Score Linking and Equating in R The intercept b can also be defined using the slope a and any pair of X and Y coordinates ( where j = 1, 2, . . . , J indexes the points on scale X and k = 1, 2, . . . , K indexes the points on scale Y . The identity linking function is then expressed as When the scales of X and Y are the same, a = 1 and b = 0, and Equation 7 reduces to the identity equating function:

Mean functions
In mean linking and equating, form difficulty differences are estimated by the mean difference µ Y − µ X . Equation 7 is used to define a line that passes through the means of X and Y , rather than the point (x j , y k ). The intercept from Equation 6 is expressed as The mean linking function is then where a is found using Equation 4. When the scales of X and Y are the same, the slope a is 1, which leads to the mean equating function: In mean equating, coordinates for the line are based on deviation scores: In mean linking, coordinates are based on deviation scores relative to the scales of X and Y :

Linear functions
The linear linking and equating functions also assume that the difficulty difference between X and Y changes by a constant amount a across the score scale. However, in linear equating the slope is estimated using the standard deviations of X and Y as The linear linking and equating functions are defined as In both the linear linking and equating functions, coordinates for the line are based on standardized deviation scores:

General linear functions
The identity, mean, and linear linking and equating functions presented above can all be obtained as variations of a general linear function glin Y (x) (Albano 2015). The general linear function is defined based on Equation 1 as Here, α is a general scaling parameter that can be estimated using σ, ∆, another fixed value, or weighted combinations of these values. β is a general centrality parameter that can be estimated using µ, x j or y k , other values, or weighted combinations of these values. Applications of the general linear function are discussed below and in Albano (2015).

Equipercentile functions
Equipercentile linking and equating define a nonlinear relationship between score scales by setting equal the cumulative distribution functions for X and Y : F (x) = G(y). Solving for y produces the equipercentile linking function: which is also the equipercentile equating function equipe Y (x). When the score scales are discrete, which is often the case, the cumulative distribution function can be approximated using percentile ranks. This is a simple approach to continuizing the discrete score distributions (for details, see Kolen and Brennan 2014, Chapter 2). Kernel equating, using Gaussian kernels, offers a more flexible approach to continuization (von Davier, Holland, and Thayer 2004), but differences between the methods tend to be negligible. The percentile rank method is currently used in the equate package. The equipercentile equivalent of a form-X score on the Y scale is calculated by finding the percentile rank in X of a particular score, and then finding the form-Y score associated with that form-Y percentile rank.
Equipercentile equating is appropriate when X and Y differ nonlinearly in difficulty, that is, when difficulty differences fluctuate across the score scale, potentially at each score point. Each coordinate on the equipercentile curve is estimated using information from the distributions of X and Y . Thus, compared to identity, mean, and linear equating, equipercentile equating is more susceptible to sampling error because it involves the estimation of as many parameters as there are unique score points on X.
Smoothing methods are typically used to reduce irregularities due to sampling error in either the score distributions or the equipercentile equating function itself. Two commonly used smoothing methods include polynomial loglinear presmoothing (Holland and Thayer 2000) and cubic-spline postsmoothing (Kolen 1984). The equate package currently supports loglinear presmoothing via the glm function. Details are provided below.

Circle-arc functions
Circle-arc linking and equating (Livingston and Kim 2009) also define a nonlinear relationship between score scales; however, they utilize only three score points in X and Y to do so: the low and high points, as defined above for the identity function, and a midpoint (x j , y k ). On their own, the low and high points define the identity linking function id Y (x), a straight line.
When (x j , y k ) does not fall on the identity linking line, it can be connected to (x 1 , y 1 ) and (x 2 , y 2 ) by the circumference of a circle with center (x c , y c ) and radius r.
There are multiple ways of solving for (x c , y c ) and r based on the three known points (x 1 , y 1 ), (x j , y k ), and (x 2 , y 2 ). For example, the center coordinates can be found by solving the following system of equations: Subtracting Equation 23 from (21) and (22) and rearranging terms leads to the following linear system: The center coordinates can then be obtained by plugging in the known values for (x 1 , y 1 ), (x j , y k ), and (x 2 , y 2 ) and again combining equations. The center and any other coordinate pair, e.g., (x 1 , y 1 ), are then used to find the radius: Finally, solving Equation 26 for y results in the circle-arc linking function: where the second quantity, under the square root, is added to y c when y k > id Y (x j ) and subtracted when y k < id Y (x j ). The circle-arc equating function circe Y (x) is obtained by using ide Y (x j ) in place of id Y (x j ) above. Livingston and Kim (2010) refer to the circle connecting (x 1 , y 1 ), (x j , y k ), and (x 2 , y 2 ) as symmetric circle-arc equating. They also present a simplified approach, where the circle-arc function is decomposed into the linear component defined by (x 1 , y 1 ) and (x 2 , y 2 ), which is the identity function, and the circle defined by the points (x 1 , y 1 −id Y (x 1 )), (x j , y k −id Y (x j )), and (x 2 , y 2 − id Y (x 2 )). These low and high points reduce to (x 1 , 0) and (x 2 , 0), and the center coordinates can then be found as and where y * k = y k − id Y (x j ). Equation 26 is used to find the radius. Then, the simplified circlearc function is the combination of the resulting circle-arc circ * Y (x) and the identity function:

Composite functions
The circle-arc linking and equating functions involve a curvilinear combination of the identity and mean functions, where the circle-arc overlaps with the identity function at the low and high points, and with the mean function at the midpoint (µ X , µ Y ). A circle then defines the coordinates that connect these three points. This is a unique example of what is referred to here as a composite function.
The composite linking function is the weighted combination of any linear and/or nonlinear linking or equating functions: where w h is a weight specifying the influence of function link hY (x) in determining the composite.
Equation 31 is referred to as a linking function, rather than as equating function, because it will typically not meet the symmetry requirement of equating. For symmetry to hold, the inverse of the function that links X to Y must be the same as the function that links Y to X, that is, comp −1 Y (x) = comp X (y), which is generally not true when using Equation 31. Holland and Strawderman (2011) show how symmetry can be maintained for any combination of two or more linear functions. The weighting system must be adjusted by the slopes for the linear functions being combined, where the adjusted weight W h is found as Here, a h is the slope for a given linear function link h , and p specifies the type of L p -circle with which symmetry is defined. For details, see Holland and Strawderman (2011).

Equating methods
The linking and equating functions presented above are defined in terms of a single target population T , and they are assumed to generalize to this population. A subscript, e.g., X T , is omitted for simplicity; it is presumed that X = X T and Y = Y T . In the nonequivalent-groups design, scores come from two distinct populations, referred to here as populations P and Q.
Because the target population is not directly sampled, the linking and equating functions are redefined in terms of a weighted combination of P and Q, where T = w P P + w Q Q and w P and w Q are proportions that sum to 1. This mixture of P and Q is referred to as the synthetic population (Braun and Holland 1982).
Linear equating is presented for the synthetic population first. All of the means and standard deviations in Equation 15 are estimated as weighted combinations of estimates from P and Q, where Because X is not administered to Q and Y is not administered to P , the terms µ X Q , µ Y P , σ 2 X Q , and σ 2 X Q in Equations 33 through 36 are obtained using available information for X, Y , and the anchor test V . This results in the following synthetic parameter estimates (for details, see Kolen and Brennan 2014): The γ terms in Equations 37 through 40 represent the relationship between total scores on X and Y and the respective anchor scores on V . γ P and γ Q are used along with the weights to adjust the observed µ and σ 2 for X and Y in order to obtain corresponding estimates for the synthetic population. For example, when w P = 0 and w Q = 1, µ Y = µ Y Q , and conversely µ X Q will be adjusted the maximum amount when obtaining µ X . The same would occur with the estimation of synthetic variances. Furthermore, the adjustments would be completely removed if populations P and Q did not differ in ability, where µ V P = µ V Q and σ 2 V P = σ 2 V Q . A handful of techniques have been developed for estimating the linear γ terms required by Equations 37 through 40, and the terms required for equipercentile equating, as described below. These techniques all make certain assumptions about the relationships between total scores and anchor scores for populations P and Q. The techniques are referred to here as equating methods. The equate package supports the Tucker, nominal weights, Levine observed-score, Levine true-score, Braun/Holland, frequency estimation, and chained equating methods (although chained equating does not rely on γ, it does make assumptions about the relationship between total and anchor scores). The Tucker, nominal weights, Braun/Holland, and frequency estimation methods are also available for use with multiple anchor tests; see Appendix A. Table 1 shows the supported methods that apply to each equating type.

Tucker method
In Tucker equating (Levine 1955) the relationship between total and anchor test scores is Method Type nominal tucker levine braun frequency chained mean defined in terms of regression slopes, where γ P is the slope resulting from the regression of X on V for population P , and γ Q the slope from a regression of Y on V for population Q: The Tucker method assumes that across populations: (1) the coefficients resulting from a regression of X on V are the same, and (2) the conditional variance of X given V is the same. These assumptions apply to the regression of Y on V and the covariance of Y given V as well.
They also apply to the regression of X or Y on multiple anchor tests (e.g., Angoff 1984); see Appendix A.1 for details.

Nominal weights method
Nominal weights equating is a simplified version of the Tucker method where the total and anchor tests are assumed to have similar statistical properties and to correlate perfectly within populations P and Q. In this case the γ terms can be approximated by the ratios where K is the number of items on the test. See Babcock, Albano, and Raymond (2012) for a description and examples with a single anchor. When using multiple anchor tests, a γ term is again estimated for each anchor test, as in the multi-anchor Tucker method; see Appendix A.2.

Levine method
Assumptions for the Levine (Levine 1955) observed-score method are stated in terms of true scores (though only observed scores are used), where, across both populations: (1) the correlation between true scores on X and V is 1, as is the correlation between true scores on Y and V ; (2) the coefficients resulting from a linear regression of true scores for X on V are the same, as with true scores for Y on V ; and (3) measurement error variance is the same (across populations) for X, Y , and V . These assumptions make possible the estimation of γ as which are the inverses of the respective regression slopes for V on X and V on Y . The Levine true-score method is based on the same assumptions as the observed-score method; however, it uses a slightly different linear equating function in place of Equation 15: with γ defined by Equation 43. Hanson (1991) and Kolen and Brennan (2014) provide justifications for using the Levine true-score method.

Frequency estimation method
The frequency estimation or poststratification method is used in equipercentile equating under the nonequivalent-groups design. It is similar to the methods described above in that it involves a synthetic population. However, in this method full score distributions for the synthetic population taking forms X and Y are required. When the assumptions are made that (1) the conditional distribution of total scores on X for a given score point in V is the same across populations, and (2) the conditional distribution of total scores on Y for a given score point in V is the same across populations, the synthetic distributions can be obtained as and Here, P(x), P(y), and P(v) denote the distribution functions for forms X, Y , and V respectively. Percentile ranks can be taken for the cumulative versions of Equations 45 and 46 to obtain Equation 20. The frequency estimation method can also accommodate multiple anchor tests and covariates; see Appendix A.3.

Braun/Holland method
As a kind of extension of the frequency estimation method, the Braun/Holland method (Braun and Holland 1982) defines a linear function relating X and Y that is based on the means and standard deviations of the synthetic distributions obtained via frequency estimation. Thus the full synthetic distributions are estimated, as with frequency estimation, but only in order to obtain their means and standard deviations.

Chained method
Finally, chained equating (Livingston, Dorans, and Wright 1990) can be applied to both linear and equipercentile equating under the nonequivalent-groups with anchor test design. The chained method differs from all other methods discussed here in that it does not explicitly reference a synthetic population. Instead, it introduces an additional equating function in the process of estimating score equivalents; see Appendix B for details. For both linear and equipercentile equating the steps are as follows: 1. Define the function relating X to V for population P, link V P (x).

Define the function relating
Chained methods are based on the assumptions that (1) the equating of X to V is the same for P and Q, and (2) the equating of V to Y is the same for P and Q.

Methods for general linear functions
The general linear equating function can be utilized with any combination of weighted means and standard deviations from Equations 33 through 36. Thus, any methods for nonequivalent groups that estimate means or means and standard deviations for the synthetic population can be implemented within the general linear function. In the equate package, these methods are currently Tucker, nominal-weights, Levine observed-score, and Braun/Holland. See Albano (2015) for examples. Composites of these methods can also be obtained.

Methods for circle-arc functions
As discussed above, the circle-arc equating function combines a linear with a curvilinear component based on three points in the X and Y score distributions. Although all three points can be obtained using the general linear function, the first and third of these points are typically determined by the score scale whereas the midpoint is estimated. Equating methods used with circle-arc equating in the nonequivalent-groups design apply only to estimation of this midpoint. Livingston and Kim (2009) demonstrated chained linear equating of means, under a nonequivalent-groups design. The midpoint could also be estimated using other linear methods, such as Tucker or Levine.
Note that circle-arc equating is defined here as an equating type, and equating methods are used to estimate the midpoint. When groups are considered equivalent (i.e., an anchor test is not used) equating at the midpoint is simply mean equating, as mentioned above (replace x with µ X in Equation 15 to see why this is the case). With scores on an anchor test, both Tucker and Levine equating at the midpoint also reduce to mean equating. However, chained linear equating at the midpoint differs from chained mean (see Appendix B).

Sample data
The equate package includes three sample data sets. The first, ACTmath, comes from two administrations of the ACT mathematics test, and is used throughout Kolen and Brennan (2014). The test scores are based on an equivalent-groups design and are contained in a three-column data frame where column one is the 40-point score scale and columns two and three are the number of examinees for X and Y obtaining each score point.
The second data set, KBneat, is also used in Kolen and Brennan (2014). It contains scores for two forms of a 36-item test administered under a nonequivalent-groups design. A 12-item anchor test is internal to the total test, where anchor scores contribute to an examinee's total score. The number of non-anchor items that are unique to each form is 24, and the highest possible score is 36. KBneat contains a separate total and anchor score for each examinee. It is a list of length two where the list elements x and y each contain a two-column data frame of scores on the total test and scores on the anchor test.
The third data set, PISA, contains scored cognitive item response data from the 2009 administration of the Programme for International Assessment (PISA). Four data frames are included in PISA: PISA$students contains scores on the cognitive assessment items in math, reading, and science for all 5233 students in the USA cohort; PISA$booklets contains information about the structure of the test design, where multiple item sets, or clusters, were administered across 13 test booklets; PISA$items contains the cluster, subject, maximum possible score, item format, and number of response options for each item; and PISA$totals contains a list of cluster total scores for each booklet, calculated using PISA$students and PISA$booklets. For additional details, see the PISA help file which includes references to technical documentation.

Preparing score distributions
The equate package analyzes score distributions primarily as frequency table objects with class 'freqtab'. For example, to equate the ACTmath forms, they must first be converted to frequency tables as follows. The 'freqtab' class stores frequency distributions as table arrays, with a dimension for each of the variables included in the distribution. The function as.freqtab is used above because ACTmath already contains tabulated values; this code simply restructures the scales and counts for the two test forms and gives them the appropriate attributes. When printed to the console, 'freqtab' objects are converted to data.frames. They are summarized with the summary method. The constructor freqtab creates a frequency table from a vector or data frame of observed scores. With an anchor test, this becomes a bivariate frequency table. Bivariate distributions contain counts for all possible score combinations for the total and anchor scores. Multivariate distributions, e.g., containing scores on multiple anchor tests and external covariates, are also supported.
R> neat.x <-freqtab(KBneat$x, scales = list(0:36, 0:12)) R> neat.y <-freqtab(KBneat$y, scales = list(0:36, 0:12)) Finally, freqtab can also be used to sum scored item responses, and then tabulate the resulting total scores. In this case, the list items must contain vectors of the columns over which total scores should be calculated. For example, the following syntax creates a frequency   A basic plot method is provided for the 'freqtab' class. Univariate frequencies are plotted as vertical lines for the argument x, similar to a bar chart, and as superimposed curves for the argument y. When y is a matrix, each column of frequencies is added to the plot as a separate line. This feature is useful when examining smoothed frequency distributions, as demonstrated below. When x is a bivariate frequency table, a scatter plot with marginal frequency distributions is produced. See Figure 1 for an example of a univariate plot, and Figure 2 for an example of a bivariate plot.

Presmoothing
The distributions in Figures 1 and 2 contain irregularities in their shapes that likely result, to some extent, from sampling error. The population distributions that these samples estimate are expected to be smoother, with less jaggedness between adjacent scores. Three methods are available for smoothing sample distributions with the objective of reducing these irregularities. The first, frequency averaging (Moses and Holland 2008) replaces scores falling below jmin with averages based on adjacent scores. This is implemented with smoothmethod = "average" in the presmoothing function. The second, adds a small relative frequency (again, jmin) to each score point while adjusting the probabilities to sum to one (as described by Kolen and Brennan 2014, p. 46). This is implemented using smoothmethod = "bump" in the presmoothing function.
The third method of reducing irregularities in sample distributions is polynomial loglinear smoothing. Appendix C contains details on the model formula itself. In the equate package, loglinear models are fit using the presmoothing function with smoothmethod = "loglinear", which calls the generalized linear model (glm) function in R. Models can be fit in three different ways. The first way is with a formula object, as follows.
R> presmoothing(~poly(total, 3, raw = TRUE) + poly(anchor, 3, raw = TRUE) + + total:anchor, data = neat.x) This is similar to the approach used in other modeling functions in R, but with two restrictions: (1) the dependent variable, to the left of the~, is set to be the score frequencies contained in data, and it does not need to be specified in the formula; and (2) the intercept is required and will be added if it has been explicitly removed in the formula.
The formula smoothing method is efficient to write for simpler models, but it can be cumbersome for more complex models containing multiple interaction terms. The second way to specify the model is with a matrix of score functions (scorefun) similar to a model.matrix in R, where each column is a predictor variable in the model, as follows.
R> neat.xsf <-with(as.data.frame(neat.x), cbind(total, total^2, + total^3, anchor, anchor^2, anchor^3, total * anchor)) R> presmoothing(neat.x, smooth = "loglinear", scorefun = neat.xsf) The object neat.xsf is a matrix containing the total and anchor score scales to the first, second, and third powers, and the interaction between the two. The presmoothing results based on this score function are the same as those for the formula method above. One benefit of creating the score function externally is that it can be easily modified and used with other models. It can also include variables not contained within the raw frequency table, neat.x, in this example. The formula interface is limited in this respect, as the data argument must be a frequency table, and it cannot include variables besides the total and anchor scores.
The most efficient approach to specify a loglinear model in the presmoothing function is by including the degrees of the highest polynomial terms for each variable at each level of interaction. For example, in the formula and score function methods above, terms are included for both the total and anchor tests up to the third power and for the two-way interaction to the first power. This is specified compactly using degrees = list(c(3, 3), c(1, 1)), which can be reduced to degrees = list(3, 1).
R> neat.xs <-presmoothing(neat.x, smooth = "log", degrees = list(3, 1)) This functionality extends to multivariate distributions. For example, three way interactions would be specified by including a third vector in the degrees list.
For the bivariate example, the smoothed distributions in Figure 3 can be compared to the unsmoothed ones in Figure 2. Figure 4 superimposes the smoothed frequencies on the unsmoothed marginal distributions for a more detailed comparison of the different models. Descriptive statistics show that the smoothed distributions match the unsmoothed in the first three moments.
R> neat.xsmat <-presmoothing(neat.x, smooth = "loglinear", + degrees = list(3, 1), stepup = TRUE) R> plot(neat.xs) R> plot(neat.x, neat.xsmat, ylty = 1:4) R> round(rbind(x = summary(neat.x), xs = summary(neat.xs)), 2) The presmoothing function is used above to compare results from a sequence of nested models. The argument stepup = TRUE returns a matrix of fitted frequencies for models based on subsets of columns in scorefun, where the columns for each model can be specified with the argument models. The presmoothing function can also infer nested models when the degrees argument is used. In this case, terms are added sequentially for all variables within each level of interaction. For the example above, the first model in neat.xsmat includes the total and anchor scales to the first power, the second additionally includes both scales to the second power, and the third includes both to the third power. A fourth model contains the interaction term. The smoothed curves in the marginal distributions of Figure 4 show the loglinear smoothing results for each nested model that is fit in neat.xsmat. The legend text defines each smoothing curve using two numbers, the first for the level of interaction (1 for the first three models, and 2 for the fourth), and the second for the highest power included in a model (1, 2, and 3 for the main effects, and 1 for the interaction).
Using the argument compare = TRUE in presmoothing, an ANOVA table of deviance statistics is returned for sequentially nested models. Model fit is compared using functions from the base packages in R, which provide AIC (Akaike's information criterion), BIC (Bayesian information criterion), and likelihood ratio χ 2 tests. In the output below, AIC and BIC are smallest for the most complex model, labeled "Model 4", which also results in the largest decrease in deviance. The models being compared are the same as those contained in neat.xsmat.
R> presmoothing(neat.x, smooth = "loglinear", + degrees = list(c(3, 3), c(1, 1)), compare = TRUE) equate: Observed-Score Linking and Equating in R Finally, with choose = TRUE, the presmoothing function will automatically select the best fitting model and return a smoothed frequency distribution based on that model. The statistic for selection is indicated in the argument choosemethod, with options χ 2 ("chi"), AIC ("aic"), and BIC ("bic"). For "aic" and "bic", the model with the smallest value is chosen. For "chi", the most complex model with p value less than the argument chip is chosen, with default of 1 − (1 − 0.05) ( 1/(#models − 1)). This automatic model selection is useful in simulation and resampling studies where unique presmoothing models must be fit at each replication.

The equate function
Most of the functionality of the equate package can be accessed via the function equate, which integrates all of the equating types and methods introduced above. The equivalentgroups design provides a simple example: Besides the X and Y frequency tables, only the equating type, i.e., the requested equating function, is required.
R> equate(act.x, act.y, type = "mean") The nonequivalent-groups design is specified with an equating method, and smoothing with a smoothmethod.
An equating object such as neat.ef contains basic information about the type, method, design, smoothing, and synthetic population weighting for the equating, in addition to the conversion Here, the argument y passed to equate is the frequency estimation equipercentile equating object from above, which is an object of class 'equate'. Since the equating function from neat.ef relates scores on X to the scale of Y , anchor test scores are not needed for the examinees in newx. This feature provides a quick way to convert a score vector of any size from X to Y . Because this feature does not rely on the discrete concordance  Kim, von Davier, and Haberman (2008) refer to as synthetic linear linking. The argument symmetric = TRUE is used to adjust the weighting system so that the resulting function is symmetric. Figure 5 shows the composite line in relation to the identity and linear components.  The two forms differ in length and item type. R3R6 contains 30 items, one of which has a maximum possible score of 2, and the remainder of which are scored dichotomously. This results in a score scale ranging from 0 to 31. However, 14 of the 30 items in R3R6 were multiple-choice, mostly with four response options. The remaining items were either constructed-response or complex multiple-choice, where examinees were unlikely to guess the correct response. Thus, the lowest score expected by chance for R3R6 is 14/4 = 3.5. R5R7 contains 29 items, all of which are scored dichotomously. Eight of these items are multiple-choice with four response options and the remainder are constructed-response or complex multiple-choice, resulting in a lowest expected chance score of 8/4 = 2. The summary statistics above show that, despite having a slightly smaller score scale, the mean for R5R7 is slightly higher than for R3R5.

Linking with different scale lengths and item types
Results for linking R3R6 to R5R7 are compared here for five linking types: identity, mean, linear, circle-arc, and equipercentile with loglinear presmoothing (using the default degrees). By default, the identity linking component of each linear function is based on the minimum and maximum possible points for each scale, that is, (0, 0) and (31, 29). The low points were modified to be (3.5, 2) to reflect the lowest scores expected by chance.

Linking with multiple anchors and covariates
The PISA data are used here to demonstrate linking with multiple anchor tests. As noted above, the PISA data come from a cluster rotation design, where different groups of students, organized by test booklet, saw different clusters of math, reading, and science items. Data from booklet 4 are used here to create two pseudo forms for a NEAT design with an external covariate. The unique reading items for each form come from clusters R3 and R4, the first anchor comes from reading cluster R2, and the covariate from science cluster S2.
If multiple anchor tests and/or covariates are contained within a frequency table, they are automatically used by the equate function for any methods that support multi-anchor equating, and they are ignored for methods that do not support them. The following code conducts linking with covariates using the nominal weights, Tucker, and frequency estimation methods.

Bootstrapping
All but the identity linking and equating functions estimate a statistical relationship between score scales. Like any statistical estimate, equated scores are susceptible to bias and random sampling error, for example, as defined in Appendix D. Standard error (SE), bias, and root mean square error (RMSE) can be estimated in the equate package using empirical and parametric bootstrapping.
With the argument boot = TRUE, the equate function will return bootstrap standard errors based on sample sizes of xn and yn taken across reps = 100 replications from x and y. Individuals are sampled with replacement, and the default sample sizes xn and yn will match those observed in x and y. Equating is performed at each replication, and the estimated equating functions are saved. bias and RMSE can be obtained by including a vector of criterion equating scores via crit. Finally, the matrix of estimated equatings at each replication can be obtained with eqs = TRUE.
Parametric bootstrapping involves resampling as described above, but from a smoothed score distribution that is assumed to produce more reliable results with small samples (Kolen and Brennan 2014). In simulation studies this smoothed distribution is sometimes treated as a pseudo-population. Parametric bootstrapping is performed within the equate function by providing the optional frequency distributions xp and yp. These simply replace the sample distributions x and y when the bootstrap resampling is performed. Additionally, the bootstrap function can be used directly to perform multiple equatings at each bootstrap replication. SE, bias, and RMSE can then be obtained for each equating function using the same bootstrap data.
Note that the number of bootstrap replications, specified via the reps argument, can impact the stability of the results, with error estimates varying noticeably for replications below 100.
Bootstrapping studies vary widely in the number of replications utilized. It is recommended that no fewer than 100 be used. For more stable results, 500 to 1000 replications may be necessary, as computing time permits.
Parametric bootstrapping using the bootstrap function is demonstrated here for eight equatings of form X to Y in KBneat: Tucker and chained mean, Tucker and chained linear, frequency estimation and chained equipercentile, and Tucker and chained-linear circle-arc. Identity equating is also included. Smoothed population distributions are first created. Based on model fit comparisons, loglinear models were chosen to preserve 4 univariate and 2 bivariate moments in the smoothed distributions of X and Y . Plots are shown in Figures 8 and 9.
R> set.seed(131031) R> reps <-100 R> xn <-100 R> yn <-100 R> crit <-equate(neat.xp, neat.yp, "e", "c")$conc$yx Finally, to run multiple equatings in a single bootstrapping study, the arguments for each equating must be combined into a single object. Here, each element in neat.args is a named list of arguments for each equating. This object is then used in the bootstrap function, which carries out the bootstrapping.
R> neat.args <-list(i = list(type = "i"), + mt = list(type = "mean", method = "t"), + mc = list(type = "mean", method = "c"), + lt = list(type = "lin", method = "t"), + lc = list(type = "lin", method = "c"), + ef = list(type = "equip", method = "f", smooth = "log"), + ec = list(type = "equip", method = "c", smooth = "log"), + ct = list(type = "circ", method = "t"), + cc = list(type = "circ", method = "c", chainmidp = "lin")) R> bootout <-bootstrap(x = neat.xp, y = neat.yp, xn = xn, yn = yn, + reps = reps, crit = crit, args = neat.args) A plot method is available for visualizing output from the bootstrap function, as demonstrated below. Figures 10 through 13 contain the mean equated scores across replications for each method, the SE, bias, and RMSE. In Figure 10, the mean equated scores appear to be similar across much of the scale. Chained mean equating (the light orange line) consistently produces the highest mean equated scores. Mean equated scores for the remaining methods fall below those of chained mean and above those of identity equating (the black line). In Figure 11, standard errors tend to be highest for the equipercentile methods, especially chained equipercentile (the dark blue line), followed by the linear methods (green lines). SE are lowest for the circle-arc methods (purple and pink), especially in the tails of the score scale where the identity function has more of an influence. In Figure 12, bias is highest for chained mean equating, and is negative for the identity function; otherwise, bias for the remaining methods falls roughly between −0.5 and 0.5. Finally, in Figure 13, RMSE tends to be highest for chained mean and the linear and equipercentile methods. RMSE for Tucker mean and the circle-arc methods tended to fall at or below 0.5.

Summary
This paper presents some basic concepts and procedures for observed-score linking and equating of measurement scales. Linear and nonlinear functions are discussed, and various methods for applying them to nonequivalent groups are reviewed. Finally, the equate package is introduced, and its basic functionality is demonstrated using three data sets.
The equate package is designed to be a resource for teaching, research, and applying observedscore linking and equating procedures. A simple interface, via the equate function, can be used to control most of the necessary functionality, including data preparation, presmoothing, linking and equating, and managing output. Summary and plot methods facilitate the comparison of results. Future versions of the equate package will be extended to support additional procedures, for example, postsmoothing (e.g., Kolen 1991), nonlinear continuization (von Davier et al. 2004), additional asymptotic standard errors of equating, and new linking and equating functions as they are developed.

A. Equating methods for multiple anchors
The assumptions presented above for the Tucker, nominal weights, and frequency estimation methods are extended here to the total score distributions X P and Y Q and two or more anchor tests V 1 , V 2 , . . . , V m .

A.1. Tucker method
Tucker equating with multiple anchor tests involves a linear regression of total scores X or Y on a matrix of anchor scores V. For population P taking X, the regression model can be expressed in matrix notation as where x P is an n × 1 column vector of total scores on X for n individuals, the scalar δ X is the intercept, V P is an n × m matrix containing scores across m anchor tests, γ X is an m × 1 column vector of regression slopes, and e X P is the n × 1 column vector of residuals. Equation 47 is used here to derive the unobserved components of Equations 37 and 39 as functions of γ X and the observed means and standard deviations on X and V. Similar procedures are used to derive the unobserved values for Y in P .

Note that in Equation 47
there are no population subscripts on the regression coefficients δ X and γ X . The Tucker method assumes that these coefficients are the same in P as in Q. As a result, the total score means for X can be expressed as and where the means of e X P and e X Q are 0. The unobserved mean on X for Q can then be obtained through substitution and some rearranging as The unobserved mean on Y for P is obtained through a similar process as These values can then be used in Equations 33 and 34 to estimate µ X and µ Y .
Next, it is assumed that the residual variance for each regression model is equal across populations. With σ 2 E X representing the residual variance for both P and Q, assuming these variances are equal, the total score variances on X can be expressed as and where Σ V P and Σ V Q are the m × m variance-covariance matrices for the anchor tests, and γ X is the row vector of m anchor test slopes. The unobserved variance on X for Q can then be obtained through substitution and some rearranging as The unobserved variance on Y for P can be obtained similarly as Equations 54 and 55 can then be used within Equations 39 and 40 to obtain σ 2 X and σ 2 Y . These derivations for the Tucker method are simple extensions of equations presented for a single anchor (e.g., Kolen and Brennan 2014). Tucker equating with multiple anchors was first discussed by Gulliksen (1950).

A.2. Nominal weights method
Nominal weights mean equating with a single anchor (Babcock et al. 2012) can also be extended to accommodate multiple anchor tests. Test X is presented as an example. The column vector of total scores x P is first expressed as where x P represents an item response on X, µ x P is a mean item score or proportion correct taken across all items for an individual, and µ x P is a column vector of mean proportion correct for n individuals. Total scores on the anchor tests are expressed similarly as where v P represents an item response on V , µ v P is an n × m matrix of mean proportion corrects for n individuals on m anchor tests, and K V is a diagonal matrix containing the numbers of items on each anchor test.
Next, the mean proportion correct versions of x P and V P in Equations 56 and 57 are used to obtain the regression slopes in Equation 47. With x P and V P expressed in deviation scores, solving Equation 47 for γ X results in the following familiar equation: Substituting Equations 56 and 57 into 58 then results in The nominal weights method for a single anchor test assumes that, after accounting for the number of items on the total and anchor tests, mean performance on each test is equal in terms of proportion correct. This assumption is simply extended here to multiple anchors, where mean performance in terms of proportion correct is assumed to be equal for each individual across all anchors and the total test.
Postmultiplying µ x P by a row vector of 1s with length m, denoted by 1 , produces an n × m matrix where each column contains the mean proportion correct for each individual on X.
Under the assumption of equivalent mean proportion correct, µ x P 1 = µ v P , which leads to Postmultiplying by K V K −1 V then reduces Equation 60 to which can be simplified further by postmultiplying by a column vector of length m containing only 1 m , i.e., 1 1 m , as Equation 62, and the equivalent for Y , provide the nominal weights γ terms needed to estimate the unobserved means in (50) and (51).

A.3. Frequency estimation method
When using frequency estimation with multiple anchors, the multivariate score distribution of X and V 1 , V 2 , . . . , V m for the synthetic population is obtained as and the distribution of Y and V is obtained similarly as P(y, v 1 , v 2 , . . . , v m ) = w P P P (y, v 1 , v 2 , . . . , v m ) + w Q P Q (y, v 1 , v 2 , . . . , v m ).
Like in the single anchor case, the distributions P Q (x, v 1 , v 2 , . . . , v m ) and P Q (y, v 1 , v 2 , . . . , v m ) are unavailable, and must be estimated indirectly. This estimation is made possible by assuming that the conditional distributions of X given the anchors and Y given the anchors are population invariant. This results in and P P (y, v 1 , v 2 , . . . , v m ) = P Q (y, v 1 , v 2 , . . . , v m ) which are used in Equations 63 and 64 to obtain the multivariate synthetic distributions. Note that Equations 65 and 66 are simple generalizations of equations for equating with one anchor (e.g., Kolen and Brennan 2014;von Davier et al. 2004) and equating with two (Moses, Deng, and Zhang 2010) and three anchors (Angoff 1984).

C. Loglinear presmoothing
Polynomial loglinear modeling is a flexible procedure for smoothing distributions of various shapes to varying degrees. The structure of a distribution can either be maintained or ignored depending on the complexity of the model, where the degree of the polynomial term included determines the moment of the raw score distribution to be preserved. For example, a model with terms to the first, second, and third powers would create a smoothed distribution which matches the raw in mean, variance, and skewness. As shown below, the log of the expected relative frequency p for score point x is modeled as a function of a normalizing constant (the intercept β 0 ) and the observed-score value to the first, second, and third powers: Indicator variables may also be included to preserve specific moments for subsets of score points. In the next model, the mean and variance of a sub-distribution are preserved, in addition to the first three moments of the full distribution. When S = 1, score point x is included in this sub-distribution, and when S = 0, it is ignored: An acceptable degree of smoothing is typically achieved by comparing multiple models with different numbers of polynomial terms based on their fit to the data (Kolen and Brennan 2014). The loglinear function in equate is a wrapper for the glm function in the stats package. It can be used to fit and compare nested models up to specified maximum polynomial terms. For additional details, see the presmoothing help file.

D. Error in equating
In simulation and resampling studies, equating functions are typically compared based on both random and systematic error (or "differences"), where the first is estimated by the standard error of equating SE and the second by the bias. In the equate package, error is estimated in terms of the criterion equating function link Y (x) and estimate link Y r (x) for samples r = 1, 2, . . . , R. Systematic error is estimated as where is the average estimated equated score over R samples. The random error is estimated as Combining both systematic error and random error, the root mean squared error is estimated as RMSE = bias 2 + SE 2 .