Colombian Women's Life Patterns: A Multivariate Density Regression Approach

TL;DR

This paper introduces a Bayesian multivariate density regression model to analyze Colombian women's socio-demographic factors and their effects on fertility, partnership, and work, revealing the negative impact of family violence.

Contribution

The paper develops a flexible Bayesian density regression model capable of handling mixed responses and complex error structures, applied to study women's life patterns in Colombia.

Findings

Family violence negatively influences women's life choices.

The model effectively captures nonlinear relationships and complex error features.

Insights into socio-demographic factors affecting Colombian women.

Abstract

Women in Colombia face difficulties related to the patriarchal traits of their societies and well-known conflict afflicting the country since 1948. In this critical context, our aim is to study the relationship between baseline socio-demographic factors and variables associated to fertility, partnership patterns, and work activity. To best exploit the explanatory structure, we propose a Bayesian multivariate density regression model, which can accommodate mixed responses with censored, constrained, and binary traits. The flexible nature of the models allows for nonlinear regression functions and non-standard features in the errors, such as asymmetry or multi-modality. The model has interpretable covariate-dependent weights constructed through normalization, allowing for combinations of categorical and continuous covariates. Computational difficulties for inference are overcome through an adaptive truncation algorithm combining adaptive Metropolis-Hastings and sequential Monte Carlo to create a sequence of automatically truncated posterior mixtures. For our study on Colombian women's life patterns, a variety of quantities are visualised and described, and in particular, our findings highlight the detrimental impact of family violence on women's choices and behaviors.

Tables

Table 1. Table B.1: Simulated data set. Summaries of the performance: computational burden, mixing, goodness of fit, and predictive errors in mean and density obtained with the parametric model (first row) and the nonparametric model for different values of J 0 subscript 𝐽 0 J_{0} . Results are reported for the adaptive truncation algorithm based on the ESS and CESS stopping rules.

	$J_0$	$J^{\ast }$	CPU	${\text{ESS}}_{\text{MCMC}}$	${\text{ESS}}_{J^{\ast }}$	LPML (${10}^3$)			${\text{ERR}}_{\text{Mean}}$			${\text{ERR}}_{\text{Dens}}$
						$Z_1$	$Z_2$	$Z_3$	$Z_1$	$Z_2$	$Z_3$	$Z_1$	$Z_2$	$Z_3$
Parametric	1	1	0.66	767.3		-1.17	-0.82	-0.34	5.82	3.26	15.39	163.79	115.32	13.86
${\text{ESS}}_J$	2	13	1.90	495.0	1,125.5	-1.01	-0.76	-0.34	3.93	5.38	6.99	132.55	119.69	7.23
	5	14	3.93	519.6	1,966.7	-0.91	-0.71	-0.34	2.94	4.89	7.41	97.80	104.74	6.59
	10	19	5.29	537.0	1,918.9	-0.94	-0.71	-0.34	2.41	3.12	8.42	90.47	98.56	7.99
	15	26	5.77	543.5	1,266.0	-0.93	-0.73	-0.35	2.94	3.75	9.30	82.99	109.44	8.32
	20	24	5.43	567.8	1,990.3	-0.86	-0.69	-0.34	2.31	3.56	7.92	78.89	95.45	7.88
	30	34	9.04	550.8	2,000.0	-0.85	-0.69	-0.34	2.59	3.67	8.37	78.06	97.60	7.96
${\text{CESS}}_J$	2	14	3.67	495.0	1,989.8	-1.01	-0.76	-0.34	3.89	5.26	6.82	129.35	116.34	7.09
	5	14	3.90	519.6	1,978	-0.91	-0.71	-0.34	2.94	4.89	7.41	97.80	104.79	6.59
	10	17	5.18	537.0	1,905.9	-0.92	-0.71	-0.34	2.35	2.92	8.21	89.79	98.81	7.96
	15	23	6.13	543.5	1,974.9	-0.93	-0.73	-0.35	3.00	3.75	9.35	84.40	109.41	8.38
	20	24	5.51	567.8	1,994.1	-0.86	-0.69	-0.34	2.31	3.56	7.92	78.89	95.45	7.88
	30	34	11.17	550.8	2,000	-0.85	-0.69	-0.34	2.59	3.67	8.37	78.06	97.60	7.96

Table 2. Table B.2: Simulated data set. Alternative scenarios: (i) longer MCMC chain; (ii) n = 1 , 400 𝑛 1 400 n=1,400 ; (iii) no censoring. Summaries of the performance: computational burden, mixing, goodness of fit, and predictive errors in mean and density obtained with the algorithm initialized at J 0 = 15 subscript 𝐽 0 15 J_{0}=15 .

	$J^{\ast }$	CPU	${\text{ESS}}_{\text{MCMC}}$	${\text{ESS}}_{J^{\ast }}$	LPML (${10}^3$)			${\text{ERR}}_{\text{Mean}}$			${\text{ERR}}_{\text{Dens}}$
					$Z_1$	$Z_2$	$Z_3$	$Z_1$	$Z_2$	$Z_3$	$Z_1$	$Z_2$	$Z_3$
(i)	23	10.24	619.1	1,352	-0.93	-0.74	-0.35	2.60	3.85	9.41	80.83	105.49	8.28
(ii)	23	12.15	516.0	1,117.2	-1.75	-1.4	-0.73	2.60	3.86	9.71	87.85	103.60	9.32
(iii)	23	4.42	560.5	1,070.9	-0.88	-0.96	-0.34	2.33	2.13	8.59	75.34	79.74	7.66

Equations

