Understanding Reliability from a Regression Perspective
Yang Liu Jolynn Pek Alberto Maydeu-Olivares




Reliability
[1,2,3]Yang Liu, Jolynn Pek, Alberto Maydeu-Olivares

  Department of Human Development and Quantitative Methodology
University of Maryland, College Park,
  Department of Psychology
The Ohio State University,
  Department of Psychology
University of South Carolina

Correspondence should be made to Yang Liu at
3304R Benjamin Bldg, 3942 Campus Dr, University of Maryland, College Park, MD
20742. Email: yliu87@umd.edu.



Reliability is an important quantification of measurement precision based on a latent variable measurement model. Inspired by mcdonald.2011, we present a regression framework of reliability, placing emphasis on whether latent or observed scores serve as the regression outcome. Our theory unifies two extant perspectives of reliability: (a) classical test theory (measurement decomposition), and (b) optimal prediction of latent scores (prediction decomposition). Importantly, reliability should be treated as a property of the observed score under a measurement decomposition, but a property of the latent score under a prediction decomposition. To facilitate the evaluation and interpretation of distinct reliability coefficients for complex measurement models, we introduce a Monte Carlo approach for approximate calculation of reliability. We illustrate the proposed computational procedure with an empirical data analysis, which concerns measuring susceptibility and severity of depressive symptoms using a two-dimensional item response theory model. We conclude with a discussion on computing reliability coefficients and outline future avenues of research.
  
  

Understanding Reliability from a Regression Perspective
    Samuel Mallick, Azita Dabiri, and Bart De Schutter, Fellow, IEEE
	
	This paper is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon
		2020 research and innovation programme (Grant agreement No. 101018826
		- CLariNet).
	All authors are affiliated with Delft Center for Systems and Control, Delft
		University of Technology, Delft, The Netherlands (e-mail: {s.h.mallick,a.dabiri,b.deschutter}@tudelft.nl).
    April 29, 2024
==============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================



  Constructs (e.g., personality, intelligence, attitudes), which are not directly observable, are often treated as building blocks in psychological theories (e.g., cronbach   meehl.1955). Often, constructs are operationalized as latent variables (LVs), which are indicated by manifest variables (MVs; also commonly referred to as observed variables, indicator variables, or items). For example, MVs can be participants' responses to items in standardized tests and survey questionnaires, which are designed to measure LVs such as cognitive ability and personality. MVs can also be subtest or test scores.

  Observed scores (e.g., summed scores and estimated factor scores) are functions of MVs and are often employed in downstream analysis including scoring, classification, and fitting explanatory models as proxies of latent scores, which are functions of LVs (e.g., see liu   pek.inreview). Because observed scores contain measurement error and measurement error might result in biased inferences (bollen.1989, Chapter 5; cole   preacher.2014), it is pertinent to quantify this error and adjust our interpretation of inferential results accordingly (e.g., see schmidt   hunter.1999). Measurement error can be assessed by two distinct approaches. One approach is to select an observed score and evaluate the precision with which the observed score can be estimated by LVs. An alternative approach is to select a latent score and to investigate the extent to which it can be predicted by MVs. These two approaches result in the two popular definitions of reliability: (a) classical test theory (CTT) reliability and (b) proportional reduction in mean squared error (PRMSE), respectively.

  

  The contribution of our paper is two-fold. First, we provide a road map for organizing different types of common reliability coefficients, clarifying how they relate to one another, and what distinct information they quantify.  Methodological developments in reliability are numerous and nuanced (e.g., revelle   condon.2019), making it challenging for substantive researchers to recognize subtle differences among extant reliability coefficients and employ appropriate coefficients for their research. Thus, our systematization of reliability coefficients, based on McDonald's  regression perspective, aims to facilitate the proper application of reliability coefficients in substantive research. Drawing on this reliability-as-regression perspective, our second contribution is to introduce a novel and straightforward Monte Carlo (MC) procedure that can be used to approximate reliability coefficients for any measurement model. In linear measurement models (e.g., common factor model; thurstone.1935) or in nonlinear measurement models (e.g., IRT; thissen   steinberg.2009) that measure a single LV, reliability estimation is relatively straightforward. We show how our MC procedure provides comparable results to these classical approaches of reliability estimation. More important, we also show that our MC procedure is invaluable when the computation of reliability coefficients becomes intractable in more complex measurement models (e.g., multidimensional IRT models).

  

  

  The paper is structured as follows. We begin by describing reliability from a regression framework, which stems from mcdonald.2011. The framework encompasses two types of decompositions called measurement and prediction decompositions. The measurement decomposition maps onto CTT reliability, and the prediction decomposition maps onto PRMSE. We describe each decomposition using theoretical examples linked to simple measurement models. Next, we introduce an MC procedure for obtaining approximate calculations of reliability coefficients. To illustrate and contrast various types of reliability coefficients, we step through an empirical example on the depressive symptom scale within the Collaborative Psychiatric Epidemiological Surveys (CPES; magnus   liu.2022), making use of the novel MC procedure for computing reliability coefficients. Finally, we discuss implications of our work on reporting reliability coefficients in substantive research and outline future directions of research.

Reliability From a Regression Framework

Assumptions and Notations


We will restrict our discussion of reliability at the level of the population, assuming that the specified measurement model is (or at least closely approximates) the data generating mechanism and that all the parameters in the measurement model are known. In practice, however, parameters in the measurement model must be estimated from a finite sample of MVs. The model can be either linear (e.g,. factor analysis) or nonlinear (e.g., IRT). Let _i be the m× 1 vector of MVs associated with person i, and _i be the d× 1 vector of LVs for person i. Here, m denotes the number of MVs (i.e., the length of the measurement instrument) and d denotes the number of LVs (i.e., dimensionality of the measurement model).  To distinguish between random variables (i.e., variables before they are observed) and their empirical realizations (i.e., variables when they take on fixed values), we underline the random variables.

