Bi-factor and Second-Order Copula Models for Item Response Data

Sayed H. Kadhem; Aristidis K. Nikoloulopoulos

doi:10.1007/s11336-022-09894-2

Bi-factor and Second-Order Copula Models for Item Response Data

Published online by Cambridge University Press: 01 January 2025

Sayed H. Kadhem and

Aristidis K. Nikoloulopoulos

Show author details

Sayed H. Kadhem: Affiliation:
University of East Anglia
Aristidis K. Nikoloulopoulos*: Affiliation:
University of East Anglia
*: Correspondence should be made to Aristidis K. Nikoloulopoulos, School of Computing Sciences, University of East Anglia, Norwich NR4 7TJ, UK. Email: a.nikoloulopoulos@uea.ac.uk

Article contents

Abstract
Bi-factor and Second-Order Copula Models
Estimation and Computational Details
Bivariate Copula Selection
Goodness of Fit
Simulations
Application
Discussion
Software
Footnotes
References

Rights & Permissions

Abstract

Bi-factor and second-order models based on copulas are proposed for item response data, where the items are sampled from identified subdomains of some larger domain such that there is a homogeneous dependence within each domain. Our general models include the Gaussian bi-factor and second-order models as special cases and can lead to more probability in the joint upper or lower tail compared with the Gaussian bi-factor and second-order models. Details on maximum likelihood estimation of parameters for the bi-factor and second-order copula models are given, as well as model selection and goodness-of-fit techniques. Our general methodology is demonstrated with an extensive simulation study and illustrated for the Toronto Alexithymia Scale. Our studies suggest that there can be a substantial improvement over the Gaussian bi-factor and second-order models both conceptually, as the items can have interpretations of discretized maxima/minima or mixtures of discretized means in comparison with discretized means, and in fit to data.

Keywords

Bi-factor model conditional independence limited information second-order model tail dependence/asymmetry truncated vines

Information

Type: Theory and Methods
Information: Psychometrika , Volume 88 , Issue 1 , March 2023 , pp. 132 - 157

DOI: https://doi.org/10.1007/s11336-022-09894-2 [Opens in a new window]
Creative Commons: This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Copyright: Copyright © 2022 The Author(s) under exclusive licence to The Psychometric Society

Psychological scales and educational tests are developed to measure a particular construct by selecting items from several identified domains (Gibbons et al. Reference Gibbons, Bock, Hedeker, Weiss, Segawa, Bhaumik, Kupfer, Frank, Grochocinski and Stover2007). For example, a questionnaire or instrument, used in psychometrics to assess abstract concepts, such as the well-being has a large number of items or questions that are sampled from several subdomains such as depression, anxiety and stress. This special classification of items in educational assessments is termed ‘testlets’ (Wainer and Kiely Reference Wainer and Kiely1987). It is essential to investigate the factorial structure, as implementing unstructured factor models on testlet-based items could result in biased estimates and a poor fit (Wang and Wilson Reference Wang and Wilson2005; DeMars Reference DeMars2006; Zenisky et al. Reference Zenisky, Hambleton and Sireci2002; Sireci et al. Reference Sireci, Thissen and Wainer1991; Lee and Frisbie Reference Lee and Frisbie1999; Wainer and Thissen Reference Wainer and Thissen1996).

To account for the homogeneous dependence in several subdomains of some larger domain, Gibbons and Hedeker (Reference Gibbons and Hedeker1992) and Gibbons et al. (Reference Gibbons, Bock, Hedeker, Weiss, Segawa, Bhaumik, Kupfer, Frank, Grochocinski and Stover2007) proposed bi-factor models for binary and ordinal items, respectively. They consist of a common factor, that is linked to all items, and non-overlapping group-specific factors. The common factor explains the dependence between all the items, while the group-specific factors explain the dependence amongst items within each domain or 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 domains is via the second-order model (e.g., de la Torre and Song Reference de la Torre and Song2009; Rijmen Reference Rijmen2010), where items are indirectly mapped to an overall (second-order) factor via non-overlapping group-specific (first-order) factors. Second-order models are suitable when the first-order factors are associated with each other, and there is a second-order factor that accounts for the relations among the first-order factors. The second-order model can be described as an independent clusters factor model (McDonald Reference McDonald1999) with a single second-order factor.

The bi-factor and the second-order models are not generally equivalent (Yung et al. Reference Yung, Thissen and McLeod1999; Gustafsson and Balke Reference Gustafsson and Balke1993; Mulaik and Quartetti Reference Mulaik and Quartetti1997; Rijmen Reference Rijmen2010), unless proportionality constraints are imposed by using the Schmid–Leiman transformation method (Schmid and Leiman Reference Schmid and Leiman1957). More importantly, both models are restricted to the MVN assumption for the latent variables, which might not be valid. Nikoloulopoulos and Joe (Reference Nikoloulopoulos and Joe2015) emphasized that if the ordinal variables in item response can be thought of as discretization of latent random variables that are maxima/minima or mixtures of means, then the use of factor models based on the MVN 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 a respondent considers specific items. An item might make the respondent think about M past events which, say, have values $W_{1}, \dots, W_{M}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$W_1,\ldots ,W_M$$\end{document} . In answering the item, the subject might take the average, maximum or minimum of $W_{1}, \dots, W_{M}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$W_1,\ldots ,W_M$$\end{document} and then convert to the ordinal scale depending on the magnitude. The case of a latent maxima/minima can occur if the response is based on a best or worst case.

Nikoloulopoulos and Joe (Reference Nikoloulopoulos and Joe2015) have studied factor copula models for item response where for the first factor there are bivariate copulas that couple each item to the first latent variable, and for the second factor there are copulas that link each item to the second latent variable conditioned on the first factor (leading to conditional dependence parameters), etc. They have shown that there is an improvement on the factor models based on the MVN assumption for the latent variables both conceptually and in fit to data. This improvement relies on the aforementioned reasons, i.e., items can have more probability in joint upper or lower tail than would be expected with a discretized MVN or items can be considered as discretized maxima/minima or mixtures of discretized means rather than discretized means. When all the bivariate copulas are bivariate normal (BVN), then the resulting model is the same as the discretized MVN model with a p-factor correlation matrix (Maydeu-Olivares Reference Maydeu-Olivares2006), also known as the p-dimensional normal ogive model (Jöreskog and Moustaki Reference Jöreskog and Moustaki2001). For example, the 1-factor copula model with BVN copulas is the same as the variant of Samejima’s (Reference Samejima1969) graded response IRT model, known as normal ogive model (McDonald Reference McDonald, van der Linden and Hambleton1997) with an 1-factor correlation matrix. We refer to Nikoloulopoulos and Joe (Reference Nikoloulopoulos and Joe2015, Section 2.3) for further details and explanations on the normal ogive models as special cases of factor copula models.

In this paper, we propose copula extensions for bi-factor and second-order models. The construction of the bi-factor copula model exploits the use of bivariate copulas that link the items to the common and group-specific factors. Note that if there is only one group of items, then the bi-factor model reduces to the two-factor copula model in Nikoloulopoulos and Joe (Reference Nikoloulopoulos and Joe2015). Similarly with the bi-factor copula model, we also use bivariate copulas to construct the second-order copula model. In this case, there are bivariate copulas that link the items 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 items and group-specific factors, each group of variables in fact is modelled using the one-factor copula model proposed by Nikoloulopoulos and Joe (Reference Nikoloulopoulos and Joe2015). In addition, if there is only one group of items, then the second-order copula model reduces to the one-factor copula model. Hence, the proposed models contain the one- and two-factor copula models in Nikoloulopoulos and Joe (Reference Nikoloulopoulos and Joe2015) as special cases, while allowing flexible dependence structure for 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 copula constructions are truncated vine copula models (Brechmann et al. Reference Brechmann, Czado and Aas2012) that involve both observed and latent variables. Joe et al. (Reference Joe, Li and Nikoloulopoulos2010) have shown that by choosing bivariate linking copulas appropriately, truncated vine copula models can have a wide range of asymmetric dependence as well as tail dependence (dependence among extreme values) and different lower/upper tail dependence parameters for each bivariate margin. Hence, the bi-factor and second-order copula models will be useful when the items have more probability in joint upper or lower tail than would be expected with a discretized MVN. If the bivariate linking copulas are BVN, then the Gaussian bi-factor and second-order models are special cases of our constructions which are the discrete counterparts of the structured factor copula models Krupskii and Joe (Reference Krupskii and Joe2015) where dependence and tail properties are obtained.

The remainder of the paper proceeds as follows: Section 1 introduces the bi-factor and second-order copula models for item response and discusses their relationship with the existing models. Estimation techniques and computational details are provided in Sect. 2. Section 3 proposes simple diagnostics based on semi-correlations and a heuristic method to select suitable bivariate copulas and build plausible bi-factor and second-order copula models. Section 4 summarizes the assessment of goodness of fit of these models using the $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistic of Maydeu-Olivares and Joe (Reference Maydeu-Olivares and Joe2006), which is based on a quadratic form of the deviations of sample and model-based proportions over all bivariate margins. Section 5 contains an extensive simulation study to gauge the small-sample efficiency of the proposed estimation, investigate the misspecification of the bivariate copulas, and examine the reliability of the model selection and goodness-of-fit techniques. Section 6 presents an application 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 discretized latent minima or mixtures of discretized means, provide better fit than the Gaussian bi-factor and second-order models. We conclude with some discussion in Sect. 7.

1. Bi-factor and Second-Order Copula Models

Let ${\underset{⏟}{Y_{11}, \dots, Y_{d_{1} 1}}}_{1}, \dots, {\underset{⏟}{Y_{1 g}, \dots, Y_{d_{g} g}}}_{g}, \dots, {\underset{⏟}{Y_{1 G}, \dots, Y_{d_{G} G}}}_{G}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\underbrace{Y_{11}, \ldots ,Y_{d_11}}_1,\ldots , \underbrace{Y_{1g},\ldots , Y_{d_gg}}_g, \ldots , \underbrace{Y_{1G}, \ldots , Y_{d_GG}}_G$$\end{document} denote the item response variables classified into the G non-overlapping groups. There are $d_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d_g$$\end{document} items in group g; $g = 1, \dots, G, j = 1, \dots, d_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g=1,\ldots ,G,\,j=1,\ldots ,d_g$$\end{document} and collectively there are $d = \sum_{g = 1}^{G} d_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d=\sum _{g=1}^{G} d_g$$\end{document} items, which are all measured on an ordinal scale; $Y_{jg} \in {0, \dots, K - 1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg}\in \{0,\ldots ,K-1\}$$\end{document} . Let the cutpoints in the uniform U(0, 1) scale for the jg’th item be $a_{j g, k}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$a_{jg,k}$$\end{document} , $k = 1, \dots, K - 1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$k=1,\ldots ,K-1$$\end{document} , with $a_{j g, 0} = 0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$a_{jg,0}=0$$\end{document} and $a_{j g, K} = 1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$a_{jg,K}=1$$\end{document} . These correspond to $a_{j g, k} = Φ (α_{j g, k})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$a_{jg,k}=\Phi (\alpha _{jg,k})$$\end{document} , where $α_{j g, k}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha _{jg,k}$$\end{document} are cutpoints in the normal N(0, 1) scale.

The bi-factor and second-order factor copula models are presented in Sects. 1.1 and 1.2, respectively. Section 1.3 discusses their relationship with the existing Gaussian bi-factor and second-order models, and Sect. 1.4 provides the bivariate linking copulas we consider along with their properties.

1.1. Bi-factor Copula Model

Consider a common factor $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} and G group-specific factors $V_{1}, \dots, V_{G}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_1,\ldots ,V_G$$\end{document} , where $V_{0}, V_{1}, \dots, V_{G}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0,V_1,\ldots ,V_G$$\end{document} are independent and standard uniformly distributed. Let $Y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg}$$\end{document} be the jth observed variable in group g, with $y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$y_{jg}$$\end{document} being the realization. The bi-factor model assumes that $Y_{1 g}, \dots, Y_{d_{g} g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{1g},\ldots , Y_{d_gg}$$\end{document} are conditionally independent given $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} and $V_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_g$$\end{document} , and that $Y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg}$$\end{document} in group g does not depend on $V_{g^{'}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_{g'}$$\end{document} for $g \neq g^{'}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g\ne g'$$\end{document} . Figure 1 depicts a graphical representation of the model.

Figure 1. Graphical representation of the bi-factor copula model with G group-specific factors and a common factor $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} .

The joint probability mass function (pmf) is given by:

\begin{matrix} π (y) = & Pr (Y_{jg} = y_{jg}; j = 1, \dots, d_{g}, g = 1, \dots, G) \\ = & \int_{{[0, 1]}^{G + 1}} \prod_{g = 1}^{G} \prod_{j = 1}^{d_{g}} Pr (Y_{jg} = y_{jg} | V_{0} = v_{0}, V_{g} = v_{g}) d v_{1} \dots d v_{G} d v_{0} . \end{matrix}

According to Sklar’s theorem (Reference Sklar1959), there exists a bivariate copula $C_{Y_{jg}, V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_0}$$\end{document} such that $Pr (Y_{jg} \leq y_{jg}, V_{0} \leq v_{0}) = C_{Y_{jg}, V_{0}} (F_{Y_{jg}} (y_{jg}), v_{0})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Pr (Y_{jg} \le y_{jg}, V_0 \le v_0)= C_{Y_{jg},V_0}\bigl (F_{{Y_{jg}}}(y_{jg}), v_0\bigr )$$\end{document} , for $v_{0} \in [0, 1]$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$v_0 \in [0,1]$$\end{document} , where $C_{Y_{jg}, V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_0}$$\end{document} is the copula that links the item $Y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg}$$\end{document} with the common factor $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} , $F_{Y_{jg}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$F_{Y_{jg}}$$\end{document} is the cumulative distribution function (cdf) of $Y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg}$$\end{document} ; note that $F_{Y_{jg}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$F_{Y_{jg}}$$\end{document} is a step function with jumps at $0, \dots, K - 1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$0,\ldots ,K-1$$\end{document} , i.e., $F_{Y_{jg}} (y_{jg}) = a_{j g, y_{jg} + 1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$F_{Y_{jg}}(y_{jg})=a_{jg,y_{jg}+1}$$\end{document} . Then, it follows that,

\begin{matrix} F_{Y_{jg} | V_{0}} (y_{jg} | v_{0}) : = Pr (Y_{jg} \leq y_{jg} | V_{0} = v_{0}) = \frac{\partial}{\partial v_{0}} C_{Y_{jg}, V_{0}} (a_{j g, y_{jg} + 1}, v_{0}) . \end{matrix}

For shorthand notation, we let $C_{Y_{jg} | V_{0}} (a_{j g, y_{jg} + 1} | v_{0}) = \frac{\partial}{\partial v_{0}} C_{Y_{jg}, V_{0}} (a_{j g, y_{jg} + 1}, v_{0})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg}|V_0}\bigl (a_{jg,y_{jg}+1}|v_0\bigr ) = \frac{\partial }{\partial v_0} C_{Y_{jg},V_0}\bigl ( a_{jg,y_{jg}+1}, v_0\bigr )$$\end{document} .