$$\mathbb{P}(Z_{\ast ,\ell }=1|\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })=\mathbb{P}(Y_{\ast ,\ell }>0|\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })\approx \sum _{m=1}^M{\vartheta }^m\sum _{j=1}^J{w}_j^m({\mathbf{x}}_{\ast })\Phi \left(\frac{{\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,\ell )}^m}{\sqrt{{\mathbf{\Sigma }}_{j,(\ell ,\ell )}^m}}\right).$$

$$f({\tilde{z}}_{\ast ,\ell }|\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })$$

$$=\sum _{m=1}^M{\vartheta }^m\sum _{j=1}^J{w}_j^m({\mathbf{x}}_{\ast })\mathrm{logN}({\tilde{z}}_{\ast ,\ell }|{\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,\ell )}^m,{\mathbf{\Sigma }}_{j,(\ell ,\ell )}^m),$$

$$\mathbb{P}({\tilde{Z}}_{\ast ,\ell }\ge (x_{\ast ,1}+1)$$

$$\approx \sum _{m=1}^M{\vartheta }^m\sum _{j=1}^J{w}_j^m({\mathbf{x}}_{\ast })\left(1-\Phi \left(\frac{\log (x_{\ast ,1}+1)-{\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,\ell )}^m}{\sqrt{{\mathbf{\Sigma }}_{j,(\ell ,\ell )}^m}}\right)\right).$$

$$\mathbb{P}(Z_{\ast ,\ell }=1|{\tilde{z}}_{\ast ,{\ell }^{\prime }},\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })\approx \sum _{m=1}^M{\vartheta }^m\sum _{j=1}^J{w}_j^m({\mathbf{x}}_{\ast })\Phi \left(\frac{{\mu }_{j,\ell |{\ell }^{\prime }}^m}{\sqrt{{\sigma }_{j,\ell |{\ell }^{\prime }}^{2m}}}\right)\frac{\mathrm{logN}({\tilde{z}}_{\ast ,{\ell }^{\prime }}|{\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,{\ell }^{\prime })}^m,{\mathbf{\Sigma }}_{j,({\ell }^{\prime },{\ell }^{\prime })}^m)}{f({\tilde{z}}_{\ast ,{\ell }^{\prime }}|\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })},$$

$${\mu }_{j,\ell |{\ell }^{\prime }}^m={\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,\ell )}^m+{\mathbf{\Sigma }}_{j,(\ell ,{\ell }^{\prime })}^m({\mathbf{\Sigma }}_{j,({\ell }^{\prime },{\ell }^{\prime })}^m{)}^{-1}(\log ({\tilde{z}}_{\ast ,{\ell }^{\prime }})-{\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,{\ell }^{\prime })}^m),$$

$${\sigma }_{j,\ell |{\ell }^{\prime }}^{2m}={\mathbf{\Sigma }}_{j,(\ell ,\ell )}^m-({\mathbf{\Sigma }}_{j,(\ell ,{\ell }^{\prime })}^m{)}^2({\mathbf{\Sigma }}_{j,({\ell }^{\prime },{\ell }^{\prime })}^m{)}^{-1},$$

$$f({\tilde{z}}_{\ast ,\ell }|{\tilde{z}}_{\ast ,{\ell }^{\prime }},\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })=$$

$$\sum _{m=1}^M{\vartheta }^m\sum _{j=1}^J{w}_j^m({\mathbf{x}}_{\ast })\mathrm{logN}({\tilde{z}}_{\ast ,\ell }|{\mu }_{j,\ell |{\ell }^{\prime }}^m,{\sigma }_{j,\ell |{\ell }^{\prime }}^{2m})\frac{\mathrm{logN}({\tilde{z}}_{\ast ,{\ell }^{\prime }}|{\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,{\ell }^{\prime })}^m,{\mathbf{\Sigma }}_{j,({\ell }^{\prime },{\ell }^{\prime })}^m)}{f({\tilde{z}}_{\ast ,{\ell }^{\prime }}|\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })}.$$

$$f({\breve{z}}_{\ast ,3}|\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })=\int f({\tilde{z}}_{\ast ,3}|{\tilde{z}}_{\ast ,1},\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })f({\tilde{z}}_{\ast ,1}|\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })d{\tilde{z}}_{\ast ,1}$$

$$\approx \sum _{m=1}^M{\vartheta }^m\sum _{j=1}^J{w}_j^m({\mathbf{x}}_{\ast }){\int }_{-\infty }^{{\breve{z}}_{\ast ,3}}\mathrm{logN}({\tilde{z}}_{\ast ,3}|{\mu }_{j,3|1}^m,{\sigma }_{j,3|1}^{2m})\mathrm{logN}({\tilde{z}}_{\ast ,1}|{\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,1)}^m,{\mathbf{\Sigma }}_{j,(1,1)}^m)d{\tilde{z}}_{\ast ,1},$$

$$\mathbb{P}({\breve{Z}}_{\ast ,3}>x_{\ast ,1}|\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })=\mathbb{P}({\tilde{Z}}_{\ast ,3}+{\tilde{Z}}_{\ast ,1}\ge x_{\ast ,1}+1|\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })$$

$$\approx \sum _{m=1}^M{\vartheta }^m\sum _{j=1}^J{w}_j^m({\mathbf{x}}_{\ast })\int \left(1-\Phi \left(\frac{l(x_{\ast ,1}+1)-{\mu }_{j,3|1}^m}{\sqrt{{\sigma }_{j,3|1}^{2m}}}\right)\right)\mathrm{logN}({\tilde{z}}_{\ast ,1}|{\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,1)}^m,{\mathbf{\Sigma }}_{j,(1,1)}^m)d{\tilde{z}}_{\ast ,1},$$

