\section{INTRODUCTION}\label{sec:intro}

The influence of agent heterogeneity on cooperation in social learning has been a recent focus of research in many fields, including ecology, sociology, and decision theory \cite{li2019influence}. Studies on evolutionary human behavior provide evidence that individual differences can be leveraged to enhance collective prosperity \cite{amaral2016evolutionary}. Motivated by applications such as  social foraging and multi-robot coordination tasks, we study and design cooperative strategies for a group of agents making sequential explore-exploit decisions in an uncertain environment. The strategies we design incorporate agent heterogeneity to optimize the performance of the group through collective learning.

Consider a group of agents, each making a sequence of choices among options in an uncertain environment in order to maximize collective payoff. At each time step in the sequence, each agent chooses an option depending on the knowledge it has acquired about the environment up to that time step. Maximizing payoff  necessitates striking a balance between making choices that yield high immediate payoff, i.e., exploiting, and making choices that yield high information content and possibly high future payoffs, i.e., exploring. When an agent fails to acquire sufficient information about the environment to make optimal decisions, it must sacrifice exploitation potential in order to explore. However, in the group setting, agents can recover exploitation potential by gaining information through cooperation i.e., through collective learning.

Sequential decision making in uncertain environments that requires trading off exploitation and exploration is modeled mathematically by the bandit framework \cite{robbins1952some}. In the multi-armed bandit (MAB) problem, an agent is repeatedly faced with the task of choosing an option from a given set of options. At each time step the agent receives a stochastic reward drawn from a fixed probability distribution associated with the chosen option. The agent's goal is to maximize the  cumulative reward by the end of the decision-making process. This requires choosing frequently enough the optimal option i.e., the option with highest expected reward. In order to meet this requirement, the agent must simultaneously choose options that are known to provide high rewards (exploit) and choose lesser known options (explore) that might potentially provide even higher rewards \cite{lai1985asymptotically, auer2002finite}.

Maximizing   cumulative reward is equivalent to minimizing   cumulative regret, defined as the loss incurred by an agent  choosing a sub-optimal option instead of the optimal option. Since the probability distribution associated with each option is fixed, cumulative regret can be minimized by reducing the number of times sub-optimal options are chosen. Performance of the proposed algorithms for this problem is measured using expected cumulative regret. The paper \cite{lai1985asymptotically} establishes that any efficient policy chooses subotimal options asymptotically logarithmically in time. The paper \cite{auer2002finite} proposes an Upper Confidence Bound (UCB) based sampling rule that achieves a logarithmic expected cumulative regret uniformly in time.

The papers \cite{landgren2016distributed,landgren2016distributedCDC,landgren2020distributed,martinez2019decentralized,szorenyi2013gossip,madhushani2020dynamic,madhushani2020distributed,chakraborty2017coordinated,madhushani2019heterogeneous,landgren2018social,wang2020optimal} extend to the multi-agent setting and capture different aspects of collective learning.
In \cite{landgren2016distributed,landgren2016distributedCDC,landgren2020distributed,martinez2019decentralized,szorenyi2013gossip}, agents share their estimates of the expected reward of options with neighbors according to fixed communication structures.  The papers \cite{landgren2016distributed,landgren2016distributedCDC,landgren2020distributed}  use a running consensus algorithm to update estimates and provide graph-structure-dependent performance measures that predict the relative performance of agents and networks. The paper \cite{landgren2020distributed} also addresses the case of a constrained reward model in which agents that choose the same option at the same time step receive no reward.  The paper \cite{martinez2019decentralized}
proposes an accelerated consensus procedure assuming agents know the spectral gap of the communication graph and designs a decentralized UCB algorithm based on delayed rewards. The paper \cite{szorenyi2013gossip} considers a P2P communication where an agent is only allowed to communicate with two other agents at each time step.

In \cite{madhushani2020dynamic,madhushani2020distributed,chakraborty2017coordinated},  agents  share reward values and choices with neighbors according to stochastic communication structures that depend on the decision-making process. In \cite{madhushani2020dynamic}, each agent observes rewards and actions of its neighbors when it is exploring. In \cite{madhushani2020distributed}, each agent instead broadcasts its rewards and actions to its neighbors when it is exploring. In \cite{chakraborty2017coordinated}, at each time step, agents decide either to sample an option or to broadcast the last obtained reward to the entire group.

The setup in our earlier paper \cite{madhushani2019heterogeneous} is closest to that in the present paper: agents observe reward values and actions of their neighbors defined by a network graph that changes in time according to probabilistic edge weights. An underlying fixed network graph is given, and each agent $k$ observes its neighbors with probability $p_k$.  The communication structure is independent of the decision-making process.

A communication strategy where agents observe  rewards and choices of their neighbors in a leader-follower setting is considered in \cite{landgren2018social}. In \cite{wang2020optimal}
only the leader explores and estimates the mean reward of options, while other agents choose the option with highest estimated mean per the leader.

 All of the referenced papers consider homogeneous explore-exploit strategies for the group of agents, except for \cite{landgren2018social,wang2020optimal} that consider leaders and followers. However, when communication among agents is defined by an irregular network graph, e.g., some agents serve as information hubs, group performance can be improved by using heterogeneous explore-exploit strategies.
To understand this, consider an environment with unconstrained resources.  Then, agents can only influence the decisions of one another through the information they share, and the structure of interactions that defines neighbors, i.e., who is sharing information with whom, strongly affects the quality and quantity of information received by each individual.

In this paper we consider the case where all the agent broadcast their rewards and actions to their neighbors with the same probability. This probability captures the cooperative effort of agents in sharing information. In irregular and centralized networks like the multi-star, center agents have more neighbors and thereby  receive more information than peripheral agents. This leads to an imbalanced exploitation potential across the group \cite{sueur2012social,li2019influence}, and group performance degrades with increasing number of peripheral agents. We investigate improving group performance by leveraging heterogeneity in the exploitation potential of agents. To do so we propose  heterogeneous explore-exploit strategies that require center agents to explore more and thus increase the exploitation potential of peripheral agents.

The paper is organized as follows. In Section \ref{sec:problem} we provide the problem formulation and notation. Section \ref{sec:algo} presents the proposed algorithm and intuition. We analyze  performance of the proposed algorithm in Section \ref{sec:performance} and provide improved theoretical bounds for the expected cumulative group regret.  In Section \ref{sec:simulation} we show numerical simulations to illustrate and validate the theoretical results. We conclude in Section \ref{sec:conl}.

\section{PROBLEM FORMULATION}\label{sec:problem}
In this section we present the problem formulation and relevant mathematical notations. Consider a group of $K$ agents, each faced with the same $N$-armed bandit problem for $T$ time steps. At each time step $t \in \{1, \ldots, T\}$, each agent chooses an option and receives a stochastic reward associated with the chosen option. Let $X_i$ be a sub-Gaussian random variable that denotes the {\em reward associated with option} $i\in \{1,\ldots,N\}.$ Define $\mu_i=\mathbb{E}(X_i)$ and $\sigma^2_i$ as the {\em expected reward} and {\em variance proxy} associated with option $i$, respectively. Let $i^*=\arg_{i} \max \{\mu_1,\ldots, \mu_N\}$ be the {\em optimal option} with highest expected reward. Define $\Delta_i=\mu_{i^*}-\mu_i$ as the {\em expected reward gap} between option $i^*$ and option $i.$

Let $G(\mathcal{V},\mathcal{E})$ be a fixed undirected network graph that defines the structure of the interactions between agents. This captures the inherent hard communication constraints of the system. Here $\mathcal{V}$ is a set of $K$ vertices such that each vertex in the graph corresponds to an agent. Each edge $e(k,j)\in \mathcal{E}$ in the graph denotes that agent $k$ and agent $j$ are {\em neighbors}. At each time step, each agent broadcasts its reward value and action to its neighbors with {\em broadcasting probability} $p$. Let $\mathbb{I}^t_{\{k,j\}}$ be the indicator random variable that takes value 1 if agent $k$ receives information from agent $j$ at time $t$ and 0 otherwise. Then, for every time $t$,  $\mathbb{E}(\mathbb{I}^t_{\{k,j\}})=p,\forall k,j$ such that $e(k,j)\in \mathcal{E}$, and $\mathbb{E}(\mathbb{I}^t_{\{k,j\}})=0$ otherwise. We define $\mathbb{I}_{\{k,k\}}^t=1,\forall k,t.$

Let $d_k$ be the {\em degree} (number of neighbors) of agent $k$ and $d_{avg}=\frac{1}{K}\sum_{k=1}^K d_k$ be the {\em average degree of the network}. Let $d^{avg}_k$ be the {\em average degree of neighbors of agent $k$}:
\begin{align*}
d^{avg}_k=\frac{1}{d_k}\sum^K_{j=1,e(k,j)\in\mathcal{E}} d_j.
\end{align*}

We focus on {\em mult-star} graphs defined as follows. Let there be $m$ center agents and $K-m$ peripheral agents.  Without loss of generality let each agent $k$, $k \leq m$, be a {\em center agent}.  All center agents are neighbors of one another, i.e., $e(k,j)\in \mathcal{E}, \forall k,j \leq m$, and a center agent's degree $d_k$ is at least $m-1$. Each {\em peripheral agent} $k$, $k > m$, has exactly one neighbor ($d_k=1$), and the neighbor is a center agent. To reduce complexity, we assume the graph is symmetric, which implies that all the center agents have the same number of neighbors. Thus $K-m$ is an integer multiple of $m$. As long as $K > 2$ and $m<K$, the multi-star graph is {\em irregular}, i.e., the degree of center agents differs from the degree of peripheral agents. Let $d_{cen}$ be the degree of each center agent. Then, $d_{cen}=\frac{K-m}{m}+m-1.$ When $m=1$ the graph is a star, the most irregular multi-star graph. When $m=K$, there are no peripheral agents and the graph is all-to-all and therefore regular.