For ease of presenting theoretical examples and without loss of generality, we focus on unidimensional scores. For person i, let s(_i) be an observed score ( a function of MVs _i) and ξ(_i) be a latent score (a function of LVs _i). Both observed and latent scores are unidimensional as implied by non-bolded expressions of s() and ξ(), respectively. Observed scores are often used to rank or classify individuals as well as get incorporated in downstream statistical analyses. An observed score s(_i) is a function of MVs _i and depends stochastically on LVs _i through the conditional distribution s(_i)|_i. A latent score ξ(_i) can never be observed but can be inferred from the conditional distribution of ξ(_i)|_i. Both conditional distributions can be obtained from the joint distribution  of _i and _i, with probability density function (pdf) f(_i, _i). Stated differently, these conditional distributions are fully determined by the measurement model.
  

  Below, we describe two simple measurement models at the level of the population to illustrate different types of observed and latent scores. After describing these examples, we introduce our regression framework of reliability which encompass CTT reliability and PRMSE.

  Examples and Observed Scores

  Example 1: One-Factor Model


  Consider a common factor model with a single LV η_i (d=1), whose mean equals to 0 and variance equals to φ. The m MVs in _i are associated with the LV by the following equation:
    
    _i =  + η_i + _i,

    in which  is an m× 1 vector of MV intercepts and  is an m× 1 vector of factor loadings. The m× 1 vector of unique factors, denoted by _i, are uncorrelated with η_i. The covariance among the unique factors, (_i, _i) is denoted by , which is an m× m diagonal matrix. As an example, consider a one-factor model for three MVs (m = 3) with parameters
    
    = [ 0; 0; 0 ], = [ 0.3; 0.5; 0.7 ], φ = 1, and  = [ 0.91    0    0;    0 0.75    0;    0    0 0.51;      ].

     Equation <ref> implies that the LV η_i and the three MVs in _i = [y_i1, y_i2, y_i3] are standardized.

   Due to linearity of the one-factor model, commonly used observed scores are weighted or unweighted sums of MVs. Let us first consider the regression factor score <cit.>,
    
    s_1(_i) = ^-1_i/^-1 + φ^-1,

    which is a weighted sum of the elements in _i. Let us further assume that the LV follows a normal distribution, η_i∼ N(0, φ), and that the unique factors follow (independent to the LV) a normal distribution as well, _i∼ N( 0_m, ), in which 0_m denotes an m× 1 vector of zeros. Taken together, the regression factor score s_1(_i) coincides with the conditional mean of η_i given _i, 𝔼(η_i|_i). The conditional mean is also known as the expected a posteriori [EAP] score of the LV η_i <cit.>. As such, we use the term “regression factor scores” interchangeably with “EAP scores” in the linear one-factor model example. The second observed score under consideration is the (unweighted) summed score of the m MVs,
    
    s_2(_i) = _m_i,

    in which _m is an m× 1 vector of ones.
    Given the model parameter values in Equation <ref>, the regression factor score for individual i is s_1(_i) = 0.14y_i1 + 0.28y_i2 + 0.57y_i3 (by Equation <ref>) and the summed score is s_2(_i) = y_i1 + y_i2 + y_i3 (by Equation <ref>).

    Under a unidimensional measurement model, the latent scores of interest can be the LV itself; i.e., ξ_1(η_i) = η_i. Alternatively, it might be desirable to align the scale of the latent score to the scale of an observed score (e.g., the summed score). For instance, let the expected summed score be
    
    ξ_2(η_i) = [s_2(_i)|η_i] = _m( + η_i),

    in which [·|η_i] is the conditional expectation over the MVs given the LV η_i. Because ξ_2(η_i) is a linear transformation of η_i, ξ_2(η_i) is perfectly correlated with η_i. Substituting  and  with population parameter values from Equation <ref> results in Equation <ref> becoming ξ_2(η_i) = 1.5η_i. 

Example 2: Two-Parameter Logistic Model


   In our second example, the MVs are dichotomous item responses. Conditional on a single LV η_i, assume that the MVs y_ij, j = 1,…, m, are independent (i.e., local independence; mcdonald.1994; stout.2002). The conditional probability of y_ij = k given η_i, where k ∈{0, 1}, is parameterized by a two-parameter logistic (2PL) model <cit.>:

    f_j(k|η_i) = {y_ij = k|η_i} = exp[k(α_j + β_jη_i)]/1 + exp(α_j + β_jη_i),

  in which α_j and β_j are the intercept and slope parameters for the jth MV, respectively. Equation <ref> is commonly referred to as the item response function (IRF). We further assume that the LV η_i follows a standard normal distribution: η_i∼ N(0, 1). As a numerical example, consider a test composed of m = 3 dichotomous MVs with

    = [  1;  0; -2 ], and  = [   1; 1.5;   2 ],

  in which  and  contain the intercept and slope parameters, respectively.



   Similar to the linear one-factor model example, we consider two types of observed scores. The first type of score is the EAP score of η_i based on the 2PL model:

    s_1(_i) = ( η_i|_i ) = ∫η_i f(_i|η_i)ϕ(η_i)dη_i/f(_i),

in which

    f(_i|η_i) = ∏_j=1^m f_j(y_ij|η_i)

  denotes the conditional pdf of _i given η_i and ϕ(η) denotes the standard normal pdf of η_i. Additionally, 
    f(_i) = ∫ f(_i|η_i)ϕ(η_i)dη_i

denotes the marginal pdf of _i. The integrals in Equations <ref> and <ref> do not have closed-form expressions and are often approximated by numerical quadrature in practice.[Consistent with the default configuration in the R package  <cit.>, we always use 61 equally-spaced quadrature nodes ranging from -6 to 6 to approximate intractable integrals in the 2PL example.] With three binary MVs, there are 2^3 = 8 response patterns. Using item parameter values specified in Equation <ref>, we computed the marginal probabilities and EAP scores for these 8 response patterns, which are summarized in Table <ref>. The second type of score is the summed score s_2(_i) = _m_i (presented in the fourth column of Table <ref>). Note that there are only four distinct values for the summed score for the eight response patterns in a three-item test (see Table <ref>).



  In addition to latent scores that represent the LV, ξ_1(η_i) = η_i, the expected summed score (also known as the test characteristic curve; <e.g.,>[pp. 159-160]thissen   wainer.2001) is another popular latent score in IRT applications. The expected summed score is expressed as
    
    ξ_2(η_i) = [s_2(_i)|η_i] = ∑_j=1^mf_j(1|η_i).

    In contrast to Equation <ref> for the linear one-factor model example, Equation <ref> is often a nonlinear function of the LV, η. Figure <ref> is a visualization of Equation <ref> using values of the slope and intercept parameters specified in Equation <ref>. When all the slope parameters are positive, ξ_2(η_i) is a nonlinear and strictly increasing function of η_i, which is bounded between 0 and the total number of MVs (i.e., m = 3 in the current example).

    In the next section, we review two approaches to defining reliability coefficients, mirroring the measurement and prediction decompositions by mcdonald.2011. Both decompositions stem from regression equations but differ in whether a latent score or an observed score serves as the outcome variable. The measurement decomposition yields the well-known definition of reliability coefficients in CTT, while the prediction decomposition introduces PRMSE as a measure of reliability that is more frequently seen in the IRT literature.