$$f({\breve{z}}_{\ast ,3}|{\tilde{z}}_{\ast ,2},\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })\approx \sum _{m=1}^M{\vartheta }^m\sum _{j=1}^J{w}_j^m({\mathbf{x}}_{\ast })f({\breve{z}}_{\ast ,3}|{\tilde{z}}_{\ast ,2},{\mathbf{\theta }}_j^m,{\mathbf{x}}_{\ast })\frac{\mathrm{logN}({\tilde{z}}_{\ast ,2}|{\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,2)}^m,{\mathbf{\Sigma }}_{j,(2,2)}^m)}{f({\tilde{z}}_{\ast ,2}|\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })}.$$

$$f({\breve{z}}_{\ast ,3}|{\tilde{z}}_{\ast ,2},{\mathbf{\theta }}_j^m,{\mathbf{x}}_{\ast })={\int }_{-\infty }^{{\breve{z}}_{\ast ,3}}\mathrm{logN}({\breve{z}}_{\ast ,3}-{\tilde{z}}_{\ast ,1}|{\mu }_{j,3|(1,2)}^m,{\sigma }_{j,3|(1,2)}^{2m})\mathrm{logN}({\tilde{z}}_{\ast ,1}|{\mu }_{j,1|2}^m,{\sigma }_{j,1|2}^{2m})d{\tilde{z}}_{\ast ,1},$$

$$\begin{array}{rl}{\mu }_{j,3|(1,2)}^m& ={\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,3)}^m+{\mathbf{\Sigma }}_{j,(3,1:2)}^m{\mathbf{\Sigma }}_{j,(1:2,1:2)}^{-1m}(\log ({\tilde{z}}_{\ast ,1:2})-{\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,1:2)}^m),\\ {\sigma }_{j,3|(1,2)}^{2m}& ={\mathbf{\Sigma }}_{j,(3,3)}^m-{\mathbf{\Sigma }}_{j,(3,1:2)}^m{\mathbf{\Sigma }}_{j,(1:2,1:2)}^{-1m}{\mathbf{\Sigma }}_{j,(1:2,3)}^m.\end{array}$$

$$\mathbb{P}(Z_{\ast ,4}=1|{\breve{z}}_{\ast ,3},\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })\approx \sum _{m=1}^M{\vartheta }^m\sum _{j=1}^J{w}_j^m({\mathbf{x}}_{\ast })\mathbb{P}(Y_{\ast ,4}>0|{\breve{z}}_{\ast ,3},{\mathbf{\theta }}_j^m,{\mathbf{x}}_{\ast })\frac{f({\breve{z}}_{\ast ,3}|{\mathbf{\theta }}_j^m,{\mathbf{x}}_{\ast })}{f({\breve{z}}_{\ast ,3}|\mathbf{x},\mathbf{z},{\mathbf{x}}_{\ast })},$$

$$\mathbb{P}(Y_{\ast ,4}>0|{\breve{z}}_{\ast ,3},{\mathbf{\theta }}_j^m,{\mathbf{x}}_{\ast })f({\breve{z}}_{\ast ,3}|{\mathbf{\theta }}_j^m,{\mathbf{x}}_{\ast })$$

$$={\int }_{-\infty }^{\log ({\breve{z}}_{\ast ,3})}\Phi \left(\frac{{\mu }_{j,4|(1,3)}^m}{\sqrt{{\sigma }_{j,4|(1,3)}^{2m}}}\right)\mathrm{logN}({\breve{z}}_{\ast ,3}-{\tilde{z}}_{\ast ,1}|{\mu }_{j,3|1}^m,{\sigma }_{j,3|1}^{2m})\mathrm{logN}({\tilde{z}}_{\ast ,1}|{\mathbf{x}}_{\ast }{\mathbf{\beta }}_{j,(\cdot ,1)}^m,{\mathbf{\Sigma }}_{j,(1,1)}^m)d{\tilde{z}}_{\ast ,1},$$

$${\tilde{Z}}_{i,1}={\mu }_1^t({\tilde{x}}_{i,1},x_{i,2}^{\ast },x_{i,3}^{\ast })+{ϵ}_{i,1},\text{mbox}fori=1,\dots ,n,$$

