Factor copula models for mixed data

TL;DR

This paper introduces factor copula models for analyzing complex dependence structures in mixed continuous and discrete data, offering improved flexibility over traditional models through the use of vine copulas with latent variables.

Contribution

It develops a novel factor copula modeling framework that captures tail, asymmetric, and non-linear dependence, with methods for model selection and goodness-of-fit.

Findings

Substantial improvement over standard factor models for mixed data.

Demonstrated effectiveness through simulation studies.

Applied to real datasets showing practical utility.

Abstract

We develop factor copula models for analysing the dependence among mixed continuous and discrete responses. Factor copula models are canonical vine copulas that involve both observed and latent variables, hence they allow tail, asymmetric and non-linear dependence. They can be explained as conditional independence models with latent variables that don't necessarily have an additive latent structure. We focus on important issues that would interest the social data analyst, such as model selection and goodness-of-fit. Our general methodology is demonstrated with an extensive simulation study and illustrated by re-analysing three mixed response datasets. Our study suggests that there can be a substantial improvement over the standard factor model for mixed data and makes the argument for moving to factor copula models.

Tables

Table 1. Table 1 : Lower semi-correlations ρ N − superscript subscript 𝜌 𝑁 \rho_{N}^{-} , upper semi-correlations ρ N + superscript subscript 𝜌 𝑁 \rho_{N}^{+} , lower tail dependence λ L subscript 𝜆 𝐿 \lambda_{L} , and upper tail dependence λ U subscript 𝜆 𝑈 \lambda_{U} , with τ = { 0.3 , 0.5 , 0.7 } 𝜏 0.3 0.5 0.7 \tau=\{0.3,0.5,0.7\} for 1-parameter and 2-parameter bivariate copulas.

Bivariate copula	$\tau$		$\theta$	$\delta$		${\rho }_N^{-}$	${\rho }_N^{+}$		${\lambda }_L$	${\lambda }_U$
BVN	0.3		0.45			0.23	0.23		0.00	0.00
	0.5		0.71			0.47	0.47		0.00	0.00
	0.7		0.89			0.75	0.75		0.00	0.00
$t_3$	0.3		0.45			0.45	0.45		0.29	0.29
	0.5		0.71			0.61	0.61		0.45	0.45
	0.7		0.89			0.80	0.80		0.66	0.66
Frank	0.3		2.92			0.15	0.15		0.00	0.00
	0.5		5.74			0.32	0.32		0.00	0.00
	0.7		11.41			0.60	0.60		0.00	0.00
Joe	0.3		1.77			0.05	0.58		0.00	0.52
	0.5		2.86			0.14	0.78		0.00	0.73
	0.7		5.46			0.37	0.92		0.00	0.86
Gumbel	0.3		1.43			0.16	0.46		0.00	0.38
	0.5		2.00			0.36	0.67		0.00	0.59
	0.7		3.33			0.64	0.85		0.00	0.77
BB1	0.3		0.50	1.14		0.43	0.25		0.30	0.17
	0.5		0.35	1.71		0.52	0.59		0.31	0.50
	0.7		1.33	2.00		0.85	0.72		0.77	0.59
BB7	0.3		1.40	0.40		0.28	0.37		0.18	0.36
	0.5		1.50	1.57		0.66	0.42		0.64	0.41
	0.7		4.00	2.00		0.73	0.85		0.71	0.81
BB8	0.3		3.92	0.60		0.10	0.22		0.00	0.00
	0.5		4.51	0.80		0.20	0.52		0.00	0.00
	0.7		6.89	0.90		0.41	0.84		0.00	0.00
BB10	0.3		1.60	0.83		0.18	0.09		0.00	0.00
	0.5		2.50	0.98		0.43	0.19		0.00	0.00
	0.7		10.00	1.00		0.25	0.66		0.00	0.00

Table 2. Table 2 : The sample correlation ρ N subscript 𝜌 𝑁 \rho_{N} , lower semi-correlation ρ N − superscript subscript 𝜌 𝑁 \rho_{N}^{-} , and upper semi-correlation ρ N + superscript subscript 𝜌 𝑁 \rho_{N}^{+} for each pair of variables, along with the measures of discrepancy between the sample and the resulting correlation matrix of linear factor analysis with 1 and 2 factors for the political-economic risk data.

pairs of variables		${\rho }_N$	${\rho }_N^{-}$	${\rho }_N^{+}$
BM	GDP	0.53	-0.04	0.57
BM	IJ	0.61	-	-
BM	XPR	0.67	0.88	0.63
BM	CRP	0.62	0.16	0.55
GDP	IJ	0.78	-	-
GDP	XPR	0.55	0.11	0.75
GDP	CRP	0.77	0.24	0.63
IJ	XPR	0.91	-	-
IJ	CRP	0.87	-	-
XPR	CRP	0.76	0.71	0.71

Table 3. Table 3 : Estimated parameters, their standard errors (SE) in Kendall’s τ 𝜏 \tau scale, joint log-likelihoods, the 95 % percent 95 95\% CIs of Vuong’s statistics, and the M 2 subscript 𝑀 2 M_{2} statistics for the one-factor copula models for the political-economic risk data.

1-factor	BVN¶			$t_5$			Frank			Selected
	$\hat{\tau }$	SE		$\hat{\tau }$	SE		$\hat{\tau }$	SE		copulas	$\hat{\tau }$	SE
BM	0.50	0.06		0.51	0.07		0.49	0.06		Joe	0.51	0.05
GDP	0.57	0.05		0.57	0.06		0.58	0.06		Joe	0.58	0.05
IJ	0.80	0.09		0.81	0.09		0.75	0.09		reflected Joe	0.80	0.07
XPR	0.66	0.06		0.68	0.07		0.66	0.06		Joe	0.69	0.06
CRP	0.71	0.06		0.70	0.06		0.72	0.06		Gumbel	0.74	0.06
$\ell$	-165.15			-166.25			-164.89			-151.98
Vuong $95\%$CI				(-0.051,0.015)			(-0.077,0.085)			(0.073,0.352)
$M_2$	179.2			187.4			177.6			129.2
df	134			134			134			134
$p$-value										0.60

Table 4. Table 4 : The sample correlation ρ N subscript 𝜌 𝑁 \rho_{N} , lower semi-correlation ρ N − superscript subscript 𝜌 𝑁 \rho_{N}^{-} , and upper semi-correlation ρ N + superscript subscript 𝜌 𝑁 \rho_{N}^{+} for each pair of variables, along with the measures of discrepancy between the sample and the resulting correlation matrix of linear factor analysis with 1, 2 and 3 factors for the general social survey dataset.

pairs of variables		${\rho }_N$	${\rho }_N^{-}$	${\rho }_N^{+}$
income	age	0.29	0.48	0.23
income	degree	0.52	0.24	0.33
income	pincome	0.14	0.02	0.28
income	pdegree	0.24	0.04	0.08
income	child	0.22	0.23	0.01
income	pchild	-0.09	0.06	0.00
age	degree	0.06	0.22	-0.04
age	pincome	-0.11	-0.02	0.12
age	pdegree	-0.14	-0.42	0.44
age	child	0.58	0.36	0.26
age	pchild	0.12	0.18	0.07
degree	pincome	0.21	0.17	-0.05
degree	pdegree	0.46	0.46	0.41
degree	child	-0.11	-0.10	-0.09
degree	pchild	-0.25	-0.14	-0.30
pincome	pdegree	0.44	0.44	0.34
pincome	child	-0.16	-0.15	0.11
pincome	pchild	-0.23	0.13	-0.30
pdegree	child	-0.21	0.08	0.10
pdegree	pchild	-0.34	0.19	-0.32
child	pchild	0.20	-0.11	-0.06

Table 5. Table 5 : Estimated parameters, their standard errors (SE) in Kendall’s τ 𝜏 \tau scale, joint log-likelihoods, the 95 % percent 95 95\% CIs of Vuong’s statistics, and the M 2 subscript 𝑀 2 M_{2} statistics for the 1- and 2-factor copula models for the general social survey dataset.

1-factor	BVN¶			$t_9$			Frank			Selected
	$\hat{\tau }$	SE		$\hat{\tau }$	SE		$\hat{\tau }$	SE		copulas	$\hat{\tau }$	SE
income	0.20	0.04		0.20	0.04		0.20	0.04		Joe	0.29	0.04
age	-0.14	0.04		-0.14	0.04		-0.14	0.04		2-reflected Joe	-0.14	0.03
degree	0.40	0.04		0.39	0.04		0.38	0.04		$t_3$	0.45	0.04
pincome	0.33	0.03		0.34	0.04		0.35	0.04		$t_3$	0.33	0.05
pdegree	0.62	0.05		0.65	0.05		0.68	0.06		reflected Gumbel	0.56	0.05
child	-0.20	0.04		-0.19	0.04		-0.19	0.04		2-reflected Joe	-0.14	0.03
pchild	-0.32	0.03		-0.31	0.04		-0.32	0.04		2-reflected Gumbel	-0.27	0.03
$\ell$	-3425.39			-3420.56			-3433.83			-3397.79
Vuong $95\%$CI				(-0.005,-0.025)			(-0.037,0.001)			(0.022,0.097)
$M_2$	743.74			715.45			738.76			660.47
df	348			348			348			348
$p$-value

Table 6. Table 6 : The sample correlation ρ N subscript 𝜌 𝑁 \rho_{N} , lower semi-correlation ρ N − superscript subscript 𝜌 𝑁 \rho_{N}^{-} , and upper semi-correlation ρ N + superscript subscript 𝜌 𝑁 \rho_{N}^{+} for each pair of variables, along with the measures of discrepancy between the sample and the resulting correlation matrix of linear factor analysis with 1, 2 and 3 factors for the Swiss consumption survey dataset.

pairs of variables		${\rho }_N$	${\rho }_N^{-}$	${\rho }_N^{+}$
food	clothes	0.65	0.21	0.76
food	leisure	0.60	0.18	0.76
food	dishwasher	0.31	-	-
food	car	0.38	-	-
food	motorcycle	0.11	-	-
food	bicycles	0.21	0.22	0.02
clothes	leisure	0.52	0.02	0.63
clothes	dishwasher	0.23	-	-
clothes	car	0.25	-	-
clothes	motorcycle	0.07	-	-
clothes	bicycles	0.18	0.15	0.02
leisure	dishwasher	0.24	-	-
leisure	car	0.18	-	-
leisure	motorcycle	0.01	-	-
leisure	bicycles	0.08	0.04	0.08
dishwasher	car	0.43	-	-
dishwasher	motorcycle	0.03	-	-
dishwasher	bicycles	0.24	-	-
car	motorcycle	0.18	-	-
car	bicycles	0.26	-	-
motorcycle	bicycles	0.21	-	-

Table 7. Table 7 : Estimated parameters, their standard errors (SE) in Kendall’s τ 𝜏 \tau scale, joint log-likelihoods, the 95 % percent 95 95\% CIs of Vuong’s statistics, and the M 2 subscript 𝑀 2 M_{2} statistics for the 1- and 2-factor copula models for the Swiss consumption survey dataset.

1-factor		BVN¶			$t_5$			Frank			Selected
		$\hat{\tau }$	SE		$\hat{\tau }$	SE		$\hat{\tau }$	SE		Copulas	$\hat{\tau }$	SE