The observed variables also load on the group-specific factors; hence to account for this dependence, we let $C_{Y_{jg}, V_{g} | V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_g|V_0 }$$\end{document} be a bivariate copula that links the item $Y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg}$$\end{document} with the group-specific factor $V_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_g$$\end{document} given the common factor $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} . Hence,

\begin{matrix} Pr (Y_{jg} \leq y_{jg} | V_{0} = v_{0}, V_{g} = v_{g}) = \frac{\partial}{\partial v_{g}} Pr (Y_{jg} \leq y_{jg}, V_{g} \leq v_{g} | V_{0} = v_{0}) \\ = \frac{\partial}{\partial v_{g}} C_{Y_{jg}, V_{g} | V_{0}} (F_{Y_{j} g | V_{0}} (y_{jg} | v_{0}), v_{g}) = C_{Y_{jg} | V_{g}; V_{0}} (F_{Y_{j} g | V_{0}} (y_{jg} | v_{0}) | v_{g}) . \end{matrix}

To this end, the pmf of the bi-factor copula model takes the form

(1)

\begin{matrix} π (y) & = \int_{{[0, 1]}^{G + 1}} \prod_{g = 1}^{G} \prod_{j = 1}^{d_{g}} {C_{Y_{jg} | V_{g}; V_{0}} (F_{Y_{jg} | V_{0}} (y_{jg} | v_{0}) | v_{g}) \\ - C_{Y_{jg} | V_{g}; V_{0}} (F_{Y_{jg} | V_{0}} (y_{jg} - 1 | v_{0}) | v_{g})} d v_{1} \dots d v_{G} d v_{0} \\ = \int_{0}^{1} \prod_{g = 1}^{G} {\int_{0}^{1} \prod_{j = 1}^{d_{g}} [C_{Y_{jg} | V_{g}; V_{0}} (F_{Y_{jg} | V_{0}} (y_{jg} | v_{0}) | v_{g}) \\ - C_{Y_{jg} | V_{g}; V_{0}} (F_{Y_{jg} | V_{0}} (y_{jg} - 1 | v_{0}) | v_{g})] d v_{g}} d v_{0} \\ = \int_{0}^{1} \prod_{g = 1}^{G} {\int_{0}^{1} \prod_{j = 1}^{d_{g}} [C_{Y_{jg} | V_{g}; V_{0}} (C_{Y_{jg} | V_{0}} (a_{j g, y_{jg} + 1} | v_{0}) | v_{g}) \\ - C_{Y_{jg} | V_{g}; V_{0}} (C_{Y_{jg} | V_{0}} (a_{j g, y_{jg}} | v_{0}) | v_{g})] d v_{g}} d v_{0} \\ = \int_{0}^{1} \prod_{g = 1}^{G} {\int_{0}^{1} \prod_{j = 1}^{d_{g}} f_{Y_{jg} | V_{g}; V_{0}} (y_{jg} | v_{g}, v_{0}) d v_{g}} d v_{0} . \end{matrix}

\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \pi (\textbf{y})&= \int _{[0,1]^{G+1}} \prod _{g=1}^{G} \prod _{j=1}^{d_g} \biggr \{C_{Y_{jg}|V_g;V_0}\bigl (F_{Y_{jg}|V_0}(y_{jg}|v_0)| v_g\bigr ) \nonumber \\ {}&\quad -C_{Y_{jg}|V_g;V_0}\bigr (F_{Y_{jg}|V_0}(y_{jg}-1|v_0)| v_g\bigl ) \biggl \} d{v_1}\cdots d{v_G} \textrm{d}{v_0} \nonumber \\&= \int _0^1 \prod _{g=1}^{G} \Bigg \{ \int _0^1 \prod _{j=1}^{d_g} \biggl [ C_{Y_{jg}|V_g;V_0}\bigl (F_{Y_{jg}|V_0}(y_{jg}|v_0)| v_g\bigr )\nonumber \\ {}&\quad -C_{Y_{jg}|V_g;V_0}\bigr (F_{Y_{jg}|V_0}(y_{jg}-1|v_0)| v_g\bigl ) \biggr ] \textrm{d}{v_g} \Bigg \} \textrm{d}{v_0} \nonumber \\&= \int _0^1 \prod _{g=1}^{G} \Bigg \{ \int _0^1 \prod _{j=1}^{d_g} \biggr [ C_{Y_{jg}|V_g;V_0}\bigl ( C_{Y_{jg}|V_0}(a_{jg,y_{jg}+1}|v_0) | v_g\bigr ) \nonumber \\ {}&\quad -C_{Y_{jg}|V_g;V_0}\bigr ( C_{Y_{jg}|V_0}(a_{jg,y_{jg}}|v_0) | v_g\bigr ) \biggl ] \textrm{d}{v_g} \Bigg \} \textrm{d}{v_0}\nonumber \\&= \int _0^1 \prod _{g=1}^{G} \Bigg \{ \int _0^1 \prod _{j=1}^{d_g} f_{Y_{jg}|V_{g};V_0}(y_{jg}|v_g,v_0) \textrm{d}{v_g} \Bigg \} \textrm{d}{v_0}. \end{aligned}$$\end{document}

It is shown that the pmf is represented as an one-dimensional integral of a function which is in turn a product of G one-dimensional integrals. Thus, we avoid $(G + 1)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(G+1)$$\end{document} -dimensional numerical integration.

In addition to the computational advancements the proposed model offers, it can provide, with appropriately chosen linking copulas, more probability in joint upper or lower tail than would be expected with a discretized MVN. The bi-factor copula can be explained as a 2-truncated vine. d-dimensional vine copulas can cover flexible dependence structures through the specification of d bivariate marginal copulas at level 1 and $d (d - 1) / 2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d(d-1)/2$$\end{document} bivariate conditional copulas at higher levels (Nikoloulopoulos et al. Reference Nikoloulopoulos, Joe and Li2012). For the d-dimensional bi-factor copula, the pairs at level 1 are $Y_{jg}, V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg},V_0$$\end{document} for $g = 1, \dots, G, j = 1, \dots, d_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g=1,\ldots ,G,\,j=1,\ldots ,d_g$$\end{document} , the pairs at level 2 are $Y_{jg}, V_{g} | V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg},V_g|V_0$$\end{document} for $g = 1, \dots, G, j = 1, \dots, d_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g=1,\ldots ,G,\,j=1,\ldots ,d_g$$\end{document} , and for higher levels the (conditional) copula pairs are set to independence. That is the bi-factor copula has d bivariate copulas $C_{Y_{jg}, V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_0}$$\end{document} that link $Y_{jg}, g = 1, \dots, G, j = 1, \dots, d_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg},\,g=1,\ldots ,G,\,j=1,\ldots ,d_g$$\end{document} with $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} in the 1st level of the vine, d bivariate copulas $C_{Y_{jg}, V_{g} | V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_g|V_0}$$\end{document} that link $Y_{jg}, g = 1, \dots, G, j = 1, \dots, d_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg},\,g=1,\ldots ,G,\,j=1,\ldots ,d_g$$\end{document} with $V_{g}, g = 1, \dots, G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_g,\,g=1,\ldots ,G$$\end{document} given $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} in the 2nd level of the vine, and independence copulas in all the remaining levels of the vine (truncated after the 2nd level). From results in Joe et al. (Reference Joe, Li and Nikoloulopoulos2010) and Krupskii and Joe (Reference Krupskii and Joe2015), upper or lower tail dependent copulas in levels 1 and 2 will lead to items that have more probability in joint upper or lower tail than would be expected with a discretized MVN.

For the parametric version of the bi-factor copula model, we let $C_{Y_{jg}, V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_0}$$\end{document} and $C_{Y_{jg}, V_{g} | V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_g|V_0}$$\end{document} be parametric copulas with dependence parameters $θ_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{jg}$$\end{document} and $δ_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\delta _{jg}$$\end{document} , respectively.

1.2. Second-Order Copula Model

Assume that for a fixed $g = 1, \dots, G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g=1,\ldots ,G$$\end{document} , the items $Y_{1 g}, \dots, Y_{d_{g} g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{1g},\ldots , Y_{d_gg}$$\end{document} are conditionally independent given the first-order factors $V_{g} \sim U (0, 1), g = 1, \dots, G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_g \sim U(0,1),\,g=1,\ldots ,G$$\end{document} and that $V = (V_{1}, \dots, V_{G})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{V}=(V_1,\ldots ,V_G)$$\end{document} are conditionally independent given the second-order factor $V_{0} \sim U (0, 1)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0 \sim U(0,1)$$\end{document} . That is the joint distribution of $V$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{V}$$\end{document} has an one-factor structure. We also assume that $Y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg}$$\end{document} in group g does not depend on $V_{g^{'}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_{g'}$$\end{document} for $g \neq g^{'}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g\ne g'$$\end{document} . Figure 2 depicts the graphical representation of the model.

The joint pmf takes the form

\begin{matrix} π (y) = \int_{{[0, 1]}^{G}} {\prod_{g = 1}^{G} \prod_{j = 1}^{d_{g}} Pr (Y_{jg} = y_{jg} | V_{g} = v_{g})} c_{V} (v_{1}, \dots, v_{G}) d v_{1} \dots d v_{G}; \end{matrix}

$c_{V}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c_\textbf{V}$$\end{document} is the one-factor copula density (Krupskii and Joe Reference Krupskii and Joe2013) of $V$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{V}$$\end{document} , viz.

\begin{matrix} c_{V} (v_{1}, \dots, v_{G}) = \int_{0}^{1} \prod_{g = 1}^{G} c_{V_{g}, V_{0}} (v_{g}, v_{0}) d v_{0}, \end{matrix}

where $c_{V_{g}, V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c_{V_g,V_0}$$\end{document} is the bivariate copula density of the copula $C_{V_{g}, V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{V_g,V_0}$$\end{document} linking $V_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_g$$\end{document} and $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} .

Figure 2. Graphical representation of the second-order copula model with G first-order factors and one second-order factor $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} .

Letting $C_{Y_{jg}, V_{g}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_g}$$\end{document} be a bivariate copula that joins the item $Y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg}$$\end{document} and the group-specific factor $V_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_g$$\end{document} such that

\begin{matrix} F_{Y_{jg} | V_{g}} (y_{jg} | v_{g}) : = Pr (Y_{jg} \leq y_{jg} | V_{g} = v_{g}) = \frac{\partial}{\partial v_{g}} C_{Y_{jg}, V_{g}} (a_{j g, y_{jg} + 1}, v_{g}) = C_{Y_{jg} | V_{g}} (a_{j g, y_{jg} + 1} | v_{g}), \end{matrix}

the pmf of the second-order copula model becomes

(2)

\begin{matrix} π (y) & = \int_{0}^{1} \int_{{[0, 1]}^{G}} {\prod_{g = 1}^{G} \prod_{j = 1}^{d_{g}} (C_{Y_{jg} | V_{g}} (a_{j g, y_{jg} + 1} | v_{g}) - C_{Y_{jg} | V_{g}} (a_{j g, y_{jg}} | v_{g}))} \\ \times {\prod_{g = 1}^{G} c_{V_{g}, V_{0}} (v_{g}, v_{0})} d v_{1} \dots d v_{G} d v_{0} \\ = \int_{0}^{1} {\prod_{g = 1}^{G} \int_{0}^{1} [\prod_{j = 1}^{d_{g}} (C_{Y_{jg} | V_{g}} (a_{j g, y_{jg} + 1} | v_{g}) - C_{Y_{jg} | V_{g}} (a_{j g, y_{jg}} | v_{g}))] c_{V_{g}, V_{0}} (v_{g}, v_{0}) d v_{g}} d v_{0} \\ = \int_{0}^{1} {\prod_{g = 1}^{G} \int_{0}^{1} [\prod_{j = 1}^{d_{g}} f_{Y_{jg} | V_{g}} (y_{jg} | v_{g})] c_{V_{g}, V_{0}} (v_{g}, v_{0}) d v_{g}} d v_{0} . \end{matrix}

\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \pi (\textbf{y})&= \int _0^1 \int _{[0,1]^G} \Bigg \{ \prod _{g=1}^{G} \prod _{j=1}^{d_g} \Big ( C_{Y_{jg}|V_g}\bigl (a_{jg,y_{jg}+1} | v_g\bigr ) - C_{Y_{jg}|V_g}\bigl (a_{jg,y_{jg}} | v_g \bigr ) \Big ) \Bigg \}\nonumber \\ {}&\quad \times \Bigg \{ \prod _{g=1}^{G} c_{V_g,V_0}\bigl (v_g,v_0\bigr ) \Bigg \} \textrm{d}{v_1} \cdots \textrm{d}{v_G} \textrm{d}{v_0}\nonumber \\&= \int _0^1 \Bigg \{ \prod _{g=1}^{G} \int _0^1 \Bigg [ \prod _{j=1}^{d_g} \Big ( C_{Y_{jg}|V_g}\bigl (a_{jg,y_{jg}+1} | v_g\bigr ) - C_{Y_{jg}|V_g}\bigl (a_{jg,y_{jg}} | v_g \bigr ) \Big ) \Bigg ] c_{V_g,V_0}\bigl (v_g,v_0\bigr ) \textrm{d}{v_g} \Bigg \} \textrm{d}{v_0}\nonumber \\&=\int _0^1 \Bigg \{ \prod _{g=1}^{G} \int _0^1 \Big [ \prod _{j=1}^{d_g} f_{Y_{jg}|V_g}(y_{jg}|v_g) \Big ] c_{V_g,V_0}\bigl (v_g,v_0\bigr ) \textrm{d}{v_g} \Bigg \} \textrm{d}{v_0}. \end{aligned}$$\end{document}

Similarly with the bi-factor copula model, the pmf is represented as an one-dimensional integral of a function, which is in turn a product of G one-dimensional integrals.

In addition to the computational advancements, the second-order model offers, it can provide, with appropriately chosen linking copulas, more probability in joint upper or lower tail than would be expected with a discretized MVN. The second-order copula can be explained as an 1-truncated vine. For the d-dimensional second-order copula, the pairs at level 1 are $Y_{jg}, V_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg},V_g$$\end{document} for $g = 1, \dots, G, j = 1, \dots, d_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g=1,\ldots ,G,\,j=1,\ldots ,d_g$$\end{document} and $V_{0}, V_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0,V_g$$\end{document} for $g = 1, \dots, G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g=1,\ldots ,G$$\end{document} , and for higher levels the (conditional) copula pairs are set to independence. That is the second copula has d bivariate copulas $C_{Y_{jg}, V_{g}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_g}$$\end{document} that link $Y_{jg}, g = 1, \dots, G, j = 1, \dots, d_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg},\,g=1,\ldots ,G,\,j=1,\ldots ,d_g$$\end{document} with $V_{g}, g = 1, \dots, G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_g,\,g=1,\ldots ,G$$\end{document} and G bivariate copulas $C_{V_{g}, V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{V_g,V_0}$$\end{document} that link $V_{g}, g = 1, \dots, G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_{g},\,g=1,\ldots ,G$$\end{document} with $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} in the 1st level of the vine, and independence copulas in all the remaining levels of the vine (truncated after the 1st level). Joe et al. (Reference Joe, Li and Nikoloulopoulos2010) have shown that in order for a vine copula to have tail dependence for all bivariate margins, it is only necessary for the bivariate copulas in level 1 to have tail dependence and it is not necessary for the conditional bivariate copulas in levels $2, \dots, d$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$2,\ldots ,d$$\end{document} to have tail dependence. Hence, upper or lower tail dependent copulas in level 1 will lead to will lead to items that have more probability in joint upper or lower tail than would be expected with a discretized MVN.