$$\mu^{t}_{1}(\tilde{x}_{i,1},x^{*}_{i,2},x^{*}_{i,3})=\left\{\begin{array}[]{ll}-0.057\tilde{x}_{i,1}^{2}+3.08\tilde{x}_{i,1}-21.247&\mbox{if }x^{*}_{i,2}\neq 1,x^{*}_{i,3}=2\\ \frac{1}{3}\tilde{x}_{i,1}+10&\mbox{if }x^{*}_{i,2}\neq 1,x^{*}_{i,3}=1\\ 0.0001\tilde{x}_{i,1}^{3}-0.0695\tilde{x}_{i,1}^{2}+3.83\tilde{x}_{i,1}-30.584&\mbox{if }x^{*}_{i,2}=1,x^{*}_{i,3}=2\\ \frac{8}{15}\tilde{x}_{i,1}+7&\mbox{if }x^{*}_{i,2}=1,x^{*}_{i,3}=1\end{array}\right..$$

$$\tilde{Z}_{i,2}=\left\{\begin{array}[]{ll}-0.056\tilde{x}_{i,1}^{2}+3.08\tilde{x}_{i,1}-18+0.75\left[\tilde{z}_{i,1}-\mu^{t}_{1}(\tilde{x}_{i,1},x^{*}_{i,2},x^{*}_{i,3})\right]+\epsilon_{i,2}&\mbox{if }x^{*}_{i,3}=2\\ 0.5\tilde{x}_{i,1}+8+0.75\left[\tilde{z}_{i,1}-\mu^{t}_{1}(\tilde{x}_{i,1},x^{*}_{i,2},x^{*}_{i,3})\right]+\epsilon_{i,2}&\mbox{if }x^{*}_{i,3}=1\end{array}\right.,$$

$$\epsilon_{i,2}\sim\left\{\begin{array}[]{ll}0.9{\rm N}(-\frac{1}{6},0.4^{2})+0.1{\rm N}(1.5,0.75^{2})&\mbox{if }x^{*}_{i,3}=2\\ 0.9{\rm N}\left(-\frac{1}{6},\left(\frac{7.5}{\tilde{x}_{i,1}}\right)^{2}\right)+0.1{\rm N}\left(1.5,\left(\frac{7.5}{\tilde{x}_{i,1}}\right)^{2}\right)&\mbox{if }x^{*}_{i,3}=1\end{array}\right..$$

$$Z_{3,i}\sim \mathrm{Bern}\left(\Phi \left(\frac{{\tilde{x}}_{1,i}-18}{6}\right)\right).$$

$$\begin{array}{rl}{z}_{\ell }& =h_{\ell }(\mathbf{y},\mathbf{x})=c_{\ell }(\mathbf{y},\mathbf{x})\lfloor \exp (y_{\ell })\rfloor ,\text{for }\ell =1,2,\\ z_3& =h_3(\mathbf{y},\mathbf{x})=\text{mathbbm}{1}_{[0,\infty )}(y_3),\end{array}$$

$$(l_{\ell },u_{\ell })$$

$$(l_3,u_3)$$

$$\mathbb{E}[{\mathbf{\beta }}_j]={\mathbf{\beta }}_0=\hat{\mathbf{\beta }}\text{and}\mathbb{E}[{\mathbf{\Sigma }}_j]=\frac{1}{\nu -b-1}{\mathbf{\Sigma }}_0=\hat{\mathbf{\Sigma }}.$$

$${\text{LPML}}^{\ell }=\sum _{i=1}^n\log ({\text{CPO}}_i^{\ell })\text{with}{\text{CPO}}_i^{\ell }={\left(\frac{1}{M}\sum _{m=1}^M\frac{1}{f(z_{i,\ell }|{\mathbf{w}}^m,{\mathbf{\psi }}^m,{\mathbf{\theta }}^m,{\mathbf{x}}_i)}\right)}^{-1},$$

$${\text{ERR}}_{\text{Mean}}^{\ell }=\frac{100}{n^{\ast }}\sum _{i=1}^{n^{\ast }}\frac{|{\mu }_{\ell }^t({\mathbf{x}}_i^{\ast })-{\hat{\mu }}_{\ell }({\mathbf{x}}_i^{\ast })|}{|{\mu }_{\ell }^t({\mathbf{x}}_i^{\ast })|},$$

$${\text{ERR}}_{\text{Dens}}^{\ell }=\frac{100}{n^{\ast }}\sum _{i=1}^{n^{\ast }}\frac{\int |f^t(z_{\ell }^{\ast }|{\mathbf{x}}_i^{\ast })-\hat{f}(z_{\ell }^{\ast }|{\mathbf{x}}_i^{\ast })|d{z}_{\ell }^{\ast }}{\int |f^t(z_{\ell }^{\ast }|{\mathbf{x}}_i^{\ast })|d{z}_{\ell }^{\ast }}\approx \frac{100}{n^{\ast }}\sum _{i=1}^{n^{\ast }}\sum _{g=1}^G|f^t(z_{g,\ell }^{\ast }|{\mathbf{x}}_i^{\ast })-\hat{f}(z_{g,\ell }^{\ast }|{\mathbf{x}}_i^{\ast })|\Delta ,$$

Full text

Supplementary Material for “Colombian Women’s Life Patterns: A Multivariate Density Regression Approach”

Sara Wadelabel=e1]sara.wade@ed.ac.uk [   

Raffaella Piccarretalabel=e2]raffaella.piccarreta@unibocconi.it [   

Andrea Cremaschilabel=e3]andrea.cremaschi@yale-nus.edu.sg [   

Isadora Antoniano-Villalobos label=e4]isadora.antoniano@unive.it [ School of Mathematics, University of Edinburgh, Edinburgh , UK

Department of Decision Sciences, Dondena Research Centre, Bocconi Institute for Data Science and Analytics (BIDSA), Bocconi University, Milan, Italy

Yale-NUS College, Singapore

Department of Environmental Sciences, Informatics and Statistics, Ca’ Foscari University of Venice and BIDSA, Italy

\startlocaldefs\startlocaldefs

, , , and

\addappheadtotoc

In this Supplementary Material, we describe how to compute predictive quantities of interest from the MCMC output in Section A, and in Sections B and C, we report additional results for the simulated data example and the application to the Colombian women dataset described in the paper.

Appendix A Predictions

The weighted posterior samples obtained with the adaptive truncation algorithm can be used to produce various posterior and predictive quantities of interest. Here, we describe how to compute the predictive densities, medians, probability of censoring for the (undiscretized) ages at event and probability of success for binary variables. Full details and implementation for other quantities are provided in the accompanying software and documentation: https://github.com/sarawade/BNPDensityRegression_AdaptiveTruncation.

Focusing on the application in Section 6, for $\ell=1,2,3$, we denote by $\tilde{Z}_{\ell}$ the (undiscretized) age at sexual debut, the (undiscretized) age at union, and the time from sexual debut to first child, respectively. These are linked to our model by the relation $\tilde{Z}_{\ell}=\exp(Y_{\ell})$, and the corresponding ages are obtained through discretization. The (undiscretized) age at first child is denoted as $\breve{Z}_{3}=\tilde{Z}_{1}+\tilde{Z}_{3}$. For Work Status, we have $Z_{4}=\mathbbm{1}_{(0,\infty)}(Y_{4})$. In the following, let $J$ denote the final truncation level, with corresponding weighted particles $({\bf w}_{1:J}^{m},\bm{\theta}_{1:J}^{m},\bm{\psi}_{1:J}^{m},{\bf y}_{1:n}^{m})$ and unnormalized particle weights $\tilde{\vartheta}^{m}$, for $m=1,\ldots,M$ (without loss of generality, we drop the subscript $J$). We indicate with $\vartheta^{m}$, for $m=1,\ldots,M$, the normalized particle weights.

We begin with marginal predictive quantities of interest. First, the predictive probability of success for a binary response given ${\bf x}_{*}$ (e.g. for $\ell=4$, shown in Figure 5, bottom row) is:

[TABLE]

For $\ell=1,2,3$, the marginal predictive density of $\tilde{Z}_{*,\ell}$ given ${\bf x}_{*}$, shown in Figure 6 for some values of ${\bf x}_{*}$, is given by:

[TABLE]

for $\tilde{z}_{*,\ell}>0$, where $\bm{\beta}_{j,(\cdot,\ell)}^{m}$ denotes the $\ell$-th column of $\bm{\beta}$ in component $j$ and particle $m$; $\bm{\Sigma}_{j,(\ell,\ell)}^{m}$ denotes element $(\ell,\ell)$ of the matrix $\bm{\Sigma}$ in component $j$ and particle $m$; and ${\rm logN}(\cdot|\mu,\sigma^{2})$ denotes the log-normal density with parameters $\mu$ and $\sigma^{2}$. Due to the skewness of the predictive densities in our application (Section 6), we focus on the predictive median over the predictive mean to better represent the central tendencies and summarize the predictive densities. The marginal predictive median (Figure 5) can be computed numerically by evaluating the marginal predictive density (1) on a sufficiently dense grid of $\tilde{z}_{*,\ell}$ values.

For $\ell=1,2$ corresponding to age at sexual debut and union, an interesting quantity is the predictive probability that the indexed event has not yet occurred for a new individual with $x_{*,1}$ years of age (Figure C.8), computed as:

[TABLE]

This can be interpreted as the predictive probability of censoring of the event for a new individual and corresponds to the mass above the dashed line of Figure 6, given $x_{*,1}$.

Our model also recovers the joint relationship between responses, which allows inference on conditional properties. Specifically, when $\ell$ indexes a binary response and $\ell^{\prime}$ indexes an age at event response, the conditional predictive probability of success given $\tilde{z}_{*,\ell^{\prime}}$ and ${\bf x}_{*}$ (Figure 8) is:

[TABLE]

where

[TABLE]

[TABLE]

and the density in the denominator is the marginal predictive of equation (1). For $\ell\neq\ell^{\prime}$ both indexing ages at event, the conditional predictive density of $\tilde{Z}_{*,\ell}$ given $\tilde{z}_{*,\ell^{\prime}}$ and ${\bf x}_{*}$ takes the form:

[TABLE]

Figure C.11 shows the conditional predictive density of $\tilde{Z}_{*,2}-\tilde{z}_{*,1}$ given $\tilde{z}_{*,1}$ and ${\bf x}_{*}$, which can be easily computed from (3). The corresponding predictive medians (Figure 7) can be obtained numerically from evaluations of this density on an adequate, dense grid of values. The conditional density plot for $\tilde{Z}_{*,3}$ given $\tilde{z}_{*,1}$ (Figure C.12) and the corresponding median (Figure C.9) can be obtained directly from equation (3).

Finally, we note that for the (undiscretized) age at first child, $\breve{Z}_{*,3}=\tilde{Z}_{*,1}+\tilde{Z}_{*,3}$, and more generally constrained responses, the corresponding marginal and conditional predictive quantities may require integration over $\tilde{Z}_{*,1}$. For example, the conditional predictive density of $\breve{Z}_{*,3}$ given $\tilde{z}_{*,1}$ is simply the conditional predictive density of equation (3), evaluated at $\breve{z}_{*,3}-\tilde{z}_{*,1}$. While the marginal predictive density of $\breve{Z}_{*,3}$ given ${\bf x}_{*}$ (Figure 6) is obtained as:

[TABLE]

where $\tilde{z}_{*,3}=\breve{z}_{*,3}-\tilde{z}_{*,1}$. We evaluate the integral stochastically, via a Monte Carlo approximation, and compute the marginal predictive median of the undiscretized age at first child (Figure 5) numerically from the marginal predictive density in (4). Also, the predictive probability that the woman has not yet had a child at $x_{*,1}$ years of age (Figure C.8) takes the form:

[TABLE]

where $l(z)=\log(\max(0,z-\tilde{z}_{*,1}))$. The conditional predictive density of the (undiscretized) age at first child $\breve{Z}_{*,3}$ given the (undiscretized) age at union $\tilde{z}_{*,2}$ and $x_{*,1}$ is:

[TABLE]

Notice that this expression differs from equation (3) in that

[TABLE]

where

[TABLE]

Figure C.13 shows the conditional predictive density of $\breve{Z}_{*,3}-\tilde{z}_{*,2}$ given $\tilde{z}_{*,2}$ and ${\bf x}_{*}$, which can be easily computed from (5), with the corresponding predictive medians in Figure C.10. Lastly, the conditional predictive probability of success for a binary response, e.g. $\ell=4$ in our application, given $\breve{z}_{*,3}$ and ${\bf x}_{*}$ is:

[TABLE]

where

[TABLE]

where $\mu_{j,4|(1,3)}$ and $\sigma_{j,4|(1,3)}^{2}$ are calculated analogously to expression (6).

Appendix B Simulation study

We assess the performance of the proposed procedure on a simulated data set including $q^{*}=3$ covariates and $d=3$ responses. The first covariate mimics Age and, as such, is assumed to be registered at a discrete level: $x_{1}=\lfloor\tilde{x}_{1}\rfloor$, where $\tilde{x}_{1}\sim{\rm U}(15,30)$. The remaining covariates are categorical: $x_{2}^{*}$ has three levels with probabilities $(0.5,0.3,0.2)$ while $x_{3}^{*}$ has two levels with probabilities $(0.4,0.6)$.

We generate two positive discretized responses and one binary response. The first response, $Z_{1}$, is a discretized noisy observation of a nonlinear function of $x_{1}$. To build $Z_{1}$, we first generate:

[TABLE]

where $\epsilon_{1,1},\dots,\epsilon_{n,1}\stackrel{{\scriptstyle{\scriptstyle i.i.d.}}}{{\sim}}0.9{\rm N}(-15/90,0.5^{2})+0.1{\rm N}(1.5,0.75^{2})$, and

[TABLE]

Similarly, $Z_{2}$ is a discretized noisy observation of a nonlinear function of $x_{1}$ and the realized $z_{1}$ and it is built by generating:

[TABLE]

where the errors are assumed to depend also on $\tilde{x}_{1}$ and $x^{*}_{3}$:

[TABLE]

Note that the response curves are the same for $x_{2}^{*}=2,3$ and differ for other categorical combinations, while the errors are not normal but right skewed, additionally depending on $x_{1}$ and $x_{3}^{*}$ for the second response. Observed responses are set to missing for censored observations, defined as individuals with $\tilde{z}_{1,i}>\tilde{x}_{1,i}$ or $\tilde{z}_{2,i}>\tilde{x}_{1,i}$. Since the age-related variables in our motivating application are registered at a discrete level, the observed responses were rounded down to the nearest integer, i.e. $z_{1}=\lfloor\tilde{z}_{1}\rfloor$, $z_{2}=\lfloor\tilde{z}_{2}\rfloor$. The true curves and densities are depicted in Figure B.1 (top row) for selected combinations of the covariates. The behavior of $Z_{1}$ and $Z_{2}$ may at a first sight seem a simple function of $x_{1}$ alone, however it depends also on the interactions between $x_{1}$ and the categorical covariates. Moreover, this consideration holds both for the regression function and for the variance. Combined together these aspects, which are of the same nature as those present in our motivating data, pose challenges for parametric and semi-parametric models. Indeed, the relationship between each response and the $Age$ at interview is relatively smooth, but interactions with the categorical covariates and changes in variability increase complexity. Finally, a binary response is simulated from a linear probit model depending only on $x_{1}$ (Figure B.2):

[TABLE]

We seek to recover the conditional distribution of the response variables given the covariates using our proposed model, from a sample of size $n=700$. We define the link functions $h_{\ell}({\bf y},{\bf x})$ as:

[TABLE]

where $c_{\ell}({\bf y},{\bf x})=\mathbbm{1}_{(0,x_{1}+1)}(\exp(y_{\ell}))$. In this case, the bounds required in the adaptive MCMC are obtained from inverting $z_{\ell}=h_{\ell}({\bf y},{\bf x})$; concretely,

[TABLE]

Prior specification and computational details.

Prior parameters for the linear coefficients and covariance matrix of each component are specified empirically based on multivariate linear regression fit to the data. Specifically, for $\ell=1,2$, we set $y_{i,\ell}=(l_{i,\ell}+u_{i,\ell})/2$ and $y_{i,\ell}=\log(x_{i,1}+2)$ for uncensored and censored observations, respectively, where the bounds $l_{i,\ell}$ and $u_{i,\ell}$ are defined in Section 4. Additionally, we let $y_{i,3}=-1$ for $z_{i,3}=0$ and $y_{i,3}=1$ for $z_{i,3}=1$. A multivariate linear regression fit on these auxiliary responses gives estimates $\widehat{\bm{\beta}}$ of the linear coefficients and $\widehat{\bm{\Sigma}}$ of the covariance matrix. We then define

[TABLE]

Together, ${\bf U}$ and $\bm{\Sigma}_{j}$ reflect the variability of $\bm{\beta}_{j}$ across components, and we set the matrix ${\bf U}$ such that \min({\rm diag}(\widehat{\bm{\Sigma}}))\,{\bf U}=10({\bf X}^{\mathchoice{\raisebox{0.0pt}{\displaystyle\intercal}}{\raisebox{0.0pt}{\textstyle\intercal}}{\raisebox{0.0pt}{\scriptstyle\intercal}}{\raisebox{0.0pt}{\scriptscriptstyle\intercal}}}{\bf X})^{-1}. The factor of this g-prior was selected to ensure reasonable ages (i.e. mostly lower than 100) in prior simulations. We explored more uninformative and vague prior choices but found that this could lead to quite large and unreasonable imputed ages for censored data. We further set $\nu=b+3$, to ensure the existence of the first and second moments of $\bm{\Sigma}_{j}$ a-priori. Other specified hyperparameters include $\mu_{0,1}=\overline{x}_{1}$, $u_{1}=1/2$, $\alpha_{1}=2$, $\gamma_{1}=u_{1}(\text{range}(x_{1:n,1})/4)^{2}$, $\bm{\varrho}_{k}=(1,1)$ for $k=p+1,\ldots,q$, and the parameters of the stick-breaking prior are $\zeta_{j,1}=1$ and $\zeta_{j,2}=1$. Here $\overline{x}_{1}$ and $\text{range}(x_{1:n,1})$ denote the sample mean and range of $(x_{1,1},\dots,x_{n,1})$.

The MCMC stage of the adaptive truncation algorithm, with $J_{0}=15$ components, is run for 20,000 iterations after discarding the first 10,000 as burn-in. Every 10-th iteration is saved to produce $M$= 2,000 initial values for the particles in the SMC stage.

Results

In Appendix A, we describe various posterior and predictive quantities that can be computed from the weighted particles to describe the relationship between the observed response ${\bf z}$ and covariates ${\bf x}$. Here, we focus on the marginal predictive mean and density functions for (undiscretized) $Z_{1}$ and $Z_{2}$, as well as on the marginal predictive probability of success for $Z_{3}$, and compare them with the true data-generating functions in Figures B.1 and B.2, for a selected combinations of the categorical covariates. We also show the conditional mean of $Z_{2}$ as a function of $x_{1}$ given different values of $z_{1}$ in Figure B.3. Overall, the model is able to recover the underlying structure present in the data, despite the heavy censoring of $Z_{2}$ for lower levels of $x_{1}$, particularly when $x_{3}^{*}=2$. The true conditional structure is well recovered in areas where data (crosses) is available, i.e. in the left plot of Figure B.3, for $15\leq x_{1}\leq 18$ given $z_{1}=15$ (magenta) and for $27\leq x_{1}\leq 29$ given $z_{1}=23$ (black). However, as can be expected, the model struggles when predicting at values far from the observed data, i.e in the left plot of Figure B.3, for $25\leq x_{1}\leq 30$ given $z_{1}=15$ (magenta) and for $15\leq x_{1}\leq 20$ given $z_{1}=23$ (black). We highlight that interpretation of the conditional dependence structure in the latent scale as well as the latent covariance matrices of the mixture components and its relation to the dependence structure on the observed scale is an open and interesting direction of research, which would expand the work of García-Zattera et al. (2007) in the parametric setting.

Robustness analysis.

We perform a robustness analysis comparing several initialization specifications, namely by setting $J_{0}=2,3,5,10,15,20,30$. We also compare with a parametric version of the model that is similar in nature to the parametric model of Korsgaard et al. (2003), i.e. a multivariate Gaussian regression model with the link functions $h_{\ell}({\bf y},{\bf x})$ in (7) and a prior given by the base measure ${\bf P}_{0}$. For the sake of comparison, we use the Metropolis-within-Gibbs scheme for inference. In all scenarios, the adaptive MCMC algorithm is run for 30,000 iterations, discarding the first 10,000 as burn-in, and saving only every 10th iteration for a total of $M=$ 2,000 particles to be used in the SMC step. A summary of the analysis is reported in Table B.1. Besides the number of components inferred by the model ($J^{*}$), and the elapsed CPU time (in hours), the table reports the ESS for the log-likelihood in the MCMC stage ($\text{ESS}_{\text{MCMC}}$) computed with the mcmcse package in R (Flegal et al., 2017), and the $\text{ESS}_{J^{*}}$ of the final iteration of the SMC, that can be used to compare mixing. Results are reported for two different discrepancy measures used to define the stopping rule of SMC, namely the ESS and the CESS (Zhou et al., 2016). In order to assess the fit of the model, we compute for each $Z_{\ell}$ the log-pseudo marginal likelihood (LPML, Geisser and Eddy, 1979), and the percentage absolute errors with respect to the true mean and true density at a set of new test covariates, ${\bf x}_{i}^{*}$, for $i=1,\ldots,n^{*}$, denoted by $\text{ERR}_{\text{Mean}}$ and $\text{ERR}_{\text{Dens}}$:

[TABLE]

where for each response $\ell=1,\ldots,d$, $\mu_{\ell}^{t}({\bf x}_{i}^{*})$ and $\widehat{\mu}_{\ell}({\bf x}_{i}^{*})$ indicate the true and estimated mean functions, and $f^{t}(\cdot|{\bf x}_{i}^{*})$ and $\widehat{f}(\cdot|{\bf x}_{i}^{*})$ indicate the true and estimated densities. For each response, densities are evaluated on a grid of values, $z_{1,\ell}^{*},\ldots,z_{G,\ell}^{*}$, with grid size $\Delta$. The results show robustness with respect to the choice of the discrepancy measure.

We observe that for $J_{0}\geq 20$ only a moderate number of components are added, suggesting that a sufficient approximation is obtained with around 20 components. Recall that the SMC is run for at least $I=4$ cycles, i.e. at least four new components are added to the initial model. Therefore, if $J_{0}$ is large enough, we have $J^{*}=J_{0}+I$. Generally, the computational time increases with $J_{0}$, although this is not always the case, especially if $\text{ESS}_{J}$ becomes too low so that resampling and rejuvenation are required. Despite the increased number of parameters for large $J_{0}$, the mixing of the MCMC, reflected in the $\text{ESS}_{\text{MCMC}}$, does not deteriorate; however, note the improved mixing for the parametric model, which has the least number of parameters, due to the absence of the covariate-dependent weights. Focusing on the SMC, a larger $J_{0}$ generally results in less degeneracy of the particles, reflected in a higher $\text{ESS}_{J^{*}}$. Finally the LPML, measuring the goodness of fit of the model, increases with $J_{0}$, while the errors in predictive mean and density both decrease. This is particularly true for $Z_{1}$, the most nonlinear response, while there is little improvement in the binary response $Z_{3}$, which is indeed simulated from a linear probit model. Similar results are obtained when substituting the ESS with the CESS in the discrepancy measure of the SMC, confirming robustness to the choice of the stopping rule. To conclude, initializing the algorithm with a conservative number of components provides a good compromise between computational time, mixing, and accuracy.

Focusing on the model flexibility and its ability to recover the correct structure present in the data, we explore three additional scenarios: (i) longer number of iterations (25,000 burn-in; 10 thinning; $M=2,500$); (ii) doubled sample size ($n=1,400$); (iii) omitted censoring. In all cases, $J_{0}=15$. The results are reported in Table B.2, showing the same goodness-of-fit and error indicators used in Table B.1. As expected, mixing improves for scenario (i), due to the higher number of iterations, and modest improvements in the predictive power are observed. Increasing the sample size (scenario (ii)) without adjusting the algorithm settings (e.g., length of the chain, $J_{0}$), does not lead to consistent improvements. Note that this scenario yields approximately doubled LPML values, due to the definition of this quantity. The removal of censoring (scenario (iii)) yields faster and more precise computations. Estimates of the predictive densities for the undiscretized variable $Z_{1}$ as function of $x_{1}$ are reported in Figure B.4.

Appendix C Application: life patterns of Colombian women

In the following sections, we provide some further discussion on the DHS data characteristics and report some additional figures, complementing and enriching the results reported in the main paper. We begin by providing a map of the regions of Colombia used in the study in Figure C.5. Different territorial divisions of Colombia are used in different contexts. We follow Ojeda et al. (2011), considering six regions: Atlantica, Oriental, Central, Pacifica, Bogota, Territorios Nacionales. We point out that Bogota is in fact a part of the Oriental region but was treated separately because of its peculiar features in terms of social and economic development. Territorios Nacionales could be divided into smaller, more homogeneous regions (e.g., Orinoquia and Amazon), but this is the definition employed by the DHS for the 2010 survey.

C.1 The data: some considerations on sample weights

The data arise from a complex survey design, and thus have associated weights with a rather complicated structure due both to survey design and to additional post-stratification carried out to adjust for various factors (e.g., the total number of women interviewed in each municipality, non-response, etc.; see Ojeda et al., 2011, for details). This creates difficulties in further adjusting the weights based on the filtering considered (see Section 2 in the main text). Specifically, we do not know the population proportions for each filter within each municipality (e.g. of ethnicity, or of residents that have lived in the same region from at least the age of 6). Thus, the weights for the filtered data would need to be post-adjusted based on the observed sample proportions, leading to uncertainty in the weights. In addition, how this correction should relate to the post-stratification is unclear. We emphasize that even if we were able to obtain accurate weights for the filtered data, the use of survey weights in regression analyses is debated within the literature. From a frequentist perspective, Winship and Radbill (1994) point out that when ordinary least squares (OLS) and weighted OLS lead to similar estimates the former is preferred, because it leads to improved efficiency, power, coverage. In fact, the authors state that “where the sampling weights are solely a function of the independent variables” (as in our case), “unweighted OLS estimates are preferred because they are unbiased, consistent, and have smaller standard errors.” Moreover, the authors suggest that when OLS and weighted OLS estimates are substantially different one should investigate more flexible models, with e.g. non-linearity or interactions. Thus, since our model is nonparametric and can flexibly recover non-linearity, interactions, multi-modality (that may be due to missing or unmeasured covariates), we prefer the unweighted median regression and density regression estimates.

From a Bayesian perspective, although our model and algorithm could be adjusted to include the sample weights, for example through a pseudo (i.e. weighted) likelihood or via data augmentation, as outlined by Gunawan et al. (2020), this may lead to Bayesian interval estimates that do not have the correct frequentist coverage. Indeed, the authors show that this is the case with the pseudo likelihood approach, and provide empirical evidence that data augmentation can perform better. In summary, due to the issues involved in the determination of accurate weights for the filtered data and to the poor coverage and efficiency of the weighted estimates that has been identified in literature, we carry out inference based on the unweighted data, which is further supported by the flexibility of our model.

C.2 Additional figures

We display additional figures, enriching the results reported in the main text, and for convenience, comments on possibly relevant findings are reported in the figures’ captions. Figure C.6 complements Figure 5 by reporting median ages at events for women who grew up in violent environments with only physical punishment or only parental domestic violence, i.e (${\mathbf{P}},{\mathbf{\bar{B}}}$) or (${\mathbf{\bar{P}}},{\mathbf{B}}$). Figure C.7 completes Figure 6 by displaying the predictive density of the age at sexual debut. Figure C.8 reports the predictive probability of censoring, that is the probability that a woman will experience the event after the given Age, as a function of Age. Turning to the conditional analysis, the conditional predictive medians for the time from sexual debut to first child given the age at sexual debut is shown in Figure C.9 and for the time from union to first child given the age at union is shown in Figure C.10. The underlying conditional predictive densities for selected covariate combinations are visualized in Figures C.11, C.12, and C.13. Finally, to explore the possible relation between anticipation of union on work activity, Figure C.14 reports the conditional predictive probability of working as function of $Age$ given different ages at union.