Let $\varphi_k^t$ be a random variable that denotes the {\em option chosen by agent $k\in \{1,\ldots,K\}$ at time $t\in\{1,\ldots,T\}$}. Let $\mathbb{I}_{\{\varphi_t^k=i\}}$ be an indicator random variable that takes value 1 if agent $k$ chooses option $i$ at time $t$ and 0 otherwise. Let $n_i^k(t)$  be the {\em total number of times agent $k$ chooses option $i$ until time $t$} and let $N_i^k(t)$ be the {\em total number of times agent $k$ observes option $i$ until time $t$}. The total number of observations is the sum of the number of samples taken from option $i$ by agent $k$ and the number of broadcasts on option $i$ by its neighbors:
\begin{align}
n_i^k(t)=\sum_{\tau=1}^t\mathbb{I}_{\{\varphi_{\tau}^k=i\}},\:\:\:N_i^k(t)=\sum_{\tau=1}^t\sum_{ j=1}^K\mathbb{I}_{\{\varphi_{\tau}^j=i\}}\mathbb{I}^{\tau}_{\{k,j\}}.
\label{eq:nandN}
\end{align}
Let $\widehat{\mu}_i^k(t)$ denote the estimate of expected reward of agent $k$ for option $i$ at time $t.$ Then, $
\widehat{\mu}^k_i(t)=\frac{S_i^k(t)}{N^k_i(t)}$, where $S_i^k(t)=\sum_{\tau=1}^t\sum_{j=1}^K X_i\mathbb{I}_{\{\varphi^j_{\tau}=i\}}\mathbb{I}^{\tau}_{\{k,j\}}.$

{\em Expected regret} is defined as the expected loss suffered by agents by sampling sub-optimal options. Let $R(t)$ be the {\em cumulative group regret} at time $t.$ Then {\em expected cumulative group regret} can be computed as
\begin{align}
\mathbb{E}\left(R(t)\right)=\sum_{i=1}^N\sum_{k=1}^K\Delta_i\mathbb{E}\left(n_i^k(t)\right)\label{eq:regret}.
\end{align}

\section{ALGORITHM}\label{sec:algo}
To realize the goal of maximizing cumulative group reward, agents should minimize the number of times they sample sub-optimal options. Each agent employs an agent-based strategy that captures the trade-off between exploring and exploiting by constructing an objective function that strikes a balance between the estimation of the expected reward and the uncertainty associated with the estimate \cite{auer2002finite}.

Since center agents have more neighbors they are more likely to obtain a high number of observations. This reduces the uncertainty associated with their estimate of the expected reward of options. Thus,  identifying the optimal option requires less exploring, which increases their exploitation potential. Since peripheral agents only have one neighbor they are more likely to obtain a low number of observations. Thus, identifying the optimal option requires more exploring, which decreases their exploitation potential. Further, since center agents do less exploring, the usefulness of the information they broadcast is reduced, also decreasing the peripheral agents' exploitation potential. Accordingly, homogeneous sampling rules in irregular, multi-star networks lead to imbalanced exploitation potential across the group and thus degraded group performance.

To improve group performance, we propose heterogeneous explore-exploit strategies that regulate  exploitation potential across the network. When center agents are more exploratory their performance degrades, but the usefulness of the information they broadcast increases and so the performance of peripheral agents improves. When there are more peripheral agents than center agents, and broadcasting probability $p$ is sufficiently high, the performance improvement obtained by peripheral agents outweighs the performance degradation incurred by center agents, and group performance increases.  If $p$ is too small, for example, when broadcasting is costly or risky, center agents do not broadcast enough information to benefit peripheral agents. Thus it doesn't pay for center agents to increase their exploration. Indeed, when $p=0$ all  agents have the same exploitation potential.

