Efficient Bayesian estimation for flexible panel models for multivariate outcomes: Impact of life events on mental health and excessive alcohol consumption

TL;DR

This paper introduces a novel Bayesian estimation method using particle Metropolis within Gibbs and Hamiltonian Monte Carlo for complex multivariate panel data with mixed outcomes, applied to study life events' effects on mental health and alcohol use.

Contribution

It develops a new Bayesian inference approach for multivariate panel models with mixed outcomes and applies it to analyze the impact of life events on health and behavior.

Findings

Life events significantly affect mental health and alcohol consumption.

Dependence between outcomes is stronger with poor mental health and excessive drinking.

Method improves inference for complex panel data models.

Abstract

The problem we consider considers estimating a multivariate longitudinal panel data model whose outcomes can be a combination of discrete and continuous variables. This problem is challenging because the likelihood is usually analytically intractable. Our article makes both a methodological contribution and also a substantive contribution to the application. The methodological contribution is to introduce into the panel data literature a particle Metropolis within Gibbs method to carry out Bayesian inference, using a Hamiltonian Monte Carlo (Neal 2011} proposal for sampling the vector of unknown parameters. We note that in panel data models the Our second contribution is to apply our method to carry out a serious analysis of the impact of serious life events on mental health and excessive alcohol consumption. The dependence between these two outcomes may be more pronounced when consumption of alcohol is excessive and mental health poor, which in turn has implications for how life events impact the joint distribution of the outcomes.

Tables

Table 1. Table 1: TNV comparison of different sampling schemes (PG, data augmentation, MCMC-MH) for the bivariate probit regression simulation with random effects with P = 1000 𝑃 1000 P=1000 and T = 4 𝑇 4 T=4

	PG	Data Aug.	MH1	MH10	MH20	MH50
Time	$0.62$	$0.13$	$0.48$	$3.21$	$6.56$	$15.13$
$IAC{T}_{mean}$	$4.42$	$44.23$	$42.88$	$7.13$	$5.94$	$4.60$
TNV	$2.74$	$5.75$	$20.58$	$22.89$	$38.97$	$69.60$
Rel. TNV	$1$	$2.09$	$7.51$	$8.35$	$14.22$	$25.40$

Table 2. Table 2: Average Partial Effects on the probability of binge drinking and mental health for major life events variables for male. Symbol ∗ denotes statistical significance

Variables	Probit	IACT	Gaussian	IACT
gave birth	$\underset{(-0.02,0.00)}{-0.01}$	1.17	$\underset{(-0.03,0.00)}{-0.01}$	1.22
death of a friend	$\underset{(-0.00,0.01)}{0.00}$	1.28	$\underset{(-0.01,0.01)}{0.00}$	1.38
death of a spouse/child	$\underset{(-0.01,0.08)}{0.03}$	1.05	$\underset{(-0.02,0.08)}{0.02}$	1.12
personal injury	$\underset{{(0.01,0.04)}^{\ast }}{0.02}$	1.19	$\underset{{(0.01,0.04)}^{\ast }}{0.02}$	1.45
getting married	$\underset{(-0.02,0.01)}{-0.01}$	1.15	$\underset{(-0.03,0.01)}{-0.01}$	1.14
changed residence	$\underset{(-0.01,0.01)}{-0.00}$	1.24	$\underset{(-0.01,0.01)}{0.00}$	1.22
victim of crime	$\underset{(-0.01,0.02)}{0.00}$	1.09	$\underset{(-0.01,0.02)}{0.01}$	1.16
promoted at work	$\underset{(-0.01,0.02)}{0.00}$	1.29	$\underset{(-0.01,0.01)}{0.00}$	1.26
back with spouse	$\underset{(-0.02,0.04)}{0.00}$	1.26	$\underset{(-0.03,0.03)}{-0.00}$	1.38
separated from spouse	$\underset{(-0.00,0.03)}{0.01}$	1.14	$\underset{(-0.00,0.03)}{0.01}$	1.22
improvement in financial	$\underset{(-0.01,0.01)}{-0.00}$	1.04	$\underset{(-0.02,0.02)}{-0.00}$	1.10
worsening in financial	$\underset{(-0.01,0.02)}{0.01}$	1.26	$\underset{(-0.01,0.03)}{0.01}$	1.29
Unconditional Prob.	0.065

Table 3. Table 3: Average Partial Effects for the probability of binge drinking and mental health for major life events variables for Female. Symbol ∗ denotes statistical significance.

Variables	Probit	IACT	Gaussian	IACT
gave birth	$\underset{(-0.03,-0.01)}{-0.02}$	1.60	$\underset{(-0.04,-0.01)}{-0.03}$	1.56
death of a friend	$\underset{(-0.00,0.02)}{0.01}$	1.53	$\underset{(0.00,0.02)}{0.01}$	1.66
death of a spouse/child	$\underset{(-0.02,0.05)}{0.01}$	1.32	$\underset{(-0.03,0.04)}{0.01}$	1.51
personal injury	$\underset{(-0.00,0.02)}{0.01}$	1.49	$\underset{(-0.00,0.02)}{0.01}$	1.59
getting married	$\underset{(-0.01,0.02)}{0.00}$	1.36	$\underset{(-0.01,0.03)}{0.01}$	1.55
changed residence	$\underset{{(0.00,0.01)}^{\ast }}{0.01}$	1.47	$\underset{{(0.00,0.02)}^{\ast }}{0.01}$	1.75
victim of crime	$\underset{(-0.01,0.01)}{-0.00}$	1.46	$\underset{(-0.02,0.01)}{-0.00}$	1.56
promoted at work	$\underset{{(0.01,0.03)}^{\ast }}{0.02}$	1.45	$\underset{{(0.01,0.03)}^{\ast }}{0.02}$	1.59
back with spouse	$\underset{(-0.02,0.03)}{-0.00}$	1.42	$\underset{(-0.03,0.03)}{-0.00}$	1.64
separated from spouse	$\underset{{(0.01,0.04)}^{\ast }}{0.03}$	1.36	$\underset{{(0.01,0.04)}^{\ast }}{0.03}$	1.60
improvement in financial	$\underset{(-0.01,0.02)}{0.01}$	1.49	$\underset{(-0.01,0.02)}{0.01}$	1.44
worsening in financial	$\underset{{(0.01,0.05)}^{\ast }}{0.03}$	1.61	$\underset{{(0.01,0.06)}^{\ast }}{0.03}$	1.74
Unconditional Prob.	0.058

Table 4. Table 4: Estimation results for male Individual Effects and Dependence Parameters. Posterior mean estimates with 95% credible intervals (in brackets).

	Gaussian	Clayton	Gumbel	Probit
${\rho }_{\alpha }$	$\underset{(-0.03,0.10)}{0.03}$	$\underset{(-0.04,0.09)}{0.02}$	$\underset{(-0.03,0.10)}{0.03}$	$\underset{(0.03,0.15)}{0.09}$
Kendall tau	$\underset{(-0.02,0.07)}{0.02}$	$\underset{(-0.03,0.06)}{0.01}$	$\underset{(-0.02,0.06)}{0.02}$	$\underset{(0.02,0.09)}{0.06}$
${\theta }_{dep}$	$\underset{(-0.07,0.11)}{0.02}$	$\underset{(0.11,0.30)}{0.19}$	$\underset{(1.00,1.06)}{1.02}$	$\underset{(-0.10,0.09)}{-0.00}$
Kendall tau	$\underset{(-0.04,0.07)}{0.01}$	$\underset{(0.05,0.13)}{0.09}$	$\underset{(0.00,0.06)}{0.02}$	$\underset{(-0.06,0.06)}{-0.00}$
Lower/Upper Tail	NA	$\underset{(0.00,0.10)}{0.03}$	$\underset{(0.00,0.08)}{0.03}$	NA
${\tau }_1^2$	$\underset{(3.63,5.04)}{4.31}$	$\underset{(3.60,5.09)}{4.29}$	$\underset{(3.71,5.15)}{4.37}$	$\underset{(3.66,5.07)}{4.34}$
${\tau }_2^2$	$\underset{(0.40,0.47)}{0.43}$	$\underset{(0.40,0.47)}{0.43}$	$\underset{(0.40,0.47)}{0.44}$	$\underset{(2.45,3.09)}{2.76}$

Table 5. Table 5: Estimation results for female Individual Effects and Dependence Parameters. Posterior mean estimates with 95% credible intervals (in brackets).

	Gaussian	Clayton	Gumbel	Probit
${\rho }_{\alpha }$	$\underset{(-0.06,0.08)}{0.01}$	$\underset{(-0.08,0.06)}{-0.01}$	$\underset{(-0.06,0.07)}{0.00}$	$\underset{(-0.04,0.08)}{0.02}$
Kendall tau	$\underset{(-0.04,0.05)}{0.00}$	$\underset{(-0.05,0.04)}{-0.01}$	$\underset{(-0.04,0.04)}{0.00}$	$\underset{(-0.02,0.05)}{0.01}$
${\theta }_{dep}$	$\underset{(-0.07,0.11)}{0.02}$	$\underset{(0.11,0.29)}{0.19}$	$\underset{(1.00,1.06)}{1.02}$	$\underset{(-0.08,0.08)}{-0.00}$
Kendall tau	$\underset{(-0.05,0.07)}{0.01}$	$\underset{(0.05,0.13)}{0.09}$	$\underset{(0.00,0.05)}{0.02}$	$\underset{(-0.05,0.05)}{-0.00}$
Lower/Upper tail	NA	$\underset{(0.00,0.09)}{0.03}$	$\underset{(0.00,0.07)}{0.03}$	NA
${\tau }_1^2$	$\underset{(3.35,4.62)}{3.93}$	$\underset{(3.31,4.58)}{3.90}$	$\underset{(3.37,4.69)}{3.99}$	$\underset{(3.38,4.63)}{3.96}$
${\tau }_2^2$	$\underset{(0.34,0.40)}{0.37}$	$\underset{(0.34,0.40)}{0.37}$	$\underset{(0.34,0.40)}{0.37}$	$\underset{(1.91,2.36)}{2.13}$

Table 6. Table 6: Estimation results for male Binge/Excessive Drinking ( y 1 ) subscript 𝑦 1 \left(y_{1}\right) (balanced panel). Posterior mean estimates with 95% credible intervals (in brackets).

	Gaussian	Probit	Clayton	Gumbel
