Functional Mixed Membership Models

To aid in computation, we assume that $f^{(k)}$ is a smooth function and is in the $P$-dimensional subspace $\boldsymbol{\mathcal{S}}\subset L^{2}(\mathcal{T})$, spanned by a set of linearly independent square-integrable basis functions $\{b_{1},\dots,b_{P}\}$ ($b_{p}:\mathcal{T}\rightarrow\mathbb{R}$). While the choice of basis functions is user-defined, in practice B-splines (or the tensor product of B-splines for $\mathcal{T}\subset\mathbb{R}^{d}$ for $d\geq 2$) are a common choice due to their flexibility. De Boor (2000) has shown that for any smooth function, the distance between the true function and the approximation constructed using B-splines with uniform knots goes to zero as the number of knots goes to infinity. In practice, the number of basis functions used is often proportional to the average number of time points observed for each functional observation (Section E of the web-based appendix contains a detailed discussion on the choice of $P$). Alternative basis systems can be selected in relation to application-specific considerations.

The multivariate Karhunen-Loève (KL) theorem, proposed by Ramsay and Silverman (2005), can be used to jointly decompose our $K$ stochastic processes. To allow our model to handle higher-dimensional functional data, such as images, we will use the extension of the multivariate KL theorem for different dimensional domains proposed by Happ and Greven (2018). Although they state the theorem in full generality, we will only consider the case where $f^{(k)}\in\boldsymbol{\mathcal{S}}$. In this section, we will show that the multivariate KL theorem for different dimensional domains still holds under our assumption that $f^{(k)}\in\boldsymbol{\mathcal{S}}$, given that the conditions in Lemma 2.2 are satisfied.

We will start by defining the multivariate function $\mathbf{f}(\mathbf{t})$ in the following way:

where $\mathbf{t}$ is a $K$-tuple such that $\mathbf{t}=\left(t^{(1)},t^{(2)},\dots,t^{(K)}\right)$, where $t^{(1)},t^{(2)},\dots,t^{(K)}\in\mathcal{T}$. This construction allows for the joint representation of $K$ stochastic processes, each at different points in their domain. Since $f^{(k)}\in\boldsymbol{\mathcal{S}}$, we have $\mathbf{f}\in\boldsymbol{\mathcal{H}}:=\boldsymbol{\mathcal{S}}\times\boldsymbol{\mathcal{S}}\times\dots\times\boldsymbol{\mathcal{S}}:=\boldsymbol{\mathcal{S}}^{K}$. We define the corresponding mean of $\mathbf{f}(\mathbf{t})$, such that $\boldsymbol{\mu}(\mathbf{t}):=\left(\mu^{(1)}\left(t^{(1)}\right),\mu^{(2)}\left(t^{(2)}\right),\dots,\mu^{(K)}\left(t^{(K)}\right)\right)$, where $\boldsymbol{\mu}(\mathbf{t})\in\boldsymbol{\mathcal{H}}$, and the mean-centered cross-covariance functions as $C^{(k,k^{\prime})}(s,t):=\text{Cov}\left(f^{(k)}(s)-\mu^{(k)}(s),f^{(k^{\prime})}(t)-\mu^{(k^{\prime})}(t)\right)$, where $s,t\in\mathcal{T}$.

Proofs for all results are given in the Supplementary Materials. From Lemma 2.1, since $\boldsymbol{\mathcal{S}}$ is a closed linear subspace of $L^{2}(\mathcal{T})$, we have that $\boldsymbol{\mathcal{S}}$ is a Hilbert space with respect to the inner product $\langle f,g\rangle_{\boldsymbol{\mathcal{S}}}=\int_{\mathcal{T}}f(t)g(t)\text{d}t$ where $f,g\in\boldsymbol{\mathcal{S}}$. By defining the inner product on $\boldsymbol{\mathcal{H}}$ as