For the parametric version of the second-order copula model, we let $C_{Y_{jg}, V_{g}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_g}$$\end{document} and $C_{V_{g}, V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{V_g,V_0}$$\end{document} be parametric copulas with dependence parameters $θ_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{jg}$$\end{document} and $δ_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\delta _{g}$$\end{document} , respectively.

1.3. Special Cases

In this subsection, we show what happens when all bivariate copulas are BVN. Let $Z_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_{jg}$$\end{document} be the underlying continuous variable of the ordinal variable $Y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg}$$\end{document} , i.e., $Y_{jg} = y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg} = y_{jg}$$\end{document} if $α_{j g, y_{jg}} \leq Z_{jg} \leq α_{j g, y_{jg} + 1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha _{jg,y_{jg}} \le Z_{jg} \le \alpha _{jg,y_{jg}+1}$$\end{document} with $α_{j g, K} = \infty$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha _{jg,K} = \infty $$\end{document} and $α_{j g, 0} = - \infty$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha _{jg,0} = -\infty $$\end{document} .

For the bi-factor model, if $C_{Y_{jg}, V_{0}} (\cdot; θ_{jg})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_0}(\cdot ;\theta _{jg})$$\end{document} and $C_{Y_{jg}, V_{g} | V_{0}} (\cdot; δ_{jg})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_g|V_0}(\cdot ;\delta _{jg})$$\end{document} are BVN copulas,

\begin{matrix} C_{Y_{jg} | V_{g}; V_{0}} (C_{Y_{jg} | V_{0}} (a_{j g, y_{jg} + 1} | v_{0}) | v_{g}) = Φ (\frac{α_{j g, y_{jg} + 1} - θ_{jg} Φ^{- 1} (v_{0}) - δ_{jg} \sqrt{1 - θ_{jg}^{2}} Φ^{- 1} (v_{g})}{\sqrt{(1 - θ_{jg}^{2}) (1 - δ_{jg}^{2})}}) . \end{matrix}

Hence, the pmf for the bi-factor copula model in (1) becomes:

\begin{matrix} π (y) = & \int_{0}^{1} \prod_{g = 1}^{G} {\int_{0}^{1} \prod_{j = 1}^{d_{g}} [Φ (\frac{α_{j g, y_{jg} + 1} - θ_{jg} Φ^{- 1} (v_{0}) - δ_{jg} \sqrt{1 - θ_{jg}^{2}} Φ^{- 1} (v_{g})}{\sqrt{(1 - θ_{jg}^{2}) (1 - δ_{jg}^{2})}}) \\ - Φ (\frac{α_{j g, y_{jg}} - θ_{jg} Φ^{- 1} (v_{0}) - δ_{jg} \sqrt{1 - θ_{jg}^{2}} Φ^{- 1} (v_{g})}{\sqrt{(1 - θ_{jg}^{2}) (1 - δ_{jg}^{2})}})] d v_{g}} d v_{0} \\ = & \int_{- \infty}^{\infty} \prod_{g = 1}^{G} {\int_{- \infty}^{\infty} \prod_{j = 1}^{d_{g}} [Φ (\frac{α_{j g, y_{jg} + 1} - θ_{jg} z_{0} - δ_{jg} \sqrt{1 - θ_{jg}^{2}} z_{g}}{\sqrt{(1 - θ_{jg}^{2}) (1 - δ_{jg}^{2})}}) \\ - Φ (\frac{α_{j g, y_{jg}} - θ_{jg} z_{0} - δ_{jg} \sqrt{1 - θ_{jg}^{2}} z_{g}}{\sqrt{(1 - θ_{jg}^{2}) (1 - δ_{jg}^{2})}})] ϕ (z_{g}) d z_{g}} ϕ (z_{0}) d z_{0} . \end{matrix}

\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \pi (\textbf{y})= & {} \int _0^1 \prod _{g=1}^G \Biggl \{ \int _0^1 \prod _{j=1}^{d_g} \Biggl [ \Phi \left( \frac{ \alpha _{jg,y_{jg}+1} - \theta _{jg} \Phi ^{-1}(v_0) - \delta _{jg} \sqrt{1-\theta _{jg}^2} \Phi ^{-1}(v_g)}{\sqrt{(1 - \theta _{jg}^2) (1-\delta _{jg}^2) }} \right) \\{} & {} -\Phi \left( \frac{ \alpha _{jg,y_{jg}} - \theta _{jg} \Phi ^{-1}(v_0) - \delta _{jg} \sqrt{1-\theta _{jg}^2} \Phi ^{-1}(v_g)}{\sqrt{(1 - \theta _{jg}^2) (1-\delta _{jg}^2) }} \right) \Biggr ] dv_g \Biggr \} \textrm{d}v_0 \\= & {} \int _{-\infty }^{\infty } \prod _{g=1}^G \Biggl \{ \int _{-\infty }^{\infty } \prod _{j=1}^{d_g} \Biggl [ \Phi \left( \frac{ \alpha _{jg,y_{jg}+1} - \theta _{jg} z_{0} - \delta _{jg} \sqrt{1-\theta _{jg}^2} z_{g}}{\sqrt{(1 - \theta _{jg}^2) (1-\delta _{jg}^2) }} \right) \\{} & {} - \Phi \left( \frac{ \alpha _{jg,y_{jg}} - \theta _{jg} z_{0} - \delta _{jg} \sqrt{1-\theta _{jg}^2} z_{g}}{\sqrt{(1 - \theta _{jg}^2) (1-\delta _{jg}^2) }} \right) \Biggr ] \phi (z_{g}) dz_{g} \Biggr \} \phi (z_{0}) \textrm{d}z_{0}. \end{aligned}$$\end{document}

This model is the same as the Gaussian bi-factor model (Gibbons and Hedeker Reference Gibbons and Hedeker1992; Gibbons et al. Reference Gibbons, Bock, Hedeker, Weiss, Segawa, Bhaumik, Kupfer, Frank, Grochocinski and Stover2007) with stochastic representation

(3)

\begin{matrix} Z_{jg} = θ_{jg} Z_{0} + γ_{jg} Z_{g} + \sqrt{1 - θ_{jg}^{2} - γ_{jg}^{2}} ϵ_{jg}, g = 1, \dots, G, j = 1, \dots, d_{g}, \end{matrix}