Measurement Decomposition


Classical Test Theory


  We adapt classical definitions and results in Lord and Novick (lord   novick.1968, p. 34) to the context of LV measurement models using our notation to describe the measurement decomposition of reliability. Recall that an underlying LV model determines the joint distribution of MVs _i and LVs _i. It follows that the conditional distribution of an observed score s(_i) given the LVs _i, denoted by s(_i)|_i, is uniquely identified for every _i. This conditional distribution is called the propensity distribution of the observed score s(_i) (cf. lord   novick.1968, p. 29–38), which reflects the variability of observed scores across independent and identically distributed (i.i.d.) measurement instances for the same person i.[In a different view, the conditional distribution s(_i)|_i captures the variability of observed scores across persons with the same LV value. Correspondingly, the conditional expectation [s(_i)|_i] is the average of observed scores within the subpopulation at the LV level _i. See holland.1990 for relevant discussions.] The true score corresponding to s(_i) is defined as the expectation of the propensity distribution; i.e., [ s(_i)|_i ] (cf. Lord & Novick, 1968, Definition 2.3.1). By definition, the true score is the long-run average of observed scores across repeated i.i.d. measurements of the same person, which can be conceived as a version of the observed score that is free of measurement error. The error score for person i, ε_i (cf. lord   novick.1968, Definition 2.4.1.), is defined as the difference between the observed score and the true score

    
    ε_i = s(_i) - [ s(_i)|_i].


    Over the population of persons, the error score has a mean of zero and is uncorrelated with the true score. That is, ε_i = 0 and  (ε_i, [s(_i)|_i]) = 0. Note that _i is underlined in [s(_i)|_i] in the covariance expression because the conditional expectation [s(_i)|_i] is a function of _i and _i is treated as random. CTT reliability of s(_i), denoted as Rel_ CTT(s(_i)), can be defined in two main ways described below.

    First, CTT reliability quantifies the proportion of observed score variance that can be explained by true scores. Mathematically, CTT reliability is the ratio of true score variance over observed score variance, or one minus the ratio of error score variance over observed score variance:

    
    Rel_ CTT(s(_i)) = ([s(_i)|_i])/[s(_i)] = 1 - (ε_i)/[s(_i)].

   The second equality in Equation <ref> follows from the fact that the true score is uncorrelated with the error score. Second, assuming that the true score variance ([s(_i)|_i]) is positive, CTT reliability also equals to the squared correlation between the observed score s(_i) and the true score [s(_i)|_i]:
    
    Rel_ CTT(s(_i)) = (s(_i), [s(_i)|_i])^2.

    

    In sum, CTT reliability can be interpreted in two exchangeable ways: a ratio of true- and observed-score variances (Equation <ref>) and a squared-correlation between true and observed scores (Equation <ref>).[A third definition of CTT reliability can be defined using the concept of a parallel score, which is less relevant to the present paper. Let s̃_i be another observable score produced by person i. Then, s̃_i is parallel to s(_i) if (a) (s̃_i|_i) = [s(_i)|_i] (i.e., the true scores for s̃ and s(_i) are equal), (b) (ε_i) = (ε̃_i) (i.e., the variance of error scores are equal), where ε̃_i = s̃_i - (s̃_i|_i) denotes the error score of s̃_i, and (c) (ε_i, ε̃_i) = 0 (i.e., the error scores are uncorrelated; lord   novick.1968, p. 47). Thus, CTT reliability equals to the correlation between two parallel test scores (cf. lord   novick.1968, p. 58):
    
      Rel_ CTT(s(_i)) = [s(_i), s̃_i].
      
    
    The parallel score s̃_i, often regarded as a concept, might not be necessarily computed from MVs. A more intuitive approach to construct parallel scores, which requires stronger assumptions, was described by Lord (lord.1983; see also kim.2012). Let _i denote an i.i.d. copy of the MV vector _i conditional on _i; we may think of _i as independent responses to an equivalent form of the same measurement instrument produced by the same person. lord.1983 directly set s̃_i = s(_i), which satisfies all the three requirements of a parallel score; the independence requirement, however, is stronger than having uncorrelated error scores.]
 

    A Regression Formulation


    Rearranging Equation <ref> leads to the widely known true score formula:
 
    
    s(_i) = [ s(_i)|_i ] + ε_i,

    or equivalently in words,
    
    observed score = true score + measurement error.

 
    
    Equation <ref> is hereinafter termed the measurement decomposition of reliability <cit.>. Recall that a regression traces the conditional expectation of an outcome variable given one or more explanatory variables (e.g., fox.2015, p. 15). Thus, Equation <ref> can be viewed as the regression of the observed score s(_i) on all the LVs _i. This regression is in general nonlinear because the predicted value of the model, which coincides with the true score [s(_i)|_i] might not be a linear function of _i (e.g., in a 2PL model). If we regard the true score as the explanatory variable, then Equation <ref> can be alternatively understood as a unit-weight linear regression (i.e., with an intercept of zero and a slope of one) of the observed score s(_i) on the true score [s(_i)|_i]. This second regression interpretation follows from the towering property of the conditional expectation: {s(_i)|[s(_i)|_i]} = {[s(_i)|_i]|[s(_i)|_i]} = [s(_i)|_i]. Thus, when the observed score s(_i) is regressed onto the true score [s(_i)|_i], the predicted value is simply the true score itself, leading to a unit-weight regression.

    Let ϱ^2(u, ) denote the population coefficient of determination when regressing an outcome variable u on explanatory variables : The regression is implied by the joint distribution of u and  and is potentially nonlinear. Both views of Equation <ref> (i.e., regression of the observed score on all the LVs and regression of the observed score on the true score) yield the same coefficient of determination that equals to CTT reliability (Equation <ref>):
    
    ϱ^2(s(_i), _i) = ϱ^2(s(_i), [s(_i)|_i]) = Rel_ CTT(s(_i)).

    In words, reliability from the measurement decomposition perspective is the amount of variance in the observed score s(_i) explained by the latent score _i (i.e., treating _i as explanatory variables) or by the true score [s(_i)|_i] (i.e., treating [s(_i)|_i] as the explanatory variable). Implicit in this interpretation is the assumption that the chosen observed score (which is a summary statistic of the MVs) reflects underlying LVs (constructs). Stated differently, the LVs give rise to the MVs that are eventually combined to form an observed score (cf. borsboom   mellenbergh.2002).

    

    Below, we illustrate how to calculate CTT reliability from the regression framework for the examples of a one-factor model (Equation <ref> with parameter values in Equation <ref>) and a 2PL model (Equation <ref> with parameter values Equation <ref>). Within each model, CTT reliability coefficients are computed for the EAP score s_1(_i) and the summed score s_2(_i).