where $\mathbf{f},\mathbf{g}\in\boldsymbol{\mathcal{H}}$, we have that $\boldsymbol{\mathcal{H}}$ is the direct sum of Hilbert spaces. Since $\boldsymbol{\mathcal{H}}$ is the direct sum of Hilbert spaces, $\boldsymbol{\mathcal{H}}$ is a Hilbert space (Reed and Simon, 1972). Thus, we have constructed a subspace $\boldsymbol{\mathcal{H}}$ using our assumption that $f^{(k)}\in\boldsymbol{\mathcal{S}}$, which satisfies Proposition 1 in Happ and Greven (2018). Letting $\mathbf{f}\in\boldsymbol{\mathcal{H}}$, we can define the covariance operator $\boldsymbol{\mathcal{K}}:\boldsymbol{\mathcal{H}}\rightarrow\boldsymbol{\mathcal{H}}$ element-wise in the following way:

Since we made the assumption $f^{(k)}\in\boldsymbol{\mathcal{S}}$, we can simplify Equation 7. Since $\boldsymbol{\mathcal{S}}$ is the span of the basis functions, we have $f^{(k)}(t)-\mu^{(k)}(t)=\sum_{p=1}^{P}\theta_{(k,p)}b_{p}(t)=\mathbf{B}^{\prime}(t)\boldsymbol{\theta}_{k},$ for some $\boldsymbol{\theta}_{k}\in\mathbb{R}^{P}$ and $\mathbf{B}^{\prime}(t)=[b_{1}(t)\cdots b_{P}(t)]$. Thus, we can see that

and we can rewrite Equation 7 as $\left(\boldsymbol{\mathcal{K}}\mathbf{f}\right)^{(k)}(\mathbf{t})=\sum_{k^{\prime}=1}^{K}\int_{\mathcal{T}}\mathbf{B}^{\prime}(s)\text{Cov}\left(\boldsymbol{\theta}_{k^{\prime}},\boldsymbol{\theta}_{k}\right)\mathbf{B}\left(t^{(k)}\right)f^{(k^{\prime})}(s)\text{d}s,$ for some $\mathbf{f}\in\mathcal{H}$. The following lemma establishes conditions under which $\boldsymbol{\mathcal{K}}$ is a bounded and compact operator, which is a necessary condition for the multivariate KL decomposition to exist.

Assuming the conditions specified in lemma 2.2 hold, by the Hilbert-Schmidt theorem, since $\boldsymbol{\mathcal{K}}$ is a bounded, compact, and self-adjoint operator, we know that there exists real eigenvalues $\lambda_{1},\dots,\lambda_{KP}$ and a complete set of eigenfunctions $\boldsymbol{\Psi}_{1},\dots,\boldsymbol{\Psi}_{KP}\in\boldsymbol{\mathcal{H}}$ such that $\boldsymbol{\mathcal{K}}\boldsymbol{\Psi}_{p}=\lambda_{p}\boldsymbol{\Psi}_{p},$ for $p=1,\dots,KP.$ Since $\boldsymbol{\mathcal{K}}$ is a non-negative operator, we know that $\lambda_{p}\geq 0$ for $p=1,\dots,KP$ (Happ and Greven (2018), Proposition 2). From Theorem VI.17 of Reed and Simon (1972), we have that the positive, bounded, self-adjoint, and compact operator $\boldsymbol{\mathcal{K}}$ can be written as $\boldsymbol{\mathcal{K}}\mathbf{f}=\sum_{p=1}^{KP}\lambda_{p}\left\langle\boldsymbol{\Psi}_{p},\mathbf{f}\right\rangle_{\boldsymbol{\mathcal{H}}}\boldsymbol{\Psi}_{p}.$ Thus, from Equation 7, we have that

where $\Psi_{p}^{(k)}(t^{(k)})$ is the $k^{th}$ element of $\boldsymbol{\Psi}_{p}(\mathbf{t})$. Thus, we can see that the covariance kernel can be written as the finite sum of eigenfunctions and eigenvalues,

Since we are working under the assumption that the random function, $\mathbf{f}\in\boldsymbol{\mathcal{H}}$, we can use a modified version of the Multivariate KL theorem. Considering that $\{\boldsymbol{\Psi}_{1},\dots,\boldsymbol{\Psi}_{KP}\}$ form a complete basis for $\boldsymbol{\mathcal{H}}$ and $\mathbf{f}(\mathbf{t})-\boldsymbol{\mu}(\mathbf{t})\in\boldsymbol{\mathcal{H}}$, we have

