\section{Introduction}

The regression tree is a nonparametric method for estimating a regression function. It has the advantages of being flexible, easily interpretable, well-suited to uncovering interactions among predictors, and having the ability to produce well-defined subgroups of the data via membership in the terminal nodes of the constructed tree. Tree methods have been developed primarily for situations with a univariate response variable \citep[see, e.g.,][]{breiman1984classification}, but they can be applied as well to multivariate response data. A multivariate regression tree provides two additional potential benefits to those for a univariate response: the ability to allow for correlation between the response variables, and the ability to impose the restriction of a single tree structure that is driven by the multiple responses simultaneously. Consider, for example, data with a bivariate blood pressure response (systolic and diastolic): a method using various predictors that determines a single set of subgroups of patients with ``similar blood pressure" based on both responses is likely to be much more useful than a method that produces two different sets of subgroups for the two responses. \cite{de2002multivariate} proposed such an algorithm called Multivariate Regression Trees (MRT) that can construct a single tree for data with multiple response variables, implemented in R package $\mathtt{mvpart}$. MRT straightforwardly extends the univariate regression tree by defining a new impurity measure for multiple response variables. It can then split the tree nodes in the same way as the univariate regression tree based on this new measure.

The traditional tree-based methods, including MRT, assume independence over the observations. This assumption is violated in longitudinal or clustered data, as longitudinal structure often implies correlation among the observations within an object. Ignoring the structure can result in overstating the correlation in the population and negatively affecting model selection, often leading to models that are too complex. Therefore, the tree-based algorithms have to be specifically modified to tackle this issue. \cite{segal1992tree} first extended the regression tree to a multiple response version that is designed for longitudinal data. However, the ``multiple response" in this method only refers to the outcomes of a single response variable at different time points. It considers a covariance structure among those outcomes instead of assuming them as independent, but does not allow for different response variables such as the blood pressure variables mentioned above. \cite{hajjemetal} proposed a mixed effects regression tree (MERT) algorithm for clustered data, which divides the model into two parts: a fixed (population-averaged) effect and a random (cluster-specific) component. They use a standard regression tree (CART) for the fixed effect and a linear structure for the random effects. \cite{sela2012re} independently proposed a similar estimation method called RE-EM that also uses CART to estimate the fixed effects and a linear model for the random effects. They alternate between estimating the regression tree, assuming that estimates of the random effects are correct, and estimating the random effects, assuming that the regression tree for the fixed effects is correct. Both MERT and RE-EM are designed for univariate response data, which motivates us to extend their approach to a multivariate version by replacing CART with MRT. 

In this paper we extend the RE-EM tree method of \cite{sela2012re} to multivariate response data. In Section \ref{sec:method} we outline the algorithm for producing the tree. Section \ref{sec:MCsimul} provides discussion of Monte Carlo simulations investigating the performance of the method. In Section \ref{sec:poverty} we examine a real data example that explores the multivariate nature of poverty, and Section \ref{sec:concl} concludes the paper.

\section{The multivariate RE-EM tree}
\label{sec:method}

To formally define the algorithm, we partially follow the notation in \cite{sela2012re}. Suppose there are $I$ objects $i=1,2,.., I$, and each object is observed at times $t=1,2,\cdots,T_i$. Let $\bm{Y}_{it}=\left(Y^{(1)}_{it},\cdots,Y^{(J)}_{it}\right)^{\top}$ denote the $J$ response variables of object $i$ at time $t$, and $\bm{X}_{it} = \left(X_{it1},\cdots,X_{itk}\right)^{\top}$ denote the $k$ corresponding covariates. Given a design matrix $\bm{Z}_{it} \in \mathbb{R}^{q\times 1}$ at each observation, the multivariate model can be formalized as 
\begin{equation}
	\begin{aligned}
		\bm{Y}_{it} &= \bm{f}(\bm{X}_{it}) + \bm{B}_{i}\bm{Z}_{it} +  \bm{\varepsilon}_{it},
	\end{aligned}
	\label{model}
\end{equation}
where $\bm{f}(\bm{X})$ is a multivariate regression tree,  $\bm{B}_{i}$ is the random effect matrix of object $i$, and $\bm{\varepsilon}_{it}$ is the noise of observation $(i, t)$. Concretely, if we let $g_1,\cdots,g_L$ denote the $L$ terminal nodes (i.e. leaves) of $\bm{f}$, and $\mu^{(j)}_{1},\cdots,\mu^{(j)}_{L}$ denote the predicted response of $Y^{(j)}$ on each leaf, we can rewrite $\bm{f}(\bm{X})$ as  
\begin{equation}
	\bm{f}(\bm{X})= \left(\sum_{l=1}^L \mathbb{I} (\bm{X} \in g_l)\mu^{(1)}_{l},\, \cdots \, ,\sum_{l=1}^L \mathbb{I} (\bm{X} \in g_l)\mu^{(J)}_{l}\right)^{\top},
	\label{model_tree}
