Individual differences in paired comparison data

British Journal of Mathematical and Statistical Psychology (2001), 54, 265±277 © 2001 The British Psychological Society

Printed in Great Britain

265

Individual differences in paired comparison data Ulf Bo¨ckenholt* and Rung-Ching Tsai University of Illinois at Urbana-Champaign, USA

Thurstonian models provide a ¯exible framework for the analysis of multiple paired comparison judgments because they allow a wide range of hypotheses about the judgments’ mean and covariance structures to be tested. However, applications have been limited to a large extent by the computational intractability involved in ®tting this class of models. This paper demonstrates that the Monte Carlo EM algorithm facilitates maximum likelihood estimation of Thurstonian paired comparison models even when the number of items is large. A paired comparison study is presented in detail to illustrate the estimation approach.

1. Introduction Probabilistic models of attitudes and preferences originated with Thurstone’s (1927) seminal work. By invoking the psychophysical concept of a sensory continuum, Thurstone explained the observation that choices by the same person may vary under seemingly identical conditions. He argued that choice options should be represented along this sensory continuum by random variables that describe the options’ effects on a person’s sensory apparatus. The comparison process is then reduced to selecting the option with the greatest realized value. This conceptualization proved to be particularly attractive in the modelling of paired comparison data where each judge is asked to compare two items at a time. Thurstone assumed, mainly for estimation reasons, that choice options are judged independently of each other; this led to the famous Case V model. However, despite the simplicity of this proposal, more general forms of the dependence structure are worthy of investigation. In particular, when a judge is asked to compare all possible pairs of choice options, it seems likely that the notion of independent evaluations as implied by Thurstone’s Case V model is too restrictive and that more general formulations are desirable (Takane, 1987). Initial empirical studies of this issue by BoÈckenholt & Dillon (1997) were restricted to three items for computational reasons. Here we investigate the usefulness of Thurstonian models with an unrestricted covariance matrix and show that the Monte Carlo expectation±maximization (MCEM) algorithm (Wei & Tanner, 1990) provides an attractive estimation approach even when the number of items to be compared is large. When judges compare several pairs of items, it is important to distinguish sources of within- and between-person variability. For this purpose, Bock (1958) introduced a `threecomponent model’ that includes two random effects representing individual preference This research was partially supported by NSF grant SBR97-30197. * Requests for reprints should be addressed to Ulf BoÈckenholt, Department of Psychology, University of Illinois, Champaign, Il 61820, USA, e-mail: ([email protected]).

266

U. Bo¨ckenholt and R.-C. Tsai

differences for the options and a response error speci®c to each pair of options. Takane (1987) extended Bock’s (1958) model by allowing for different covariance structures for the individual difference components. These extensions are particularly useful when the investigation of similarity relationships among the choice options is of interest. However, empirical applications of Takane’s (1987) approach have not been reported in the literature because of computational dif®culties in estimating the model parameters. High-dimensional integration is required when analysing paired comparison data obtained in a multiplejudgment experiment. Fortunately, recent advances in the implementation of both the MCEM algorithm (Tanner, 1996) and the Gibbs sampler for the truncated multivariate normal distribution (Hajivassiliou, McFadden & Ruud, 1996) have overcome these computational problems. Using these methods, we can avoid numerical integration whose accuracy is known to diminish rapidly with the dimensionality of the integration. In contrast, Monte Carlo simulations are more stable in high dimensions (Meng & Schilling, 1996), as recently demonstrated by Yao & BoÈckenholt (1999) for the estimation of Thurstonian ranking models, albeit in a Bayesian context. The next section presents a general formulation of Thurstonian paired comparison models. Section 3 discusses estimation by MCEM and model testing issues. The remaining part of the paper focuses on simulation studies to investigate small-sample properties of the estimation approach and an application designed to test the presence of similarity relationships in paired comparison data. 2. The Thurstonian paired comparison model In a complete paired comparison experiment a judge makes r(r 2 1)/2 comparisons of r items. However, it is not necessary to assume that each individual compares all items. Instead, item pairs may be selected according to an incomplete paired comparison design to reduce the number of judgments made by each individual (Bock & Jones, 1968). This approach is particularly desirable when the number of items to be compared exceeds nine or ten. To simplify notation we discuss the case of a complete paired comparison design. Let W 5 (w1 , . . . , wn ) represent the observed individual response vectors, where wi 5 (wi12, wi13 , . . . , wi(r2 1)r ) are the binary paired comparison outcomes of person i. According to Thurstone (1927), each binary response is a result of a latent difference judgment. Thus, in a comparison between items j and k, the latent judgment yi j k is represented as a difference between two continuous random variables, yi j k 5 yi j (k) 2 yi k ( j ) , where yi j (k) 5

mj 1

ni j 1

ei j (k) ,

yi k ( j ) 5

mk 1

ni k 1

ei k ( j ) .

(1)

The ®xed effects for items j and k are denoted by mj and mk , and the corresponding personspeci®c random effects are given by ni j and ni k . The terms ei j (k) and ei k ( j ) capture random variations in the judgmental process and are assumed to be identically and independently distributed for all item pairs. Brackets around j and k indicate that these subscripts serve merely as labels. The choice between the two items is determined by the sign of the difference yi j k . Since no assumptions are made about the origin of the ®xed or random effects, (1) is applicable both in preference and in discrimination studies. Thus, an item may be chosen because it obtained a

Individual differences in paired comparison data

267

higher value on some hedonic or on some sensory continuum. In the former case, it seems likely that one of the primary mental activities involves the reduction of multiple option attributes to a single one through the use of trade-offs. The notion of maximizing a unidimensional criterion is consistent with alternative representations of preferential choice that involve a person-speci®c ideal point (Andrich, 1989; Luo, 1998). For example, the random effect vector ni may be a result of a comparison between the ideal point of person i and the options’ values on some continuum. We will return to this point in the discussion section. Let ei 5 (ei12 , ei13, . . . , ei (r2 1)r) , with ei j k 5 ei j (k) 2 ei k ( j ) , and ni 5 (ni1 , . . . , nir ) be normally distributed with ei ~ N(0, j2 I), ni ~ N(0, Sn ), and ei \ni . Then the joint probability of 2r paired comparisons yi 5 ( yi12 , yi13, . . . , yi(r2 1)r ) is given by an 2r -dimensional normal distribution function F(2r) with mean vector A*m 5 (m12 , m13, . . . , m(r2 1)r ) and covariance matrix V 5 A*Sn A*9 1 j2 I, where mj k 5 mj 2 mk and A* is an 2r 3 r paired comparison design matrix. Each column of A* corresponds to one of the items, and each row to one of the r 3, 2 paired comparisons. For example, when r 5 A 5

1

2 1

1

0

2 1 .

0

1

2 1

0

The marginal probability that item j is selected in a comparison of items j and k is P(wj k 5

1 | mj , mk ) 5

F

j2 1

mj 2

mk

j 2nj 2

2j nj k 1

j 2nk

,

(2)

where F denotes the standard normal cumulative distribution function (Thurstone, 1927). Thurstone’s (1927) well-known Case V is obtained when the n effects are normally distributed with covariance matrix Sn 5 j2n I. In this case, the marginal paired comparison probabilities reduce to P(wj k 5

1 | mj , mk ) 5

F

mj 2 j 1 2

mk 2j2n

,

(3)

(see Bock & Jones, 1968). This parsimonious representation is useful when items are evaluated independently of each other. It is well known that because of the binary nature of the data only the correlation matrix corresponding to V, denoted by Corr(V), can be identi®ed. From this result it follows that the pair-speci®c variance j 2 cannot be identi®ed, and for convenience we set j2 5 1. In addition, because of the difference structure of the judgments only b 5 Cm, u 5 Cn and Su 5 CSn C9 are identi®able, where C is an (r 2 1) 3 r contrast matrix such as C 5 [Ir2 1 2 1] (Tsai, 2000). As a result, when reporting the item-speci®c scale values, we arbitrarily ®x the scale of the mean values by setting mr 5 0. An identi®ed structure of Sn is obtained by constraining the diagonal elements of Sn equal to 1 and one additional off-diagonal element to zero (Dansie, 1986). In the reported application we apply Dansie’s (1986) parameterization, instead of reported Sn which is dif®cult to interpret. However, we note that although the selection of the off-diagonal cell is arbitrary as long as Sn remains positive de®nite, it affects the values of

U. Bo¨ckenholt and R.-C. Tsai

268

the remaining off-diagonal elements. The relationship between the elements of two (2) (1) (2) covariance matrices S(1) n and Sn with different ®xed values j nj k and j nj k is linear and (1) (2) (1) È ckenholt, given by j ni j 5 (1 2 f) 1 f j ni j , where f 5 (1 2 j njk )/(1 2 j nj(2) k ) (Yao & Bo 1999). 3. Model estimation and testing 3.1. Estimation The simplicity and stability of the EM algorithm makes it an attractive choice for estimating the parameters of the Thurstonian paired comparison model. Because the high dimensionality of y prevents direct numerical calculations (except when r 5 3), we use Gibbs sampling in the E-step of the algorithm to obtain the conditional mean vector and covariance matrix of y given W (Tanner, 1996). McCulloch (1994) and Chan & Kuk (1997) apply the same approach for the estimation of mixed effects models with binary data. The implementation of MCEM is described below in detail. Let Y 5 (y1 , . . . , yn ) represent the unobserved latent response utilities, where yi 5 ( yi12 , yi13, . . . , yi(r 2 1)r )9 are the latent difference judgments of person i. Denote v º (b, Su ) and A as the full-column-rank version of A*, by omitting the last column of A*. We treat {Y, u, W} as the complete data and the observed data W as the incomplete data. à Sà u ) are Given the complete data {Y, u, W}, the maximum likelihood estimators of v º (b, bà 5

1 n

n

(A9 A)2 A9 (yi 2

Aui ),

(4)

i5 1

and SÃ u 5

1 n

n

ui u9i ,

(5)

i5 1

where (A9 A)2 is the generalized inverse of A9 A. M-step. By replacing yi , ui and ui u9i in (4) and (5) with their conditional expectations evaluated at the current estimates v(k) º (b(k) , S(k) u ) from the E-step, we update our parameter estimates such that b(k1

1)

5

1 n

n

(A9 A)2 A9 [E(yi | W, v(k) ) 2

A(E(ui | W, v(k) ))]

(6)

i5 1

and 1 S(k u

1)

5

1 n

n

E(ui u9i | W, v(k) ).

(7)

i5 1

E-step. Unfortunately, the conditional expectations cannot be expressed in closed form for a given W and parameter estimate v. Instead, they are approximated by simulations in the MCEM approach (Wei & Tanner, 1990; McCulloch, 1994). However, both E(ui | W, v) and

Individual differences in paired comparison data

269

E(ui u9i | W, v) can be expressed entirely in terms of Cov(yi | W, v) and E(yi | W, v) such that E(ui | W, v(k) ) 5

E[E(ui | yi )| W, v(k) ] 5

(k) 9 2 1 S(k) u A V [E(yi | W, v ) 2

Ab(k) ]

(8)

and E(ui u9i | W, v(k) ) 5 5

E[E(ui u9i | yi )| W, v(k) ] (k) 2 1 9 2 1 S(k) u A V DV ASu 1

(k) 9 2 1 S(k) u A V ASu ,

S(k) u 2

(9)

where D 5 Cov(yi | W, v(k) ) 1 [E(yi | W, v(k) ) 2 Ab(k) ][E(yi | W, v(k) ) 2 Ab(k) ]9 . Thus, it is suf®cient to simulate only y to obtain approximations of the conditional expectations required for the M-step. Because of the high dimensionality of y, we use Gibbs sampling (Geman & Geman, 1984) to simulate draws from the distribution of f (y| W, v) (Tanner, 1996). In other words, we draw yi from yi | W, v ~ TMVN(Ab, ASu A9 1 I), where TMVN denotes a truncated multivariate normal distribution. To incorporate the information of W, a random draw yi from the TMVN distribution can be obtained by drawing from N(2r ) (Ab, ASu A9 1 I) and then discarding those yi with values inconsistent in sign with wi . But as r gets large, this method becomes inef®cient since a large number of draws are rejected before the sign condition is satis®ed. Fortunately, a more ef®cient process was introduced by Hajivassiliou (1993) that generate samples from a TMVN distribution by the Gibbs sampler. Instead of drawing directly from a TMVN distribution, the TMVN Gibbs sampler cycles successively through 2r fully conditional truncated univariate normal distributions in which the range of each yi j k of yi is determined by wi j k such that (T 1 1) wi j k 5 1 implies yi j k > 0. The sequence of draws {y(T) , . . .} converges to f (yi | W, v) as i , yi a limiting distribution. Hajivassiliou et al. (1996) provide a detailed discussion of this approach.1 Vijverberg (1997) reviews several comparative studies which demonstrate that this Gibbs sampler simulator is among the best of the currently available simulation algorithms for the (truncated) multivariate normal distribution. Finally, the conditional expectation and covariance of y needed for the M-step can be approximated by averages of their sample values. To be precise, we have E(yi | W, v) 5

1 T1 S m t5

yÅ i | W, v <

m (t) T (yi |

W, v)

and Cov(yi | W, v) <

1 T1 S m t5

m (t) T (yi

2

yÅ i | W, v)(y(t) i 2

yÅ i | W, v)9 .

From the initial starting values v(0) º (b(0) , S(0) u ), the MCEM algorithm substitutes the expectation and covariance of y obtained from Gibbs sampling into (6) and (7) to update the estimates v in the M-step and continues to carry out the EM iterations until convergence is attained to obtain the ®nal estimate of và º (m, à Sà u ). Only minor modi®cations are required for 2 the estimation of the Case V model Sn 5 jn I, which are outlined in the Appendix. 1

FORTRAN and GAUSS computer code for the Gibbs sampler are available via anonymous ftp at the Internet site: http://econ.lse.ac.uk/staff/vassilis.

270

U. Bo¨ckenholt and R.-C. Tsai

As noted by Tanner (1996), two important considerations are the speci®cation of the number of Monte Carlo draws, m, and the monitoring of convergence. In the reported applications, we found that m 5 100 led to a smooth convergence behaviour of the algorithm. We assessed convergence by plotting v(t) versus the iteration number (t). After a certain number of iterations, the plot will reveal random ¯uctuation about the estimated population parameters. This approach is illustrated in the reported application. 3.2. Model testing Large-sample tests of ®t are available based on the likelihood ratio (LR) x2 statistic (G2 ) and/ or Pearson’s goodness-of-®t test statistic (P2 ). Asymptotically, if a paired comparison model provides an adequate description of the data, then both statistics follow a x2 distribution with 2r(r2 1)/2 2 s 2 1 degrees of freedom, where s refers to the number of parameters to be estimated. The LR test is most useful when the number of items to be compared is small. Otherwise, only a small subset of the possible paired comparison patterns may be observed. In this case, it is doubtful that the test statistics will follow approximately a x2 distribution. However, useful information about a model ®t may be obtained by inspecting standardized differences between the observed and ®tted paired comparison probabilities for subsets (e.g., all triplets) of the items. Although these residuals are not independent, a careful inspection of their direction and size is useful in identifying systematic mis®t. Nested paired comparison models can be compared by computing twice the difference between their log-likelihood values. This difference is asymptotically distributed as a x2 statistic with the degrees of freedom equal to the difference between the number of parameters in the unrestricted and restricted models. Values of the log-likelihood functions are computed based on the GHK smooth recursive stimulator (Geweke, 1989; Hajivassiliou, 1993; Keane, 1994). 4. Simulation studies In this section, we present the results of several small simulation studies to evaluate the performance of the MCEM algorithm for estimating the paired comparison parameters. We ®rst compare estimates obtained by numerical integration and the MCEM algorithm for r 5 3 and n 5 100 judges. The mean vector was speci®ed as m 5 [2 1, 2 0.5, 0], with covariance matrices for the ®rst and second study given by Sn 5 I and Sn 5

1

0

0.6

0

1

0.6

0.6

0.6 , 1

(10)

respectively. The within-pair variance j2 was set equal to 2. For both sets of population parameters, 100 frequency tables were generated. Table 1 contains the identi®ed parameters as well as the means of the parameter estimates and their respective standard deviations estimated by numerical integration and MCEM. The MCEM estimates were obtained by setting m 5 100 and the number of iterations equal to 200. Because convergence was obtained in less than 50 iterations with random starting values, the MCEM estimates were computed as averages of the last 50 iterations. We note that the

Individual differences in paired comparison data

271

Table 1. Two sets of identi®ed population parameter values and maximum likelihood estimates obtained by numerical integration and the MCEM algorithm (r 5 3, n 5 100) Study 1 Parameters m1 2 m3 m2 2 m3 j2u1 ju12 j2u2

Study 2

True

Num. Int

MCEM

True

Num. Int.

MCEM

2 0.71

2 0.72 (0.21)

2 0.72 (0.21)

2 0.71

2 0.71 (0.17)

2 0.72 (0.17)

1.00 0.50 1.00

1.22 (0.94) 0.64 (0.56) 1.29 (0.97)

1.20 (0.91) 0.63 (0.55) 1.28 (0.93)

0.40 2 0.10 0.40

0.46 (0.44) 2 0.06 (0.29) 0.48 (0.42)

2 0.04 (0.26)

2 0.35

2 0.37 (0.20)

2 0.37 (0.19)

2 0.35

2 0.38 (0.12)

2 0.38 (0.12)

0.48 (0.41) 0.50 (0.40)

Note: The numbers in parentheses are the standard deviations of the maximum likelihood estimates obtained in 100 replications.

estimates obtained by numerical integration and MCEM are very close. As expected, the mean estimates exhibit much less bias than the estimates of the variance components. The standard deviations of the variances and covariances are large, indicating that with 23 data cells, it is dif®cult to obtain precise estimates of the random effects structure. Numerical integration is cumbersome with current methods for r > 3. In Table 2 we report the results of two additional simulation studies with r 5 4 for the two sample sizes n 5 100 and n 5 300. As before, convergence was rapid with m 5 100, and Table 2 reports the averages of the last 50 iterations. We note that the positive bias of the random effects estimates is much reduced when n 5 300. Equally important, the standard deviations of the parameter estimates decrease substantially in comparison to r 5 3, which demonstrates that the larger number of cells facilitates a more precise estimation of the random effects. Overall, these results indicate that the MCEM algorithm provides a promising approach for the estimation of Thurstonian paired comparison models. 5. Similarity relationships in paired comparison data Adopting the design of the Rumelhart & Greeno (1971) study, Kroeger (1992) presented 164 students (68 males and 96 females) with pairs taken from a set of three actors, three political Table 2. Identi®ed population parameter values and maximum likelihood estimates obtained by the MCEM algorithm for (r 5 4, n 5 100) and (r 5 4, n 5 300) Parameters m1 2 m4 m2 2 m4 m3 2 m4 j2u1 ju12 j2u2 ju13 ju23 j2u3

True

n5

100

n5

300

2 1.20

2 1.22 (0.20)

2 1.22 (0.09)

2 0.40

2 0.40 (0.20)

2 0.42 (0.11)

2 0.80

0.80 0.10 1.40 0.70 0.70 2.00

2 0.81 (0.17)

0.87 (0.45) 0.18 (0.32) 1.54 (0.63) 0.75 (0.40) 0.82 (0.45) 2.16 (0.75)

2 0.78 (0.10)

0.83 (0.31) 0.10 (0.21) 1.42 (0.36) 0.69 (0.29) 0.71 (0.28) 2.02 (0.42)

Note: The numbers in parentheses are the standard deviations of the maximum likelihood estimates obtained in 100 replications.

272

U. Bo¨ckenholt and R.-C. Tsai

Individual differences in paired comparison data Table 3. Observed and expected frequencies for three athletes Comparisons

273

a

Females

Males

yj k

yj h

yk h

obs.

pred.

obs.

pred.

0 0 0 0 1 1 1 1

0 0 1 1 0 0 1 1

0 1 0 1 0 1 0 1

19 23 2 16 12 2 7 15

19.5 20.7 1.9 19.5 13.0 1.9 8.2 11.2

21 17 0 5 19 0 6 0

22.1 16.0 1.1 3.8 18.1 1.1 4.3 1.5

a

j, Bonnie Blair; k, Jackie Joyner-Kersee; h, Jennifer Capriati.

Table 4. Observed and expected frequencies for three actorsa Comparisons

Females

Males

yj k

yj h

yk h

obs.

pred.

obs.

pred.

0 0 0 0 1 1 1 1

0 0 1 1 0 0 1 1

0 1 0 1 0 1 0 1

21 18 1 17 11 1 18 9

17.9 19.2 3.2 15.1 12.5 3.2 11.8 13.3

7 15 0 7 20 0 17 2

8.6 10.0 3.0 7.6 18.8 3.0 12.5 4.4

a

j, Robert De Niro; k, Dustin Hoffman; h, Jack Nicholson.

®gures and three athletes. For each of the 36 pairs students were asked to indicate which person they would prefer to interview. Because of the small sample size the following analyses are restricted to the 15 paired comparison outcomes based on the actors (Robert De Niro, Dustin Hoffman and Jack Nicholson) and athletes (Bonnie Blair, Jackie Joyner-Kersee and Jennifer Capriati). The purpose of the following analyses is to investigate whether the similarity relationships among the choice options (male actors versus female athletes) is re¯ected by the covariance matrix of the person-speci®c random effects. To test this hypothesis we ®tted two models: model (a) with independently distributed random effects and model (b) with an unconstrained covariance matrix. The four panels of Figure 1 illustrate the MCEM convergence behaviour when both models are ®tted to the paired comparison data of the male respondents and m 5 100. Figures 1(a) and 1(b) contain plots of the Case V parameter deviations from the ®nal parameter estimates separately for the mean and variance (j 2n ) effects. Figures 1(c) and 1(d) contain the

Figure 1. Parameter changes for paired comparison data set of male respondents. (a) Main effects versus EM iterations, Case V; (b) variance effect versus EM iterations, Case V; (c) main effects versus EM iterations, general case; (d) (co)variance effects versus EM iterations, general case.

U. Bo¨ckenholt and R.-C. Tsai

274

Table 5. Fit statistics of binary paired comparison data Model

No. parameters

(a) (b)

6 20

females

males

2 Log-lik.

av.P2

2 Log-lik.

av.P2

695.6 670.0

8.9 4.6

447.8 418.8

9.2 4.5

corresponding results for the unrestricted covariance matrix. Although the Case V estimates converge in a slightly smaller number of iterations than the general case estimates, Fig. 1 demonstrates that the deviations display a smooth convergence pattern. Figure 1(d) is illustrative of the larger variation and autocorrelation in the variance components, but considering the size of the data set (n 5 68), these variations are to be expected. The MCEM results for the data set of the female respondents are similar but exhibit less variability in the random effects parameters. Because the data table is too sparse for an overall goodness-of-®t test to be used, we report the log-likelihood function values and determine the ®t of a model by computing residual statistics based on collapsed cells of the data matrix. In particular, we computed the expected frequencies for all 63 5 20 triplets that can be formed from the six items. The observed and expected frequencies under a model with a general covariance matrix for the actor and athlete triplets are reported in Tables 3 and 4, respectively. Note that the number of intransitive cycles (i N j N k N i and i N k N j N i) is quite small in the data set. Only four female respondents exhibit intransitive judgments when comparing the athletes, and two female respondents are intransitive in their comparisons of the actors. Table 5 summarizes the ®t statistics obtained under the two speci®cations (a) and (b). Using an LR test to compare both models, we note that the model with an unrestricted covariance matrix ®ts signi®cantly better for both the female and male respondents (G2 5 51.2 (df 5 14) for females; G2 5 58 (df 5 14) for males). To assess the overall ®t of the models, we computed Pearson’s x2 statistics (P2 ) based on the observed and expected frequencies of all possible triplets. The fourth and sixth columns of Table 5 contain the average values of the statistics obtained for the 20 subtables. With average values of 4.6 and 4.5 for the data of the female and males respondents, respectively, the ®t of model (b) appears Table 6. Estimated covariance matrix and mean values RN DH JN BB JK JC

RN

DH

JN

BB

JK

JC

mà F

mà M

1.00 0.56 0.54 2 0.19 2 0.12 0*

0.64 1.00 0.62 0.07 0.11 0.28

0.89 0.61 1.00 2 0.16 0.02 2 0.23

0.26 0.49 0.31 1.00 0.36 0.22

0.39 0.28 0.42 0.29 1.00 0.14

0* 0.38 0.10 0.63 0.29 1.00

0.77 0.98 0.91 2 0.26 0.20 0*

0.26 0.08 0.43 2 1.02 2 0.57 0*

Estimated values below and above the main diagonal are obtained from the female and male respondents, respectively. RN, Robert De Niro; DH, Dustin Hoffman; JN, Jack Nicholson; BB, Bonnie Blair; JK, Jackie Joyner-Kersee; JC, Jennifer Capriati. 0* 5 ®xed; F 5 females; M 5 males.

Individual differences in paired comparison data

275

to be satisfactory. The worst-®tting set of triplets is obtained for the three actors (P2 5 8.7 for females, P2 5 11.9 for males). As can be seen from Table 4, the paired comparison model does not capture well the fact that the male respondents are perfectly consistent in their comparison of the actors. The estimated covariance parameters and mean values are reported in Table 6. The lower diagonal of the covariance matrix contains the results obtained from the female respondents. We note that associations within actors and athletes are higher than between the two sets, which supports our initial expectation about the similarity relationships among the choice options. A similar but less pronounced pattern is obtained for the male respondents, whose covariances are listed in the upper triangular part of the matrix. Although both male and female respondents prefer to interview actors rather than athletes, there are strong gender differences within both sets. For example, while female students prefer Dustin Hoffman, males prefer Jack Nicholson. 6. Conclusion Over the years the method of paired comparisons has proven to be an important tool for preference and value measurements in the social sciences (David, 1988). By asking judges to compare two items with respect to some criterion, this technique typically leads to sharper discriminations than obtainable by other methods. Thurstonian paired comparison models provide an attractive framework for the analysis of paired comparison judgments because they facilitate the testing of a wide range of hypotheses about the mean and covariance structure of the data. This paper has shown that the MCEM algorithm facilitates the estimation of Thurstonian paired comparison models. Both the simulation studies and the application demonstrate that the MCEM algorithm is easy to implement and converges quickly even when the paired comparison data are sparse. In our presentation of the algorithm we have focused on the identity and unrestricted covariance structures of the random effects. In some applications it may be desirable to specify models between these two covariance forms to obtain a parsimonious description of the random effects. Such an approach would also be useful in view of the small-sample bias found in the simulation studies. Perhaps factor-analytic models are most promising in this context because they facilitate the identi®cation of individual difference and similarity effects (Brady, 1989; Takane, 1987). In this case, yi j (k) can be written as S

yi j (k) 5

mj 1

zis lsj 1

di j 1

ei j (k) ,

(11)

s5 1

where zis denotes the weight person i places on attribute s (s 5 1, . . . , S), lsj the factor coef®cient of attribute s on item j, and di j the unique characteristics of item j speci®c to person i. Alternatively, we may postulate an ideal point model by assuming that items can be represented as points in an S-dimensional perceptual space and person i’s preference judgment for item j is determined by the squared Euclidean distance between person i’s ideal point and item j. As a result, yi j (k) is written as yi j (k) 5

mj 2

1 2

S

(zis 2 s5 1

lsj )2 1

di j 1

ei j (k) ,

(12)

U. Bo¨ckenholt and R.-C. Tsai

276

where zis and lsj denote the sth coordinates of the ideal point for person i and the locations of item j in an S dimensional space, respectively, and di j denotes the unique characteristic of item j speci®c to person i. It is important to note that although both representations are quite different, their difference structures in the form of yi j k are identical. Thus, for both (11) and (12), we obtain yi j k 5

(mj 2

mk ) 1

S

zis (lsj 2

lsk ) 1

(di j 2

di k ) 1

ei j k .

(13)

s5 1

This result demonstrates that the Thurstonian approach may be equally useful for the analysis of individual difference effects in preference and in discrimination data since only (the signs of the) difference judgments are observable. Based on our results reported here, we expect that the MCEM algorithm is well suited to estimating this class of models and, more generally, that it will become one of the standard estimation approaches for Thurstonian paired comparison models. References Andrich, D. (1989). A probabilistic item response theory model for unfolding preference data. Applied Psychological Measurement, 13, 193±216. Bock, R. D. (1958). Remarks on the test of signi®cance for the method of paired comparisons. Psychometrika, 23, 323±334. Bock, R. D., & Jones, L. V. (1968). The measurement and prediction of judgment and choice. San Francisco: Holden-Day. BoÈckenholt, U., & Dillon, W. R. (1997). Modeling within-subject dependencies in ordinal paired comparison data. Psychometrika, 62, 414±434. Brady, H. E. (1989). Factor and ideal point analysis for interpersonally incomparable data. Psychometrika, 54, 181±202. Chan, S. K., & Kuk, A. Y. C. (1997). Maximum likelihood estimation for probit±linear mixed models with correlated random effects. Biometrics, 53, 86±97. Dansie, B. R. (1986). Normal order statistics as permutation probability models. Applied Statistics, 35, 269±275. David, H. A. (1988). The method of paired comparisons. London: Grif®n. Geman, D., & Geman, D. (1984). Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6, 721±741. Geweke, J. (1989). Bayesian inference in econometric models using Monte Carlo integration. Econometrica, 57, 1317±1340. Hajivassiliou, V. A. (1993). Simulation estimation methods for limited dependent variable models. In G. S. Maddala, C. R. Rao & H. D. Vinod (Eds.), Handbook of Statistics, Vol. 11. Amsterdam: Elsevier. Hajivassiliou, V. A., McFadden, D., & Ruud, P. (1996). Simulation of multivariate normal rectangle probabilities and their derivatives. Theoretical and computational results. Journal of Econometrics, 72, 85±134. Keane, M. (1994). A computationally practical simulation estimator for panel data. Econometrica, 62, 95±116. Kroeger, K. (1992). {Modeling intransitivity through log-bilinear models.} Unpublished M.A. thesis. Department of Psychology, University of Illinois, Urbana-Champaign. Luo, G. (1998). A general formulation of unidimensional unfolding and pairwise preference models: Making explicit the latitude of acceptance. Journal of Mathematical Psychology, 42, 400±417. McCulloch, C. E. (1994). Maximum likelihood variance components estimation for binary data. Journal of the American Statistical Association, 89, 330±335.

Individual differences in paired comparison data

277

Meng, X.-L., & Schilling, S. (1996). Fitting full-information item factor models and an empirical investigation of bridge sampling. Journal of the American Statistical Association, 91, 1254± 1267. Rumelhart, D. L., & Greeno, J. G. (1971). Similarity between stimuli. An experimental test of the Luce and Restle choice models. Journal of Mathematical Psychology, 8, 370±381. Takane, Y. (1987). Analysis of covariance structures and probabilistic binary choice data. Cognition and Communication, 20, 45±62. Tanner, M. A. (1996). Tools for statistical inference. New York: Springer-Verlag. Thurstone, L. L. (1927). A law of comparative judgment. Psychological Review, 34, 273±286. Tsai, R. (2000). Remarks on the identi®ability of Thurstonian ranking models: Case V, Case III, or Neither? Psychometrika, 65, 233±240. Vijverberg, W. P. M. (1997). Monte Carlo evaluation of multivariate normal probabilities. Journal of Econometrics, 76, 281±307. Wei, G. C. G., & Tanner, M. A. (1990). A Monte Carlo implementation of the EM algorithm and the poor man’s data augmentation algorithms. Journal of the American Statistical Association, 85, 699±704. Yao, G., & BoÈckenholt, U. (1999). Bayesian estimation of Thurstonian ranking models based on the Gibbs sampler. British Journal of Mathematical and Statistical Psychology, 52, 79±92. Received 17 December 1999; revised version received 11 August 2000

Appendix: EM steps for Case V covariance structure Let v º (A m, j2n A A 9 1

(m, j 2n ).

Using Gibbs sampling, we obtain draws of yi | W, v(k) from TMVN I).

E-step: E(ni | W, v(k) ) 5

E[E(ni | yi )| W, v(k) ] 5

(k)

j 2n A 9 V2 1 [E(yi | W, v(k) ) 2

A m(k) ]

and E(ni n9i | W, v(k) ) 5

E[E(ni n9i | yi )| W, v(k) ] (k)

j4n A 9 V2 1 DV2 1 A 1 5

where D 5

Cov(yi | W, v(k) ) 1

(E[yi | W, v(k) ] 2

(k)

j 2n (I 2

(k)

j 2n A 9 V2 1 A ),

A m(k) )(E[yi | W, v(k) ] 2

A m(k) )9 .

M-step: m (k1

1)

5

1 n

n

(A 9 A )2 A 9 [E(yi | W, v(k) ) 2

A (E(ni | W, v(k) ))]

i5 1

and j 2n

(k1 1)

5

1 9 1 diag nr

n

E(ni n9i | W, v(k) ) ,

i5 1

where diag(M) is a vector containing the diagonal elements of matrix M.