where $\mathbf{P}_{\boldsymbol{\mathcal{H}}}$ is the projection operator onto $\boldsymbol{\mathcal{H}}$. Letting $\rho_{p}=\left\langle\boldsymbol{\Psi}_{p},\mathbf{f}-\boldsymbol{\mu}\right\rangle_{\boldsymbol{\mathcal{H}}}$, we have $\mathbf{f}(\mathbf{t})-\boldsymbol{\mu}(\mathbf{t})=\sum_{p=1}^{KP}\rho_{p}\boldsymbol{\Psi}_{p}(\mathbf{t})$, where $\mathbb{E}(\rho_{p})=0$ and $\text{Cov}(\rho_{p},\rho_{p^{\prime}})=\lambda_{p}\delta_{pp^{\prime}}$ (Happ and Greven, 2018).

Since $\boldsymbol{\Psi}_{p}\in\boldsymbol{\mathcal{H}}$ and $\boldsymbol{\mu}\in\boldsymbol{\mathcal{H}}$, there exist $\boldsymbol{\nu}_{k}\in\mathbb{R}^{P}$ and $\boldsymbol{\phi}_{kp}\in\mathbb{R}^{P}$ such that $\boldsymbol{\mu}^{(k)}(\mathbf{t})=\boldsymbol{\nu}_{k}^{\prime}\mathbf{B}\left(t^{(k)}\right)$ and $\sqrt{\lambda_{p}}\boldsymbol{\Psi}^{(k)}_{p}(\mathbf{t})=\boldsymbol{\phi}_{kp}^{\prime}\mathbf{B}\left(t^{(k)}\right):=\boldsymbol{\Phi}_{p}^{(k)}(\mathbf{t})$. These scaled eigenfunctions, $\boldsymbol{\Phi}_{p}^{(k)}(\mathbf{t})$, are used over $\boldsymbol{\Psi}^{(k)}_{p}(\mathbf{t})$ because they fully specify the covariance structure of the latent features, as described in Section 2.2. From a modeling prospective, the scaled eigenfunction parameterization is advantageous, as it admits a prior model based on the multiplicative gamma process shrinkage prior proposed by Bhattacharya and Dunson (2011), allowing for adaptive regularized estimation (Shamshoian et al., 2022). Thus we have

where $\chi_{p}=\left\langle(\lambda_{p})^{-1/2}\boldsymbol{\Psi}_{p},\mathbf{f}-\boldsymbol{\mu}\right\rangle_{\boldsymbol{\mathcal{H}}}$, $\mathbb{E}(\chi_{p})=0$, and $\text{Cov}(\chi_{p},\chi_{p^{\prime}})=\delta_{pp^{\prime}}$. Equation 10 gives us a way to jointly decompose any realization of the $K$ stochastic processes, $\mathbf{f}\in\boldsymbol{\mathcal{H}}$, as a finite weighted sum of basis functions.

In this section, we will describe a Bayesian additive model that allows for the constructive representation of mixed memberships for functional data. Our model allows for direct inference on the mean and covariance structures of the $K$ stochastic processes, which are often of scientific interest. To aid computational tractability, we will use the joint decomposition described in Section 2.1. In this section, we will assume that $f^{(k)}\in\boldsymbol{\mathcal{S}}$, and that the conditions of Lemma 2.2 hold.

We aim to model the mixed membership of $N$ observed sample paths, $\{\mathbf{Y}_{i}(.)\}_{i=1}^{N}$, to $K$ latent functional features ${\bf f}=\left(f^{(1)},f^{(2)},\ldots,f^{(K)}\right)$. Assuming that each path is observed at $n_{i}$ evaluation points ${\bf t}_{i}=[t_{i1}\cdots t_{in_{i}}]^{\prime}$, without loss of generality, we define a sampling model for the finite-dimensional marginals of ${\bf Y}_{i}({\bf t}_{i})$. To allow for path-specific partial membership, we extend the finite mixture model in Equation 1 and introduce path-specific mixing proportions $Z_{ik}\in(0,1)$, s.t. $\sum_{k=1}^{K}Z_{ik}=1$, for $i=1,2,\ldots,N$.

Let $\mathbf{S}(\mathbf{t}_{i})=[\mathbf{B}(t_{1})\cdots\mathbf{B}(t_{n_{i}})]\in\mathbb{R}^{P\times n_{i}}$, $\chi_{im}\sim\mathcal{N}(0,1)$, and $\boldsymbol{\Theta}$ denote the collection of model parameters. Let $M\leq KP$ be the number of eigenfunctions used to approximate the covariance structure of the $K$ stochastic processes. Using the decomposition of $f^{(k)}(t)$ in Equation 10 and assuming a normal distribution on ${\bf Y}_{i}({\bf t}_{i})$, we obtain:

If we integrate out the latent $\chi_{im}$ variables, we obtain a more transparent form for the proposed functional mixed membership process. Specifically, we have

where $\boldsymbol{\Theta}_{-\chi}$ is the collection of our model parameters excluding the $\chi_{im}$ variables, and $\mathbf{V}_{i}=\sum_{k=1}^{K}\sum_{k^{\prime}=1}^{K}Z_{ik}Z_{ik^{\prime}}\left\{\mathbf{S}^{\prime}(\mathbf{t}_{i})\sum_{m=1}^{M}\left(\boldsymbol{\phi}_{km}\boldsymbol{\phi}^{\prime}_{k^{\prime}m}\right)\mathbf{S}(\mathbf{t}_{i})\right\}.$ Thus, for a sample path, we have that the mixed membership mean is a convex combination of the functional feature means, and the mixed membership covariance, is a weighted sum of the covariance and cross-covariance functions between different functional features, following from the multivariate KL characterization in Section 2.1. Furthermore, from Equation 9 it is easy to show that, for large enough $M$, we have $\mathbf{S}^{\prime}(\mathbf{t}_{i})\sum_{m=1}^{M}\left(\boldsymbol{\phi}_{kp}\boldsymbol{\phi}^{\prime}_{k^{\prime}p}\right)\mathbf{S}(\mathbf{t}_{i})\approx C^{(k,k^{\prime})}(\mathbf{t}_{i},\mathbf{t}_{i}),$with equality when $M=KP$.

Mixed membership models can be thought of as a generalization of clustering. As such, these stochastic schemes are characterized by an inherent lack of likelihood identifiability. A typical source of non-indentifiability is the common label switching problem. To deal with the label switching problem, a relabelling algorithm can be derived for this model directly from the work of Stephens (2000). A second source of non-identifiabilty stems from allowing $Z_{ik}$ to be continuous random variables. Specifically, consider a model with 2 features and let $\boldsymbol{\Theta}_{0}$ be the set of “true” parameters. Let $Z_{i1}^{*}=0.5(Z_{i1})_{0}$ and $Z_{i2}^{*}=(Z_{i2})_{0}+0.5(Z_{i1})_{0}$. If we let $\boldsymbol{\nu}_{1}^{*}=2(\boldsymbol{\nu}_{1})_{0}-(\boldsymbol{\nu}_{2})_{0}$, $\boldsymbol{\nu}_{2}^{*}=(\boldsymbol{\nu}_{2})_{0}$, $\boldsymbol{\phi}_{1m}^{*}=2(\boldsymbol{\phi}_{1m})_{0}-(\boldsymbol{\phi}_{2m})_{0}$, $\boldsymbol{\phi}_{2m}^{*}=(\boldsymbol{\phi}_{2m})_{0}$, $\chi_{im}^{*}=(\chi_{im})_{0}$, and $(\sigma^{2})^{*}=\sigma^{2}_{0}$, we have $P(Y_{i}(t)|\boldsymbol{\Theta}_{0})=P(Y_{i}(t)|\boldsymbol{\Theta}^{*})$ (Equation 11). Thus, we can see that our model is not identifiable, and we will refer to this problem as the rescaling problem. To address the rescaling problem, additional assumptions such as the separability condition (Papadimitriou et al., 1998; McSherry, 2001; Azar et al., 2001; Chen et al., 2022) or the sufficiently scattered condition (Huang et al., 2016; Jang and Hero, 2019; Chen et al., 2022) are often made on the allocation parameters. The separability condition assumes that at least one observation belongs entirely to each feature, making it relatively easy to incorporate into a sampling model. To ensure that the separability condition holds in a two-feature model, we use the Membership Rescale Algorithm (Algorithm 2 in the Supplementary Materials). Compared to the separability condition, the sufficiently scattered condition makes weaker geometric assumptions on the sampling model when $K\geq 3$, but it is more difficult to implement in a sampling scheme. Section C.4 in the Supplementary Materials has an in-depth review of the two common assumptions made to address the rescaling problem.