\end{equation}
where $\mathbb{I}(\cdot)$ is the indicator function. $\bm{B}_{i}=\left(\bm{b}^{(1)}_{i},\cdots,\bm{b}^{(J)}_{i}\right)^{\top} \in \mathbb{R}^{J\times q}$, where $\bm{b}^{(j)}_{i}$ stands for the time-constant, object-specific random effects of object $i$ on the $j$-th response variable. In particular, if we set $\bm{Z_{it}}=1$ and $\bm{B}_{i}=\left(b^{(1)}_{i},\cdots,b^{(J)}_{i}\right)^{\top} \in \mathbb{R}^{J\times 1}$, only a random intercept for each object for each response variable will be included in the model. We assume \[\mathrm{vec}(\bm{B}_i)=\begin{pmatrix}\bm{b}^{(1)}_{i}\\ \vdots \\ \bm{b}^{(J)}_{i}\end{pmatrix} \sim \mathcal{N} \left(\bm{0},\bm{D}\right),\] 
which allows correlation between $\bm{b}^{(j)}_{i}$s. The noise $\bm{\varepsilon}_{it}=\left(\varepsilon^{(1)}_{it},\cdots,\varepsilon^{(J)}_{it}\right)^{\top}$, and we suppose for each $i$,
\[ \begin{pmatrix}\bm{\varepsilon}_{i1}\\
	\bm{\varepsilon}_{i2}\\
	\vdots\\
	\bm{\varepsilon}_{iT_i}\end{pmatrix} \sim \mathcal{N} (\bm{0},\bm{R}_i).\]
In other words, we assume the noise random variables are independent across different objects, but within one object they are allowed to be correlated by setting non-diagonal $\bm{R}_i$.

We are now ready to define the multivariate RE-EM algorithm:

\begin{enumerate}
	\item[\textbf{Step 1}] Initialize $\widehat{\bm{B}}_i=\left(\widehat{\bm{b}_i^{(1)}}, \cdots,\widehat{\bm{b}_i^{(J)}}\right)=\bm{0}$ for each $i$.
	
	
	\item[\textbf{Step 2}] Iterate the following (a) and (b) until convergence (the increase of the likelihood function in step (b) being less than some tolerance
	value).
	
	(a) Fit a Multivariate Regression Tree with target variables $\left\lbrace\bm{Y}_{it}-\widehat{\bm{B}}_i\bm{Z}_{it}\right\rbrace$ 
	and predictors $\left\lbrace\bm{X}_{it}\right\rbrace$, $i\in\{1,2,...,I\}, t\in\{1,..., T_i\}$, using the method proposed in \cite{de2002multivariate}. Record the terminal nodes of this tree, denoted by $\widehat{g}_1,\cdots,\widehat{g}_L$.
	
	(b) Fit a multivariate mixed linear model as follows,
	\[ \bm{Y}_{it}= \begin{pmatrix}
		\mu_1^{(1)}&\mu_2^{(1)}&
		\cdots&\mu_L^{(1)}\\
		\vdots & \vdots  & \ddots   & \vdots \\
		\mu_1^{(J)}&\mu_2^{(J)}&\cdots&\mu_L^{(J)}
	\end{pmatrix}\begin{pmatrix}
		\mathbb{I} (\bm{X}_{it} \in \widehat{g}_1)\\
		\vdots\\ 
		\mathbb{I} (\bm{X}_{it} \in \widehat{g}_L)
	\end{pmatrix} + \bm{B}_i\bm{Z}_{it}  + 
	\bm{\varepsilon}_{it}.\]
	Extract the estimated random effects $\widehat{\bm{B}}_{i}$ from this model.
	
	\item[\textbf{Step 3}] Take the mixed linear model in the last iteration to be the final result, with parameters $\widehat{g}_l, \hat{\mu}^{(j)}_{l}, \widehat{\bm{B}}_{i}$.
	
\end{enumerate}

The Multivariate Regression Tree proposed in \cite{de2002multivariate} uses Euclidean distance for measuring the diversity within a tree node. Since the response variables may have different scales, we should potentially standardize them at the beginning of the algorithm. For pruning, the RE-EM tree uses $k$-fold cross-validation in each iteration to select the tree, using the selecting method ``minimum" or ``one-standard-error". More details about standardization and pruning will be given in Section \ref{sec:methods}. 