where $γ_{jg} = δ_{jg} \sqrt{1 - θ_{jg}^{2}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _{jg}=\delta _{jg}\sqrt{1 - \theta _{jg}^2}$$\end{document} and $Z_{0}, Z_{g}, ϵ_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_{0},Z_{g},\epsilon _{jg}$$\end{document} are iid N(0, 1) random variables. The parameter $θ_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{jg}$$\end{document} of $C_{Y_{jg}, V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_0}$$\end{document} is the correlation of $Z_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_{jg}$$\end{document} and $Z_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_0$$\end{document} , and the parameter $δ_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\delta _{jg}$$\end{document} of $C_{Y_{jg}, V_{g} | V_{0}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{Y_{jg},V_g|V_0}$$\end{document} is the partial correlation between $Z_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_{jg}$$\end{document} and $Z_{g} = Φ^{- 1} (V_{g})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_{g}=\Phi ^{-1}(V_g)$$\end{document} given $Z_{0} = Φ^{- 1} (V_{0})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_0=\Phi ^{-1}(V_0)$$\end{document} .

It implies that the underlying random variables $Z_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_{jg}$$\end{document} ’s have a multivariate Gaussian distribution where the off-diagonal entries of the correlation matrix have the form $θ_{j_{1} g} θ_{j_{2} g} + γ_{j_{1} g} γ_{j_{2} g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{j_1g}\theta _{j_2g} + \gamma _{j_1g}\gamma _{j_2g}$$\end{document} and $θ_{j_{1} g_{1}} θ_{j_{2} g_{2}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{j_1g_1}\theta _{j_2g_2}$$\end{document} for $j_{1} \neq j_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$ j_1 \ne j_2$$\end{document} and $g_{1} \neq g_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g_1 \ne g_2$$\end{document} , respectively. For the Gaussian bi-factor model to be identifiable, the number of dependence parameters has to be $2 d - N_{1} - N_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$2d - N_{1}-N_{2}$$\end{document} , where $N_{1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$N_{1}$$\end{document} and $N_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$N_{2}$$\end{document} 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}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_g$$\end{document} is absorbed with $ϵ_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\epsilon _{jg}$$\end{document} because $γ_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _{jg}$$\end{document} would not be identifiable. For a group g of size 2 with variable indices $j_{1}, j_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$j_1,j_2$$\end{document} , the parameters $γ_{j_{1} g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _{j_1g}$$\end{document} and $γ_{j_{2} g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _{j_2g}$$\end{document} appear only in one correlation; hence, one of $γ_{j_{1} g}, γ_{j_{2} g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _{j_1g},\gamma _{j_2g}$$\end{document} can be taken as 1 without loss of generality. For the bi-factor copula with non-Gaussian linking copulas, near non-identifiability can occur when there are groups of size 2; in this case, one of the linking copulas to the group latent variable can be fixed at comonotonicity.

For the Gaussian second-order model, let $Z_{0}, Z_{1}^{'}, \dots, Z_{G}^{'}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_0,Z_1',\ldots ,Z_G'$$\end{document} be the dependent latent N(0, 1) variables, where $Z_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_0$$\end{document} is the second-order factor and $Z_{g}^{'} = β_{g} Z_{0} + \sqrt{1 - β_{g}^{2}} Z_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_g'=\beta _{g}Z_0+\sqrt{1-\beta _{g}^2}Z_g$$\end{document} is the first-order factor for group g. That is, there is an one second-order factor $Z_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_0$$\end{document} , and the first-order factors $Z_{1}^{'}, \dots, Z_{G}^{'}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_1',\ldots ,Z_G'$$\end{document} are linear combinations of the second-order factor, plus a unique variable $Z_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_g$$\end{document} for each first-order factor. The stochastic representation is (Krupskii and Joe Reference Krupskii and Joe2015):

\begin{matrix} Z_{jg} = & β_{jg} Z_{g}^{'} + \sqrt{1 - β_{jg}^{2}} ϵ_{jg} \\ Z_{g}^{'} = & β_{g} Z_{0} + \sqrt{1 - β_{g}^{2}} Z_{g}, g = 1, \dots, G, j = 1, \dots, d_{g}, \end{matrix}

(4)

\begin{matrix} Z_{jg} = β_{jg} β_{g} Z_{0} + β_{jg} \sqrt{1 - β_{g}^{2}} Z_{g} + \sqrt{1 - β_{jg}^{2}} ϵ_{jg}, g = 1, \dots, G, j = 1, \dots, d_{g} . \end{matrix}

Hence, this is a special case of (3) where $θ_{jg} = β_{jg} β_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{jg} = \beta _{jg} \beta _{g}$$\end{document} and $γ_{jg} = β_{jg} \sqrt{1 - β_{g}^{2}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _{jg} = \beta _{jg} \sqrt{1-\beta _{g}^2}$$\end{document} .

1.4. Other Choices of Parametric Bivariate Copulas

In line with Nikoloulopoulos and Joe (Reference Nikoloulopoulos and Joe2015), we use bivariate parametric copulas that can be used when considering latent maxima, minima or mixtures of means. For different dependent items based on latent maxima or minima, multivariate extreme value and copula theory (e.g., Joe Reference Joe1997) can be used to select suitable copulas that link observed to latent variables. Copulas that arise from extreme value theory have more probability in one joint tail (upper or lower) than expected with a discretized MVN distribution or a MVN copula with discrete margins. If item responses are based on discretizations of latent variables that are means, then it is possible that there can be more probability in both the joint upper and joint lower tail, compared with discretized MVN models. This happens if the respondents consist of a ‘mixture’ population (e.g., different locations or genders). From the theory of elliptical distributions and copulas (e.g., McNeil et al. Reference McNeil, Frey and Embrechts2005), it is known that the multivariate Student-t distribution as a scale mixture of MVN has more dependence in the tails. Extreme value and elliptical copulas can model item response data that have reflection asymmetric and symmetric dependence, respectively.

A bivariate copula C is reflection symmetric if its density satisfies $c (u_{1}, u_{2}) = c (1 - u_{1}, 1 - u_{2})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c(u_1,u_2)=c(1-u_1,1-u_2)$$\end{document} for all $0 \leq u_{1}, u_{2} \leq 1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$0\le u_1,u_2\le 1$$\end{document} . Otherwise, it is reflection asymmetric often with more probability in the joint upper tail or joint lower tail. Upper tail dependence means that $c (1 - u, 1 - u) = O (u^{- 1})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c(1-u,1-u)=O(u^{-1})$$\end{document} as $u \to 0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$u\rightarrow 0$$\end{document} and lower tail dependence means that $c (u, u) = O (u^{- 1})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c(u,u)=O(u^{-1})$$\end{document} as $u \to 0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$u\rightarrow 0$$\end{document} . If $(U_{1}, U_{2}) \sim C$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(U_1,U_2)\sim C$$\end{document} for a bivariate copula C, then $(1 - U_{1}, 1 - U_{2}) \sim \hat{C}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(1-U_1,1-U_2)\sim \widehat{C}$$\end{document} , where $\hat{C} (u_{1}, u_{2}) = u_{1} + u_{2} - 1 + C (1 - u_{1}, 1 - u_{2})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\widehat{C}(u_1,u_2)=u_1+u_2-1+C(1-u_1,1-u_2)$$\end{document} is the survival or reflected copula of C; this “reflection” of each uniform U(0, 1) random variable about 1/2 changes the direction of tail asymmetry.

After briefly providing definitions of tail dependence and reflection symmetry/asymmetry, we provide below the bivariate copula choices we consider:

• The extreme value Gumbel copula with cdf
$\begin{matrix} C (u_{1}, u_{2}; θ) = exp [- {{(- log u_{1})}^{θ} + {(- log u_{2})}^{θ}}^{1 / θ}], θ \geq 1 . \end{matrix}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} C(u_1,u_2;\theta )=\exp \Bigl [-\Bigl \{(-\log u_1)^{\theta } +(-\log u_2)^{\theta }\Bigr \}^{1/\theta }\Bigr ], \theta \ge 1. \end{aligned}$$\end{document}
A model with bivariate Gumbel copulas has latent (ordinal) variables that can be considered as (discretized) maxima, and there is more probability in the joint upper tail as the Gumbel copula has reflection asymmetry and upper tail dependence.
• The survival Gumbel (s.Gumbel) copula with cdf
$\begin{matrix} C (u_{1}, u_{2}; θ) = & u_{1} + u_{2} - 1 \\ + exp [- {(- log (1 - u_{1}))^{θ} + (- log (1 - u_{2}))^{θ}}^{1 / θ}], θ \geq 1 . \end{matrix}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} C(u_1,u_2;\theta )= & {} u_1+u_2-1 \\ {}{} & {} + \exp \Bigl [-\Bigl \{\bigl (-\log (1-u_1)\bigr )^{\theta } +\bigl (-\log (1-u_2)\bigr )^{\theta }\Bigr \}^{1/\theta }\Bigr ], \theta \ge 1. \end{aligned}$$\end{document}
A model with bivariate s.Gumbel copulas has latent (ordinal) variables that can be considered as (discretized) minima, and there is more probability in the joint lower tail as the s.Gumbel copula has reflection asymmetry and lower tail dependence.
• The elliptical bivariate $t_{ν}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_\nu $$\end{document} copula with cdf
$\begin{matrix} C (u_{1}, u_{2}; θ) = T_{2} (T^{- 1} (u_{1}; ν), T^{- 1} (u_{2}; ν); θ, ν), - 1 \leq θ \leq 1, \end{matrix}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} C(u_1,u_2;\theta )=\mathcal {T}_2\Bigl (\mathcal {T}^{-1}(u_1;\nu ),\mathcal {T}^{-1}(u_2;\nu );\theta ,\nu \Bigr ),-1\le \theta \le 1, \end{aligned}$$\end{document}
where $T (; ν)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\mathcal {T}(;\nu )$$\end{document} is the univariate Student t cdf with (non-integer) $ν$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\nu $$\end{document} degrees of freedom, and $T_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\mathcal {T}_2$$\end{document} is the cdf of a bivariate Student-t distribution with $ν$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\nu $$\end{document} degrees of freedom and correlation parameter $θ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta $$\end{document} . A model with bivariate $t_{ν}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_\nu $$\end{document} copulas 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 $ν$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\nu $$\end{document} , such as $1 \leq ν \leq 5$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$1 \le \nu \le 5$$\end{document} , leads to a model with more probabilities in the joint upper and joint lower tails compared with the BVN copula as the $t_{ν}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_\nu $$\end{document} copula has reflection symmetric upper and lower tail dependence.

The BVN and $t_{ν}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_\nu $$\end{document} are comprehensive copulas, i.e., they interpolate between countermonotonicity (perfect negative dependence) to comonotonicity (perfect positive dependence), but the Gumbel copulas interpolates between independence and perfect positive dependence. Nevertheless, negative dependence or interpolation between perfect negative dependence and independence can be obtained from the Gumbel copulas by considering reflection of one of the uniform random variables on (0, 1). If $(U_{1}, U_{2}) \sim C$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(U_1,U_2)\sim C$$\end{document} for a bivariate copula C with positive dependence, then

• $(1 - U_{1}, U_{2}) \sim {\hat{C}}^{(1)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(1-U_1,U_2)\sim \widehat{C}^{(1)}$$\end{document} , where ${\hat{C}}^{(1)} (u_{1}, u_{2}) = u_{2} - C (1 - u_{1}, u_{2})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\widehat{C}^{(1)}(u_1,u_2)=u_2-C(1-u_1,u_2)$$\end{document} is the 1-reflected copula of C with negative lower-upper tail dependence;
• $(U_{1}, 1 - U_{2}) \sim {\hat{C}}^{(2)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(U_1,1-U_2)\sim \widehat{C}^{(2)}$$\end{document} , where ${\hat{C}}^{(2)} (u_{1}, u_{2}) = u_{1} - C (u_{1}, 1 - u_{2})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\widehat{C}^{(2)}(u_1,u_2)=u_1-C(u_1,1-u_2)$$\end{document} is the 2-reflected copula of C with negative upper-lower dependence.

Negative upper-lower tail dependence means that $c (1 - u, u) = O (u^{- 1})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c(1-u,u)=O(u^{-1})$$\end{document} as $u \to 0^{+}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$u\rightarrow 0^+$$\end{document} and negative lower-upper tail dependence means that $c (u, 1 - u) = O (u^{- 1})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c(u,1-u)=O(u^{-1})$$\end{document} as $u \to 0^{+}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$u\rightarrow 0^+$$\end{document} (Joe Reference Joe, Kurowicka and Joe2011). The proposed models can provide, with linking copulas that allow for negative (tail) dependence, negative (tail) dependence between the observed variables as the dependence between the observed and latent variables is inherited to the dependence amongst the observed variables.

In Fig. 3, to depict the concepts of refection symmetric or asymmetric tail dependence, we show contour plots of the corresponding copula densities with standard normal margins and dependence parameters corresponding to Kendall’s $τ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} value of 0.6 in absolute value. To make the models comparable, we convert the BVN/ $t_{ν}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_\nu $$\end{document} and (reflected) Gumbel copula parameters to Kendall’s $τ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} ’s via

(5)

\begin{matrix} τ (θ) = \frac{2}{π} arcsin (θ) \end{matrix}

and

(6)

\begin{matrix} τ (θ) = (\begin{matrix} 1 - θ^{- 1}, & C, \hat{C} \\ θ^{- 1} - 1, & {\hat{C}}^{(1)}, {\hat{C}}^{(2)} \end{matrix}), \end{matrix}

respectively. Sharper corners (relative to ellipse) in Fig. 3 indicate tail dependence.

Figure 3. Contour plots of bivariate copulas with standard normal margins and dependence parameters corresponding to Kendall’s $τ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} value of 0.6 in absolute value.

The Kendall’s tau parameters, which are strictly increasing functions of the copula parameters, account for the dependence dominated by the middle of the data, and they are expected to be similar among different families of bivariate copulas. However, the tail dependence varies and it is a property to consider when choosing among different families of bivariate copulas (Nikoloulopoulos and Karlis Reference Nikoloulopoulos and Karlis2008).

2. Estimation and Computational Details

For the set of all parameters, let $θ = (a, θ_{g}, δ_{g})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\theta }=(\textbf{a},\varvec{\theta }_g,\varvec{\delta }_{g})$$\end{document} for the bi-factor copula model and $θ = (a, θ_{g}, δ)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\theta }=(\textbf{a},\varvec{\theta }_g,\varvec{\delta })$$\end{document} for the second-order copula model, where $a = (a_{j g, k} : j = 1, \dots, d_{g}, g = 1, \dots, G, k = 1, \dots, K - 1)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{a}=(a_{jg,k}: j=1,\ldots ,d_g, g=1,\ldots ,G, k=1,\ldots ,K-1)$$\end{document} , $θ_{g} = (θ_{1 g}, \dots, θ_{jg}, \dots, θ_{d_{g} g} : g = 1, \dots, G)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\theta }_g=(\theta _{1g}, \ldots , \theta _{jg}, \ldots ,\theta _{d_gg}: g=1,\ldots ,G)$$\end{document} , $δ_{g} = (δ_{1 g}, \dots, δ_{jg}, \dots, δ_{d_{g} g} : g = 1, \dots, G)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\delta }_{g}=(\delta _{1g}, \ldots , \delta _{jg}, \ldots ,\delta _{d_gg}: g=1,\ldots ,G)$$\end{document} and $δ = (δ_{1}, \dots, δ_{G})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\delta }=(\delta _{1}, \ldots ,\delta _{G})$$\end{document} . The dimension of $a$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{a}$$\end{document} , $θ_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\theta }_g$$\end{document} , $δ_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\delta }_{g}$$\end{document} and $δ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\delta }$$\end{document} are $d (K - 1)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d(K-1)$$\end{document} , d, d and G, respectively. Hence, the dimension q of $θ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\theta }$$\end{document} is $d (K + 1)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d(K+1)$$\end{document} and $d K + G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$dK +G$$\end{document} for the bi-factor and second-order copula model, respectively.

With sample size n and data $y_{1}, \dots, y_{n}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{y}_1,\ldots ,\textbf{y}_n$$\end{document} , the joint log-likelihood of the bi-factor and second-order copula is

(7)

\begin{matrix} ℓ (θ; y_{1}, \dots, y_{n}) = \sum_{i = 1}^{n} log π (y_{i}; θ) . \end{matrix}

with $π (y_{i}; θ)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\pi (\textbf{y}_i;\varvec{\theta })$$\end{document} as in (1) and (2), respectively. Maximum likelihood (ML) estimation, i.e., maximization of (7), is numerically possible but time-consuming for large d because the large number of univariate cutpoints and dependence parameters. Hence, we approach estimation using the two-step Inference Function of Margins (IFM) method (Joe and Xu Reference Joe and Xu1996; Joe Reference Joe1997). Joe (Reference Joe2005) has established its asymptotic efficiency and has shown that can efficiently, in the sense of computing time and asymptotic variance, estimate the univariate and dependence parameters.

In the first step of the IFM, the univariate parameters, i.e., the cutpoints, are estimated using the univariate sample proportions. The univariate cutpoints for the jth item in group g are estimated as ${\hat{a}}_{j g, k} = \sum_{y = 0}^{k} p_{j g, y}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\hat{a}_{jg,k} = \sum _{y=0}^{k} p_{jg,y}$$\end{document} , where $p_{j g, y}, y = 0, \dots, K - 1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$p_{jg,y}\,,y=0,\ldots ,K-1$$\end{document} for $g = 1, \dots, G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g=1,\ldots ,G$$\end{document} and $j = 1, \dots, d_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$j=1,\ldots ,d_g$$\end{document} are the univariate sample proportions. In the second step of the IFM method, the joint log-likelihood in (7) is maximized over the copula parameters with the cutpoints fixed as estimated at the first step. That is for $i = 1, \dots, n$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i=1,\ldots ,n$$\end{document} we start from a d-variate sample $y_{i 11}, \dots, y_{i d_{1} 1}, \dots, y_{i 1 G}, \dots, y_{i d_{G} G}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$y_{i11}, \ldots ,y_{id_11},\ldots , y_{i1G}, \ldots , y_{id_GG}$$\end{document} from which d estimators $F_{11} (y_{i 11}), \dots, F_{d_{1} 1} (y_{i d_{1} 1}), \dots, F_{1 G} (y_{i 1 G}), \dots,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$F_{11}(y_{i11}), \ldots ,F_{d_11}(y_{id_11}), \ldots , F_{1G}(y_{i1G}), \ldots ,$$\end{document} $F_{d_{G} G} (y_{i d_{G} G})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$F_{d_GG}(y_{id_GG})$$\end{document} are obtained. We use these estimators, i.e., the cutpoints, to transform the $y_{i 11}, \dots, y_{i d_{1} 1}, \dots, y_{i 1 G},$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$y_{i11}, \ldots ,y_{id_11},\ldots , y_{i1G},$$\end{document} $\dots, y_{i d_{G} G}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\ldots , y_{id_GG}$$\end{document} sample to a uniform sample ${\hat{α}}_{11, y_{i 11} + 1}, \dots, {\hat{α}}_{d_{1} 1, y_{i d_{1} 1} + 1}, \dots {\hat{α}}_{1 G, y_{i 1 G} + 1}, \dots, {\hat{α}}_{d_{G} G, y_{i d_{G} G} + 1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{\alpha }}_{11,y_{i11}+1},\ldots ,{\hat{\alpha }}_{d_11,y_{id_11}+1},\ldots {\hat{\alpha }}_{1G,y_{i1G}+1}, \ldots , {\hat{\alpha }}_{d_GG,y_{id_GG}+1}$$\end{document} on ${[0, 1]}^{d}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$[0, 1]^d$$\end{document} and then fit the factor copula model at the second step. Hence, the IFM approach can be regarded as a two-step approach on the original data or simply as the standard one-step ML method on the transformed (copula) data. Note also in passing that compared to the ML, the IFM method is not as punishing for misspecification of the dependence structure (Joe and Xu Reference Joe and Xu1996; Xu Reference Xu1996).