food		0.69	0.01		0.73	0.01		0.74	0.01		reflected BB10	0.79	0.00
clothes		0.53	0.01		0.53	0.01		0.53	0.01		BB10	0.38	0.00
leisure		0.47	0.01		0.50	0.01		0.50	0.01		BB10	0.39	0.00
dishwasher		0.24	0.01		0.25	0.01		0.23	0.01		reflected Joe	0.28	0.01
car		0.27	0.01		0.30	0.01		0.28	0.01		reflected Joe	0.23	0.01
motorcycle		0.07	0.01		0.06	0.01		0.08	0.01		reflected Joe	0.13	0.01
bicycles		0.15	0.01		0.15	0.01		0.16	0.01		reflected Joe	0.17	0.01
AIC		55004.24			54221.36			55105.88			48932.32
Vuong $95\%$ CI					(0.032,0.046)			(-0.015,0.005)			(0.286,0.324)
$M_2$		2775.73			2734.05			2808.53			1626.54
df		71			71			71			68
$p-$value

Table 8. Table 8 : Maximum deviations D j 1 ​ j 2 subscript 𝐷 subscript 𝑗 1 subscript 𝑗 2 D_{j_{1}j_{2}} of observed and expected counts for each bivariate margin ( j 1 , j 2 ) subscript 𝑗 1 subscript 𝑗 2 (j_{1},j_{2}) for the 1- and 2-factor copula models for the Swiss consumption survey dataset.

	1-factor model					2-factor model
$D_{j_1,j_2}$	BVN	$t_5$	Frank	Selected		BVN	$t_7$	Frank	Selected
$D_{1,2}$	347	317	303	167		349	311	270	40
$D_{1,3}$	511	468	456	183		509	460	428	133
$D_{1,4}$	158	177	163	70		159	185	161	56
$D_{1,5}$	231	189	223	119		233	181	230	60
$D_{1,6}$	87	117	88	60		87	130	72	12
$D_{1,7}$	78	92	79	88		78	110	89	81
$D_{2,3}$	442	418	431	69		433	403	393	54
$D_{2,4}$	59	80	84	145		38	56	64	86
$D_{2,5}$	96	107	107	201		60	47	93	36
$D_{2,6}$	18	3	18	27		19	15	29	39
$D_{2,7}$	51	76	60	83		49	91	52	61
$D_{3,4}$	182	146	141	196		253	216	168	83
$D_{3,5}$	82	105	106	191		59	13	83	61
$D_{3,6}$	59	58	69	71		13	23	27	45
$D_{3,7}$	62	54	64	103		65	67	69	59
$D_{4,5}$	289	276	286	223		66	74	207	2
$D_{4,6}$	9	5	11	29		133	138	100	96
$D_{4,7}$	82	81	81	88		28	20	46	54
$D_{5,6}$	111	123	111	77		15	22	19	20
$D_{5,7}$	101	96	95	68		33	25	40	64
$D_{6,7}$	70	74	70	61		80	96	87	52

Table 9. Table 9 : Small sample of sizes n = { 100 , 300 , 500 } 𝑛 100 300 500 n=\{100,300,500\} simulations ( 10 4 superscript 10 4 10^{4} replications) from the selected factor copula models in Section 6 to assess the measures of discrepancy D 1 subscript 𝐷 1 D_{1} , D 2 subscript 𝐷 2 D_{2} , and D 3 subscript 𝐷 3 D_{3} between the observed and the resulting correlation matrix of linear factor analysis for 1, 2 and 3 factors, with resultant means and standard deviations (SD).

Political-economic dataset – 1-factor model
			$D_1$			$D_2$			$D_3$
$n$	$\#$ factors		mean	SD		mean	SD		mean	SD
100	1		0.061	0.027		0.016	0.006		0.101	0.071
	2		0.022	0.016		0.004	0.003		0.014	0.023
300	1		0.038	0.017		0.010	0.004		0.036	0.023
	2		0.011	0.008		0.002	0.002		0.004	0.005
500	1		0.033	0.014		0.009	0.003		0.024	0.015
	2		0.009	0.006		0.002	0.001		0.002	0.003

Table 10. Table 10 : Small sample of sizes n = { 100 , 300 , 500 } 𝑛 100 300 500 n=\{100,300,500\} simulations ( 10 4 superscript 10 4 10^{4} replications) from the Gumbel copula with Kendall’s τ = { 0.3 , 0.5 , 0.7 } 𝜏 0.3 0.5 0.7 \tau=\{0.3,0.5,0.7\} for mixed continuous, ordinal, and count data with resultant biases, root mean square errors (RMSE) and standard deviations (SD), scaled by n 𝑛 n , for the estimated correlation ρ N subscript 𝜌 𝑁 \rho_{N} , lower semi-correlation ρ N − superscript subscript 𝜌 𝑁 \rho_{N}^{-} , and upper semi-correlation ρ N + superscript subscript 𝜌 𝑁 \rho_{N}^{+} .

			(continuous, continuous)				(continuous, ordinal)				(continuous, count)				(ordinal, ordinal)				(ordinal, count)				(count, count)
$n$	$\tau$		${\rho }_N$	${\rho }_N^{-}$	${\rho }_N^{+}$		${\rho }_N$	${\rho }_N^{-}$	${\rho }_N^{+}$		${\rho }_N$	${\rho }_N^{-}$	${\rho }_N^{+}$		${\rho }_N$	${\rho }_N^{-}$	${\rho }_N^{+}$		${\rho }_N$	${\rho }_N^{-}$	${\rho }_N^{+}$		${\rho }_N$	${\rho }_N^{-}$	${\rho }_N^{+}$
100	0.3	True values	0.46	0.16	0.46		0.46	0.16	0.46		0.46	0.16	0.46		0.46	0.16	0.46		0.46	0.16	0.46		0.46	0.16	0.46
		$n$Bias	-1.09	-0.79	-4.18		-1.05	0.16	-9.25		1.62	1.40	-4.41		-0.55	2.35	-9.86		0.62	4.88	-8.25		2.03	9.54	-5.58
		$n$SD	8.57	18.10	16.83		9.03	18.64	16.71		9.23	16.54	20.46		9.31	18.95	18.03		9.37	17.08	21.78		9.55	14.73	24.24
		$n$RMSE	8.64	18.12	17.34		9.09	18.64	19.10		9.37	16.60	20.93		9.32	19.10	20.55		9.39	17.76	23.29		9.76	17.55	24.87
	0.5	True values	0.70	0.36	0.67		0.70	0.36	0.67		0.70	0.36	0.67		0.70	0.36	0.67		0.70	0.36	0.67		0.70	0.36	0.67
		$n$Bias	-0.99	-1.65	-3.98		-0.16	-0.38	-10.21		2.43	0.10	-3.78		0.26	3.93	-8.58		1.41	7.42	-8.18		2.80	14.88	-5.57
		$n$SD	5.77	15.73	11.72		6.26	15.67	12.41		6.19	14.49	14.93		6.30	15.91	13.52		6.35	14.71	16.73		6.34	12.18	17.35
		$n$RMSE	5.85	15.82	12.37		6.26	15.68	16.07		6.65	14.49	15.40		6.31	16.39	16.01		6.51	16.48	18.62		6.94	19.23	18.22
	0.7	True values	0.88	0.64	0.85		0.88	0.64	0.85		0.88	0.64	0.85		0.88	0.64	0.85		0.88	0.64	0.85		0.88	0.64	0.85
		$n$Bias	-0.71	-2.12	-2.74		0.78	-2.21	-8.77		2.16	-5.21	-0.28		0.55	4.34	-4.46		1.23	6.29	-4.75		2.12	13.76	-2.19
		$n$SD	2.71	10.76	5.99		3.02	10.84	7.29		2.80	10.74	8.40		3.07	10.48	7.77		3.03	10.24	10.91		2.94	7.26	9.42
		$n$RMSE	2.80	10.96	6.59		3.12	11.06	11.40		3.53	11.94	8.40		3.12	11.35	8.96		3.27	12.02	11.90		3.63	15.56	9.67
300	0.3	True values	0.46	0.16	0.46		0.46	0.16	0.46		0.46	0.16	0.46		0.46	0.16	0.46		0.46	0.16	0.46		0.46	0.16	0.46
		$n$Bias	-1.44	-1.48	-5.96		-2.88	1.14	-26.54		5.52	4.43	-12.35		-1.56	7.52	-28.59		2.04	14.61	-25.56		6.36	28.76	-16.44
		$n$SD	15.04	30.94	28.32		15.75	31.26	27.83		16.11	28.02	33.34		16.32	32.42	30.34		16.45	28.98	36.06		16.65	25.39	40.28
		$n$RMSE	15.11	30.97	28.94		16.01	31.28	38.46		17.03	28.37	35.55		16.39	33.29	41.69		16.58	32.46	44.20		17.82	38.36	43.50
	0.5	True values	0.70	0.36	0.67		0.70	0.36	0.67		0.70	0.36	0.67		0.70	0.36	0.67		0.70	0.36	0.67		0.70	0.36	0.67
		$n$Bias	-1.23	-2.48	-5.34		-0.77	-1.16	-30.78		7.39	-0.39	-11.03		0.64	11.81	-25.37		4.11	21.74	-25.71		8.39	44.61	-16.40
		$n$SD	9.99	26.98	19.09		10.87	26.60	20.40		10.62	24.72	24.25		11.08	27.48	22.59		11.09	25.22	27.56		10.96	20.84	28.88
		$n$RMSE	10.06	27.09	19.82		10.90	26.63	36.93		12.94	24.72	26.64		11.10	29.91	33.97		11.82	33.30	37.69		13.80	49.24	33.22
	0.7	True values	0.88	0.64	0.85		0.88	0.64	0.85		0.88	0.64	0.85		0.88	0.64	0.85		0.88	0.64	0.85		0.88	0.64	0.85
		$n$Bias	-0.83	-2.93	-3.43		1.42	-7.89	-28.35		5.84	-18.52	-1.43		1.31	12.56	-13.92		3.27	17.32	-16.87		5.97	40.56	-7.09
		$n$SD	4.60	18.37	9.35		5.16	18.37	11.94		4.71	18.05	13.58		5.34	18.16	13.05		5.26	17.59	18.02		5.04	12.35	15.54
		$n$RMSE	4.68	18.61	9.96		5.35	19.99	30.76		7.50	25.86	13.66		5.50	22.08	19.08		6.20	24.68	24.69		7.81	42.40	17.08
500	0.3	True values	0.46	0.16	0.46		0.46	0.16	0.46		0.46	0.16	0.46		0.46	0.16	0.46		0.46	0.16	0.46		0.46	0.16	0.46
		$n$Bias	-1.37	-1.08	-7.04		-4.45	2.42	-44.04		9.65	7.93	-20.25		-2.25	12.60	-47.33		3.75	24.71	-42.32		10.96	48.14	-27.35
		$n$SD	19.06	39.98	36.96		19.95	39.89	35.47		20.49	35.93	42.97		20.75	41.68	39.18		21.00	37.65	46.91		21.35	32.71	52.79
		$n$RMSE	19.11	40.00	37.63		20.44	39.97	56.55		22.64	36.79	47.51		20.87	43.54	61.45		21.33	45.04	63.17		24.00	58.20	59.46
	0.5	True values	0.70	0.36	0.67		0.70	0.36	0.67		0.70	0.36	0.67		0.70	0.36	0.67		0.70	0.36	0.67		0.70	0.36	0.67
		$n$Bias	-1.11	-2.31	-6.27		-1.16	-1.38	-51.40		12.58	-0.39	-18.48		1.34	19.78	-41.78		7.16	36.42	-42.79		14.31	74.64	-27.23
		$n$SD	12.58	35.08	24.56		13.67	34.07	26.18		13.31	31.87	31.24		14.04	35.55	29.22		13.99	32.68	36.00		13.91	26.57	37.57
		$n$RMSE	12.63	35.16	25.35		13.72	34.10	57.68		18.31	31.88	36.29		14.10	40.69	50.99		15.71	48.93	55.92		19.95	79.23	46.40
	0.7	True values	0.88	0.64	0.85		0.88	0.64	0.85		0.88	0.64	0.85		0.88	0.64	0.85		0.88	0.64	0.85		0.88	0.64	0.85
		$n$Bias	-0.76	-2.83	-3.63		2.03	-13.45	-47.97		9.47	-31.71	-2.93		2.21	21.01	-22.91		5.43	28.29	-28.60		9.95	67.47	-11.82
		$n$SD	5.81	23.66	11.69		6.48	23.48	15.30		5.95	23.22	17.23		6.75	23.36	16.73		6.68	22.72	23.21		6.42	15.65	20.03
		$n$RMSE	5.86	23.82	12.24		6.79	27.06	50.35		11.18	39.30	17.48		7.10	31.42	28.37		8.61	36.29	36.84		11.84	69.26	23.26