\section{Simulation results}
\label{sec:MCsimul}

\subsection{Bivariate RE-EM with simple tree structure}\label{sec:bi-simple}

We first set $\bm{f}(\bm{X})$ to be a relatively simple structure and $J=2$ to compare the multivariate RE-EM tree to other competitor methods listed in section \ref{sec:methods}. 

\subsubsection{Settings} \label{sec:setting}

We generate simulated data sets with $I=50, 100, 200, 400, 800$ and $T_i=5, 10, 25, 50$ for all $i$. Following the procedure in \cite{fu2015unbiased}, we first generate $X_1,X_2,X_3,X_4,X_5,X_6 \sim Unif(0,10)$ (i.i.d), then  round $X_6$ to the nearest integer. Additionally, we generate another predictor $X_7 \sim Unif\{0,5.77\}$, which has the same variance as $X_1,..,X_5$. The tree structure is set as follows (using the notation defined in \eqref{model}, also shown in Figure \ref{fig:simpletree}): 
\[ g_1=(X_1 \leq 5) \wedge (X_2 \leq 5), \mu^{(1)}_{1} = 10, \mu^{(2)}_{1} = 6.\]
\[ g_2=(X_1 \leq 5) \wedge (X_2 > 5), \mu^{(1)}_{2} = 11, \mu^{(2)}_{2} =7.\]
\[ g_3=(X_1 > 5) \wedge (X_3 \leq 5), \mu^{(1)}_{3} = 12, \mu^{(2)}_{3} = 8.\]
\[ g_4=(X_1 > 5) \wedge (X_3 > 5), \mu^{(1)}_{4} = 13, \mu^{(2)}_{4} = 9.\]

\begin{figure} 
	\centering   
	\begin{forest} 
		[$X_1<5$, circle,draw
		[$X_2<5$, circle, draw
		[{$\mu^{(1)}_1=10$}\\ {$\mu^{(2)}_1=6$}, align=center, base=bottom]
		[{$\mu^{(1)}_2=11$}\\ {$\mu^{(1)}_2=7$}, align=center, base=bottom]
		]
		[$X_3<5$, circle, draw
		[{$\mu^{(1)}_3=12$}\\ {$\mu^{(2)}_3=8$}, align=center, base=bottom]
		[{$\mu^{(1)}_4=13$}\\ {$\mu^{(2)}_4=9$}, align=center, base=bottom]
		]
		]
	\end{forest}
	\caption{Simple tree structure, as discussed in section \ref{sec:setting}}\label{fig:simpletree}
\end{figure}
For the random effect, we take the design matrix $\bm Z_{it}=1$, with
$\begin{pmatrix}\bm{b}^{(1)}_{i}\\ \bm{b}^{(2)}_{i}\end{pmatrix} \sim \mathcal{N} \left(\bm{0},\bm{D}\right)$ and $\bm{D}=\begin{pmatrix}
	1& \sigma_{12}\\\sigma_{12} & 1
\end{pmatrix}$, where  $\sigma_{12}=0$, $.25$, $.5$, $.75$. The covariance matrix of noise $\bm{R}_i$ is set to be $ \begin{pmatrix}
	\sigma_\varepsilon& 0\\0 & \sigma_\varepsilon
\end{pmatrix}$ with $\sigma_\varepsilon=.5$, $1$, $1.5$, $2$. Based on this setting, the response values of $\bm{Y}$ are generated by \eqref{model} and \eqref{model_tree}. 

For each $(I, \sigma_{12}, \sigma_\varepsilon)$, we generate a test set containing the same $I$ objects as in the training set, each with 20 new observations. Concretely, we keep $\bm{b}^{(1)}_{i}, \bm{b}^{(2)}_{i}$ to be the same as the training set for all $i$, then randomly sample $\bm{X}_{it_{\mathrm{test}}}$ and $\bm{\varepsilon}_{it_{\mathrm{test}}}$ for $t_{\mathrm{test}}=1,2,...,20$ from the same distribution as above. 



\subsubsection{Methods for Comparison}\label{sec:methods}

For comparison, we apply the following 4 competitive methods on the simulation sets:

\begin{enumerate}
	\item[1.] Separate univariate linear mixed models, referred to as ``unilme". This method only fits a linear mixed model of each response on the covariates, and does not account for correlation within $\bm{B}_i$.
	
	\item[2.] One multivariate linear mixed model, referred to as ``multilme". This method still imposes linearity, but allows correlation between the responses.
	
	\item[3.] One multivariate tree that doesn't account for the longitudinal structure, referred to as ``multitree".
	
	
	\item[4.] Separate univariate RE-EM trees, referred to as ``uniREEM". This is similar to method 1, but removes the linear model restriction. Note that this doesn't impose a restriction that there is a single tree structure driving the response variables simultaneously, which potentially adds noise if there actually is a single tree structure driving both responses.