For the bi-factor copula model, numerical evaluation of the joint pmf can be achieved with the following steps:

(a) Calculate Gauss–Legendre quadrature (Stroud and Secrest Reference Stroud and Secrest1966) points ${v_{q} : q = 1, \dots, n_{q}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{\mathcal {v}_q: q=1,\ldots ,n_q\}$$\end{document} and weights ${w_{q} : q = 1, \dots, n_{q}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{w_q: q=1,\ldots ,n_q\}$$\end{document} in terms of standard uniform.
(b) Numerically evaluate the joint pmf
$\begin{matrix} \int_{0}^{1} \prod_{g = 1}^{G} {\int_{0}^{1} \prod_{j = 1}^{d_{g}} f_{Y_{jg} | V_{jg}; V_{0}} (y_{jg} | v_{g}, v_{0}) d v_{g}} d v_{0} \end{matrix}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \int _0^1 \prod _{g=1}^{G} \Bigg \{ \int _0^1 \prod _{j=1}^{d_g} f_{Y_{jg}|V_{jg};V_0}(y_{jg}|v_g,v_0) \textrm{d}{v_g} \Bigg \} \textrm{d}{v_0} \end{aligned}$$\end{document}
in a double sum
$\begin{matrix} \sum_{q_{1} = 1}^{n_{q}} w_{q_{1}} \prod_{g = 1}^{G} {\sum_{q_{2} = 1}^{n_{q}} w_{q_{2}} \prod_{j = 1}^{d_{g}} f_{Y_{jg} | V_{jg}; V_{0}} (y_{jg} | v_{q_{2}}, v_{q_{1}})} . \end{matrix}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \sum _{q_1=1}^{n_q} w_{q_1} \prod _{g=1}^{G} \Bigg \{ \sum _{q_2=1}^{n_q} w_{q_2} \prod _{j=1}^{d_g} f_{Y_{jg}|V_{jg};V_0}(y_{jg}|\mathcal {v}_{q_2},\mathcal {v}_{q_1}) \Bigg \}. \end{aligned}$$\end{document}

For the second-order copula model, numerical evaluation of the joint pmf can be achieved with the following steps:

1. Calculate Gauss–Legendre quadrature points ${v_{q} : q = 1, \dots, n_{q}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{\mathcal {v}_q : q=1,\ldots , n_q \}$$\end{document} and weights ${w_{q} : q = 1, \dots, n_{q}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{w_q : q=1,\ldots , n_q \}$$\end{document} in terms of standard uniform.
2. Numerically evaluate the joint pmf
$\begin{matrix} \int_{0}^{1} {\prod_{g = 1}^{G} \int_{0}^{1} [\prod_{j = 1}^{d_{g}} f_{Y_{jg} | V_{g}} (y_{jg} | v_{g}; θ_{jg})] c_{V_{g}, V_{0}} (v_{g}, v_{0}; δ_{g}) d v_{g}} d v_{0} \end{matrix}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \int _0^1 \Bigg \{ \prod _{g=1}^{G} \int _0^1 \Big [ \prod _{j=1}^{d_g} f_{Y_{jg}|V_g}(y_{jg}|v_g;\theta _{jg}) \Big ] c_{V_g,V_0}\bigl (v_g,v_0;\delta _{g}\bigr ) \textrm{d}{v_g} \Bigg \} \textrm{d}{v_0} \end{aligned}$$\end{document}
in a double sum
$\begin{matrix} \sum_{q_{1} = 1}^{n_{q}} w_{q_{1}} {\prod_{g = 1}^{G} \sum_{q_{2} = 1}^{n_{q}} w_{q_{2}} [\prod_{j = 1}^{d_{g}} f_{Y_{jg} | V_{g}} (y_{jg} | v_{q_{2} | q_{1}}; θ_{jg})]}, \end{matrix}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \sum _{q_1=1}^{n_q} w_{q_1} \Bigg \{ \prod _{g=1}^{G} \sum _{q_2=1}^{n_q} w_{q_2} \Big [ \prod _{j=1}^{d_g} f_{Y_{jg}|V_g}(y_{jg}|\mathcal {v}_{q_2|q_1};\theta _{jg}) \Big ] \Bigg \}, \end{aligned}$$\end{document}
where $v_{q_{2} | q_{1}} = C_{Y_{jg} | V_{g}; V_{0}}^{- 1} (v_{q_{2}} | v_{q_{1}}; δ_{g})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\mathcal {v}_{q_2|q_1} = C^{-1}_{Y_{jg}|V_g;V_0}( {v}_{q_2} | {v}_{q_1};\delta _{g})$$\end{document} . Note that the independent quadrature points ${v_{q_{1}} : q_{1} = 1, \dots, n_{q}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{\mathcal {v}_{q_1}: q_1 = 1,\ldots ,n_q\}$$\end{document} and ${v_{q_{2}} : q_{2} = 1, \dots, n_{q}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{\mathcal {v}_{q_2}: q_2 = 1,\ldots ,n_q\}$$\end{document} have been converted to dependent quadrature points that have an one-factor copula distribution $C_{X} (\cdot; δ)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{X}(\cdot ;\varvec{\delta })$$\end{document} .

The estimated copula parameters can be obtained by using a quasi-Newton (Nash Reference Nash1990) method applied to the logarithm of the joint likelihood. With Gauss–Legendre quadrature, the same nodes and weights are used for different functions; this helps in yielding smooth numerical derivatives for numerical optimization via quasi-Newton. Our comparisons show that $n_{q} = 25$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n_q=25$$\end{document} quadrature points are adequate with good precision.

3. Bivariate Copula Selection

In the following subsections, we describe simple diagnostics based on semi-correlations and a heuristic method that automatically selects the bivariate parametric copula families that build either the bi-factor or the second-order copula model.

Choices of copulas with upper or lower tail dependence are better if the items have more probability in joint lower or upper tail than would be expected with the BVN copula. This can be shown with summaries of polychoric correlations in the upper and lower joint tail (Kadhem and Nikoloulopoulos Reference Kadhem and Nikoloulopoulos2021). In the context of items that can be split into G non-overlapping groups, such that there is homogeneous dependence within each group, it is sufficient to (a) summarize the average of the polychoric semi-correlations for all pairs within each of the G groups and for all pairs of items, and (b) not mix bivariate copulas for a single factor; hence, for both the bi-factor and second-order copula models we allow $G + 1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$G+1$$\end{document} different copula families, one for each group specific factor $V_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_g$$\end{document} and one for $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} .

We distinguish the simple descriptives, i.e., the semi-correlations, from the heuristic algorithm and model fitting. On the one hand, the descriptive statistics can suggest more probability to the lower or upper tail for many pairs of items, but bi-factor and second-order copula models with asymmetrical dependence can be used to check whether the two tails are significantly different.

3.1. Simple Diagnostics Based on Semi-correlations

Consider again the underlying N(0, 1) latent variables $Z_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_{jg}$$\end{document} ’s of the ordinal variables $Y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg}$$\end{document} ’s. The correlations of $Z_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Z_{jg}$$\end{document} ’s in the upper and lower tail, hereafter semi-correlations, are defined as (Joe Reference Joe2014, p. 71):

(8)

\begin{matrix} ρ_{N}^{+} = & Cor (Z_{jg}, Z_{j^{'} g^{'}} | Z_{jg} > 0, Z_{j^{'} g^{'}} > 0) \\ = & \frac{\int_{0}^{\infty} \int_{0}^{\infty} z_{1} z_{2} ϕ (z_{1}) ϕ (z_{2}) c (Φ (z_{1}), Φ (z_{2})) d z_{1} d z_{2} - (\int_{0}^{\infty} z ϕ (z) (1 - C_{2 | 1} (0.5 | Φ (z))) d z)^{2} / C (0.5, 0.5)}{\int_{0}^{\infty} z^{2} ϕ (z) (1 - C_{2 | 1} (0.5 | Φ (z))) d z - (\int_{0}^{\infty} z ϕ (z) (1 - C_{2 | 1} (0.5 | Φ (z))) d z)^{2} / C (0.5, 0.5)}; \\ ρ_{N}^{-} = & Cor (Z_{j_{1} g}, Z_{j_{2} g} | Z_{j_{1} g} < 0, Z_{j_{2} g} < 0) \\ = & \frac{\int_{- \infty}^{0} \int_{- \infty}^{0} z_{1} z_{2} ϕ (z_{1}) ϕ (z_{2}) c (Φ (z_{1}), Φ (z_{2})) d z_{1} d z_{2} - (\int_{- \infty}^{0} z ϕ (z) C_{2 | 1} (0.5 | Φ (z)) d z)^{2} / C (0.5, 0.5)}{\int_{- \infty}^{0} z^{2} ϕ (z) C_{2 | 1} (0.5 | Φ (z)) d z - (\int_{- \infty}^{0} z ϕ (z) C_{2 | 1} (0.5 | Φ (z)) d z)^{2} / C (0.5, 0.5)} . \end{matrix}

\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \rho _N^+= & {} \text{ Cor }\Bigl (Z_{jg},Z_{j'g'}|Z_{jg}>0,Z_{j'g'}>0\Bigr )\\= & {} \frac{\int _{0}^{\infty }\int _{0}^{\infty } z_1z_2\phi (z_1)\phi (z_2)c\bigl (\Phi (z_1),\Phi (z_2)\bigr )dz_1dz_2-{\biggr (\int _0^\infty z\phi (z)\Bigr (1- C_{2|1}\bigr (0.5|\Phi (z)\bigl )\Bigr )dz\biggr )^2}/{C(0.5,0.5)}}{\int _0^\infty z^2\phi (z)\Bigr (1- C_{2|1}\bigr (0.5|\Phi (z)\bigl )\Bigr )\textrm{d}z-{\biggr (\int _0^\infty z\phi (z)\Bigr (1- C_{2|1}\bigr (0.5|\Phi (z)\bigl )\Bigr )\textrm{d}z\biggr )^2}/{C(0.5,0.5)}};\nonumber \\ \rho _N^-= & {} \text{ Cor }\Bigl (Z_{j_1g},Z_{j_2g}|Z_{j_1g}<0,Z_{j_2g}<0\Bigr )\nonumber \\= & {} \frac{\int _{-\infty }^{0}\int _{-\infty }^{0} z_1z_2\phi (z_1)\phi (z_2)c\bigl (\Phi (z_1),\Phi (z_2)\bigr )\textrm{d}z_1\textrm{d}z_2-{\biggr (\int _{-\infty }^0 z\phi (z)C_{2|1}\bigr (0.5|\Phi (z)\bigl )\textrm{d}z\biggr )^2}/{C(0.5,0.5)}}{\int _{-\infty }^0 z^2\phi (z) C_{2|1}\bigr (0.5|\Phi (z)\bigl )dz-{\biggr (\int _{-\infty }^0 z\phi (z)C_{2|1}\bigr (0.5|\Phi (z)\bigl )\textrm{d}z\biggr )^2}/{C(0.5,0.5)}}.\nonumber \end{aligned}$$\end{document}

From the above expressions, it is clear that the semi-correlations depend only on the copula C of $(Φ (Z_{jg}),$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Bigl (\Phi (Z_{jg}),$$\end{document} $Φ (Z_{j^{'} g^{'}}))$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Phi (Z_{j'g'})\Bigr )$$\end{document} ; $C_{2 | 1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C_{2|1}$$\end{document} is the conditional copula cdf. For the BVN and $t_{ν}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_\nu $$\end{document} copulas $ρ_{N}^{-} = ρ_{N}^{+}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\rho _N^{-}=\rho _N^{+}$$\end{document} , while for the Gumbel and s.Gumbel copulas $ρ_{N}^{-} < ρ_{N}^{+}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\rho _N^{-}<\rho _N^{+}$$\end{document} and $ρ_{N}^{-} > ρ_{N}^{+}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\rho _N^{-}>\rho _N^{+}$$\end{document} , respectively. The sample versions of $ρ_{N}^{+}, ρ_{N}^{-}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\rho _N^{+},\rho _N^{-}$$\end{document} for item response data are the polychoric correlations in the joint lower and upper quadrants of $Y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg}$$\end{document} and $Y_{j^{'} g^{'}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{j'g'}$$\end{document} (Kadhem and Nikoloulopoulos Reference Kadhem and Nikoloulopoulos2021).

3.2. Selection Algorithm

We propose a heuristic method that selects appropriate bivariate copulas for each factor of the bi-factor and second-order copula models. It starts with an initial assumption, that all bivariate linking copulas are BVN copulas, i.e. the starting model is either the Gaussian bi-factor or second-order model, and then, sequentially other copulas with lower or upper tail dependence are assigned to the factors where necessary to account for more 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.
2. Fit all the possible bi-factor or second-order copula models, iterating over all the copula candidates that link all items $Y_{jg}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{jg}$$\end{document} ’s in group g or each group-specific factor $V_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_g$$\end{document} , respectively, to $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} .
3. Select the copula family that corresponds to the lowest Akaike information criterion (AIC), that is, $AIC = - 2 \times ℓ + 2 \times # copula parameters$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\text {AIC}= -2 \times \ell +2 \times \#\text {copula parameters}$$\end{document} .
4. Fix the selected copula family that links the observed (bi-factor model) or latent (second-order model) variables to $V_{0}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document} .
5. For $g = 1, \dots, G$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g=1,\dots ,G$$\end{document} :
- (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 $V_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_g$$\end{document} .
- (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 $V_{g}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_g$$\end{document} .

For vine copula models (bi-factor and second-order copula models are vine copula models that involve both observed and latent variables), Dissmann et al. (Reference Dissmann, Brechmann, Czado and Kurowicka2013) also found that bivariate copula selection based on AIC seems to be better than even using bivariate goodness-of-fit tests. The goodness-of-fit procedures involve a global distance measure between the model-based and empirical distribution; hence, they might not be sensitive to tail behaviours and are not diagnostics in the sense of suggesting improved parametric models in the case of small p-values (Joe Reference Joe2014, p. 254). A smaller AIC indicates a model that better approximates both the dependence structure of the data and the strength of dependence in the tails.

4. Goodness of Fit

We will use the limited information $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistic proposed by Maydeu-Olivares and Joe (Reference Maydeu-Olivares and Joe2006) to evaluate the overall fit of the proposed bi-factor and second-order copula models. The $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistic is based on a quadratic form of the deviations of sample and model-based proportions over all bivariate margins. For the bi-factor and second-order copula models with parameter vector $θ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\theta }$$\end{document} of dimension q, let $π_{2} (θ) = ({\dot{π}}_{1} {(θ)}^{⊤}, {\dot{π}}_{2} {(θ)}^{⊤})^{⊤}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\pi }_2(\varvec{\theta })=\bigl (\dot{\varvec{\pi }}_1(\varvec{\theta })^\top ,\dot{\varvec{\pi }}_2(\varvec{\theta })^\top \bigr )^\top $$\end{document} be the column vector of the univariate and bivariate model-based marginal probabilities that do not include category 0 with sample counterpart $p_{2} = {({\dot{p}}_{1}^{⊤}, {\dot{p}}_{2}^{⊤})}^{⊤}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{p}_2=(\dot{\textbf{p}}_1^\top ,\dot{\textbf{p}}_2^\top )^\top $$\end{document} . The total number of the univariate and bivariate residuals $(p_{2} - π_{2} (\hat{θ}))^{⊤}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\bigl (\textbf{p}_2-\varvec{\pi }_2({\hat{\varvec{\theta }}})\bigr )^\top $$\end{document} is

\begin{matrix} s = d (K - 1) + (\begin{matrix} d \\ 2 \end{matrix}) {(K - 1)}^{2}, \end{matrix}

where $d (K - 1)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d (K-1)$$\end{document} is the dimension of the univariate residuals and $(\begin{matrix} d \\ 2 \end{matrix}) {(K - 1)}^{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\left( {\begin{array}{c}d\\ 2\end{array}}\right) (K - 1)^2$$\end{document} is the dimension of the bivariate residuals excluding category 0.

With a sample size n, the limited-information $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistic is given by

\begin{matrix} M_{2} = M_{2} (\hat{θ}) = n (p_{2} - π_{2} (\hat{θ}))^{⊤} C_{2} (\hat{θ}) (p_{2} - π_{2} (\hat{θ})), \end{matrix}

with

\begin{matrix} C_{2} (θ) = Ξ_{2}^{- 1} - Ξ_{2}^{- 1} Δ_{2} {(Δ_{2}^{⊤} Ξ_{2}^{- 1} Δ_{2})}^{- 1} Δ_{2}^{⊤} Ξ_{2}^{- 1} = Δ_{2}^{(c)} ({[Δ_{2}^{(c)}]}^{⊤} Ξ_{2} Δ_{2}^{(c)})^{- 1} {[Δ_{2}^{(c)}]}^{⊤}, \end{matrix}

where $Δ_{2} = \partial π_{2} (θ) / \partial θ^{⊤}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Delta }_2=\partial \varvec{\pi }_2(\varvec{\theta })/\partial \varvec{\theta }^\top $$\end{document} is an $s \times q$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$s\times q$$\end{document} matrix of full column rank q with the first-order derivatives of the univariate and bivariate marginal probabilities with respect to the estimated model parameters (in the Supplementary Tables 1–5 of the electronic supplementary material, we provide details on the calculation of these derivatives), $Δ_{2}^{(c)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Delta }_2^{(c)}$$\end{document} is an $s \times (s - q)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$s\times (s-q)$$\end{document} orthogonal complement to $Δ_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Delta }_2$$\end{document} of full column rank $s - q$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$s-q$$\end{document} , such that ${[Δ_{2}^{(c)}]}^{⊤} Δ_{2} = 0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$[\varvec{\Delta }_2^{(c)}]^\top \varvec{\Delta }_2=\textbf{0}$$\end{document} and $Ξ_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Xi }_2$$\end{document} is the asymptotic $s \times s$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$s \times s$$\end{document} covariance matrix of $\sqrt{n} (p_{2} - π_{2} (\hat{θ}))^{⊤}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sqrt{n}\bigl (\textbf{p}_2-\varvec{\pi }_2({\hat{\varvec{\theta }}})\bigr )^\top $$\end{document} .

The asymptotic covariance matrix $Ξ_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Xi }_2$$\end{document} can be partitioned according to the partitioning of $p_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\textbf{p}_2$$\end{document} into $Ξ_{11} = \sqrt{n} Acov ({\dot{p}}_{1})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Xi }_{11}=\sqrt{n}\text{ Acov }(\dot{\textbf{p}}_1)$$\end{document} , $Ξ_{21} = \sqrt{n} Acov ({\dot{p}}_{2}, {\dot{p}}_{1})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Xi }_{21}=\sqrt{n}\text{ Acov }(\dot{\textbf{p}}_2,\dot{\textbf{p}}_1)$$\end{document} and $Ξ_{22} = \sqrt{n} Acov ({\dot{p}}_{2})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Xi }_{22}=\sqrt{n}\text{ Acov }(\dot{\textbf{p}}_2)$$\end{document} , where $Acov (\cdot)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\text{ Acov }(\cdot )$$\end{document} denotes asymptotic covariance matrix. The elements of $Ξ_{11}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Xi }_{11}$$\end{document} , $Ξ_{21}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Xi }_{21}$$\end{document} and $Ξ_{22}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Xi }_{22}$$\end{document} involve up to the 4-dimensional probabilities as shown below:

\begin{matrix} \sqrt{n} Acov (p_{j_{1}, y_{1}}, p_{j_{2}, y_{2}}) = π_{j_{1} j_{2}, y_{1} y_{2}} - π_{j_{1}, y_{1}} π_{j_{2}, y_{2}} \\ \sqrt{n} Acov (p_{j_{1} j_{2}, y_{1} y_{2}}, p_{j_{3}, y_{3}}) = π_{j_{1} j_{2} j_{3}, y_{1} y_{2} y_{3}} - π_{j_{1} j_{2}, y_{1} y_{2}} π_{j_{3}, y_{3}} \\ \sqrt{n} Acov (p_{j_{1} j_{2}, y_{1} y_{2}}, p_{j_{3} j_{4}, y_{3} y_{4}}) = π_{j_{1} j_{2} j_{3} j_{4}, y_{1} y_{2} y_{3} y_{4}} - π_{j_{1} j_{2}, y_{1} y_{2}} π_{j_{3} j_{4}, y_{3} y_{4}}, \end{matrix}

where $π_{j, y} = Pr (Y_{j} = y)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$ \pi _{j,y} = \Pr (Y_j = y)$$\end{document} , $π_{j_{1} j_{2}, y_{1} y_{2}} = Pr (Y_{j_{1}} = y_{1}, Y_{j_{2}} = y_{2})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\pi _{j_1j_2,y_1y_2}= \Pr (Y_{j_1} = y_1,Y_{j_2} = y_2)$$\end{document} , $π_{j_{1} j_{2} j_{3}, y_{1} y_{2} y_{3}} = Pr (Y_{j_{1}} = y_{1}, Y_{j_{2}} = y_{2}, Y_{j_{3}} = y_{3})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\pi _{j_1j_2j_3,y_1y_2y_3}=\Pr (Y_{j_1} = y_1,Y_{j_2} = y_2,Y_{j_3} = y_3)$$\end{document} , and $π_{j_{1} j_{2} j_{3} j_{4}, y_{1} y_{2} y_{3} y_{4}} = Pr (Y_{j_{1}} = y_{1}, Y_{j_{2}} = y_{2}, Y_{j_{3}} = y_{3}, Y_{j_{4}} = y_{4})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\pi _{j_1j_2j_3j_4,y_1y_2y_3y_4}=\Pr (Y_{j_1} = y_1,Y_{j_2} = y_2,Y_{j_3} = y_3,Y_{j_4} = y_4)$$\end{document} .

The limited information statistic $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} under the null hypothesis has an asymptotic distribution that is $χ^{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\chi ^2$$\end{document} with $s - q$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$s-q$$\end{document} degrees of freedom as the IFM estimate $\hat{θ}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{\varvec{\theta }}}$$\end{document} is $\sqrt{n}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sqrt{n}$$\end{document} -consistent.

5. Simulations

An extensive simulation study is conducted to (a) gauge the small-sample efficiency of the IFM estimation method and investigate the misspecification of the bivariate pair-copulas, (b) examine the reliability of using the heuristic algorithm to select the true (simulated) bivariate linking copulas, and (c) study the small-sample performance of the $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistic.

We randomly generate 1000 datasets with samples of size $n = 500$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n=500$$\end{document} or 1000 and $d = 16$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d=16$$\end{document} items, with $K = 3$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$K=3$$\end{document} or $K = 5$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$K=5$$\end{document} equally weighted categories, that are equally separated into $G = 4$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$G=4$$\end{document} non-overlapping groups from the bi-factor and second-order copula model. In each simulated model, we use different linking copulas to cover different types of dependence. To make the models comparable, we convert the BVN/ $t_{ν}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_\nu $$\end{document} and Gumbel/s.Gumbel copula parameters to Kendall’s $τ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} ’s via the relations in (5) and (6), respectively. For the bi-factor copula models, we set $τ (θ_{g}) = (0.45, 0.55, 0.65, 0.75)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau (\varvec{\theta }_g)=(0.45,0.55,0.65,0.75)$$\end{document} and $τ (δ_{g}) = (0.30, 0.35,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau (\varvec{\delta }_g)=(0.30,0.35,$$\end{document} 0.40, 0.50) for $g = 1, \dots, 4$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g=1,\ldots ,4$$\end{document} . For the second-order copula models, we set $τ (θ_{g}) = (0.4, 0.5, 0.6, 0.7)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau (\varvec{\theta }_g)=(0.4,0.5,0.6,0.7)$$\end{document} for $g = 1, \dots, 4$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g=1,\ldots ,4$$\end{document} and $τ (δ) = (0.30, 0.35, 0.40, 0.45)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau (\varvec{\delta })=(0.30,0.35,0.40,0.45)$$\end{document} .

The Kendall’s tau parameters $τ (θ_{g})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau (\varvec{\theta }_g)$$\end{document} and $τ (δ_{g})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau (\varvec{\delta }_g)$$\end{document} as described above are common for each group; hence, Table 1 (Table 2) contains the group estimated average biases, root mean square errors (RMSE), and standard deviations (SD), scaled by n, for the IFM (ML) estimates under different pair-copulas from the bi-factor and second-order copula models. In the true (simulated) models, the linking copulas are Gumbel copulas. Given the large number of cutpoints as the number of categories K increases, for ML estimation we restrict ourselves to $K = 3$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$K=3$$\end{document} categories.

Table 1. Small sample of size $n = 500$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n = 500$$\end{document} simulations (10 $^{3}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$^3$$\end{document} replications) from the bi-factor and second-order copula 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.

Table 2. Small sample of size $n = 500$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n = 500$$\end{document} simulations (10 $^{3}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$^3$$\end{document} replications) from the bi-factor and second-order copula models with Gumbel copulas and group estimated average biases, root mean square errors (RMSE), and standard deviations (SD), scaled by n, for the ML estimates under different pair-copulas from the bi-factor and second-order copula models.

Conclusions from the values in the tables 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 $τ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} ’s are not robust under bivariate copula misspecification and their biases increase when the 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 linking copulas are the s.Gumbel copulas.
• IFM is not as punishing for bivariate copula misspecification as ML estimation. For example, the scaled biases for the ML estimates are even larger when the bivariate linking copulas are misspecified to the s.Gumbel copulas.

To examine the reliability of using the heuristic algorithm to select the true (simulated) bivariate linking copulas, samples of size 500 were generated from various bi-factor and second-order copula models. Table 3 presents the number of times that the true (simulated) bivariate linking copulas were chosen over 1000 simulation runs. It is revealed that the model selection algorithm performs extremely well for various bi-factor and second-order copulas models with different choices of linking copulas as the number of categories K increases. For a small K dependence in the tails cannot be easily quantified. Hence, for example, when the true copula is the $t_{5}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_5$$\end{document} which has the same upper and lower tail dependence, the algorithm selected either $t_{5}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_5$$\end{document} or BVN which has zero lower and upper tail dependence, because both copulas provide reflection symmetric dependence.

Table 3. Small sample of size $n = 500$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n = 500$$\end{document} simulations ( $10^{3}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$10^3$$\end{document} replications) from the bi-factor and second-order copula models with various linking copulas and frequencies of the true bivariate copula identified using the model selection algorithm.

To check whether the $χ_{s - q}^{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\chi ^2_{s-q}$$\end{document} is a good approximation for the distribution of the $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistic under the null hypothesis, samples of sizes 500 and 1000 were generated from various bi-factor second-order copula models. Table 4 contains four common nominal levels of the $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistic under the bi-factor and second-order copula models with different bivariate copulas. As can be seen in the table, the observed levels of $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} are close to the nominal $α$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha $$\end{document} levels and remain accurate even for extremely sparse tables ( $d = 16$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d=16$$\end{document} and $K = 5$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$K=5$$\end{document} ).

