genDA

class: split-40 hide-slide-number with-thick-border border-white
background-image: url("bkg/bg1.png")
background-size: cover

]]
.column.shade_main[.content.vmiddle[

# .large[**genDA**]
## A generalised discriminant analysis algorithm
## for mixed-response type data

### Sarah Romanes sarah_romanes

### Young Stats Conference, 2019

### &nbsp;bit.ly/SR-YoungStats

]]

---

class: split-70 hide-slide-number
background-image: url("bkg/bg3.png")
background-size: cover

.column.slide-in-left[
.sliderbox.vmiddle.shade_main.center[
.font5[Discriminant Analysis]]]
.column[
]

---
class: split-two white

.column.bg-main2[.content[

# **What is Discriminant Analysis?**

### Discriminant Analysis (Fisher, 1936) is a ML technique that seeks to find a linear combination of features that separates classes of objects.

### It *strictly* assumes the conditional distribution of the data, given class grouping, is .orange[multivariate normal].

### Available through `MASS` package in with functions `lda` (common covariance) and `qda`.

]]
.column.bg-main5[.content.vmiddle.center[

]]

---

.column.bg-main2[.content[

# **Advantages**

## Intuitive, and easy to use.

## Describes data generating process as well as provide a classifier for new points.

# .orange[**but...**]

]]
.column.bg-main5[.content[

# **Disadvantages**

## Does not work when `$p > n$` due to MLE covariance matrix estimates being singular.

## Does not work for non-Normal data types.

]]

---

class: middle center hide-slide-number
background-image: url("bkg/bg2.png")
background-size: cover

# .black[]

---

class: split-60 hide-slide-number
background-image: url("bkg/bg2.png")
background-size: cover

.column.slide-in-right[
.sliderbox.vmiddle.shade_main.center[
.font5[ for GLLVMs ]]]
.column[
]

.column.slide-in-left[
.sliderbox.shade_main.center[
.font5[Swap multivariate Normals ]]]
.column[
]

---

.column.bg-main2[.content[

## ** Generalised Linear Latent Variable Models**

### Let `$\mathbf{Y}$` denote a `$n \times m$` response matrix. 
### Let `$\mathbf{X}$` denote a design matrix of class information for each observation.

### In a GLLVM, the mean responses are regressed against a .orange[vector of] `$\color{orange}{d \ll m}$` .orange[latent variables] `$\color{orange}{\mathbf{u}_i}$` along with the class covariates.

]]
.column[.content.vmiddle.center[

`$$g(\mu_{ij}) = \eta_{ij} = \beta_0 + \mathbf{x}_i^T\boldsymbol{\beta}_j  \color{orange}{ + \mathbf{u}_i^T\boldsymbol{\lambda}_j }$$`

]]

### .black[*where*]

]]

---
class: split-two

.column.bg-main2[.content[

## ** Generalised Linear Latent Variable Models**

### To complete the formulation we assume that conditional on the latent variables `$\mathbf{u}_i$` and parameter vector `$\boldsymbol{\Psi}$`, the responses are independent observations from the exponential family of distributions.

### .orange[**We can capture differing response types through altering the exponential family representation for each column as needed**]

]]
.column[.content.vmiddle.center[

`$$p(y_{ij}| \mathbf{u}_i, \boldsymbol{\Psi}) = \left\{ \frac{y_{ij} a_j(\eta_{ij}) - b_j(\eta_{ij}) }{\phi_j} + c_j(y_{ij}, \phi_j)  \right\}$$`

]]

---

.column.bg-main2[.content[

# ** GLLVMs for DA **

## .orange[*Common Covariance*]

## In the common covariance model, class information is captured in `$\mathbf{x}_i$`.

]]
.column[.content.vmiddle.center[

]]
---

.column.bg-main2[.content[

# ** GLLVMs for DA **

## .orange[*Separate Covariance*]

## We can model .orange[differing] correlation structures by fitting multiple GLLVMs for different classes.

]]
.column[.content.vmiddle.center[

]]

---