Table 11. Table 11 : Small sample of sizes n = { 100 , 300 , 500 } 𝑛 100 300 500 n=\{100,300,500\} simulations ( 10 4 superscript 10 4 10^{4} replications) from the selected factor copula models in Section 6 with resultant biases, root mean square errors (RMSE) and standard deviations (SD), scaled by n 𝑛 n , for the estimated parameters.

Political-economic dataset – 1-factor model
$\tau$	0.51				0.58				0.80				0.68				0.74
$n$	100	300	500		100	300	500		100	300	500		100	300	500		100	300	500
$n$Bias	0.88	2.30	3.17		-1.36	-3.39	-4.87		0.75	-0.55	0.64		0.21	-0.27	0.19		0.29	2.57	0.73
$n$SD	4.28	7.60	9.63		4.19	7.50	9.08		5.41	10.91	11.98		4.58	8.43	9.84		4.46	14.92	11.78
$n$RMSE	4.37	7.95	10.13		4.40	8.23	10.31		5.47	10.92	12.00		4.59	8.44	9.84		4.47	15.13	11.80

Table 12. Table 12 : Small sample of sizes n = { 100 , 300 , 500 } 𝑛 100 300 500 n=\{100,300,500\} simulations ( 10 4 superscript 10 4 10^{4} replications) from the selected 1-factor copula models in Section 6 with resultant biases, root mean square errors (RMSE) and standard deviations (SD), scaled by n 𝑛 n , for the estimated parameters under an 1-factor copula model with BVN copulas, i.e. the standard factor model.

Political-economic dataset – 1-factor model
$\tau$	0.51				0.58				0.80				0.68				0.74
$n$	100	300	500		100	300	500		100	300	500		100	300	500		100	300	500
$n$Bias	-0.35	-0.96	-1.56		-1.40	-3.90	-6.41		-1.29	-6.19	-10.54		0.51	0.89	1.40		-0.18	-1.08	-2.15
$n$SD	5.24	9.16	11.57		4.95	8.57	11.13		6.03	9.52	12.14		4.60	7.91	10.01		4.42	7.49	9.69
$n$RMSE	5.25	9.21	11.68		5.15	9.42	12.85		6.17	11.35	16.07		4.63	7.96	10.11		4.42	7.57	9.92

Table 13. Table 13 : Small sample of sizes n = { 100 , 300 , 500 } 𝑛 100 300 500 n=\{100,300,500\} distribution for M 2 subscript 𝑀 2 M_{2} ( 10 4 superscript 10 4 10^{4} replications). Empirical rejection levels at α = { 0.20 , 0.10 , 0.05 , 0.01 } 𝛼 0.20 0.10 0.05 0.01 \alpha=\{0.20,0.10,0.05,0.01\} , degrees of freedom (df), and mean under the factor copula models. Continuous and count variables are transformed to ordinal with K = { 3 , 4 , 5 } 𝐾 3 4 5 K=\{3,4,5\} and K = { 3 , 4 } 𝐾 3 4 K=\{3,4\} categories, respectively, using the general strategies proposed in Section 5.2 . Count variables area also transformed to ordinal with K = 5 𝐾 5 K=5 categories by treating them as ordinal where the 5th category contained all the counts greater than 3.

	$n=100$				$n=300$				$n=500$
	$K=3$	$K=4$	$K=5$		$K=3$	$K=4$	$K=5$		$K=3$	$K=4$	$K=5$
Political-economic dataset – 1-factor model
df	92	121	152		92	121	152		92	121	152
mean	89.3	118.3	148.4		91.0	119.7	152.6		91.0	119.6	152.3
$\alpha =0.20$	0.183	0.192	0.197		0.196	0.194	0.195		0.196	0.189	0.190
$\alpha =0.10$	0.121	0.125	0.134		0.122	0.121	0.119		0.114	0.109	0.109
$\alpha =0.05$	0.083	0.089	0.098		0.076	0.077	0.077		0.072	0.070	0.067
$\alpha =0.01$	0.044	0.046	0.055		0.036	0.034	0.037		0.027	0.030	0.026
General social survey – 1-factor model
df	161	239	329		161	239	329		161	239	329
mean	161.5	240.0	333.0		160.7	239.4	329.7		161.3	240.2	329.6
$\alpha =0.20$	0.213	0.220	0.240		0.202	0.216	0.203		0.211	0.228	0.212
$\alpha =0.10$	0.110	0.121	0.122		0.106	0.118	0.102		0.118	0.127	0.108
$\alpha =0.05$	0.058	0.070	0.061		0.054	0.067	0.051		0.065	0.073	0.056
$\alpha =0.01$	0.013	0.018	0.014		0.014	0.019	0.012		0.016	0.023	0.011
Swiss consumption survey – 1-factor model
df	74	128	194		74	128	194		74	128	194
mean	75.4	130.1	197.8		74.6	128.5	195.1		74.5	128.0	194.4
$\alpha =0.20$	0.229	0.239	0.254		0.214	0.209	0.221		0.210	0.202	0.207
$\alpha =0.10$	0.121	0.135	0.147		0.111	0.104	0.113		0.105	0.099	0.103
$\alpha =0.05$	0.067	0.076	0.086		0.056	0.055	0.060		0.051	0.053	0.053
$\alpha =0.01$	0.016	0.024	0.030		0.011	0.013	0.013		0.012	0.011	0.012

Table 14. Table 14 : Frequencies of the true bivariate copula identified using the model selection algorithm from 100 simulation runs.

Political-economic dataset – 1-factor model
	Continuous			Ordinal
$n$	1-reflected Joe	Joe		reflected Joe	Joe	Gumbel
100	88	81		45	82	34
300	88	93		54	83	60
500	91	100		66	100	79

Equations

$$F_{j|X_{1}}(y|x):=\Pr(Y_{j}\leq y|X_{1}=x)={\partial C_{X_{1}j}\bigl{(}x,F_{j}(y)\bigr{)}\over\partial x}.$$

$$f_{\mathbf{Y}}(\mathbf{y})={\int }_0^1\prod _{j=1}^d{f}_{j|X_1}(y_j|x)dx,$$