\end{enumerate}

For the multivariate RE-EM tree algorithm, we also compare four versions, with two standardization methods and two criteria for selecting the tree after the $k$-fold cross-validation (we choose $k=10$ as the default value). The standardization methods are 
\begin{enumerate}
	\item ``marg": Marginalizing each response variable $Y^{(j)}$ by 
	\[Y_{it, \mathrm{scaled}}^{(j)}\leftarrow \frac{Y_{it}^{(j)}-\mathrm{mean}_{i,t}\left(Y_{it}^{(j)}\right)}{\mathrm{sd}_{i,t}\left(Y_{it}^{(j)}\right)}.\]
	\item ``cov": Estimate the covariance matrix $\mathrm{cov}(\bm Y)$ at the beginning, where $\bm Y=\left(Y^{(1)},...,Y^{(J)}\right)$. Let
	\[\bm Y_{\mathrm{scaled}} \leftarrow \bm Y \left(\mathrm{cov}(\bm Y)\right)^{-1/2}.\]
	In this way, the Euclidean distance computed by the Multivariate Regression Tree will be actually the Mahalanobis distance.
\end{enumerate}

And the choosing criteria are
\begin{enumerate}
	
	\item ``min": Choosing the tree with minimum estimated predictive error based on $k$-fold
	cross-validation. 
	
	\item ``1se": Choosing the simplest tree with cross-validated estimated predictive error within one standard error of the minimized estimated predictive error.
	
\end{enumerate}

The main measure to evaluate the performance of the models is the object-level Predicted Mean Squared Error (PMSE), which is defined by
\[ \widehat{PMSE} = \frac{1}{20I} \sum_{i, t_{\mathrm{test}},j}\left(\widehat{Y}^{(j)}_{it_{\mathrm{test}}}-Y^{(j)}_{it_{\mathrm{test}}}\right)^2.\]
We also report the EMSE (Estimation Mean Squared Error) of the population-level fixed effect for all the methods. Specifically, for all methods except ``multi-tree", given a fitted linear mixed model $\bm{Y}_{it}=\widehat{\bm{f}}(\bm{X}_{it})+\widehat{\bm{B}}_i\bm{Z}_{it}$ and a new data point $(\bm{X}_0,\bm{Z}_{i_0t_0})$ $(i_0 \in \{1,2,..,I\})$, the estimated fixed effect at the new point is $\widehat{\bm{f}}(\bm{X}_{0})$ and the predicted value $\widehat{\bm Y}_{it_{\mathrm{test}}} = \widehat{\bm{f}}(\bm{X}_0) + \widehat{\bm{B}}_{i_0}\bm{Z}_{i_0t_0}$. For ``multi-tree", the estimated fixed effect and the predicted value are the same since there is no random effect. Further, we also record the tree-recovering rates of the tree-based methods, i.e., the rate that a tree method can capture the tree structure (the shape of the tree and the variable on which each split is based) correctly. In the ``uni-REEM" method, we record the recovering rates of the $J$ univariate RE-EM trees separately.  The final results are averaged over 400 repeats.



\subsubsection{Results for simpler bivariate tree structure}

The results for the estimation in population and prediction in object for all of the methods when $T_i=5$ are shown in Figures \ref{fig:ti5-str1-fixed} and \ref{fig:ti5-str1-group-all}, respectively. In both figures, it is obvious that the tree-based methods performed better than the linear methods. Figure \ref{fig:ti5-str1-fixed} indicates that all of the tree methods are competitive for population-level estimation. This illustrates that the dominant factor in estimating the population-level expected response is to recognize it as a tree; beyond that, doing so with a tree that accounts for the longitudinal structure provides an advantage over one that does not account for it, but this effect is virtually gone when the number of objects $I$ exceeds 200 to 400. For this reason, we will focus on object-level prediction here. 

For the object-level prediction in Figure \ref{fig:ti5-str1-group-all}, both of the versions of the proposed algorithm ``multiREEM" with ``marg" standardization outperform the other methods, including ``multitree" and ``uniREEM'', especially when the correlation $\sigma_{12}$ or the noise $\sigma_\varepsilon$ is large. As would be expected, the trees that account for the random effects clearly outperform the tree that ignores the longitudinal structure when attempting prediction at the level of the individual, so the latter will not be considered further.

\newpage