.column.bg-main2[.content[

# ** Estimating GLLVMs **

## To obtain estimates for `$\boldsymbol{\Psi}$`, we must first marginalise over the latent variables `$\mathbf{u}_i$` as they are unobserved.

]]
.column[.content.vmiddle.center[

`$$\begin{align}
      \ell(\boldsymbol{\Psi}) & = \sum_{i=1}^n \log \left\{p(\mathbf{y}_i, \boldsymbol{\Psi}) \right\} \\
                              & = \sum_{i=1}^n \log \left( \int \prod_{i=1}^{m} p(y_{ij}| \mathbf{u}_i,\boldsymbol{\Psi}) p(\mathbf{u}_i) d \mathbf{u}_i \right)
   \end{align}$$`
]]

]]

---

.column.bg-main2[.content[

# ** Estimating GLLVMs **

## To obtain estimates for `$\boldsymbol{\Psi}$`, we must first marginalise over the latent variables `$\mathbf{u}_i$` as they are unobserved.

## However, the integral needed to obtain the marginal likelihood proves to be intractable.

]]
.column[.content.vmiddle.center[

]]

---

.column.bg-main2[.content[

# ** Estimating GLLVMs **

## To obtain estimates for `$\boldsymbol{\Psi}$`, we must first marginalise over the latent variables `$\mathbf{u}_i$` as they are unobserved.

## However, the integral needed to obtain the marginal likelihood proves to be intractable. .orange[How do we overcome this?]

]]
.column[.content.vmiddle.center[

]]
---
class: center bg-main2

# Variational Approximations to the rescue!

<img src="images/superman.gif", width="50%">
---

class: middle center hide-slide-number
background-image: url("images/VA.png")
background-size: cover

---
class: split-two

.column.bg-main2[.content[

# **VA for GLLVMs **

### Hui et. al. (2017) show that it is optimal to pick `$q(\mathbf{u}_i) \sim \mathcal{N}(\boldsymbol{\mu}_i, \boldsymbol{\Sigma}_i)$`, leading to the following expression for the approximate marginal log likelihood.

]]
.column[.content.vmiddle.center[

`$$\begin{align}
      \underline{\ell}(\boldsymbol{\Psi}, \boldsymbol{\xi}) & = \sum_{i=1}^n \sum_{j=1}^{m}  \left\{ \frac{y_{ij} a(\widetilde{\eta}_{ij}) - E_q\{b(\eta_{ij}) \}}{\phi_j} + c(y_{ij}, \phi_j)  \right\} \\
                              & + \frac{1}{2} \sum_{i=1}^{n} ( \log \det(\boldsymbol{\Sigma}_i) - \mbox{tr}(\boldsymbol{\Sigma}_i) - \boldsymbol{\mu}_i^T\boldsymbol{\mu}_i ),
   \end{align}$$`
]]

### .black[*where*]

.medium[.black[
`$$\widetilde{\eta}_{ij} = \beta_{0j} + \mathbf{x}_i^T\boldsymbol{\beta}_j  + \boldsymbol{\mu}_i^T\boldsymbol{\lambda}_j$$`  ]]

]]

---

.column.bg-main2[.content[

# **VA for GLLVMs **

### Hui et al. (2017) show that it is optimal to pick `$q(\mathbf{u}_i) \sim \mathcal{N}(\boldsymbol{\mu}_i, \boldsymbol{\Sigma}_i)$`, leading to the following expression for the variational lower bound for the approximate marginal log likelihood.

### .orange[Maximising the approximate lower bound is equivalent to minimising KL divergence between approximation and true marginal likelihood.]

]]
.column[.content.vmiddle.center[

### .black[*where*]

.medium[.black[
`$$\widetilde{\eta}_{ij} = \beta_{0j} + \mathbf{x}_i^T\boldsymbol{\beta}_j  + \boldsymbol{\mu}_i^T\boldsymbol{\lambda}_j$$`  ]]

]]

---