$$f_{j|X_{1}}(y|x)=\left\{\begin{array}[]{ccc}C_{j|X_{1}}\bigl{(}F_{j}(y)|x\bigr{)}-C_{j|X_{1}}\bigl{(}F_{j}(y-1)|x\bigr{)}&\mbox{if}&Y_{j}\quad\mbox{is discrete};\\ c_{X_{1}j}\bigl{(}x,F_{j}(y)\bigr{)}f_{j}(y)&\mbox{if}&Y_{j}\quad\mbox{is continuous},\end{array}\right.$$

$$\Pr(X_{2}\leq x_{2},Y_{j}\leq y|X_{1}=x_{1})=C_{X_{2}j}\bigl{(}x_{2},F_{j|X_{1}}(y|x_{1})\bigr{)},$$

$$\mathrm{Pr}(Y_j\le y|X_1=x_1,X_2=x_2)=\frac{\partial }{\partial x_2}\mathrm{Pr}(X_2\le x_2,Y_j\le y|X_1=x_1)$$

$$={\partial\over\partial x_{2}}C_{X_{2}j}\bigl{(}x_{2},F_{j|X_{1}}(y|x_{1})\bigr{)}=C_{j|X_{2}}\bigl{(}F_{j|X_{1}}(y|x_{1})|x_{2}\bigr{)}.$$

$$f_{\mathbf{Y}}(\mathbf{y})=\int_{0}^{1}\int_{0}^{1}\prod_{j=1}^{d}f_{X_{2}j|X_{1}}\bigl{(}x_{2},y_{j}|x_{1}\bigr{)}\,dx_{1}dx_{2},$$

$$f_{X_{2}j|X_{1}}(x_{2},y|x_{1})=\left\{\begin{array}[]{ccc}C_{j|X_{2}}\bigl{(}F_{j|X_{1}}(y|x_{1})|x_{2}\bigr{)}-C_{j|X_{2}}\bigl{(}F_{j|X_{1}}(y-1|x_{1})|x_{2}\bigr{)}&\mbox{if}&Y_{j}\quad\mbox{is discrete};\\ c_{jX_{2};X_{1}}\bigl{(}F_{j|X_{1}}(y|x_{1}),x_{2}\bigr{)}c_{X_{1}j}\bigl{(}x_{1},F_{j}(y)\bigr{)}f_{j}(y)&\mbox{if}&Y_{j}\quad\mbox{is continuous}.\end{array}\right.$$

$${\rho }_N^{+}$$

$${\rho }_N^{-}$$

$$Y_j=y_j\text{mbox}if{\alpha }_{y_j-1,j}\le Z_j\le {\alpha }_{y_{j}j},y_j=1,\dots ,K_j,$$

$$\hat{\rho}_{N}=\mbox{argmax}_{\rho}\sum_{i=1}^{n}\log\Bigl{(}\Phi_{2}(\alpha_{y_{i1}},\alpha_{y_{i2}})-\Phi_{2}(\alpha_{y_{i1}-1},\alpha_{y_{i2}})-\Phi_{2}(\alpha_{y_{i1}},\alpha_{y_{i2}-1})+\Phi_{2}(\alpha_{y_{i1}-1},\alpha_{y_{i2}-1})\Bigr{)}$$

$$\hat{\rho}_{N}=\mbox{argmax}_{\rho}\sum_{i=1}^{n}\log\biggl{\{}\phi(z_{i1})\biggl{(}\Phi\Bigl{(}\frac{\alpha_{y_{i2}}-\rho z_{i1}}{(1-\rho^{2})^{1/2}}\Bigr{)}-\Phi\Bigl{(}\frac{\alpha_{y_{i2}-1}-\rho z_{i1}}{(1-\rho^{2})^{1/2}}\Bigr{)}\biggr{)}\biggr{\}}$$

$$F_j(y_{ij})=\frac{1}{n+1}\sum _{i=1}^{n}1(Y_{ij}\le y_{ij})=R_{ij}/(n+1),$$

$$f_j(y_j;{\mathbf{\gamma }}_j)=\Phi ({\alpha }_{y_{j}j})-\Phi ({\alpha }_{y_j-1,j}),$$

$$f_j(y_j;{\mathbf{\gamma }}_j)=\frac{\Gamma ({\xi }_j^{-1}+y_j)}{\Gamma ({\xi }_j^{-1})y_j!}\frac{{\mu }_j^y{\xi }_j^y}{(1+{\xi }_j^{-1}{)}^{{\xi }_j^{-1}+y_j}},y_j=0,1,2,\dots ,{\mu }_j>0,{\xi }_j>0,$$

$${\ell }_j({\mathbf{\gamma }}_j)=\sum _{i=1}^n\log f_j(y_{ij};{\mathbf{\gamma }}_j)$$

$${\ell }_{\mathbf{Y}}(\mathbf{\theta })=\sum _{i=1}^n\log f_{\mathbf{Y}}(y_{i1},\dots ,y_{id};\mathbf{\theta }).$$

$$f_{\mathbf{Y}}(\mathbf{y})={\int }_0^1\prod _{j=1}^d{f}_{j|X_1}(y_j|x)dx\approx \sum _{q=1}^{n_q}{w}_q\prod _{j=1}^d{f}_{j|X_1}(y_j|x_q),$$

$$f_{\mathbf{Y}}(\mathbf{y})$$

$$\displaystyle\Delta_{1f_{\mathbf{Y}}}=n^{-1}\Big{[}\sum_{i=1}^{n}\Bigl{\{}E_{f_{\mathbf{Y}}}\log f_{\mathbf{Y}}(\mathbf{y}_{i})-E_{f_{\mathbf{Y}}}\log f^{(1)}_{\mathbf{Y}}(\mathbf{y}_{i};\boldsymbol{\theta}_{1})\Bigr{\}}\Bigr{]},$$

$$\displaystyle\Delta_{2f_{\mathbf{Y}}}=n^{-1}\Bigl{[}\sum_{i=1}^{n}\Bigl{\{}E_{f_{\mathbf{Y}}}\log f_{\mathbf{Y}}(\mathbf{y}_{i})-E_{f_{\mathbf{Y}}}\log f^{(2)}_{\mathbf{Y}}(\mathbf{y}_{i};\boldsymbol{\theta}_{2})\Bigr{\}}\Bigr{]}.$$

$$\bar{D}=\sum _{i=1}^n{D}_i/n,$$

$$\sqrt{n}\bar{D}/s\sim N(0,1),$$

$$s=\sum _{j=1}^d(K_j-1)+\sum _{1\le j_1<j_2\le d}(K_{j_1}-1)(K_{j_2}-1).$$

$$M_{2}=M_{2}(\hat{\boldsymbol{\theta}})=n\bigl{(}\mathbf{p}_{2}-\boldsymbol{\pi}_{2}(\hat{\boldsymbol{\theta}})\bigr{)}^{\top}\mathbf{C}_{2}(\hat{\boldsymbol{\theta}})\bigl{(}\mathbf{p}_{2}-\boldsymbol{\pi}_{2}\bigl{(}\hat{\boldsymbol{\theta}})\bigr{)},$$

$$\mathbf{C}_{2}(\boldsymbol{\theta})=\boldsymbol{\Xi}_{2}^{-1}-\boldsymbol{\Xi}_{2}^{-1}\boldsymbol{\Delta}_{2}(\boldsymbol{\Delta}_{2}^{\top}\boldsymbol{\Xi}_{2}^{-1}\boldsymbol{\Delta}_{2})^{-1}\boldsymbol{\Delta}_{2}^{\top}\boldsymbol{\Xi}_{2}^{-1}=\boldsymbol{\Delta}_{2}^{(c)}\bigl{(}[\boldsymbol{\Delta}_{2}^{(c)}]^{\top}\boldsymbol{\Xi}_{2}\boldsymbol{\Delta}_{2}^{(c)}\bigr{)}^{-1}[\boldsymbol{\Delta}_{2}^{(c)}]^{\top},$$

$${\theta }_j={\beta }_{j1},{\delta }_j=\frac{{\beta }_{j2}}{(1-{\beta }_{j1}^2{)}^{1/2}},$$

Full text

Key Words: Canonical vines; Conditional independence; Goodness-of-fit; Latent variable models; Model selection; Tail dependence/asymmetry.

1 Introduction

It is very common in social science, e.g., in surveys, to deal with datasets that have mixed continuous and discrete responses. In the literature, two broad frameworks have been considered to model the dependence among such mixed continuous and discrete responses, namely the latent variable and copula framework.

There are two approaches for modelling multivariate mixed data with latent variables: the underlying variable approach that treats all variables as continuous by assuming the discrete responses as a manifestation of underlying continuous variables that usually follow the normal distribution (e.g., Muthén, 1984; Lee et al., 1992; Quinn, 2004), and the response function approach that postulates distributions on the observed variables conditional on the latent variables usually from the exponential family (e.g., Moustaki, 1996; Moustaki and Knott, 2000; Wedel and Kamakura, 2001; Huber et al., 2004; Moustaki and Victoria-Feser, 2006). The former method almost invariably assumes that the underlying variables (linked to the observed variables via a threshold process to yield ordinal data and an identity process to yield continuous data) follow a multivariate normal (MVN) distribution, while the latter assumes that the observed variables are conditionally independent usually given MVN distributed latent variables. They are equivalent when in the underlying and the response function approach the MVN distribution has a factor and an independence correlation structure, respectively (Takane and de Leeuw,, 1987).

The underlying variable approach calls the MVN distribution as a latent model for the discrete responses, and therefore maximum likelihood (ML) estimation requires multidimensional integrations (Nikoloulopoulos,, 2013, 2016); their dimension is equal with the number of observed discrete variables. This is why alternative estimation methods such as the three stage weighted least squares and composite likelihood have been proposed; see e.g., Katsikatsou et al., (2012). The response function approach, with the dependence coming from $p$ latent (unobservable) variables/factors where $p<<d$ (the number of observed variables), requires $p$- rather than $d$- dimensional integration. Hence, ML estimation is feasible, especially when the number of latent variables is small.

Nevertheless, both approaches are restricted to MVN assumption for the observed or the latent variables that is not valid if tail asymmetry or tail dependence exists in the mixed data which is a realistic scenario. Ma and Genton, (2010), Montanari and Viroli, (2010), and Irincheeva et al., 2012a stress that the MVN assumption might not be adequate, and acknowledge that the effect of misspecifying the distribution of the latent variables could lead to biased model estimates and poor fit. To this end, Irincheeva et al., 2012b proposed a more flexible response function approach by strategically multiplying the MVN density of the latent variables by a polynomial function to achieve departures from normality.

As we have discussed, the underlying variable approach exploits the use of the MVN assumption to model the joint distribution of mixed data. The univariate margins are transformed to normality and then the MVN distribution is fitted to the transformed data. This construction is apparently the MVN copula applied to mixed data (Shen and Weissfeld,, 2006; Hoff,, 2007; Song et al.,, 2009; He et al.,, 2012; Jiryaie et al.,, 2016), but previous papers (e.g., Quinn, 2004) do not refer to copulas as the approach can be explained without copulas.

Smith and Khaled, (2012), Stöber et al., (2015) and Zilko and Kurowicka, (2016) called vine copulas to model mixed data. Vine copulas have two major advantages over the MVN copula as emphasized in Panagiotelis et al., (2017). The first is that the computational complexity of computing the joint probability distribution function grows quadratically with $d$, whereas for the MVN copula the computational complexity grows exponentially with $d$. The second is that vine copulas are highly flexible through their specification from bivariate parametric copulas with different tail dependence or asymmetry properties. They have as special case the MVN copula, if all the bivariate parametric copulas are bivariate normal (BVN).

In this paper, we extend the factor copula models in Krupskii and Joe, (2013) and Nikoloulopoulos and Joe, (2015) to the case of mixed continuous and discrete responses. Factor copula are vine copula models that involve both observed and latent variables. Hence, they are highly flexible through their specification from bivariate parametric copulas with different tail dependence or asymmetry properties. The underlying variable approach where the MVN distribution has a $p$-factor correlation structure or its equivalent, i.e., the response function approach where the MVN distribution has an independence correlation structure is a special case of factor copula models when all the bivariate parametric copulas are BVN (hereafter standard factor model). Factor copula models are more interpretable and fit better than vine copula models, when dependence can be explained through latent variables. Furthermore, they are closed under margins, i.e., lower order marginals belong to the same parametric family of copulas and a different permutation of the observed variables has exactly the same distribution. This is not the case for vine copulas without latent variables, where a different permutation of the observed variables could lead to a different distribution.

We tackle issues that particularly interest the social data analyst as model selection and goodness-of-fit. Model selection in previous papers on factor copula models (Krupskii and Joe,, 2013; Nikoloulopoulos and Joe,, 2015) was mainly based on simple diagnostics. In addition to simple diagnostics based on semi-correlations, we propose an heuristic method that automatically selects the bivariate parametric copula families. As regard as to the issue of goodness-of-fit testing, we propose a technique that is based on the $M_{2}$ goodness-of-fit statistic (Maydeu-Olivares and Joe,, 2006) in multidimensional contingency tables to overcome the shortage of goodness-of-fit statistics for mixed continuous and discrete response data (e.g., Moustaki and Knott, 2000).

The remainder of the paper proceeds as follows. Section 2 introduces the factor copula models for mixed data and provides choices of parametric bivariate copulas with latent variables. Estimation techniques and computational details are provided in Section 3. Sections 4 and 5 propose methods for model selection and goodness-of-fit, respectively. Section 6 presents applications of our methodology to three mixed response data sets. Section 7 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. We conclude with some discussion in Section 8, followed by a brief section with software details.

2 The factor copula model for mixed responses

Although the factor copula models can be explained as truncated canonical vines rooted at the latent variables, we derive the models as conditional independence models, i.e., a response function approach with dependence coming from latent (unobservable) variables/factors. The $p$-factor model assumes that the mixed continuous and discrete responses $\mathbf{Y}=(Y_{1},\ldots,Y_{d})$ are conditionally independent given $p$ latent variables $X_{1},\ldots,X_{p}$. In line with Krupskii and Joe, (2013) and Nikoloulopoulos and Joe, (2015), we use a general copula construction, based on a set of bivariate copulas that link observed to latent variables, to specify the factor copula models for mixed continuous and discrete variables. The idea in the derivation of this $p$-factor model will be shown below for the 1-factor and 2-factor case. It can be extended to $p\geq 3$ factors or latent variables in a similar manner. The evaluation of a $p$-dimensional integral can be successively performed as we strategically assume that the factors or latent variables are independent.

For the 1-factor model, let $X_{1}$ be a latent variable, which we assume to be standard uniform (without loss of generality). From Sklar, (1959), there is a bivariate copula $C_{X_{1}j}$ such that \Pr(X_{1}\leq x,Y_{j}\leq y)=C_{X_{1}j}\bigl{(}x,F_{j}(y)\bigr{)} for $0\leq x\leq 1$ where $F_{j}$ is the cumulative distribution function (cdf) of $Y_{j}$. Then it follows that

[TABLE]

Letting $C_{j|X_{1}}(F_{j}(y)|x)=\partial C_{X_{1}j}(x,F_{j}(y))/\partial x$ for shorthand notation and $\mathbf{y}=(y_{1},\ldots,y_{d})$ be realizations of $\mathbf{Y}$, the density333We mean the density of $\mathbf{Y}$ w.r.t. the product measure on the respective supports of the marginal variables. For discrete margins with integer values this is the counting measure on the set of possible outcomes, for continuous margins we consider the Lebesgue measure in $\mathbb{R}$. of the observed data in the 1-factor model case is

[TABLE]

where

[TABLE]

is the density of $Y_{j}=y$ conditional on $X_{1}=x$; $c_{X_{1}j}$ is the bivariate copula density of $X_{1}$ and $Y_{j}$ and $f_{j}$ is the univariate density of $Y_{j}$.

For the 2-factor model, consider two latent variables $X_{1},X_{2}$ that are, without loss of generality, independent uniform $U(0,1)$ random variables. Let $C_{X_{1}j}$ be defined as in the 1-factor model, and let $C_{X_{2}j}$ be a bivariate copula such that

[TABLE]

where $F_{j|X_{1}}$ is given in (1). Then for $0\leq x_{1},x_{2}\leq 1$,

[TABLE]

[TABLE]

The density of the observed data in the 2-factor model case is

[TABLE]

where

[TABLE]

Note that the copula $C_{X_{1}j}$ links the $j$th response to the first latent variable $X_{1}$, and the copula $C_{X_{2}j}$ links the $j$th response to the second latent variable $X_{2}$ conditional on $X_{1}$. In our general statistical model there are no constraints in the choices of the parametric marginal $F_{j}$ or copula $\{C_{X_{1}j},C_{X_{2}j}\}$ distributions.

2.1 Choices of bivariate copulas with latent variables

We provide choices of parametric bivariate copulas that can be used to link the latent to the observed variables. We will consider copula families that have different tail dependence (Joe,, 1993) or tail order (Hua and Joe,, 2011).

A bivariate copula $C$ is reflection symmetric if its density satisfies $c(u_{1},u_{2})=c(1-u_{1},1-u_{2})$ for all $0\leq u_{1},u_{2}\leq 1$. 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})$ as $u\to 0$ and lower tail dependence means that $c(u,u)=O(u^{-1})$ as $u\to 0$. If $(U_{1},U_{2})\sim C$ for a bivariate copula $C$, then $(1-U_{1},1-U_{2})\sim\widehat{C}$, where $\widehat{C}(u_{1},u_{2})=u_{1}+u_{2}-1+C(1-u_{1},1-u_{2})$ 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. Under some regularity conditions (e.g., existing finite density in the interior of the unit square, ultimately monotone in the tail), if there exists $\kappa_{L}(C)>0$ and some $L(u)$ that is slowly varying at $0^{+}$ (i.e., $\frac{L(ut)}{L(u)}\sim 1,$ as $u\to 0^{+}$ for all $t>0$), then $\kappa_{L}(C)$ is the lower tail order of $C$. The upper tail order $\kappa_{U}(C)$ can be defined by the reflection of $(U_{1},U_{2})$, i.e., ${\overline{C}}(1-u,1-u)\sim u^{\kappa_{U}(C)}L^{*}(u)$ as $u\to 0^{+}$, where ${\overline{C}}$ is the survival function of the copula and $L^{*}(u)$ is a slowly varying function. With $\kappa=\kappa_{L}$ or $\kappa_{U}$, a bivariate copula has intermediate tail dependence if $\kappa\in(1,2)$, tail dependence if $\kappa=1$, and tail quadrant independence if $\kappa=2$ with $L(u)$ being asymptomatically a constant.