\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/all/Ti=5-str=1-fixed}}
	\caption{EMSE in population for $T_i=5$ for the simpler bivariate tree. Note that the lines for separate univariate linear mixed effects models (``unilme") and a single multivariate linear mixed effects model (``multilme") are almost the same.}
	\label{fig:ti5-str1-fixed}
\end{figure}

\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/all/Ti=5-str=1-group}}
	\caption{PMSE in object for $T_i=5$ for the simpler bivariate tree structure.}
	\label{fig:ti5-str1-group-all}
\end{figure}

To see the differences among the tree-based methods more clearly, we discarded linear methods ``unilme" and ``multilme" in the following figures. The object-level PMSE for $T_i=5, 10, 25, 50$ are shown in Figures \ref{fig:ti5-str1-group-part1}--\ref{fig:ti50-str1-group-part1}. In almost all of the cases, the multiREEM tree ``min\_marg" and ``1se\_marg" perform the best, with no meaningful difference between them. The other standardizing method ``cov" does not perform very well, and ``1se\_cov" is always worse than ``min\_cov". For $T_i= 5, 10, 25$, the multiREEMtree ``1se\_cov" is never better than the ``uniREEM" either for small $I$, small $\sigma_{12}$ or small $\sigma_\varepsilon$. Only when $I, \sigma_{12}, \sigma_\varepsilon$ are all large enough, or $T_i=50$, can the ``cov" standardization show advantage over ``uniREEM", but it is still worse than the ``marg" standardization. 

The better performance of the ``min\_cov" pruning versus the ``1se\_cov" pruning suggests that the covariance standardization tends to result in fitted trees that are too simple. This is further supported by simulations based on a more complex fixed effects structure, in which the ``cov" standardization performs even more poorly; for this reason, it will not be considered further in the results presented here.

Using the ``min" criterion, the advantage of fitting a multivariate tree versus two separate univariate trees when the true underlying structure is a single tree is evident from the figures, as the univariate RE-EM trees always lag behind in performance, sometime by a considerable amount. Clearly, avoiding the additional noise of fitting two separate longitudinal trees provides the potential of meaningful gains in predictive accuracy.

As would be expected, predictive performance improves with increasing number of objects $I$, and increasing number of time points $T_i$, reflecting both better estimation of the fixed effects and better prediction of the random effects.

The results are fairly insensitive to the presence of correlation in the random effects. Separate univariate RE-EM trees do not use that information, of course, and their performance suffers accordingly. Interestingly, the ``cov" standardization is also negatively affected by correlation, presumably because correlation of random effects is not the same as correlation between responses. 

\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/part1/Ti=5-str=1-group}}
	\caption{PMSE in object for $T_i=5$ for the simpler bivariate tree structure.}
	\label{fig:ti5-str1-group-part1}
\end{figure}
\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/part1/Ti=10-str=1-group}}
	\caption{PMSE in object for $T_i=10$ for the simpler bivariate tree structure.}
	\label{fig:ti10-str1-group-part1}
\end{figure}
\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/part1/Ti=25-str=1-group}}
	\caption{PMSE in object for $T_i=25$ for the simpler bivariate tree structure.}
	\label{fig:ti25-str1-group-part1}
\end{figure}
\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/part1/Ti=50-str=1-group}}
	\caption{PMSE in object for $T_i=50$ for the simpler bivariate tree structure.}
	\label{fig:ti50-str1-group-part1}
\end{figure}

The tree recovering rates for the Multivariate Regression Tree, Univariate RE-EM trees and Multivariate RE-EM tree with marginal normalization are shown in Figures \ref{fig:ti5-str1-rec}--\ref{fig:ti50-str1-rec}. When $I$ is large enough, the multivariate RE-EM tree methods with marginal normalization and ``1se" criterion always have the highest recovering rate, while the ``min" criterion does not achieve this even when $I=800$. This implies that choosing the tree based on the minimum CV error may lead to overfitting, although this does not impact the Prediction Error, and the ``one-standard-error" choosing strategy overcomes this difficulty. However, when the noise is large and the sample size $I$ is small, the ``min" criterion achieves the best recovering rate compared to other competitors. It is striking that the univariate trees achieve perfect recovery of the underlying tree structure for a large enough number of objects and/or time points, but still typically lag behind the multivariate trees in terms of PMSE. This is because of both poorer estimation of fixed effects and poorer prediction of random effects; see Figures \ref{fig:SM:ti5-str1-fixed}--\ref{fig:SM:ti50-str1-RE} in Section \ref{sec:SM:RE} of the Appendix. 

\newpage

\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/part2/Ti=5-str=1-rec}}
	\caption{Recovering rate for $T_i=5$ for the simpler bivariate tree structure.}
	\label{fig:ti5-str1-rec}