Table 4. Small sample of size $n = {500, 1000}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n =\{500,1000\}$$\end{document} simulations (10 $^{3}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$^3$$\end{document} replications) from bi-factor and second-order copula models and the empirical rejection levels at $α = {0.20, 0.10, 0.05, 0.01}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha = \{0.20, 0.10, 0.05, 0.01\}$$\end{document} , degrees of freedom (df), mean and variance.

6. Application

The Toronto Alexithymia Scale (TAS) is the most utilized measure of alexithymia in empirical research (e.g., Bagby et al. Reference Bagby, Parker and Taylor1994; Taylor et al. Reference Taylor, Bagby and Parker2003; Parker et al. Reference Parker, Taylor and Bagby2003; Gignac et al. Reference Gignac, Palmer and Stough2007; Reise et al. Reference Reise, Bonifay and Haviland2013; Tuliao et al. Reference Tuliao, Klanecky, Landoy and McChargue2020; Carnovale et al. Reference Carnovale, Taylor, Parker, Sanches and Bagby2021) and is composed of $d = 20$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d=20$$\end{document} items. The aforementioned works suggest that the items measure 3 facets of alexithymia, namely Difficulty Identifying Feelings (DIF; $d_{1} = 7$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d_1=7$$\end{document} items), Difficulty Describing Feelings (DDF; $d_{2} = 5$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d_2=5$$\end{document} items), Externally Oriented Thinking (EOT; $d_{3} = 8$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d_3=8$$\end{document} items), where each facet represents different non-overlapping items. Therefore, a 3-factor model was initially called, but given the hypothesized association among the three facets of the alexithymia construct, oblique (correlated) factor or bi-factor models have been used (e.g., Bagby et al. Reference Bagby, Parker and Taylor1994; Parker et al. Reference Parker, Taylor and Bagby2003; Gignac et al. Reference Gignac, Palmer and Stough2007). Tuliao et al. (Reference Tuliao, Klanecky, Landoy and McChargue2020) has demonstrated that a bi-factor model outperforms any other competing factor model for the TAS scale. To this end, recent studies adopted the bi-factor structure for the TAS scale (e.g., Carnovale et al. Reference Carnovale, Taylor, Parker, Sanches and Bagby2021) and further supported a general alexithymia factor and group-specific factors (DIF, DDf and EOT) that account for homogeneous dependence amongst the non-overlapping items. Note also in passing that estimating a 3-factor or an oblique 3-factor model is computationally demanding as it requires 3-dimensional integration. This is not case for the bi-factor or second-order (an oblique 3-factor model where the group-specific factors are linked to another latent variable via a 1-factor model) models. In spite of the fact they involve $G + 1 = 4$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$G+1=4$$\end{document} latent variables, only require one-dimensional integrals of a function which in turn is a product of 3 one-dimensional integrals.

We use a dataset of 1925 university students from the French-speaking region of Belgium (Briganti and Linkowski Reference Briganti and Linkowski2020). Students were 17 to 25 years old, and 58% of them were female and 42% were male. They were asked to respond to each item using one of $K = 5$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$K=5$$\end{document} categories: “ $1 =$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$1=$$\end{document} completely disagree”, “ $2 =$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$2=$$\end{document} disagree”, “ $3 =$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$3=$$\end{document} neutral, “ $4 =$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$4=$$\end{document} agree”, “ $5 =$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$5=$$\end{document} completely agree”. The dataset and full description of the items can be found in the R package BGGM (Williams and Mulder Reference Williams and Mulder2020).

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 discretized maxima or minima. Since the sample is a mixture (male and female students), we can expect a priori that a bi-factor or second-order copula model with $t_{ν}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_\nu $$\end{document} copulas might be plausible, as in this case the items can be considered as mixtures of discretized means.

In Table 5, we summarize the averages of polychoric semi-correlations for all pairs within each facet of alexithymia and for all pairs of items along with the theoretical semi-correlations in (8) under different choices of copulas. For a BVN/ $t_{ν}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_\nu $$\end{document} copula, the copula parameter is the sample polychoric correlation, while for a Gumbel/s.Gumbel copula the polychoric correlation was converted to Kendall’s tau with the relation in (5) and then from Kendall’s $τ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} to Gumbel/s.Gumbel copula parameter via the functional inverse in (6). The summary of findings from the diagnostics in the table show that the items appear to be a mixed selection between discretized means and minima. For the indicators of the DIF factor (items in group 1) there is more correlation in the joint lower tail, i.e., the items are based on discretizations of latent variables that are minima and have more probability in the joint lower tail, suggesting the use of s.Gumbel linking copulas when modelling the responses to these items. All the other items have stronger correlation in both the joint upper and joint lower tail than with the BVN, i.e., the items are based on discretizations of latent variables that are means and have more probability in both the joint lower and upper tail, suggesting the use of $t_{ν}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_\nu $$\end{document} bivariate linking copulas as the respondents consist of a “mixture” population (different genders), Hence, a bi-factor or second-order copula model with the aforementioned linking copulas might provide a better fit than the (Gaussian) models with BVN copulas.

Table 5. Average observed polychoric correlations and semi-correlations for all pairs within each group and for all pairs of items for the Toronto Alexithymia Scale (TAS), along with the corresponding theoretical semi-correlations for BVN, $t_{5}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_5$$\end{document} , Frank, Gumbel , and survival Gumbel (s.Gumbel) copulas.

Then, we fit the bi-factor and second-order models with the bivariate copulas selected by the heuristic algorithm in Sect. 3.2. For a baseline comparison, we also fit their special cases; these are the one- and two-factor copula models where we have also selected the bivariate copulas using the heuristic algorithm proposed by Kadhem and Nikoloulopoulos (Reference Kadhem and Nikoloulopoulos2021). To show the improvement of the copula models over their Gaussian analogues, we have also fitted all the classes of copula models with BVN copulas. The fitted models are compared via the AIC, since the number of parameters is not the same between the models. In addition, we use Vuong’s test (Vuong Reference Vuong1989) to show if (a) the best fitted model according to the AICs provides better fit than the other fitted models and (b) a model with the selected copulas provides better fit than the one with BVN copulas. The Vuong’s test is the sample version of the difference in Kullback–Leibler divergence between two models and can be used to differentiate two parametric models which could be non-nested. For the Vuong’s test we provide the 95% confidence interval of the test statistic (Joe Reference Joe2014, p. 258). If the interval does not 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_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistic, along with the Root Mean Square Error of Approximation based on $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} (Maydeu-Olivares and Joe Reference Maydeu-Olivares and Joe2014), viz.

\begin{matrix} {RMSEA}_{2} = \sqrt{Max (\frac{M_{2} - d f}{n \times d f}, 0)} . \end{matrix}

Table 6. AICs, Vuong’s 95% CIs, and $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistics for the 1-factor, 2-factor, bi-factor and second-order copula models with BVN copulas and selected copulas, along with the maximum deviations of observed and expected counts for all pairs within each group and for all pairs of items for the Toronto Alexithymia Scale.

$^{a} Selected$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$^\textrm{a}\hbox {Selected}$$\end{document} factor copula model versus its Gaussian special case.

$^{b} Selected$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$^\textrm{b}\hbox {Selected}$$\end{document} Bi-factor copula model versus any other fitted model.

Table 6 gives the AICs, the 95% CIs of Vuong’s tests and the $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistics for all the fitted models. The best 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_{3}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_3$$\end{document} for both the DDF and EOT factors and $t_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_2$$\end{document} for the common factor (alexithymia). This is in line with the preliminary analyses based on the interpretations of items as mixtures of means and the diagnostics in Table 5. The proposed model selection algorithm has selected the $t_{ν}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_\nu $$\end{document} copula that has the same lower and upper tail dependence for the common and all the group specific factors except the group specific factor in group 1 for which the survival Gumbel copula has been selected. For the items in group 1, the largest difference ${\hat{ρ}}_{N}^{-} - {\hat{ρ}}_{N}^{+} = 0.07$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{\rho }}_N^--{\hat{\rho }}_N^+=0.07$$\end{document} is found showing more dependence in the lower tail and the fact that for this group the survival Gumbel copula is selected it shows that this difference is statistically significant. No other difference is statistically significant. It is revealed that the DIF 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_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistics, it is shown that factor copula models provide a big improvement over their Gaussian analogues and that the selected bi-factor copula model outperforms all the fitted models.

The highly statistical significant $M_{2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistics are not surprising since one should expect discrepancies between the postulated parametric model and the population probabilities, when the sample size or dimension is sufficiently large (Maydeu-Olivares and Joe Reference Maydeu-Olivares and Joe2014) as in our case; none should expect a model with 3000 df to fit exactly. To further show that the fit has been improved we have calculated the maximum deviations of observed and model-based counts for each bivariate margin, that is, $D_{j_{1} j_{2}} = n {max}_{y_{1}, y_{2}} | p_{j_{1}, j_{2}, y_{1}, y_{2}} - π_{j_{1}, j_{2}, y_{1}, y_{2}} (\hat{θ}) |$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$D_{j_1j_2}=n\max _{y_1,y_2}|p_{j_1,j_2,y_1,y_2}-\pi _{j_1,j_2,y_1,y_2}({\hat{\varvec{\theta }}})|$$\end{document} . In Table 6, we summarize the averages of these deviations for all pairs within each group and for all pairs of items. Overall, the maximum discrepancies have been sufficiently reduced in the selected bi-factor model.

Table 7 gives the copula parameter estimates in Kendall’s $τ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} scale and their standard errors (SE) for the selected bi-factor copula model and the Gaussian bi-factor model as the benchmark model. The SEs of the estimated 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 estimation of cutpoints can be obtained by jackknifing the two-stage estimation procedure. The loading parameters ( $\hat{τ}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{\tau }}$$\end{document} ’s converted to BVN copula parameters via the functional inverse in (5) and then to loadings using the relations in Section 2.3) show that the common alexithymia factor is mostly loaded on DIF and DDF items, suggesting that items in the domains DIF and DDF are good indicators for alexithymia. The items in the EOT although they loaded on the EOT latent factor, they had poor loadings in the common alexithymia factor.

Table 7. Estimated copula parameters and their standard errors (SE) in Kendall’s $τ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} scale for the bi-factor copula models with BVN copulas and selected copulas for the Toronto Alexithymia Scale.

7. Discussion

For items from several domains, we have proposed bi-factor and second-order copula models where we replace BVN distributions, between observed and latent variables, with bivariate copulas. Our copula constructions include the Gaussian bi-factor and second-order models as special cases and can provide a substantial improvement over the Gaussian models based on AIC, Vuong’s and goodness-of-fit statistics. Hence, superior statistical inference for the loadings can be achieved. We have demonstrated that the Kendall’s $τ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} ’s or loading parameters are not robust under bivariate linking copula misspecification and their biases increase when the assumed bivariate copula has tail dependence of opposite direction from the true bivariate copula.

The improvement relies on the fact that when we use appropriate bivariate copulas other than BVN copulas in the construction, there is an interpretation of latent variables that can be maxima/minima or mixture of means instead of means. The bi-factor and second-order copula models, if other than BVN copulas are called, have a latent structure that is not additive as in (3) and (4), respectively. The bi-factor copula (dependence) parameters are interpretable as dependence of an observed variable with the common factor, or conditional dependence of an observed variable with the group-specific latent variable given the common factor, i.e., the bi-factor copula model permits conditional dependence within identified subsets of items.

Both the bi-factor and second-order copula models lead to substantial simplification of the joint likelihood. The joint pmfs in (1) and (2) reduce to one-dimensional integrals of a function which in turn is a product of G one-dimensional integrals. Hence, the evaluation of the joint likelihood requires only low-dimensional integration, as in the one- and two-factor copula models, regardless of the dimension $G + 1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$G+1$$\end{document} of the factors. This is an advantage over the p-factor ( $p > 2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$p>2$$\end{document} ) copula models where the joint pmf requires p-dimensional integration and becomes intractable as the number of factors increases. Hence, the proposed structured multidimensional factor models provide parsimonious factor solutions without any computational deficiencies as in the p-factor copula models when p increases.

We have proposed a numerically stable likelihood estimation technique based on Gauss–Legendre quadrature. The use of independent Gauss–Legendre quadrature points for this kind of models has been proposed in Krupskii and Joe (Reference Krupskii and Joe2015). For the bi-factor copula models for item response the integrants are bounded and thus independent Gauss–Legendre points are fine. For the second-order copula models for item response the integrants can be unbounded as copula densities can be unbounded, hence we have proposed the novel use of dependent Gauss–Legendre quadrature points that have an 1-factor copula distribution.

Building on the models proposed in this paper, there are several extensions that can be implemented. The adoption of the structure of the Gaussian tri-factor and third-order models (e.g., Rijmen et al. Reference Rijmen, Jeon, von Davier and Rabe-Hesketh2014), to account for any additional layer of dependence, is feasible using the notion of truncated vine copulas that involve both observed and latent variables.

Software

R functions for estimation, simulation, model selection and goodness-of-fit of the bi-factor and second-order copula models are part of the R package FactorCopula (Kadhem and Nikoloulopoulos 2022). All the analyses presented in Sect. 6 are given as code examples in the package.

Acknowledgements

We would like to thank the associate editor and referees for their careful reading and insightful comments that led to an improved presentation. The simulations presented in this paper were carried out on the High Performance Computing Cluster supported by the Research and Specialist Computing Support service at the University of East Anglia.

Footnotes

Supplementary Information The online version contains supplementary material available at https://doi.org/10.1007/S0033312300006050.

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

Bagby, R., Parker, J. D., Taylor, G. J., (1994). The twenty-item Toronto Alexithymia Scale-I. Item selection and cross-validation of the factor structure. Journal of Psychosomatic Research 38(1) 23–32 8126686 10.1016/0022-3999(94)90005-1CrossRef Google Scholar PubMed

