Bi-factor and second-order copula models for item response data
aa r X i v : . [ s t a t . M E ] F e b Bi-factor and second-order copula models for item response data
Sayed H. Kadhem Aristidis K. Nikoloulopoulos * Abstract
Bi-factor and second-order models based on copulas are proposed for item response data, where the itemscan be split into non-overlapping groups such that there is a homogeneous dependence within each group.Our general models include the Gaussian bi-factor and second-order models as special cases and can leadto more probability in the joint upper or lower tail compared with the Gaussian bi-factor and second-ordermodels. Details on maximum likelihood estimation of parameters for the bi-factor and second-order copulamodels are given, as well as model selection and goodness-of-fit techniques. Our general methodology isdemonstrated with an extensive simulation study and illustrated for the Toronto Alexithymia Scale. Ourstudies suggest that there can be a substantial improvement over the Gaussian bi-factor and second-ordermodels both conceptually, as the items can have interpretations of latent maxima/minima or mixtures ofmeans in comparison with latent means, and in fit to data.
Key Words:
Bi-factor model; Conditional independence; Second-order model; Truncated vines; Limitedinformation; Tail dependence/asymmetry.
Factor models are statistical frameworks that analyse multivariate observed variables using few latent variablesor factors (Bartholomew et al., 2011). They are increasingly popular in psychometrics applications, where con-structs such as general well-being could be assessed via survey questions (referred to as items). Datasets withlarge number of items are naturally divided into subgroups, in such, each group of items has homogeneous de-pendence. For example, the well-being (common factor) of patients is usually assessed via items that arise fromseveral sub-domains to assess several group-specific factors such as the depression, anxiety and stress. This spe-cial classification of items is also common in educational assessments and termed as testlets (Wainer and Kiely,1987). It is essential to investigate the items structure, as implementing factor models on testlet-based itemscould result in biased estimates and poor fit (Wang and Wilson, 2005; DeMars, 2006; Zenisky et al., 2002;Sireci et al., 1991; Lee and Frisbie, 1999; Wainer and Thissen, 1996). * Correspondence to: [email protected], Aristidis K. Nikoloulopoulos, School of Computing Sciences, University ofEast Anglia, Norwich NR4 7TJ, U.K.
1o account for the homogeneous dependence in each group of items, Gibbons and Hedeker (1992) andGibbons et al. (2007) proposed bi-factor models for binary and ordinal response data, respectively. The bi-factor models have become omnipresent in analysing survey items that arise from several sub-domains orgroups. They consist of a common factor that is linked to all items, and non-overlapping group-specific factors.The common factor explains dependence between items for all groups, while the group-specific factors explaindependence amongst items within each group. The items are assumed to be independent given the group-specific and common factors.An alternative way of modelling items that are split into several groups is via the second-order model (e.g.,de la Torre and Song 2009; Rijmen 2010), where items are indirectly mapped to an overall (second-order) factorvia non-overlapping group-specific (first-order) factors. Second-order models are suitable when the first-orderfactors are associated with each other, and there is a second-order factor that accounts for the relations amongthe first-order factors.The bi-factor and the second-order models are not generally equivalent (Yung et al., 1999; Gustafsson and Balke,1993; Mulaik and Quartetti, 1997; Rijmen, 2010), unless proportionality constraints are imposed by using theSchmid-Leiman transformation method (Schmid and Leiman, 1957). More importantly, both models are re-stricted to the MVN assumption for the latent variables, which might not be valid. Nikoloulopoulos and Joe(2015) emphasized that if the ordinal variables in item response can be thought of as discretization of latentrandom variables that are maxima/minima or mixtures of means then the use of factor models based on theMVN assumption for the latent variables could provide poor fit.In the context of item response data, latent maxima, minima and means can arise depending on how arespondent considers specific items. An item might make the respondent think about M past events which,say, have values W , . . . , W M . In answering the item, the subject might take the average, maximum or mini-mum of W , . . . , W M and then convert to the ordinal scale depending on the magnitude. The case of a latentmaxima/minima can occur if the response is based on a best or worst case. For different dependent itemsbased on latent maxima or minima, multivariate extreme value and copula theory can be used to select suitabledistributions for the latent variables. Copulas that arise from extreme value theory have more probability inone joint tail (upper or lower) than expected with a MVN distribution where and have latent variables that aremaxima/minima instead of means. Even, in the case where the item responses are based on discretizations oflatent variables that are means, then it is possible that there can be more probability in both the joint upper andjoint lower tail, compared with MVN distributed latent variables. This happens if the respondents consist of a“mixture” population (e.g., different locations or genders). From the theory of elliptical distributions and cop-2las, it is known that the multivariate Student- t distribution as a scale mixture of MVN has more dependencein the tails.Nikoloulopoulos and Joe (2015) have studied factor copula models for item response and have shown thatthere is an improvement on the factor models based on the MVN assumption for the latent variables bothconceptually and in fit to data. This improvement relies on the aforementioned reasons, i.e., items can havemore probability in joint upper or lower tail than would be expected with a MVN or items can be considered asdiscretized maxima/minima or mixtures of discretized means rather than discretized means.In this paper, we propose copula extensions for bi-factor and second-order models. The construction ofthe bi-factor copula model exploits the use of bivariate copulas that link the observed variables to the commonand group-specific factors. Note that if there is only one group of items, then the bi-factor model reducesto the 2-factor copula model in Nikoloulopoulos and Joe (2015). Similarly with the bi-factor copula model,we also use bivariate copulas to construct the second-order copula model. In this case, there are bivariatecopulas that link the observed to the group-specific factors, and also bivariate copulas that link the group-specific to the second-order factor. To account for the dependence between the observed variables and group-specific factors, each group of variables in fact is modelled using the one-factor copula model proposed byNikoloulopoulos and Joe (2015). In addition, if there is only one group of items, then the second-order copulamodel reduces to the one-factor copula model. Hence, the proposed models contain the one- and two-factorcopula models in Nikoloulopoulos and Joe (2015) as special cases, while allowing flexible dependence structurefor both within and between group dependence. As a result, the models are suitable for modelling a high-dimensional item response classified into non-overlapping groups.The proposed models copula models are truncated vine copulas (Brechmann et al., 2012) that involve bothobserved and latent variables. They provide flexible dependence by selecting arbitrary bivariate linking copulas(Joe et al., 2010) to link the items to latent factors. If the bivariate linking copulas are BVN, then the Gaussianbi-factor and second-order models are special cases of our constructions which are the discrete counterparts ofthe structured factor copula models introduced by Krupskii and Joe (2015).The remainder of the paper proceeds as follows. Section 2 introduces the bi-factor and second-order copulamodels for item response and discusses their relationship with the existing models. Estimation techniquesand computational details are provided in Section 3. Section 4 proposes simple diagnostics based on semi-correlations and an heuristic method to select suitable bivariate copulas and build plausible bi-factor and second-order copula models. Section 5 summarizes the assessment of goodness-of-fit of these models using the M statistic of Maydeu-Olivares and Joe (2006), which is based on a quadratic form of the deviations of sample3nd model-based proportions over all bivariate margins. Section 6 contains an extensive simulation study togauge the small-sample efficiency of the proposed estimation, investigate the misspecification of the bivariatecopulas, and examine the reliability of the model selection and goodness-of-fit techniques. Section 7 presents anapplication of our methodology to the Toronto Alexithymia Scale. In this example, it turns out that our models,with linking copulas selected according to the items being latent minima or mixtures of means, provide better fitthan the Gaussian bi-factor and second-order models. We conclude with some discussion in Section 8, followedby a technical Appendix. Let Y , . . . , Y d | {z } , . . . , Y g , . . . , Y d g g | {z } g , . . . , Y G , . . . , Y d G G | {z } G denote the item response variables classified into the G non-overlapping groups. There are d g items in group g ; g = 1 , . . . , G, j = 1 , . . . , d g and collectively thereare d = P Gg =1 d g items, which are all measured on an ordinal scale; Y jg ∈ { , . . . , K jg − } . Let the cutpointsin the uniform U (0 , scale for the jg ’th item be a jg,k , k = 1 , . . . , K − , with a jg, = 0 and a jg,K = 1 .These correspond to a jg,k = Φ( α jg,k ) , where α jg,k are cutpoints in the normal N (0 , scale.The bi-factor and second-order factor copula models are presented in Subsections 2.1 and 2.2, respectively.Subsection 2.3 discusses their relationship with the existing Gaussian bi-factor and second-order models. Consider a common factor X and G group-specific factors X , . . . , X G , where X , X , . . . , X G are indepen-dent and standard uniformly distributed. Let Y jg be the j th observed variable in group g , with y jg being therealization. The bi-factor model assumes that Y g , . . . , Y d g g are conditionally independent given X and X g ,and that Y jg in group g does not depend on X g ′ for g = g ′ . Figure 1 depicts a graphical representation of themodel.The joint probability mass function (pmf) is given by π ( y ) = Pr( Y jg = y jg ; j = 1 , . . . , d g , g = 1 , . . . , G )= Z [0 , G +1 G Y g =1 d g Y j =1 Pr( Y jg = y jg | X = x , X g = x g ) dx · · · dx G dx . According to Sklar’s theorem (Sklar, 1959) there exists a bivariate copula C Y jg ,X such that Pr( Y jg ≤ y jg , X ≤ x ) = C Y jg ,X (cid:0) F Y jg ( y jg ) , x (cid:1) , for x ∈ [0 , , where C Y jg ,X is the copula that links observed variable with4 X X g X G Y · · · Y j · · · Y d Y X Y j X Y d X · · · Y g · · · Y jg · · · Y d g g Y g X Y j g X Y d g g X · · · Y G · · · Y jG · · · Y d G G Y G X Y j G X Y d G G X Y X | X Y j X | X Y d X | X Y g X g | X Y j g X g | X Y d g g X g | X Y G X G | X Y j G X G | X Y d G G X G | X Figure 1: Graphical representation of the bi-factor copula model with G group-specific factors and a common factor X . the common factor X , F Y jg is the cumulative distribution function (cdf) of Y jg ; note that F Y jg is a step functionwith jumps at , . . . , K − , i.e., F Y jg ( y jg ) = a jg,y jg +1 . Then it follows that, F Y jg | X ( y jg | x ) := Pr( Y jg ≤ y jg | X = x ) = ∂∂x C Y jg ,X (cid:0) F Y jg ( y jg ) , x (cid:1) . For shorthand notation, we let C Y jg | X (cid:0) F Y jg ( y jg ) | x (cid:1) = ∂∂x C Y jg ,X (cid:0) F Y jg ( y jg ) , x (cid:1) .The observed variables also load on the group-specific factors, hence to account for this dependence, welet C Y jg ,X g | X be a bivariate copula that links the observed variable Y jg with the group-specific factor X g giventhe common factor X . Hence, Pr( Y jg ≤ y jg | X = x , X g = x g ) = ∂∂x g Pr( Y jg ≤ y jg , X g ≤ x g | X = x )= ∂∂x g C Y jg ,X g | X (cid:0) F Y j g | X ( y jg | x ) , x g (cid:1) = C Y jg | X g ; X (cid:0) F Y j g | X ( y jg | x ) | x g (cid:1) . To this end, the pmf of the bi-factor copula model takes the form π ( y ) = Z [0 , G +1 G Y g =1 d g Y j =1 (cid:26) C Y jg | X g ; X (cid:0) F Y jg | X ( y jg | x ) | x g (cid:1) − C Y jg | X g ; X (cid:0) F Y jg | X ( y jg − | x ) | x g (cid:1)(cid:27) dx · · · dx G dx = Z G Y g =1 ( Z d g Y j =1 (cid:20) C Y jg | X g ; X (cid:0) F Y jg | X ( y jg | x ) | x g (cid:1) − C Y jg | X g ; X (cid:0) F Y jg | X ( y jg − | x ) | x g (cid:1)(cid:21) dx g ) dx = Z G Y g =1 ( Z d g Y j =1 (cid:20) C Y jg | X g ; X (cid:0) C Y jg | X ( a jg,y jg +1 | x ) | x g (cid:1) − C Y jg | X g ; X (cid:0) C Y jg | X ( a jg,y jg | x ) | x g (cid:1)(cid:21) dx g ) dx = Z G Y g =1 ( Z d g Y j =1 f Y jg | X g ; X ( y jg | x g , x ) dx g ) dx . (1)It is shown that the pmf is represented as an one-dimensional integral of a function which is in turn is a productof G one-dimensional integrals. Thus we avoid ( G + 1) -dimensional numerical integration.5or the parametric version of the bi-factor copula model, we let C Y jg ,X and C Y jg ,X g | X be parametriccopulas with dependence parameters θ jg and δ jg , respectively. Assume that for a fixed g = 1 , . . . , G , the items Y g , . . . , Y d g g are conditionally independent given the first-order factors X g ∼ U (0 , , g = 1 , . . . , G and that X = ( X , · · · , X G ) are conditionally independent giventhe second-order factor X ∼ U (0 , . That is the joint distribution of X has an one-factor structure. We alsoassume that Y jg in group g does not depend on X g ′ for g = g ′ . Figure 2 depicts the graphical representation ofthe model. X X X g X G Y · · · Y j · · · Y d · · · Y g · · · Y jg · · · Y d g g · · · Y G · · · Y jG · · · Y d G G X X X X g X X G Y X Y j X Y d X Y g X g Y j g X g Y d g g X g Y G X G Y j G X G Y d G G X G Figure 2: Graphical representation of the second-order copula model with G first-order factors and one second-order factor X . The joint pmf takes the form π ( y ) = Z [0 , G ( G Y g =1 d g Y j =1 Pr( Y jg = y jg | X g = x g ) ) c X ( x , . . . , x G ) dx · · · dx G ; c X is the one-factor copula density (Krupskii and Joe, 2013) of X = ( X , . . . , X G ) , viz. c X ( x , . . . , x G ) = Z G Y g =1 c X g ,X ( x g , x ) dx , where c X g ,X is the bivariate copula density of the copula C X g ,X linking X g and X .Letting C Y jg ,X g be a bivariate copula that joins the observed variable Y jg and the group-specific factor X g such that F Y jg | X g ( y jg | x g ) := Pr( Y jg ≤ y jg | X g = x g ) = ∂∂x g C Y jg ,X g (cid:0) F Y jg ( y jg ) , x g (cid:1) = C Y jg | X g (cid:0) F Y jg ( y jg ) | x g (cid:1) , the pmf of the second-order copula model becomes 6 ( y ) = Z Z [0 , G ( G Y g =1 d g Y j =1 (cid:16) C Y jg | X g (cid:0) F Y jg ( y jg ) | x g (cid:1) − C Y jg | X g (cid:0) F Y jg ( y jg − | x g (cid:1)(cid:17))( G Y g =1 c X g ,X (cid:0) x g , x (cid:1)) dx · · · dx G dx = Z ( G Y g =1 Z " d g Y j =1 (cid:16) C Y jg | X g (cid:0) F Y jg ( y jg ) | x g (cid:1) − C Y jg | X g (cid:0) F Y jg ( y jg − | x g (cid:1)(cid:17) c X g ,X (cid:0) x g , x (cid:1) dx g ) dx = Z ( G Y g =1 Z " d g Y j =1 (cid:16) C Y jg | X g (cid:0) a jg,y jg +1 | x g (cid:1) − C Y jg | X g (cid:0) a jg,y jg | x g (cid:1)(cid:17) c X g ,X (cid:0) x g , x (cid:1) dx g ) dx = Z ( G Y g =1 Z h d g Y j =1 f Y jg | X g ( y jg | x g ) i c X g ,X (cid:0) x g , x (cid:1) dx g ) dx . (2) Similarly with the bi-factor copula model, the pmf is represented as an one-dimensional integral of a functionwhich is in turn is a product of G one-dimensional integrals.For the parametric version of the second-order copula model, we let C Y jg ,X g and C X g ,X be parametriccopulas with dependence parameters θ jg and δ g , respectively. In this subsection we show what happens when all bivariate copulas are BVN. Let Z jg be the underlyingcontinuous variable of the ordinal variable Y jg , i.e., Y jg = y jg if α jg,y jg ≤ Z jg ≤ α jg,y jg +1 with α jg,K = ∞ and α jg, = −∞ .For the bi-factor model, if C Y jg ,X ( · ; θ jg ) and C Y jg ,X g | X ( · ; δ jg ) are BVN copulas, C Y jg | X g ; X ( C Y jg | X ( F jg ( y jg ) | x ) | x g ) = Φ α jg,y jg +1 − θ jg Φ − ( x ) − δ jg q − θ jg Φ − ( x g ) q (1 − θ jg )(1 − δ jg ) . Hence, the pmf for the bi-factor copula model in (1) becomes π ( y ) = Z G Y g =1 (Z d g Y j =1 " Φ α jg,y jg +1 − θ jg Φ − ( x ) − δ jg q − θ jg Φ − ( x g ) q (1 − θ jg )(1 − δ jg ) − Φ α jg,y jg − θ jg Φ − ( x ) − δ jg q − θ jg Φ − ( x g ) q (1 − θ jg )(1 − δ jg ) dx g ) dx = Z ∞−∞ G Y g =1 (Z ∞−∞ d g Y j =1 " Φ α jg,y jg +1 − θ jg z − δ jg q − θ jg z g q (1 − θ jg )(1 − δ jg ) − Φ α jg,y jg − θ jg z − δ jg q − θ jg z g q (1 − θ jg )(1 − δ jg ) φ ( z g ) dz g ) φ ( z ) dz . Z jg = θ jg Z + γ jg Z g + q − θ jg − γ jg ǫ jg , g = 1 , . . . , G, j = 1 , · · · , d g , (3)where γ jg = δ jg q − θ jg and Z , Z g , ǫ jg are iid N (0 , random variables. The parameter θ jg of C Y jg ,X isthe correlation of Z jg and Z , and the parameter δ jg of C Y jg ,X g | X is the partial correlation between Z jg and Z g = Φ − ( X g ) given Z = Φ − ( X ) .It implies that the underlying random variables Z jg ’s have a multivariate Gaussian distribution where theoff-diagonal entries of the correlation matrix have the form θ j g θ j g + γ j g γ j g and θ j g θ j g for j = j and g = g , respectively. For the Gaussian bi-factor model to be identifiable, the number of dependence parametershas to be d − N − N , where N and N is the number of groups that consist of 1 and 2 items, respectively.For a group g of size 1 with variable j , Z g is absorbed with ǫ jg because γ jg would not be identifiable. Fora group g of size 2 with variable indices j , j , the parameters γ j g and γ j g appear only in one correlation,hence one of γ j g , γ j g can be taken as 1 without loss of generality. For the bi-factor copula with non-Gaussianlinking copulas, near non-identifiability can occur when there are groups of size 2; in this case, one of thelinking copulas to the group latent variable can be fixed at comonotonicity.For the Gaussian second-order model let Z , Z ′ , . . . , Z ′ G be the dependent latent N (0 , variables, where Z is the second-order factor and Z ′ g = β g Z +(1 − β g ) Z g is the first-order factor for group g . That is, there is anone second-order factor Z , and the first-order factors Z ′ , . . . , Z ′ G are linear combinations of the second-orderfactor, plus a unique variable Z g for each first-order factor. The stochastic representation is (Krupskii and Joe,2015): Z jg = β jg Z ′ g + q − β jg ǫ jg Z ′ g = β g Z + q − β g Z g , g = 1 , . . . , G, j = 1 , · · · , d g , or Z jg = β jg β g Z + β jg q − β g Z g + q − β jg ǫ jg , j = 1 , · · · , d g . (4)Hence, this is a special case of (3) where θ jg = β jg β g and γ jg = β jg q − β g .8 Estimation and computational details
For the set of all parameters, let θ = ( a , θ g , δ g ) for the bi-factor copula model and θ = ( a , θ g , δ ) for thesecond-order copula model, where a = ( a jg,k : j = 1 , . . . , d g , g = 1 , . . . , G, k = 1 , . . . , K − , θ g =( θ g , . . . , θ jg , . . . , θ d g g : g = 1 , . . . , G ) , δ g = ( δ g , . . . , δ jg , . . . , δ d g g : g = 1 , . . . , G ) and δ = ( δ , . . . , δ G ) .With sample size n and data y , . . . , y n , the joint log-likelihood of the bi-factor and second-order copula is ℓ ( θ ; y , . . . , y n ) = n X i =1 log π ( y i ; θ ) . (5)with π ( y i ; θ ) as in (1) and (2), respectively. Maximization of (5) is numerically possible but is time-consumingfor large d because of many univariate cutpoints and dependence parameters. Hence, we approach estimationusing the two-step IFM method proposed by Joe (2005) that can efficiently, in the sense of computing time andasymptotic variance, estimate the model parameters.In the first step, the cutpoints are estimated using the univariate sample proportions. The univariate cutpointsfor the j th item in group g are estimated as ˆ a jg,k = P ky =0 p jg,y , where p jg,y , y = 0 , . . . , K − for g = 1 , . . . , G and j = 1 , . . . , d g are the univariate sample proportions. In the second step of the IFM method, the joint log-likelihood in (5) is maximized over the copula parameters with the cutpoints fixed as estimated at the first step.The estimated copula parameters can be obtained by using a quasi-Newton (Nash, 1990) method applied to thelogarithm of the joint likelihood.For the bi-factor copula model numerical evaluation of the joint pmf can be achieved with the followingsteps:1. Calculate Gauss-Legendre quadrature (Stroud and Secrest, 1966) points { x q : q = 1 , . . . , n q } andweights { w q : q = 1 , . . . , n q } in terms of standard uniform.2. Numerically evaluate the joint pmf Z G Y g =1 ( Z d g Y j =1 f Y jg | X jg ; X ( y jg | x g , x ) dx g ) dx in a double sum n q X q =1 w q G Y g =1 ( n q X q =1 w q d g Y j =1 f Y jg | X jg ; X ( y jg | x q , x q ) ) For the second-order copula model numerical evaluation of the joint pmf can be achieved with the followingsteps: 9. Calculate Gauss-Legendre quadrature points { x q : q = 1 , . . . , n q } and weights { w q : q = 1 , . . . , n q } interms of stand uniform.2. Numerically evaluate the joint pmf Z ( G Y g =1 Z h d g Y j =1 f Y jg | X g ( y jg | x g ; θ jg ) i c X g ,X (cid:0) x g , x ; δ g (cid:1) dx g ) dx in a double sum n q X q =1 w q ( G Y g =1 n q X q =1 w q h d g Y j =1 f Y jg | X g ( y jg | x q | q ; θ jg ) i) where x q | q = C − Y jg | X g ; X ( x q | x q ; δ g ) . Note that the independent quadrature points { x q : q =1 , . . . , n q } and { x q : q = 1 , . . . , n q } have converted to dependent quadrature points that have anone-factor copula distribution C X ( · ; δ ) .With Gauss-Legendre quadrature, the same nodes and weights are used for different functions; this helps inyielding smooth numerical derivatives for numerical optimization via quasi-Newton. Our comparisons showthat n q = 25 quadrature points are adequate with good precision. In line with Nikoloulopoulos and Joe (2015), we use bivariate parametric copulas that can be used when con-sidering latent maxima, minima or mixtures of means, namely the Gumbel, survival Gumbel (s.Gumbel) andStudent t ν copulas, respectively. A model with bivariate Gumbel copulas that possess upper tail dependencehas latent (ordinal) variables that can be considered as (discretized) maxima, and there is more probability inthe joint upper tail. A model with bivariate s.Gumbel copulas that possess lower tail dependence has latent(ordinal) variables that can be considered as (discretized) minima, and there is more probability in the jointlower tail. A model with bivariate t ν copulas that possess the same lower and upper tail dependence has latent(ordinal) variables that can be considered as mixtures of (discretized) means, since the bivariate Student t ν distribution arises as a scale mixture of bivariate normals. A small value of ν , such as ≤ ν ≤ , leads to amodel with more probabilities in the joint upper and joint lower tails compared with the BVN copula.In the following subsections we describe simple diagnostics based on semi-correlations and an heuristicmethod that automatically selects the bivariate parametric copula families that build either the bi-factor or thesecond-order copula model. In the context of items that can be split into G non-overlapping groups, suchthat there is homogeneous dependence within each group, it is sufficient to (a) summarize the average of the10olychoric semi-correlations for all pairs within each of the G groups and for all pairs of items, and (b) not mixbivariate copulas for a single factor; hence, for both the bi-factor and second-order copula models we allow G + 1 different copula families, one for each group specific factor X g and one for X . Choices of copulas with upper or lower tail dependence are better if the items have more probability in jointlower or upper tail than would be expected with the BVN copula. This can be shown with summaries ofcorrelations in the upper joint tail and lower joint tail.Consider again the underlying N (0 , latent variables Z jg ’s of the ordinal variables Y jg ’s. The correlationsof Z jg ’s in the upper and lower tail, hereafter semi-correlations, are defined as (Joe, 2014, page 71): ρ + N = Cor (cid:16) Z j g , Z j g | Z j g > , Z j g > (cid:17) (6) = R ∞ R ∞ z z φ ( z ) φ ( z ) c (cid:0) Φ( z ) , Φ( z ) (cid:1) dz dz − (cid:18) R ∞ zφ ( z ) (cid:16) − C | (cid:0) . | Φ( z ) (cid:1)(cid:17) dz (cid:19) /C (0 . , . R ∞ z φ ( z ) (cid:16) − C | (cid:0) . | Φ( z ) (cid:1)(cid:17) dz − (cid:18) R ∞ zφ ( z ) (cid:16) − C | (cid:0) . | Φ( z ) (cid:1)(cid:17) dz (cid:19) /C (0 . , .
5) ; ρ − N = Cor (cid:16) Z j g , Z j g | Z j g < , Z j g < (cid:17) = R −∞ R −∞ z z φ ( z ) φ ( z ) c (cid:0) Φ( z ) , Φ( z ) (cid:1) dz dz − (cid:18) R −∞ zφ ( z ) C | (cid:0) . | Φ( z ) (cid:1) dz (cid:19) /C (0 . , . R −∞ z φ ( z ) C | (cid:0) . | Φ( z ) (cid:1) dz − (cid:18) R −∞ zφ ( z ) C | (cid:0) . | Φ( z ) (cid:1) dz (cid:19) /C (0 . , . . From the above expressions, it is clear that the semi-correlations depend only on the copula C of (cid:16) Φ( Z j g ) , Φ( Z j g ) (cid:17) ; C | is the conditional copula cdf. For the BVN and t ν copulas ρ − N = ρ + N , while for the Gum-bel and s.Gumbel copulas ρ − N < ρ + N and ρ − N > ρ + N , respectively. The sample versions of ρ + N , ρ − N foritem response data are the polychoric correlations in the joint lower and upper quadrants of Y j g and Y j g (Kadhem and Nikoloulopoulos, 2021). We propose an heuristic method that selects appropriate bivariate copulas for each factor of the bi-factor andsecond-order copula models. It starts with an initial assumption, that all bivariate linking copulas are BVNcopulas, i.e. the starting model is either the Gaussian bi-factor or second-order model, and then sequentiallyother copulas with lower or upper tail dependence are assigned to the factors where necessary to account formore probability in one or both joint tails. The selection algorithm involves the following steps:1. Fit the bi-factor or second-order copula model with BVN copulas.11. Fit all the possible bi-factor or second-order copula models, iterating over all the copula candidates thatlink all items Y jg ’s in group g or each group-specific factor X g , respectively, to X .3. Select the copula family that corresponds to the lowest Akaike information criterion (AIC), that is, AIC = − × ℓ + 2 × copula parameters.4. Fix the selected copula family that links the observed (bi-factor model) or latent (second-order model)variables to X .5. For g = 1 , . . . , G :(a) Fit all the possible models, iterating over all the copula candidates that link all the items in group g to the group-specific factor X g .(b) Select the copula family that corresponds to the lowest AIC.(c) Fix the selected linking copula family for all the items in group g with X g . We will use the limited information M statistic proposed by Maydeu-Olivares and Joe (2006) to evaluate theoverall fit of the proposed bi-factor and second-order copula models. The M statistic is based on a quadraticform of the deviations of sample and model-based proportions over all bivariate margins. For our parametricmodels with parameter vector θ of dimension q , let π ( θ ) = (cid:0) ˙ π ( θ ) ⊤ , ˙ π ( θ ) ⊤ (cid:1) ⊤ be the column vector ofthe univariate and bivariate model-based marginal probabilities that do not include category 0 with samplecounterpart p = ( ˙ p ⊤ , ˙ p ⊤ ) ⊤ . The total number of the univariate and bivariate residuals (cid:0) p − π (ˆ θ ) (cid:1) ⊤ is s = d ( K −
1) + (cid:18) d (cid:19) ( K − , where d ( K − is the dimension of the univariate residuals and (cid:0) d (cid:1) ( K − is the dimension of the bivariateresiduals excluding category 0.With a sample size n , the limited-information M statistic is given by M = M (ˆ θ ) = n (cid:0) p − π (ˆ θ ) (cid:1) ⊤ C (ˆ θ ) (cid:0) p − π (cid:0) ˆ θ ) (cid:1) , (7)with C ( θ ) = Ξ − − Ξ − ∆ ( ∆ ⊤ Ξ − ∆ ) − ∆ ⊤ Ξ − = ∆ ( c )2 (cid:0) [ ∆ ( c )2 ] ⊤ Ξ ∆ ( c )2 (cid:1) − [ ∆ ( c )2 ] ⊤ , (8)12here ∆ = ∂ π ( θ ) /∂ θ ⊤ is an s × q matrix with the first order derivatives of the univariate and bivariatemarginal probabilities with respect to the estimated model parameters (in the Appendix, we provide detailson the calculation of these derivatives), ∆ ( c )2 is an s × ( s − q ) orthogonal complement to ∆ , such that [ ∆ ( c )2 ] ⊤ ∆ = , and Ξ is the asymptotic s × s covariance matrix of √ n (cid:0) p − π (ˆ θ ) (cid:1) ⊤ . The limited in-formation statistic M under the null hypothesis has an asymptotic distribution that is χ with s − q degrees offreedom when the estimate ˆ θ is √ n -consistent.The asymptotic covariance matrix Ξ can be partitioned according to the portioning of p into Ξ = √ n Acov ( ˙ p ) , Ξ = √ n Acov ( ˙ p , ˙ p ) and Ξ = √ n Acov ( ˙ p ) , where Acov ( · ) denotes asymptotic covari-ance matrix. The elements of Ξ , Ξ and Ξ involve up to the 4-dimensional probabilities as shown below: √ n Acov ( p j ,y , p j ,y ) = π j j ,y y − π j ,y π j ,y √ n Acov ( p j j ,y y , p j ,y ) = π j j j ,y y y − π j j ,y y π j ,y √ n Acov ( p j j ,y y , p j j ,y y ) = π j j j j ,y y y y − π j j ,y y π j j ,y y , where π j,y = Pr( Y j = y ) , π j j ,y y = Pr( Y j = y , Y j = y ) , π j j j ,y y y = Pr( Y j = y , Y j = y , Y j = y ) , and π j j j j ,y y y y = Pr( Y j = y , Y j = y , Y j = y , Y j = y ) . An extensive simulation study is conducted to (a) gauge the small-sample efficiency of the IFM estimationmethod and investigate the misspecification of the bivariate pair-copulas, (b) examine the reliability of usingthe heuristic algorithm to select the true (simulated) bivariate linking copulas, and (c) study the small-sampleperformance of the M statistic.We randomly generate 1,000 datasets with samples of size n = 500 or 1000 and d = 16 items, with K = 3 or K = 5 equally weighted categories, that are equally separated into G = 4 non-overlapping groups from thebi-factor and second-order copula model. In each simulated model, we use different linking copulas to coverdifferent types of dependence. To make the models comparable, we convert the BVN/ t ν and Gumbel/s.Gumbelcopula parameters to Kendall’s τ ’s via τ ( θ ) = 2 π arcsin( θ ) (9)and τ ( θ ) = 1 − θ − , (10)respectively. For the bi-factor copula models we set τ ( θ g ) = (0 . , . , . , . and τ ( δ g ) = (0 . , . , able 1: Small sample of size n = 500 simulations (10 replications) from the bi-factor and second-order factor models with Gumbel copulas and group estimated average biases, root mean square errors(RMSE), and standard deviations (SD), scaled by n , for the IFM estimates under different pair-copulas from the bi-factor and second-order copula models. Bi-factor copula model Second-order copula model τ ( θ g ) , g = 1 , . . . , τ ( δ g ) , g = 1 , . . . , τ ( δ ) τ ( θ g ) , g = 1 , . . . , Fitted model K 0.45 0.55 0.65 0.75 0.30 0.35 0.40 0.50 0.30 0.35 0.40 0.45 0.40 0.50 0.60 0.70 n bias BVN 3 2.65 2.54 2.66 2.16 6.60 7.81 6.99 6.39 5.58 5.34 5.33 5.60 0.41 0.86 0.62 0.275 1.98 2.27 2.54 2.53 5.99 6.27 5.42 2.31 8.71 8.36 7.94 8.52 0.93 0.51 0.58 2.52Gumbel 3 0.39 0.35 0.28 0.34 0.89 1.02 1.62 3.40 -0.18 0.18 0.18 1.88 0.22 0.67 1.14 2.375 0.23 0.23 0.07 0.20 0.84 0.85 1.02 1.98 0.22 0.13 -0.25 1.15 0.23 0.43 0.63 0.60s.Gumbel 3 3.59 3.03 1.51 0.31 4.86 4.52 4.21 1.19 18.43 18.29 18.54 18.68 6.32 6.18 5.47 3.675 0.79 2.25 3.80 5.30 15.89 15.82 13.89 14.52 25.65 24.80 23.58 22.59 3.77 2.54 1.24 2.74 t n SE BVN 3 15.03 13.42 12.37 11.06 30.77 31.20 33.07 39.93 22.80 24.94 24.97 27.03 16.82 16.41 17.06 21.325 13.68 11.89 10.63 8.95 24.58 25.33 25.70 29.86 21.28 23.04 22.45 24.72 15.09 14.27 14.01 15.41Gumbel 3 15.10 13.81 12.33 10.97 29.61 31.34 32.82 42.17 22.58 24.73 25.35 27.87 16.99 16.73 17.66 22.025 13.67 12.29 10.55 8.76 23.60 24.72 25.39 31.13 20.75 22.75 22.69 24.86 15.31 14.62 14.33 15.72s.Gumbel 3 15.58 13.76 12.60 11.27 33.77 34.80 38.18 51.31 25.34 26.80 27.19 29.36 17.40 16.49 16.59 18.465 14.11 12.30 11.16 9.66 27.08 28.44 30.18 40.10 22.61 24.13 23.36 25.46 15.90 14.57 14.38 16.89 t n RMSE BVN 3 15.28 13.66 12.66 11.27 31.48 32.19 33.81 40.45 23.47 25.50 25.53 27.60 16.83 16.44 17.08 21.335 13.83 12.11 10.93 9.30 25.31 26.10 26.27 29.96 22.99 24.51 23.81 26.14 15.12 14.28 14.03 15.62Gumbel 3 15.10 13.81 12.34 10.98 29.63 31.37 32.87 42.31 22.58 24.73 25.35 27.94 16.99 16.75 17.71 22.155 13.67 12.30 10.55 8.77 23.62 24.74 25.42 31.20 20.75 22.75 22.69 24.88 15.31 14.63 14.35 15.73s.Gumbel 3 16.00 14.09 12.69 11.27 34.13 35.13 38.42 51.33 31.33 32.45 32.91 34.80 18.52 17.65 17.49 18.825 14.14 12.51 11.79 11.02 31.41 32.55 33.22 42.67 34.19 34.60 33.19 34.04 16.35 14.82 14.44 17.13 t . , . for g = 1 , . . . , . For the second-order copula models we set τ ( θ g ) = (0 . , . , . , . for g = 1 , . . . , and τ ( δ ) = (0 . , . , . , . .The Kendall’s tau parameters τ ( θ g ) and τ ( δ g ) as described above are common for each group, hence Table1 contains the group estimated average biases, root mean square errors (RMSE), and standard deviations (SD),scaled by n , for the IFM estimates under different pair-copulas from the bi-factor and second-order copulamodels. In the true (simulated) models the linking copulas are Gumbel copulas.Conclusions from the values in the table are the following:• IFM with the true bi-factor or second-order model is highly efficient according to the simulated biases,SDs and RMSEs.• The IFM estimates of τ ’s are not robust under copula misspecification and their biases increase whenthe assumed bivariate copula has tail dependence of opposite direction from the true bivariate copula.For example, in Table 1 the scaled biases for the IFM estimates increase substantially when the linkingcopulas are the s.Gumbel copulas. Table 2: Small sample of size n = 500 simulations ( replications) from the bi-factor and second-order factor models with variouslinking copulas and frequencies of the true bivariate copula identified using the model selection algorithm. Bi-factor Model 1 Model 2 Model 3 Model 4Copula K = 3 K = 5 Copula K = 3 K = 5 Copula K = 3 K = 5 Copula K = 3 K = 5 X Gumbel 992 1000 t
984 1000 Gumbel 996 1000 t
975 1000 X Gumbel 858 956 t
597 806 t
585 789 Gumbel 888 958 X Gumbel 870 951 t
588 799 t
569 775 Gumbel 894 969 X Gumbel 846 950 t
546 777 s.Gumbel 844 945 s.Gumbel 865 947 X Gumbel 844 942 t
589 805 s.Gumbel 878 949 s.Gumbel 900 956Second-order Model 1 Model 2 Model 3 Model 4Copula K = 3 K = 5 Copula K = 3 K = 5 Copula K = 3 K = 5 Copula K = 3 K = 5 X Gumbel 901 848 t
664 819 Gumbel 892 987 t
648 765 X Gumbel 895 975 t
735 939 t
756 933 Gumbel 918 990 X Gumbel 892 962 t
686 911 t
705 910 Gumbel 918 991 X Gumbel 891 981 t
711 915 s.Gumbel 901 980 s.Gumbel 902 982 X Gumbel 900 984 t
743 926 s.Gumbel 904 984 s.Gumbel 919 980
To examine the reliability of using the heuristic algorithm to select the true (simulated) bivariate linkingcopulas, samples of size 500 were generated from various bi-factor and second-order copula models. Table 2presents the number of times that the true (simulated) linking copulas were chosen over 1,000 simulation runs.It is revealed that the model selection algorithm performs extremely well for various bi-factor and second-ordercopulas models with different choices of linking copulas as the number of categories K increases. For a small15 able 3: Small sample of size n = { , } simulations (10 replications) from bi-factor and second-order copula models and theempirical rejection levels at α = { . , . , . , . } , degrees of freedom (df), mean and variance. Bi-factor copula model M Copula n K df Mean Variance α =0.20 α =0.10 α =0.05 α =0.01BVN 500 3 448 449.0 912.8 0.206 0.100 0.060 0.0165 1888 1885.5 4858.3 0.210 0.117 0.065 0.0241000 3 448 448.7 879.0 0.192 0.097 0.051 0.0205 1888 1886.5 4332.5 0.202 0.108 0.064 0.015Gumbel 500 3 448 449.9 887.3 0.216 0.111 0.053 0.0115 1888 1886.6 4709.7 0.225 0.126 0.070 0.0151000 3 448 448.9 864.0 0.201 0.102 0.050 0.0155 1888 1888.6 4332.1 0.226 0.107 0.069 0.014 t
500 3 448 448.7 907.3 0.202 0.088 0.048 0.0185 1888 1886.6 4479.4 0.204 0.107 0.053 0.0171000 3 448 448.6 834.9 0.184 0.090 0.050 0.0145 1888 1890.3 4008.5 0.218 0.103 0.052 0.015Second-order copula model M Copula n K df Mean Variance α =0.20 α =0.10 α =0.05 α =0.01BVN 500 3 460 462.2 1001.2 0.220 0.113 0.055 0.0165 1900 1903.5 3736.2 0.214 0.112 0.052 0.0101000 3 460 461.3 1023.9 0.220 0.109 0.064 0.0135 1900 1906.5 3918.2 0.230 0.130 0.068 0.012Gumbel 500 3 460 464.5 1011.3 0.233 0.117 0.073 0.0245 1900 1909.2 5099.8 0.245 0.129 0.064 0.0081000 3 460 461.9 871.2 0.203 0.106 0.049 0.0095 1900 1908.5 3977.0 0.239 0.129 0.067 0.015 t
500 3 460 465.3 1362.4 0.247 0.145 0.091 0.0395 1900 1904.7 3740.6 0.226 0.113 0.050 0.0101000 3 460 461.8 900.1 0.214 0.108 0.055 0.0105 1900 1908.1 3864.9 0.229 0.131 0.072 0.015 K dependence in the tails cannot be easily quantified. Hence, for example, when the true copula is the t whichhas the same upper and lower tail dependence, the algorithm selected either t or BVN which has zero lowerand upper tail dependence, because both copulas provide reflection symmetric dependence.To check whether the χ s − q is a good approximation for the distribution of the M statistic under the nullhypothesis, samples of sizes 500 and 1000 were generated from various bi-factor second-order copula models.Table 3 contains four common nominal levels of the M statistic under the bi-factor and second-order copulamodels with different bivariate copulas. As can be seen in the table the observed levels of M are close to thenominal α levels and remain accurate even for extremely sparse tables ( d = 16 and K = 5 ).16 Application
The Toronto Alexithymia Scale is the most utilized measure of alexithymia in empirical research (Bagby et al.1994; Gignac et al. 2007; Tuliao et al. 2020) and is composed of d = 20 items that can be subdivided into G = 3 non-overlapping groups: d = 7 items to assess difficulty identifying feelings (DIF), d = 5 items toassess difficulty describing feelings (DDF) and d = 8 items to assess externally oriented thinking (EOT). Weuse a dataset of 1925 university students from the French-speaking region of Belgium (Briganti and Linkowski,2020). Students were 17 to 25 years old and 58% of them were female and 42% were male. They were askedto respond to each item using one of K = 5 categories: “ completely disagree”, “ disagree”, “ neutral, “ agree”, “ completely agree”. The dataset and full description of the items can be found in the R package BGGM (Williams and Mulder, 2020).For these items, a respondent might be thinking about the average “sensation” of many past relevant events,leading to latent means. That is, based on the item descriptions, this seems more natural than a discretizedmaxima or minima. Since the sample is a mixture (male and female students) we can expect a priori that abi-factor or second-order copula model with t ν copulas might be plausible, as in this case the items can beconsidered as mixtures of discretized means. Table 4: Average observed polychoric correlations and semi-correlations for all pairs within each group and for all pairs of items forthe Toronto Alexithymia Scale (TAS), along with the corresponding theoretical semi-correlations for BVN, t , Frank, Gumbel , andsurvival Gumbel (s.Gumbel) copulas. All items Items in group 1 Items in group 2 Items in group 3 ρ N ρ − N ρ + N ρ N ρ − N ρ + N ρ N ρ − N ρ + N ρ N ρ − N ρ + N Observed 0.17 0.21 0.20 0.34 0.36 0.29 0.42 0.37 0.40 0.19 0.26 0.29BVN 0.17 0.07 0.07 0.34 0.16 0.16 0.42 0.21 0.21 0.19 0.08 0.08 t In Table 4 we summarize the averages of polychoric semi-correlations for all pairs within each group andfor all pairs of items along with the theoretical semi-correlations in (6) under different choices of copulas.For a BVN/ t ν copula the copula parameter is the sample polychoric correlation, while for a Gumbel/s.Gumbelcopula the polychoric correlation was converted to Kendall’s tau with the relation in (9) and then from Kendall’s τ to Gumbel/s.Gumbel copula parameter via the functional inverse in (10). The summary of findings from thediagnostics in the table show that• for the first group of items there is more probability in the joint lower tail suggesting s.Gumbel linking17opulas to join each item in this group with the DIF factor;• for the second group of items there is more probability in the joint lower and upper tail suggesting t ν linking copulas to join each item in this group with the DDF factor;• for the third group of items there is more probability in the joint lower and upper tail suggesting t ν linkingcopulas to join each item in this group with the EOT factor;• for the items overall there is more probability in the joint lower and upper tail suggesting t ν linkingcopulas to join each item or group specific factor (second-order model) with the common factor.Hence, a bi-factor or second-order copula model with the aforementioned linking copulas might provide a betterfit that the (Gaussian) models with BVN copulas.Then, we fit the bi-factor and second-order models with the copulas selected by the heuristic algorithmin Section 4.2. For a baseline comparison, we also fit their special cases; these are the one- and two-factorcopula models where we have also selected the bivariate copulas using the heuristic algorithm proposed byKadhem and Nikoloulopoulos (2021). To show the improvement of the copula models over their Gaussiananalogues, we have also fitted all the classes of copula models with BVN copulas. The fitted models arecompared via the AIC, since the number of parameters is not the same between the models. In addition, weuse the Vuong’s test (Vuong, 1989) to show if (a) the best fitted model according to the AICs provides better fitthan the other fitted models and (b) a model with the selected copulas provides better fit than the one with BVNcopulas. The Vuong’s test is the sample version of the difference in Kullback-Leibler divergence between twomodels and can be used to differentiate two parametric models which could be non-nested. For the Vuong’s testwe provide the 95% confidence interval of the Vuong’s test statistic (Joe, 2014, page 258). If the interval doesnot contain 0, then the best fitted model according to the AICs is better if the interval is completely above 0.To assess the overall goodness-of-fit of the bi-factor and second-order copula models, we use the M statistic(Maydeu-Olivares and Joe, 2006).Table 5 gives the AICs, the 95% CIs of Vuong’s tests and the M statistics for all the fitted models. Thebest fitted model, based on AIC values, is the bi-factor copula model obtained from the selection algorithm.The best fitted bi-factor copula model results when we use s.Gumbel for the DIF factor, t for both the DDFand EOT factors and t for the common factor (alexithymia). This is in line with the preliminary analyses basedon the interpretations of items as mixtures of means and the diagnostics in Table 4. It is revealed that the DIF18 able 5: AICs, Vuong’s 95% CIs, and M statistics for the 1-factor, 2-factor, bi-factor and second-order copula models with BVNcopulas and selected copulas, along with the maximum deviations of observed and expected counts for all pairs within each group andfor all pairs of items for the Toronto Alexithymia Scale. a (0.35,0.50) (0.53,0.69) (0.51,0.69) (0.38,0.52)Vuong’s 95% CI b (0.93,1.13) (0.55,0.67) (0.69,0.88) (0.13,0.23) (0.51,0.69) (0.61,0.80) (0.21,0.29) M p -value < . < . < . < . < . < . < . < . Maximum discrepancyItems in Group 1 71 63 71 60 69 55 70 61Items in Group 2 112 98 113 83 77 48 84 55Items in Group 3 87 74 81 52 80 45 82 53All items 112 98 113 83 80 55 84 61 a Selected factor copula model versus its Gaussian special case. b Selected Bi-factor copula model versus any other fitted model. items and DIF factor are discretized and latent minima, respectively, as the participants seem to reflect that they“disagree” or “completely disagree” having difficulty identifying feelings. From the Vuong’s 95% Cls and M statistics it is shown that factor copula models provide a big improvement over their Gaussian analogues andthat the selected bi-factor copula model outperforms all the fitted models.Although the selected bi-factor copula model shows substantial improvement over the Gaussian bi-factormodel or any other fitted model, it is not so clear from the goodness-of-fit p -values that the response patterns aresatisfactorily explained by using the linking copulas selected by the heuristic algorithm. This is not surprisingsince one should expect discrepancies between the postulated parametric model and the population probabili-ties, when the sample size or dimension is sufficiently large (Maydeu-Olivares and Joe, 2014). To further showthat the fit has been improved we have calculated the maximum deviations of observed and model-based countsfor each bivariate margin, that is, D j j = n max y ,y | p j ,j ,y ,y − π j ,j ,y ,y (ˆ θ ) | . In Table 5 we summarizethe averages of these deviations for all pairs within each group and for all pairs of items. Overall, the maximumdiscrepancies have been sufficiently reduced in the selected bi-factor model.Table 6 gives the copula parameter estimates in Kendall’s τ scale and their standard errors (SE) for theselected bi-factor copula model and the Gaussian bi-factor model as the benchmark model. The SEs of theestimated parameters are obtained by the inversion of the Hessian matrix at the second step of the IFM method.These SEs are adequate to assess the flatness of the log-likelihood. Proper SEs that account for the estimationof cutpoints can be obtained by jackknifing the two-stage estimation procedure. The loading parameters ( ˆ τ ’s19 able 6: Estimated copula parameters and their standard errors (SE) in Kendall’s τ scale for the Bi-factor copula models with BVNcopulas and selected copulas for the Toronto Alexithymia Scale. Bi-factor copula model with BVN copulas Bi-factor copula model with selected copulasCommon factor Group-specific factors Common factor Group-specific factorsItems Est SE Est SE Copulas Est SE Copulas Est SE1 0.42 0.01 0.23 0.02 t t t t t t t t t t t t t t t t t t t t t t t t t t t t t t t t t For item response data that can be split into non-overlapping groups, we have proposed bi-factor and second-order copula models where we replace BVN distributions, between observed and latent variables, with bivariatecopulas. Our copula constructions include the Gaussian bi-factor and second-order models as special cases andcan provide a substantial improvement over the Gaussian models based on AIC, Vuong’s and goodness-of-fitstatistics. Hence, superior statistical inference for the loading parameters of interest can be achieved. Theimprovement relies on the fact that when we use appropriate bivariate copulas other than BVN copulas in theconstruction, there is an interpretation of latent variables that can be maxima/minima or mixture of meansinstead of means. 20ur constructions have a latent structure that is not additive as in (3) and (4) if other than BVN copulasare called and the bi-factor copula (dependence) parameters are interpretable as dependence of an observedvariable with the common factor, or conditional dependence of an observed variable with the group-specificlatent variable given the common factor.We have proposed a fast and efficient likelihood estimation technique based on Gauss-Legendre quadraturepoints. The joint pmfs in (1) and (2) reduce to one-dimensional integrals of a function which in turn is a productof G one-dimensional integrals. Hence, the evaluation of the the joint likelihood requires only low-dimensionalintegration regardless of the dimension of the latent variables.Building on the models proposed in this paper, there are several extensions that can be implemented. Theadoption of the structure of the Gaussian tri-factor and the third-order models (e.g., Rijmen et al. 2014), toaccount for any additional layer of dependence, is feasible using the notion of truncated vine copulas thatinvolve both observed and latent variables. Software R functions for estimation, model selection and goodness-of-fit of the bi-factor and second-order copula modelswill be part of the next major release of the R package FactorCopula (Kadhem and Nikoloulopoulos, 2020).
Acknowledgements
The simulations presented in this paper were carried out on the High Performance Computing Cluster supportedby the Research and Specialist Computing Support service at the University of East Anglia.
Appendix
We provide the form of the derivatives of the univariate and bivariate marginal probabilities with respect to theestimated model parameters in the Appendix Tables 1–5.
Appendix Table 1: Derivatives of the univariate probability π jg,y = Φ( α jg,y +1 ) − Φ( α jg,y ) with respect to the cutpoint α jg,k for g = 1 . . . , G, j = 1 , . . . , d g , y = 1 , . . . , K − , and k = 1 , . . . , K − . ∂π jg,y /∂α jg,k If φ ( α jg,y +1 ) k = y + 1 − φ ( α jg,y ) k = y ppendix Table 2: Derivatives of the bivariate probability π j j g,y ,y = Pr( Y j g = y , Y j g = y ) with respect to the cutpoint α jg,k , the copula parameter θ jg for thecommon factor X , and the copula parameter δ jg for the group-specific factor X g for the bi-factor copula model for g = 1 . . . , G, j, j , j = 1 , . . . , d g , y, y , y =1 , . . . , K − , and k = 1 , . . . , K − . Note that f Y jg | X g ; X ( y jg | x g ; x ) = (cid:16) C Y jg | X g ; X (cid:0) C Y jg | X ( a jg,y +1 | x ; θ jg ) | x g ; δ jg (cid:1) − C Y jg | X g ; X (cid:0) C Y jg | X ( a jg,y | x ; θ jg ) | x g ; δ jg (cid:1)(cid:17) where a jg,k = Φ( α jg,k ) , c X Y jg ( x , a ) = ∂ C X Y jg ( x , a ) /∂x ∂a , ˙ C jg | X ( · ; θ jg ) = ∂C jg | X ( · ; θ jg ) /∂θ jg , ˙ C Y jg | X g ; X ( · ; δ jg ) = ∂C Y jg | X g ; X ( · ; δ jg ) /∂δ jg , ˙ f Y jg | X jg ; X ( y jg | x g ; x ) = ∂f Y jg | X jg ; X ( y jg | x g ; x ) /∂δ jg = ˙ C Y jg | X g ; X (cid:0) C Y jg | X ( a jg,y +1 | x ) | x g (cid:1) − ˙ C Y jg | X g ; X (cid:0) C Y jg | X ( a jg,y | x ) | x g (cid:1) , ¯ f Y jg | X jg ; X ( y jg | x g ; x ) = ∂f Y jg | X jg ; X ( y jg | x g ; x ) /∂θ jg = c X g Y jg (cid:0) x g , C Y jg | X ( a jg,y +1 | x ) (cid:1) ˙ C Y jg | X ( a jg,y +1 | x ) − c X g Y jg (cid:0) x g , C Y jg | X ( a jg,y | x ) (cid:1) ˙ C Y jg | X ( a jg,y | x ) . ∂π j j g,y ,y /∂α jg,k If φ ( α j g,y +1 ) R R f Y j g | X g ; X ( y j g | x g ; x ) c X g Y j g (cid:0) x g , C Y j g | X ( a j g,y +1 | x ) (cid:1) c X Y j g ( x , a j g,y +1 ) dx g dx j = j , k = y + 1 − φ ( α j g,y ) R R f Y j g | X g ; X ( y j g | x g ; x ) c X g Y j g (cid:0) x g , C Y j g | X ( a j g,y | x ) (cid:1) c X Y j g ( x , a j g,y ) dx g dx j = j , k = y φ ( α j g,y +1 ) R R f Y j g | X g ; X ( y j g | x g ; x ) c X g Y j g (cid:0) x g , C Y j g | X ( a j g,y +1 | x ) (cid:1) c X Y j g ( x , a j g,y +1 ) dx g dx j = j , k = y + 1 − φ ( α j g,y ) R R f Y j g | X g ; X ( y j g | x g ; x ) c X g Y j g (cid:0) x g , C Y j g | X ( a j g,y | x ) (cid:1) c X Y j g ( x , a j g,y ) dx g dx j = j , k = y ∂π j j g,y ,y /∂θ jg If R R f Y j g | X g ; X ( y j g | x g ; x ) ¯ f Y j g | X j g ; X ( y j g | x g ; x ) dx g dx j = j R R f Y j g | X g ; X ( y j g | x g ; x ) ¯ f Y j g | X j g ; X ( y j g | x g ; x ) dx g dx j = j ∂π j j g,y ,y /∂δ jg If R R f Y j g | X g ; X ( y j g | x g ; x ) ˙ f Y j g | X j g ; X ( y j g | x g ; x ) dx g dx j = j R R f Y j g | X g ; X ( y j g | x g ; x ) ˙ f Y j g | X j g ; X ( y j g | x g ; x ) dx g dx j = j ppendix Table 3: Derivatives of the bivariate probability π j g j g ,y ,y = Pr( Y j g = y , Y j g = y ) with respect to the cutpoint α jg,k , the copula parameter θ jg forthe common factor X , and the copula parameter δ jg for the group-specific factor X g for the bi-factor copula model for g = 1 . . . , G, j, j , j = 1 , . . . , d g , y, y , y =1 , . . . , K − , and k = 1 , . . . , K − . Note that f Y jg | X g ; X ( y jg | x g ; x ) = (cid:16) C Y jg | X g ; X (cid:0) C Y jg | X ( a jg,y +1 | x ; θ jg ) | x g ; δ jg (cid:1) − C Y jg | X g ; X (cid:0) C Y jg | X ( a jg,y | x ; θ jg ) | x g ; δ jg (cid:1)(cid:17) where a jg,k = Φ( α jg,k ) , c X Y jg ( x , a ) = ∂ C X Y jg ( x , a ) /∂x ∂a , ˙ C jg | X ( · ; θ jg ) = ∂C jg | X ( · ; θ jg ) /∂θ jg , ˙ C Y jg | X g ; X ( · ; δ jg ) = ∂C Y jg | X g ; X ( · ; δ jg ) /∂δ jg , ˙ f Y jg | X jg ; X ( y jg | x g ; x ) = ∂f Y jg | X jg ; X ( y jg | x g ; x ) /∂δ jg = ˙ C Y jg | X g ; X (cid:0) C Y jg | X ( a jg,y +1 | x ) | x g (cid:1) − ˙ C Y jg | X g ; X (cid:0) C Y jg | X ( a jg,y | x ) | x g (cid:1) , ¯ f Y jg | X jg ; X ( y jg | x g ; x ) = ∂f Y jg | X jg ; X ( y jg | x g ; x ) /∂θ jg = c X g Y jg (cid:0) x g , C Y jg | X ( a jg,y +1 | x ) (cid:1) ˙ C Y jg | X ( a jg,y +1 | x ) − c X g Y jg (cid:0) x g , C Y jg | X ( a jg,y | x ) (cid:1) ˙ C Y jg | X ( a jg,y | x ) . ∂π j g j g ,y ,y /∂α jg,k If φ ( α j g ,y +1 ) R R f Y j g | X g ; X ( y j g | x g ; x ) dx g R c X g Y j g (cid:0) x g , C Y j g | X ( a j g ,y +1 | x ) (cid:1) c X Y j g ( x , a j g ,y +1 ) dx g dx j = j , g = g , k = y + 1 − φ ( α j g ,y ) R R f Y j g | X g ; X ( y j g | x g ; x ) dx g R c X g Y j g (cid:0) x g , C Y j g | X ( a j g ,y | x ) (cid:1) c X Y j g ( x , a j g ,y ) dx g dx j = j , g = g , k = y φ ( α j g ,y +1 ) R R f Y j g | X g ; X ( y j g | x g ; x ) dx g R c X g Y j g (cid:0) x g , C Y j g | X ( a j g ,y +1 | x ) (cid:1) c X Y j g ( x , a j g ,y +1 ) dx g dx j = j , g = g , k = y + 1 − φ ( α j g ,y ) R R f Y j g | X g ; X ( y j g | x g ; x ) dx g R c X g Y j g (cid:0) x g , C Y j g | X ( a j g ,y | x ) (cid:1) c X Y j g ( x , a j g ,y ) dx g dx j = j , g = g , k = y ∂π j g j g ,y ,y /∂θ jg If R R f Y j g | X g ; X ( y j g | x g ; x ) dx g R ¯ f Y j g | X j g ; X ( y j g | x g ; x ) dx g dx j = j , g = g R R f Y j g | X g ; X ( y j g | x g ; x ) dx g R ¯ f Y j g | X j g ; X ( y j g | x g ; x ) dx g dx j = j , g = g ∂π j g j g ,y ,y /∂δ jg If R R f Y j g | X g ; X ( y j g | x g ; x ) dx g R ˙ f Y j g | X j g ; X ( y j g | x g ; x ) dx g dx j = j , g = g R R f Y j g | X g ; X ( y j g | x g ; x ) dx g R ˙ f Y j g | X j g ; X ( y j g | x g ; x ) dx g dx j = j , g = g ppendix Table 4: Derivatives of the bivariate probabilities π j j g,y ,y = Pr( Y j g = y , Y j g = y ) with respect to the cutpoint α jg,k , the copula parameter θ jg for the first-order factor X g ,and the copula parameter δ g for the the second-order factor X for the second-order copula model for g = 1 . . . , G, j, j , j = 1 , . . . , d g , y, y , y = 1 , . . . , K − , and k = 1 , . . . , K − .Note that f Y jg | X g ( y jg | x g ) = C Y jg | X g (cid:0) a jg,y +1 | x g ; θ jg (cid:1) − C Y jg | X g (cid:0) a jg,y | x g ; θ jg (cid:1) , c X g Y jg ( x g , a ) = ∂ C X g Y jg ( x g , a ) /∂x g ∂a, ˙ C Y jg | X g ( · ; θ jg ) = ∂C Y jg | X g ( · ; θ jg ) /∂θ jg , ˙ f Y jg | X jg ( y jg | x g ) = ∂f Y jg | X jg ( y jg | x g ) /∂θ jg = ˙ C Y jg | X g (cid:0) a jg,y +1 | x g (cid:1) − ˙ C Y jg | X g (cid:0) a jg,y | x g (cid:1) , ˙ c X g X ( x g , x ; δ g ) = ∂c X g X ( x g , x ; δ g ) /∂δ g . ∂π j j g,y ,y /∂α jg,k If φ ( α j g,y +1 ) R R f Y j g | X g ( y j g | x g ) c X g Y j g ( x g , a j g,y +1 ) c X g X ( x g , x ) dx g dx j = j , k = y + 1 − φ ( α j g,y ) R R f Y j g | X g ( y j g | x g ) c X g Y j g ( x g , a j g,y ) c X g X ( x g , x ) dx g dx j = j , k = y φ ( α j g,y +1 ) R R f Y j g | X g ( y j g | x g ) c X g Y j g ( x g , a j g,y +1 ) c X g X ( x g , x ) dx g dx j = j , k = y + 1 − φ ( α j g,y ) R R f Y j g | X g ( y j g | x g ) c X g Y j g ( x g , a j g,y ) c X g X ( x g , x ) dx g dx j = j , k = y ∂π j j g,y ,y /∂θ jg If R R f Y j g | X g ( y j g | x g ) ˙ f Y j g | X g ( y j g | x g ) c X g X ( x g , x ) dx g dx j = j R R f Y j g | X g ( y j g | x g ) ˙ f Y j g | X g ( y j g | x g ) c X g X ( x g , x ) dx g dx j = j ∂π j j g,y ,y /∂δ g R R f Y j g | X g ( y j g | x g ) f Y j g | X g ( y j g | x g ) ˙ c X g X ( x g , x ) dx g dx ppendix Table 5: Derivatives of the bivariate probability π j g j g ,y ,y = Pr( Y j g = y , Y j g = y ) with respect to the cutpoint α jg,k , the copula parameter θ jg for the first-order factor X g , and the copula parameter δ g for the second-order factor X for the second-order copula model for g = 1 . . . , G, j, j , j = 1 , . . . , d g , y, y , y = 1 , . . . , K − , and k = 1 , . . . , K − .Note that f Y jg | X g ( y jg | x g ) = C Y jg | X g (cid:0) a jg,y +1 | x g ; θ jg (cid:1) − C Y jg | X g (cid:0) a jg,y | x g ; θ jg (cid:1) , c X g Y jg ( x g , a ) = ∂ C X g Y jg ( x g , a ) /∂x g ∂a, ˙ C Y jg | X g ( · ; θ jg ) = ∂C Y jg | X g ( · ; θ jg ) /∂θ jg , ˙ f Y jg | X jg ( y jg | x g ) = ∂f Y jg | X jg ( y jg | x g ) /∂θ jg = ˙ C Y jg | X g (cid:0) a jg,y +1 | x g (cid:1) − ˙ C Y jg | X g (cid:0) a jg,y | x g (cid:1) , ˙ c X g X ( x g , x ; δ g ) = ∂c X g X ( x g , x ; δ g ) /∂δ g . ∂π j g j g ,y ,y /∂α jg,k If φ ( α j g ,y +1 ) R R f Y j g | X g ( y j g | x g ) c X g X ( x g , x ) dx g R c X g Y j g ( x g , a j g ,y +1 ) c X g X ( x g , x ) dx g dx j = j , g = g , k = y + 1 − φ ( α j g ,y ) R R f Y j g | X g ( y j g | x g ) c X g X ( x g , x ) dx g R c X g Y j g ( x g , a j g ,y ) c X g X ( x g , x ) dx g dx j = j , g = g , k = y φ ( α j g ,y +1 ) R R f Y j g | X g ( y j g | x g ) c X g X ( x g , x ) dx g R c X g Y j g ( x g , a j g ,y +1 ) c X g X ( x g , x ) dx g dx j = j , g = g , k = y + 1 − φ ( α j g ,y ) R R f Y j g | X g ( y j g | x g ) c X g X ( x g , x ) dx g R c X g Y j g ( x g , a j g ,y ) c X g X ( x g , x ) dx g dx j = j , g = g , k = y ∂π j g j g ,y ,y /∂θ jg If R R f Y j g | X g ( y j g | x g ) c X g X ( x g , x ) dx g R ˙ f Y j g | X g ( y j g | x g ) c X g X ( x g , x ) dx g dx j = j , g = g R R f Y j g | X g ( y j g | x g ) c X g X ( x g , x ) dx g R ˙ f Y j g | X g ( y j g | x g ) c X g X ( x g , x ) dx g dx j = j , g = g ∂π j g j g ,y ,y /∂δ g If R R f Y j g | X g ( y j g | x g ) c X g X ( x g , x ) dx g R f Y j g | X g ( y j g | x g ) ˙ c X g X ( x g , x ) dx g dx g = g R R f Y j g | X g ( y j g | x g ) c X g X ( x g , x ) dx g R f Y j g | X g ( y j g | x g ) ˙ c X g X ( x g , x ) dx g dx g = g eferences Bagby, R., Parker, J. D., and Taylor, G. J. (1994). The twenty-item toronto alexithymia scale–i. item selectionand cross-validation of the factor structure.
Journal of Psychosomatic Research , 38(1):23–32.Bartholomew, D. J., Knott, M., and Moustaki, I. (2011).
Latent Variable Models and Factor Analysis: A UnifiedApproach . Wiley Series in Probability and Statistics. Wiley.Brechmann, E. C., Czado, C., and Aas, K. (2012). Truncated regular vines in high dimensions with applicationsto financial data.
Canadian Journal of Statistics , 40(1):68–85.Briganti, G. and Linkowski, P. (2020). Network approach to items and domains from the toronto alexithymiascale.
Psychological Reports , 123(5):2038–2052.de la Torre, J. and Song, H. (2009). Simultaneous estimation of overall and domain abilities: A higher-orderIRT model approach.
Applied Psychological Measurement , 33(8):620–639.DeMars, C. E. (2006). Application of the bi-factor multidimensional item response theory model to testlet-based tests.
Journal of Educational Measurement , 43(2):145–168.Gibbons, R. D., Bock, R. D., Hedeker, D., Weiss, D. J., Segawa, E., Bhaumik, D. K., Kupfer, D. J., Frank, E.,Grochocinski, V. J., and Stover, A. (2007). Full-information item bifactor analysis of graded response data.
Applied Psychological Measurement , 31(1):4–19.Gibbons, R. D. and Hedeker, D. R. (1992). Full-information item bi-factor analysis.
Psychometrika , 57(3):423–436.Gignac, G. E., Palmer, B. R., and Stough, C. (2007). A confirmatory factor analytic investigation of the tas–20:Corroboration of a five-factor model and suggestions for improvement.
Journal of Personality Assessment ,89(3):247–257.Gustafsson, J.-E. and Balke, G. (1993). General and specific abilities as predictors of school achievement.
Multivariate Behavioral Research , 28(4):407–434.Joe, H. (2005). Asymptotic efficiency of the two-stage estimation method for copula-based models.
Journal ofMultivariate Analysis , 94:401–419. 26oe, H. (2014).
Dependence Modelling with Copulas . Chapman and Hall/CRC.Joe, H., Li, H., and Nikoloulopoulos, A. K. (2010). Tail dependence functions and vine copulas.
Journal ofMultivariate Analysis , 101(1):252–270.Kadhem, S. H. and Nikoloulopoulos, A. K. (2020).
FactorCopula: Factor copula models for mixed data . Rpackage version 0.5. URL: http://CRAN.R-project.org/package=FactorCopula.Kadhem, S. H. and Nikoloulopoulos, A. K. (2021). Factor copula models for mixed data.
Bristish Journal ofMathematical and Statsitical Psycology . DOI:10.1111/bmsp.12231.Krupskii, P. and Joe, H. (2013). Factor copula models for multivariate data.
Journal of Multivariate Analysis ,120:85–101.Krupskii, P. and Joe, H. (2015). Structured factor copula models: Theory, inference and computation.
Journalof Multivariate Analysis , 138:53–73.Lee, G. and Frisbie, D. A. (1999). Estimating reliability under a generalizability theory model for test scorescomposed of testlets.
Applied Measurement in Education , 12(3):237–255.Maydeu-Olivares, A. and Joe, H. (2006). Limited information goodness-of-fit testing in multidimensionalcontingency tables.
Psychometrika , 71(4):713–732.Maydeu-Olivares, A. and Joe, H. (2014). Assessing approximate fit in categorical data analysis.
MultivariateBehavioral Research , 49(4):305–328.Mulaik, S. A. and Quartetti, D. A. (1997). First order or higher order general factor?
Structural EquationModeling: A Multidisciplinary Journal , 4(3):193–211.Nash, J. (1990).
Compact Numerical Methods for Computers: Linear Algebra and Function Minimisation .Hilger, New York. 2nd edition.Nikoloulopoulos, A. K. and Joe, H. (2015). Factor copula models for item response data.
Psychometrika ,80(1):126–150.Rijmen, F. (2010). Formal relations and an empirical comparison among the bi-factor, the testlet, and a second-order multidimensional irt model.
Journal of Educational Measurement , 47(3):361–372.27ijmen, F., Jeon, M., von Davier, M., and Rabe-Hesketh, S. (2014). A third-order item response theory modelfor modeling the effects of domains and subdomains in large-scale educational assessment surveys.
Journalof Educational and Behavioral Statistics , 39(4):235–256.Schmid, J. and Leiman, J. M. (1957). The development of hierarchical factor solutions.
Psychometrika ,22(1):53–61.Sireci, S. G., Thissen, D., and Wainer, H. (1991). On the reliability of testlet-based tests.
Journal of EducationalMeasurement , 28(3):237–247.Sklar, A. (1959). Fonctions de r´epartition `a n dimensions et leurs marges. Publications de l’Institut de Statis-tique de l’Universit´e de Paris , 8:229–231.Stroud, A. and Secrest, D. (1966).
Gaussian Quadrature Formulas . Prentice-Hall, Englewood Cliffs, NJ.Tuliao, A. P., Klanecky, A. K., Landoy, B. V. N., and McChargue, D. E. (2020). Toronto alexithymia scale–20:Examining 18 competing factor structure solutions in a u.s. sample and a philippines sample.
Assessment ,27(7):1515–1531.Vuong, Q. H. (1989). Likelihood ratio tests for model selection and non-nested hypotheses.
Econometrica ,57(2):307–333.Wainer, H. and Kiely, G. L. (1987). Item clusters and computerized adaptive testing: A case for testlets.
Journalof Educational Measurement , 24(3):185–201.Wainer, H. and Thissen, D. (1996). How is reliability related to the quality of test scores? what is the effect oflocal dependence on reliability?
Educational Measurement: Issues and Practice , 15(1):22–29.Wang, W.-C. and Wilson, M. (2005). The rasch testlet model.
Applied Psychological Measurement , 29(2):126–149.Williams, D. and Mulder, J. (2020).
BGGM: Bayesian Gaussian Graphical Models . R package version 1.0.0.Yung, Y.-F., Thissen, D., and McLeod, L. D. (1999). On the relationship between the higher-order factor modeland the hierarchical factor model.
Psychometrika , 64(2):113–128.Zenisky, A. L., Hambleton, R. K., and Sireci, S. G. (2002). Identification and evaluation of local item depen-dencies in the medical college admissions test.