Example 1: One-Factor Model


  Consider first the regression factor score, which is also the EAP score under the normality assumption (i.e., s_1(_i); Equation <ref>). By Equations <ref> and <ref>, the true score of s_1(_i) is

    [s_1(_i)|η_i] =
      ^-1( + η_i)/^-1 + φ^-1.

  That is, the regression of the EAP score s_1(_i) on the LV η_i is linear. With covariance algebra, it can be shown that the observed score variance of the EAP score s_1(_i) is
    
    [s_1(_i)] = φ^-1/^-1 + φ^-1,

  and the corresponding true score variance (over the population distribution of η_i) is

    ([s_1(_i)|η_i]) =
    φ[^-1/^-1 + φ^-1]^2.

  Taking the ratio of Equation <ref> over Equation <ref> yields CTT reliability for the EAP score:

    ϱ^2(s_1(_i), _i) = ([s_1(_i)|η_i])/[s_1(_i)] = ^-1/^-1 + φ^-1.

Because the true score underlying the EAP score s_1(_i) is a linear function of the LV η_i (see Equation <ref>), Equation <ref> is not only the squared correlation between s_1(_i) and its true score (i.e., the squared-correlation interpretation of CTT reliability) but also the squared correlation between s_1(_i) and η_i. When evaluated at the model parameter values in Equation <ref>, CTT reliability for the regression factor score is 0.58.

  The regression of the summed score s_2(_i) on the latent score η_i is also linear. Observe that the true score underlying s_2(_i), given by Equation <ref>, is a linear function of η_i. The observed score variance and the true score variance associated with summed score s_2(_i) are
    
    [s_2(_i)] = φ(_m)^2 + tr(),

    and

    ([s_2(_i)|η_i]) =
      φ(_m)^2,

   respectively. CTT reliability for the summed score, expressed as the ratio of Equation <ref> over Equation <ref> is
    
    ϱ^2(s_2(_i), η_i) = ([s_2(_i)|η_i])/[s_2(_i)] = φ(_m)^2/φ(_m)^2 + tr().

    Equation <ref> is widely known as the coefficient omega <cit.>.   Due to the linearity of the regression, Equation <ref> gives not only the squared correlation between the summed score and its underlying true score but also the squared correlation between the summed score and the LV. With values specified for φ, , and  in Equation <ref>, CTT reliability for the summed score (Equation <ref>) is equal to 0.51, which is noticeably lower than that of the EAP score (i.e., 0.58).

    
    

Example 2: Two-Parameter Logistic Model

   Reliability calculations are typically less tractable for nonlinear measurement models. As will be shown, calculations for the simple 2PL example become rather involved. Thus, we introduce a more intuitive and widely applicable numerical recipe for approximate reliability calculation for nonlinear measurement models in the section “A Monte Carlo Procedure."

  Consider the EAP score of η_1, denoted s_1(_i). The true score underlying s_1(_i) is

    [s_1(_i)|η_i] = ∑__is_1(_i) f(_i|η_i),

  in which the summation is taken over all 2^m possible response patterns. The regression of the EAP score s_1(_i) on the LV η_i is nonlinear due to the logistic IRFs on the right-hand side of Equation <ref>. We proceed by computing the  error score and observed score variances. The error score variance can be expressed as
    
    (ε_i) = ([s_1(_i)|η_i]) = ∫[s_1(_i)|η_i]ϕ(η_i)dη_i,

    in which the conditional variance of the EAP score s_1(_i) given the LV η_i is
    
    [s_1(_i)|η_i] = ∑__is_1(_i)^2f(_i|η_i) - [ ∑__is_1(_i) f(_i|η_i) ]^2.

    Meanwhile, the observed variance of the EAP score s_1(_i) can be computed as
    
    [s_1(_i)] = ∑__is_1(_i)^2f(_i) - [  ∑__is_1(_i)f(_i) ]^2.

    CTT reliability for the EAP score s_1(_i) is then equal to one minus the ratio of Equation <ref> over Equation <ref>:

    ϱ^2(s_1(_i), η_i) = 1 - ([s_1(_i)|η_i])/[s_1(_i)] = 1 - ∫[s_1(_i)|η_i]ϕ(η_i)dη_i/∑__is_1(_i)^2f(_i) - [  ∑__is_1(_i)f(_i) ]^2.

In practice, Equation <ref> can only be evaluated approximately as integrals involved in Equation <ref> and the marginal pdf of the MVs f(_i) defined by Equation <ref>, and are computationally intractable. Using the item parameters in Equation <ref> and numerical quadrature to approximate integrals, we obtain 0.51 as CTT reliability of the EAP score s_1(_i) in the three-item example. When the EAP score s_1(_i) is the observed score of interest, its underlying true score (Equation <ref>) is a nonlinear function of the LV η_i. Hence, the resulting CTT reliability (Equation <ref>), while remaining to be the squared correlation between the observed score and its underlying true score, is no longer the same as the squared correlation between the observed score and the LV. As we will see in the “Prediction Decomposition” section, the squared correlation between the observed score and the LV coincides with PRMSE of the LV η_i.

  CTT reliability for the summed score s_2(_i) can be calculated in the same fashion.

    ϱ^2(s_2(_i), η_i) = 1 - ([s_2(_i)|η_i])/[s_2(_i)] = 1 - ∫[s_2(_i)|η_i]ϕ(η_i)dη_i/∑__is_2(_i)^2f(_i) - [  ∑__is_2(_i)f(_i) ]^2.


   The conditional variance of the summed score (cf. Equation <ref> for the EAP score) has the following simpler expression:

    [s_2(_i)|η_i] = ∑_j = 1^mf_j(0|η_i)f_j(1|η_i),

which is derived from assuming local independence and using the variance formula for a Bernoulli random variable. In contrast to the one-factor example, CTT reliability of the summed score s_2(_i) is different from the squared correlation between the summed score and the LV η_i, for the reason that the true score underlying the summed score (Equation <ref>) is a nonlinear transformation of the LV. Under the 2PL model with item parameters in Equation <ref>, CTT reliability of the summed score is 0.5, which is only slightly lower than that of the EAP score (i.e., 0.51).

Prediction Decomposition