A third non-identifiability problem may arise numerically as a form of concurvity, i.e. when $\boldsymbol{\nu}_{k^{\prime}}\propto\boldsymbol{\phi}_{km}$ in Equation 11. Typically, an overestimation of the magnitude of $\boldsymbol{\phi}_{km}$, may result in small variance estimates for the allocation parameters (smaller credible intervals). To address this problem, we can assume the separability condition; however, under weaker geometric assumptions this problem may persist, but is typically of little practical relevance.

The sampling model in Section 2.2, allows a practitioner to select how many eigenfunctions are to be used in the approximation of the covariance function. In the case of $M=KP$, we have a fully saturated model and can represent any realization $\mathbf{f}\in\boldsymbol{\mathcal{H}}$. In Equation 10, $\boldsymbol{\Phi}$ parameters are mutually orthogonal (where orthogonality is defined by the inner product defined in Equation 6), and have a magnitude proportional to the square root of the corresponding eigenvalue, $\lambda_{p}$. Thus, a modified version of the multiplicative gamma process shrinkage prior proposed by Bhattacharya and Dunson (2011) will be used as our prior for $\boldsymbol{\phi}_{km}$. By using this prior, we promote shrinkage across the $\boldsymbol{\phi}_{km}$ coefficient vectors, with increasing prior shrinkage towards zero as $m$ increases.

To facilitate MCMC sampling from the posterior target we will remove the assumption that $\boldsymbol{\Phi}$ parameters are mutually orthogonal. Even though $\boldsymbol{\Phi}$ parameters can no longer be thought of as scaled eigenfunctions, posterior inference can still be conducted on the eigenfunctions by post-processing posterior samples of $\boldsymbol{\phi}_{km}$ (given that we can recover the true covariance operator). Specifically, given posterior samples from $\boldsymbol{\phi}_{km}$, we obtain posterior realizations for the covariance function and then calculate the eigenpairs of the posterior samples of the covariance operator using the method described in Happ and Greven (2018). Thus, letting $\phi_{kpm}$ be the $p^{th}$ element of $\boldsymbol{\phi}_{km}$, we have

where $1\leq k\leq K$, $1\leq p\leq P$, $1\leq m\leq M$, and $2\leq j\leq M$. To promote shrinkage across the $M$ matrices, we need $\delta_{jk}>1$. Thus, letting $\alpha_{2}>\beta_{2}$, we have $\mathbb{E}(\delta_{jk})>1$, which will promote shrinkage. In this construction, we allow for different rates of shrinkage across different functional features, which is particularly important in cases where the covariance functions of different features have different magnitudes. In cases where the magnitudes of the covariance functions are different, we would expect $\delta_{mk}$ to be relatively smaller in the $k$ associated with the functional feature with a large covariance function.

To promote adaptive smoothing in the mean function, we will use a first-order random walk penalty proposed by Lang and Brezger (2004). The first-order random walk penalty penalizes differences in adjacent B-spline coefficients. In the case where $\mathcal{T}\subset\mathbb{R}$, we have that $P(\boldsymbol{\nu}_{k}|\tau_{k})\propto exp\left(-\frac{\tau_{k}}{2}\sum_{p=1}^{P-1}\left(\nu_{pk}^{\prime}-{\nu}_{(p+1)k}\right)^{2}\right),$ for $k=1,\dots,K$, where $\tau_{k}\sim\Gamma(\alpha,\beta)$ and $\nu_{pk}$ is the $p^{th}$ element of $\boldsymbol{\nu}_{k}$. Since we have $Z_{ik}\in(0,1)$ and $\sum_{k=1}^{K}Z_{ik}=1$, it is natural to consider prior Dirichlet sampling for $\mathbf{z}_{i}=[Z_{i1}\cdots Z_{iK}]$. Therefore, following Heller et al. (2008), we have

for $i=1,\dots,N$. Lastly, we will use a conjugate prior for our random error component of the model, such that $\sigma^{2}\sim IG(\alpha_{0},\beta_{0})$. Although we relax the assumption of orthogonality, in Section C.1 of the Supplementary Materials, we outline an alternative sampling scheme where we impose the condition that the $\boldsymbol{\Phi}$ parameters form orthogonal eigenfunctions using the work of Kowal et al. (2017).