\end{figure}
\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/part2/Ti=10-str=1-rec}}
	\caption{Recovering rate for $T_i=10$ for the simpler bivariate tree structure.}
	\label{fig:ti10-str1-rec}
\end{figure}
\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/part2/Ti=25-str=1-rec}}
	\caption{Recovering rate for $T_i=25$ for the simpler bivariate tree structure.}
	\label{fig:ti25-str1-rec}
\end{figure}
\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/part2/Ti=50-str=1-rec}}
	\caption{Recovering rate for $T_i=50$ for the simpler bivariate tree structure.}
	\label{fig:ti50-str1-rec}
\end{figure}


\subsection{A more complex tree structure}\label{sec:bi-complex}

In this section, we consider a more complex tree structure as shown in Figure \ref{fig:complextree}, with seven terminal nodes. As compared to the simpler tree in Figure \ref{fig:simpletree}, this tree is deeper and has an asymmetrical structure, which makes it more difficult to estimate the fixed effects and recover the tree structure.  

\newpage

\begin{figure}[H] 
	\centering   
	\begin{forest} 
		[$X_1<5$, circle,draw
		[$X_2<5$, circle, draw
		[$X_4<5$, circle, draw
		[{$\mu^{(1)}_1=6$}\\ {$\mu^{(2)}_1=4.5$}, align=center, base=bottom]
		[{$\mu^{(1)}_2=8$}\\ {$\mu^{(1)}_2=6.5$}, align=center, base=bottom]
		]
		[$X_5<5$, circle, draw
		[{$\mu^{(1)}_3=10$}\\ {$\mu^{(2)}_3=8.5$}, align=center, base=bottom]
		[{$\mu^{(1)}_4=12$}\\ {$\mu^{(1)}_4=10.5$}, align=center, base=bottom]
		]
		]
		[$X_3<5$, circle, draw
		[{$\mu^{(1)}_5=14$}\\ {$\mu^{(2)}_5=10.5$}, align=center, base=bottom]
		[$X_6<5$, circle, draw
		[{$\mu^{(1)}_6=16$}\\ {$\mu^{(2)}_6=12.5$}, align=center, base=bottom]
		[{$\mu^{(1)}_7=18$}\\ {$\mu^{(2)}_7=14.5$}, align=center, base=bottom]
		]
		]
		]
	\end{forest}
	\caption{A deeper and asymmetrical tree structure for section \ref{sec:bi-complex}. The true value of the $j$-th response variable on the $l$-th terminal node are given by $\mu_l^{(j)}$.}\label{fig:complextree}
\end{figure}

We present the object-level PMSE and the tree recovering rate for $T_i=5$ in Figure \ref{fig:ti5-str2-group-part1} and Figure \ref{fig:ti5-str2-rec}, respectively. The results for $T_i>5$ are shown in Section \ref{sec:SM:complex-bivariate} of the Appendix. The patterns are similar to those seen for the simpler model, with the multivariate RE-EM tree superior in PMSE to two separate RE-EM trees. It is also clear that the simpler trees that can come from using the ``1se" rule are sometimes too simple (that is, the final trees are underfit), and result in poorer performance than using the ``min" rule, particularly for smaller sample sizes. The ``min" rule can lead to overfitting, and hence a lower recovering rate, but this clearly has relatively little effect on the PMSE compared to that of underfitting.

\newpage

\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/part2/Ti=5-str=2-group}}
	\caption{PMSE in object for $T_i=5$ for the complex bivariate tree}
	\label{fig:ti5-str2-group-part1}
\end{figure}



\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/part2/Ti=5-str=2-rec}}
	\caption{Recovering rate for $T_i=5$ for the complex bivariate tree}
	\label{fig:ti5-str2-rec}
\end{figure}


\subsection{More Response Variables}

\label{sec:multi-simple}

