In alexanderjwhite/hthmixture: High Dimensional Multivariate Mixture Regression

library(dplyr)
knitr::opts_chunk$set(echo = FALSE)

Motivation

Understanding how DNA methylation affects/predicts gene expression in heterogeneous patient populations.

Suppose for some patient population we have:

• DNA methylation levels for some set of CpG sites
• Gene expression levels for some clinically relevant sites.

We are interested in:

1. Clustering patients into the correct sub-group
2. Within sub-group, finding relationship between CpG site & gene expression

Idea

Goal: Correctly cluster observations & regress in high dimensional $X$ & $Y$.

• $Y_{n \times q}$ Matrix of responses
• $X_{n \times p}$ Design matrix
• $A_{p \times q}$ Coefficient matrix (sparse in p)
• $k$ clusters
• Assume (for now) response covariance ($\Sigma_k$) is diagonal
• Strong yet common assumption due to lack of degrees of freedom for correct variance specification

$$ f\left(\boldsymbol{y}{i} \mid \boldsymbol{x}{\boldsymbol{i}} ; \boldsymbol{\theta}\right)=\Sigma_{k=1}^{K} \pi_{k} \mathcal{N}{q}\left(\boldsymbol{y}{\boldsymbol{i}} ; \boldsymbol{x}{\boldsymbol{i}} A{k}, \Sigma_{k}\right) $$

• Parameter space $\theta={\pi_k,A_k,\Sigma_k;k=1\ldots K }$ solved by general EM using SARRS to compute $A_k$.

Multivariate Regression

• Used when there are multiple response variables


Mixed Regression

• Used when there are multiple response variables


Reduced Rank Regression

• Coefficient matrix sparse in q dimension
• Only few genes are associated with methylation


HTH Mixture Algorithm

• Initialize: $\pi_k^{(0)}=\frac{n_k^{(0)}}{n}$
• Randomly initialize observations into k clusters

While not converged ($m=1,\ldots, M$) do:

• for $k=1,\ldots,K$ apply SARRS on all observations in $C_k^{(m-1)}$ to obtain $A_k^{(m)}$, $\Sigma_k^{(m)}$
• compute $\mu_{ik}^{(m)} = \mathcal{N}{p}\left(\boldsymbol{y}{\boldsymbol{i}} ; A_{k}^{(m)} \boldsymbol{x}{\boldsymbol{i}}, \Sigma{k}^{(m)}\right)$
• $C_k^{(m)}={i|\text{ML component k} }$
• $\pi_k^{(m)}=\frac{n_k^{(m)}}{n}$

HTH Mixture Algorithm

• Empirically, HTH Mixture reaches local maximum quickly
• Need to determine method for efficiently exploring the likelihood space with different random initialization.

Likelihood Space Exploration Idea 1

• Mimic the idea of chains from MCMC
• In parallel, run large number of chains and select the one which reaches the largest likelihood.
• This method works well, although not the most efficient way to explore the space.

Likelihood Space Exploration Idea 2

• To more efficiently explore the space, when a local maximum is reached, perturb the initialized state by a proportion $p$ (e.g. $20\%$) and rerun.

Data Simulation

• $X_k$ consists of iid random vectors sample from $MVN(\mathbf{0},\Sigma_k)$
• $\Sigma_k$ diagonal
• Noise matrix $Z_k \in \mathbb{R}^{n \times q}$ has iid $N(0,\sigma^2)$ entries
• $A_k=\left(\begin{array}{c}b_k B_{0_k} B_{1_k} \0\end{array}\right)$
• with $b>0$, $B_0 \in \mathbb{R}^{s \times r}$, $B_1 \in \mathbb{R}^{r \times q}$
• $Y_k = X_kA_k+Z_k$

Finally, combine $X$ & $Y$

Simulation

• N: $50 - 400$
• s: $5 - 20$
• rank: $1 - 2$
• B: $1$
• $p=m=1000$
• Perturb: $20\%$
• Cluster with equal probability
• K: $2$
• $\sigma$: $1$
• CHains: 500

Results

• HTH algorithm consistently reaches local maximum

Results

Results

• Idea 1 $\rightarrow$ random initialization each time

Results

• Idea 2 $\rightarrow$ perturb local maximum
• Appears to be more efficient

Results

data <- readr::read_rds("G:/My Drive/Dissertation/HTH Mixture/hthmixture/results/20210415_slide_results.rds")


data %>% 
  slice(1:12) %>% 
  knitr::kable(format = "latex", digits = 2)

Results

data %>% 
  slice(13:24) %>% 
  knitr::kable(format = "latex", digits = 2)

Performance

• In simulated data, current algorithm clusters well (perfectly in most cases given enough chains)

• Challenges:

• Likelihood space is explored inefficiently, sometimes requires large number of chains
• Computationally cumbersome/inefficient

alexanderjwhite/hthmixture documentation built on Sept. 4, 2022, 4:48 a.m.