## How to Calculate Expected A Posteriori (EAP) Scores for Unidimensional and Multidimensional Latent Trait Models

### 1. Introduction

Sometimes when one performs a latent trait analysis the goal is only to characterize the latent factor structure of a set of binary or ordered-category variables. That is, one mainly wishes to know the loadings (correlations) of latent response variables on the latent trait or latent traits.

Definitions: a latent response variable Y is the unobserved (latent) continuous precursor variable associated with an observed (manifest) binary or ordered-category variable X. If there are k manifest variables X1, X2, ..., Xk being anaylyzed, there are k associated latent response variables, Y1, Y2, ... Yk. The latent response variables are distinct from the latent trait(s); latent traits are to latent response variables what factors are to continuous measures in a standard factor analysis model.
In some applications the locations of the discretizing thresholds may also be of interest.

Other times one may also wish to score cases relative to the latent trait(s)--this is the equivalent of calculating factor scores with continuous-data factor analysis. With a multidimensional latent trait model, this can also be understood as estimating the location of an object in the space defined by the latent traits.

Such scores can be obtained with different methods. One method both practical and theoretically appealing is Expected A Posteriori (EAP) estimation. Bock and Aitkin (1981; pp. 452-454) give a general but brief description of the method. Uebersax (1993, Section 2.5) and Clogg (1988) also give a simple formulas for EAP estimation that apply to unidimensional discrete latent trait (located latent class) models. However, an introductory discussion of the subject, such as might suit one who is not already experienced with latent trait models, appears lacking in the literature; we will attempt to remedy this here.

The present plan: Section 2 briefly reviews the premise of EAP estimation; Section 3 describes its use to estimate latent trait scores for a unidimensional latent trait model; and Section 4 extend this to multidimensional latent traits.

Attention is restricted to fully Gaussian latent trait models--that is, here (1) the latent trait is assumed to have a normal (or, for multiple latent traits, a multivariate normal) distribution; and (2) response functions are modeled as normal-ogives (normal cumulative distribution functions; this assumption follows from the assumption of normally distributed measurement error). The methods described here, though, are applicable to other types of latent trait models, such as those based on logistic ogive assumptions.

We consider here only dichotomous data, but the principles easily generalize to EAP estimation of latent trait scores for ordered-category data (for an example of this in the unidimensional case, see Uebersax, 1993).

The methods here apply regardless of how the latent trait model is estimated. The latent trait model can be estimated by full-information maximum likelihood (FIML) estimation (e.g., using the TESTFACT program, or, for unidimensional models, some of the programs available on these web pages), but that is not necessary. As already described, one can also estimate a fully Gaussian latent trait model by the heuristic method of factoring tetrachoric correlations; EAP estimation can be used in this case as well.

### 2. EAP Estimation Defined

EAP or Expected A Posteriori estimation derives from Bayesian statistical principles. The concept is not remote or difficult to understand. Rather, EAP estimation is a very common-sense and easily accessible method. We will present its derivation in Bayesian terms. However, we will see that the result is rather obvious.

The term "a posteriori" derives from the Bayesian concept of a posterior probability. In this context it refers to a posterior probability distribution of latent trait scores--specifically, the predicted distribution of scores for a given case given (a) the response pattern of that case, and (b) the estimated model parameters. The term "expected" derives from the concept of an expected value. Thus an "expected a posteriori" estimate refers to the expected value of the posterior probability distribution of latent trait scores for a given case.

### 3. EAP Estimation with a Unidimensional Latent Trait

We begin with a simple formula that applies Bayes' rule:

```
Pr(T) Pr(X|T, pi)
Pr(T|X, pi) =  ---------------------         (1)
SUM Pr(T) Pr(X|T, pi)

where:  T     is the score of a case relative to the latent trait;
Pr(T) is the prior probability that the case has latent
trait level T; this corresponds to our "belief"
about the distribution of T *before* we consider
any information about the particular case (e.g.,
the response pattern of the case).
X     (a vector) is the pattern of responses of the case
across the manifest variables;
pi    is the vector of estimated model parameters (i.e.,
latent correlation and threshold parameters);

```
and summation takes place across all values of T. In theory, this would be shown by an integral sign, since T is continuous; however, in practice, we will inevitably perform integrations approximately--by summing--so it is not misleading to show SUM here

The meaning of Equation (1) is as follows. The term on the left side is the probability that a case has latent trait level T, given the case's particular observed pattern of responses, X, and the parameter estimates for the model, pi. There is a value of this probability associated with every possible value of T (often one considers every possible value of T from -3 to 3, assuming that the latent trait has a standard deviation of 1 and a mean of 0).

The set of all values for T and the associated probabilities Pr(T|X, pi) define the posterior probability distribution of T for a case. Such a distribution is illustrated by Figure 1.