After briefly providing definitions of tail dependence and tail order we provide below a list of bivariate parametric copulas with varying tail behaviour:

• •

Reflection symmetric copulas with intermediate tail dependence such as the BVN copula with $\kappa_{L}=\kappa_{U}=2/(1+\theta)$, where $\theta$ is the copula (correlation) parameter.

• •

Reflection symmetric copulas with tail quadrant independence ($\kappa_{L}=\kappa_{U}=2$), such as the Frank copula.

• •

Reflection asymmetric copulas with upper tail dependence only such as

• –

the Gumbel copula with $\kappa_{L}=2^{1/\theta}$ and $\kappa_{U}=1$, where $\theta$ is the copula parameter;

• –

the Joe copula with $\kappa_{L}=2$ and $\kappa_{U}=1$.

• •

Reflection symmetric copulas with tail dependence, such as the $t_{\nu}$ copula with $\nu$ the degrees of freedom and $\kappa_{L}=\kappa_{U}=1$.

• •

Reflection asymmetric copulas with upper and lower tail dependence that can range independently from 0 to 1, such as the BB1 and BB7 copulas with $\kappa_{L}=1$ and $\kappa_{U}=1$.

• •

Reflection asymmetric copulas with tail quadrant independence, such as the the BB8 and BB10 copulas.

The BVN, Frank, and $t_{\nu}$ are comprehensive copulas, i.e., they interpolate between countermonotonicity (perfect negative dependence) to comonotonicity (perfect positive dependence). The other aforementioned parametric families of copulas, namely Gumbel, Joe, BB1, BB7, BB8 and BB10 interpolate between independence and perfect positive dependence. Nevertheless, negative dependence can be obtained from these copulas by considering reflection of one of the uniform random variables on (0, 1). If $(U_{1},U_{2})\sim C$ for a bivariate copula $C$ with positive dependence, then

• •

$(1-U_{1},U_{2})\sim\widehat{C}^{(1)}$, where $\widehat{C}^{(1)}(u_{1},u_{2})=u_{2}-C(1-u_{1},u_{2})$ is the 1-reflected copula of $C$ with negative lower-upper tail dependence;

• •

$(U_{1},1-U_{2})\sim\widehat{C}^{(2)}$, where $\widehat{C}^{(2)}(u_{1},u_{2})=u_{1}-C(u_{1},1-u_{2})$ 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})$ as $u\to 0^{+}$ and negative lower-upper tail dependence means that $c(u,1-u)=O(u^{-1})$ as $u\to 0^{+}$ (Joe,, 2011).

In Figure 1, to depict the concepts of refection symmetric or asymmetric tail dependence or quadrant tail independence, we show contour plots of the corresponding copula densities with standard normal margins and dependence parameters corresponding to Kendall’s $\tau$ value of 0.5. Sharper corners (relative to ellipse) indicate tail dependence.

2.2 **Semi-correlations to detect tail dependence or tail asymmetry **

Choices of copulas with upper or lower tail dependence are better if the observed variables have more probability in joint upper or lower tail than would be expected with the standard factor model. This can be shown with summaries of correlations in the upper joint tail and lower joint tail.

For continuous variables, although copula theory uses transforms to standard uniform margins $U_{j}=F_{j}(Y_{j})$, we convert to normal scores $Z_{j}=\Phi^{-1}(U_{j})$ to check deviations from the elliptical shape that would be expected with the BVN copula (Nikoloulopoulos et al.,, 2012). The correlations of normal scores in the upper and lower tail, hereafter semi-correlations, are defined as (Joe,, 2014, page 71):

[TABLE]

Note in passing that for the BVN copula $\rho_{N}^{+}=\rho_{N}^{-}$ and has a closed form; see (Joe,, 2014, page 71).

From the above expressions, it is apparent that the normal scores semi-correlations depend only on the copula $C$ of $(U_{j_{1}},U_{j_{2}})$. Table 1 has semi-correlations for all the aforementioned bivariate parametric copulas with $\tau=\{0.3,0.5,0.7\}$. From the table we can see that $\rho_{N}^{+}=\rho_{N}^{-}$ for any reflection symmetric copula, while they are different for any reflection asymmetric one. If there is stronger upper (lower) tail dependence than with the BVN, then the upper (lower) semi-correlation is larger.

The population versions $\rho_{N}^{+},\rho_{N}^{-}$ also apply when the variables $Y_{j}$ are ordinal. Under the univariate probit model (Agresti,, 2010, Section 3.3.2) $Z_{j}$ are standard normal underlying latent variables, such that

[TABLE]

where $K_{j}$ is the number of categories of $Y_{j}$ and $a_{1j},\ldots,a_{K_{j}-1,j}$ are the univariate cutpoints (without loss of generality, we assume $\alpha_{0j}=-\infty$ and $\alpha_{K_{j}j}=\infty$). Note in passing that for binary variables ($K_{j}=2$) the calculation of the semi-correlations is meaningless as the binary variables have no tail asymmetries.

The sample versions of $\rho_{N}^{+},\rho_{N}^{-}$ are sample linear (when both variables are continuous), polychoric (when both variables are ordinal), and polyserial (when one variable is continuous and the other is ordinal) correlations in the joint lower and upper quadrants of the two variables. The sample polychoric and polyserial correlation is defined as

[TABLE]

and

[TABLE]

with z_{ij}=\Phi\Bigl{(}\frac{1}{n+1}\sum_{i=1}^{n}\mathbf{1}(Y_{ij}\leq y_{ij})\Bigr{)}, respectively.

3 Estimation

We use a two-stage copula modelling approach toward the estimation of a multivariate model that borrows the strengths of the semi-parametric and IFM approach in Genest et al., (1995) and Joe, (2005), respectively. Suppose that the data are $y_{ij},\,j=1,\ldots,d,\,i=1,\ldots,n$, where $i$ is an index for individuals or clusters and $j$ is an index for the within-cluster measurements. For $i=1,\ldots,n$ we start from a $d$-variate sample $y_{i1},\ldots,y_{id}$ from which $d$ estimators $F_{1}(y_{i1}),\ldots,F_{d}(y_{id})$ can be obtained. We use these to transform the $y_{i1},\ldots,y_{id}$ sample to a uniform sample $u_{i1}=F_{1}(y_{i1}),\ldots,u_{id}=F_{d}(y_{id})$ on $[0,1]^{d}$ and then fit the factor copula model at the second step. For continuous and discrete data $y_{ij}$, we use non-parametric and parametric univariate distributions, respectively, to transform the data $y_{ij}$ to copula data $u_{ij}=F_{j}(y_{ij})$, i.e., data on the uniform scale. Hence our proposed approach, in line with the approaches in Genest et al., (1995) and Joe, (2005), 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.

3.1 Univariate modelling

For continuous random variables, we estimate each marginal distribution non-parametrically by the empirical distribution function of $Y_{j}$, viz.

[TABLE]

where $R_{ij}$ denotes the rank of $Y_{ij}$ as in the semi-parametric estimation of Genest et al., (1995) and Shih and Louis, (1995). Hence we allow the distribution of the continuous margins to be quite free and not restricted by parametric families.

Nevertheless, rank-based methods cannot be used for discrete variables with copulas (Genest and Nešlehová,, 2007). Hence, for both ordinal and count variables we have chosen realistic parametric models:

• •

For an ordinal response variable $Y_{j}$ we use the univariate probit model in (4). The ordinal response $Y_{j}$ is assumed to have density

[TABLE]

where $\boldsymbol{\gamma}_{j}=(a_{1j},\ldots,a_{K_{j}-1,j})$ is the vector of the univariate cutpoints.

• •

For a count response variable $Y_{j}$ we use the negative binomial distribution (Lawless,, 1987). It allows for over-dispersion and its probability mass function is

[TABLE]

where $\boldsymbol{\gamma}_{j}=\{\mu_{j},\xi_{j}\}$ is the vector with the mean and dispersion parameters. In the limit $\xi\to 0$ the negative binomial reduces to Poisson, which belongs to the exponential family of distributions and it is the only distribution for count data that existing latent variable models for mixed data can accommodate.

To this end, for a discrete random variable $Y_{j}$, we approach estimation by maximizing the univariate log-likelihoods

[TABLE]

over the vector of the univariate parameters $\boldsymbol{\gamma}_{j}$. That is equivalent with the first step of the IFM method in Joe, (1997, 2005).

In line with the IFM method, if one uses a misspecified univariate model for the discrete responses at the first step, then the estimation of the copula parameters at the second step deteriorates as demostrated in Kim et al., (2007). Nevertheless, there is no “correct specification” of the margins or copula for data analysis. If one does a proper analysis of the univariate margins for goodness-of-fit, then the proposed two-stage (or the IFM) method should be fine. Kim et al., (2007) have “true univariate distributions for simulations” and “specified univariate distributions for estimation” that were very far apart and unrealistic, because the difference of the two is easily detected without too much data.

3.2 Copula modelling

After estimating the univariate marginal distributions we proceed to estimation of the dependence parameters. For the 1-factor and 2-factor models, we let $C_{X_{1}j}$ and $C_{X_{2}j}$ be parametric bivariate copulas, say with dependence parameters $\theta_{j}$ and $\delta_{j}$, respectively. Let also $\boldsymbol{\theta}=\{\boldsymbol{\gamma}_{j},\theta_{j}:j=1,\ldots,d\}$ and $\boldsymbol{\theta}=\{\boldsymbol{\gamma}_{j},\theta_{j},\delta_{j}:j=1,\ldots,d\}$ to denote the set of all parameters for the 1- and 2-factor model, respectively. Estimation can be achieved by maximizing the joint log-likelihood