In this section, we describe simulations with $J=5$ to evaluate the method with more response variables. The settings are similar to those in Section \ref{sec:setting}, only with the following tree structure: 
\[ g_1=(X_1 \leq 5) \wedge (X_2 \leq 5), \bm{\mu}_1 = (10, 9, 8, 4, 6)^{\top}.\]
\[ g_2=(X_1 \leq 5) \wedge (X_2 > 5), \bm{\mu}_2 = (11, 10, 9, 5, 7)^{\top}.\]
\[ g_3=(X_1 > 5) \wedge (X_3 \leq 5), \bm{\mu}_3 = (12, 11, 10, 6, 8)^{\top}.\]
\[ g_4=(X_1 > 5) \wedge (X_3 > 5), \bm{\mu}_4 = (13, 12, 11, 7, 9)^{\top}.\]
where $\bm{\mu}_l=\left({\mu}^{(1)}_l,\cdots,{\mu}^{(5)}_l\right)^{\top}$. For the random effect, we again take the design matrix $\bm{Z}_{it}=1$ and adjust
the covariance matrix $\bm{D}=(\sigma_{j_1j_2})_{5\times 5}$ as: $\sigma_{j_1j_2}=\sigma_B$ if $j_1\neq j_2$ and $\sigma_{j_1j_2}=1$ if $j_1=j_2$. In particular, we take $\sigma_{B}=0, .25, .5, .75$. The object-level PMSE is shown in Figure \ref{fig:ti5-str3-group-part1}, the tree recovering rate for $T_i=5$ is shown in Figure \ref{fig:ti5-str3-rec}, and results for $T_i>5$ are given in Section \ref{sec:SM:multivariate} of the Appendix. 





The results of the multivariate RE-EM tree with five response variables are similar to those of the bivariate tree. Moreover, with more response variables, the advantage of the multivariate RE-EM tree methods over ``uniREEM" is more pronounced, especially, the correlation is more apparent for higher correlation. As we hypothesized earlier, fitting a single multivariate tree apparently avoids the additional noise of fitting separate univariate longitudinal trees, and accounting for the correlation between the random effects further improves prediction performance. The tree recovering rates again verify our findings in the bivariate-tree simulations. When the sample size is small, the Multivariate RE-EM tree methods always achieve significantly higher recovering accuracy than each separate univariate RE-EM tree, and the ``1se" criterion performs worse in tree recovering compared to the ``min" criterion, which leads to a poorer prediction. However, when the sample size is large enough, the ``min" criterion may overfit the tree structure, but this does not affect the prediction. 


\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/part2/Ti=5-str=3-group}}
	\caption{PMSE in object for $T_i=5$ for multivariate tree with 5 response variables}
	\label{fig:ti5-str3-group-part1}
\end{figure}

\begin{figure}[H]
	\centerline{
		\includegraphics[width=1.1\linewidth, height=0.5\textheight]{figs/part2/Ti=5-str=3-rec}}
	\caption{Recovering rate for $T_i=5$ multivariate tree with five response variables}
	\label{fig:ti5-str3-rec}
\end{figure}



\section{Real data example}
\label{sec:poverty}

Poverty is typically measured based on income, focusing on the expenditure necessary to buy a basket of essential goods and services. Such a monetary measure does not necessarily reflect other, non-financial, aspects of well-being. The World Bank, in its 2021 document {\it A Roadmap for Countries Measuring Multidimensional Poverty\/} \citep{WorldBank} argues for the concept of multidimensional poverty, which reflects that human deprivation encompasses characteristics beyond income and consumption expenditures. Such characteristics include access to services (electricity and clean water), health care, and education. 

