Software design, selection and estimation for latent variable models

background-color: #98ebee
class: middle center

<div class="shade_black" style="width:calc(45%);right:0;bottom:0;padding-left:5px;border: dashed 4px white;margin: auto;">
 These slides are viewed best by Chrome and occasionally need to be refreshed if elements did not load properly. See here for <a href=https://www.dropbox.com/s/y2lwp798th6t9y1/slides_ofal.pdf?dl=0 style="color:black!important"/>PDF </a>. 
</div>

---

background-image: url(images/bg3.jpg)
background-size: cover
class: hide-slide-number
count: false

:::::::::: { .grid grid: 1fr / 7fr 3fr;}

::: {.item .shade_black border-right-style: solid; border-right-color: white;}

# Software design, selection and estimation for latent variable models

##

.bottom_abs.width100[

Presented by Emi Tanaka

School of Mathematics and Statistics 
 dr.emi.tanaka@gmail.com
 @statsgen

29th Nov 2019 @ WOMBAT2019 | Melbourne, Australia
]

::: item

:::

::::::::::

---

# Motivating Data

---

# .circle[1] Marketing: Product testing survey

---

# .circle[2] Psychometrics: Human intelligence test

---

# .circle[3] Bioinformatics: Gene expression data

---

# .circle[4] Ecology: Multi-abundance species data

---

# .circle[5] Agriculture: Multi-environmental field trial

---

# Multi-variate data .font_small[(where the response variables span multiple columns)]

::: paddings5

* The unit .pink[(row)] response is measured for different variables/traits/attributes .purple[(column)]
* This data structure is common in many disciplines

:::

::: { .grid grid: 1fr / 1fr 1fr; }

::: { .item font-size: 0.7em; .bg-grey }

:::

::: { .item font-size: 0.8em; .bg-grey }

.circle-list[1] **Marketing**: Product testing survey   
&emsp; .pink[Person] vs .purple[product attributes]  
.circle-list[2] **Psychometrics**: Human intelligence test   
&emsp; .pink[Student] vs .purple[intellect test]  
.circle-list[3] **Bioinformatics**: Gene expression data   
&emsp; .pink[Sample] vs .purple[genes]  
.circle-list[4] **Ecology**: Multi-abundance species data   
&emsp; .pink[Site] vs .purple[species]  
.circle-list[5] **Agriculture**: Multi-environmental field trial   
&emsp; .pink[Variety] vs .purple[site]

:::

---

# 🎯 Aim and assumptions

::: { .paddings margin-top: 10px; }

* The typical aim of analysis is to fit a reasonable model to data for inference or prediction
* The observations from a .pink[row] unit would likely be **correlated** across .purple[column]
* Assume responses across .pink[row] are independent for this talk (but in theory does not have to be)

:::

---

# Latent variable models

can model this reasonably well

---

# Latent variable models .font_small[(or **factor analytic models**)]

::: { .paddings margin-top: 10px; }

::: grid

::: { .info-box margin-bottom: 10px; width: 370px; margin-right: 10px; }

Non-Gaussian Distribution

`$$g(\mu_{ij}) = \eta_{ij} = \boldsymbol{x}_i^\top\boldsymbol{\beta}_j + \boldsymbol{u}_i^\top\color{blue}{\boldsymbol{\lambda}_j}$$`

:::

::: { .info-box  width: 370px; }

Gaussian Distribution

`$$y_{ij} =  \boldsymbol{x}_i^\top\boldsymbol{\beta}_j + \boldsymbol{u}_i^\top  \color{blue}{ \boldsymbol{\lambda}_j} + \epsilon_{ij}$$`

:::