[TABLE]

over the copula parameters $\theta_{j}$ or $\delta_{j},\,j=1,\ldots,d$ with the univariate parameters/distributions fixed as estimated at the first step of the proposed two-step estimation approach. The estimated parameters can be obtained by using a quasi-Newton Nash, (1990) method applied to the logarithm of the joint likelihood. This numerical method requires only the objective function, i.e., the logarithm of the joint likelihood, while the gradients are computed numerically and the Hessian matrix of the second order derivatives is updated in each iteration. The standard errors (SE) of the estimates can be obtained via the gradients and the Hessian computed numerically during the maximization process. These SEs are adequate to assess the flatness of the log-likelihood. Proper SEs that account for the estimation of univariate parameters can be obtained by maximizing the joint likelihood in (5) at one step over $\boldsymbol{\theta}$.

For factor copula models numerical evaluation of the joint density $f_{\mathbf{Y}}(\mathbf{y;\boldsymbol{\theta}})$ can be easily done using Gauss-Legendre quadrature (Stroud and Secrest,, 1966). To compute one-dimensional integrals for the 1-factor model, we use the following approximation:

[TABLE]

where $\{x_{q}:q=1,\ldots,n_{q}\}$ are the quadrature points and $\{w_{q}:q=1,\ldots,n_{q}\}$ are the quadrature weights. To compute two-dimensional integrals for the 2-factor model, the approximation uses Gauss-Legendre quadrature points in a double sum:

[TABLE]

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 (Nash,, 1990). Our comparisons show that $n_{q}=25$ is adequate with good precision.

4 Model selection

In this section we propose an heuristic method that automatically selects the bivariate parametric copula families that link the observed to the latent variables. This is very useful when the direction to the tail asymmetry based on semi-correlations is not consistent or clear. For multivariate mixed data, it is infeasible to estimate all possible combinations of bivariate parametric copula families, and compare them on the basis of information criteria. We develop an algorithm that can quickly select a factor copula model that accurately captures the (tail) dependence features in the data on hand. The linking copulas at each factor are selected with a sequential algorithm under the initial assumption that linking copulas are Frank, and then sequentially copulas with non-tail quadrant independence are assigned to any of pairs where necessary to account for tail asymmetry (discrete data) or tail dependence (continuous data).

For the 1-factor model, the proposed model selection algorithm is summarized in the following steps:

For $j=1,\ldots,d$ estimate the marginal distributions $F_{j}(y)$. 2. 2.

Fit the 1-factor copula model with Frank copulas to link each of the $d$ observed variables with the latent variable, i.e., maximise the log-likelihood function of the factor copula model in (5) over the vector of copula parameters $(\theta_{1},\ldots,\theta_{d})$. 3. 3.

If the $j$th linking copula has $\hat{\theta}_{j}>0$, then select a set of copula candidates with ability to interpolate between independence and comonotonicity, otherwise select a set of copula candidates with ability to interpolate between countermonotonicity and independence. 4. 4.

For $j=1,\ldots,d$,

1. (a)

Fit all the possible 1-factor copula models iterating over all the copula candidates for the $j$th variable. 2. (b)

Select the copula family that corresponds to the lowest information criterion, say the Akaike, that is $\text{AIC}=-2\times\ell+2\times\#\text{copula parameters}$. 3. (c)

Fix the selected linking copula family for the $j$th variable. 4. (d)

Iterate through step (a) – (c) to select the copulas that link all the observed variables to the 1st factor.

For more than one factor we can select the appropriate linking copulas accordingly. We first select copula families in the first factor, and then we proceed to the next factor and apply exactly the same algorithm.

5 Techniques for parametric model comparison and goodness-of-fit

Factor copula models with different bivariate linking copulas can be compared via the log-likelihood or AIC at the maximum likelihood estimate. In addition, we will use the Vuong’s test (Vuong,, 1989) to show if a factor copula model provides better fit than the standard factor model with a latent additive structure, that is a factor copula model with BVN bivariate linking copulas (Krupskii and Joe,, 2013; Nikoloulopoulos and Joe,, 2015). 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. This test has been used extensively in the copula literature to compare vine copula models (e.g., Brechmann et al., 2012; Joe, 2014; Nikoloulopoulos, 2017). We provide specific details in Section 5.1.

Furthermore, to assess the overall goodness-of-fit of the factor copula models for mixed data, we will use appropriately the limited information $M_{2}$ statistic (Maydeu-Olivares and Joe,, 2006). The $M_{2}$ statistic has been developed for goodness-of-fit testing in multidimensional contingency tables. Nikoloulopoulos and Joe, (2015) has used the $M_{2}$ statistic to assess the goodness-of-fit of factor copula models for ordinal data. We build on the aforementioned papers and propose a methodology to assess the overall goodness-of-fit of factor copula models for mixed continuous and discrete responses. We provide the specifics for the $M_{2}$ statistic in Section 5.2.

5.1 Vuong’s test for parametric model comparison

In this subsection, we summarize the Vuong’s test for comparing parametric models (Vuong,, 1989). Assume that we have Models 1 and 2 with parametric densities $f^{(1)}_{\mathbf{Y}}$ and $f^{(2)}_{\mathbf{Y}}$, respectively. We can compare

[TABLE]

where $\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{2}$ are the parameters in Models 1 and 2, respectively, that lead to the closest Kullback-Leibler divergence to the true $f_{\mathbf{Y}}$; equivalently, they are the limits in probability of the MLEs based on Models 1 and 2, respectively.

Model 1 is closer to the true $f_{\mathbf{Y}}$, i.e., is the better fitting model if $\Delta=\Delta_{1f_{\mathbf{Y}}}-\Delta_{2f_{\mathbf{Y}}}<0$, and Model 2 is the better fitting model if $\Delta>0$. The sample version of $\Delta$ with MLEs $\widehat{\boldsymbol{\theta}}_{1},\widehat{\boldsymbol{\theta}}_{2}$ is

[TABLE]

where $D_{i}=\log\left[\frac{f^{(2)}_{\mathbf{Y}}(\mathbf{y}_{i};\widehat{\boldsymbol{\theta}}_{2})}{f^{(1)}_{\mathbf{Y}}(\mathbf{y}_{i};\widehat{\boldsymbol{\theta}}_{1})}\right]$. Vuong, (1989) has shown that asymptotically that

[TABLE]

where $s^{2}=\frac{1}{n-1}\sum_{i=1}^{n}(D_{i}-\bar{D})^{2}$. Hence, its $95\%$ confidence interval (CI) is $\bar{D}\pm 1.96\times\frac{1}{\sqrt{n}}\sigma$.

5.2 M2 goodness-of-fit statistic

Since the $M_{2}$ statistic has been developed for multivariate ordinal data (Maydeu-Olivares and Joe,, 2006), we propose to first transform the continuous and count variables to ordinal and then calculate the $M_{2}$ statistic at the maximum likelihood estimate before transformation.

Continuous variables can be transformed to ordinal with categories that are meaningful both practically and scientifically. If this is not the case, we propose an unsupervised strategy of transforming a continuous into an ordinal variable:

Set the number of ordinal categories $K_{j}$. 2. 2.

Transform $Y_{j}$ to a standard uniform random variable $U_{j}$ using its empirical distribution function. 3. 3.

Set the ordinal cutpoints on the uniform scale by generating a regular sequence from 1 to $K_{j}-1$ and then dividing over $K_{j}$. 4. 4.

Divide the range of $U_{j}$ into intervals with breaks the ordinal cutpoints. 5. 5.

Transform $U_{j}$ to an ordinal variable $Y_{j}$ according to which interval its values fall.

Count variables that contain very high counts or very low counts, can be treated as ordinal where the first or the last category contains all the low or high counts, respectively, and their other values remain as they are. We further propose an unsupervised strategy of categorising a count into an ordinal variable:

Set the number of ordinal categories $K_{j}$. 2. 2.

Divide the range of $Y_{j}$ into intervals with breaks a regular sequence of length $K_{j}+1$ from $\min(Y_{j})$ to $\max(Y_{j})$. 3. 3.

Transform $Y_{j}$ to an ordinal variable according to which interval its values fall.

After applying the transformations as above for each continuous or count variable, we have $d$ ordinal variables $Y_{1},\ldots,Y_{d}$ (both the original and the transformed ones) where the $j$th $(1\leq j\leq d)$ variable consists of $K_{j}\geq 2$ categories labelled as $0,1,\ldots,K_{j}-1$. Consider the set of univariate and bivariate residuals that do not include category 0. This is a residual vector of dimension

[TABLE]

For a factor copula model with parameter vector $\boldsymbol{\theta}$ of dimension $q$, let \boldsymbol{\pi}_{2}(\boldsymbol{\theta})=\bigl{(}\dot{\boldsymbol{\pi}}_{1}(\boldsymbol{\theta})^{\top},\dot{\boldsymbol{\pi}}_{2}(\boldsymbol{\theta})^{\top}\bigr{)}^{\top} be the column vector of the model-based marginal probabilities with $\dot{\boldsymbol{\pi}}_{1}(\boldsymbol{\theta})$ the vector of univariate marginal probabilities, and $\dot{\boldsymbol{\pi}}_{2}(\boldsymbol{\theta})$ the vector of bivariate marginal probabilities. Also, let $\mathbf{p}_{2}=(\dot{\mathbf{p}}_{1}^{\top},\dot{\mathbf{p}}_{2}^{\top})^{\top}$ be the vector of the observed sample proportions, with $\dot{\mathbf{p}}_{1}$ the vector of univariate marginal proportions, and $\dot{\mathbf{p}}_{2}$ the vector of the bivariate marginal proportions.

With a sample size $n$, the limited information statistic $M_{2}$ is given by

[TABLE]

with

[TABLE]

where $\boldsymbol{\Delta}_{2}=\partial\boldsymbol{\pi}_{2}(\boldsymbol{\theta})/\partial\boldsymbol{\theta}^{\top}$ is an $s\times q$ matrix with the derivatives of all the univariate and bivariate marginal probabilities with respect to the model parameters, $\boldsymbol{\Delta}_{2}^{(c)}$ is an $s\times(s-q)$ orthogonal complement to $\boldsymbol{\Delta}_{2}$, such that $[\boldsymbol{\Delta}_{2}^{(c)}]^{\top}\boldsymbol{\Delta}_{2}=\mathbf{0}$, and $\boldsymbol{\Xi}_{2}=\mbox{diag}(\boldsymbol{\pi}_{2}(\boldsymbol{\theta}))-\boldsymbol{\pi}_{2}(\boldsymbol{\theta})\boldsymbol{\pi}_{2}(\boldsymbol{\theta})^{\top}$ is the $s\times s$ covariance matrix of all the univariate and bivariate marginal sample proportions, excluding category 0. Due to equality in (7), $\mathbf{C}_{2}$ is invariant to the choice of orthogonal complement. The limited information statistic $M_{2}$ has a null asymptotic distribution that is $\chi^{2}$ with $s-q$ degrees of freedom when the estimate $\hat{\boldsymbol{\theta}}$ is $\sqrt{n}$-consistent. For details on the computation of $\boldsymbol{\Xi}_{2}$ and $\boldsymbol{\Delta}_{2}$ for factor copula models we refer the interested reader to Nikoloulopoulos and Joe, (2015).

6 Applications

In this section we illustrate the proposed methodology by re-analysing three mixed response datasets.