Optimal Prediction of a Latent Score

  
The prediction decomposition of reliability concerns using MVs to predict a latent score, ξ(_i), which is a function of the LVs _i. The uncertainty in prediction can be quantified by the mean squared error (MSE). Generally speaking, a predictor is (a possibly constant) function of MVs and is optimal when it minimizes the MSE. Let us first consider making predictions without MVs (i.e., using only a constant predictor g∈). A well-known result in elementary statistics is that the optimal predicting constant of the latent score ξ(_i) is its unconditional mean ξ(_i) <cit.>e.g.,>[Example 2.2.6]casella   berger.2002, and the minimized MSE is
    
    min_g∈[ξ(_i) - g]^2 = [ξ(_i) -  ξ(_i)]^2 = [ξ(_i)].

    Now suppose that we employ the MVs _i to predict the latent score ξ(_i) by considering predictors of the form g(_i), where g denotes a real-valued function of the MVs. Similar to Equation <ref>, it can be shown that the MSE-minimizing predictor is [ξ(_i)|_i], i.e., the EAP score of ξ(_i) (e.g., casella   berger.2002, Exercise 4.13). The minimized MSE then equals to
    
    min_g[ξ(_i) - g(_i)]^2 = {ξ(_i) - [ξ(_i)|_i]}^2 = ([ξ(_i)|_i]).

    

    When predicting ξ(_i) with a constant, the minimized MSE (Equation <ref>) coincides with the unconditional variance of the latent score, [ξ(_i)]. This unconditional variance can be interpreted as the total amount of variability (uncertainty) due to individual differences in ξ(_i). When the MVs _i are used for predicting the latent score ξ(_i), the minimized MSE, which can be understood as the remaining uncertainty, is reduced to ([ξ(_i)|_i]). The amount of uncertainty reduction, expressed as [ξ(_i)] - ([ξ(_i)|_i]) = ([ξ(_i)|_i]), thus reflects the degree to which we can precisely recover the latent score using MVs. From this perspective, measurement error can be understood as the remaining uncertainty in the latent score ξ(_i) that cannot be captured by the MVs _i. When the remaining uncertainty in prediction is zero, the latent score can be reproduced from the MVs without error. Taken together, define 
    
    PRMSE(ξ(_i)) = 1 - ([ξ(_i)|_i])/[ξ(_i)] = ([ξ(_i)|_i])/[ξ(_i)].

    
   PRMSE quantifies the proportion of latent score variance that can be explained by the MVs. 
   Assuming that the EAP score has a positive variance, ([ξ(_i)|_i]) > 0, it can also be shown that PRMSE(ξ(_i)) = (ξ(_i), [ ξ(_i)|_i ])^2 (e.g., kim.2012); that is, PRMSE equals to the squared correlation between the latent score ξ(_i) and its EAP score [ξ(_i)|_i].
  
  
  Note further that PRMSE for predicting ξ(_i) cannot exceed the CTT reliability of the EAP score [ξ(_i)|_i], which follows from Equation 31 of  (kim.2012, p.160).

    A Regression Formulation


    PRMSE can also be approached from a regression perspective. In contrast to the measurement decomposition in which an observed score serves as the regression outcome (Equation <ref>), the prediction decomposition treats the latent score as the outcome <cit.> in
    
    ξ(_i) = [ ξ(_i)|_i ] + δ_i,

    or equivalently in words,
    
    latent score = EAP score + prediction error.

    The prediction error δ_i is defined as the difference between the latent score ξ(_i) and its EAP score [ ξ(_i)|_i ]. Equation <ref> implies that the prediction error has mean zero (i.e., δ_i=0) and is uncorrelated with the EAP score (i.e., (δ_i, [ξ(_i)|_i])=0) in the population; these properties are analogous to the facts that the error score has mean zero (i.e., ε_i = 0) and is uncorrelated with the true score (i.e., (ε_i, [s(_i)|_i])=0). The orthogonality of prediction error with the EAP score allows Equation <ref> to be interpreted as a regression of the latent score ξ(_i) onto all the MVs _i. Treating the EAP score as the explanatory variable in the prediction decomposition, Equation <ref> can be alternatively viewed as a unit-weight linear regression of the latent score ξ(_i) on the EAP score [ ξ(_i)|_i ].

    The prediction decomposition of the latent score ξ(_i), whether viewed as a potentially nonlinear regression on _i or a unit-weight linear regression on the EAP score, [ ξ(_i)|_i ], yields a coefficient of determination that equals to the PRMSE:
    
    ϱ^2(ξ(_i), _i) =  ϱ^2(ξ(_i), [ ξ(_i)|_i ]) =  PRMSE(ξ(_i)).

    In words, reliability from the prediction decomposition is interpreted as the amount of variance in the latent score ξ(_i) explained by the MVs _i or equivalently the EAP estimator of the latent score, denoted [ ξ(_i)|_i ]. Different from the measurement decomposition, the prediction decomposition does not emphasize the reflective nature of the measurement model (i.e., MVs are manifestations of underlying LVs). 
    

    

    Illustrations on calculating the PRMSE reliability coefficient for the one-factor and 2PL models are presented in the next section. For each model, we calculate PRMSE for the LV ξ_1(η_i) = η_i and the expected summed score ξ_2(η_i) (Equation <ref>).