.font_small[
<details>
<ul>
<li> $g(\cdot)$ is a known link function</li>
<li> $\boldsymbol{\beta}_j$ is a $q$-vector of fixed effects associated with covariates $\boldsymbol{x}_i$</li>
<li> $\boldsymbol{\lambda}_j$ is a $k$-vector of loadings for response $j$ </li> 
</ul>
Let $\boldsymbol{\eta}_i = (\eta_{i1}, ..., \eta_{ip})^\top$ and $\boldsymbol{\epsilon}_i = (\epsilon_{i1}, ..., \epsilon_{ip})^\top$
<ul>
<li>$\mathbf{\Lambda} = \begin{bmatrix}\boldsymbol{\lambda}_1\cdots \boldsymbol{\lambda}_n\end{bmatrix}^\top$</li>
<li>$\mathbf{\Psi}$ is a diagonal matrix with diagonal entries $\psi_i$</li>
<li> $\boldsymbol{u}_i \sim N(\boldsymbol{0}_k, \mathbf{I}_k)$, $Cov(\boldsymbol{\eta}_i|\boldsymbol{x}_i) = \mathbf{\Lambda}\mathbf{\Lambda}^\top$, $\boldsymbol{\epsilon}_i \sim N(\boldsymbol{0}_p, \mathbf{\Psi})$ and $Cov(\boldsymbol{y}_i | \boldsymbol{x}_i) = \mathbf{\Lambda}\mathbf{\Lambda}^\top + \mathbf{\Psi}$ </li>
</ul>
</details>
]

::: { .grid grid: 1fr / 2fr 1fr; }

::: item

**Main idea**

* `$p$`-vector response of .pink[row] unit is a result of `$k$` (independent) latent variables where `$k < p$`
* this results in **dimension reduction**

:::

::: item center

:::

---

# But how many latent variables?

(in another words what should `$k$` be?)

---

# Estimating variance parameters

::: paddings

* Note that under our model for given `$k$`

`$$Cov(\boldsymbol{\eta}_i|~\boldsymbol{x}_i) = \mathbf{\Lambda}\mathbf{\Lambda}^\top\quad\text{or}\quad Cov(\boldsymbol{y}_i|~\boldsymbol{x}_i) = \mathbf{\Lambda}\mathbf{\Lambda}^\top + \mathbf{\Psi}$$`

* So `$\mathbf{\Lambda}$` is a set of variance parameters

* Let `$\boldsymbol{\theta}$` denote a vector of all variance parameters in a model then

::: { .info-box width: 60%; margin: auto; font-size: 1.2em; }

`$$\DeclareMathOperator*{\argmax}{arg\,max}\hat{\boldsymbol{\theta}}_{\text{ML/REML}} = \argmax_{\boldsymbol{\theta}}~ \ell(\boldsymbol{\theta}| ~\cdot)$$`

:::

* `$k$` is chosen using likelihood ratio tests or some ad hoc approaches

:::

---

# Ordered Factor LASSO (OFAL)

::: paddings

* Assume that `$\mathbf{\Lambda}_0$` is a `$p\times d$` pseudo factor loading matrix where `$k \leq d \leq p$`

* .blue2[**OFAL estimate**]: our approach via penalised likelihood

:::

`$$\displaystyle\hat{\boldsymbol{\theta}}_{\text{OFAL}}  = \argmax_{\boldsymbol{\theta}} \left\{\ell(\boldsymbol{\theta}) - s \sum _{l=1}^d \omega_{g,l} \sqrt{ \sum_{i=1}^p\sum_{j=l}^d \lambda_{0,ij}^2} - s \sum_{i=1}^p\sum_{j=1}^d \omega_{e,ij}|\lambda_{0,ij}|\right\}$$`

::: { .paddings font-size: 0.8em; }

where 
* `$s$` is a **tuning parameter**,
* `$\omega_{g,l}$` is a **group-wise adaptive weight** for `$l$`th column of `$\mathbf{\Lambda}_0$`, and
* `$\omega_{e,ij}$` is an **element-wise adaptive weight** for `$i,j$`th entry of `$\mathbf{\Lambda}_0$`.

:::

::: { .footnote font-size: 0.63em; }

Hui, Tanaka & Warton (2018) Order Selection and Sparsity in Latent Variable Models via the Ordered Factor LASSO. Biometrics