```
Pr(T|X, pi))
|
|
|                   ..
|                 .    .
|               .        .
|              .          .
|             .            .
|            .              .
|          .                  .
|        .                      .
|  . . .                          . . . . . . . . . . . .
+----+-----+-----+-----+-----+-----+-----+-----+-----+-----+----+----+
-3         -2          -1           0           1           2         3
Latent Trait Level (T)

Figure 1 (draft).  Distribution of the posterior probability of
latent trait level T given estimated model parameters pi and
a case with response pattern X.

```
Recall that in Eq. (1), Pr(T) refers to the prior probability distribution for the latent trait level of any given case. We will adopt a "subjective Bayes" view here (the particular version of "Bayesianism" one adopts makes little practical difference here). That is, we view the prior distribution as representing our prior belief about the distribution of a given case's level of T, before we consider the case's responses.

The obvious choice for the prior distribution of T is the distribution chosen in the process of estimating model parameters. If parameters are estimated via FIML estimation with an assumed normal (or multivariate normal) distribution for T, then that is the logical form of the prior distribution of T. If parameters are estimated by factoring tetrachoric correlations, that again supposes a normal (for a 1-factor model; or multivariate normal for a multifactor model) distribution, which is again the logical choice for the prior distribution of T.

We will shortly consider, step by step, what is involved in calculating the posterior probability distribution P(T|X, pi) and, from that, the EAP estimate of a case's latent trait score.

First, though, we briefly consider what it means that we are seeking an expected value. As noted above, Pr(T|X, pi) defines the posterior probability distribution of T for a given case. What we seek, of course, is not a distribution of scores, but a single best score for that case. With EAP estimation, that single best score is taken as the expected value of the posterior probability distribution. The expected value is merely a weighted average of all possible values of T, the weight being the probability of that level of given X and pi, or:

```
S
E(T|X, pi) = SUM T(s) Pr[T(s)|X, pi]   (2)
s=1

where:  E(T|X, pi) is the EAP estimate of a case's latent trait level;

s indexes one of S intervals into which we have discretized
T for computational purposes;

T(s) is the value of T in the s'th interval.

```
To calculate Pr(T(s)|X, pi), let us first rewrite Eq. 1 as:
```                         Pr[T(s)] Pr[X|T(s), pi]
Pr[T(s)|X, pi] =  ---------------------------    (3)
SUM Pr[T(s)] Pr[X|T(s), pi]

```
Suppose, then, that we have divided T into, say, 20 equally-spaced intervals from T = -3 to T = +3 (we could use more intervals, but that would not make much difference). To apply this equation we must, for a given case, know the values of Pr[T(s)] and Pr[X|T(s), pi] at each value level of T(s).

Pr[T(s)] is known from the prior distribution. Assuming a normal prior with mean of 0 and standard deviation of 1, Pr[(T(s)] is merely the ordinate of the standard normal curve for z = T(s), multiplied times the interval size, here 6/20 = .3.

We need therefore calculate only Pr[X|T(s), pi)]. The formula for this is closely related to the basic definitional formula of the fully Gaussian latent trait model:

```
J
Pr[X|T(s), pi] = PROD  Pr[X(j)|T(s)]      (4)
j=1

```
where j indexes the J dichotomous manifest variables and X(j) is the observed response (0 or 1) of the case on manifest variable j.
```
phi[z(j, s)],      if X(j) = 1
Pr[X(j)|T(s)]   =                                    (5)
1 - phi[z(j, s)],  if X(j) = 0

```
phi() is the cumulative distribution function of the standard normal curve, and z(j, s) is a z score associated with manifest variable j and latent trait level interval s, calculated as:
```
T(s) - tau(j)
z(j, s) =   --------------    (6)
e(j)

where:

tau(j) is the response threshold parameter for manifest variable j;

e(j)   is the standard deviation of measurement error for manifest
variable j, equal to the square root of 1 - r(j)^2, where
r(j) is the correlation parameter for manifest variable j
the latent trait/factor).

```
We now consider the calculation of EAP scores step by step:
1. Divide the latent trait into S intervals (for a unidimensional latent trait, S = 20 is about right; we will assume that here).

2. Calculate the prior probability distribution, Pr[T(s)]. Assuming the normal prior, divide T into 20 intervals in the range -3 to +3. This gives 21 values of T(s), from -3 to +3, inclusive. For each value of T(s), calculate Pr[T(s)] as the ordinate of the standard normal curve at z = T(s), multiplied times the interval size of .3. Save the results in an array of length 21.

3. For each case, do the following:

3a. For each value T(s), calculate Pr[X|T(s), pi)]. Remember that pi is simply the estimated model parameters--i.e., the threshold and correlation parameters of the manifest variables. The calculation is done via Eqs. (4), (5) and (6).