Example 1: One-Factor Model


  Assuming that the common factor η_i and the unique factors _i (Equation <ref>) are normally distributed, the EAP estimate of ξ_1(η_i) = η_i is identical to the regression factor score s_1(_i) <cit.>. The expression of s_1(_i) under a one-factor model (i.e., Equation <ref>) implies that the regression of the LV η_i on the MVs _i is linear. Note that the variance of the regression factor score [s_1(_i)] has been derived in Equation <ref> and the variance of the latent score is [ξ_1(η_i)] = (η_i) = φ. Thus, the PRMSE of ξ_1(η_i) is the ratio of these two variances:

    ϱ^2(ξ_1(η_i), _i)  = [s_1(_i)]/[ξ_1(η_i)] = ([s_1(_i)|η_i])/[s_1(_i)] = ^-1/^-1 + φ^-1.

  Note that PRMSE for ξ_1(η_i) = η_i under a linear normal one-factor model (Equation <ref>) is equivalent to CTT reliability for the regression factor score s_1(_i) ((Equation <ref>). This is to be expected because the regression of s_1(_i) on the LV η_i is linear (expressed by Equation <ref>) and the regression of the LV η_i on s_1(_i) is also linear (with unit weight; see Equation <ref>). Taken together, the two coefficients of determination from the measurement and prediction decompositions coincide and are equal to the squared Pearson correlation between the regression factor score s_1(_i) and the LV η_i.[An alternative explanation is that s_1(_i) and η_i follow a bivariate normal distribution.] With specific model parameter values in Equation <ref>, PRMSE of the latent score ξ_1(η_i) = η_i is 0.58 and identical to CTT reliability of the regression factor score s_1(_i).

    Recall that the expected summed score ξ_2(η_i) for the linear one-factor model is also a linear transformation of the LV η_i (Equation <ref>). Thus, PRMSE(ξ_2(η_i)) = PRMSE(ξ_1(η_i)) = 0.58 for this numerical example. Note that PRMSE of ξ_2(η_i) is usually distinct from (and higher than) CTT reliability of the summed score s_2(_i) (Equation <ref>; cf. PRMSE of ξ_1(η_i) is identical to CTT reliability of s_1(_i)).  To see why, recall that the summed score s_2(_i) and the expected summed score ξ_2(η_i) follow a bivariate normal distribution jointly. Regressing s_2(_i) on ξ_2(η_i) and regressing ξ_2(η_i) on s_2(_i) result in the same coefficient of determination. While the formal coefficient of determination is CTT reliability of s_2(η_i), the latter coefficient is in general smaller than PRMSE of ξ_2(η_i). This is because the summed score s_2(_i) is a summary statistic of the MVs _i (i.e., summed scores have less information than _i), and a regression onto the summed score s_2(_i) usually produces a smaller coefficient of determination than a regression onto the MVs _i.

Example 2: Two-Parameter Logistic Model


  For the first latent score that is the LV itself (i.e., ξ_1(η_i) = η_i), its EAP estimator is s_1(_i), which we have defined in Equation <ref>. The variance of s_1(_i) has been expressed in Equation <ref>, and the variance of the LV η_i is one. As such, PRMSE of the first latent score ξ_1(η_i) is

    ϱ^2(ξ_1(η_i), _i) = [s_1(η_i)|_i]/(η_i) = ∑__is_1(_i)^2f(_i) - [  ∑__is_1(_i)f(_i) ]^2.

Table <ref> presents the response-pattern probabilities of the MVs (i.e., f(_i)) and the EAP estimates of ξ_1(η_i) (i.e., s_1(_i)) for the three-item 2PL example. Plugging in the specified values of the population parameters (Equation <ref>) into Equation <ref>, we have PRMSE(ξ_1(η_i)) = 0.5.

  In the prediction decomposition (Equation <ref>) of the expected summed score ξ_2(η_i) (Equation <ref>), the predicted value is its EAP estimate:

    s_3(_i) = [ξ_2(η_i)|_i] = ∫[ ∑_j=1^mf_j(1|η_i) ] f(_i|η_i)ϕ(η_i)dη_i/f(_i),

which we momentarily denote by s_3(_i) for ease of reference. Similar to Equation <ref>, the variance of s_3(_i) over the population distribution of _i is given by

    [s_3(_i)] = ∑__is_3(_i)^2f(_i) - [  ∑__is_3(_i)f(_i) ]^2.

  In addition, the variance of the expected summed score ξ_2(η_i) can be calculated as

    [ξ_2(η_i)] = ∫ξ_2(η_i)^2ϕ(η_i)dη_i - [ ∫ξ_2(η_i)ϕ(η_i)dη_i ]^2,

  which does not possess a closed-form expression. PRMSE for the expected summed score ξ_2(η_i) is then the ratio of Equation <ref> to Equation <ref>:

    ϱ^2(ξ_2(η_i), _i) = [s_3(_i)]/[ξ_2(_i)] = ∑__is_3(_i)^2f(_i) - [  ∑__is_3(_i)f(_i) ]^2/∫ξ_2(η_i)^2ϕ(η_i)dη_i - [ ∫ξ_2(η_i)ϕ(η_i)dη_i ]^2.

  For the 2PL example, we first compute Equation <ref> (i.e., s_3(_i), the EAP estimates of expected summed scores) for all eight response patterns as ordered in Table <ref>:

    (0.80, 1.11, 1.28, 1.45, 1.63, 1.80, 1.97, 2.29).

Then Equation <ref> (i.e., the variance of s_3(_i)) and Equation <ref> (i.e., the variance of ξ_2(η_i)) are estimated as 0.24 and 0.46, respectively. Plugging these values into Equation <ref>, we obtain PRMSE of the expected summed score ξ_2(η_i) as 0.52.

A Monte Carlo Procedure

Estimating Reliability by Simulation

In principle, reliability coefficients compatible with the regression framework can be directly evaluated as functions of model parameters (as we have done for the two simple illustrative examples). When the measurement model is complex (e.g., a multidimensional IRT model), however, computational difficulties in evaluating expectations over the spaces of LVs and MVs ensue. For example, expectations with respect to LVs can be intractable and numerical integration can become inefficient when the latent dimensionality is high. In addition, expectations with respect to discrete MVs often translate to a finite sum across all possible response patterns whose number grows exponentially as the number of MVs increases. As a workaround to such computational challenges, reliability coefficients can be approximated by simulation.


  
  * Simulate latent scores and MVs. Given a measurement model with known parameters (e.g., using estimated parameters from data), simulate a large sample of i.i.d. LV vectors _i, i = 1,…, M, in which M is the Monte Carlo sample size. Then for each i, simulate an MV vector _i conditional on _i.
  
  * Nonparametric regression. For a measurement decomposition, compute the observed score s(_i) for each i = 1, …, M and fit a nonparametric regression predicting the observed scores by all the LVs. For a prediction decomposition, compute the latent score ξ(_i) for each i = 1, …, M and then fit a nonparametric regression predicting the latent scores by all the MVs.
  
  Step 2'. Unit weight regression. For a measurement decomposition, compute the observed score s(_i) and its true score [s(_i)|_i] for each i = 1,…, M, followed by fitting a unit-weight regression predicting the observed scores by the true scores. For a prediction decomposition, compute the latent score ξ(_i) and its EAP estimate [ξ(_i)|_i] for each i = 1,…, M and then fit a unit-weight regression predicting the latent scores by the EAP estimates.
  
  * Estimate reliability coefficient. Obtain the sample coefficient of determination (e.g., the R^2 statistic) from the fitted regression, which is an estimate of reliability.


  Step 1 produces a large sample of LVs and MVs based on a specified measurement model, which can be utilized to approximate the joint distribution of MVs _i (and thus observed scores s(_i)) and LVs _i (and thus latent scores ξ(_i)). Step 2 and Step 2' involve fitting regressions to the MC samples. In practice, only one of the two steps needs to be performed. Step 2 is generally applicable when a nonparametric regression can be estimated with high precision. As an alternative, Step 2' is useful when the computation of true scores or EAP estimates is viable, avoiding the fitting of nonparametric regression. The explanatory variables in the regression may be discrete; for example, in measurement decompositions with discrete LVs (e.g., latent profile analysis) and prediction decompositions with discrete MVs (e.g., binary item responses). In cases with discrete MVs, Step 2 reduces to fitting a linear regression to dummy-coded patterns for explanatory variables.

Throughout the rest of the paper, all MC-based reliability estimates are computed from M = 10^6 MC samples in Step 1, which ensures negligible Monte Carlo error. In Step 2, the R package  <cit.> was used to fit nonparametric regressions and the default thin-plate spline smoother <cit.> was adopted. R code for numerical illustrations in the present paper will be provided as Supplemental Material when the paper is accepted for publication. 

  

  

  Tables <ref> and <ref> present MC-based reliability estimates for the linear one-factor model and the 2PL model, respectively. MC estimates (i.e., R^2) are juxtaposed against their population values that have been presented in earlier sections. Observe that the maximum absolute deviation between the MC-based estimates and population values is less than 0.001, demonstrating accuracy of approximation for these examples. Moreover, Step 2 (i.e., nonparametric regression) and Step 2' (unit weight regression) yield estimates that coincide up to the fourth digit. 



Empirical Example

  The two simple examples described above were designed for two conceptual purposes. First, these examples illustrated the distinction between the measurement and prediction decompositions of reliability, emphasizing the regression framework. Second, these examples demonstrated the accuracy of the novel MC algorithm. When the true parameter values are known, the MC-based reliability estimates closely approximated their population counterparts (see Tables <ref> and <ref>). Next, we apply the MC algorithm to an empirical example to illustrate calculating reliability coefficients from the measurement and prediction decompositions in a more realistic and complex measurement model.

Data and Measurement Model


  The data are from the National Comorbility Survey Replication (NCS-R), which is part of  the Collaborative Psychiatric Epidemiological Surveys <cit.>CPES;>cpes. We focus on a 14-item scale that measures depressive symptoms in the past 30 days. Each item was rated on a scale of 0 = “never,” 1 = “rarely,” 2 = “sometimes,” and 3 = “often.” In magnus   liu.2022, a random subset of 3000 complete responses was analyzed to illustrate the multidimensional hurdle graded response model (MH-GRM). Here we perform reliability calculations on the same data subset (i.e., n = 3000 and m = 14), which was retrieved from https://osf.io/frjm6/https://osf.io/frjm6/.

  The MH-GRM was developed to account for an excessive amount of “never” endorsements in the data set <cit.>, and is a special case of the item response tree model <cit.>e.g.,>jeon   deboeck.2016. Let each MV y_ij be the outcome of a two-stage decision process, in which y_ij∈{0, 1, 2, 3}, j = 1,…, 14. The first stage represents the presence or absence of the symptom. Let y_ij^(1)∈{0, 1} be a dichotomous indicator of the first stage such that y_ij^(1) = 0 maps onto y_ij = 0 (i.e., the “never” category), and y_ij^(1) = 1 maps onto  y_ij > 0 (i.e., one of the three categories indicating symptom presence). The second stage determines the frequency of the symptom, which is encoded by a trichotomous indicator y_ij^(2)∈{1, 2, 3}. Assuming that the symptom is present (i.e., y_ij^(1) = 1, set y_ij = y_ij^(2)). Conversely, y_ij^(1) and y_ij^(2) can be recoded from the item response y_ij by

    y_ij^(1) = {[       0, y_ij = 0;       1, y_ij > 0 ]. and  
      y_ij^(2) = {[      𝙽𝙰, y_ij = 0;    y_ij, y_ij > 0 ].

  If the symptom is absent (i.e., y_ij = y_ij^(1) = 0) no information about y_ij^(2) can be inferred from y_ij and thus we code y_ij^(2) as missing (i.e., ) in Equation <ref>. The two-stage decision process is illustrated in Figure <ref>A.

  Under the MH-GRM model, each individual i is characterized by two LVs, denoted by _i = (η_i1, η_i2). The susceptibility LV η_i1 is indicated by the symptom presence indicators y_i1^(1), …, y_i, 14^(1). The conditional distribution y_ij^(1)|η_i1 is characterized by a 2PL model:

    {y_ij^(1) = 1|η_i1} = exp(α_j^(1) + β_j^(1)η_i1)/1 + exp(α_j^(1) + β_j^(1)η_i1),

  in which α_j^(1) and β_j^(1) represent the intercept and slope corresponding to the susceptibility LV η_i1. The severity LV η_i2 is indicated only by the symptom frequency indicators y_i1^(2), …, y_i,14^(2) through a graded response model:

    {y_ij^(2) = k|η_i2} = exp(α_j,k-1^(2) + β_j^(2)η_i2)/1 + exp(α_j,k-1^(2) + β_j^(2)η_i2) - exp(α_jk^(2) + β_j^(2)η_i2)/1 + exp(α_jk^(2) + β_j^(2)η_i2), k = 1, 2, or 3.

  In Equation <ref>, α_jk^(2)'s and β_j^(2) represent the intercepts and slope corresponding to the severity LV η_i2, in which the intercepts must be in a descending order; i.e., ∞ = α_j0 > α_j1 > α_j2 > α_j3 = -∞. Additionally, η_i1 and η_i2 are allowed to covary (see Figure <ref>B). The estimated correlation between susceptibility and severity is 0.58, implying that they are related but distinctive constructs. magnus   liu.2022 also reported that EAP scores for susceptibility and severity exhibited different patterns in predicting health-related outcomes. For instance, when both sets of scores were used in a logistic regression to predict attempted suicide, only severity scores had a significant partial effect.



  To illustrate various reliability coefficients organized under the regression framework, we treat estimated parameters as known and ignore sampling variability because of the large sample size (n = 3000). Slightly different from magnus   liu.2022, in which model parameters are estimated by directly maximizing the marginal likelihood, we recoded the raw item response data using Equation <ref> and fit the independent cluster model to the recoded data using the R package  <cit.>. The software package implements an Expectation-Maximization algorithm for parameter estimation <cit.>, and the default configuration of numerical quadrature and convergence criteria were used.

Measurement Decomposition


  We first compute CTT reliability coefficients for three types of observed scores to illustrate measurement decompositions. Recall that under a measurement decomposition, we regress an observed score onto all LVs. The observed scores being considered here are: (a) the EAP score for the susceptibility factor (η_i1|_i), (b) the EAP score for the severity factor (η_i2|_i), and (c) the summed score _m_i. The two LVs in the MH-GRM model are the susceptibility LV η_i1 and the severity LV η_i2. The three nonparametric regressions obtained by Step 2 of our MC algorithm are visualized in separate panels of Figure <ref>. For a less dense visualization, we randomly selected a subset of 10^4 MC samples and generated three-dimensional scatter plots of the scores (i.e., the gray dots in Figure <ref>), in which the x- and y-axes represent the latent scores (i.e., the susceptibility and severity LVs) and the z-axis represents the observed score. The fitted values of the nonparametric regressions are depicted by mesh surfaces. The corresponding R^2 statistics, which are estimates of CTT reliability, are presented on top of the graphics.

  The R^2 statistic obtained for the susceptibility EAP score is higher than that of the severity EAP score (i.e., 0.91 versus 0.85, respectively), indicating that the EAP score is more reliable than the CTT score. Moreover, the mesh surface in Figure <ref>A shows that the association between the susceptibility LV and its EAP score is strong, positive, and nonlinear, whereas the association between severity LV and its EAP score is close to zero. In contrast, the severity EAP score is positively related to both susceptibility and severity except when both LVs have low values (see Figure <ref>B). Taken together, the susceptibility EAP score almost exclusively reflects the susceptibility LV. However, the severity EAP score reflects strong relations with the susceptibility and severity LVs. These observations can be gleaned from values of the estimated reliability coefficients and the wire frame plots made available via the MC algorithm.

  The estimated CTT reliability for the summed score is 0.93, which is higher than the two EAP scores. This is expected because the fitted regression surface in Figure <ref>C suggests that both the susceptibility and severity LVs have strong partial effects on the summed score. Stated differently, the high reliability of the summed score is attributed to its strong, positive dependencies on both the susceptibility and severity LVs (relative to the two EAP scores). As a corollary of Theorem 4.4.3 in lord   novick.1968, CTT reliability of the summed score is bounded from below by coefficient alpha under the MH-GRM. In this example, the estimated coefficient alpha is 0.93, which is a very tight lower bound for our MC estimate of 0.93.

Prediction Decomposition



  The two latent scores of interest are η_i1 and η_i2, representing the constructs of susceptbility and severity, respectively. Recall that all the MVs in _i serve as the observed scores in a prediction decomposition. Because the MVs are discrete, the prediction decomposition is a linear regression with 4^14 coefficients in which MVs are dummy coded and we include lower order terms up to the 14-way interaction among the dummy variables. Because it is difficult to estimate and visualize this regression due to the large number of explanatory variables, we computed the response-pattern EAP scores using the  package and applied Step 2' of the MC algorithm. We then plot the latent scores against these EAP scores which are composed of a random subset of 10^4 MC samples. In theory, the regression of a latent score on the associated EAP estimate should be linear with intercept zero and slope one (i.e., a unit-weight linear regression; see Equation <ref>).

  The estimated PRMSE for the susceptibility and severity LVs are 0.88 and 0.72, respectively, implying that the MVs predict the construct of symptom presence better than the construct of symptom frequency under the MH-GRM. Observe also that the PRMSEs for the two LVs (0.88 for susceptibility and 0.72 for severity) are lower than their corresponding CTT reliabilities of the two EAP scores (0.91 for susceptibility and 0.85 for severity), which follows from Equation 31 of kim.2012. The fitted linear regression has a zero intercept and a unit slope, presenting as diagonal straight lines (see Figure <ref>).

  In sum, we presented an empirical example using a more complex measurement model for NCS-R depressive symptom data to illustrate approximate computations of CTT reliability and PRMSE with the proposed MC algorithm. Recall that these reliability coefficients can be obtained from regressing observed scores onto LVs (i.e., measurement decomposition) or regressing latent scores onto MVs (i.e., prediction decomposition). Estimated coefficients of determination (namely, R^2s from the corresponding regressions) approximate CTT reliability and PRMSE. We emphasize that reliability is a property of the observed score under a measurement decomposition, but a property of the latent score under a prediction decomposition. Therefore, different reliability coefficients communicate distinct information and direct comparisons of reliability coefficients from different decompositions are discouraged.

Discussion

  Consistent with the regression framework originated from mcdonald.2011, reliability coefficients are coefficients of determination. Given a population measurement model that specifies the joint distribution of MVs and LVs, we can regress an observed score of interest onto all the LVs, leading to the measurement decomposition of the observed score (Equation <ref>). Alternatively, we can regress a latent score of interest onto all the MVs, leading to the prediction decomposition of the latent score (Equation <ref>). In a measurement decomposition, the predicted value of an observed score amounts to the underlying true score, and the corresponding coefficient of determination coincides with CTT reliability of the observed score <cit.>. In a prediction decomposition, the predicted value of a latent score amounts to its EAP estimate, and the corresponding coefficient of determination coincides with PRMSE of the latent score <cit.>. Hence, we can alternatively view the measurement decomposition as a unit-weight linear regression of the observed score on its true score; similarly, the prediction decomposition can be viewed as a unit-weight linear regression of the latent score on its EAP estimate.

  To compute approximate measures of reliability for computationally intractable situations, we introduced an MC procedure. With a large sample of simulated LVs and MVs, we can fit (nonparametric) regressions predicting either observed or latent scores and approximate the reliability coefficients by R^2 statistics. The MC algorithm allows one to obtain estimates for various types of reliability coefficients without directly evaluating integrals or large finite sums. For the linear one-factor and 2PL examples, we observed that MC-based reliability estimates are nearly identical to their theoretical values. Furthermore, we demonstrated the adaptability of the proposed algorithm to more complex IRT models with an illustration using the NCS-R depressive symptom data.  We also highlighted why different types of reliability coefficients (measurement versus prediction decomposition) should be expected to be different. Thus, when reporting reliability coefficients, it is essential to specify the underlying measurement model and the type of score (i.e., observed versus latent) used as an outcome in the regression framework.

  We have re-introduced and expanded on McDonald's  view on reliability with illustrations using a novel MC approach, leaving a number of related topics to be studied in the future.   First, the performance of the MC procedure should be broadly examined under more extensive simulation settings (e.g., under different measurement models, sample sizes, and model error) as well as a broader set of empirical examples. It is anticipated that the optimal tuning of the MC algorithm (e.g., the number of MC samples M and the choice of nonparametric regression estimators) would differ from case to case. Second, in our empirical example, we treated estimated model parameters as population values, ignoring sampling variability. In practice, however, calibration studies of measurement models tend to have limited sample sizes. Thus, it would be important to quantify sampling variability inherent in estimated reliability coefficients based on our MC algorithm. Finally, to better understand the characteristics of various reliability coefficients, methodologists would need to develop benchmarks or recommendations on how these distinct measures of reliability might be qualitatively interpreted.