In the previous section, we saw that the $\boldsymbol{\Phi}$ parameters are not assumed to be mutually orthogonal. By removing this constraint, we facilitate MCMC sampling from the posterior target and can perform inference on the eigenpairs of the covariance operator, as long as we can recover the covariance structure. Due to the complex identifiability issues mentioned in Section 2.2, establishing weak posterior consistency with unknown allocation parameters unattainable, even if we include the orthogonality constraint on the $\boldsymbol{\Phi}$ parameters. However, the model becomes identifiable when we condition on the allocation parameters. In this section, we show that we can achieve conditional weak posterior consistency without enforcing the orthogonality constraint of the $\boldsymbol{\Phi}$ parameter. Therefore, we will show that conditional on the allocation parameters, relaxing the orthogonality constraint does not affect our ability to recover the mean and covariance structure of the $K$ underlying stochastic processes from Equation 12. Let $\boldsymbol{\Pi}$ be the prior distribution on $\boldsymbol{\omega}:=\left\{\boldsymbol{\nu}_{1},\dots,\boldsymbol{\nu}_{K},\boldsymbol{\Sigma}_{11},\dots,\boldsymbol{\Sigma}_{1K},\dots,\boldsymbol{\Sigma}_{KK},\sigma^{2}\right\}$, where $\boldsymbol{\Sigma}_{kk^{\prime}}:=\sum_{p=1}^{KP}\left(\boldsymbol{\phi}_{kp}\boldsymbol{\phi}^{\prime}_{k^{\prime}p}\right)$. We will be proving weak posterior consistency with respect to the parameters $\boldsymbol{\Sigma}_{kk^{\prime}}$ because the parameters $\boldsymbol{\phi}_{kp}$ are non-identifiable. Since the covariance and cross-covariance structure are completely specified by the $\boldsymbol{\Sigma}_{kk^{\prime}}$ parameters, the lack of identifiability of the $\boldsymbol{\phi}_{kp}$ parameters bears no importance on inferential considerations. We will denote the true set of parameter values as $\boldsymbol{\omega}_{0}=\left\{(\boldsymbol{\nu}_{1})_{0},\dots,(\boldsymbol{\nu}_{K})_{0},(\boldsymbol{\Sigma}_{11})_{0},\dots,(\boldsymbol{\Sigma}_{1K})_{0},\dots,(\boldsymbol{\Sigma}_{KK})_{0},\sigma^{2}_{0}\right\}.$ To prove weak posterior consistency, we will make the following assumptions:

In order to prove weak posterior consistency, we will first specify the following quantities related to the Kullback–Leibler (KL) divergence. Following the notation of Choi and Schervish (2007), we define the following quantities:

where $f_{i}(\mathbf{Y}_{i};\boldsymbol{\omega}_{0})$ is the likelihood under $\boldsymbol{\omega}_{0}$. To simplify the notation, we define the following two quantities:

where $\mathbf{U}_{i}^{\prime}\mathbf{D}_{i}\mathbf{U}_{i}$ is the corresponding spectral decomposition. Let $d_{il}$ be the $l^{th}$ diagonal element of $\mathbf{D}_{i}$. Let $\boldsymbol{\Omega}_{\epsilon}(\boldsymbol{\omega}_{0})$ be the set of parameters such that the KL divergence is less than some $\epsilon>0$ ($\boldsymbol{\Omega}_{\epsilon}(\boldsymbol{\omega}_{0}):=\{\boldsymbol{\omega}:K_{i}(\boldsymbol{\omega}_{0},\boldsymbol{\omega})<\epsilon\text{ for all i}\}$). Let $a,b\in\mathbb{R}$ be such that $a>1$ and $b>0$, and define $\mathcal{B}(\boldsymbol{\omega}_{0}):=\left\{\boldsymbol{\omega}:\frac{1}{a}\left((d_{il})_{0}+\sigma^{2}_{0}\right)\leq d_{il}+\sigma^{2}\leq a\left((d_{il})_{0}+\sigma^{2}_{0}\right),\|\left(\boldsymbol{\mu}_{i}\right)_{0}-\boldsymbol{\mu}_{i}\|\leq b\right\}.$