university/degree	$\underset{(-0.97,-0.41)}{-0.69}$	$\underset{(-0.97,-0.40)}{-0.69}$	$\underset{(-0.96,-0.40)}{-0.68}$	$\underset{(-0.97,-0.40)}{-0.68}$
diploma/certificate	$\underset{(-0.22,0.22)}{-0.01}$	$\underset{(-0.23,0.22)}{0.00}$	$\underset{(-0.21,0.23)}{0.01}$	$\underset{(-0.22,0.23)}{0.00}$
married	$\underset{(-0.61,-0.21)}{-0.41}$	$\underset{(-0.61,-0.21)}{-0.41}$	$\underset{(-0.62,-0.22)}{-0.41}$	$\underset{(-0.62,-0.22)}{-0.42}$
income	$\underset{(-0.11,0.02)}{-0.05}$	$\underset{(-0.11,0.02)}{-0.05}$	$\underset{(-0.11,0.01)}{-0.05}$	$\underset{(-0.11,0.02)}{-0.05}$
num. child	$\underset{(-0.27,-0.05)}{-0.16}$	$\underset{(-0.27,-0.05)}{-0.16}$	$\underset{(-0.27,-0.05)}{-0.16}$	$\underset{(-0.27,-0.05)}{-0.16}$
gave birth	$\underset{(-0.66,-0.05)}{-0.35}$	$\underset{(-0.66,-0.06)}{-0.35}$	$\underset{(-0.66,-0.04)}{-0.34}$	$\underset{(-0.67,-0.05)}{-0.36}$
death of a friend	$\underset{(-0.12,0.26)}{0.07}$	$\underset{(-0.12,0.25)}{0.07}$	$\underset{(-0.12,0.25)}{0.07}$	$\underset{(-0.12,0.26)}{0.07}$
death of a spouse/child	$\underset{(-0.62,0.98)}{0.20}$	$\underset{(-0.63,1.00)}{0.21}$	$\underset{(-0.74,1.00)}{0.17}$	$\underset{(-0.74,1.01)}{0.17}$
personal injury	$\underset{(-0.08,0.32)}{0.12}$	$\underset{(-0.08,0.32)}{0.12}$	$\underset{(-0.06,0.34)}{0.14}$	$\underset{(-0.08,0.33)}{0.13}$
getting married	$\underset{(-0.61,0.15)}{-0.22}$	$\underset{(-0.60,0.14)}{-0.22}$	$\underset{(-0.61,0.13)}{-0.23}$	$\underset{(-0.61,0.14)}{-0.23}$
changed residence	$\underset{(-0.11,0.20)}{0.04}$	$\underset{(-0.12,0.19)}{0.04}$	$\underset{(-0.10,0.20)}{0.05}$	$\underset{(-0.12,0.20)}{0.04}$
victim of crime	$\underset{(-0.31,0.24)}{-0.04}$	$\underset{(-0.32,0.22)}{-0.04}$	$\underset{(-0.30,0.24)}{-0.03}$	$\underset{(-0.33,0.23)}{-0.05}$
promoted at work	$\underset{(-0.14,0.30)}{0.08}$	$\underset{(-0.14,0.30)}{0.08}$	$\underset{(-0.15,0.30)}{0.08}$	$\underset{(-0.14,0.30)}{0.08}$
back with spouse	$\underset{(-0.84,0.40)}{-0.21}$	$\underset{(-0.84,0.41)}{-0.21}$	$\underset{(-0.84,0.42)}{-0.20}$	$\underset{(-0.85,0.39)}{-0.21}$
separated from spouse	$\underset{(-0.26,0.33)}{0.03}$	$\underset{(-0.27,0.31)}{0.02}$	$\underset{(-0.27,0.32)}{0.03}$	$\underset{(-0.27,0.31)}{0.02}$
improvement in financial	$\underset{(-0.31,0.30)}{-0.00}$	$\underset{(-0.31,0.32)}{-0.00}$	$\underset{(-0.32,0.31)}{-0.00}$	$\underset{(-0.32,0.32)}{-0.00}$
worsening in financial	$\underset{(-0.43,0.21)}{-0.11}$	$\underset{(-0.45,0.20)}{-0.12}$	$\underset{(-0.43,0.21)}{-0.10}$	$\underset{(-0.44,0.20)}{-0.12}$
${\min }_{1:23}IACT\left({\beta }_{1i}\right)$	1.00	1.00	1.00	1.13
${\max }_{1:23}IACT\left({\beta }_{1i}\right)$	7.58	8.19	8.67	8.64
${\mathrm{mean}}_{1:23}\left({\beta }_{1i}\right)$	2.55	2.55	2.59	2.60

Table 7. Table 7: Estimation results for male Mental Health score ( y 2 ) subscript 𝑦 2 \left(y_{2}\right) (balanced panel). Posterior mean estimates with 95% credible intervals (in brackets).

	Gaussian	Probit	Clayton	Gumbel
university/degree	$\underset{(-0.14,0.03)}{-0.05}$	$\underset{(-0.29,0.07)}{-0.11}$	$\underset{(-0.13,0.03)}{-0.05}$	$\underset{(-0.14,0.03)}{-0.05}$
diploma/certificate	$\underset{(-0.10,0.05)}{-0.02}$	$\underset{(-0.22,0.10)}{-0.06}$	$\underset{(-0.09,0.05)}{-0.02}$	$\underset{(-0.10,0.05)}{-0.03}$
married	$\underset{(-0.12,0.01)}{-0.05}$	$\underset{(-0.28,0.01)}{-0.13}$	$\underset{(-0.12,0.01)}{-0.05}$	$\underset{(-0.12,0.01)}{-0.05}$
income	$\underset{(-0.02,0.03)}{0.01}$	$\underset{(-0.07,0.02)}{-0.02}$	$\underset{(-0.02,0.03)}{0.01}$	$\underset{(-0.02,0.04)}{0.01}$
num. child	$\underset{(-0.03,0.06)}{0.02}$	$\underset{(-0.05,0.10)}{0.03}$	$\underset{(-0.03,0.06)}{0.02}$	$\underset{(-0.03,0.06)}{0.02}$
gave birth	$\underset{(-0.05,0.19)}{0.07}$	$\underset{(-0.11,0.30)}{0.09}$	$\underset{(-0.05,0.19)}{0.07}$	$\underset{(-0.05,0.20)}{0.07}$
death of a friend	$\underset{(-0.05,0.10)}{0.02}$	$\underset{(-0.07,0.19)}{0.06}$	$\underset{(-0.05,0.10)}{0.02}$	$\underset{(-0.05,0.09)}{0.02}$
death of a spouse/child	$\underset{(-0.14,0.49)}{0.18}$	$\underset{(0.02,1.11)}{0.56}$	$\underset{(-0.13,0.50)}{0.18}$	$\underset{(-0.14,0.49)}{0.18}$
personal injury	$\underset{(0.30,0.46)}{0.38}$	$\underset{(0.48,0.75)}{0.62}$	$\underset{(0.30,0.46)}{0.38}$	$\underset{(0.30,0.46)}{0.38}$
getting married	$\underset{(-0.18,0.15)}{-0.02}$	$\underset{(-0.32,0.24)}{-0.04}$	$\underset{(-0.18,0.15)}{-0.01}$	$\underset{(-0.18,0.15)}{-0.01}$
changed residence	$\underset{(-0.06,0.07)}{0.00}$	$\underset{(-0.18,0.05)}{-0.07}$	$\underset{(-0.06,0.07)}{0.00}$	$\underset{(-0.07,0.07)}{0.00}$
victim of crime	$\underset{(0.04,0.27)}{0.16}$	$\underset{(-0.01,0.40)}{0.19}$	$\underset{(0.04,0.27)}{0.16}$	$\underset{(0.04,0.27)}{0.15}$
promoted at work	$\underset{(-0.12,0.06)}{-0.03}$	$\underset{(-0.12,0.19)}{0.04}$	$\underset{(-0.12,0.06)}{-0.03}$	$\underset{(-0.12,0.06)}{-0.03}$
back with spouse	$\underset{(-0.12,0.40)}{0.14}$	$\underset{(-0.11,0.82)}{0.36}$	$\underset{(-0.12,0.41)}{0.14}$	$\underset{(-0.13,0.40)}{0.13}$
separated from spouse	$\underset{(0.10,0.37)}{0.24}$	$\underset{(0.13,0.60)}{0.36}$	$\underset{(0.10,0.37)}{0.23}$	$\underset{(0.10,0.37)}{0.24}$
improvement in financial	$\underset{(-0.13,0.11)}{-0.01}$	$\underset{(-0.27,0.16)}{-0.05}$	$\underset{(-0.14,0.11)}{-0.01}$	$\underset{(-0.13,0.12)}{-0.01}$
worsening in financial	$\underset{(0.24,0.49)}{0.37}$	$\underset{(0.24,0.69)}{0.47}$	$\underset{(0.24,0.49)}{0.37}$	$\underset{(0.24,0.49)}{0.37}$
${\min }_{1:23}IACT\left({\beta }_{2i}\right)$	1.00	1.00	1.00	1.01
${\max }_{1:23}IACT\left({\beta }_{2i}\right)$	3.31	9.92	3.60	2.73
${\mathrm{mean}}_{1:23}\left({\beta }_{2i}\right)$	1.41	2.43	1.30	1.51

Table 8. Table 8: Estimation results for female Binge/Excessive Drinking ( y 1 ) subscript 𝑦 1 \left(y_{1}\right) . Posterior mean estimates with 95% credible intervals (in brackets).

	Gaussian	Probit	Clayton	Gumbel
university/degree	$\underset{(-0.61,-0.16)}{-0.38}$	$\underset{(-0.61,-0.16)}{-0.39}$	$\underset{(-0.62,-0.18)}{-0.40}$	$\underset{(-0.63,-0.17)}{-0.40}$
diploma/certificate	$\underset{(-0.30,0.13)}{-0.09}$	$\underset{(-0.30,0.13)}{-0.09}$	$\underset{(-0.31,0.12)}{-0.09}$	$\underset{(-0.31,0.12)}{-0.09}$
married	$\underset{(-0.68,-0.35)}{-0.51}$	$\underset{(-0.69,-0.35)}{-0.52}$	$\underset{(-0.70,-0.36)}{-0.53}$	$\underset{(-0.70,-0.36)}{-0.53}$
income	$\underset{(0.05,0.19)}{0.12}$	$\underset{(0.05,0.20)}{0.12}$	$\underset{(0.05,0.20)}{0.13}$	$\underset{(0.05,0.20)}{0.13}$
num. child	$\underset{(-0.22,0.00)}{-0.11}$	$\underset{(-0.22,0.00)}{-0.11}$	$\underset{(-0.22,0.01)}{-0.10}$	$\underset{(-0.22,0.00)}{-0.11}$
gave birth	$\underset{(-1.16,-0.49)}{-0.81}$	$\underset{(-1.16,-0.48)}{-0.81}$	$\underset{(-1.15,-0.47)}{-0.81}$	$\underset{(-1.18,-0.48)}{-0.83}$
death of a friend	$\underset{(-0.00,0.37)}{0.19}$	$\underset{(-0.00,0.36)}{0.18}$	$\underset{(-0.01,0.36)}{0.17}$	$\underset{(-0.01,0.37)}{0.18}$
death of a spouse/child	$\underset{(-0.90,0.56)}{-0.16}$	$\underset{(-0.92,0.56)}{-0.16}$	$\underset{(-0.95,0.56)}{-0.17}$	$\underset{(-0.93,0.59)}{-0.16}$
personal injury	$\underset{(-0.43,0.03)}{-0.20}$	$\underset{(-0.43,0.03)}{-0.20}$	$\underset{(-0.42,0.04)}{-0.18}$	$\underset{(-0.44,0.03)}{-0.20}$
getting married	$\underset{(-0.17,0.55)}{0.19}$	$\underset{(-0.17,0.55)}{0.20}$	$\underset{(-0.16,0.55)}{0.20}$	$\underset{(-0.17,0.56)}{0.20}$
changed residence	$\underset{(0.04,0.32)}{0.18}$	$\underset{(0.04,0.31)}{0.17}$	$\underset{(0.04,0.32)}{0.18}$	$\underset{(0.04,0.32)}{0.18}$
victim of crime	$\underset{(-0.43,0.15)}{-0.14}$	$\underset{(-0.44,0.15)}{-0.14}$	$\underset{(-0.43,0.15)}{-0.13}$	$\underset{(-0.44,0.14)}{-0.15}$
promoted at work	$\underset{(0.12,0.53)}{0.33}$	$\underset{(0.12,0.53)}{0.33}$	$\underset{(0.11,0.53)}{0.32}$	$\underset{(0.12,0.53)}{0.33}$
back with spouse	$\underset{(-0.91,0.29)}{-0.30}$	$\underset{(-0.88,0.29)}{-0.31}$	$\underset{(-0.91,0.30)}{-0.30}$	$\underset{(-0.93,0.27)}{-0.32}$
separated from spouse	$\underset{(0.11,0.65)}{0.38}$	$\underset{(0.11,0.65)}{0.38}$	$\underset{(0.11,0.65)}{0.38}$	$\underset{(0.10,0.65)}{0.38}$
improvement in financial	$\underset{(-0.15,0.44)}{0.14}$	$\underset{(-0.17,0.44)}{0.14}$	$\underset{(-0.16,0.45)}{0.14}$	$\underset{(-0.15,0.45)}{0.15}$
worsening in financial	$\underset{(-0.10,0.56)}{0.24}$	$\underset{(-0.10,0.56)}{0.24}$	$\underset{(-0.06,0.61)}{0.27}$	$\underset{(-0.08,0.57)}{0.25}$
${\min }_{1:23}IACT\left({\beta }_{1i}\right)$	1.25	1.01	1.00	1.00
${\max }_{1:23}IACT\left({\beta }_{1i}\right)$	6.55	7.56	6.78	7.58
${\mathrm{mean}}_{1:23}\left({\beta }_{1i}\right)$	2.35	2.41	2.41	2.27

Table 9. Table 9: Estimation results for female Mental Health score ( y 2 ) subscript 𝑦 2 \left(y_{2}\right) (balanced panel) Posterior mean estimates with 95% credible intervals (in brackets).

	Gaussian	Probit	Clayton	Gumbel