:::

---

# OFAL concept demo: element-wise sparsity

::: { .grid grid: 1fr / 1fr 900px; }

::: item

`$$\mathbf{\Lambda}_0$$`

:::

::: item

Say `$s\omega_{e,15} \rightarrow \infty$`, then you would expect `$|\lambda_{15}| \rightarrow 0$`.

:::

---

# OFAL concept demo: group-wise sparsity

::: { .grid grid: 1fr / 1fr 900px; }

::: item

`$$\mathbf{\Lambda}_0$$`

:::

::: {.item font-size: 0.8em; }

* Say `$s\omega_{g,5} \rightarrow \infty$`, then you would expect `$\left(\sum_{l=1}^p(\lambda_{l5}^2 + \lambda^2_{l6})\right)^{1/2} \rightarrow 0$`.

* Sum of squares is zero only if each element is zero, so `$\lambda_{l5}\rightarrow 0$` and `$\lambda_{l6}\rightarrow 0$` for `$l=1, ..., p$`.

:::

---

# Algorithm in action 🎞

::: grid

::: item

&emsp; `$\mathbf{\Lambda}_0$`

:::

::: item

* Our procedure uses EM-algorithm
* At the M-step, we reformulate the OFAL penalised likelihood into an elastic-net type regularisation problem
* Then employ coordinate-wise optimisation to obtain loading estimates
* Full algorithm (including code), choice of adapative weights and tuning parameter selection are in Hui, Tanaka & Warton (2018)

:::

::: { .footnote font-size: 0.63em; }

Hui, Tanaka & Warton (2018) Order Selection and Sparsity in Latent Variable Models via the Ordered Factor LASSO. Biometrics

:::

---

Well now how to get people to use your method?

# Software designing for practitioners

---

# Specifying the model 🎬 .circle2[1]

::: { .paddings margin-top: 10px; }

::: grid

::: { .info-box margin-bottom: 10px; width: 370px; margin-right: 10px; }

Non-Gaussian Distribution

`$$\color{green}{g}(\mu_{ij}) = \eta_{ij} = \color{purple}{\boldsymbol{x}_i}^\top\boldsymbol{\beta}_j + \boldsymbol{u}_i^\top\boldsymbol{\lambda}_j$$`

:::

::: { .info-box  width: 370px; }

Gaussian Distribution

`$$y_{ij} =  \color{purple}{\boldsymbol{x}_i}^\top\boldsymbol{\beta}_j + \boldsymbol{u}_i^\top  \boldsymbol{\lambda}_j + \epsilon_{ij}$$`

:::

:::

---

# Unbalanced data

What if the cells contain `$>1$` observations or no observation?

---

Users: "Let's make it rectangular [to fit the software]"

::: { .info-box width: 50%; margin: auto; }

**Bottomline**: software design encourages certain user behaviour whether that be good or not

:::

---

# Multivariate to univariate data

---

# Specifying the model 🎬 .circle2[2]

---

# Specifying the model 🎬 .circle2[3]

<center>
Symbolic model formulae

.footnote[
Tanaka and Hui (2019) Symbolic formulae for linear mixed models. https://arxiv.org/abs/1911.08628
]

---

# Software extensibility and generalisability

::: paddings

* Practitioners normally need to fit many kinds of models 
* If software is restrictive, user beahviour is restrictive (for good or bad)

:::

.footnote[
Tanaka and Hui (2019) Symbolic formulae for linear mixed models. https://arxiv.org/abs/1911.08628
]

---

background-image: url("images/bg3.jpg")
background-size: cover
class: hide-slide-number animated pulse fast

:::::::::: { .grid  .white grid: 1fr / 3fr 1fr;}

::: {.item .shade_black  border-right-style: solid; border-right-color: white;}

<h1>Thanks!</h1>

These slides are made using `xaringan` R-package and can be found at

Emi Tanaka

dr.emi.tanaka@gmail.com 
 @statsgen