Initially, we use the diagnostic method in Joe, (2014, pages 245-246) to show that each dataset (or more precisely the correlation matrix of the observed variables for each dataset) has a factor structure based on linear factor analysis. The correlation matrix $\mathbf{R}_{\mathrm{observed}}$ has been obtained based on the sample correlations from the bivariate pairs of the observed variables. These are the linear (when both variables are continuous), polychoric (when both variables are discrete), and polyserial (when one variable is continuous and the other is discrete) sample correlations among the observed variables. The resulting $\mathbf{R}_{\mathrm{observed}}$ is generally positive definite if the sample size is not small enough; if not one has to convert it to positive definite. We calculate various measures of discrepancy between $\mathbf{R}_{\mathrm{observed}}$ and $\mathbf{R}_{\mathrm{model}}$ (the resulting correlation matrix of linear factor analysis), such as the maximum absolute correlation difference $D_{1}=\max|\mathbf{R}_{\mathrm{model}}-\mathbf{R}_{\mathrm{observed}}|$, the average absolute correlation difference $D_{2}=\mathrm{avg}|\mathbf{R}_{\mathrm{model}}-\mathbf{R}_{\mathrm{observed}}|$, and the correlation matrix discrepancy measure D_{3}=\log\bigl{(}\det(\mathbf{R}_{\mathrm{model}})\bigr{)}-\log\bigl{(}\det(\mathbf{R}_{\mathrm{observed}})\bigr{)}+\mathrm{tr}(\mathbf{R}^{-1}_{\mathrm{model}}\mathbf{R}_{\mathrm{observed}})-d.

After confirming that a factor model with a parsimonious correlation structure is reasonable, we calculate the semi-correlations for each pair of observed variables to check if there is tail asymmetry. This will be a useful information for choosing potential parametric bivariate copulas other than the BVN copulas that lead to the standard factor model. Note that when the variables are negatively associated we calculate the sample semi-correlations in the lower-upper and upper-lower quadrant.

After motivating why more flexible dependencies are needed in cases of mixed data and how those dependencies in the data can be captured by suitable bivariate copulas, we proceed with factor copula models and construct a plausible factor copula model, to capture any type of reflection asymmetric dependence, by using the proposed algorithm in Section 4. For a baseline comparison, we first fit the factor copula models with the comprehensive bivariate parametric copula families that allow for reflection symmetric dependence; these are the BVN, Frank, and $t_{\nu}$ copulas. For $t_{\nu}$ copulas, we summarize the choice of integer $\nu$ with the largest log-likelihood. For the standard 2-factor model, to obtain a unique solution we must impose sufficient constraints. One parameter for the second factor can be set to zero and the likelihood can be maximized with respect to other $2d-1$ parameters. We report the varimax transform (Kaiser,, 1958) of the loadings (a reparametrization of $2d$ parameters), converted to factor copula parameters via the relations

[TABLE]

where $\beta_{j1}$ and $\beta_{j2}$ are the loadings at the first and second factor, respectively (Krupskii and Joe,, 2013; Nikoloulopoulos and Joe,, 2015).

If the number of parameters is not the same between the models, we use the AIC as a rough diagnostic measure for goodness-of-fit between the models, otherwise we use the likelihood at the maximum likelihood estimates. We further compute the Vuong’s tests with Model 1 being the factor copula model with BVN copulas, that is the standard factor model, to reveal if any other factor copula model provides better fit than the standard factor model. To make it easier to compare strengths of dependence, we convert the estimated parameters to Kendall’s $\tau$’s in $(-1,1)$ via the relations in Joe, (2014, Chapter 4); SEs are also converted via the delta method. For the model that provides the best fit, we provide the estimates and SEs that are obtained by maximizing the joint likelihood in (5) at one step over $\boldsymbol{\theta}$. Although, the two-stage estimation approach in Section 3 is a convenient way to quickly compare candidate factor copula models, the full likelihood is applied for the best fitting factor copula model. The overall fit of the factor copula models is evaluated using the $M_{2}$ statistic. Note that the $M_{2}$ statistic in the case with $2d-1$ copulas (one set to independence for the second factor) is computed with $\boldsymbol{\Delta}_{2}$ having one less column.

6.1 Political-economic dataset

Quinn, (2004) considered measuring the (latent) political-economic risk of 62 countries, for the year 1987. The political-economic risk is defined as the country’s risk in manipulating economic rules for its own and constituents’ advantages (see e.g., North and Weingast, 1989). Quinn, (2004) used 5 mixed variables, namely the continuous variable black-market premium in each country (used as a proxy for illegal economic activity), the continuous variable productivity as measured by real gross domestic product per worker in 1985 international prices, the binary variable independence of the national judiciary (1 if the judiciary is judged to be independent and 0 otherwise), and the ordinal variables measuring the lack of expropriation risk and lack of corruption. The dataset and its complete description can be found in Quinn, (2004) or in the R package MCMCpack (Martin et al.,, 2011). Note that since the continuous variable black-market premium is negatively associated with the remaining variables (from the context), we re-orient it leading to positive dependence among all the observed variables.

Table 2 shows that the sample correlation matrix of the mixed responses has an 1-factor structure based on linear factor analysis (large $D_{3}$ is due to the small sample size as demonstrated using simulated data in Section 7). The sample semi-correlations in Table 2 show that there is more probability in the upper tail or lower tail compared with a discretized MVN, suggesting that a factor model with bivariate parametric copulas with upper or lower tail dependence might provide a better fit. Table 3 gives the estimated parameters, their standard errors (SE) in Kendall’s $\tau$ scale, joint log-likelihoods, the $95\%$ CIs of Vuong’s tests, and the $M_{2}$ statistics for the 1-factor copula models. Table 3 also indicates the parametric copula family chosen for each pair using the proposed heuristic algorithm. Copulas with asymmetric dependence are selected for all the copulas that link the latent variable to each of the observed variables. Hence, it is revealed that there are features in the data such as tail dependence and asymmetry which cannot be captured by copulas with reflection symmetric dependence such as BVN, Frank and $t_{\nu}$ copulas.

In all the fitted models the estimated Kendall’s $\tau$’s are similar. Kendall’s $\tau$ only accounts for the dependence dominated by the middle of the data, and it is expected to be similar amongst different families of copulas. However, the tail dependence and tail order vary, as explained in Section 2.1, and they are properties to consider when choosing amongst different families of copulas (Nikoloulopoulos and Karlis,, 2008).

The table shows that the selected model using the proposed algorithm provides the best fit and there is a substantial improvement over the standard factor model as indicated by the Vuong’s and $M_{2}$ statistics. To compute the $M_{2}$ statistics we transformed the continuous variables to ordinal with 5 categories using the unsupervised strategy in Section 5.2; similar inference was drawn, when we transformed them to ordinal with 3, 4, or 6 categories. The factor copula parameter of $0.51$ on negative black market premium indicates a negative association between the illegal economic activity and the latent variable. All the other estimated factor copula parameters indicate a positive association between each of the other observed variables (independent judiciary, productivity, lack of expropriation, and lack of corruption) with the latent variable. Hence, we can interpret the latent variable to be the political economical certainty.

6.2 General social survey

Hoff, (2007) analysed seven demographic variables of 464 male respondents to the 1994 General Social Survey. Of these seven, two were continuous (income and age of the respondents), three were ordinal with $5$ categories (highest degree of the survey respondent, income and highest degree of respondent’s parents), and two were count variables (number of children of the survey respondent and respondent’s parents). The data are available in Hoff, (2007, Supplemental materials).

Table 4 shows that the sample correlation matrix of the mixed responses has a 2- or even a 3-factor structure based on linear factor analysis. The direction to the tail asymmetry based on sample semi-correlations in Table 4 is not consistent, and this emerges the usefulness of the proposed model selection technique. Table 5 gives the estimated parameters, their standard errors (SE) in Kendall’s $\tau$ scale, the joint log-likelihoods, the $95\%$ CIs of Vuong’s tests, and the $M_{2}$ statistics for the 1-factor and 2-factor copula models. The best fit for the 1-factor model is based on the bivariate copulas selected by the proposed algorithm, where there is improvement over the factor copula model with BVN copulas according to the Vuong’s statistic. However, assessing the overall goodness-of-fit via the $M_{2}$ statistic, it is revealed that one latent variable is not adequate to explain the dependencies among the mixed responses. To apply the $M_{2}$ statistic, age and income were transformed to ordinal with $4$ (18–24, 25–44, 45–64, and $65+$) and $5$ (0–10, 11–19, 20–29, 30–40, and $41+$) categories, respectively, and number of children of the survey respondent and respondent’s parents were treated as ordinal where the 4th (more than 3 children) and 8th (more than 7 children) category, respectively, contained all the high counts.

The 2-factor copula models with BVN, $t_{\nu}$, and Frank copulas provide some improvement over the 1-factor copula models but according to the $M_{2}$ statistic they still have a poor fit. Note that the factor copula model with $t_{9}$ copulas was not identifiable (large SEs) in line with Nikoloulopoulos and Joe, (2015), hence one parameter for the second factor was set to zero and the likelihood was maximized with respect to the remaining parameters. We report the varimax transform (Kaiser,, 1958) of the loadings, converted to factor copula parameters via the relations in (8).

The selected 2-factor copula model using the algorithm in Section 4 shows improvement over the standard factor model according to the Vuong’s statistic and better fit according to the $M_{2}$ statistic; it changes a $p$-value $<0.001$ to one $>0.10$. For the 2-factor model based on the proposed algorithm for model selection, note that, without the need for a varimax rotation, the unique loading parameters ($\hat{\tau}$’s converted to normal copula parameters $\hat{\theta}_{j}$’s and $\hat{\delta}_{j}$’s and then to loadings using the relations in (8)) show that one factor is loaded only on the demographic variables of the respondent’s parents.

6.3 Swiss consumption survey

Irincheeva et al., 2012b considered measuring the latent variable ‘financial wealth of the household’ in its different realizations by analysing seven household variables of $n=9960$ respondents to the Swiss consumption survey. Out of these seven, three were continuous (food, clothing and leisure expenses), three were binary (dishwasher, car, and motorcycle), and one was count variable (the number of bicycles in possession of the household).

Irincheeva et al., 2012b , with simple descriptive statistics such as scatter plots of the original data, have shown that these mixed responses have reflection asymmetric dependence, and fitted their latent variable approach with one and two latent variables. In Figure 2 we depict the bivariate normal scores plots for the continuous data along with their correlations and semi-correlations. With a bivariate normal scores plot one can check for deviations from the elliptical shape that would be expected with the BVN copula, and hence assess if tail asymmetry and tail dependence exists on the data. For all the pairs the upper semi-correlation is larger, and interestingly, contrasting the bivariate normal scores plots in Figure 2 with the contour plots in Figure 1, it is apparent that for the continuous variables the linking copulas might be the BB10 copulas.

Table 6 shows that the sample correlation matrix of the mixed responses has a 2-factor structure based on linear factor analysis. The sample semi-correlations in Table 6 show that there is more probability in the upper tail and lower tail among the continuous variables and between each of the continuous variables with the count variable, respectively, suggesting that a factor model with bivariate parametric copulas with asymmetric tail dependence might provide a better fit. Table 7 gives the estimated parameters, their standard errors (SE) in Kendall’s tau scale, the joint log-likelihoods, the $95\%$ CIs of Vuong’s tests, and the $M_{2}$ statistics for the 1-factor and 2-factor copula models. The best fitted 1- and 2-factor models result when we use BB10 copulas with asymmetric quadrant tail independence to link the latent variable to each of the continuous observed variables and copulas with lower tail dependence to link the latent variables to the discrete observed variables. Once again the one-factor copula model is not adequate to explain the dependence amongst the mixed responses based on the $M_{2}$ statistic (Table 7, 1-factor). To apply the $M_{2}$ statistic, we transformed the continuous to ordinal variables with 3 categories using the unsupervised strategy in section 5 and the count variable bicycle was treated as ordinal where the 6th category contained all the high counts ($5$ bicycles or more).