Lemma 3.1 shows that the prior probability of our parameters being arbitrarily close (defined by KL divergence) to the true parameters is positive. Since $\mathbf{Y}_{1},\dots,\mathbf{Y}_{N}$ are not identically distributed, condition (1) in Lemma 3.1 is necessary to prove Lemma 3.2.

All proofs are provided in the Supplementary Materials. Lemma 3.2 shows us that, conditional on the allocation parameters, we are able to recover the covariance structure of the $K$ stochastic processes. Thus, relaxation of the orthogonality constraint on the $\boldsymbol{\phi}_{km}$ parameters does not affect our ability to perform posterior inference on the main functions of scientific interest. Inference on the eigenstructure can still be performed by calculating the eigenvalues and eigenfunctions of the covariance operator using the MCMC samples of the $\boldsymbol{\phi}_{km}$ parameter. Finally, we point out that, in most cases, the parameters $Z_{ik}$ are unknown. Although theoretical guarantees for consistent estimation of the latent mixed allocation parameters are still elusive, we provide some empirical evidence of convergence in Section 4.1.

In this simulation study, we examine how well our model can recover the mean and covariance functions when we vary the number of observed functions. The model used in this section is a mixed membership model with $2$ functional features ($K=2$), which can be represented by a basis constructed of $8$ B-splines ($P=8$), and uses $3$ eigenfunctions ($M=3$). We consider the case where we have $40$, $80$, and $160$ observed functional observations, where each observation is uniformly observed at $100$ time points. We then simulated 50 data sets for each of the three cases ($N=40,80,160$). To measure how well we recovered the the mean, covariance, and cross-covariance functions, we estimate the integrated relative mean square error (R-MISE), where $\text{R-MISE}=\frac{\int\{f(t)-\hat{f}(t)\}^{2}\text{d}t}{\int f(t)^{2}\text{d}t}\times 100\%$ and $\hat{f}$ is the estimated posterior median of the function $f$. To measure how well we recovered the allocation parameters, $Z_{ik}$, we calculated the root-mean-square error (RMSE). Section B.1 of the Supplementary Materials goes into more detail of how the model parameters were drawn, how estimates of all quantities of interest were calculated in this simulation, and includes visualizations of some of the recovered covariance structures.

From Figure 2, we can see that we have a good recovery of the mean structure with as few as $40$ functional observations. While the R-MISE of the mean functions improves as we increase the number of functional observations ($N$), this improvement will likely have little practical impact. However, when looking at the recovery of the covariance and cross-covariance functions, we can see that the R-MISE noticeably decreases as more functional observations are added. As the recovery of the mean and covariance structures improves, the recovery of the allocation structure $(\mathbf{Z})$ improves. A supplementary simulation study examining how well our model can recover the mean, covariance, and allocation structures in a three-feature model ($K=3$) can be found in Section B.1 of the Supplementary Materials.

Choosing the number of latent features can be challenging, especially when no prior knowledge is available for this quantity. Information criteria, such as the AIC, BIC, or DIC, are often used to aid practitioners in the selection of a data-supported value for $K$. In this simulation, we evaluate how various types of information criteria, as well as simple heuristics such as the “elbow-method”, perform in recovering the true number of latent features. To do this, we simulate 10 different data sets, each with 200 functional observations, from a mixed membership model with three functional features. We then calculate these information criteria on the 10 data sets for mixed membership models where $K=2,3,4,5$. Additional information on how the simulation study was conducted, as well as definitions of the information criterion, can be found in Section B.2 of the Supplementary Materials. As specified, the optimal model should have the largest BIC, the smallest AIC, and the smallest DIC.

From the simulation results presented in Figure 3, we can see that on average each information criterion overestimated the number of functional features in the model. While the three information criteria seem to greatly penalize models that do not have enough features, they do not seem to punish overfit models to a great enough extent. As expected, the log-likelihood increases as we add more features; however, we can see that there is often an elbow at $K=3$ (Figure 3). Using the “elbow-method” led to selecting the correct number of latent functional features 8 times out of 10, while BIC picked the correct number of latent functional features twice. DIC and AIC were found to be the least reliable information criteria, only choosing the correct number of functional features once. Thus, through empirical consideration, we suggest using the “elbow-method” along with the BIC to aid in the final selection of the number of latent features to be interpreted in analyses. Although formal considerations of model-selection consistency are out of scope for the current contribution, we maintain that some of these techniques are best interpreted in the context of data exploration, with a potential for great improvement in interpretation if semi-supervised considerations allow for an a-priori informed choice of $K$.