university/degree	$\underset{(-0.18,-0.04)}{-0.11}$	$\underset{(-0.39,-0.10)}{-0.24}$	$\underset{(-0.18,-0.04)}{-0.11}$	$\underset{(-0.18,-0.04)}{-0.11}$
diploma/certificate	$\underset{(-0.10,0.03)}{-0.04}$	$\underset{(-0.21,0.07)}{-0.07}$	$\underset{(-0.10,0.03)}{-0.04}$	$\underset{(-0.10,0.03)}{-0.04}$
married	$\underset{(-0.21,-0.09)}{-0.15}$	$\underset{(-0.36,-0.14)}{-0.25}$	$\underset{(-0.21,-0.09)}{-0.15}$	$\underset{(-0.20,-0.09)}{-0.15}$
income	$\underset{(-0.02,0.03)}{0.01}$	$\underset{(-0.04,0.04)}{0.00}$	$\underset{(-0.02,0.03)}{0.01}$	$\underset{(-0.02,0.03)}{0.01}$
num. child	$\underset{(-0.04,0.05)}{0.01}$	$\underset{(-0.05,0.09)}{0.02}$	$\underset{(-0.03,0.06)}{0.01}$	$\underset{(-0.03,0.06)}{0.01}$
gave birth	$\underset{(0.02,0.26)}{0.14}$	$\underset{(0.09,0.47)}{0.28}$	$\underset{(0.03,0.26)}{0.14}$	$\underset{(0.02,0.25)}{0.14}$
death of a friend	$\underset{(-0.04,0.10)}{0.03}$	$\underset{(-0.09,0.14)}{0.02}$	$\underset{(-0.04,0.10)}{0.03}$	$\underset{(-0.04,0.10)}{0.03}$
death of a spouse/child	$\underset{(0.01,0.52)}{0.27}$	$\underset{(0.08,0.97)}{0.52}$	$\underset{(0.01,0.52)}{0.26}$	$\underset{(0.01,0.52)}{0.27}$
personal injury	$\underset{(0.36,0.51)}{0.44}$	$\underset{(0.58,0.85)}{0.72}$	$\underset{(0.36,0.51)}{0.44}$	$\underset{(0.36,0.52)}{0.44}$
getting married	$\underset{(-0.19,0.11)}{-0.04}$	$\underset{(-0.41,0.09)}{-0.16}$	$\underset{(-0.19,0.12)}{-0.04}$	$\underset{(-0.19,0.11)}{-0.04}$
changed residence	$\underset{(-0.04,0.09)}{0.03}$	$\underset{(-0.11,0.09)}{-0.01}$	$\underset{(-0.03,0.09)}{0.03}$	$\underset{(-0.04,0.08)}{0.02}$
victim of crime	$\underset{(-0.04,0.18)}{0.07}$	$\underset{(-0.08,0.30)}{0.11}$	$\underset{(-0.04,0.18)}{0.07}$	$\underset{(-0.05,0.18)}{0.07}$
promoted at work	$\underset{(-0.08,0.10)}{0.01}$	$\underset{(0.02,0.32)}{0.17}$	$\underset{(-0.08,0.11)}{0.01}$	$\underset{(-0.08,0.10)}{0.01}$
back with spouse	$\underset{(-0.01,0.51)}{0.25}$	$\underset{(-0.10,0.81)}{0.35}$	$\underset{(-0.01,0.52)}{0.25}$	$\underset{(-0.01,0.51)}{0.26}$
separated from spouse	$\underset{(0.04,0.30)}{0.17}$	$\underset{(0.11,0.52)}{0.32}$	$\underset{(0.04,0.30)}{0.17}$	$\underset{(0.04,0.30)}{0.17}$
improvement in financial	$\underset{(-0.14,0.09)}{-0.02}$	$\underset{(-0.11,0.26)}{0.07}$	$\underset{(-0.14,0.09)}{-0.02}$	$\underset{(-0.14,0.08)}{-0.02}$
worsening in financial	$\underset{(0.35,0.60)}{0.48}$	$\underset{(0.40,0.86)}{0.63}$	$\underset{(0.35,0.61)}{0.48}$	$\underset{(0.34,0.60)}{0.47}$
${\min }_{1:23}IACT\left({\beta }_{2i}\right)$	1.23	1.00	1.00	1.00
${\max }_{1:23}IACT\left({\beta }_{2i}\right)$	2.48	5.59	3.17	2.94
${\mathrm{mean}}_{1:23}\left({\beta }_{2i}\right)$	1.58	2.09	1.41	1.18

Table 10. Table S1: Comparison of Inefficiency Factors (IACT) for the Parameters with Different Sampling Schemes (PG, data augmentation, MCMC-MH) of Mixed Discrete-Linear Gaussian regression Simulation with random effects P = 1000 𝑃 1000 P=1000 and T = 4 𝑇 4 T=4 .

Param.	PG	MH1	MH10	MH20	MH50
${\beta }_{11}$	$1.58$	$3.64$	$1.10$	$1.34$	$1.00$
${\beta }_{21}$	$1.71$	$4.15$	$1.11$	$1.55$	$1.00$
${\beta }_{31}$	$1.65$	$4.49$	$1.00$	$1.40$	$1.00$
${\beta }_{41}$	$1.82$	$3.34$	$1.10$	$1.53$	$1.00$
${\beta }_{51}$	$1.55$	$3.13$	$1.17$	$1.48$	$1.00$
${\beta }_{61}$	$1.55$	$3.67$	$1.06$	$1.40$	$1.00$
${\beta }_{71}$	$1.63$	$4.48$	$1.00$	$1.39$	$1.01$
${\beta }_{81}$	$1.63$	$3.86$	$1.15$	$1.46$	$1.00$
${\beta }_{91}$	$1.56$	$4.30$	$1.03$	$1.32$	$1.00$
${\beta }_{101}$	$1.55$	$3.69$	$1.07$	$1.36$	$1.00$
${\tau }_1^2$	$30.54$	$398.18$	$59.13$	$46.50$	$28.69$
${\rho }_{\alpha }$	$6.78$	$83.99$	$14.06$	$10.29$	$6.71$

Table 11. Table S2: Comparison of Inefficiency Factors (IACT) for the Parameters with Different Sampling Schemes (PG, data augmentation, MCMC-MH) of Mixed Discrete-Linear Gaussian regression Simulation with random effects P = 1000 𝑃 1000 P=1000 and T = 4 𝑇 4 T=4 .

Param.	PG	MH1	MH10	MH20	MH50
${\beta }_{12}$	$1.51$	$3.50$	$1.15$	$1.21$	$1.00$
${\beta }_{22}$	$1.59$	$3.52$	$1.07$	$1.30$	$1.00$
${\beta }_{32}$	$1.53$	$3.35$	$1.01$	$1.40$	$1.00$
${\beta }_{42}$	$1.59$	$4.71$	$1.00$	$1.45$	$1.00$
${\beta }_{52}$	$1.53$	$4.66$	$1.00$	$1.34$	$1.00$
${\beta }_{62}$	$1.61$	$3.99$	$1.01$	$1.32$	$1.00$
${\beta }_{72}$	$1.52$	$4.06$	$1.00$	$1.27$	$1.00$
${\beta }_{82}$	$1.50$	$3.29$	$1.00$	$1.40$	$1.00$
${\beta }_{92}$	$1.50$	$4.15$	$1.01$	$1.35$	$1.00$
${\beta }_{102}$	$1.58$	$4.82$	$1.01$	$1.27$	$1.00$
${\tau }_2^2$	$1.42$	$12.62$	$1.63$	$1.60$	$1.49$
$\rho$	$7.24$	$13.35$	$9.85$	$7.72$	$7.56$

Table 12. Table S3: TNV comparison of different sampling schemes (PG, data augmentation, MCMC-MH) of Mixed Discrete-Linear Gaussian regression Simulation with random effects P = 1000 𝑃 1000 P=1000 and T = 4 𝑇 4 T=4 .

	PG	MH1	MH10	MH20	MH50
Time	$0.21$	$0.18$	$0.67$	$1.23$	$2.89$
$IAC{T}_{mean}$	$3.24$	$24.46$	$4.40$	$3.90$	$2.65$
TNV	$0.68$	$4.40$	$2.95$	$4.80$	$7.66$
rel. TNV	$1$	$6.47$	$4.34$	$7.06$	$11.26$

Table 13. Table S4: Comparison of Inefficiency Factors (IACT) for the Parameters with Different Sampling Schemes (PG, data augmentation, MCMC-MH) of bivariate Probit regression Simulation with random effects P = 1000 𝑃 1000 P=1000 and T = 4 𝑇 4 T=4

Param.	PG	Data Aug.	MH1	MH10	MH20	MH50
${\beta }_{11}$	$1.00$	$9.73$	$2.13$	$1.00$	$1.00$	$1.17$
${\beta }_{21}$	$1.00$	$11.59$	$2.15$	$1.00$	$1.00$	$1.19$
${\beta }_{31}$	$1.00$	$10.46$	$1.00$	$1.00$	$1.00$	$1.12$
${\beta }_{41}$	$1.00$	$11.01$	$2.42$	$1.01$	$1.00$	$1.19$
${\beta }_{51}$	$1.00$	$11.49$	$2.82$	$1.00$	$1.00$	$1.17$
${\beta }_{61}$	$1.03$	$9.91$	$1.01$	$1.00$	$1.00$	$1.16$
${\beta }_{71}$	$1.00$	$10.39$	$2.49$	$1.00$	$1.00$	$1.15$
${\beta }_{81}$	$1.00$	$11.40$	$2.78$	$1.00$	$1.00$	$1.20$
${\beta }_{91}$	$1.01$	$14.43$	$2.24$	$1.00$	$1.00$	$1.25$
${\beta }_{101}$	$1.00$	$10.44$	$3.89$	$1.00$	$1.01$	$1.20$
${\tau }_1^2$	$18.20$	$79.95$	$178.23$	$19.18$	$20.17$	$17.17$
${\rho }_{\alpha }$	$13.91$	$62.36$	$92.64$	$20.25$	$15.44$	$13.41$

Table 14. Table S5: Comparison of Inefficiency Factors (IACT) for the Parameters with Different Sampling Schemes (PG, data augmentation, MCMC-MH) of bivariate Probit regression Simulation with random effects P = 1000 𝑃 1000 P=1000 and T = 4 𝑇 4 T=4

Param.	PG	Data Aug.	MH1	MH10	MH20	MH50
${\beta }_{12}$	$1.00$	$16.49$	$1.02$	$1.00$	$1.00$	$1.12$
${\beta }_{22}$	$1.00$	$17.03$	$1.00$	$1.09$	$1.00$	$1.07$
${\beta }_{32}$	$1.00$	$19.18$	$4.32$	$1.00$	$1.04$	$1.17$
${\beta }_{42}$	$1.00$	$23.78$	$1.02$	$1.00$	$1.00$	$1.15$
${\beta }_{52}$	$1.00$	$16.50$	$5.52$	$1.01$	$1.03$	$1.19$
${\beta }_{62}$	$1.00$	$21.58$	$1.02$	$1.00$	$1.00$	$1.15$
${\beta }_{72}$	$1.00$	$15.53$	$1.00$	$1.00$	$1.00$	$1.06$
${\beta }_{82}$	$1.00$	$17.08$	1.03	$1.05$	$1.00$	$1.13$
${\beta }_{92}$	$1.00$	$15.37$	$1.00$	$1.00$	$1.00$	$1.11$
${\beta }_{102}$	$1.01$	$17.17$	$2.79$	$1.00$	$1.00$	$1.19$
${\tau }_2^2$	$45.11$	$207.69$	$709.23$	$102.65$	$78.98$	$47.92$
$\rho$	$9.69$	$420.90$	$14.29$	$10.38$	$9.09$	$8.64$

Equations

$$y_{1,it}^{\ast }$$

$${\mathbf{\alpha }}_i$$

$${\mathbf{\Sigma }}_{\mathbf{\alpha }}$$

$${\overline{\mathbf{x}}}_i$$

$$y_{1,it}:=I\left(y_{1,it}^{\ast }>0\right)\text{and}{y}_{2,it}:=I\left(y_{2,it}^{\ast }>0\right).$$

$$p\left(\mathbf{y}|\mathbf{\theta },\mathbf{\alpha }\right)$$

$${\mathbf{\mu }}_{it}$$

$$p\left(\mathbf{y}|\mathbf{\theta },\mathbf{\alpha }\right)$$

$${\left[\left(1-\Phi \left(\frac{{\mu }_{1|2}}{{\sigma }_{1|2}}\right)\right)\varphi \left(y_{2,it}-\left({\mathbf{x}}_{2,it}^{\prime }{\mathbf{\beta }}_{21}+{\overline{\mathbf{x}}}_{2,i}^{\prime }{\mathbf{\beta }}_{22}+{\alpha }_{2,i}\right)\right)\right]}^{1-y_{1,it}}$$