Using this intuition, we propose the following heterogeneous sampling rules.
Assume that variance proxy  $\sigma^2_i$ for each option $i$ is known to all agents.
\begin{definition} {\bf{(Heterogeneous Sampling Rules)}}\label{def:samplerule}
The sampling rule $\{\varphi^k_t\}_1^{T}$ of agent $k$ at time $t \in \{1, \ldots, T\}$ is
    \begin{align*}
    \mathbb{I}_{\{\varphi^k_{t+1}=i\}}=\left\{
    \begin{array}{cl} 1 &, \:\:\:i=\arg \max \{Q^k_1(t),\cdots,Q^k_N(t)\}\\
      0 &, \:\:\: {\mathrm{o.w.}}\end{array}\right.
    \end{align*}
    with
    \begin{align}
    Q^k_i(t)&= \widehat{\mu}^k_{i}(t)+C^k_i(t)\label{eq:UCBQ}\\
    C^k_i(t)&=\sigma_{i}\sqrt{\frac{2(1+\alpha_k)(\xi+1)\log t}{N^k_{i}(t)}}\label{eq:Uncertainity}
    \end{align}
    where $\xi>1$ and
    \begin{align}
\alpha_k=\left\{
    \begin{array}{cl}
     \frac{ p^{1-p}(d_k-d^{avg}_k)}{d_k} &, \:\:\: k\leq m \\
     0 &, \:\:\:k>m.
     \end{array}\right.
     \label{eq:alpha}
    \end{align}
\end{definition}

$C_i^k(t)$ in \eqref{eq:Uncertainity} represents {\em agent $k$'s uncertainty in its estimated mean of option $i$}, and Definition~\ref{def:samplerule} implies that for any agent $k$,  when $C_i^k(t)$ is high, agent $k$ will more likely explore. By \eqref{eq:Uncertainity}, $C_i^k(t)$ can be high when $N_i^k$, the number of agent $k$'s observations of option $i$, is low, i.e., when option $i$ is under-sampled.  $C_i^k(t)$ can also be high when {\em agent  $k$'s exploration bias} $\alpha_k>0$ is high.

By \eqref{eq:alpha}, $\alpha_k \neq 0$ only for center agents. Since peripheral agents have one center agent neighbor, $d_k^{avg} \leq d_k$ and thus $\alpha_k \geq 0$ for every center agent $k \leq m$. In fact, $\alpha_k\geq 0$ is designed to grow with increasing irregularity: in the regular case (all-to-all) when $m=K$, $\alpha_k = 0$, and in the most irregular case (star) when $m=1$, $d_1 = K-1$ and $d_1^{avg} = 1$ so $(d_1-d_1^{avg})/d_1 = (K-2)/(K-1)$. Further, $\alpha_k$ grows with $p$ according to the factor $p^{1-p}$, which is large only when center agents are broadcasting their reward values and actions with sufficiently high probability $p$.

\begin{remark}\label{rem:homrule}
The corresponding homogeneous sampling rules can be obtained by letting $\alpha_k=0,\forall k$, in Definition~\ref{def:samplerule}. The heterogeneous sampling rule and the homogeneous sampling rule for peripheral agents are the same.
\end{remark}

By design, the heterogeneous  rules of Definition \ref{def:samplerule} drive center agents to explore more than  the corresponding homogeneous  rules and only when it benefits group performance.

\section{PERFORMANCE ANALYSIS}\label{sec:performance}
In this section we analyze the performance of the heterogeneous sampling rules of Definition~\ref{def:samplerule}.  Using an  approach similar to \cite{auer2002finite} with a few key modifications, we upper bound the expected cumulative group regret $\mathbb{E}(R(T))$. We  show that the bound is lower than the upper bound in the case of the corresponding homogeneous sampling rules, and so we can conclude that the designed heterogeneous strategies provide better group performance than the homogeneous strategies.

By (\ref{eq:regret}), we  upper bound  $\mathbb{E}(R(T))$ if we   upper bound  $\sum_{k=1}^{K} \mathbb{E}(n_i^k(T))$, where  $n_i^k(T)$ is the number of times   agent $k$ samples sub-optimal option $i$ until time $T$. By Definition~\ref{def:samplerule}, agent $k$ chooses sub-optimal option $i$ at time  $t$ if $Q_i^k(t)\geq Q_{i^*}^k(t)$.  Let $\mathbb{I}_{\{Q_i^k(\tau)\geq Q_{i^*}^k(\tau)\}}$ be an indicator random variable that takes value 1 if $Q_i^k(\tau)\geq Q_{i^*}^k(\tau)$ and 0 otherwise. Then,
\begin{align*}
n_i^k(t)=\sum_{\tau=1}^t\mathbb{I}_{\{\varphi_\tau^k=i\}}\leq \sum_{\tau=1}^t\mathbb{I}_{\{Q_i^k(\tau)\geq Q_{i^*}^k(\tau)\}}.
\end{align*}
For each option $i$ and agent $k$ let $\{\eta_i^k(t)\}_1^T$ be a sequence of nonnegative nondecreasing functions. Then,
\begin{align}
\sum_{k=1}^K\mathbb{E}&\left(n_i^k(T)\right)
\leq\sum_{k=1}^K\sum_{t=1}^T\mathbb{E}\left(\mathbb{I}_{\{\varphi_t^k=i\}},N_i^k(t)\leq\eta_i^k(t)\right)\nonumber\\
+&\sum_{k=1}^K\sum_{t=1}^T\mathbb{P}\left(Q_i^k(t)\geq Q_{i^*}^k(t),N_i^k(t)>\eta_i^k(t)\right). \label{neq:step1regret}
\end{align}

It remains to upper bound the right hand side of  (\ref{neq:step1regret}) and we do so in two steps. First, we  upper bound the second summation term of  (\ref{neq:step1regret}) as follows. From  (\ref{eq:UCBQ}) we have
\begin{align}
&\left\{Q_{i}^{k}(t)\geq Q_{i^{*}}^{k}(t)\right\}\subseteq \left\{\mu_{i^{*}}<\mu_{i}+2C_{i}^{k}(t)\right\}\nonumber\\
&\cup\left\{\widehat{\mu}_{i^{*}}^{k}(t)\leq \mu_{i^{*}}-C_{i^{*}}^{k}(t)\right\}\cup \left\{\widehat{\mu}_{i}^{k}(t)\geq \mu_{i}+C_{i}^{k}(t)\right\}\label{eq:seteq}.
\end{align}
For all $k$ let
\begin{align}
\eta_i^k(t)=(1+\alpha_k)\eta_i(t), \;\;\; \eta_i(t)=\frac{8\sigma_i^2(\xi+1)\log t}{\Delta_i^2}.
\label{eq:eta}
\end{align}
Then, by (\ref{eq:Uncertainity}),  $\{\mu_{i^{*}}<\mu_{i}+2C_{i}^{k}(t)\}\cap\{N_i^k(t)> \eta_i^k(t)\}=\emptyset$ where $\emptyset$ is the empty set.
Using (\ref{eq:seteq}) we obtain
\begin{align}
\mathbb{P}&\left(Q_i^k(t)\geq Q_{i^*}^k(t),N_i^k(t)>\eta_i^k(t)\right)\nonumber\\
&\leq \mathbb{P}\left(\widehat{\mu}_{i^{*}}^{k}(t)
\leq \mu_{i^{*}}-C_{i^{*}}^{k}(t)\right)
+\mathbb{P}\left(\widehat{\mu}_{i}^{k}(t)
\geq \mu_{i}+C_{i}^{k}(t)\right).\label{eq:twotail}
\end{align}
To upper bound the right hand side of \eqref{eq:twotail} we use the tail probability bound provided in the following lemma.
\begin{lemma} \label{lem:tailProb}
For any $\xi>1,$ some $\zeta>1$ and for $\sigma_i>0$ in the uncertainty $C_i^k(t)$ given by \eqref{eq:Uncertainity}, we get
\begin{align*}
    \mathbb{P}\left(\big |\widehat{\mu}_i^k(t)-{\mu}_i\big |>C_i^k(t)\right) \leq \frac{1}{\log \zeta}\frac{\log \left((1+d_k)t\right)}{t^{(\xi+1)(1+\alpha_k)}}.
\end{align*}
\end{lemma}
\begin{proof}
From Theorem 1 in the paper \cite{madhushani2020dynamic} we have
for some $\zeta>1$ and for $\sigma_i>0$ there exists a $\vartheta_k>0$ such that
\begin{align*}
    \mathbb{P}\left(\widehat{\mu}_i^k(T)-{\mu}_i>\sqrt{\frac{\vartheta}{N_{i}^{k}(T)}}\right)\leq \frac{\nu\log ((d_k+1)T)}{\exp(2\kappa \vartheta_k)}
\end{align*}
where, $\nu=\frac{1}{\log \zeta},\:\:\:\: \kappa=\frac{1}{\sigma_i^2\left(\zeta^{\frac{1}{4}}+\zeta^{-\frac{1}{4}}\right)^2}.$ Since $\alpha_k\geq 0,\forall k$, we can use $\vartheta_k=2\sigma_{i}^{2}(1+\alpha_k)(\xi+1)\log t$ to get the statement of the lemma.
\end{proof}
Using the statement of Lemma \ref{lem:tailProb} in \eqref{eq:twotail},
\begin{align}
\mathbb{P}\left(Q_i^k(t)\geq Q_{i^*}^k(t),N_i^k(t)>\eta_i^k(t)\right)\nonumber \\
\leq
\frac{2}{\log \zeta} \frac{\log \left((1+d_k)t\right)}{t^{(\xi+1)(1+\alpha_k)}}.
\label{inequality2}
\end{align}
Summing the right hand side of \eqref{inequality2} over $t$ we get
\begin{align}
&\sum_{t=1}^T\frac{\log \left((1+d_k)t\right)}{t^{(\xi+1)(1+\alpha_k)}} \leq \log (1+d_k) \nonumber \\ & \! + \frac{\log (2(1+d_k))(\xi\alpha_k+\xi+\alpha_k)+1}{(\xi\alpha_k+\xi+\alpha_k)^22^{\xi\alpha_k+\xi+\alpha_k}}.
\label{expand}
\end{align}
Since $\log $ is  concave, substituting \eqref{expand} into \eqref{inequality2} we get
\begin{align}
\sum_{k=1}^K\sum_{t=1}^T\mathbb{P}\left(Q_i^k(t)\geq Q_{i^*}^k(t),N_i^k(t)>\eta_i^k(t)\right)\nonumber\\
\leq \frac{2K}{\log \zeta}\log (1+d_{avg})
\nonumber\\+  \frac{2}{\log \zeta}\sum_{k=1}^K\frac{\log (2(1+d_k))(\xi\alpha_k+\xi+\alpha_k)+1}{(\xi\alpha_k+\xi+\alpha_k)^22^{\xi\alpha_k+\xi+\alpha_k}}.\label{eq:tail}
\end{align}
This concludes upper bounding the second summation of (\ref{neq:step1regret}).

As our second step, we upper bound the first summation term of  (\ref{neq:step1regret}) as follows. Since we restrict to  symmetric graphs where all center agents have the same number and type of neighbors, it follows that $\alpha_k=\alpha, \forall k\leq m$. Then, by \eqref{eq:eta} we have $\eta_i^k(t)=(1+\alpha)\eta_i(t),\forall k\leq m$, and $\eta_i^k(t)=\eta_i(t) ,\forall k>m.$ Let $[x]^+=\max\{x,0\}.$

\begin{lemma} \label{lem:Bound}
Let $G$ be a symmetric multi-star graph with $m$ center agents and $K-m$ peripheral agents. Let $\{\eta_i^k(t)\}_1^T$  be the sequence of nonnegative nondecreasing functions given by \eqref{eq:eta}. Then
\begin{align*}
\sum_{k=1}^K\sum_{t=1}^T\mathbb{P}\left(\mathbb{I}_{\{\varphi_t^k=i\}},N_i^k(t)\leq \eta_i^k(t)\right)
\leq (K-m)\eta_i(T)\\
+\frac{m}{1+p(m-1)}\left[]1-p\frac{K-m}{m}\right]^+(1+\alpha)\eta_i(T).
\end{align*}
\end{lemma}
\begin{proof}
\normalfont
Recall the definitions of $n_i^k(t)$ and $N_i^k(t)$ in \eqref{eq:nandN}.
Since the communication structure is independent of the decision-making process $\forall k$,
\begin{align}
\mathbb{E}\left(n_i^k(t)\right)+p\sum_{  \stackrel{j=1}{e(k,j)\in\mathcal{E}}}^K \mathbb{E}\left(n_i^j(t)\right)=\mathbb{E}\left(N_i^k(t)\right).\label{eq:probcond}
\end{align}
Since $N_i^k(t)$ is a nonnegative random variable, $N_i^k(t)\leq \eta_i^k(t)\implies \mathbb{E}\left(N_i^k(t)\right)\leq \eta_i^k(t).$ Thus, from (\ref{eq:probcond}), for all $k$,
\begin{align}
\mathbb{E}&\left(n_i^k(t),N_i^k(t)\leq \eta_i^k(t)\right)\nonumber\\
&+p\sum_{  \stackrel{j=1}{e(k,j)\in\mathcal{E}}}^K \mathbb{E}\left(n_i^j(t), N_i^k(t)\leq \eta_i^k(t)\right)\leq \eta_i^k(t). \label{eq:expNum}
\end{align}
To upper bound
\begin{align*}
\sum_{k=1}^K\sum_{t=1}^T\mathbb{P}\left(\mathbb{I}_{\{\varphi_t^k=i\}},N_i^k(t)\leq \eta_i^k(t)\right)
\end{align*}
we maximize $\sum_{k=1}^K\mathbb{E}(n_i^k(t))$ subject to the constraint given by (\ref{eq:expNum}). This can be formulated as a linear programming optimization problem. Since $\eta_{i}(t)$ is a nondecreasing function,
\begin{align*}
\sum_{k=1}^K\sum_{t=1}^T\mathbb{P}\left(\mathbb{I}_{\{\varphi_t^k=i\}},N_i^k(t)\leq \eta_i^k(t)\right)
\leq \sum_{k=m+1}^K\eta_i^k(T)\\
+\sum_{k=1}^m\frac{1}{1+p(m-1)}\left[1-p\frac{K-m}{m}\right]^+\eta_i^k(T),
\end{align*}
and the statement of the lemma follows.
\end{proof}

This concludes upper bounding the first summation of \eqref{neq:step1regret}, and
we are prepared to state the main theorem of the paper.

\begin{theorem} \label{thm:mainresut}
Consider a distributed stochastic bandit problem with $N$ options, $K$ agents, and $T$ time steps. Let communication graph $G$ be a symmetric multi-star graph with $m$ center agents and $K-m$ peripheral agents. If all agents sample according to the heterogeneous sampling rules defined in Definition \ref{def:samplerule}, the expected cumulative group regret satisfies
\begin{align*}
&\mathbb{E}\left(R(T)\right) \leq c_1(K,m,\alpha,p)\sum_{i=1}^N\frac{8\sigma_i^2(\xi+1)\log T}{\Delta_i}\\
&
+\frac{2}{\log \zeta}\sum_{i=1}^N\Delta_i \Bigg( K\log (1+d_{avg}) +(K-m)\frac{\xi\log 4 + 1}{\xi^22^{\xi}} \\
& \quad \quad\quad \quad \quad  + m  \frac{\log (2(1+d_{cen}))(\xi\alpha+\xi+\alpha) + 1}{(\xi\alpha+\xi+\alpha)^22^{\xi\alpha+\xi+\alpha}} \Bigg), \\
& c_1(K,m,\alpha,p) = K-m + \frac{m(1+\alpha)}{1+p(m-1)}\left[1-p\frac{K-m}{m}\right]^+,
\end{align*}
where $\Delta_i$ is the expected reward gap between the options $i^*$ and  $i$, $\sigma_i^2$ is the variance proxy, and $\xi,\zeta>1$.
\end{theorem}

\begin{proof}
\normalfont
The result follows from \eqref{eq:regret}, (\ref{neq:step1regret}), (\ref{eq:tail}) and Lemma \ref{lem:Bound}.
\end{proof}

\begin{remark}\label{rem:homregret}
Recall that under the corresponding homogeneous sampling rules we have $\alpha_k=0, \forall k.$ Thus, we can recover the expected cumulative group regret bound for the homogeneous sampling rules as follows:
\begin{align*}
&\mathbb{E}\left(R(T)\right) \leq c_2(K,m,p)\sum_{i=1}^N\frac{8\sigma_i^2(\xi+1)\log T}{\Delta_i}\\
&
+\frac{2}{\log \zeta}\sum_{i=1}^N\Delta_i \Bigg( K\log (1+d_{avg})  +(K-m)\frac{\xi\log 4 + 1}{\xi^22^{\xi}}
\\ & \quad \quad\quad \quad \quad \quad + m  \frac{\log (2(1+d_{cen}))\xi + 1}{\xi^22^{\xi}} \Bigg), \\
& c_2(K,m,p) = K-m+\frac{m}{1+p(m-1)}\left[1-p\frac{K-m}{m}\right]^+.
\end{align*}
\end{remark}
When the network graph has a large enough ratio of peripheral agents to center agents and a sufficiently high broadcasting probability $p$, i.e. $p\frac{K-m}{m}>1$, we have $
\left[1-p\frac{K-m}{m}\right]^+=0$,
which implies $c_1=c_2=K-m.$ And since $\alpha>0$ we have
\begin{align*}
\frac{\log (2(1+d_{cen}))(\xi\alpha+\xi+\alpha) + 1}{(\xi\alpha+\xi+\alpha)^22^{\xi\alpha+\xi+\alpha}} \\
< \frac{\log (2(1+d_{cen}))\xi + 1}{\xi^22^{\xi}}.
\end{align*}
Plugging these results into the bounds of Theorem \ref{thm:mainresut} and Remark \ref{rem:homregret}, we see that the heterogeneous sampling rules provide a lower theoretical regret bound than the corresponding homogeneous sampling rules, which implies that the heterogeneous sampling rules provide better group performance than the corresponding homogeneous sampling rules.

\section{SIMULATION RESULTS}\label{sec:simulation}
In this section we provide numerical simulations to illustrate results and validate theoretical bounds. For all simulations, we consider 10  options ($N=10$) with Gaussian reward distributions. Expected reward for the optimal option  is $\mu_{i^*}=11$ and for all sub-optimal options $i\neq i^*$ is $\mu_{i}=10$. We let variance associated with all options $i$ be $\sigma_i^2 =1$. Because the expected reward gaps  $\Delta_i=1$, $i\neq i^*$, are equal to the variances $\sigma_i^2 = 1$, it is a challenging problem to distinguish the optimal option from the sub-optimal options. For all simulations, we consider 1000 time steps ($T=1000$) and use 1000 Monte Carlo simulations with $\xi=1.01$.

We show simulation results for the performance of a group of $K=36$ agents  that communicate over two different symmetric multi-star graphs and use the heterogeneous sampling rules of Definition~\ref{def:samplerule}. We compare to the case when agents use the corresponding homogeneous sampling rules of Remark~\ref{rem:homrule}. The first multi-star graph has $m=2$ center agents and $K-m =34$ peripheral agents, with each center agent communicating with 17 peripheral agents and the other center agent. The second multi-star graph has $m=3$ center agents and $K-m=33$ peripheral agents, with each center agent communicating with 11 peripheral agents and the other center agents. In each case, the center agents are interchangeable and the peripheral agents are interchangeable, so the average performance of a center (peripheral) agent is the same as the individual performance of a center (peripheral) agent.

Figure \ref{fig:netprob} shows how average expected cumulative group regret varies with broadcasting probability $p$ for agents using the heterogeneous rules (dotted) and homogeneous rules (solid). Regret is inversely related to  performance: lower group regret implies higher group performance. Results are plotted on the left for the graph with 2 center agents and on the right for the graph with 3 center agents. When $p=0$ there is no communication at all.  So when $p$ becomes even just a little positive and agents learn about options from their neighbors,  regret falls, i.e.,  group performance rises.
\vspace{-10pt}
\begin{figure}[h]
    \centering
    \includegraphics[width=0.48\textwidth]{prob_net.pdf}
    \caption{Average expected  cumulative group regret for $K=36$ agents at time $t=1000$ as a function of broadcasting probability $p$ with communication over a symmetric multi-star graph. Left  corresponds to 2 center and 34 peripheral agents. Right  corresponds to 3 center and 33 peripheral agents. Dotted lines and solid line shows average regret when agents use  heterogeneous and homogeneous sampling rules, respectively.
    }
    \label{fig:netprob}
    \vspace{-10pt}
\end{figure}

In the case of the homogeneous rules, as $p$ increases through intermediate values, center agents do less and less exploring and the usefulness of the information received by peripheral agents decreases. This leads to increased regret for peripheral agents, and the group overall, and thus degraded group performance. When $p$ approaches 1,  center agents receive sufficient  information from their peripheral neighbors such that their improved performance outweighs the degraded performance of peripheral agents.
This leads to a final decrease in group regret and increase in group performance.

The improvement in performance provided by the heterogeneous rules relative to the homogeneous rules, as predicted by Theorem~\ref{thm:mainresut} and Remark~\ref{rem:homregret}, can be clearly seen in Figure~\ref{fig:netprob} by observing how much lower the dotted regret curve is than the solid regret curve.   The growth in regret in the  homogeneous case, as $p$ increases through intermediate values, is avoided in the heterogeneous case. This is because, by design, center agents are biased toward more exploring, which improves the information that peripheral agents receive. The group performance increase that comes, as $p$ increases further, occurs in the heterogeneous case well before $p$ approaches 1.

The influence of irregularity of the graph can be observed in Figure~\ref{fig:netprob} by comparing the left plot (2 center agents and more irregular) to the right plot (3 center agents and less irregular). The results suggest that performance is higher with more center agents, i.e., with greater regularity in the graph.

Figure \ref{fig:nettime} shows expected cumulative regret as a function of time $t$ for  center (blue), peripheral (pink), and average (black) agents,  when $p=0.8$ and agents use the heterogeneous rules (dotted) and homogeneous rules (solid). Results are plotted on the left for the graph with 2 center agents and on the right for the graph with 3 center agents.  It can be observed that, as predicted for the heterogeneous rules, the peripheral agent performance increases and the center agent performance decreases, such that group performance (as represented by the average agent) improves.  Further, a comparison of left and right plots suggests that group performance improves with more center agents (more regularity).
\vspace{-10pt}
\begin{figure}[h]
    \centering
    \includegraphics[width=0.48\textwidth]{Images/time_net.pdf}
    \caption{Expected cumulative regret of center agent, peripheral agent, and average agent for $K=36$ agents as a function of time $t$ for $p=0.8$ and the same two symmetric multi-star graphs as in Figure~\ref{fig:netprob} (on the left with 2 center agents and on the right with 3 center agents) where agents use heterogeneous (dotted) and homogeneous (solid) sampling rules.}
    \label{fig:nettime}
    \vspace{-10pt}
\end{figure}

\section{CONCLUSIONS}\label{sec:conl}
We have designed and analyzed new heterogeneous rules for how a group of agents that share information over a network should sample an uncertain environment to maximize group reward. We consider communication networks defined by symmetric multi-star graphs, since these exemplify realistic settings in which some agents may have many neighbors while others may have very few.  Using the multi-armed bandit problem as the explore-exploit framework, we show how sampling rules for center agents that favor exploring over exploiting make the information that center agents broadcast to their neighbors more useful, thereby increasing the total reward accumulated by the group.

Our analysis and design advance  understanding of the role that heterogeneity does and can play in collective decision making.  And our demonstration that heterogeneity can be leveraged to improve the performance of a cooperative multi-agent system suggests that further investigation is warranted.

\bibliographystyle{IEEEtran}


\section{Introduction}

Meta-analysis is a common approach used to combine quantitative evidence across studies to generate a global exposure/treatment effect that is more precise than estimates obtainable from individual studies alone.  Traditionally, this has been achieved by obtaining summary statistics from published studies and then combining those estimates using methods of meta-analyses. Although it is cost effective and easy to implement, meta-analyses tend to be underpowered and prone to ecological and confounding bias \citep{doi:10.1002/sim.2768, doi:10.1002/jrsm.4}. Individual patient data (IPD) meta-analysis provides a potential solution to avoid such biases, along with modeling flexibilities such as accounting for the correlation between endpoints and dealing with missing data at the patient level \citep{Riley2010MetaanalysisOI}. If one has access to the original patient level data, there may be a choice between a fully specified multivariate IPD  and a two-stage IPD approach. The full multivariate approach generally uses mixed-effects multilevel regressions to model between and within heterogeneity and quantify the effect of interest in a single model. Although this approach is considered flexible, it may be challenging for conducting and communicating the findings, especially regarding visualization with the hallmark forest plot. The alternative IPD approach involves modeling the data in two stages. In the first step, study-level estimates are obtained using separate regression models. In the second step, standard methods of meta-analysis are used to obtain an overall estimate. An important disadvantage of standard methods of meta-analysis is the assumption that effect size estimates being combined are independent. This assumption is violated when multiple outcomes are interest for synthesis across studies. To avoid the dependence of the effect sizes, several ad-hoc methods have been proposed including averaging the effect sizes and selecting one effect size per study. These ad-hoc approaches may lead to missed opportunities to utilize all available data to address the relevant research questions \citep{cheung-NR2019}.

More principled approaches have been proposed to deal with correlated effects when conducting IPD meta-analysis.  These advances include multivariate meta-analysis which has been used to jointly synthesize the outcomes observed across studies to estimate multiple pooled effects simultaneously \citep{cite-key}. Another approach is  the three-level meta analytic model \citep{cheung2013, konstantopoulos-RSM2011,van-den-noortgate-BRM2013}, which has been used to account for dependence of effect sizes within clusters. This approach considers participants to contribute to only one effect size, so the non-independence is primarily introduced due to the nested structure of the effects \citep{cheung-NR2019}.
An additional approach is based on a two-stage meta-analysis that uses summary measures. In this approach, dependency among effect sizes is handled via robust variance estimation (RVE) in which the dependence between the endpoints is not explicitly modeled, but instead the standard errors for the overall treatment effect or meta-regression coefficients are adjusted \citep{hedges-RSM2010}. This approach may require making a reasonable guess about the between outcome correlation to estimate the between-study variance and to approximate the optimal weights.

In this paper, we propose an innovative approach: a hierarchical meta-analysis for the settings in which each cohort study provides multiple endpoints,  resulting in correlated sampling errors of the estimated effects. The work is motivated by a project that involves the integration of data from six longitudinal cohorts, each one using using multiple inter-related tests and assessment tools to measure child cognition. We assume that all endpoints are measuring an effect of exposure on a common underlying construct  A major strength of our approach is that it will enable us to pool data from diverse endpoints within each cohort and then assess the degree to which consistent patterns emerge across cohorts.

In the proposed approach, we first derive the estimates and the standard errors by fitting univariate regression models for each separate outcome of interest. In the second stage, we  combine the summary measures  using a random-effects model. In this stage, unlike the existing methods of two-stage IPD analysis, we account for the correlated effect sizes within each cohort. First, some of the endpoints within a cohort were highly correlated causing unstable results. To overcome this problem, we obtain robust covariance estimates for effect size estimates within cohorts which we use at the second synthesis stage of the hierarchical meta-analytic approach which not all endpoints was observed for all children within a cohort. To address this issue, we derived a formula for the pairwise correlation between estimated effects using an adjustment, reflecting the fact that not all endpoints are observed for all children. To facilitate the creation of forest plots, downloading the results, and view information about the endpoints included in the analysis in an efficient way, we developed a web-based application which is written using the Shiny library and is hosted on a server using the open-source Shiny Server software (Rstudio and Inc., 2014).  In the last stage of our hierarchical meta analytic approach, we combine the independent effect size estimates obtained for each cohort in a random-effects model to obtain a global measure of the effect size across cohorts \citep{whitehead-book2002,10.1093/biomet/asq006}.

We compare and contrast the findings from proposed approach and a full multivariate analysis to answer a question that naturally arises when results of these two models coincide. Previous studies evaluated this question in different settings. Olkin and Simpsin \citep{olkin-sampson-BIOMETRICS1998} showed that in the case of comparing multiple treatments and a control with respect to a continuous outcome, the traditional meta-analysis based on estimated treatment contrasts is equivalent to the least squares regression analysis of individual patient data if there are no study-by-treatment interactions and the error variances are constant across trials. Mathew and Nordstorm \citep{mathew-nordstrom-BIOMETRICS1999} claimed that the equivalence holds even if the error variances are different across trials. Empirically, meta-analysis using original data has been found to be generally similar but not identical to meta-analysis using summary statistics. Whitehead and Lin and Zeng \citep{whitehead-book2002,10.1093/biomet/asq006}, showed that for all commonly used parametric and semi-parametric models, there is no asymptotic efficiency gain by analyzing original data if the parameter of main interest has a common value across studies, the nuisance parameters have distinct values among studies, and the summary statistics are based on maximum likelihood. More recently Kontopantelis \citep{doi:10.1002/jrsm.1303} conducted a comprehensive simulation study to compare one-stage and two-stage IPD analysis and concluded that a fully specified one-stage model is preferable especially when investigating interactions. We extend the results from these existing studies to the setting where there are correlated endpoints across multiple cohorts.

The remainder of the article is organized as follows. In Section 2, we introduce our motivating application based on a meta-analysis of correlated endpoints measures the effect of prenatal alcohol exposure (PAE) on cognition in six large cohort studies. In Section 3, we introduce notation and describe the two-stage analysis and modeling framework used to combine multiple correlated endpoints within a single cohort. In Section 4, we present the modeling framework used to combine pooled effect size estimates across cohorts. In section 5, we compare/contrast the proposed approach with the corresponding one-stage approach via simulation studies. In section 6, we illustrate our method using data from our motivating application.  Finally, in section 7, we discuss the strengths and limitations of our method.

\section{Motivating Application}
Evidence from animal models and epidemiological studies has linked PAE to a broad range of cognitive and behavioral deficits, growth impairment, and physical anomalies, which are known collectively as fetal alcohol syndrome disorders (FASD) \citep{doi:10.1097/01.ALC.0000145691.81233.FA,doi:10.1111/j.1530-0277.2007.00585.x,Cartere20160775,doi:10.1111/acer.14040}. Fetal alcohol syndrome (FAS), the most severe of the FASD, is characterized by distinctive craniofacial dysmorphology (small palpebral fissures, flat philtrum, thin vermillion), small head circumference, and growth retardation \citep{stratton-book1996}; \citep{hoyme-PEDIATRICS2005} while partial FAS (PFAS) is diagnosed in the presence of facial dysmorphology, a history of PAE and growth retardation, microcephaly, or central nervous system (CNS) impairment. Individuals with PAE who lack the characteristic pattern of dysmorphic features but exhibit cognitive and/or behavioral impairment are often given diagnosis of alcohol-related neurodevelopmental disorder (ARND), which is the most prevalent FASD.  Although the diagnosis of ARND requires a confirmed history of maternal alcohol consumption during pregnancy, there is virtually no information in the scientific literature regarding the levels of PAE associated with an increased risk of clinically significant adverse effects. In a recently updated set of clinical guidelines \citep{hoyme-PEDIATRICS2005}, risky drinking during pregnancy was defined as 6 or more drinks/week for a 2-week period or 3 or more drinks per occassion consumed on at least 2 occasions. However, in human studies, there is little empirical evidence of adverse effects at these levels of exposure.

Between 1975 and 1993, the National Institutes of Health (NIH) funded six longitudinal cohort studies in four U.S. cities (Seattle, Atlanta (two cohorts), Pittsburgh (two cohorts), and Detroit); these are described briefly in Section 1.2. In all but one of these studies, mothers were recruited and interviewed prospectively about their alcohol use during pregnancy and their children were followed longitudinally from infancy through young adulthood; the Atlanta 2 study \citep{doi:10.1002/dev.20215} recruited the mothers shortly following delivery, interviewed them about their drinking during pregnancy, and followed the children as in the other studies. All six studies had measured a variety of neurodevelopmental endpoints that assess five specific measures of cognitive function: IQ, executive function, academic achievement, and learning and memory.

\subsection{The Six Cohort Studies}

\subsubsection{Seattle Cohort\normalfont{\citep{Streissguthetal1981}}  \hfill birth years: 1975-76}
All women who were enrolled in prenatal care by the 5th month of pregnancy at two large Seattle hospitals were eligible to participate. To  ascertain PAE,  participating mothers (N = 1529) were administered a Quantity-Frequency-Variability interview \citep{doi:10.15288/qjsa.1968.29.130} regarding alcohol, tobacco, and drug use for two time periods:  during pregnancy and just prior to pregnancy recognition. 462 newborns were selected based on an algorithm derived from maternal absolute alcohol (AA)/day, alcohol use/occasion, volume variability, and frequency of intoxication constructed to over-represent infants born to heavier drinkers. Controls included both abstainers and light drinkers. Infants were followed up at 8 and 18 months and 4,  7,  11,  14, 21,  25,  and 30 years. Although cohort retention was high (e.g., 82\% at 14 years), other children not initially selected whose mothers had been interviewed during pregnancy were added at  follow-
up assessments to keep the sample size close to 500 at each examination.

\subsubsection{Atlanta Cohort 1\normalfont{\citep{article}} \hfill birth years: 1980-86}
527 low socioeconomic status (SES), pregnant women were recruited at their first prenatal  visit at an urban Atlanta hospital serving a primarily African American, low income population. Women who reported drinking at least 1 oz AA/week during pregnancy were recruited. Nondrinkers, who were similar in demographic background, were recruited at the same time to serve as controls. Women were interviewed at recruitment about their alcohol and drug use; the majority reported drinking on weekends in a ‘binge‘ pattern. Infants were evaluated following birth. Sub-samples were followed up at 6 and 12 months and 7, 14, and 22 years.

\subsubsection{Atlanta Cohort 2\normalfont{\citep{Brownarticle}} \hfill birth years: 1992-1994}
306 mothers and their infants were recruited shortly after delivery at an urban Atlanta hospital.  111 reported having drunk alcohol during pregnancy, 71 of whom also had used cocaine (based on self-report or urine screen); 44 used cocaine but no alcohol; 151 did not drink alcohol or use cocaine. All participants were English speaking, 19 years or older, and had singleton births; most were African American and low SES. The infants were assessed at 2 and 8 years.

\subsubsection{Pittsburgh Cohort 1\normalfont{\citep{DAY1991329}}  \hfill birth years: 1983-86}
Participants were recruited from the prenatal clinic at a maternity hospital if they were English-speaking, age 18 or older,  and in their 4th  or 5th gestational month. A total of 1360 women were interviewed regarding their usual, maximum, and minimum consumption of beer, wine, liquor, tobacco, marijuana, and other illicit drugs prior to pregnancy and during the first trimester. Two cohorts were selected: (1) women who drank 3 or more drinks/week in the first trimester and a random sample of women who drank less or abstained, and (2) women who used marijuana at least twice/month in the first trimester and a random sample of those who used less or abstained.  Because women could be selected for either or both cohorts, there was a 48\% overlap. The cohorts were combined for this analysis: The birth sample consisted of 763 live singleton infants. The alcohol, tobacco, and drug use interview was repeated in the 7th gestational month and at delivery, when second and third trimester substance use information was obtained. The cohort consisted of women who were predominantly low income and of fairly equal numbers of Caucasian and African American women. Participants were followed up at 8 and 18 months, and 3, 6, 10, 14, 16, and 22 years.

\subsubsection{Pittsburgh Cohort 2\normalfont{\citep{Galearticle}} \hfill birth years: 1988-93}
English-speaking women in their 4th or 5th month of pregnancy attending the prenatal clinic at a large inner-city hospital who were 18 years old or older were interviewed regarding their usual, maximum, and minimum consumption of cocaine, alcohol, marijuana, tobacco, and other drugs prior to pregnancy and during the 1st trimester. Every woman who reported any cocaine/crack use during the first trimester was enrolled in the study cohort, as was the next woman interviewed who reported  no cocaine or crack use during both the year prior to pregnancy and the first trimester. Although crack/cocaine use was the criterion for recruitment, a large proportion of these women also drank moderate-to-heavy levels of alcohol; some of the adverse effects subsequently seen on development were independently attributable to the alcohol exposure. The alcohol and drug use interview was repeated at the  end of the second and third trimesters, and offspring were assessed at 1, 3, 7, 10, 15, and 21 years.   The birth cohort consisted of 295 women and infants; the women were predominantly of low socioeconomic status and were roughly equally divided by Caucasian and African American race.

\subsubsection{Detroit Cohort\normalfont{\citep{doi:10.1111/j.1530-0277.1993.tb00744.x}} \hfill birth years: 1986-89}
All women (N > 6400) enrolling in the antenatal maternity clinic at a large, inner-city hospital were interviewed regarding their alcohol use at their first prenatal visit (M = 23.4 wk gestation; SD=7.9), using a timeline follow-back interview \citep{jacobson-PREDIATRICS2002}. The mother provided detailed information about her alcohol consumption on a day-by-day basis during the previous 2 weeks with recall linked to specific times of day and activities and was also asked to recall her day-by-day drinking during a typical week around time of conception. Moderate and heavy drinking women were overrepresented in the sample by including all women reporting at least 0.5 oz AA at conception and a random sample of approximately 5\% of the lower level drinkers and abstainers. The 2-week timeline follow-back interview was repeated at each prenatal clinic visit (M = 5.4 visits). To reduce the risk that alcohol might be confounded with cocaine exposure, 78 heavy cocaine (<2 days/wk), light alcohol (<7 drinks/wk) users were also included in the final sample, which consisted of 480 pregnant women and their children. Participants were followed up at 6.5, 12, and 13 months and 7, 14, and 19 years.

\section{Notation and Model Formulation}

Let $Y_{ijk}$ be the random variable representing response $k$ for individual $j$ in cohort $i$,
$k = 1, \ldots, K$, $j = 1, \ldots, J_i$, where $J_i$ is the number of individuals in cohort $i$, $i = 1, \ldots, I$.
Let $A_{ij}$ be the exposure of interest (i.e. prenatal alcohol exposure) for individual $j$ in study $i$
and $S_{ij}$ be their corresponding propensity score.

\begin{equation} \label{eq-Y-linear-model-separate}
Y_{ijk} = \alpha_{ik} + B_{ik} \, A_{ij} + \gamma_{ik} \, S_{ij} + E_{ijk} \,,
\end{equation}
\noindent
where $B_{ik}$ is the effect of a one-unit increase in $A_{ij}$ (alcohol volume) on the mean for response $k$ in cohort $i$
given the propensity score $S_{ij}$, $j = 1, \ldots, J$, $k = 1, \ldots, K$, $i=1, \dots, I$.
The propensity score is calculated separately for each cohort in our setting
since the sets of covariates measured differ between cohorts \citep{Hwang2020}.
We assume here that conditioning on the propensity score renders the exposure variable independent of all confounders and so that it is sufficient to condition on $S_{ij}$ rather than the confounders themselves \citep{doi:10.1198/016214504000001187}.

The parameter $\gamma_{ik}$ characterizes the effect of the propensity score for a given level of alcohol exposure
and $E_{ijk}$ is the error term which has mean zero and variance $\sigma^2_{ik}$, $k = 1, \ldots, K$, $i=1, \dots,I$.

We suppose that the effects $B_{ik}$ vary about some average exposure effect, with,
\begin{equation} \label{eq-average-exposure-Bk}
B_{ik} \sim N(\beta_{i}, \phi_{i}) \,,
\end{equation}
\noindent
where $\beta_i$ is the exposure effect for cohort $i$ and $\phi_i$ represents the heterogeneity of the response-specific exposure effects for cohort $i$.

Across cohorts, we suppose that the average exposure effects $\beta_{i}$ vary about some average exposure effect, with,
\begin{equation} \label{eq-average-exposure-Bk_2}
\beta_i \sim N(\beta_{o}, \eta^2) \,,
\end{equation}
\noindent
where  $\beta_{o}$ represents the ``average effect'' of a one-unit increase in the exposure across all cohorts, and is our parameter of ultimate interest.
The variance $\eta^2$ in (\ref{eq-average-exposure-Bk_2}) reflects the extent of heterogeneity of the cohort-specific exposure effects.

In the next two subsection, we describe a two-stage approach to estimation and inference with data from a single cohort and in Section 4, we show how to synthesize cohort-specific exposure effects to obtain an estimate for the average effect of a one-unit increase in the exposure across all cohorts.

\subsection{Stage I Estimation for a Single Cohort}
In this section, we temporarily omit the subscript $i$ and describe a two-stage approach to estimate the average exposure effect for a single cohort where the effects are correlated.
Before model fitting we standardize the responses so that they have the same first two moments as the Full-Scale IQ variable
which has a mean of 100 and a standard deviation of 15.
By conducting this standardization the exposure effects can be expressed
in terms of the decrement in IQ  associated with a one-unit increase
in prenatal alcohol exposure \citep{axelrad-EHP2007}.

For the first stage, we fit separate linear models for each response, assuming
\begin{equation}
Y_{jk} = \alpha_{k} + B_{k} \, A_{j} + \gamma_{k} \, S_{j} + E_{jk} \,,
\end{equation}
\noindent
where $B_{k}$ is the effect of a one-unit increase in $A_{j}$ (alcohol volume) on the mean for response $k$,
given the propensity score $S_{j}$, $j = 1, \ldots, J$, $k = 1, \ldots, K$.
The parameter $\gamma_{k}$ characterizes the effect of the propensity score for a given level of alcohol exposure
and $E_{jk}$ is the error term which has mean zero and variance $\sigma^2_{k}$, $k = 1, \ldots, K$.

We suppose that the effects $B_{k}$, $k = 1, \ldots, K$
vary about some average exposure effect, with
\begin{equation}
B_{k} \sim N(\beta, \phi),
\end{equation}
\noindent
independently and identically distributed and $\beta$ is the average exposure effect. The variance $\phi$ reflects the extent of heterogeneity of the response-specific exposure effects for a single cohort.

If we let $X_{jk} = (1, A_{j}, S_{j})'$ be the covariate vector, we can be write
\begin{equation} \label{eq-Y-giv-BE-linear-model-compact}
Y_{jk} = X'_{jk} \, \theta_{k} + E_{jk}
\end{equation}
where $\theta_{k} = (\alpha_{k}, B_{k}, \gamma_{k})'$.
We assume $E_{jk}\independent (A_{j}, B_{k},S_{j})$ with $E_{jk} \sim N(0, \sigma_{k}^2)$ are i.i.d. for $j = 1, \ldots, J$, $k=1,2,\ldots, K$.
Note that $X_{jk}$ does not vary by response type
since the exposure variable $A_{j}$ and the propensity score $S_{j}$ are individual level covariates,
but we retain this notation for generality.

We next define $K \times 1$ vectors $Y_j = (Y_{j1}, Y_{j2}, \ldots, Y_{jK})'$,
$\alpha = (\alpha_1, \ldots, \alpha_K)'$,
$B = (B_1, \ldots, B_K)'$ and $\gamma = (\gamma_1, \ldots, \gamma_K)'$ and a $K \times 3K$ covariate matrix
\begin{equation}
X'_j = \left[
\begin{array}{ccccc}
X'_{j1} & 0          & 0 & \cdots & 0 \\
0          & X'_{j2} & 0 & \cdots & 0 \\
\vdots   &            & \ddots &   & 0 \\
0          & \cdots  & 0         & 0 & X'_{jK} \\
\end{array}
\right] \,,
\end{equation}
\noindent
where the $0$s in this matrix refer to $3\times 1$ vectors of $0$s. The model given by (\ref{eq-Y-giv-BE-linear-model-compact}) can then be represented in a unifying model
\begin{equation}
Y_j = X'_j \, \theta + E_j
\end{equation}
\noindent
where
$\theta = (\theta'_1, \ldots, \theta'_K)'$ is a $3 K \times 1$ vector of parameters,
$E_j = (E_{j1}, \ldots, E_{jK})'$ and $E_j \sim N(0, \Sigma)$,
where $\Sigma$ is a $K \times K$ covariance matrix with diagonal entries
$\Sigma_{kk} = \sigma_k^2 = {\rm var}(E_{jk})$.
The off-diagonal entries
$\Sigma_{kl} = \sigma_{kl}={\rm cov}(E_{jk}, E_{j \ell})$
accommodate a conditional dependence (given $X_{jk}$, $X_{jl}$, $B$) between the responses from the same individual.

In the second stage of estimation these estimates are pooled over
to obtain a single estimate of the global measure of the causal effect denoted by $\beta$ in (\ref{eq-average-exposure-Bk}). We begin the second stage by estimating the covariance between the errors  ${\rm cov}(E_{jk}, E_{jl}) = \sigma_{kl}$, $l \ne k = 1, \ldots, K$, characterizing the dependence between the stage I estimators $\widehat{\theta}_1, \ldots, \widehat{\theta}_K$

The challenge in estimating the covariance between the errors  ${\rm cov}(E_{jk}, E_{jl}) = \sigma_{kl}$, $l \ne k = 1, \ldots, K$ is
that not all individuals contribute data for all responses. To accommodate this fact, we introduce the indicators $R_{jk} = I(Y_{jk}~\mbox{is observed})$, $k = 1, \ldots, K$, and assume that the responses are missing at random (MAR) according to Little and Rubin (2019).

If $S_{jk}(\theta_{k})=X_{jk}(Y_{jk}-X_{jk}' \, \theta_{k})$
is the desired contribution from individual $j$ to the score function for $\theta_{k}$
given $B$,
the observed data score equation for estimating $\theta_{k}$ at stage I can be written as
\begin{equation} \label{eq-stageI-Sk-thetak}
S_{k}(\theta_{k})= \sum_{j=1}^{J} R_{jk} S_{jk}(\theta_{k}) = 0 \,,
\end{equation}
the solutions to which are
\begin{equation}
\widehat{\theta}_{k}=\sum_{j=1}^{J} R_{jk} \, (X_{jk} \, X'_{jk})^{-1} \, X_{jk} \, Y_{jk} \,, \quad k = 1, \ldots, K \,.
\end{equation}
\noindent
The maximum likelihood estimate of $\sigma^2_{k}$ is
\begin{equation} \label{eq-mle-sigma2k}
\widehat{\rm var}(E_{jk})
= \widehat{\sigma}^{2}_{k}
=\dfrac{\sum_{j=1}^{J} R_{jk} \, (Y_{jk}-X'_{jk} \, \widehat{\theta}_{k})^2}{n_k}
\end{equation}
where $n_k=\sum_{j=1}^{J} R_{jk}$
is the number of individuals contributing to the estimation of $\theta_{k}$, $k=1,2,\dots, K$.
We also obtain the maximum likelihood covariance estimate as
\begin{equation} \label{eq-mle-sigmakl}
\widehat{\rm cov}(E_{jk}, E_{jl})
=\widehat{\sigma}_{kl}
= \dfrac{\sum_{j=1}^{J} R_{jk} \, R_{jl}  \,(Y_{jk}- X'_{jk} \, \widehat{\theta}_{k}) \,  (Y_{jl}- X'_{jl} \, \widehat{\theta}_{l})}
            {n_{kl}}
\end{equation}
\noindent
where $n_{kl} = \sum_{j=1}^J R_{jk} \, R_{jl}$,
which is consistent
under a missing at random assumption \citep{little-rubin-book2019}.
We then let $\widehat{\Sigma}$ denote the estimated covariance matrix for the errors
where $\sigma_{kl}$ are adjusted for missing data as

\begin{equation}
\widehat{\rm cov}(E_{jk},E_{jl})_{adjusted}
=\dfrac{\widehat{\sigma}_{kl} \,n_{kl}}
{n_{k}\,n_{l}}
\end{equation}
\noindent

where $n_k=\sum_{j=1}^{J} R_{jk}$ and  $n_l=\sum_{j=1}^{J} R_{jl}$
are the number of individuals contributing to the estimation of $\theta_{k}$ and $\theta_{l}$ respectively.

More details regarding the derivation of the dependence between the estimates $\widehat{\theta}_{k}$ and $\widehat{\theta}_{l}$
are provided in the appendix.

\subsection{Stage II: Synthesis across Responses within a Cohort}
\label{SEC3.2}

To consider the synthesis of estimators across all responses
we note that
\begin{equation*}
E(\widehat{B}_k - \beta) = E \{ (\widehat{B}_k - B_k) + (B_k - \beta) \} = 0
\end{equation*}
\noindent
so $\widehat{B}$ is comprised of $K$ dependent unbiased estimators of $\beta$.
Thus
\begin{equation} \label{eq-Bhat-MVN}
\widehat{B} \sim {\rm MVN}(\mu(\beta), \Psi(\phi))
\end{equation}
\noindent
asymptotically,
where $\mu(\beta)$ is a $K \times 1$ vector with each element equal to $\beta$.
If we let $\widehat{\Psi}(\phi) = J^{-1} \, \widehat{\Gamma} + \Delta \, \phi$
then based on (\ref{eq-Bhat-MVN}) we may specify a pseudo-likelihood $PL(\beta, \phi)$ for $(\beta, \phi)$ given by
\begin{equation}\label{eq-PL-beta-tau2}
PL(\beta, \phi)
\propto \dfrac{1}{ (2 \pi)^{K/2} \, \sqrt{\widehat{\Psi}(\phi)|}}
\exp \left(-\frac{1}{2} \, (\widehat{B}- \mu(\beta))'  \, \widehat{\Psi}^{-1}(\phi)  \, (\widehat{B} - \mu(\beta))  \right) \,.
\end{equation}

Note that (\ref{eq-PL-beta-tau2}) could be maximized with respect to $(\beta, \phi)$,
but we proceed in a computationally convenient iterative approach based on (\ref{eq-PL-beta-tau2}).
Given an estimate $\phi^{(r)}$
we compute an estimate $\beta^{(r)}$ based on a linear combination of
$\widehat{B}_1, \ldots, \widehat{B}_K$.
The most efficient linear estimator of $\beta$ has the form
\begin{equation} \label{eq-linear-estimator-beta-hat}
\widehat{\beta} = [\mathbbm{1}' \, [\Psi(\phi)]^{-1} \, \widehat{B}] / [ \mathbbm{1}' \, [\Psi(\phi)]^{-1}  \, \mathbbm{1} ] \,,
\end{equation}
\noindent

so we replace $\Psi(\phi)$ with an estimator
$\widehat{\Psi}(\phi^{(r)}) = J^{-1} \, \widehat{\Gamma} + \Delta \, \phi^{(r)}$.
We could invert $\widehat{\Psi}(\phi^{(r)})$ but in practice it may be difficult and while a generalized inverse could be used,
the weights resulting from this approach were often found to vary greatly in magnitude and even in sign.
We therefore adopted an alternative more stable linear estimate using inverse variance weights,
whereby we replace $\Psi(\phi)$ with
$diag(J^{-1}  \, \widehat{\Gamma}_{kk} + \phi^{(r)} \,,~k = 1, \ldots, K)$ in (\ref{eq-linear-estimator-beta-hat}) to obtain $\widehat{\beta}^{(r)}$.
We then maximize $PL(\widehat{\beta}^r, \phi)$ with respect to $\phi$ to obtain $\phi^{(r+1)}$, with which we recompute $\widehat{\beta}^{(r+1)}$
and repeat iteratively until convergence;
we let $(\widehat{\beta}, \widehat{\phi})$ denote the estimates upon convergence.
A robust variance estimate is then obtained for $\widehat{\beta}$ based on
$\widehat{\Psi}(\widehat{\phi})$ which is given by
$\widehat{\rm var}(\widehat{\beta}) = [ \mathbb{I}' \, \widehat{\Psi}(\widehat{\phi}) \, \mathbb{I}]^{-1}$.
When we consider data from multiple cohorts we reintroduce the subscript $i$ and write the corresponding estimates for cohort $i$ as $\widehat{\beta}_i$ and
\begin{equation}
\widehat{V}_i(\widehat{\beta}_i) = \widehat{\rm var}(\widehat{\beta}_i) = [ \mathbb{I}' \, \widehat{\Psi}_i(\widehat{\phi_i}) \, \mathbb{I}]^{-1}
\end{equation}
respectively, $i=1,\ldots, I$.

\subsection{An Alternative One-Stage (Fully Specified Multivariate) Approach}
Note that the parameters estimated according to the hierarchical meta-analytic approach can alternatively be fitted in one-step via software for fitting hierarchical mixed effect linear models.
To do so we define $K-1$ covariates $T_{jkr} = I(k = r)$, $r = 2, \ldots, K$ which indicates the outcomes is for response $k$, $k = 2, \ldots, K$.
We may put these $K-1$ indicators in vector format and define the $(K-1)\times 1$ vector $T_{jk} = (T_{jk2}, \ldots, T_{jkK})'$.
We consider the response for verbal IQ as the reference type and let $T_{ik2} = 1$ for performance IQ,
$T_{ik3} = 1$ for freedom from distractibility (Table 2).
Then we fit the model
\begin{equation*}
Y_{jk}  = \alpha_1 + B_k \, A_j + \gamma_1 \, S_j + \tau' \, T_{jk} + \zeta' \, S_j \, T_{jk} + E_{jk} \,, \quad k = 1, \ldots, K \,.
\end{equation*}
where $\tau=(\tau_2, \ldots, \tau_K)'$ and
$\zeta = (\zeta_2,\ldots, \zeta_K)'$ are $(K-1) \times 1$ vectors with
$\tau_r = \alpha_r - \alpha_1$ and
$\zeta_r = \gamma_r - \gamma_1$, $r=2,\ldots, K$.
We assume $B_k \sim N(\beta, \phi)$ as specified in (\ref{eq-average-exposure-Bk}) where $B_k$ is the effect of a one-unit increase in $A_j$ on the
mean of response $k$ given the propensity score $S_j$,
$\beta$ is the parameter of ultimate interest representing the ``average causal effect'' of a one-unit increase in the exposure across
all responses within the cohort, and
$\phi$ characterizes the degree of heterogeneity in the effect across responses.
We also assume $E_j = (E_{j1}, \ldots, E_{jK})'$ is a $K \times 1$ error term with
$E_j \sim {\rm MVN}(0, \Sigma)$
with $\Sigma$ a $K \times K$ covariance matrix as in Section 3.1.
The one-step approach involves simultaneous estimation of all fixed effects, $\beta$, $\phi$ and $\Sigma$ at once.
This can be fitted using software for fitting hierarchical linear mixed effects models.
We let
$\tilde{\beta}_i$ denote the estimate of $\beta_i$ obtained from fitting the hierarchical model to the data from cohort $i$ and
$\tilde{V}_i(\tilde{\beta}_i)$ denote the correspondign variance estimate based on the observed information matrix, $i=1,\ldots, I$.

\section{Synthesis across Cohorts}
In the previous section, we described methods for synthesizing data across multiple outcomes to obtain estimates of the global causal effect based on a
two-stage approach and based on fitting a hierchical mixed effect model.
These methods were based on analyzing data from a single cohort. Here, we describe how to combine cohort-specific estimates to obtain a an overall estimate
of a causal effect while accommodating possible heterogeneity.
The approach described in Section 3.2 is an extension of the approach described by Viechtbauer and implemented in the metafor package \citep{metaforR} which deals with
independent estimates; Section 3.2 adapted the methods to deal with dependent effect estimates so what follows is a simplification of the approach for the
last stage of the data synthesis.
We describe it briefly as follows.

We consider $\beta_i$ as the global causal effect of exposure in cohort $i$ reflecting the impact of an increment in the volume of  prenatal alcohol exposure
on the common underlying construct; we let $\hat\beta_i$ be the corresponding estimate.
Note that the studies draw individuals from different populations and so the composition of the samples varies across cohorts.
Moreover, the methods used to measure exposure along with the precise nature of the outcomes differ between studies, even though they were measuring the same latent attributes regarding cognition.
We therefore wish to accommodate a component of variation between studies (heterogeneity) for the true effects which we accomplished by use of a random effects model of the form
\begin{eqnarray}
\hat{\beta}_{i} & = & \beta_{i}+\epsilon_{i}\\
\beta_{i} & = & \beta_\circ + u_{i}
\end{eqnarray}
where
we let $\epsilon_i \sim N(0, \widehat{V}_i(\widehat{\beta}_i))$ reflect the sampling variation of the estimator from cohort $i$ about the true effect $\beta_i$,
and $u_{i}\sim N(0,\eta^2)$ reflects the heterogeneity of the global cohort-specific causal effects across studies.
The parameter $\beta_\circ$ represents the overall global effect which is the parameter of ultimate interest.
Through this variance decomposition then upon introducing the heterogeneity between studies we have
$\widehat{\rm var}(\widehat{\beta}_i)=
\widehat{V}_i(\widehat{\beta}_i) + \eta^2$.
The synthesis is achieved in a simlar spirit to Section 3.2 whereby we consider a pseudo-likelihood of the form
\begin{equation}
\label{eq-PL-beta-eta2}
PL(\beta_\circ, \eta)
\propto
\prod_{i=1}^I \left\{ \dfrac{1}{ (2 \pi)^{I/2} \, \sqrt{(\widehat{V}_i(\widehat{\beta}_i) + \eta^2)}}
\exp \left( - \frac{(\widehat{\beta}_i- \beta_\circ)^2}{2 ( \widehat{V}_i(\widehat{\beta}_i) + \eta^2)} \right) \, \right\}.
\end{equation}
The pooled exposure effect estimate $\widehat{\beta}_\circ$ is obtained as a weighted average of the $\widehat\beta_i$ terms
with cohort weights equal to the inverse of $\widehat{V}(\widehat{\beta}_i)+\hat\eta^2$ where $\hat{\eta}^2$ is obtained as the solution to iteratively
maximizing (\ref{eq-PL-beta-eta2}).
The R package `metafor' can be used to carry out this final stage of the data synthesis.
If the linear model of Section 3.3 is used for simultaneous estimation of the overall causal effect then $\tilde{\beta}_i$ and
$\tilde{V}_i(\tilde{\beta}_i)$  can be used in a similar fashion to obtain the estimator $\tilde{\beta}_o$.

\section{Simulation Studies}
For the simulation studies, we consider $k$ correlated continuous endpoints from a single study. We generated outcomes from the following linear regression model:

\begin{equation}
Y_{jk} = \alpha_k + \beta_k \,X_j + \gamma_k \, Z_j + E_{jk} \,,
\end{equation}

where $Y_{jk}$ be the random variable representing response k for individual j, k = 1,...,K, j = 1,...,J, and  $\beta_k$ is the effect of a one-unit increase in $X_j$ on the mean for response $k$ given the covariate $Z_j$. We let $\beta_k$ vary about some average exposure effect within a study, with $\beta_k \sim N(\beta, \tau^2)$. The parameter $\gamma_k$ characterizes the effect of the covariate for a given level of exposure, $X_{j}$. We also assume $E_{j} = (E_{j1}, \dots, E_{jK} )$ is a $K \times 1$ error term with $E_{j}\sim MVN(0, \Sigma)$ with $\Sigma$ a $K \times K$ covariance matrix. The simulations were performed in different scenarios. We generated the effect size for the exposure, $\beta_k$ from normal distribution with mean 3 and variance $\tau^2$.  Scenarios were created by manipulating the number of outcomes (k) and, varying between-study heterogeneity ($\tau^2$). We consider the scenarios where the number of outcomes is equal to 3, 5, and 10 and $\tau^2$ takes the values 0.10,0.25 and 0.50.
For each combination of the simulation parameters, we generated 1000 datasets with the sample size of 500 for each endpoint. For each dataset, we performed the two types of meta-analysis, i.e. the one based on the proposed approach versus the full multivariate analysis.

We evaluated the performance of the proposed approach in simulation settings previously described, over 1000 iterations. The estimates of interest were the average exposure effect $\beta$. To allow for a comprehensive comparison, performance was assessed on a range of metrics: empirical mean bias (EBIAS), average model based standard error (ASE), empirical standard error (ESE) and coverage probability (CP).  The results are summarized in Table 1. Patterns of empirical mean bias were very low and comparable for the two methods, with the exception of larger between heterogeneity and smaller number of outcomes. In those particular scenarios, the proposed approach was the best performer. Coverage probability for the proposed method was about to the nominal 0.95 for all scenarios considered in this paper .

\section{Prenatal Alcohol Exposure and Cognitive Function in Children}
For this paper we used data from six large cohort studies to assess the effects of PAE on IQ, which is a measure of cognitive function. The proposed hierarchical meta-analytic approach is well-suited to assess the effect of PAE on IQ measure since it enables us to pool data from diverse, correlated endpoints across cohorts. Table 2 lists the endpoints by cohort considered in this paper. To yield sufficiently precise estimates of effect size, we considered a broad set of potential confounders when fitting separate linear models for each endpoint. Since each cohort provided a somewhat different set of control variables, we employed a propensity score approach to adjust for potential confounders (\citealp{AkkayaHocagil_2020}). We estimated the propensity score for each cohort separately and included the propensity score in the linear model as an additional covariate as in model (3.1).

For each outcome k, $\beta_{k}$ was estimated from model (3.1). Table 3 lists the estimated effect size and standard errors from the first stage of the hierarchical meta-analytic approach. The aim of the second stage of the proposed methods is to pool the estimates of the PAE and estimate the cohort specific overall true mean effect $\beta_{i}$ while adjusting for the fact that endpoints are correlated within a cohort and accommodating incomplete information on some endpoints. Table 4 shows the estimated effect sizes and standard errors for each cohort. To compare and contrast the results obtained from our method, we conducted a fully specified multivariate analysis to estimate a pooled effect size for each cohort study using SAS procedure `proc mixed` . Table 4 shows the estimated effect sizes for each cohort obtained from the fully specified multivariate model. Both methods provided impressively similar estimates for the effect sizes and the standard errors. Although the difference was not substantial, these two methods provided slightly different estimates for the between endpoints heterogeneity.

To combine the independent effect size estimates across cohorts and obtain a global effect size estimate of PAE on IQ at age 7 years, we used the R package ``metafor" to pool the estimates resulting from our hierarchical meta-analysis and the fully specified multivariate model. Table 5 provides the estimated effect size of the PAE on childhood IQ across cohorts. The resulting global effect sizes from the two methods were almost identical.

\section{Discussion}

	In this paper, we have proposed an extension of the standard procedure for two-stage IPD meta-analysis. Our procedure follows the same steps as conventional methods for two-stage IPD analysis for making inference about the effect size but extends these analyses by accounting for the correlation between endpoints and by accommodating incomplete data on some endpoints.

	Our approach has several advantages over the one-stage IPD meta-analysis. Firstly, it builds upon the two-stage IPD meta-analysis that practitioners are already familiar with. Secondly with our approach, one can create forest plots to visualize the estimated effect sizes for each endpoint. Thirdly, our approach is less likely to encounter convergence problems compared to the one-stage IPD meta-analysis. Finally, our approach uses the known within study variances, which helps obtain more precise estimates.

	We evaluated and compared our approach with a fully specified multivariate analysis. In simulation scenarios considered in this paper, we observed that the proposed approach can successfully reduce bias relative to the fully specified multivariate approach. Our simulation results suggest that, when the number of endpoints is small and the between endpoints variance is large, our proposed approach outperforms the multivariate analysis.

    We illustrated our approach using data on childhood IQ from six cohorts. We analyzed 18 outcomes from the six cohorts. For our data application, it was important to adjust for the fact that endpoints within a cohort are correlated and some of the endpoints have incomplete data. To obtain a global effect size estimate, we conducted a hierarchical meta-analysis. In the first stage, we obtained effect size estimates for each endpoint separately. In the second stage, we employed our proposed approach to obtain cohort-specific pooled effect size estimates while adjusting for between-endpoint correlation and incomplete data. In the last stage, we combined effect sizes across the cohorts employing a random-effects model. We compared the results from our approach with the results from the fully specified multivariate approach. In this comparison, our method performed well and thus provides an useful innovative tool for performing and interpreting meta-analyses with the correlated effect sizes.