Autism spectrum disorder (ASD) is a term used to describe individuals with a collection of social communication deficits and restricted or repetitive sensory-motor behaviors (Lord et al., 2018). While once considered a very rare disorder with very specific symptoms, today the definition is more broad and is now thought of as a spectrum. Some individuals with ASD may have minor symptoms, while others may have very severe symptoms and require lifelong support. In this case study, we will use electroencephalogram (EEG) data that was previously analyzed by Scheffler et al. (2019) in the context of regression. The data set consists of EEG recordings of 39 typically developing (TD) children and 58 children with ASD between the ages of 2 and 12 years of age, which were analyzed in the frequency domain.

Scheffler et al. (2019) observed that the data from the T8 electrode, corresponding to the right temporal region, were correlated with ASD diagnosis, so we will specifically use the data from the T8 electrode in our mixed membership model. We focus our analysis to the alpha band of frequencies (6 to 14 Hz), whose patterns at rest are thought to play a role in neural coordination and communication between distributed brain regions. As clinicians examine the sample paths for the two cohorts, shown in Figure 4 (right panel), they are often interested in the location of a single prominent peak in the spectral density located within the alpha frequency band called the peak alpha frequency (PAF). This quantity has previously been associated with neural development in TD children (Rodríguez-Martínez et al., 2017). Scheffler et al. (2019) found that as TD children age, the peak becomes more prominent and the PAF shifts to a higher frequency. Conversely, a discernible PAF pattern is attenuated for most children with ASD compared to their TD counterpart.

A visual examination of the data in Figure 4 (right panel) anticipates the potential inadequacy of cluster analysis in this application, as path-specific heterogeneity does not seem to define well-separated subpopulations. In fact, if we cluster our data using the model in Pya Arnqvist et al. (2021) with $K=4$ (BIC selection), we find cluster means of dubious interpretability and poor separation of sample paths between clusters.

In contrast to classical clustering, we use a mixed membership model with only two functional features ($K=2$), (AIC-BIC selection). We note that the enhanced flexibility of mixed membership models induces parsimony in the number of pure mixture components supported by the data. In particular, the mean function of the first feature, depicted in Figure 5, can be interpreted as $1/f$ noise or pink noise. This component noise is expected to be found in every individual to some extent, but we can see that the first feature has no discernible alpha peak. The mean function of the second feature captures a distinct alpha peak, which is typically observed in EEGs of older neurotypical individuals. These two features help to differentiate between periodic (alpha waves) and aperiodic ($1/f$ trend) neural activity patterns, which coexist in the EEG spectrum. In this context, a model of uncertain membership would be necessarily inadequate to describe the heterogeneity of the observed sample path, as we would not naturally think of subjects in our sample as belonging to one or the other cluster. Instead, assuming mixed membership between the two-feature processes is likely to represent a more realistic data generating mechanism, as we conceptualize periodic and aperiodic neural activity patterns to mix continuously.

Figure 5 shows that TD children are highly likely to heavily load feature 2 (well-defined PAF), whereas ASD children exhibit a higher level of heterogeneity. The cohort average (feature-1) membership, indicated by the triangles in Figure 5, shows that the average (feature-1)-membership for children with ASD is 0.42 compared to 0.25 for TD children. It is important to note that the allocation parameters constitute compositional data, meaning that using a Euclidean metric to measure mean membership may be inappropriate in cases where we have more than two features. Therefore, the use of Fréchet means (Fréchet, 1948) with a suitable metric for compositional data may be more appropriate to measure the centroids of the allocation parameters when we have more than two functional features. In general, these loadings suggest that clear alpha peaks take longer to emerge in children with ASD, compared to their typically developing counterparts.

In general, our findings confirm related evidence in the scientific literature on developmental neuroimaging but offer a completely novel point of view in quantifying group membership as part of a spectrum. Additional details on how the study was conducted, a posterior predictive check for a subset of individuals, and an in-depth comparison between our method and alternative methods such as FPCA and functional clustering are reported in Section B.3 of the Supplementary Materials. An extended analysis of multi-channel EEG for the same case study is reported in Section B.4 of the Supplementary Materials, investigating how spatial patterns vary across the scalp.