$${\mu }_{1|2}$$

$$c({\mathbf{u}}_{it};{\mathbf{\theta }}_{\mathrm{copula}})\varphi ({\epsilon }_{1,it})\varphi ({\epsilon }_{2,it})$$

$$p\left(\mathbf{y}|\mathbf{\theta },\mathbf{\alpha }\right)$$

$$C_{1|2}\left(u_1|u_2;{\mathbf{\theta }}_{\mathrm{cop}}\right)$$

$$\pi \left(\mathbf{\theta },{\mathbf{\alpha }}_{1:P}\right)$$

$$p(\mathbf{y}|\mathbf{\theta })$$

$$\pi \left(\mathrm{d}{\mathbf{\alpha }}_{1:\mathrm{P}}|\mathbf{\theta }\right)$$

$$=\prod _{i=1}^{P}p\left(\mathrm{d}{\mathbf{\alpha }}_{\mathrm{i}}|\mathbf{\theta },{\mathbf{y}}_{\mathrm{i}}\right),$$

$$S_i^{\mathbf{\theta }}$$

$${\psi }^{\theta }\left({\mathbf{\alpha }}_{1:P}^{1:N}\right)$$

$${\hat{\pi }}_N\left(\mathrm{d}{\mathbf{\alpha }}_{1:\mathrm{P}}|\mathbf{\theta }\right)$$

$${\hat{p}}_N(\mathbf{y}|\mathbf{\theta })$$

$${\tilde{\pi }}_N(\mathbf{k},{\mathbf{\alpha }}_{1:P}^{1:N},\mathbf{\theta })$$

$$q_N({\mathbf{k}}^{\ast },{\mathbf{\alpha }}_{1:P}^{\ast 1:N},{\mathbf{\theta }}^{\ast }|\mathbf{k},{\mathbf{\alpha }}_{1:P}^{1:N},\mathbf{\theta }):=q({\mathbf{\theta }}^{\ast }|\mathbf{\theta })\times {\psi }^{{\theta }^{\ast }}\left({\mathbf{\alpha }}_{1:P}^{\ast 1:N}\right)\times \prod _{i=1}^P{{\overline{w}}_i^{\ast }}^{k_i^{\ast }}.$$

$$\frac{{\tilde{\pi }}_N\left({\mathbf{k}}^{\ast },{\mathbf{\alpha }}_{1:P}^{\ast 1:N},{\mathbf{\theta }}^{\ast }\right)}{q_N\left({\mathbf{k}}^{\ast },{\mathbf{\alpha }}_{1:P}^{\ast 1:N},{\mathbf{\theta }}^{\ast }|\mathbf{k},{\mathbf{\alpha }}_{1:P}^{1:N},\mathbf{\theta }\right)}$$

$$=\frac{{\hat{p}}_N(\mathbf{y}|{\mathbf{\theta }}^{\ast })p({\mathbf{\theta }}^{\ast })}{q({\mathbf{\theta }}^{\ast }|\mathbf{\theta })}.$$

$$\displaystyle\min\Bigg{\{}1,\frac{\widehat{p}_{N}(\bm{y}|\bm{\theta}^{\ast})p(\bm{\theta}^{\ast})}{q(\bm{\theta}^{\ast}|\bm{\theta})}\frac{q(\bm{\theta}|\bm{\theta}^{\ast})}{\widehat{p}_{N}(\bm{y}|\bm{\theta})p(\bm{\theta})}\Bigg{\}},$$

$$||{\mathcal{L}}_N\left\{\mathbf{\theta }\left(s\right),{\mathbf{\alpha }}_{1:P}\left(s\right)\in \cdot \right\}-\pi \left(\cdot \right)|{|}_{TV}\to 0,ass\to \infty ,$$

$$\min \left\{1,\frac{{\tilde{\pi }}_N\left({\mathbf{\theta }}_r^{\ast }|\mathbf{k},{\mathbf{\alpha }}_{1:P}^{\mathbf{k}},{\mathbf{\theta }}_{-r}\right)}{{\tilde{\pi }}_N\left({\mathbf{\theta }}_r|\mathbf{k},{\mathbf{\alpha }}_{1:P}^{\mathbf{k}},{\mathbf{\theta }}_{-r}\right)}\times \frac{q_r\left({\mathbf{\theta }}_r|\mathbf{k},{\mathbf{\alpha }}_{1:P}^{\mathbf{k}},{\theta }_r^{\ast },{\mathbf{\theta }}_{-r}\right)}{q_r\left({{\mathbf{\theta }}^{\ast }}_r|\mathbf{k},{\mathbf{\alpha }}_{1:P}^{\mathbf{k}},{\theta }_r,{\mathbf{\theta }}_{-r}\right)}\right\}$$

$${\tilde{\pi }}_N\left(k_1=l_1,...,k_P=l_P|\mathbf{\theta },{\mathbf{\alpha }}_{1:P}^{1:N}\right)=\prod _{i=1}^P{\overline{w}}_i^{l_i},$$

$$\frac{{\psi }^{\mathbf{\theta }}\left({\mathbf{\alpha }}_{1:P}^{1:N}\right)}{{\prod }_{i=1}^P{m}_i\left({\alpha }_i^{k_i}\right)}$$

$$||{\mathcal{L}}_N\left\{(\mathbf{\theta }\left(s\right),{\mathbf{\alpha }}_{1:P}\left(s\right))\in \cdot \right\}-\Pi \left\{(\mathbf{\theta },{\mathbf{\alpha }}_{1:P})\in \cdot \right\}|{|}_{\mathrm{TV}}\to 0,ass\to \infty .$$

$$p\left(\mathbf{\beta },\mathbf{r}|{\mathbf{\theta }}_{-\mathbf{\beta }},\mathbf{y},{\mathbf{\alpha }}_{1:P}^{\mathbf{k}},\mathbf{k}\right)\propto \exp \left(-H\left(\mathbf{\beta },\mathbf{r}\right)\right)$$

$$H\left(\mathbf{\beta },\mathbf{r}\right):=-\mathcal{L}\left(\mathbf{\beta }\right)+\frac{1}{2}{\mathbf{r}}^{\prime }{M}^{-1}\mathbf{r}$$

Full text

Keywords: Copula; Hamiltonian Monte Carlo; Particle Gibbs; Pseudo marginal method

1 Introduction

Our article considers estimating a flexible longitudinal panel data model with multivariate outcomes that are a combination of discrete and continuous variables. In general, estimating such nonlinear and non-Gaussian longitudinal models is challenging because the likelihood is an integral over the latent individual random effects and the observations are not Gaussian. Our article makes two substantive contributions. First, we introduce into the longitudinal panel data methodology a version of particle Metropolis-within-Gibbs (PMwG) (Andrieu et al., 2010) that allows us to carry out Bayesian inference where the unknown model parameters are generated using a proposal obtained by Hamiltonian Monte Carlo (Neal, 2011). We note that the parameter vector in panel data models is often high dimensional, usually because there are many covariates, so that a Metropolis-Hastings proposal based Hamiltonian Monte Carlo (HMC) can be much more efficient than competitor proposals such as a random walk; see Section 3 for more details.

We show in a simulation study that our PG approach outperforms two other approaches for estimating such panel data models with random effects and intractable likelihoods. The first is the standard data augmentation MCMC as in Albert and Chib (1993). The second is the pseudo marginal Metropolis-Hastings (PMMH) approach in Andrieu and Roberts (2009) and Andrieu et al. (2010).

The motivation for our methodological development, and our second contribution, is to investigate the impact of life events on mental health and excessive alcohol consumption using the Household, Income, and Labour Dynamics in Australia (HILDA) panel data set. In the literature that investigates the impact of life events it is standard to consider a single outcome of interest, e.g. mental health or life satisfaction, and if multiple outcomes are considered, then they are estimated as separate models; see, for example, Lindeboom et al. (2002), Frijters et al. (2011) and Buddelmeyer and Powdthavee (2016). We contend that in many cases, we can gain additional insight by jointly estimating the models for the outcomes. Although it is unsurprising that there is an association between mental health problems and excessive alcohol consumption, establishing the nature of those links requires further research. A natural approach is to attempt to identify causal effects as in Mentzakis et al. (2015), but this relies on the availability of instruments or a natural experiment. We use a reduced form approach but argue that it has the potential to provide complementary and useful evidence; see Kleinberg et al. (2015) for similar arguments. Providing a description of the joint distribution of related outcomes is often of interest and has the potential to inform about causal links, albeit indirectly.

It is common in this literature to simplify the outcomes to allow such joint estimation. For example, Contoyannis and Jones (2004) consider a joint model of a range of lifestyle variables in their study of health and lifestyle choice. Buchmueller et al. (2013) model choices across several insurance types. Both of these studies are representative of the methodology where outcomes of interest are a mix of continuous and various types of discrete variables which are all converted into binary outcomes in order to accommodate estimation as a multivariate probit model (MVP) model using simulated maximum likelihood (SML).

The bivariate probit that accommodates the longitudinal nature of the data serves as our baseline model. This model is also be of independent interest motivated by numerous applications of the MVP model; see for example Atella et al. (2004) and Mullahy (2017). However, in general, there can be associated costs in simplifying the multivariate structure to fit into the MVP framework. Abstracting from the mix of outcomes can obscure interesting features of the relationship. Furthermore, the MVP model imposes a linear correlation structure that may lead to misspecification risk when the existence of an asymmetric or nonlinear dependence structure is plausible. Because the MVP model is limited in these ways, we also consider models that accommodate one continuous outcome and one that is categorical. Here we maintain the assumption of normality for the marginal distributions while allowing both normal and non-normal dependence of error terms, where non-normal dependence is obtained by using copulas. Problems induced by correlated individual effects are addressed with Mundlak type specifications that perform reasonably well in practice; see for example Contoyannis et al. (2004) and Woolridge (2005).

The paper is organised as follows. The panel data model is outlined in Section 2. The Bayesian estimation methodology is given in Section 3. The HILDA dataset is described in Section 4. Section 5 discusses the estimation results. Section 6 concludes. The paper has two appendices. Appendix A provides some proofs of theorem defined in Section 3. Appendix B presents some additional empirical results. The paper also has an online supplement whose sections are denoted as Sections S1, etc.

2 General Panel Data Models with Random Effects

Our motivating example is to investigate the impact of life events on excessive alcohol consumption $(y_{1})$ and mental health $(y_{2})$. To match our motivating example, we describe the following panel data models with two outcomes. The extension to three or more outcomes is conceptually straightforward, but can be much more demanding computationally. Section 2.1 defines a bivariate probit model which serves as a baseline model for comparison with the models discussed later. Section 2.2 then extends the model to accommodate one continuous outcome and one categorical variable. Both of these models impose a linear correlation structure that represents misspecification risk when the existence of asymmetric or non-linear dependence structures is plausible. Then, we extend the model such that we still maintain the assumption of normality for the marginal distributions while introducing non-normal dependence of error terms using copulas and is given in Section 2.3. Lastly, we also define and discuss briefly the Mundlak type specifications which are used for all the models in this section (see Mundlak (1978, pg 615-616); Woolridge (2010, pg 615-616)). Note that the models defined in this section can be applied more generally to any variables of interest that can be combination of discrete and continuous variables and they are not restricted to only the applications in this paper.

2.1 Bivariate Probit Model with Random Effects

We first consider the joint distribution of two binary outcomes given by a bivariate probit model. Let $y_{1,it}$ and $y_{2,it}$ be the two observed binary outcomes, for $i=1,...,P$ people and and $t=1,...,T$ time periods. The bivariate probit model is defined using the following latent variable specification

[TABLE]

In eq. 1, $\bm{x}_{j,it},j=1,2,$ are the exogenous variables which may be associated with the two outcomes (including serious/major life-events) and

[TABLE]

is the average over the sample period of the observations on a subset of the exogenous variables $\bm{x}_{it}^{(M)}$. In eqs. 1 and 2, $\bm{\alpha}_{i}$ is an individual-specific and time invariant random component comprising potentially correlated outcome-specific effects, $\bm{\varepsilon}_{it}$ is the idiosyncratic disturbance term that varies over time and individuals and is assumed to be bivariate normally distributed allowing for contemporaneous correlation but otherwise uncorrelated across individuals and time and also uncorrelated with $\bm{\alpha}_{i}$. In this model $y_{1,it}^{*}$ and $y_{2,it}^{*}$ are unobserved. The observed binary outcomes are defined as

[TABLE]