3b. For each value of T(s), multiply Pr[T(s)] times Pr[X|T(s), pi]. This product gives the numerators of Eq. (3). [We say 'numerators'-- i.e., plural, because we will calculate Eq. 3 at every value of T(s).] The denominator of Eq. (3) is the sum of this product at every level T(s).

3c. As explained above, calculate Eq. (3) for every value of T(s). Store results in an array of length 21. This set of values is the posterior probability distribution of the latent trait level for the case, Pr[T(s)|X, pi].

3d. Apply Eq. (2) to obtain the EAP estimate of the latent trait score for the case.

One can follow this method for each individual case. However, a more efficient alternative is sometimes to follow this method for each unique response pattern, then assign all cases with that response pattern the calculated latent trait score. The choice depends on whether there are fewer response patterns than cases or fewer cases than observed response patterns.

### 4. EAP estimation with a multidimensional latent trait

(Only a brief treatment of this subject is attempted at this time).

EAP estimation for a multidimensional latent trait model follows the same basic principles as described above.

For the unidimensional model, integration was approximated by dividing the latent trait continuum into a finite number of discrete levels. For a multidimensional model, one divides the latent trait space into a finite number of regions. For example, with a 2-dimensional model, one might divide each dimension into 20 equally-spaced intervals, then consider each of the 20 x 20 combinations.

Bock, Gibbons and Muraki (1988) suggest that for FIML estimation of a multidimensional latent trait model, it is the total number of regions that affects accuracy. That is, say with a 2-dimensional model, each dimension can be divided into 5 intervals, to produce 5 x 5 = 25 total combinations or regions in the space, and this is comparable in accuracy to having 25 intervals in a unidimensional latent trait model. The same principle likely applies for EAP estimation. Thus it is not necessary to have a large number of intervals for each latent trait dimension. For a multidimensional latent trait model, as few as 5 levels per latent trait dimension may suffice.

The first step, then, is to divide the space into a reasonable number of intervals by the method above.

The next step is to calculate the prior distribution. Most of the previous formulas apply, but we re-interpret the notation. We now understand T(s) to mean a vector of levels across all the latent trait dimensions, and s to be a vector or intervals across the latent trait dimensions. That is, T(s) = {T1(s1), T2(s2), ..., Tm(sm)}, where, say, T1(s1) is the level of latent trait 1 associated with the first level on that dimension, and m is the number of latent trait dimensions.

The prior distribution Pr[T(s)] is the density of the multivariate normal distribution corresponding to T(s).

Because the latent trait dimensions are assumed independent, this is calculated simply as the product of ordinates of m univariate standard normal distributions. That is,

```
m
Pr[T(s)] = PROD Pr[Tk(sk)]       (7)
k=1

```
where Pr[Tk(sk)] is the ordinate of the standard normal curve at z = Tk(sk), multiplied times the interval size for dimension k. Equations (2), (3), (4), and (5) apply as before, after making suitable adjustment for the multidimensional nature of T (for example, in Eqs. 2 and 3, summation is across all combinations of s1, s2, ..., sk.)

Equation (6) is modified however. It now has the form:

```
m
SUM w(k) Tk(sk) tau(j)
k=1
z(j, s) =   ------------------------     (8)
e(j)

```
where w(k) is a weight associated with dimension k. If one is factoring tetrachorics, w(k) is related to the eigenvalue associated with the k'th factor; I would need to think more about this to know exactly how to obtain w(k) here (perhaps Bock & Aitkin, 1981 or Bock, Gibbons & Muraki, 1988 contain clues). As a first guess, I think one might take the eigenvalues for the first k factors, rescale them so that they sum to 1, and use these as the weights; if a rotated solution is being considered, one should base weights on the eigenvalues of the rotated solution.

e(j) is also obtained differently for the multidimensional model. I believe it is calculated as 1 - R(j)^2, where R(j) is the multiple correlation of manifest variable j relative to the k latent trait dimensions. Again, I would need to think more about this (and, again, perhaps this is discussed in either of the Bock papers cited above), but since the latent trait dimensions are orthogonal, then R(j)^2 might just be the sum of r(j,k)^2 over all latent trait dimensions k = 1, ,..., m, where r(j,k)^2 is the correlation parameter of manifest variable j relative to latent trait dimension k.

### References

Bock RD, Aitkin M. Marginal maximum likelihood estimation of item parameters: Application of an EM algorithm. Psychometrika, 1981, 46, 443-459.

Bock RD, Gibbons R, Muraki E. Full-information item factor analysis. Applied Psychological Measurement, 1988, 12, 261-280.

Clogg CC. Latent class models for measuring. In Latent trait and latent class models, Langeheine R, Rost, J. (eds.), pp. 173-205. New York: Plenum, 1988.

Uebersax JS. Statistical modeling of expert ratings on medical treatment appropriateness. Journal of the American Statistical Association, 1993, 88, 421-427.