While it is revealed that the selected 2-factor copula model is the best model (lowest AIC) and there is substantial improvement over the standard 2-factor model, it is not apparent from the $M_{2}$ statistic that the response patterns are satisfactorily explained by even 2 latent variables. This is not surprising since one should expect discrepancies between the postulated parametric model and the population probabilities, when the sample size is sufficiently large (Maydeu-Olivares and Joe,, 2014). In Table 8 we list the maximum deviations of observed and expected 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}}-\pi_{j_{1},j_{2},y_{1},y_{2}}(\hat{\boldsymbol{\theta}})|$. From the table, it is revealed, that there is no misfit. The maximum discrepancy occurs between the continuous variables food and leisure. For this bivariate margin, the discrepancy of 509/9960 maximum occurs in the BVN factor copula model, while this drops to 133/9960 in the selected 2-factor copula model.

For the selected 2-factor model based on the proposed algorithm, note that, without the need for a varimax rotation, the unique loadings show that one factor is loaded only on the discrete variables (dishwasher, car, motorcycle, and bicycles), while both factors are loaded on the continuous variables (food, clothes, and leisure). This reveals that the one latent variable which is only associated with the continuous variables measures the expenses, while the the other which is associated with all the mixed variables measures the possession.

7 Simulations

An extensive simulation study is conducted to (a) examine the performance of the diagnostics to show that the correlation matrix of the simulated variables has a factor structure, (b) check the small-sample efficiency of the sample versions of $\rho_{N},\rho_{N}^{+},\rho_{N}^{-}$, (c) gauge the small-sample efficiency of the proposed estimation method and investigate the misspecification of the bivariate pair-copulas, (d) examine the reliability of using the heuristic algorithm to select the correct bivariate linking copulas, and (e) study the small-sample performance of the $M_{2}$ statistic after transforming the continuous and count variables to ordinal.

We randomly generated samples of size $n=\{100,300,500\}$ from each selected one- and two-factor copula models in the three application examples in Section 6. We set the type of the variables, the univariate margins and the bivariate linking copulas, along with their univariate and dependence parameters to mimic the real data. The binary variables don’t have tail asymmetries, hence parametric copulas are less distinguishable. Therefore instead of binary, we simulated from ordinal with 3 equally weighted categories.

Table 9 contains the simulated means and standard deviations (SD) of the discrepancy measures $D_{1}$, $D_{2}$ and $D_{3}$. The resultant summaries show that all the discrepancy measures correctly recognize both that the correlation structure has a factor structure and the number of factors. Among the discrepancy measures, $D_{2}$ has a good performance even for a small sample size ($n=100$), while this is not the case for $D_{1}$ and $D_{3}$ which require larger sample sizes to successively determine the number of adequate factors.

To check the small-sample efficiency of the sample versions of $\rho_{N}$, $\rho_{N}^{+}$ and $\rho_{N}^{-}$ we have generated $10^{4}$ random samples of size $n=\{100,300,500\}$ from all the aforementioned bivariate copulas that join the distributions of two continuous variables, two ordinal variables, one continuous and one ordinal variable, one continuous and one count variable, one ordinal and one count, and two count variables with small ($\tau=0.3$), moderate ($\tau=0.5$) and strong dependence ($\tau=0.7$). Representative results are shown in Table 10 for the Gumbel copula. Note that the count variable was treated as ordinal with 5 categories where the 5th category contained all the counts greater than 3. The resultant biases, root mean square errors (RMSE), and standard deviations (SD), scaled by $n$, show the estimation of the correlations and semi-correlations is highly efficient. Note in passing that because only part of the data are used in computing sample semi-correlations, their variability is larger than the correlations. However, if there is a consistent direction to the tail asymmetry based on semi-correlations, this is useful information for choosing potential bivariate parametric copulas.

[FIGURE:]

Table 11 contains the resultant biases, RMSEs, and SDs, scaled by $n$, for the estimates obtained using the estimation approach in Section 3. The results show that the proposed estimation approach is highly efficient according to the simulated biases, SDs and RMSEs. We further investigated the misspesification of the bivariate pair-copulas by deriving the same statistics but from 1-factor model with BVN pair copulas, i.e. the standard 1-factor model. Once again, the simulated data are based on the selected 1-factor copula models in Section 6. Table 12 contains the resultant biases, RMSEs, and SDs, scaled by $n$. The results show that the Kendall’s tau estimates are not robust to pair-copula misspesification if the true (simulated) factor copula model has different dependence in the middle of the data, e.g. when the BB10 copulas that can provide a non convex shape of dependence (see e.g., Fugure 1) are used to specify the true factor copula model (Table 12, Swiss consumption survey). As we have already mentioned the Kendall’s $\tau$ only accounts for the dependence dominated by the middle of the data, and it is expected to be similar among parametric families of copulas that provide a convex shape of dependence (Table 12, Political-economic dataset and general social survey).

Table 13 contains four common nominal levels of the $M_{2}$ statistic under the factor copula models for mixed data. We transformed the continuous and count variables to ordinal with $K=\{3,4,5\}$ and $K=\{3,4\}$ categories, respectively, using the unsupervised strategies proposed in Section 5.2. We also transformed the count variables to ordinal with $K=5$ categories by treating them as ordinal where the 5th category contained all the counts greater than 3. As the observed levels are close to nominal levels, it is demonstrated that the $M_{2}$ statistic remains reliable for mixed data and that the information loss under transformation to ordinal is minimal.

Table 14 presents the number of times that the true bivariate parametric copulas are chosen over 100 simulation runs. If the true copula has distinct dependence properties with medium to strong dependence, then the algorithm performs extremely well as the sample size increases. Low selection rates occur if the true copulas have low dependence or similar tail dependence properties, since for that case it is difficult to distinguish amongst parametric families of copulas (Nikoloulopoulos and Karlis,, 2008). For example,

• •

in the results from the 2-factor model for the general social survey, the true copula for the first continuous variable (1st factor) is the reflected Gumbel with $\tau=0.34$ and is only selected a considerable small number of times. The algorithm instead selected with a high probability the reflected Joe (results not shown here due to space constraints), because both reflected Joe and Gumbel copulas provide similar dependence properties, i.e., lower tail dependence.

• •

in the results from the 2-factor model for the Swiss consumption survey, the variables with Frank copulas have the lowest selection rates. This is due to the fact that their true Kendall’s $\tau$’s parameters are close to 0 (independence).

8 Discussion

We have extended the factor copula model proposed in Krupskii and Joe, (2013) and Nikoloulopoulos and Joe, (2015) to the case of mixed continuous and discrete responses. It is the most general factor model as (a) it has the standard factor model with an additive latent structure as a special case when the BVN copulas are used, (b) it can have a latent structure that is not additive if other than BVN copulas are called, (c) the parameters of the univariate distributions are separated from the copula (dependence) parameters which are interpretable as dependence of an observed variable with a latent variable, or conditional dependence of an observed variable with a latent variable given preceding latent variables. Other non-linear (e.g., Rizopoulos and Moustaki, 2008), semi- (e.g., Gruhl et al., 2013) or non-parametric models (e.g., Kelava et al., 2017) with latent variables have either an additive latent structure or allow polynomial and interaction terms to be added in the linear predictor, hence are not as general. Another mixed-variable model in the literature that is called factor copula model (Murray et al.,, 2013) is restricted to the MVN copula as the model proposed by Gruhl et al., (2013), hence has an additive latent structure.

We have shown that the factor copula models provide a substantial improvement over the standard factor model on the basis of the log-likelihood principle, Vuong’s and $M_{2}$ statistics. Hence, superior statistical inference for the loading parameters of interest can be achieved. This improvement relies on the fact that the latent variable distribution is expressed via factor copulas instead of the MVN distribution. The latter is restricted to linear and reflection symmetric dependence. Rizopoulos and Moustaki, (2008) stressed that the inadequacy of normally distributed latent variables can be caused by the non-linear dependence on the latent variables. The factor copula can provide flexible reflection asymmetric tail and non-linear dependence as it is a truncated canonical vine copula (Brechmann et al.,, 2012) rooted at the latent variables. Joe et al., (2010) show 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,\ldots,d-1$ to have tail dependence. The 1-factor copula has bivariate copulas with tail dependence in the 1st level and independence copulas in all the remaining levels of the vine (truncated after the 1st level). The 2-factor copula has bivariate copulas with tail dependence in the 1st and 2nd level and independence copulas in all the remaining levels (truncated after the 2nd level). Hence, the tail dependence among the latent variables and each of the observed variables is inherited to the tail dependence among the observed variables.

Even in the cases, where the effect of misspecifying the bivariate linking copula choice to build the factor copula models can be seen as minimal for the Kendall’s $\tau$ (loading) parameters, the tail dependence varies, as explained in Section 2.1, and is a property to consider when choosing amongst different families of copulas and hence affects prediction. Rabe-Hesketh et al., (2003) highlighted the importance of the correct distributional assumptions for the prediction of latent scores. The latent scores will essentially show the effect of different model assumptions, because it is an inference that depends on the joint distribution. The factor copula models have bivariate copulas that link the latent variables to each of the observed variables. If these bivariate copulas have upper or lower tail dependence, then this type of dependence is inherited to the dependence between the factor scores and each of the observed variables. Hence, factor scores are fairly different than the ones for the standard factor model if the sample size is sufficient. Figure 3 demonstrates these differences by revisiting the political-economic dataset in Section 6.1 and comparing the political-economic risk ranking obtained via our selected model, the factor copula model with BVN copulas (standard factor model), and the mixed-data factor analysis of Quinn, (2004). It is revealed that even for a small sample size ($n=62$) there are differences. Between the factor copula model with BVN copulas and the factor analysis model of Quinn, (2004), there are small to moderate differences, because while these models share the same latent variables distribution, the former model does not assume the observed variables to be normally distributed, but rather uses the empirical distribution of the continuous observed variables, i.e. allows the margins to be quite free and not restricted by normal distribution. The differences at the lower panel graph are solely due the miss-specification the latent variable distribution.

As stated by many researchers (e.g., Rabe-Hesketh and Skrondal, 2001, 2004), the major difficulty of all the models with latent variables is identifiability. For example, for the standard factor model or the more flexible model in Irincheeva et al., 2012b one of loadings in the second factor has to be set to zero, because the model with $2d$ loadings is not identifiable. The standard factor model arises as special case of our model if we use as bivariate linking copulas the BVN copulas. Hence, for the 2-factor copula model with BVN copulas, one of the BVN copulas in the second factor has to be set as an independence copula. However, using other than BVN copulas, the 2-factor copula model is near-identifiable with $2d$ bivariate linking copulas as it has been demonstrated by Krupskii and Joe, (2013) and Nikoloulopoulos and Joe, (2015).

Software

Our modelling framework is implemented in the package FactorCopula (Kadhem and Nikoloulopoulos,, 2020) within the open source statistical environment R (R Core Team,, 2020). All the analyses presented in Sections 6.1 and 6.2 are given as code examples in the package.

Acknowledgements

We would like to thank the referees and Professor Harry Joe (University of British Columbia) for their careful reading and comments that have led to an improved presentation and Dr Irina Irincheeva (University of Bern) and Professor Marc Genton (King Abdullah University of Science and Technology) for sharing the Swiss consumption survey dataset. 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.