The explanatory variables $\left(\bm{x}_{it},i=1,\dots,P,t=1,\dots,T\right)$ are also assumed to be exogenous with respect to the individual random effects $\bm{\alpha}_{i}$.

Including terms $\overline{\bm{x}}_{1,i}^{\prime}\bm{\beta}_{12}$ and $\overline{\bm{x}}_{2,i}^{\prime}\bm{\beta}_{22}$ to eq. 1 is called the Mundlak (1978) correction because we can now consider $\left(\alpha_{1i}+\overline{\bm{x}}_{1,i}^{\prime}\bm{\beta}_{12},\alpha_{2i}+\overline{\bm{x}}_{2,i}^{\prime}\bm{\beta}_{12}\right)^{\prime}$ as the $i$th composite random effect which is now potentially correlated with the exogenous covariates $\bm{x}_{it}$.

The joint density, conditional to the vector of individual random effects $\bm{\alpha}_{i}$, is

[TABLE]

where $\bm{\Sigma}_{\rm probit}$ is a $2\times 2$ covariance matrix with ones on the diagonal and $(2y_{1,it}-1)(2y_{2,it}-1)\rho_{\varepsilon}$ on the off diagonal, and

[TABLE]

2.2 Mixed marginal bivariate Model with Random Effects

We next consider an extension to the baseline model eqs. 1 to 4, where we treat one of the outcomes $y_{2,it}=y_{2,it}^{*}$ as an observed continuous variable, with $y_{1,it}$ still discrete as in eq. 4. This applies to the mental health variable where continuous values are available. The joint density conditional on the vector of individual random effects $\bm{\alpha}_{i}$ is

[TABLE]

We assume here and in Section 2.1 that the error term $\bm{\varepsilon}_{it}$ is bivariate normal and consequently impose a linear correlation structure. Section 2.3 maintains the assumption of normality for the marginal distributions of $\bm{\varepsilon}_{it}$, while introducing non-normal dependence of the error terms using bivariate copulas.

2.3 Bivariate Copula Models

Copula based models provide a flexible approach to multivariate modeling because they can: (i) capture a wide range of non-linear dependence between the marginals beyond simple linear correlation; (ii) allow the marginal distributions to come from different families of distributions; and, in particular, (iii) allow the marginal distributions to be a combination of discrete and continuous distributions as in section 2.2. There are many possible parametric copula functions proposed in the Statistics and Econometric literatures, with the choice of parametric copula determining the dependence structure of the variables being analysed. Trivedi and Zimmer (2005) discuss some of the most popular copulas. A major difference between copula distribution functions is the range of their dependence structures. Our article considers the Gaussian, Gumbel and Clayton copulas, which are three of the most commonly used bivariate copulas and they are able to capture a wide range of dependence structures. Note that the baseline model is a Gaussian copula. Appendix S1 gives some details of copula models. For a discussion of copula based models with a combination of discrete and continuous marginals, and their estimation methods see Pitt et al. (2006) and Smith and Khaled (2012). Their method augments the copula model with latent variables which are generated within an MCMC scheme. Note that in our estimation, we do not generate any copula latent variables and we work directly with the conditional density given in Equation (8) below. They also do not consider any individual random effects in their models.

We use the copula framework to obtain a more flexible joint distribution for $\bm{\varepsilon}_{it}$, while assuming that its two marginals are normally distributed. Let $c(\cdot;\bm{\theta}_{\rm copula})$ be a bivariate copula density with $\bm{\theta}_{\rm copula}$ the vector of parameters of the copula. Suppose that $\bm{u}_{it}$ has density $c(\cdot;\bm{\theta}_{\rm copula})$ and define $\varepsilon_{j,it}:=\Phi^{-1}(u_{j,it})$ for $j=1,2$, where $\Phi(\cdot)$ the standard normal cdf. Then the density of $\bm{\varepsilon}_{it}$ is

[TABLE]

The joint density of the observations conditional on the vector of individual random effects for the model in section 2.2 is

[TABLE]

where $u_{1,it}=\Phi\left(-\left(\bm{x}_{1,it}^{\prime}\bm{\beta}_{11}+\overline{\bm{x}}_{1,i}^{{}^{\prime}}\bm{\beta}_{12}+\alpha_{1,i}\right)\right)$, and $u_{2,it}=\Phi\left(y_{2,it}-\left(\bm{x}_{2,it}^{{}^{\prime}}\bm{\beta}_{21}+\overline{\bm{x}}_{2,i}^{{}^{\prime}}\bm{\beta}_{22}+\alpha_{2,i}\right)\right)$.

The conditional distribution function of $U_{1}$ given $U_{2}$ in the copula $C\left(\bm{u};\bm{\theta}_{\rm cop}\right)$ is

[TABLE]

and $c\left(\bm{u};\bm{\theta}_{\rm cop}\right)$ is the density of the copula.

Appendix S1 gives closed form expressions for the conditional copula distribution functions for the bivariate copulas used in our article.

The Pearson correlation coefficient is unsuitable for comparing the dependence structures implied by the different copula models with that of the Gaussian copula, because it only measures linear dependence. Appendix S2 discusses Kendall’s $\tau$ and the upper and lower tail dependence measures that we use in the article.

3 Bayesian Inference and particle Markov chain Monte Carlo samplers

This Section discuss efficient Bayesian inference for the random effect panel data models described in Section 2. Our approach is similar to the particle Markov chain Monte Carlo (PMCMC) approaches of Andrieu et al. (2010). However, the PMCMC approaches in Andrieu et al. (2010) are derived for state space models and the random effect panel data models we are interested in this paper have a different structures, since random effects vary across individual, but they do not change over time. This requires us to derive the PMCMC approaches for the models we are interested in from first principles. The benefit is that the simple particle structure gives straightforward derivations that make the material more accessible than the current PMCMC literature.

Let $\bm{\theta}$ be the parameters in the panel data models described in Section 2. The vector individual random effects is denoted as $\bm{\alpha}_{1:P}$, where $P$ is the number of individuals, and the vector of observations for individual $i$ is denoted by $\bm{y}_{i}$, for $i=1,\ldots,P$, with $\bm{y}_{1:P}$ denoting all the observations in the sample. In Bayesian inference, we are interested in sampling from the posterior density

[TABLE]

where $Z:=p(\bm{y}_{1:P})$ is the marginal likelihood.

The basis of our PMCMC approach is to define a target distribution on an augmented space that includes the parameters and multiple copies of the random effects, which we describe as particles. This target distribution is used to derive two samplers. The first is the Pseudo Marginal Metropolis-Hasting (PMMH) sampler and the second is the Particle Metropolis within Gibbs (PMwG) sampler. The PMMH sampler is an efficient method to sample low dimensional parameters that are highly correlated with the latent states because each PMMH step generates these parameters without conditioning on the states. However, in the panel data models defined in Section 2, the dimension of the parameter space is large so that it is difficult to implement PMMH efficiently. The reason is that it is difficult to obtain good proposals for the parameters that require derivatives of the log-posterior because they cannot be computed exactly and need to be estimated. The efficiency of PMMH then depends crucially on how accurately we can estimate the gradient of the log-posterior. If the error in the estimate of the gradient is too large, then there will be no advantage in using proposals with derivatives information over a random walk proposal (Nemeth et al., 2016). Furthermore, the random walk proposal has a step size proportional to $2.56/\sqrt{d}$, where $d$ is the number of parameters used in the random walk step (Sherlock et al., 2015). This implies that for a large number of parameters, the random walk proposal will move very slowly and so will be very inefficient.

The PMwG sampler generates the parameters conditioning on the latent random effects and the parameters of the models can be sampled in separate Gibbs or Metropolis within Gibbs step. Note that by conditioning on the states for the panel data models we are interested in, we are able to compute the gradient of conditional log-posterior analytically. Our article uses a Hamiltonian Monte Carlo proposal which requires the gradient of the log posterior to sample the high dimensional parameter $\bm{\beta}$. However, this sampler is not very efficient for the parameters that are highly correlated with the latent states. We demonstrate in Section 5 that PMwG sampler performs well for the models that we are interested in this paper.

This section is organised as follows. Section 3.1 discusses the target distribution. Section 3.2 discusses the Pseudo Marginal Metropolis-Hastings (PMMH) sampler. Section 3.3 discusses the Particle Gibbs (PG) and Particle Metropolis within Gibbs (PMwG) samplers. Section 3.4 discusses the Hamiltonian Monte Carlo proposal. Section 3.5 compares our proposed PMwG approach to some alternative approaches and shows that our method can be much more efficient.

3.1 Target Distribution

We first describe Algorithm 1 given below, which constructs a particle approximation to the distribution $\pi\left(\rm d\bm{\alpha}_{1:P}|\bm{\theta}\right)$. Note that all the models in Section 2 have the independence properties

[TABLE]

and

[TABLE]

where the independence property given in Equation 11 will be replicated in our particle appoximation.

Let $\left\{m_{i}\left(\bm{\alpha}_{i}|\bm{\theta},\bm{y}_{i}\right);i=1,...,P\right\}$ be a family of proposal densities that we will use to approximate the corresponding densities $\left\{\pi\left(\bm{\alpha}_{i}|\bm{\theta}\right);i=1,...,P\right\}$. We define

[TABLE]

The next assumption ensures that the proposal densities $m_{i}\left(\bm{\alpha}_{i}|\bm{\theta},\bm{y}_{i}\right)$ can be used to approximate approximate the corresponding densities $\left\{\pi\left(\bm{\alpha}_{i}|\bm{\theta}\right);i=1,...,P\right\}$ in Algorithm 1.

Assumption 1**.**

We assume that $S_{i}^{\bm{\theta}}\subseteq Q_{i}^{\bm{\theta}}$ for any $\bm{\theta}\in\bm{\Theta}$ and $i=1,...,P$

Note that Assumption 1 is always satisfied in our implementation because we use the prior density $p\left(\bm{\alpha}_{i}|\bm{\theta}\right)$ as a proposal density and the prior density is positive everywhere.

The generic Monte Carlo Algorithm 1 proceeds as follows.

To define the joint distribution of the particles given the parameters, let $\bm{\alpha}_{1:P}^{1:N}:=\{\bm{\alpha}_{1}^{1:N},\dots,\bm{\alpha}_{P}^{1:N}\}$. The joint distribution is

[TABLE]

Under Assumption 1, Algorithm 1 yields the approximations to $\pi\left(\rm d\bm{\alpha}_{1:P}|\bm{\theta}\right)$ and the likelihood $p(\bm{y}|\bm{\theta})$ as

[TABLE]

and

[TABLE]

It follows straightforwardly that $\widehat{p}_{N}(\bm{y}|\bm{\theta})$ is an unbiased estimator of the likelihood $p(\bm{y}|\bm{\theta})$. The proof is given in Appendix A.

To obtain particle MCMC schemes to estimate $\pi(\bm{\theta},\bm{\alpha}_{1:P})$, let $\bm{k}:=(k_{1},\dots,k_{P})$, with each $k_{i}\in\{1,\dots,N\}$, and let $\bm{\alpha}_{1:P}^{\bm{k}}:=\{\bm{\alpha}_{1}^{k_{1}},\dots,\bm{\alpha}_{P}^{k_{P}}\}$. For $N\geq 1$, we define the target density

[TABLE]

Appendix A shows that $N^{-P}\pi\left(\bm{\theta},\bm{\alpha}_{1:P}^{\bm{k}}\right)$ is the marginal probability density of $\widetilde{\pi}_{N}\left(\bm{k},\bm{\alpha}_{1:P}^{\bm{k}},\bm{\theta}\right)$. Using the target density in Equation 14, the next two sections consider two particle based methods for carrying Markov chain Monte Carlo in panel data models with an intractable likelihood.

3.2 Pseudo Marginal Metropolis Hastings (PMMH) sampler

The PMMH sampler is a Metropolis Hastings update on the extended space with target density defined in Equation (14) and the proposal density for $\bm{\theta}^{\ast},\bm{k}^{\ast},$ and ${\bm{\alpha}}_{1:P}^{\ast 1:N}$, given their current values $\bm{\theta},\bm{k}$ and $\bm{\alpha}_{1:P}^{1:N}$, as

[TABLE]

Note that this proposal density first samples $\bm{\theta}^{\ast}$ from $q(\bm{\theta}^{\ast}|\bm{\theta})$. The ${\bm{\alpha}}_{1:P}^{\ast 1:N}$ are then sampled from $\psi^{\theta^{\ast}}(\cdot)$. Finally, $\bm{k}^{\ast}$ is sampled from $\prod_{i=1}^{P}{\overline{w}^{\ast}_{i}}^{k^{\ast}_{i}}$.