.column.bg-main2[.content[

# **Implementation**

### In the optimisation process, we utilise the `TMB` package, which allows us to perform .orange[automatic differentiation], implemented in `C++`

]]
.column[.content.vmiddle.center[

<img src="images/AutomaticDifferentiationNutshell.png", width="85%">
]]

---

.column.bg-main2[.content[

# **Implementation**

### In the optimisation process, we utilise the `TMB` package, which allows us to perform .orange[automatic differentiation], implemented in `C++`.

### Automatic differentiation allows us to only define the objective function, and .orange[*automatically*] calculates the gradient function for us. It does this by breaking up the function into basic functions and applying the chain rule.

]]
.column[.content.vmiddle.center[

<img src="images/frodo.gif", width="70%">
 
 .pink[*Disclosure: unlike Mathematica it will not return a symbolic gradient formula for you.*]

]]

---

class: split-70 hide-slide-number
background-image: url("bkg/bg3.png")
background-size: cover

.column.slide-in-left[
.sliderbox.vmiddle.shade_main.center[
.font5[genDA in Action ]]]
.column[
]

---

.column.bg-main2[.content[

# **Urban Cover data**

### The study area is an urban area in Deerfield Beach, FL, USA, with a 30cm resolution colour infrared aerial orthoimagery of the study area acquired.

### Contains 9 different types of landcover.

### Data consists of .orange[n = 168] image segments to be classified with .orange[m = 147] features associated with each image segment such as Area, Brightness, etc measured at different resolutions.

]]
.column.bg-main1[.content.vmiddle.center[

##### .orange[Source: Johnson, 2013].

]]

---

.column.bg-main2[.content[

# **Urban Cover data**

# Data is .orange[*highly correlated*] with a noticible difference in correlation structure between classes.

]]
.column.bg-main1[.content.vmiddle.center[

]]

---

.column.bg-main2[.content[

# **Urban Cover data**

## EDA shows that data is .orange[not normal], with a mix of positive skew, Gaussian, and count data.

# .orange[**Let's give genDA a go on this data!**]

]]
.column.bg-main1[.content.vmiddle[

```r
# A tibble: 168 x 148
 class BrdIndx Area Round Bright Compact ShpIndx 
 <fct> <dbl> <int> <dbl> <dbl> <dbl> <dbl> 
 1 car 1.27 91 0.97 231. 1.39 1.47
 2 concrete 2.36 241 1.56 216. 2.46 2.51 
 3 concrete 2.12 266 1.47 232. 2.07 2.21
 4 concrete 2.42 399 1.28 230. 2.49 2.73 
 5 concrete 2.15 944 1.73 193. 2.28 4.1 
 6 tree 3.11 169 1.47 172. 2.49 3.35 
 7 car 1.2 44 0.79 209. 1.14 1.36 
 8 car 1 88 0.22 235. 1.11 1.12 
 9 building 1.59 1737 0.67 220. 1.3 1.64 
10 tree 2.37 153 1.3 120. 2.85 2.59
```

]]

---

.column.bg-main2[.content[

# **Urban Cover data**

```r
*fit <- genDA(y = data, class = class,
* num.lv = 2, common.covariance = FALSE)

prediction <- predict(fit, newdata = newdata) 
```

]]
.column.bg-main1[.content.vmiddle.center[

### .black[The genDA algorithm can be fit with the `genDA` function, and the generic S3 `predict` method can be used to predict new data points given a trained genDA model.]

#### .orange[(*Note: I expect the syntax will change slightly as I further improve the R package*)]

]]

---

# .black[100 Trial, 5 Fold CV]

---
class: split-two white

]]

.column.bg-main2[.content.center[

# genDA

## .white[]

### [sarahromanes/genDA](https://github.com/sarahromanes/genDA)

]]

.column.bg-main8[.content.center[

# **Acknowledgements**

### Slides are made in *rmarkdown* using *xaringan* (Yihui Xie) and *ninja* theme (Emi Tanaka).

### Thank you to SSA NSW for the YSC Travel Grant to attend this conference.
<img src="images/SSA_logo.png", width="50%">

]]