:::

:::

---

# Session Information

:::: scroll-350

```
  ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   setting  value                       
   version  R version 3.6.0 (2019-04-26)
   os       macOS Mojave 10.14.6        
   system   x86_64, darwin15.6.0        
   ui       X11                         
   language (EN)                        
   collate  en_AU.UTF-8                 
   ctype    en_AU.UTF-8                 
   tz       Australia/Melbourne         
   date     2019-11-29                  
  
  ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   package     * version    date       lib source                            
   anicon        0.1.0      2019-05-28 [1] Github (emitanaka/anicon@377aece) 
   assertthat    0.2.1      2019-03-21 [1] CRAN (R 3.6.0)                    
   backports     1.1.4      2019-04-10 [1] CRAN (R 3.6.0)                    
   broom         0.5.2      2019-04-07 [1] CRAN (R 3.6.0)                    
   callr         3.3.1      2019-07-18 [1] CRAN (R 3.6.0)                    
   cellranger    1.1.0      2016-07-27 [1] CRAN (R 3.6.0)                    
   cli           1.1.0      2019-03-19 [1] CRAN (R 3.6.0)                    
   colorspace    1.4-1      2019-03-18 [1] CRAN (R 3.6.0)                    
   crayon        1.3.4      2017-09-16 [1] CRAN (R 3.6.0)                    
   curl          4.2        2019-09-24 [1] CRAN (R 3.6.0)                    
   desc          1.2.0      2018-05-01 [1] CRAN (R 3.6.0)                    
   devtools      2.0.2      2019-04-08 [1] CRAN (R 3.6.0)                    
   digest        0.6.22     2019-10-21 [1] CRAN (R 3.6.0)                    
   dplyr       * 0.8.3      2019-07-04 [1] CRAN (R 3.6.0)                    
   emo           0.0.0.9000 2019-06-03 [1] Github (hadley/emo@02a5206)       
   evaluate      0.14       2019-05-28 [1] CRAN (R 3.6.0)                    
   forcats     * 0.4.0      2019-02-17 [1] CRAN (R 3.6.0)                    
   fs            1.3.1      2019-05-06 [1] CRAN (R 3.6.0)                    
   generics      0.0.2      2018-11-29 [1] CRAN (R 3.6.0)                    
   ggplot2     * 3.2.1      2019-08-10 [1] CRAN (R 3.6.0)                    
   glue          1.3.1.9000 2019-10-24 [1] Github (tidyverse/glue@71eeddf)   
   gtable        0.3.0      2019-03-25 [1] CRAN (R 3.6.0)                    
   haven         2.1.0      2019-02-19 [1] CRAN (R 3.6.0)                    
   hms           0.5.1      2019-08-23 [1] CRAN (R 3.6.0)                    
   htmltools     0.4.0      2019-10-04 [1] CRAN (R 3.6.0)                    
   httr          1.4.1      2019-08-05 [1] CRAN (R 3.6.0)                    
   icon          0.1.0      2019-05-28 [1] Github (ropenscilabs/icon@a510f88)
   jsonlite      1.6        2018-12-07 [1] CRAN (R 3.6.0)                    
   knitr         1.25       2019-09-18 [1] CRAN (R 3.6.0)                    
   lattice       0.20-38    2018-11-04 [1] CRAN (R 3.6.0)                    
   lazyeval      0.2.2      2019-03-15 [1] CRAN (R 3.6.0)                    
   lifecycle     0.1.0      2019-08-01 [1] CRAN (R 3.6.0)                    
   lubridate     1.7.4      2018-04-11 [1] CRAN (R 3.6.0)                    
   magrittr      1.5        2014-11-22 [1] CRAN (R 3.6.0)                    
   memoise       1.1.0      2017-04-21 [1] CRAN (R 3.6.0)                    
   modelr        0.1.4      2019-02-18 [1] CRAN (R 3.6.0)                    
   munsell       0.5.0      2018-06-12 [1] CRAN (R 3.6.0)                    
   nlme          3.1-140    2019-05-12 [1] CRAN (R 3.6.0)                    
   pillar        1.4.2      2019-06-29 [1] CRAN (R 3.6.0)                    
   pkgbuild      1.0.3      2019-03-20 [1] CRAN (R 3.6.0)                    
   pkgconfig     2.0.3      2019-09-22 [1] CRAN (R 3.6.0)                    
   pkgload       1.0.2      2018-10-29 [1] CRAN (R 3.6.0)                    
   prettyunits   1.0.2      2015-07-13 [1] CRAN (R 3.6.0)                    
   processx      3.4.1      2019-07-18 [1] CRAN (R 3.6.0)                    
   ps            1.3.0      2018-12-21 [1] CRAN (R 3.6.0)                    
   purrr       * 0.3.2      2019-03-15 [1] CRAN (R 3.6.0)                    
   quoter        0.1.0      2019-07-28 [1] local                             
   R6            2.4.0      2019-02-14 [1] CRAN (R 3.6.0)                    
   Rcpp          1.0.2      2019-07-25 [1] CRAN (R 3.6.0)                    
   readr       * 1.3.1      2018-12-21 [1] CRAN (R 3.6.0)                    
   readxl        1.3.1      2019-03-13 [1] CRAN (R 3.6.0)                    
   remotes       2.0.4      2019-04-10 [1] CRAN (R 3.6.0)                    
   rlang         0.4.0.9000 2019-08-03 [1] Github (r-lib/rlang@b0905db)      
   rmarkdown     1.16       2019-10-01 [1] CRAN (R 3.6.0)                    
   rprojroot     1.3-2      2018-01-03 [1] CRAN (R 3.6.0)                    
   rstudioapi    0.10       2019-03-19 [1] CRAN (R 3.6.0)                    
   rvest         0.3.4      2019-05-15 [1] CRAN (R 3.6.0)                    
   scales        1.0.0      2018-08-09 [1] CRAN (R 3.6.0)                    
   selectr       0.4-1      2018-04-06 [1] CRAN (R 3.6.0)                    
   sessioninfo   1.1.1      2018-11-05 [1] CRAN (R 3.6.0)                    
   stringi       1.4.3      2019-03-12 [1] CRAN (R 3.6.0)                    
   stringr     * 1.4.0      2019-02-10 [1] CRAN (R 3.6.0)                    
   testthat      2.2.1      2019-07-25 [1] CRAN (R 3.6.0)                    
   tibble      * 2.1.3      2019-06-06 [1] CRAN (R 3.6.0)                    
   tidyr       * 1.0.0      2019-09-11 [1] CRAN (R 3.6.0)                    
   tidyselect    0.2.5      2018-10-11 [1] CRAN (R 3.6.0)                    
   tidyverse   * 1.2.1      2017-11-14 [1] CRAN (R 3.6.0)                    
   usethis       1.5.0      2019-04-07 [1] CRAN (R 3.6.0)                    
   vctrs         0.2.0.9000 2019-08-03 [1] Github (r-lib/vctrs@11c34ae)      
   withr         2.1.2      2018-03-15 [1] CRAN (R 3.6.0)                    
   xaringan      0.9        2019-03-06 [1] CRAN (R 3.6.0)                    
   xfun          0.10       2019-10-01 [1] CRAN (R 3.6.0)                    
   xml2          1.2.0      2018-01-24 [1] CRAN (R 3.6.0)                    
   yaml          2.2.0      2018-07-25 [1] CRAN (R 3.6.0)                    
   zeallot       0.1.0      2018-01-28 [1] CRAN (R 3.6.0)                    
  
  [1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library
```

:::

These slides are licensed under <center><a href="https://creativecommons.org/licenses/by-sa/3.0/au/"><img src="images/cc.svg" style="height:2em;"/><img src="images/by.svg" style="height:2em;"/><img src="images/sa.svg" style="height:2em;"/></a></center>