Brechmann, E. C., Czado, C., Aas, K., (2012). Truncated regular vines in high dimensions with applications to financial data Canadian Journal of Statistics 40(1) 68–85 10.1002/cjs.10141CrossRef Google Scholar

Briganti, G., Linkowski, P., (2020). Network approach to items and domains from the Toronto Alexithymia Scale Psychological Reports 123(5) 2038–2052 31752608 10.1177/0033294119889586CrossRef Google Scholar PubMed

Carnovale, M., Taylor, G. J., Parker, J. D. A., Sanches, M., & Bagby, R. M. (2021). A bifactor analysis of the 20-item Toronto Alexithymia Scale: Further support for a general alexithymia factor. Psychological Assessment.CrossRef Google Scholar

de la Torre, J., Song, H., (2009). Simultaneous estimation of overall and domain abilities: A higher-order IRT model approach Applied Psychological Measurement 33(8) 620–639 10.1177/0146621608326423CrossRef Google Scholar

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 10.1111/j.1745-3984.2006.00010.xCrossRef Google Scholar

Dissmann, J., Brechmann, E., Czado, C., Kurowicka, D., (2013). Selecting and estimating regular vine copulae and application to financial returns Computational Statistics & Data Analysis 59 52–69 10.1016/j.csda.2012.08.010CrossRef Google Scholar

Gibbons, R. D., Bock, R. D., Hedeker, D., Weiss, D. J., Segawa, E., Bhaumik, D. K., Kupfer, D. J., Frank, E., Grochocinski, V. J., Stover, A., (2007). Full-information item bifactor analysis of graded response data Applied Psychological Measurement 31(1) 4–19 10.1177/0146621606289485CrossRef Google Scholar

Gibbons, R. D., Hedeker, D. R., (1992). Full-information item bi-factor analysis Psychometrika 57(3) 423–436 10.1007/BF02295430CrossRef Google Scholar

Gignac, G. E., Palmer, B. R., 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 18001225 10.1080/00223890701629730CrossRef Google Scholar PubMed

Gustafsson, J-E Balke, G., (1993). General and specific abilities as predictors of school achievement Multivariate Behavioral Research 28(4) 407–434 26801141 10.1207/s15327906mbr2804_2CrossRef Google Scholar PubMed

Joe, H., Multivariate models and dependence concepts London Chapman & Hall 10.1201/b13150Google Scholar

Joe, H., (1997). Asymptotic efficiency of the two-stage estimation method for copula-based models Journal of Multivariate Analysis (2005). 94 401–419 10.1016/j.jmva.2004.06.003CrossRef Google Scholar

Joe, H., Kurowicka, D., Joe, H., (2011). Tail dependence in vine copulae Dependence modeling: Vine copula handbook World Scientific 165–187Google Scholar

Joe, H., Dependence modelling with copulas Chapman & Hall 10.1201/b17116Google Scholar

Joe, H., Li, H., Nikoloulopoulos, A. K., (2014). Tail dependence functions and vine copulas Journal of Multivariate Analysis (2010). 101(1) 252–270 10.1016/j.jmva.2009.08.002CrossRef Google Scholar

Joe, H. & Xu, J. J. (1996). The estimation method of inference functions for margins for multivariate models. Technical Report #166, Department of Statistics, University of British Columbia.Google Scholar

Jöreskog, K. G., Moustaki, I., (2001). Factor analysis of ordinal variables: A comparison of three approaches Multivariate Behavioral Research 36(3) 347–387 26751181 10.1207/S15327906347-387CrossRef Google Scholar PubMed

Kadhem, S. H., Nikoloulopoulos, A. K., (2021). Factor copula models for mixed data British Journal of Mathematical and Statistical Psychology 74(3) 365–403 34626487 10.1111/bmsp.12231CrossRef Google Scholar PubMed

Kadhem and Nikoloulopoulos. (2022). FactorCopula: Factor, bi-factor and second-order copulamodels. R package version 0.8.1. http://CRAN.665R-project.org/package=FactorCopula Google Scholar

Krupskii, P., Joe, H., (2013). Factor copula models for multivariate data Journal of Multivariate Analysis 120 85–101 10.1016/j.jmva.2013.05.001CrossRef Google Scholar

Krupskii, P., Joe, H., (2015). Structured factor copula models: Theory, inference and computation Journal of Multivariate Analysis 138 53–73 10.1016/j.jmva.2014.11.002CrossRef Google Scholar

Lee, G., Frisbie, D. A., (1999). Estimating reliability under a generalizability theory model for test scores composed of testlets Applied Measurement in Education 12(3) 237–255 10.1207/S15324818AME1203_2CrossRef Google Scholar

Maydeu-Olivares, A., (2006). Limited information estimation and testing of discretised multivariate normal structural models Psychometrika 71 57–77 10.1007/s11336-005-0773-4CrossRef Google Scholar

Maydeu-Olivares, A., Joe, H., (2006). Limited information goodness-of-fit testing in multidimensional contingency tables Psychometrika 71(4) 713–732 10.1007/s11336-005-1295-9CrossRef Google Scholar

Maydeu-Olivares, A., Joe, H., (2014). Assessing approximate fit in categorical data analysis Multivariate Behavioral Research 49(4) 305–328 26765800 10.1080/00273171.2014.911075CrossRef Google Scholar

McDonald, R., Test theory: A unified treatment Lawrence Erlbaum Associates IncCrossRef Google Scholar

McDonald, R. P., van der Linden, W. J., Hambleton, R. K., (1999). Normal ogive multidimensional model Handbook of modern item response theory (1997). SpringerGoogle Scholar

McNeil, A. J., Frey, R., Embrechts, P., Quantitative risk management: Concepts, techniques and tools Princeton University PressGoogle Scholar

Mulaik, S. A., Quartetti, D. A., (2005). First order or higher order general factor? Structural Equation Modeling: A Multidisciplinary Journal (1997). 4(3) 193–211 10.1080/10705519709540071CrossRef Google Scholar

Nash, J., Compact numerical methods for computers: Linear algebra and function minimisation 2 HilgerCrossRef Google Scholar

Nikoloulopoulos, A. K., Joe, H., (1990). Factor copula models for item response data Psychometrika (2015). 80(1) 126–150 24297437 10.1007/s11336-013-9387-4CrossRef Google Scholar

Nikoloulopoulos, A. K., Joe, H., Li, H., (2012). Vine copulas with asymmetric tail dependence and applications to financial return data Computational Statistics & Data Analysis 56 3659–3673 10.1016/j.csda.2010.07.016CrossRef Google Scholar

Nikoloulopoulos, A. K., Karlis, D., (2008). Copula model evaluation based on parametric bootstrap Computational Statistics & Data Analysis 52 3342–3353 10.1016/j.csda.2007.10.028CrossRef Google Scholar

Parker, J., Taylor, G., Bagby, R., (2003). The 20-item Toronto alexithymia scale: III. Reliability and factorial validity in a community population Journal of Psychosomatic Research 55 269–275 12932802 10.1016/S0022-3999(02)00578-0CrossRef Google Scholar

Reise, S. P., Bonifay, W. E., Haviland, M. G., (2013). Scoring and modeling psychological measures in the presence of multidimensionality Journal of Personality Assessment 95(2) 129–140 23030794 10.1080/00223891.2012.725437CrossRef Google Scholar PubMed

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 10.1111/j.1745-3984.2010.00118.xCrossRef Google Scholar

Rijmen, F., Jeon, M., von Davier, M., Rabe-Hesketh, S., (2014). A third-order item response theory model for modeling the effects of domains and subdomains in large-scale educational assessment surveys Journal of Educational and Behavioral Statistics 39(4) 235–256 10.3102/1076998614531045CrossRef Google Scholar

Samejima, F., (1969). Estimation of latent ability using a response pattern of graded scores Psychometrika 34(1) 1–97 10.1007/BF03372160CrossRef Google Scholar

Schmid, J., Leiman, J. M., (1957). The development of hierarchical factor solutions Psychometrika 22(1) 53–61 10.1007/BF02289209CrossRef Google Scholar

Sireci, S. G., Thissen, D., Wainer, H., (1991). On the reliability of testlet-based tests Journal of Educational Measurement 28(3) 237–247 10.1111/j.1745-3984.1991.tb00356.xCrossRef Google Scholar

Sklar, A., Fonctions de répartition à

n

dimensions et leurs marges Publications de l’Institut de Statistique de l’Université de Paris 8 229–231Google Scholar

Stroud, A., Secrest, D., Gaussian quadrature formulas (1966). Prentice-HallGoogle Scholar

Taylor, G., Bagby, R., Parker, J., (1959). The 20-item Toronto Alexithymia Scale IV. Reliability and factorial validity in different languages and cultures Journal of Psychosomatic Research (2003). 55 277–283 12932803 10.1016/S0022-3999(02)00601-3CrossRef Google Scholar

Tuliao, A. P., Klanecky, A. K., Landoy, B. V. N., 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 30661362 10.1177/1073191118824030CrossRef Google Scholar

Vuong, Q. H., (1989). Likelihood ratio tests for model selection and non-nested hypotheses Econometrica 57(2) 307–333 10.2307/1912557CrossRef Google Scholar

Wainer, H., Kiely, G. L., (1987). Item clusters and computerized adaptive testing: A case for testlets Journal of Educational Measurement 24(3) 185–201 10.1111/j.1745-3984.1987.tb00274.xCrossRef Google Scholar

Wainer, H., Thissen, D., (1996). How is reliability related to the quality of test scores? What is the effect of local dependence on reliability? Educational Measurement: Issues and Practice 15(1) 22–29 10.1111/j.1745-3992.1996.tb00803.xCrossRef Google Scholar

Wang, W-C Wilson, M., (2005). The rasch testlet model Applied Psychological Measurement 29(2) 126–149 10.1177/0146621604271053CrossRef Google Scholar

Williams, D. & Mulder, J. (2020). BGGM: Bayesian Gaussian Graphical Models. R package version 1.0.0.CrossRef Google Scholar

Xu, J. (1996). Statistical modelling and inference for multivariate and longitudinal discrete response. The University of British Columbia, Ph.D. thesis.Google Scholar

Yung, Y-F Thissen, D., McLeod, L. D., (1999). On the relationship between the higher-order factor model and the hierarchical factor model Psychometrika 64(2) 113–128 10.1007/BF02294531CrossRef Google Scholar

Zenisky, A. L., Hambleton, R. K., Sireci, S. G., (2002). Identification and evaluation of local item dependencies in the medical college admissions test Journal of Educational Measurement 39(4) 291–309 10.1111/j.1745-3984.2002.tb01144.xCrossRef Google Scholar

Figure 1. Graphical representation of the bi-factor copula model with G group-specific factors and a common factor V0\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document}.

Figure 2. Graphical representation of the second-order copula model with G first-order factors and one second-order factor V0\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_0$$\end{document}.

Figure 3. Contour plots of bivariate copulas with standard normal margins and dependence parameters corresponding to Kendall’s τ\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} value of 0.6 in absolute value.

Table 1. Small sample of size n=500\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n = 500$$\end{document} simulations (103\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$^3$$\end{document} replications) from the bi-factor and second-order copula 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.

Table 2. Small sample of size n=500\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n = 500$$\end{document} simulations (103\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$^3$$\end{document} replications) from the bi-factor and second-order copula models with Gumbel copulas and group estimated average biases, root mean square errors (RMSE), and standard deviations (SD), scaled by n, for the ML estimates under different pair-copulas from the bi-factor and second-order copula models.

Table 3. Small sample of size n=500\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n = 500$$\end{document} simulations (103\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$10^3$$\end{document} replications) from the bi-factor and second-order copula models with various linking copulas and frequencies of the true bivariate copula identified using the model selection algorithm.

Table 4. Small sample of size n={500,1000}\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n =\{500,1000\}$$\end{document} simulations (103\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$^3$$\end{document} replications) from bi-factor and second-order copula models and the empirical rejection levels at α={0.20,0.10,0.05,0.01}\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha = \{0.20, 0.10, 0.05, 0.01\}$$\end{document}, degrees of freedom (df), mean and variance.

Table 5. Average observed polychoric correlations and semi-correlations for all pairs within each group and for all pairs of items for the Toronto Alexithymia Scale (TAS), along with the corresponding theoretical semi-correlations for BVN, t5\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t_5$$\end{document}, Frank, Gumbel , and survival Gumbel (s.Gumbel) copulas.

Table 6. AICs, Vuong’s 95% CIs, and M2\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_2$$\end{document} statistics for the 1-factor, 2-factor, bi-factor and second-order copula models with BVN copulas and selected copulas, along with the maximum deviations of observed and expected counts for all pairs within each group and for all pairs of items for the Toronto Alexithymia Scale.

Table 7. Estimated copula parameters and their standard errors (SE) in Kendall’s τ\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tau $$\end{document} scale for the bi-factor copula models with BVN copulas and selected copulas for the Toronto Alexithymia Scale.

Kadhem and Nikoloulopoulos supplementary material

Supplementary Information The online version contains supplementary material available at https://doi.org/10.1007/S0033312300006050.

File 74.5 KB

Article contents

Bi-factor and Second-Order Copula Models for Item Response Data

Abstract

Keywords

Information

1. Bi-factor and Second-Order Copula Models

1.1. Bi-factor Copula Model

1.2. Second-Order Copula Model

1.3. Special Cases

1.4. Other Choices of Parametric Bivariate Copulas

2. Estimation and Computational Details

3. Bivariate Copula Selection

3.1. Simple Diagnostics Based on Semi-correlations

3.2. Selection Algorithm

4. Goodness of Fit

5. Simulations

6. Application

7. Discussion

Software

Acknowledgements

Footnotes

References

Kadhem and Nikoloulopoulos supplementary material

Save article to Kindle

Save article to Dropbox

Save article to Google Drive

Reply to: Submit a response

Your details

You have entered the maximum number of contributors

Conflicting interests