We now consider the ratio of the extended target eq. 14 and extended proposal eq. 15 to obtain the PMMH acceptance probability for this target and proposal. This ratio is

[TABLE]

Hence, the acceptance probability is

[TABLE]

which is the acceptance probability of a PMMH scheme based on only estimating $\pi(\theta)$.

The next assumption is needed to ensure that the PMMH algorithm converges.

Assumption 2**.**

The MH sampler of the target density $\pi\left(\bm{\theta}\right)$ and proposal density $q\left(\bm{\theta}^{\ast}|\bm{\theta}\right)$ is irreducible and aperiodic.

Theorem 1**.**

Suppose Assumptions 1 and 2 hold. Then the PMMH algorithm with the expanded target density eq. 14 and expanded proposal density eq. 15 generates a sequence $\left\{\bm{\theta}\left(s\right),\bm{\alpha}_{1:P}\left(s\right)\right\}$ of iterates whose marginal distributions $\left\{\mathcal{\mathcal{L}}_{N}\left\{\bm{\theta}\left(s\right),\bm{\alpha}_{1:P}\left(s\right)\in.\right\}\right\}$ satisfy

[TABLE]

where $||\cdot||_{TV}$ is the total variation norm.

The proof is given in Appendix A

3.3 Particle Metropolis within Gibbs (PMwG) sampling

We use the same extended target distribution eq. 14 as before in order to sample from $\pi\left(\bm{\theta},\bm{\alpha}_{1:P}\right)$. Let $\bm{\theta}=\left(\bm{\theta}_{1},...,\bm{\theta}_{R}\right)$ be a partition of the parameter vector $\bm{\theta}$ into $R$ components and let $\bm{\Theta}=\left(\bm{\Theta}_{1},...,\bm{\Theta}_{R}\right)$ be the corresponding partition of the parameter space $\bm{\Theta}$. We will use the notation $\bm{\theta}_{-r}\coloneqq\left(\bm{\theta}_{1},...,\bm{\theta}_{r-1},\bm{\theta}_{r+1},...,\bm{\theta}_{R}\right).$ The Particle Gibbs (PG) and Particle Metropolis within Gibbs (PMwG) algorithm involves the following steps.

It is straightforward to implement Steps (1) and (3). We implement Step (2) using the Conditional Monte Carlo Algorithm 3 given below. Note that Step (1) might appear unusual, but it leaves the augmented target posterior density $\widetilde{\pi}^{N}\left(\bm{k},\bm{\alpha}_{1:P}^{1:N},\bm{\theta}\right)$ invariant. This is related to collapsed Gibbs sampler, see, for example Liu (2001, section 6.7).

Conditional Monte Carlo

The expression

[TABLE]

appearing in the target density eq. 14 is the density of all the variables that are generated by the Monte Carlo algorithm conditional on $\left(\bm{\alpha}_{1:P}^{\bm{k}},\bm{k}\right)$. This is a key element of the PMwG Algorithm 2. This update can be understood as updating $N-1$ Monte Carlo samples together while keeping one Monte Carlo sample fixed in $\widetilde{\pi}_{N}\left(\bm{\alpha}_{1:P}^{1:N}|\bm{\theta}\right)$.

To derive convergence results for the PMwG sampler in Algorithm 2 we require the following assumption.

Assumption 3**.**

The Metropolis within Gibbs sampler that is defined by the proposals $q_{r}\left(\cdot|\bm{\theta}_{r},\bm{\theta}_{-r},\bm{\alpha}_{1:P}\right)$, for $r=1,...,R$, and $\pi\left(\bm{\alpha}_{1:P}|\bm{\theta}\right)$ is irreducible and aperiodic.

Assumption 3 is satisfied in our applications because all the proposals and conditional distributions have strickly positive densities.

Theorem 2** (Convergence of the PMwG sampler).**

For any $N\geq 2$, the PMwG update is a transition kernel for the invariant density $\tilde{\pi}^{N}$ defined in eq. 14. If Assumptions 1 and 3 hold then the PMwG sampler generates a sequence of iterates $\left\{\bm{\theta}\left(s\right),\bm{\alpha}_{1:P}\left(s\right)\right\}$ whose marginal distributions $\left\{\mathcal{\mathcal{L}}_{N}\left\{\bm{\theta}\left(s\right),\bm{\alpha}_{1:P}\left(s\right)\in.\right\}\right\}$ satisfy

[TABLE]

The proof is given in Appendix A.

3.4 Sampling the high-Dimensional parameter vector $\bm{\beta}$ using a Hamiltonian Proposal

This section discusses the Hamiltonian Monte Carlo (HMC) proposal to sample the high dimensional parameter vector $\bm{\beta}$ from the conditional posterior density $p\left(\bm{\beta}|\bm{\theta}_{-\bm{\beta}},\bm{y},\bm{\alpha}_{1:P}^{\bm{k}},\bm{k}\right)$. It can be used to generate distant proposals for the Particle Metropolis within the Gibbs algorithm to avoid the slow exploration behaviour that results from simple random walk proposals.

Suppose we want to sample from a $d$-dimensional distribution with pdf proportional to $\exp\left(\mathcal{L\left(\bm{\beta}\right)}\right)$, where $\mathcal{L}\left(\bm{\beta}\right)=\log p\left(\bm{\beta}|\bm{\theta}_{-\bm{\beta}},\bm{y},\bm{\alpha}_{1:P}^{\bm{k}},\bm{k}\right)$ is the logarithm of the conditional posterior density of $\bm{\beta}$ (up to a normalising constant). In Hamiltonian Monte Carlo (Neal, 2011), we augment an auxiliary momentum vector $\bm{r}$ having the same dimension as the parameter vector $\bm{\beta}$ with the density $p\left(\bm{r}\right)=N\left(\bm{r}|0,\bm{M}\right)$, where $\bm{M}$ is a mass matrix of the momentum and often set to the identity matrix. We define the joint conditional density of $(\bm{\beta},\bm{r})$ as

[TABLE]

where

[TABLE]

is called the Hamiltonian.

In an idealized HMC step, the parameters $\bm{\beta}$ and the momentum variables $\bm{r}$ move continuously according to the differential equations

[TABLE]

where $\nabla_{\bm{\beta}}$ denotes the gradient with respect to $\bm{\beta}$. In a practical implementation, the continuous time HMC dynamics need to be approximated by discretizing time, using a small step size $\epsilon$. We can simulate the evolution over time of $\left(\bm{\beta},\bm{r}\right)$ via the “leapfrog”  integrator, where one step of the leapfrog update is

[TABLE]

Each leapfrog step is time reversible by negating the step size $\epsilon$. The leapfrog integrator provides a mapping $\left(\bm{\beta}^{\ast},\bm{r}^{\ast}\right)\rightarrow\left(\bm{\beta},\bm{r}\right)$ that is both time-reversible and volume preserving (Neal, 2011). It follows that the Metropolis-Hastings algorithm with acceptance probability

[TABLE]

produces an ergodic, time reversible Markov chain that satisfies detailed balance and has stationary density $p\left(\bm{\beta}|\bm{\theta}_{-\bm{\beta}},\bm{y},\bm{\alpha}_{1:P}^{\bm{k}},\bm{k}\right)$ (Liu, 2001; Neal, 1996). Algorithm 4 summarizes a single iterate of the Hamiltonian Monte Carlo method.

The performance of HMC depends strongly on choosing suitable values for $\bm{M}$, $\epsilon$, and $L$. We set $\bm{M}=\widehat{\bm{\Sigma}}^{-1}$, where $\widehat{\bm{\Sigma}}$ is an estimate of the posterior covariance matrix after some preliminary pilot runs of the HMC algorithm. The step size $\epsilon$ determines how well the leapfrog integration can approximate the Hamiltonian dynamics. If we set $\epsilon$ too large, then the simulation error is large yielding a low acceptance rate. However, if we set $\epsilon$ too small, then the computational burden is too high to obtain distant proposals. Similalry, if we set $L$ too small, the proposal will be close to the current value of the parameters, resulting in undesirable random walk behaviour and slow mixing. If $L$ is too large, HMC will generate trajectories that retrace their steps. Our article uses the No-U-Turn sampler (NUTS) with the dual averaging algorithm developed by Hoffman and Gelman (2014) and Nesterov (2009), respectively, that still leaves the target density invariant and satisfies time reversibility to adaptively select $L$ and $\epsilon$, respectively.

Appendices S5 to S8 give the derivatives required by the Hamiltonian dynamics for the panel data models given in section 2.

3.5 Comparing the performance of the PMwG with some other approaches

We now specialize the PMwG sampling scheme described in Algorithm 2, to the bivariate probit model with random effects to obtain Algorithm 5 below; see Appendix S5 for more details.

Let $\bm{\theta}=\left(\bm{\beta}_{1},\bm{\beta}_{2},\rho,\bm{\Sigma}_{\alpha}\right)$ be the set of unknown parameters of interest. We use the following prior distributions: $\rho\sim U\left(-1,1\right)$, $p\left(\bm{\Sigma}_{\alpha}^{-1}\right)\sim{\rm Wishart}\left(v_{0},\bm{R}_{0}\right)$, where $v_{0}=6,\bm{R}_{0}=400I_{2}$, and the prior distribution for the parameters of the covariates is $N\left(0,100\bm{I}_{d}\right)$. All the priors are uninformative.

In Step 2 we transform $\rho$ to $\rho_{\rm un}={\rm tanh}^{-1}(\rho)$ so that $\rho_{\rm un}$ is unconstrained and then use the adaptive random walk method of Garthwaite et al. (2015) that automatically scales univariate Gaussian random walk proposals to ensure that the acceptance rate is around 0.3. The MH acceptance probability is

[TABLE]

We can alternatively replace steps 4 and 5, and sample the latent random effects using the Metropolis-Hastings algorithm, which is denoted by MCMC-MH and is given as

• •

$4^{*})$ Sample $\bm{\alpha}_{i}\sim p\left(\bm{\alpha}_{i}|\bm{\theta}\right)=N\left(0,\bm{\Sigma}_{\alpha}\right)$ for $i=1,\dots,P.$

• •

$5^{*})$ Accept $\bm{\alpha}_{i}$ with acceptance probability given by

[TABLE]

To improve the mixing of the MCMC-MH algorithm, we can run steps $4^{*}$ and $5^{*}$ of MCMC-MH method above for a number of iterations, say $10$, $20$, or $50$ iterations, for each individual random effect. Appendix S9 describes an alternative Gibbs sampling scheme with data augmentation for the bivariate probit with random effects model.

We conducted a simulation study to compare three different approaches to estimation: PG, data augmentation, and MCMC-MH, using the base case bivariate probit model with random effects as the data generating process. To define our measure of the inefficiency of different sampling schemes that takes computing time into account, we first define the Integrated Autocorrelation Time $\left(IACT_{\bm{\theta}}\right)$. For a univariate parameter $\theta$, the IACT is estimated by

[TABLE]

where $\widehat{\rho}_{t}\left(\theta_{1:M}\right)$ denotes the empirical autocorrelation at lag $t$ of $\theta_{1:M}$ (after discarding the burnin period iterates).

A low value of the IACT estimate suggests that the chain mixes well. Here, $L$ is chosen as the first index for which the empirical autocorrelation satisfies $\left|\widehat{\rho}_{t}\left(\theta_{1:M}\right)\right|<2/\sqrt{M}$, i.e. when the empirical autocorrelation coefficient is statistically insignificant. Our measure of inefficiency of the sampling scheme is the time normalised variance

[TABLE]

where $CT$ is the computing time and $IACT_{mean}$ is the average of the IACT’s over all parameters.

For this simulation study, we generated a number of datasets with $P=1000$ people and $T=4$ time periods. The covariates are generated as $x_{1,it},...,x_{10,it}\sim U\left(0,1\right)$, and the parameters are set as

[TABLE]

with $\tau_{1}^{2}=2.5$, $\tau_{2}^{2}=1$, and, $\rho_{\epsilon}=\rho_{\alpha}=0.5$. In the simulation study, the total number of MCMC iterations was 11000, with the first 1000 discarded as burnin. The number of importance samples in the PMwG method was set as $100$.