We examine modeling of this multidimensional character of poverty based on data from 107 countries using annual data from 2012 through 2017. We measure poverty using four response variables: percentage of population with access to electricity, percentage of population using at least basic drinking water services, child (under 5) survival rate per 1000 live births, and gross percentage secondary school enrollment, respectively. We use as potential predictors national-level variables related to factors that have been proposed as potential causes of poverty (see, for example, \url{https://en.wikipedia.org/wiki/Causes\_of\_poverty}). Among these causes of global poverty, and the variables used to measure them, are inadequate food (agriculture production index), inadequate access to health care (health expenditure per capita), inequality (Gini index), poor education (per capita education expenditures), climate change (mean surface temperature change), lack of access to livelihoods or jobs (unemployment rate), poor governance and corruption (Democracy Index (DI), Corruption Perceptions Index (CPI)), limited or poor infrastructure (rural population, trade and transport-related infrastructure index, population density), conflict (violent conflict death rate), and overpopulation (population density). Data for all of the variables are available at \cite{WorldBankdata}, with the exception of the agriculture production index \citep{agprodindex}, mean surface temperature change \citep{tempchange}, and the DI \citep{DemocIndex}. Figure \ref{fig:tree} gives the multivariate RE-EM tree for these data, with pruning based on the one-SE rule and 10-fold cross-validation, and surrogate split used to handle missing values in the predictors.

\begin{figure}[H]
	\centerline{\includegraphics[width=1.3\textwidth]{fignogdptree}}
	\caption{Multivariate RE-EM tree of poverty data using all potential predictors. Each bar represents the average value of all country-years in that terminal node for each of the four response variables (Access to electricity, Access to safe water, Child Survival rate, and Secondary school enrollment).}
	\label{fig:tree}
\end{figure}

Each terminal node gives a bar plot of the average value of each response variable for all country-year pairs that fall in that node, with the bars presented in the order $\{$Access to electricity, Access to safe drinking water, Child survival rate, Secondary school enrollment$\}$. Thus, a shorter bar corresponds to a higher degree of that dimension of poverty. It should be noted, however, that those bars are in general less informative than they might appear to be, at least at the level of an individual country, since they do not reflect the predicted random effects of each response variable for that country, so we focus instead on the clustering implied by the tree: which predictors matter, in what way they matter, and how that translates to meaningful groupings of countries determined by the terminal nodes. The resultant groupings are given in the choropleth (heat) maps presented in Figure \ref{fig:mapnodes}, in which a darker color of a country corresponds to more years in which the country falls in that terminal node. 

\begin{figure}[H]
	\centerline{\includegraphics[width=1.25\textwidth]{figmapnodesnogdp}}
	\caption{Choropleth map of the number of years each country falls into each terminal node for tree using all potential predictors. Colors range from white to red as the count increases, with countries colored gray having no observations in the data set.}
	\label{fig:mapnodes}
\end{figure}


The groupings produced by the tree imply familiar geographic and socioeconomic patterns, with the tree providing nonfinancial characterizations of those patterns. The country-years with higher levels of nonfinancial poverty are separated from those with lower poverty on the basis of rural population, supporting the potential importance of infrastructure. The more urban (lower poverty) country-years are then split by CPI value. The less corrupt ones form Node 6, dominated by the wealthiest countries in the sample (Canada, western Europe, Israel, Australia, and New Zealand), but also the wealthiest countries in their geographic regions (Chile, Costa Rica, and Uruguay in Central and South America, and Slovenia, Poland, and the Baltic countries in the former Communist countries of Europe), pointing to a local relative wealth effect, rather than only an absolute one. This grouping strongly parallels an ordering based on the Human Development Index (HDI) \citep{HDI}, given that all of these countries other than Costa Rica and Uruguay rank among the top 45 countries in the world in HDI. The country-years with lower rural population and higher measured corruption form Node 5. This contains the rest of South and Central America, as well as southern Europe, the Middle East, and central Asia. These geographic groupings, perhaps pointing to cultural patterns, seem more important here, as these countries range from relatively developed Italy and Greece to less developed Honduras and El Salvador, but overall this node contains a ``second tier" level of human development.

The other branch of the tree (with rural population percentage greater than 47\%) defines four groups of country-years with higher poverty rates. The first split is on CPI, which separates the six years of data for less corrupt Bhutan into Node 4, a country with notably higher access to electricity and secondary school enrollment (and thus lower nonfinancial poverty) compared to the geographically close Bangladesh, India, and Nepal (this corresponds to a distinction not revealed by the HDI, as the four countries have similar values of that measure). The corresponding map for Node 4 is not included in Figure \ref{fig:mapnodes}, since it corresponds to only one country. The remaining country-years are next split on DI, and then population density. Node 2, with lower DI and lower population density, are the poorest countries in the sample, and among the poorest in the world, concentrating in an east-west line along the southern edge of the Sahara, southeast Africa (Mozambique and Madagascar), and Afghanistan. Node 3 splits from these countries by having higher population density, identifying slightly wealthier countries in Africa and countries in south Asia (Bangladesh, Nepal, and Pakistan). The country-years with higher DI form Node 1, and represent a middle ground of wealth, including countries such as India, Indonesia, Mauritius, and the Philippines. Thus, the multivariate tree unsurprisingly identifies many patterns related to financial poverty, but also highlights patterns related to human development, as well as geographic (and presumably cultural) effects in nonfinancial poverty.


\section{Conclusion}
\label{sec:concl}

In this paper we have described the multivariate RE-EM tree, a generalization of the univariate RE-EM tree to multivariate response data. Longitudinal trees for univariate response data have been generalized in various ways, and potential future work would be to explore such generalizations for multivariate responses. Possibilities include responses from other distributions than the normal \citep{hajjem2017generalized}, random forests \citep{hajjem2014mixed, capitaine2021random}, and more general model fits at terminal nodes than a single location summary statistic for responses \citep{fokkema2018detecting, stegmann2018recursive, wei2020precision}.

An R package to implement the multivariate RE-EM tree, {\tt multiREEMtree}, is available at \url{https://github.com/WenboJing1998/Multivariate-RE-EM-tree}. The code to implement the simulations and data analysis can also be found there. The simulations were run on R 4.0.3 using the NYU High Performance Computing cluster (NYU HPC), and the data analysis was run on R 4.1.1 under Windows 10 Pro.  