Table 1 summarises the estimation results and show that the PMwG sampler performs best. Tables S4 and S5 in Appendix S4 show the inefficiency factors (IACT) for each parameters in the bivariate probit model. In terms of TNV, PMwG is more than twice as good as the data augmentation approach and is also $7.51$, $8.35$, $14.22$, and $25.40$ times better than MCMC-MH with $1$, $10$, $20$, and $50$ iterations, respectively. This gain is mostly due to the faster computing time (CT) of the PG over MCMC-MH method. Note that with PG, the computation of importance weights in the Conditional Monte Carlo to sample each individual latent random effects can easily be paralellised. On the other hand, the MCMC-MH approach is a sequential method that may not easily be parallelised. The full Gibbs sampler with the data augmentation approach may not be available for all of the models one might want to consider. The high dimensional parameter vector $\bm{\beta}$ is sampled much more efficiently using Hamiltonian proposals compared to the data augmentation approach, which confirms the usefulness of Hamiltonian proposals for such high dimensional parameters.

We also ran a second simulation study for the mixed discrete-linear Gaussian regression. Appendix S3 reports the results.

4 The Data and their characteristics

4.1 Sample and Variable Definitions

To estimate the impact of life-shock events on the two outcomes, alcohol consumption, especially the propensity to binge drink, and the level of mental health, we use data from Release 14 of the Household, Income, and Labour Dynamics in Australia (HILDA) survey. HILDA is a nationally representative longitudinal survey, which commenced in Australia in 2001, with a survey of 13969 persons in 7682 households, and is conducted annually. Each year, all household members aged 15 years or older were interviewed and considered as part of the sample. Information is collected on education, income, health, life satisfaction, family formation, labour force dynamics, employment conditions, and other economic and subjective well-being. In our paper, the analysis is at the level of the individual, where we include those who are aged 15 years or older who have non-missing information on our two outcome variables: life shock variables, and other independent variables. We use balanced samples. 111The data used in this paper was extracted using the Add-On package PanelWhiz for Stata. PanelWhiz (http://www.PanelWhiz.eu) was written by Dr. John P. Haisken-DeNew (john@PanelWhiz.eu). See Hahn and Haisken-DeNew (2013) and Haisken-DeNew and Hahn (2010) for details.

Following Frijters et al. (2014), the data on mental health status used in this paper is generated from nine questions included in the Short-Form General Health Survey (SF-36), which is available in all the waves.

We construct a mental health score by taking the mean of the responses by the individual and then standardise so that the index has a mean zero and standard deviation one. Lower scores indicate better mental health status. Butterworth and Crosier (2004) provide evidence that the SF-36 data collected in the HILDA survey are valid and can be used as a general measure of physical and mental health status. We then categorise someone with good mental health if their score is below 0, and someone with poor mental health if their score is above 0.

The data on alcohol consumption used in this paper is generated from two questions in the HILDA survey. Subjects are asked to respond to the question: Do you drink alcohol? The second question we considered is related to the problem of binge drinking and is only available in waves 7, 9, 11, and 13. Respondents identified as drinkers from the first question are asked: how often do you have 5 or more (female) or 7 or more (male) standard drinks on any one occasion?

Similarly to Srivastava and Zhao (2010), we define the composite binary variable FREQUENT_BINGE as FREQUENT_BINGE=1 if a male respondent drinks excessive alcohol more than 1 day per week or a female respondent drinks alcohol more than 2 or 3 times a month, and is zero otherwise.

The life-shock indicators are generated from responses in a section of HILDA’s self-completion questionaire. Respondents are told ‘We now would like you to think about major events that have happened in your life over the past 12 months’ and are asked whether any of the following apply to them: (1) Separated from spouse or long-term partner, (2) Serious personal injury or illness to self, (3) Death of spouse/child, (4) Got back together with spouse or long-term partner after a separation, (5) Death of a close friend, (6) Victim of property crime (e.g. theft, housebreaking), (7) Got married, (8) Promoted at work, (9) Major improvement in financial situation, (10) Major worsening in financial situation, (11) Changed residence, (12) Partner or I gave birth to a child.

Other Variables

Other control variables included are marital status (married, single/widow/divorce) and the highest educational qualification attained (degree, diploma/certificate, high school and no qual). single/widow/divorce is the excluded category for marital status. Similarly, high school and no qual (no academic qualification) is excluded for the educational variable. We also include the total number of children below age 18 living in the household, age, and the logarithm of annualised household income. The Mundlak correction contains $\overline{age}$, $\overline{age}^{2}$, $\overline{\log\left(income\right)}$, $\overline{num.child}$.

5 Results and Discussion

5.1 Estimation Results

The PMwG sampling scheme was used to estimate the various models defined in section 2. For each panel data model, we used 11000 MCMC samples of which the first 1000 were discarded as burnin. After convergence, $M=10,000$ iterates $\left\{\bm{\theta}^{\left(m\right)}\right\}$ were collected from which we estimated the posterior means of the parameters as well as their 95% credible intervals. The Bayesian methodology provides information on the entire posterior distribution not only for the parameters of the models, but also other parameters of interest especially the partial effects. We say that the variable of interest is significant if its 95% posterior probability intervals does not cover zero.

Models for men and women were analysed separately throughout. Tables 4 and 5 in appendix B show the estimates for various model specifications for the dependence parameters. The Clayton and Gumbel copula specifications are used for the contemporaneous error terms. The overall pattern of dependence is similar for males and females. The dependence in the individual effects and the error terms of the two outcomes are weak for both male and female models as measured by the Kendall tau, which we denote as $\kappa_{\tau}$. The lower tail dependence based on the Clayton copula is very close to zero for both males and females. This indicates that there is little relationship between the unobservables who are in very good mental health and no excessive alcohol consumption in our data. Furthermore, the upper tail dependence based on the Gumbel copula is also very close to zero for both males and females. This also suggests that there is a weak relationship between having very poor mental health and excessive alcohol consumption after conditioning on the covariates. The estimates of the dependence parameters from the bivariate probit model are similar to those from the Gaussian model specification and are consistent with the expected positive correlation although only one of the correlations is significant. The estimate of $\tau_{2}^{2}$ is much bigger for the bivariate probit model than for the Gaussian model because some information is lost in going to binary variables from continuous variables. Tables 6 to 9 in appendix B give the estimates of the parameters of the main covariates. We do not report the estimates associated with the Mundlak corrections for conciseness. The patterns of the estimates for the covariates in the binge drinking equation $(y_{1})$ from the bivariate probit model and the Gaussian model specification are relatively similar across males and females. However, the estimates of the covariates in the mental health equation $(y_{2})$ are slightly different and the Gaussian model has tighter posterior probability intervals. All the copula models gave similar results. Furthermore, it can be seen that all the parameters $\left(\bm{\beta}_{1},\bm{\beta}_{2}\right)$ are estimated efficiently for all the models. For males, the mean IACT for $\bm{\beta}_{1}$ is $2.55$, $2.55$, $2.59$, and $2.60$ and the mean IACT for $\bm{\beta}_{2}$ is $1.41$, $2.43$, $1.30$, and $1.51$ for the mixed Gaussian, bivariate probit, mixed Clayton, and mixed Gumbel, respectively. For females, the mean IACT for $\bm{\beta}_{1}$ is $2.35$, $2.41$, $2.41$, and $2.27$ and the mean IACT for $\bm{\beta}_{2}$ is $1.58$, $2.09$, $1.41$, and $1.18$ for mixed Gaussian, bivariate probit, mixed Clayton, and mixed Gumbel, respectively. This confirms the usefulness of Hamiltonian Monte Carlo proposals for a high dimensional parameter $\bm{\beta}$.

Our primary interest is in the impact of the shocks on the joint outcomes of binge drinking and poor mental health. For these we compute average partial effects as described in the next section.

5.2 Average Partial Effects

We use an Average Partial Effect (APE) to study how a life event such ‘victim of a crime’ impacts on the association between the joint outcomes of binge drinking and low mental health. Let $A_{it}$ denote the event that person $i$ at time $t$ both binge drinks and has poor mental health.

We define the APE for a particular life event LE as

[TABLE]

where $\pi(\bm{\theta})$ is the posterior density of $\bm{\theta}$ and the superscript $(1)$ in $\bm{x}_{it}^{(1)}$ means that the life event of interest is set to 1, with a similar interpretation for $\bm{x}_{it}^{(0)}$. That is, $\rm{APE}_{\rm LE}$ is the average over all people and time periods of the posterior probability of both binge drinking and poor mental health given the data.

Due to the similarity of the results across copula models, we only present the results for the Gaussian copula and the bivariate probit models. Given the draws $\{\bm{\theta}^{(m)},m=1,\dots,M\}$ from the posterior $\bm{\theta}$, the estimate of the ${\rm APE}_{\rm LE}$ for the bivariate probit model is

[TABLE]

with $\left(\bm{\zeta}_{it}^{(0)}\right)^{(m)}$ defined similarly, $\Sigma^{(m)}=\bm{\Sigma}(\bm{\theta}^{(m)})$ and $\bm{\Sigma}_{\alpha}^{(m)}=\bm{\Sigma}_{\alpha}(\bm{\theta}^{(m)})$.

Tables 2 and 3 summarize the estimates of the APEs for major life events shocks for the probability of binge drinking and low mental health score for both males and females for the bivariate probit and Gaussian copula models. In this case, the sign of the APE provides a clear qualitative interpretation, with a significant positive sign implying a significant increase in the joint probability of binge drinking and having low mental health score, and vice versa. All the IACTs of the APEs, for both males and females, for all life events are very close to $1$, showing that the APEs are estimated efficiently.

Although these APE effects are small, it is necessary to compare them with the unconditional joint probability of binge drinking and poor mental health which are also small for both males and females. For example, the effect of a personal injury for males is a bit over 2 percentage points for both models but when expressed as a percentage of the unconditional probability it is 37% according to the probit estimates and 38% for the Gaussian estimates. The death of a spouse/child also has large relative effects for males but these are not significant. In fact, the only significant APE for the joint probability is for personal injury despite several of the estimates in the marginal models reported in tables 6 and 7 being significant.

In general, the results across the probit and Gaussian models are very similar for both males and females. However, unlike the results for males, several shocks have large and significant effects for females. Being separated from their spouse, changing residence, a worsening financial situation, and having a promotion at work are the shocks that are all significant for females. Each of these increases the joint probability of binge drinking and poor mental health and, using the Gaussian results, have impacts relative to the unconditional probability ranging from 16% for ,a change in residence to 58% for a worsening in the financial position. Finally, the ‘gave birth’ shock has a significant negative association which reduces the joint probability of binge drinking and poor mental health.

6 Conclusions

Based on recent advances in Particle Markov chain Monte Carlo (PMCMC), we demonstrate an approach to estimating flexible model specifications for multivariate outcomes using panel data. We propose a particle Metropolis within Gibbs (PMwG) sampling scheme for Bayesian inference of a flexible for multivariate outcomes using panel data and show that this sampler is more efficient than competing methods.

The panel data methods we develop in this paper also accommodate a mix of discrete and continuous outcomes and in doing so avoid the common approach of reducing all outcomes to binary variables so that a multivariate probit approach is possible. We demonstrate in our application that joint modelling of alcohol consumption and mental health often gave only slightly different results after discretising the outcomes. But given that more general specifications better reflect the discrete outcomes of alcohol consumption and the continuous mental health measure, there is an argument that the bivariate probit model is potentially masking important features of the relationship.

The results in the application are somewhat surprising. Specifying and comparing different copulas was motivated by the belief that the dependence structure between excessive alcohol consumption and poor mental health might potentially be very different in the tails of the distributions. The results indicate that this is not the case in that all three copulas provide qualitatively similar results. Moreover, they indicate that the relationship between alcohol consumption and mental health is weak which is a key reason why differences did not emerge across different copulas. While we have not estimated a formal model allowing two-way causality between excessive alcohol consumption and mental health, if such effects exist then we would expect them to manifest themselves in a positive relationship in our joint estimation. Not finding such a relationship is possibly evidence that the causal effects running both ways between excessive alcohol consumption and mental health are indeed weak or even non-existent. This is not inconsistent with the existing literature where evidence is mixed; see for example Boden and Ferguson (2011). Another possibility is that there are causal effects but they relate to particular subgroups of the population and our models are insufficiently rich to capture the heterogeneity in these effects. We have conducted all analyses for males and females separately and found some differences across these groups but it may be that other sources of heterogeneity may be associated with unobservable rather than observable individual features. We leave this interesting line of work for future research.

Appendices

Appendix A Proofs of Results

The first lemma shows that the estimate $\widehat{p}_{N}(\bm{y}|\bm{\theta})$ given in eq. 13 is an unbiased estimate of the likelihood $p(\bm{y}|\bm{\theta})$.

Lemma 1**.**

$E\left\{\widehat{p}_{N}(\bm{y}|\bm{\theta})\right\}=p(\bm{y}|\bm{\theta})$.

Proof.

From 1 and Steps 1 and 2 of algorithm 1, $E\left(w_{i}^{j}\right)=p(\bm{y}_{i}|\bm{\theta})$, and hence the result follows from equations (10)), (12)) and (13)). ∎

Lemma 2**.**

The marginal distribution of $\widetilde{\pi}_{N}\left(\bm{k},\bm{\alpha}_{1:P}^{\bm{k}},\bm{\theta}\right)$ is given by

[TABLE]

Proof.

We integrate the target density $\widetilde{\pi}_{N}\left(\bm{k},\bm{\alpha}_{1:P}^{1:N},\bm{\theta}\right)$ over $\bm{\alpha}_{1:P}^{\left(-\bm{k}\right)}$

[TABLE]

∎

Proof of Theorem 1.

The proof follows from Assumption 1, Lemmas 1 and 2 and Theorem 1 in Andrieu and Roberts (2009). ∎

Proof of Theorem 2.

The proof follows the approach in Theorem 5 in Andrieu et al. (2010, pg. 300). The algorithm is a Metropolis within Gibbs sampler targeting Equation 14. Hence we focus on establishing irreducibility and aperiodicity. It will be convenient to use the notation $\widetilde{\pi}_{N}\left(\bm{k},\bm{\alpha}_{1:P}^{1:N},\bm{\theta}\right)=\widetilde{\pi}_{N}\left(\bm{k},\bm{\alpha}_{1:P}^{\bm{k}},\bm{\alpha}_{1:P}^{-\bm{k}},\bm{\theta}\right)$ to partition the particles $\bm{\alpha}_{1:P}^{1:N}$ into $\bm{\alpha}_{1:P}^{\bm{k}}$ and $\bm{\alpha}_{1:P}^{-\bm{k}}$ which are the particles selected and not selected by the indices $\bm{k}$ respectively.

Let $\bm{k}\in\{1,\ldots,N\}^{P}$, $D\in{\cal B}\left({\cal R}^{P}\right)$, $E\in{\cal B}\left({\cal R}^{(N-1)P}\right)$ and $F\in{\cal B}\left(\Theta\right)$ be such that $\widetilde{\pi}_{N}\left(\{\bm{k}\}\times D\times E\times F\right)>0$.

From Assumption 1 it is possible to show that accessible sets for the Metropolis within Gibbs sampler are also marginally accessible by the particle Metropolis within Gibbs sampler. From this and Assumption 3, we deduce that there is a finite $j>0$ such that ${\cal L}_{PMwG}\left\{\left(\bm{K}(j),\bm{\alpha}_{1:P}(j),\bm{\theta}(j)\right)\in\{\bm{k}\}\times D\times F\right\}>0$.

Now because Step 2 consists of a Gibbs step using $\widetilde{\pi}_{N}(\cdot)$, we deduce that

[TABLE]

and the irreducibility of the PMwG samper follows.

Aperiodicity can be proved by contradiction since, if the PMwG sample is periodic then from Assumption 1 so is the MwG sampler, which contradicts Assumption 3. The result now follows from Theorem 1 of Tierney (1994). ∎

Appendix B Empirical Results

Online Supplement to ‘Efficient Bayesian estimation for flexible panel models for multivariate outcomes: Impact of life events on mental health and excessive alcohol consumption’

Appendix S1 Gaussian, Clayton and Gumbel Copula Models

The Gaussian copula is

[TABLE]

It can capture both positive and negative dependence and has the full range $\left(-1,1\right)$ of pairwise correlations. where $\Phi_{2}$ is the distribution function of the standard bivariate normal distribution, $\Phi$ is the distribution function for the standard univariate normal distribution, and $\theta$ is the dependence parameter. The dependence structure in the Gaussian copula is symmetric, making it unsuitable for data that exhibits strong lower tail or upper tail dependence.

The baseline model in Section 2.1 is a special case of a Gaussian copula where all the univariate marginal distributions are normally distributed.

The bivariate Clayton copula is

[TABLE]

It can only capture positive dependence, although one can reflect a Clayton copula to model the dependence between $u_{1}$ and $-u_{2}$ instead. The dependence parameter $\theta$ is defined on the interval of $\left(0,\infty\right)$. It is suitable for the data which exhibits strong lower tail dependence and weak upper tail dependence.

The bivariate Gumbel copula is

[TABLE]

For the Gumbel copula, the dependence parameter $\theta$ is defined on the interval $\left[1,\infty\right)$, where $1$ represents the independence case. The Gumbel copula only captures positive dependence. It is suitable for data which exhibits strong upper tail dependence and weak lower tail dependence.

Figure S1 plots 10, 000 draws from each of the three copula models.

The conditional copula distribution functions for the bivariate copulas used in our article can be computed in closed form and are given by

[TABLE]

[TABLE]

[TABLE]

Appendix S2 Dependence Measures

We use Kendall’s tau and upper and lower measures of tail dependence to compare the the dependence structures implied by different copula models because the Pearson (linear) correlation coefficient is not a good measure of general dependence between two random variables as it only detects linear dependence.

Kendall’s tau (Joe, 2015, p. 54), is a popular measure of the degree of concordance between two random variables. Let $\left(U_{1},V_{1}\right)$ and $\left(U_{2},V_{2}\right)$ be two draws from the joint distribution of $U$ and $V$. Kendall’s tau is defined as

[TABLE]

The value of $\kappa_{\tau}$ can vary between $-1$ to $1$ and is zero if the two random variables are independent. The Gumbel and Clayton copulas only capture positive dependence so that $0\leq\kappa_{\tau}\leq 1$ For the Gaussian copula, $-1\leq\kappa_{\tau}\leq 1$. For the copula models we consider, $\kappa_{\tau}$ can be computed in closed form as a function of its copula parameters.

[TABLE]

In many cases, the concordance between extreme (tail) values of random variables is of interest, i.e. the clustering of extreme events in the upper or lower tails. For example, suppose we are interested in the relationship between poor mental health and excessive/binge alcohol consumption or good mental health and no alcohol consumption. This requires a dependence measure for the upper and lower tails of the bivariate distribution. In this case, measures of asymmetric dependence are often based on conditional probabilities. The lower and upper tail dependence measures are defined as (Joe, 2015, p. 62),

[TABLE]

If $\lambda^{U}=0$, then the copula is said to have no upper tail dependence, and if $\lambda^{L}=0$ then the copula is said to have no lower tail dependence. The Gaussian copula has no lower or upper tail dependence. For the Clayton copula, $\lambda^{L}=2^{-1/\theta}>0$ and $\lambda^{U}=0$ so the Clayton copula has no upper tail dependence. For the Gumbel copula, $\lambda^{U}=2-2^{1/\theta}$ and $\lambda^{L}=0$. Hence the Gumbel copula has no lower tail dependence and it has upper tail dependence if and only if $\theta\neq 1$.

Appendix S3 Simulation Mixed Discrete Linear Gaussian Regression

This section provides additional simulation study to compare our particle Metropolis-within-Gibbs approach to the MCMC-MH using mixed discrete linear Gaussian regression model given in Section 2.2. The design is similar to the first experiment in Section 3.5 with $n=1000$ and $T=4$, $x_{1,it},...,x_{10,it}\sim U\left(0,1\right)$, true parameters set as follows:

$\beta_{1}=\left(-1.5,0.1,-0.2,0.2,-0.2,0.1,-0.2,0.1,-0.1,-0.2,0.2\right)^{{}^{\prime}}$,

$\beta_{2}=\left(-0.5,0.1,0.2,-0.2,0.2,0.12,0.2,-0.2,0.12,-0.12,0.12\right)^{{}^{\prime}}$, $\tau_{1}^{2}=1$, $\tau_{2}^{2}=2.5$, and, $\rho_{\epsilon}=\rho_{\alpha}=0.5$. We only compare MCMC-MH and PG methods for this simulation since they can be applied more generally to panel data models with random effects. Tables S1, S2 and S3 summarise the estimation results and show that the PG is still much better than MCMC-MH methods.

Appendix S4 Additional Results on Simulation Bivariate Probit Regression

This section provides additional simulation results on bivariate probit regression models in Section 3.5

Appendix S5 Some Further Details of the Sampling Scheme for the Bivariate Probit Model with Random Effects

Steps 3 and 4 of Algorithm 5 use a Hamiltonian Monte Carlo proposal to sample $\bm{\beta}_{1}$ and $\bm{\beta}_{2}$ conditional on the other parameters and random effects. The HMC requires the gradient of $\log p\left(\bm{y}|\bm{\theta},\bm{\alpha}\right)$ with respect to $\bm{\beta}_{1}$ and $\bm{\beta}_{2}$, where

[TABLE]

[TABLE]

and

[TABLE]

The gradient of $\log p\left(\bm{y}|\bm{\theta},\bm{\alpha}\right)$ with respect to $\bm{\beta}_{2}$ is obtained similarly.

Appendix S6 The Sampling Scheme for the Mixed Marginal Gaussian Regression

with Random Effects

The model is described in Section 2.2. The PG sampling scheme is similar to bivariate probit case and the following derivatives are needed in the HMC step. From this section onwards, we denote $\bm{x}_{j,it}^{{}^{\prime}}$ as all the covariates for the $j$th outcomes and $\eta_{j,it}=\bm{x}_{j,it}^{{}^{\prime}}\bm{\beta}_{j}+\alpha_{j,i}$ for $j=1,2$

[TABLE]

and

[TABLE]

Appendix S7 Gradients for the HMC Sampling Scheme for the Mixed Marginal Clayton Copula Regression

with Random Effects

[TABLE]

[TABLE]

where

[TABLE]

The terms $II_{it}$, $III_{it}$, and $IV_{it}$ are given by

[TABLE]

where

[TABLE]

The PG sampling scheme is similar to that for the bivariate probit case except that in step 2 we work with the unconstrained parameter $\theta_{\rm un}=\log\theta$, where $\theta>0$ for the Clayton copula.

Appendix S8 Gradients for the HMC Sampling Scheme for the Mixed Marginal Gumbel Copula Regression

with Random Effects

[TABLE]

[TABLE]

[TABLE]

The PG sampling scheme is similar to that for the bivariate probit model case except that in step 2 we reparametrize the Gumbel dependence parameter to $\theta_{\rm un}:=\log(\theta-1)$ because $\theta>1$.

Appendix S9 Data Augmentation Bivariate Probit Models with Random Effects

We will work with the augmented posterior distribution

[TABLE]

where

[TABLE]

[TABLE]

where $\bm{\mu}_{it}=\left(\bm{\eta}_{1,it},\bm{\eta}_{2,it}\right)^{T}$, $\bm{\Sigma}_{\epsilon}$ and $\bm{\Sigma}_{\alpha}$ are defined in Equations (3) in Section 2.1, and $p\left(\bm{\theta}\right)$ is the prior distribution for $\bm{\theta}$. The following complete conditional posteriors for the Gibbs sampler can then be derived.

[TABLE]

and

[TABLE]

where $\mu_{1|2}=\left(\bm{x}_{1,it}^{{}^{\prime}}\bm{\beta}_{1}+\alpha_{1,i}\right)+\rho\left(y_{2,it}^{*}-\left(\bm{x}_{2,it}^{{}^{\prime}}\bm{\beta}_{2}+\alpha_{2,i}\right)\right)$, and $\sigma_{1|2}=\sqrt{1-\rho^{2}}$ and $\mu_{2|1}$ and $\sigma_{2|1}$ are defined similarly;

[TABLE]

where $\bm{D}_{\alpha_{i}}=\left(\bm{T}\bm{\Sigma}^{-1}+\bm{\Sigma}_{\alpha}^{-1}\right)^{-1}$ and \bm{d}_{\alpha_{i}}=\bm{\Sigma}^{-1}\sum_{t}\left(\left[\begin{array}[]{c}y_{1,it}^{*}\\ y_{2,it}^{*}\end{array}\right]-\left[\begin{array}[]{cc}\bm{x}_{1,it}^{{}^{\prime}}&\bm{0}\\ \bm{0}&\bm{x}_{2,it}^{{}^{\prime}}\end{array}\right]\left[\begin{array}[]{c}\bm{\beta}_{1}\\ \bm{\beta}_{2}\end{array}\right]\right) for $i=1,...,P$;

[TABLE]

where

[TABLE]

and

[TABLE]

[TABLE]

where

[TABLE]

We use Metropolis within Gibbs to sample $\rho$.